This notebook contains material for CBE 20258 Numerical and Statistical Analysis taught at the University of Notre Dame. (c) Professors Alexander Dowling, Ryan McClarren, and Yamil Colón. This collection of notebooks cbe-xx258 is available on Github.

In [ ]:
# IMPORT DATA FILES USED BY THIS NOTEBOOK
import os,  requests

file_links = [("data/fifth_republic.txt", "https://ndcbe.github.io/cbe-xx258/data/fifth_republic.txt"),
    ("data/hats.txt", "https://ndcbe.github.io/cbe-xx258/data/hats.txt")]

# This cell has been added by nbpages. Run this cell to download data files required for this notebook.

for filepath, fileurl in file_links:
    stem, filename = os.path.split(filepath)
    if stem:
        if not os.path.exists(stem):
            os.mkdir(stem)
    if not os.path.isfile(filepath):
        with open(filepath, 'wb') as f:
            response = requests.get(fileurl)
            f.write(response.content)

1.4 Functions, Scoping, and Other Fun Stuff

Reference: Chapter 3 of Computational Nuclear Engineering and Radiological Science Using Python, R. McClarren (2018)

1.4.1 Learning Objectives

After studying this notebook, completing the activities, and asking questions in class, you should be able to:

1.4.2 Functions

Motivating Example

Why use functions? We want to write, debug, and test code once and then reuse as much as possible.

In a few class sessions, we'll formulate mass balances as linear systems and solve them using Python. But for now, let's just consider a problem you would expect to see in math class:

We want to solve the linear system,

$$\mathrm{Eqn.~1}:\quad 4.5 x + 3 y = 10.5\\\mathrm{Eqn.~2}:\quad 1.5 x + 3 y = 7.5.$$

One way to do this is using Python as calculator. In the comments below, we walk through the steps.

In [1]:
"""python code to solve 
4.5 x + 3 y = 10.5
1.5 x + 3 y = 7.5
by solving the second equation for y first,
and then solving for x"""
#step 1 solve for y, multiply equation 2 by -3, 
## and add to first equation
LHS_coefficient = -3*3 + 3 #the coefficient for y
RHS = -3*7.5 + 10.5 #the right-hand side

print('LHS_coefficient:',LHS_coefficient)
print('RHS:',RHS)

Mathematical, we started by multiplying equation 2,

$$1.5 x + 3 y = 7.5$$

by -3,

$$(-3) \times 1.5 x + (-3) \times 3 y = (-3) \times 7.5$$

and then added this scaled equation 2 to equation 1, giving:

$$(4.5 - 3 \times 1.5) x + (3 - 3 \times 3) y = 10.5 - 3 \times 7.5$$

Notice that our choice of scaling equation 2 by -3 means that the coefficient for x becomes zero after addition addition. The coefficient for $y$ is LHS_coefficient is our code. RHS is the right hand side of the new equation.

In [2]:
#now divide right-hand side by left-hand side coefficient
y = RHS / LHS_coefficient
#plug y into first equation
x = (10.5 - 3*y)/4.5 
#print the solution, note \n produces a linebreak
print("The solution to:\n4.5 x + 3 y = 10.5\n1.5 x + 3 y = 7.5\n is x =",
      x,"y=",y)
The solution to:
4.5 x + 3 y = 10.5
1.5 x + 3 y = 7.5
 is x = 1.0 y= 2.0

How to extend this code to another linear system?

Let's define a function that will solve the system for (almost) any coefficients and right-hand side.

I'll define such a function to solve $$a_1 x + b_1 y = c_1\\ a_2 x + b_2 y = c_2.$$

Home Activity: Write pseudocode to generalize the steps from the motivating example to solve the linear system with coefficients a1, a2, b1, b2, c1, and c2. Turn in your pseudocode via Gradescope. We will only look at your pseudocode for completeness and compliance with the formatting guidelines, not correctness.
Class Activity: Discuss your pseudocode with a partner. Give one compliment and one suggestion to partner's pseudocode. We'll then regroup and write it together as a class.

Below is a function that solves (most) 2x2 linear systems. Take a few minutes to study the code below. Specifically:

  • Notice the function has seven inputs. The first six are the coefficients. The seventh, LOUD, is followed by =False. This sets the default value of LOUD to false.
  • The input LOUD toggles on/off a print statement.
  • This function has a long comment string at the top. It includes a brief description, then a list of inputs (arguments) and finally a list of outputs (returns). All of the functions you write in this class must be commented in the same style.
In [3]:
def two_by_two_solver(a1,b1,c1,a2,b2,c2, LOUD=False):
    """Calculate the solution of the system 
    a1 x + b1 y = c1, 
    a2 x + b2 y = c2

    Args:
        a1: x coefficient in first equation (cannot be zero)
        b1: y coefficient in first equation
        c1: right-hand side in first equation
        a2: x coefficient in second equation 
        b2: y coefficient in second equation 
        c2: right-hand side in second equation
        LOUD: boolean that decides whether to print out the answer
        
    Returns:
        list containing the solution in the format [x,y]
    """
    #step one, eliminate x from the second equation by 
    #multiplying first equation by -a2/a1
    #and then adding it to second equation
    new_b2 = b2 - a2/a1*b1
    new_c2 = c2 - a2/a1*c1
    #solve the new equation 2
    y = new_c2/new_b2
    #plug y into original equation 1
    x = (c1-b1*y)/a1
    
    if (LOUD):
        print("The solution to:\n",a1,"x +",b1,"y =",c1,
              "\n",a2,"x +",b2,"y =",c2,"\n is x =",x,"y=",y)
    return [x,y]

We can call this function for the problem above by typing

In [4]:
two_by_two_solver(4.5,3,10.5,1.5,3,7.5,True)
The solution to:
 4.5 x + 3 y = 10.5 
 1.5 x + 3 y = 7.5 
 is x = 1.0 y= 2.0
Out[4]:
[1.0, 2.0]

We can also solve other systems, including simple ones

In [5]:
two_by_two_solver(1,0,3,0,1,2,True)
The solution to:
 1 x + 0 y = 3 
 0 x + 1 y = 2 
 is x = 3.0 y= 2.0
Out[5]:
[3.0, 2.0]

We can't solve systems where $a_1$ is zero because our function divides by $a_1$:

In [6]:
two_by_two_solver(0,1,2,1,0,3,True)
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-6-5036cb57389a> in <module>
----> 1 two_by_two_solver(0,1,2,1,0,3,True)

<ipython-input-3-bade12b85c5d> in two_by_two_solver(a1, b1, c1, a2, b2, c2, LOUD)
     19     #multiplying first equation by -a2/a1
     20     #and then adding it to second equation
---> 21     new_b2 = b2 - a2/a1*b1
     22     new_c2 = c2 - a2/a1*c1
     23     #solve the new equation 2

ZeroDivisionError: division by zero

1.4.2.1 Calling functions

Approach above: give inputs in order

We called the function two_by_two_solver by listing out the arguments in the order that it expects them a1, b1, c1, a2, b2, c2, LOUD.

Another option: use keywords

Python allows you to call them in any order, as long as you are explicit in what goes where.

In [ ]:
two_by_two_solver(a1 = 4.5, b1 = 3, 
                  a2 = 1.5, b2 = 3, 
                  c1 = 10.5, c2 = 7.5, LOUD = True)

It is often a good idea to call a function explicitly (with keywords). That way if you mess up the order of the arguments, it does not matter.

Notice that in the function definition, the argument LOUD has =False after it. This indicates that if the function is called without a value for LOUD, it assumes the caller does not what the function to "be loud".

In other words, False is the default for argument LOUD.

In [ ]:
two_by_two_solver(a1 = 4.5, b1 = 3, a2 = 1.5, 
                  b2 = 3, c1 = 10.5, c2 = 7.5)

Notice that it did not print out it's spiel about the system.

In [ ]:
two_by_two_solver(1,1,2,a2 = 1, c2 = 0, b2 = 3)

1.4.2.2 Return Values

At the end of the function we have a return statement. This tells python what the function is returning to the caller. In this case we return a list that has the solution for $x$ and $y$. We can store this in a new variable, or do whatever we like with it.

In [ ]:
answer = two_by_two_solver(a1 = 4.5, b1 = 3, a2 = 1.5, 
                           b2 = 3, c1 = 10.5, c2 = 7.5)
x = answer[0] #store in the variable x the first value in the list answer
y = answer[1] #store in the variable y the first value in the list answer
print("The list",answer,"contains",x,"and",y)
Home Activity: Solve the linear system given below and store the answers in `my_x` and `my_y`. This will be auto-graded.
$$ 2 x -1 y = 3\\ -4 x + 3 y = 0.$$
In [ ]:
# YOUR SOLUTION HERE
In [ ]:
 

We can do even fancier things, if we are so bold

In [ ]:
#just get x
x = two_by_two_solver(a1 = 4.5, b1 = 3, a2 = 1.5, 
                      b2 = 3, c1 = 10.5, c2 = 7.5)[0]
print("x =",x)

#assign variables to the output on the fly
x,y = two_by_two_solver(a1 = 4.5, b1 = 3, a2 = 1.5, 
                        b2 = 3, c1 = 10.5, c2 = 7.5)
print("x =",x,"y =",y)

These examples are more advanced and they are designed to show you some of the neat tricks you can do in python.

1.4.2.3 Docstrings and Help

Our 2x2 solver code had a long, and pretty detailed comment at the beginning of it. This is called a docstring and it is printed by a user by typing

In [ ]:
help(two_by_two_solver)

The point of this long comment is to communicate to other people:

  1. The main idea behind the function.
  2. What the function expects from you. Inputs.
  3. What the function gives you. Outputs.

In this example we can see that we need to provide at least 6 numbers, and possibly an optional boolean.

It is good programming practice to include good docstrings. It will be mandatory for any code you turn in this class.

Let's look at the docstring for some members of the math module and the random module.

In [ ]:
import math
help(math.fabs)
In [ ]:
import random
help(random.uniform)

We don't have the source code for these functions in front of us, but if we want to know how to call them (and we didn't want to Google® them), the docstrings tell us what to do.

You may wonder why those docstrings are a bit different that the one I used. The format for mine is derived from the <a, href="http://google-styleguide.googlecode.com/svn/trunk/pyguide.html#Comments"> Google coding standards for python docstrings</a>.

In [ ]:
def nothing():
    #this function does nothing
    #docstring?
    return 0
help(nothing)

1.4.3 Scope

When we call a function, it carves out in the computer's memory its own space. Variables that live in the special space, known as local scope are completely different than those in other parts of the program. Here's a simple, but illustrative example:

Home Activity: Before you run the code below, predict the output. This is good practice for exam questions, where you'll be asked to predict the output of Python code without access to a computer.

Home Activity - Your Predictions

x =

new_x =

y =

new_y =

In [ ]:
def scope_demonstration(input_variable):
    ''' A simple demonstration of scoping rules
    
    Args:
        input_variable: a string or number
        
    Returns:
        x: --REDACTED TO NOT GIVE AWAY THE ACTIVITY ANSWER --
    '''
    x = input_variable*3
    return x

#now call the function after defining some variables
x = "oui "
y = "no "

new_x = scope_demonstration(x)
new_y = scope_demonstration(y)
print("x =",x,"\nnew_x =",new_x)
print("y =",y,"\nnew_y =",new_y)

When I call scope_demonstration it creates its own memory space and any variable I create in there is different than in the rest of the program, even if the variables have the same name.

There are many subleties in scoping rules, but this example outlines just about everything that you'll need to know as a neophyte programmer.

A variable declared at the top level Python file, i.e., outside a function, is in the global scope for that file. This is easiest to see in an example:

In [ ]:
my_var = 1

def print_my_var():
    ''' Another scope demonstration
    
    Args:
        Nothing
        
    Returns:
        Nothing
        
    Other:
        prints `my_var`
    
    '''
    print(my_var)
    return None

print_my_var()

We were able to access my_var inside the function even though it was not an input to a function.

The next natural guess is can we modify my_var with a function? Let's see.

In [ ]:
def change_my_var():
    ''' A third scope demonstration
    
    Args:
        Nothing
        
    Returns:
        Nothing
        
    Other:
        Attempts to change `my_var, prints to screen
    '''
    
    my_var = 3*my_var + 2
    print("my_var = ",my_var," inside the function")
    
    return None

change_my_var()
print("my_var =",my_var," outside the function")

We got an error message! In Python, we can only access my_var inside functions. We cannot change it inside functions. If we wanted to change my_var, we would need to do the following:

In [ ]:
def change_my_var2():
    ''' A third scope demonstration
    
    Args:
        Nothing
        
    Returns:
        New value of my_var
        
    '''
    
    return 3*my_var + 2

my_var = change_my_var2()
print("Now my_var =",my_var)
Class Activity: Take 30 seconds to predict the output of the code below without running it. Below are multiple choice answers.
In [ ]:
def my_func1(x,y):
    ''' A simple function to demonstrate the nuances of scope

    Arguments:
        x: scalar real number
        y: scalar real number
        
    Returns:
        z: scalar real number

    '''
    z = x + y
    x = 3
    y = z - x + 1
    return z

# Run the function
x = 1
y = 1
z = my_func1(x,y)
# Print values of x, y, and z to screen
print("x =",x,"  y =",y,"  z =",z)

Class Activity multiple choice answers:

  1. x=1, y=1, z=2
  2. x=1, y=1, z=3
  3. x=3, y=0, z=2
  4. x=3, y=0, x=3
  5. None of these
In [ ]:
### Create your answer here.

1.4.4 Recursion

The idea behind recursion is that a function can call itself. That is really it. That capability can allow some neat tricks, but in practice there is often a faster way of doing things than using recursion.

Home Activity: Write pseudcode to calculate the factorial of integer x two ways: i) without recursion and ii) with recursion. Submit your pseudocode via Gradescope. We'll only assess for completeness (making an attempt, following guidelines) and not accuracy. Be sure to bring your pseudocode to class.
Class Activity: Explain your pseudocode to a partner. One person explains without recursion and the other explains with recurison.

Below are two functions that calculate factorial with and without recurision.

In [7]:
def factorial(n, prev=1):
    if not((n==1) or (n==0)):
        prev = factorial(n-1,prev)*n
    return prev

def factorial_no_recursion(n):
    output = 1;
    #can skip 1 because x*1 = 1
    for i in range(2,n+1):
        output = output*i
    return output
x = int(input("Enter an integer: "))
print(x,"! =",factorial(x))
print(x,"! =",factorial_no_recursion(x))
Enter an integer: 212
212 ! = 4733702182912325971198157282773458528972111665671644583063162654840567216299326200333597974632020795344694044141162288574741860330707871653991802372413420954892019572846468089404909755852192508097446724647826768577878987213960691804730882223315446309650598202756704313010742315578131345078364709758529795655446581758477730600169824143256656411069775872000000000000000000000000000000000000000000000000000
212 ! = 4733702182912325971198157282773458528972111665671644583063162654840567216299326200333597974632020795344694044141162288574741860330707871653991802372413420954892019572846468089404909755852192508097446724647826768577878987213960691804730882223315446309650598202756704313010742315578131345078364709758529795655446581758477730600169824143256656411069775872000000000000000000000000000000000000000000000000000

Let's see which is faster. This will use some ipython magic commands for timing execution. We'll compute the factorials of 0 through 20, 100,000 times.

In [8]:
%%time
for times in range(10**5):
    for n in range(21):
        factorial(n)
CPU times: user 3.19 s, sys: 26.5 ms, total: 3.22 s
Wall time: 3.22 s
In [9]:
%%time
for times in range(10**5):
    for n in range(21):
        factorial_no_recursion(n)
CPU times: user 1.38 s, sys: 12 ms, total: 1.39 s
Wall time: 1.39 s

Side-tangent: wall versus CPU user versus CPU clock time: https://stackoverflow.com/questions/7335920/what-specifically-are-wall-clock-time-user-cpu-time-and-system-cpu-time-in-uni (Side-tangents will not appear in tests or assignments, but are given to satisfy your curiosity.)

The no recursion version, while not as neat, is nearly 50% faster. It is also possible to have too many levels of recursion.

In [10]:
import sys
sys.setrecursionlimit(1000) # change this to answer the home activity question
x = 1000
#this won't work and prints ~1000 errors
#the errors are not repeated here
print(x,"! =",factorial(x))
---------------------------------------------------------------------------
RecursionError                            Traceback (most recent call last)
<ipython-input-10-f352ad03aa45> in <module>
      4 #this won't work and prints ~1000 errors
      5 #the errors are not repeated here
----> 6 print(x,"! =",factorial(x))

<ipython-input-7-330ea1369dd6> in factorial(n, prev)
      1 def factorial(n, prev=1):
      2     if not((n==1) or (n==0)):
----> 3         prev = factorial(n-1,prev)*n
      4     return prev
      5 

... last 1 frames repeated, from the frame below ...

<ipython-input-7-330ea1369dd6> in factorial(n, prev)
      1 def factorial(n, prev=1):
      2     if not((n==1) or (n==0)):
----> 3         prev = factorial(n-1,prev)*n
      4     return prev
      5 

RecursionError: maximum recursion depth exceeded in comparison
Home Activity: What is the smallest recursion limit for which factorial(1000) works?
In [11]:
x = 1000
#this works
answer = factorial_no_recursion(x)
print(x,"! =",answer)
1000 ! = 402387260077093773543702433923003985719374864210714632543799910429938512398629020592044208486969404800479988610197196058631666872994808558901323829669944590997424504087073759918823627727188732519779505950995276120874975462497043601418278094646496291056393887437886487337119181045825783647849977012476632889835955735432513185323958463075557409114262417474349347553428646576611667797396668820291207379143853719588249808126867838374559731746136085379534524221586593201928090878297308431392844403281231558611036976801357304216168747609675871348312025478589320767169132448426236131412508780208000261683151027341827977704784635868170164365024153691398281264810213092761244896359928705114964975419909342221566832572080821333186116811553615836546984046708975602900950537616475847728421889679646244945160765353408198901385442487984959953319101723355556602139450399736280750137837615307127761926849034352625200015888535147331611702103968175921510907788019393178114194545257223865541461062892187960223838971476088506276862967146674697562911234082439208160153780889893964518263243671616762179168909779911903754031274622289988005195444414282012187361745992642956581746628302955570299024324153181617210465832036786906117260158783520751516284225540265170483304226143974286933061690897968482590125458327168226458066526769958652682272807075781391858178889652208164348344825993266043367660176999612831860788386150279465955131156552036093988180612138558600301435694527224206344631797460594682573103790084024432438465657245014402821885252470935190620929023136493273497565513958720559654228749774011413346962715422845862377387538230483865688976461927383814900140767310446640259899490222221765904339901886018566526485061799702356193897017860040811889729918311021171229845901641921068884387121855646124960798722908519296819372388642614839657382291123125024186649353143970137428531926649875337218940694281434118520158014123344828015051399694290153483077644569099073152433278288269864602789864321139083506217095002597389863554277196742822248757586765752344220207573630569498825087968928162753848863396909959826280956121450994871701244516461260379029309120889086942028510640182154399457156805941872748998094254742173582401063677404595741785160829230135358081840096996372524230560855903700624271243416909004153690105933983835777939410970027753472000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

By the way, I'm surprised it is giving the correct answer.

Sometimes we have to define many functions and don't want to have one giant source file (in fact this is good programming practice). I've created a file called sphere.py in the same directory as this notebook. This file defines two functions volume and surface_area that compute the volume and surface area of a sphere. Because the file is called sphere.py we can import those functions using import sphere. The text of sphere.py is def volume(radius): """compute volume of a sphere

Args:
radius: float giving the radius of the sphere

Returns:
volume of the sphere as a float
"""
return 4.0/3.0*math.pi*radius**3

def surface_area(radius): """compute surface area of a sphere

Args:
radius: float giving the radius of the sphere

Returns:
surface area of the sphere as a float
"""
return 4.0*math.pi*radius**2

</code>

I can use the help function to tell me about the module:

In [12]:
import sphere
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-12-6ee0c90e35fb> in <module>
----> 1 import sphere

ModuleNotFoundError: No module named 'sphere'

Now that I have imported the module, we can can the help to see the docstring.

In [ ]:
help(sphere)
In [ ]:
r = 1.0
print("The volume of a sphere of radius",r,"cm is",
      sphere.volume(r),"cm**3")
print("The surface area of a sphere of radius",r,"cm is",
      sphere.surface_area(r),"cm**2")

1.4.6 Files

It is very easy to read in text files in python. The file fifth_republic.txt, lists the presidents of France's fifth republic. It is in the same folder as this notebook.

In Python, we can read it in very simply:

In [14]:
file = open('./data/fifth_republic.txt', 'r') 
#open fifth_republic.txt for reading ('r')
for line in file:
    # Repeat the first 5 characters 3 times
    print(line[0:5]*3)
file.close()
CharlCharlCharl
GeorgGeorgGeorg
ValérValérValér
FrançFrançFranç
JacquJacquJacqu
NicolNicolNicol
FrançFrançFranç
EmmanEmmanEmman

Notice how the for loop can iterate through each line of the file. You can also read a line at a time.

In [15]:
file = open('./data/fifth_republic.txt', 'r')
#open fifth_republic.txt for reading ('r')

# Read the first line
first_line = file.readline()

# Read the second line
second_line = file.readline()
print(first_line)
print(second_line)
file.close()
Charles de Gaulle

Georges Pompidou

In [16]:
help(file.readline)
Help on built-in function readline:

readline(size=-1, /) method of _io.TextIOWrapper instance
    Read until newline or EOF.
    
    Returns an empty string if EOF is hit immediately.

You can also easily write to a file.

In [17]:
writeFile = open("./data/hats.txt","w") 
#open hats.txt to write (clobber if it exists)
hats = ["fedora","trilby","porkpie",
        "tam o'shanter","Phrygian cap","Beefeaters' hat","sombrero"]
for hat in hats:
    writeFile.write(hat + "\n") #add the endline
writeFile.close()

#now open file and print
readFile = open("./data/hats.txt","r")
for line in readFile:
    print(line)
readFile.close()
fedora

trilby

porkpie

tam o'shanter

Phrygian cap

Beefeaters' hat

sombrero

We can also use enumerate with a text file:

In [18]:
import csv
file = open('./data/fifth_republic.txt', 'r')
for i,row in enumerate(file):
    print(row)
Charles de Gaulle

Georges Pompidou

Valéry Giscard d'Estaing

François Mitterrand

Jacques Chirac

Nicolas Sarkozy

François Hollande

Emmanuel Macron

1.4.7 Example: High/Low Guess My Number Game

1.4.7.1 Instructions

Write pseudocode and Python function program to interactively play the following game:

  1. Generate a random integer number between 0 and 100.
  2. Ask the player to guess a number. Capture their input from the keyboard.
  3. Check if the number is valid (between 0 and 100, integer). If not, warn the player and goto Step 2.
  4. If the player’s guess is correct, they win. Print out “Congratulations” and tell them how many guesses it took. Only count valid guesses. Terminate the program.
  5. If the guess is too low, print “Guess is too low.” and goto Step 2.
  6. If the guess is too high, print “Guess is too high.” and goto Step 2.

First write the pseudocode and then implement as a function in Python. Test your program at least 2 times.

1.4.7.2 Pseudocode

1.4.7.3 Python Implementation

In [22]:
# import libraries
import numpy as np

# function to accept guess from user and determine if guess matches secret 
# number between 0 and 100
def guessing_game(rand_num=None):
    ''' Play the guess my number interactively with keyboard input
    
        Argument:
            rand_num: the random (secret) number. If None (default), the program generates
                a random integer between 0 and 100.
            
        Returns:
            Nothing
            
        Other:
            Asks user for interactive keyboard input that are integers
    
    '''
    
    # YOUR SOLUTION HERE

1.4.7.4 Testing

Test 1: No input. The function will automatically generate a secret number.

In [20]:
# YOUR SOLUTION HERE
Guess the number: 51
Guess is too high

Guess the number: 25
Guess is too high

Guess the number: 10
Guess is too low

Guess the number: 15
Guess is too low

Guess the number: 20
Guess is too low

Guess the number: 22
Guess is too high

Guess the number: 21

Congratulations!  Your guess is correct.  Number is: 21
Number of guesses: 7

Test 2: What happens if we supply an out-of-bounds (invalid) secret number?

In [23]:
# YOUR SOLUTION HERE
Secret number if too small.
In [ ]: