1.10 Testing and Debugging in Python¶
Reference: Chapter 6 of Computational Nuclear Engineering and Radiological Science Using Python, R. McClarren (2018)
1.10.1 Learning Objectives¶
After studying this notebook, completing the activities, and asking questions in class, you should be able to:
- Create
assertstatements andtry-exceptblocks of code - Embrace the "debugging mindset" and test code you write it. We strongly recommend creating a handful of unit tests for each function you write, especially for complex codes (e.g., mini-projects)
1.10.1 Testing Your Code¶
Anytime you write code, you want to make sure that it accomplishes the desired task.
Motivating Example¶
Through an example, we will demonstrate several strategies for systematically testing code.
Phenolphthalein is a base indicator (it turns a pretty purple), but in the presence of a strong base it will fade over time. This is a reversible reaction: pseudo first order in the forward (fading) direction and first order in the reverse direction. The rate expression is:
$$ \frac{d[P]}{dt} = -k_{1}[P][OH^-]+k_{2}[POH] ~, $$where $k_{1}' = k_{1} \cdot [OH^-]$ is the pseudo first order rate. $[POH]$ is only produced from this reaction, thus we have the mass balance:
$$ [POH] = [POH]\Big|_{t=0}+([P]_{0}-[P]) ~. $$With the initial conditions:
$$ [P]\Big|_{t=0}=[P]_{0}\qquad; \qquad [POH]\Big|_{t=0}=0 ~. $$We can solve these equations to get the solution:
$$ [P] = [P]_{0}\Big(\frac{k_{2}}{k_{1}'+k_{2}}+\frac{k_{1}'}{k_{1}'+k_{2}}e^{-(k_{1}'+k_{2})t}\Big) ~. $$Let's examine a Python function to calculate the solution.
import numpy as np
def P_conc(t, P0, k1prime, k2, LOUD=False):
""" Computes the concentration of phenolphthalein as a function of time
Args:
t: time, either a scalar (float) or vector (numpy array)
P0: initial concentration of phenolphthalein
k1prime: forward reaction rate constant
k2: reverse reaction rate constant
Returns:
concentration, either a scalar (float) or vector (numpy array) depening in arg. t
"""
# compute first fraction in formula
frac1 = k2/(k1prime + k2)
# compute section fraction in formula
frac2 = k1prime/(k1prime + k2)
# compute term that goes into exponent
pre_expon = -(k1prime * k2)*t
# evaluate exponent
expon = np.exp(pre_expon)
# combine together
P = P0*(frac1 + frac2*expon)
if LOUD:
print("frac1 = ",frac1)
print("frac2 = ",frac2)
print("pre_expon = ",pre_expon)
print("expon = ",expon)
print("P = ",P)
# return final value
return P
We will now go through several tests to see if our function behaves correctly. Note: There are a few intentional mistakes above to make this exercise more interesting.
1.10.1.1 Test 0: Takes Inputs as Advertised?¶
Does the function return a single value if input t is a scalar and return a numpy array if t is a vector?
# Check with scalar t
print(P_conc(0.0,1.0,1.0,1.0))
# Check with vector t
print(P_conc(np.zeros(3),1.0,1.0,1.0))
Good. The time input works as advertised.
1.10.1.2 Test 1: Initial Condition¶
First, let's check if our function returns the initial condition at $t=0$. From the formula above, we expect:
$$ [P]\Big|_{t=0} = [P]_{0}\Big(\frac{k_{2}}{k_{1}'+k_{2}}+\frac{k_{1}'}{k_{1}'+k_{2}}e^{0}\Big) = [P]_{0}\Big(\frac{k_{2}}{k_{1}'+k_{2}}+\frac{k_{1}'}{k_{1}'+k_{2}}1\Big) = [P]_0 $$Let's start by defining a function for this test case. Notice the first input to our function is a function (see Class 4).
import math
def test1(Pfunc, P0, k1prime, k2):
''' Initial condition test
Args:
Pfunc: a Python function to calculate P concentration. Expects 4 inputs.
P0: initial condition to test
k1prime: rate constant to test
k2: rate constrant to test
Returns:
Boolean (true or false) indicated if our function passed or failed the test
Side Effect:
Prints test parameters and message to screen
'''
print("Testing initial condition with following specifications:")
print("P0 = ",P0," k1prime = ",k1prime, " k2 = ",k2)
Pcalc = Pfunc(0.0, P0, k1prime, k2)
print("Calculate value of P = ",Pcalc)
if math.fabs(Pcalc - P0) < 1E-8:
print("Test passed.\n")
return True
else:
print("Something is wrong. You need to debug your code.\n")
return False
By defining this test as a function, we can quickly test several values for the input parameters.
test1(P_conc,1.0,1.0,1.0)
test1(P_conc, 5.0, 1.0, 0.001)
1.10.1.3 Test 2: Asymptotics¶
Let's examine the solution again, but consider $t \rightarrow \infty$.
$$ [P]\Big|_{t\rightarrow\infty} = [P]_{0}\Big(\frac{k_{2}}{k_{1}'+k_{2}}+\frac{k_{1}'}{k_{1}'+k_{2}}e^{-\infty}\Big) = [P]_{0}\Big(\frac{k_{2}}{k_{1}'+k_{2}}+\frac{k_{1}'}{k_{1}'+k_{2}}0\Big) = [P]_0 \Big(\frac{k_{2}}{k_{1}'+k_{2}} \Big) $$Let's start by plotting $[P]$ as a function of time and see if it converges:
import matplotlib.pyplot as plt
#this line is only needed in iPython notebooks
%matplotlib inline
def test2(Pfunc, P0, k1prime, k2):
'''Asymptotics test
Args:
Pfunc: a Python function to calculate P concentration. Expects 4 inputs.
P0: initial condition to test
k1prime: rate constant to test
k2: rate constrant to test
Returns:
Boolean (true or false) indicated if our function passed or failed the test
Side Effect:
Creates a plot. Prints if test passed or failed.
'''
# Generate 20 points from 10^-2 to 10^4
t = np.logspace(-1,4,20)
P = P_conc(t, P0, k1prime, k2)
P_infty = P0*k2/(k1prime + k2)
# plot calculated value
plt.semilogx(t,P,'o-',label="function")
# plot calculated asymptotic value
plt.semilogx(t,P_infty*np.ones(len(t)),color='r',linestyle="--",label="$([P]_0 \cdot k_{2}) / (k_{1}'+k_{2})$")
# label plot
plt.xlabel('Time (units not given)')
plt.ylabel('$[P]$ (units not given)')
plt.title("$[P]_0$ = " + str(P0) + " $k_1'$ = " + str(k1prime) + " $k_2$ = " + str(k2))
plt.legend()
plt.show()
# check if solution converges to asymptotic value
if math.fabs(P[-1] - P_infty) < 1E-6:
# Note: P[-1] accesses the last value in array P
print("Test Passed.")
return True
else:
print("Test Failed. Either debug code or increase time in test.")
return False
test2(P_conc, 100, 10.0, 0.001)
1.10.1.4 Test 3: Check a Few Inputs¶
We can calculate with pencil, paper and a calculator values for $[P]$ for a few sets of inputs. We can then test if our Python code is close.
For simplicity, let's consider $t=1$, $P[0] = 1$, $k_1' = 1$, $k_2 = 1$.
$$-(k_{1}'+k_{2})t = - (1 + 1) \cdot (1) = -2$$$$e^{-2} \approx 0.13534$$$$\frac{k_{2}}{k_{1}'+k_{2}} = \frac{1}{1 + 1} = 0.5$$$$\frac{k_{1}}{k_{1}'+k_{2}} = \frac{1}{1 + 1} = 0.5$$Combine it all together:
$$ [P] = [P]_{0}\Big(\frac{k_{2}}{k_{1}'+k_{2}}+\frac{k_{1}'}{k_{1}'+k_{2}}e^{-(k_{1}'+k_{2})t}\Big) = (1) \left(0.5 + 0.5 \cdot e^{-2} \right) \approx 0.56766 $$Hint: Remember, do not use an exact equality with floating point numbers. Instead, we check if the answers agree within a tolerance. See Class 1 for more details.
# Define test 3 as a function
def test3(Pfunc):
'''Test function for known solution with t=1, [P0] = 1, k1prime = 1, k2 = 1
Arg:
Pfunc: a Python function to calculate P concentration. Expects 4 inputs.
Returns:
Boolean (true or false) indicated if our function passed or failed the test
Side Effect:
Prints test parameters and message to screen
'''
# Notice how using the keyword arguments makes the code very clear
P = Pfunc(t=1.0,P0=1.0,k1prime=1.0,k2=1.0)
if math.fabs(P - 0.56766764161) < 1E-6:
print("Passed Test 3.")
return True
else:
print("Failed Test 3. Recreate your calculations and debug your code.")
return False
# Run test
test3(P_conc)
1.10 Debugging¶
Now that we've identified that our code has a bug, the process of finding that error is called debugging. Given that we ran three tests, and two of them passed we can use how the tests are different to identify where the bug might be.
We'll repeat the formula and function again:
$$ [P] = [P]_{0}\Big(\frac{k_{2}}{k_{1}'+k_{2}}+\frac{k_{1}'}{k_{1}'+k_{2}}e^{-(k_{1}'+k_{2})t}\Big) ~. $$import numpy as np
def P_conc(t, P0, k1prime, k2, LOUD=False):
""" Computes the concentration of phenolphthalein as a function of time
Args:
t: time, either a scalar (float) or vector (numpy array)
P0: initial concentration of phenolphthalein
k1prime: forward reaction rate constant
k2: reverse reaction rate constant
Returns:
concentration, either a scalar (float) or vector (numpy array) depening in arg. t
"""
# compute first fraction in formula
frac1 = k2/(k1prime + k2)
# compute section fraction in formula
frac2 = k1prime/(k1prime + k2)
# compute term that goes into exponent
pre_expon = -(k1prime * k2)*t
# evaluate exponent
expon = np.exp(pre_expon)
# combine together
P = P0*(frac1 + frac2*expon)
if LOUD:
print("frac1 = ",frac1)
print("frac2 = ",frac2)
print("pre_expon = ",pre_expon)
print("expon = ",expon)
print("P = ",P)
# return final value
return P
1.10.1 Debugging is a Mindset¶
To properly identify the bug you need to think about
- How is the code failing?
- What is the code doing correctly?
- What pieces of the code are most likely to have an error?
Moreover, your best tool for debugging is the print function. When in doubt have the code print out what happens after every line and check that with your intuition/expectation.
When in doubt, print it (all) out!
Notice that we used the LOUD keyword to help us quick toggle print statements on and off. You should always use this trick in this class.
Let's use the print statements to walk through the tests.
Test 1¶
Test 1 checks that we get the initial condition at $t=0$. We expect that prexpon evaluates to 0.0 and expon evaluates to 1.0. Let's see if this is happening:
P_conc(t=0.0, P0=1.0, k1prime = 1.0, k2 = 1.0, LOUD = True)
With expon at 1.0, Test 1 checks that frac1 and frac2 sum to 1.0, which they do. Thus, it makes sense that Test 1 passes.
1.10.2 Test 2¶
Test 2 checks for the correct properties as $t \rightarrow \infty$. Let's look at the output with t set to a huge number.
P_conc(t=1E7, P0=1.0, k1prime = 1.0, k2 = 2.0, LOUD = True)
As expected, with t set to $10^7$, expon evaluates to 0.0. This kills frac2. The calculated value for P matches frac1, which we expect.
1.10.3 Test 3¶
Now let's compare the intermediate values we calculated by hand in Test 3 to those reported below:
P_conc(t=1.0, P0=1.0, k1prime = 1.0, k2 = 1.0, LOUD = True)
import numpy as np
# Find the bug and fix this function
def P_conc_fixed(t, P0, k1prime, k2, LOUD=False):
""" Computes the concentration of phenolphthalein as a function of time
Args:
t: time, either a scalar (float) or vector (numpy array)
P0: initial concentration of phenolphthalein
k1prime: forward reaction rate constant
k2: reverse reaction rate constant
Returns:
concentration, either a scalar (float) or vector (numpy array) depening in arg. t
"""
# compute first fraction in formula
frac1 = k2/(k1prime + k2)
# compute section fraction in formula
frac2 = k1prime/(k1prime + k2)
# compute term that goes into exponent
pre_expon = -(k1prime + k2)*t
# evaluate exponent
expon = np.exp(pre_expon)
# combine together
P = P0*(frac1 + frac2*expon)
if LOUD:
print("frac1 = ",frac1)
print("frac2 = ",frac2)
print("pre_expon = ",pre_expon)
print("expon = ",expon)
print("P = ",P)
# return final value
return P
# Rerun the tests
print("\nTEST 1...")
test1(P_conc_fixed,1.0,1.0,1.0)
test1(P_conc_fixed, 5.0, 1.0, 0.001)
print("\nTEST 2...")
test2(P_conc_fixed, 100, 10.0, 0.001)
print("\nTEST 3...")
test3(P_conc_fixed)
Write your answer here:
1.10 Assertions¶
In life it pays to be assertive. The same is true in programming. Let's try an example.
test_P = P_conc(t=10.0, P0=-10, k1prime=1, k2=1)
print("With negative P0, P(t=10) =", test_P)
The code worked, but the answer is non-sensical. We'd like to tell the user that P_calc was called with an improper value. We could just change the help text to indicate that each input needs to be greater than zero.
import numpy as np
def P_conc_docs(t, P0, k1prime, k2, LOUD=False):
""" Computes the concentration of phenolphthalein as a function of time
Args:
t: time, either a scalar (float) or vector (numpy array)
P0: initial concentration of phenolphthalein, MUST BE POSITIVE
k1prime: forward reaction rate constant, MUST BE POSITIVE
k2: reverse reaction rate constant, MUST BE POSITIVE
Returns:
concentration, either a scalar (float) or vector (numpy array) depening in arg. t
"""
# compute first fraction in formula
frac1 = k2/(k1prime + k2)
# compute section fraction in formula
frac2 = k1prime/(k1prime + k2)
# compute term that goes into exponent. bug has been fixed
pre_expon = -(k1prime + k2)*t
# evaluate exponent
expon = np.exp(pre_expon)
# combine together
P = P0*(frac1 + frac2*expon)
if LOUD:
print("frac1 = ",frac1)
print("frac2 = ",frac2)
print("pre_expon = ",pre_expon)
print("expon = ",expon)
print("P = ",P)
# return final value
return P
help(P_conc_docs)
#the result is the same
test_P = P_conc_docs(t=10.0, P0=-10, k1prime=1, k2=1)
print("With negative P0, P(t=10) =", test_P)
Well the result is the same: if the user doesn't respect the instruction for the function inputs, the function will behave strangely. It could be even worse:
test_P = P_conc_docs(t=10.0, P0=-10, k1prime=0.0, k2=0.0)
print("With negative k1prime and negative k2, P(t=10) =", test_P)
This is where the assert statement comes in. It can assure that the arguments to the function are what they should be.
import numpy as np
def P_conc_assert(t, P0, k1prime, k2, LOUD=False):
""" Computes the concentration of phenolphthalein as a function of time
Args:
t: time, either a scalar (float) or vector (numpy array)
P0: initial concentration of phenolphthalein, MUST BE POSITIVE
k1prime: forward reaction rate constant, MUST BE POSITIVE
k2: reverse reaction rate constant, MUST BE POSITIVE
Returns:
concentration, either a scalar (float) or vector (numpy array) depening in arg. t
"""
assert (P0 >0)
assert (k1prime > 0)
assert (k2 >0)
assert (k1prime + k2 > 1E-10)
# compute first fraction in formula
frac1 = k2/(k1prime + k2)
# compute section fraction in formula
frac2 = k1prime/(k1prime + k2)
# compute term that goes into exponent. bug has been fixed
pre_expon = -(k1prime + k2)*t
# evaluate exponent
expon = np.exp(pre_expon)
# for some problems, it makes sense to also apply assert to intermediate values
# combine together
P = P0*(frac1 + frac2*expon)
if LOUD:
print("frac1 = ",frac1)
print("frac2 = ",frac2)
print("pre_expon = ",pre_expon)
print("expon = ",expon)
print("P = ",P)
# return final value
return P
test_P = P_conc_assert(t=10.0, P0=1.0, k1prime=0.0, k2=0.0)
print("With negative k1prime and negative k2, P(t=10) =", test_P)
Notice that python tells us which assertion failed so we know that the function call had a bad value of k1prime. The program still fails, but it tells us exactly why. We can also use assertions to test that the code behaves the way we expect.
These types of assert statements can help you debug later on down the road.
1.10 Error Handling¶
There are times when you want to handle an error gracefully. When a program is running, if it encounters an error it can raise and exception. The exception will give some indication of what the error is. In Python, you can place some code in a special block of code called a try block. After the try block, the types of exception to handle are listed using except blocks. Only the exceptions that are explicitly handled are caught.
One type of exception is the ZeroDivisionError that is raised when a number is divided by zero. First, we will look at an uncaught exception:
z = 10.5/0
try:
z = 10.5/0
except ZeroDivisionError:
print("You cannot Divide by zero")
One thing that happens when you catch a raised exception, is that the program will continue on. This can be a useful feature, but often times you want to catch an exception, print a useful error message, and then have the program end. This can be done by adding a raise statement to the end of the except block. The raise statement tells Python to still fail due to the exception, despite the fact that we caught it. This changes the previous example by one line, but changes the output and forces the program to quit:
try:
z = 10.5/0
except ZeroDivisionError:
print("You cannot divide by 0, exiting")
raise # kill the program
Let's see this applied to our reaction kinetics example.
import numpy as np
def P_conc_assert2(t, P0, k1prime, k2, LOUD=False):
""" Computes the concentration of phenolphthalein as a function of time
Args:
t: time, either a scalar (float) or vector (numpy array)
P0: initial concentration of phenolphthalein, MUST BE POSITIVE
k1prime: forward reaction rate constant, MUST BE POSITIVE
k2: reverse reaction rate constant, MUST BE POSITIVE
Returns:
concentration, either a scalar (float) or vector (numpy array) depening in arg. t
"""
try:
assert (P0 >0)
assert (k1prime > 0)
assert (k2 >0)
assert (k1prime + k2 > 1E-10)
except AssertionError:
print("Input Parameters are not all positive.")
print("P0 =",P0)
print("k1prime =",k1prime)
print("k2 =",k2)
raise # kill the program
except:
print("An unexpected error occurred when",
"checking the function parameters")
raise # kill the program
# compute first fraction in formula
frac1 = k2/(k1prime + k2)
# compute section fraction in formula
frac2 = k1prime/(k1prime + k2)
# compute term that goes into exponent. bug has been fixed
pre_expon = -(k1prime + k2)*t
# evaluate exponent
expon = np.exp(pre_expon)
# for some problems, it makes sense to also apply assert to intermediate values
# combine together
P = P0*(frac1 + frac2*expon)
if LOUD:
print("frac1 = ",frac1)
print("frac2 = ",frac2)
print("pre_expon = ",pre_expon)
print("expon = ",expon)
print("P = ",P)
# return final value
return P
With this function, if it is passed a negative value, will raise an AssertionError. The code catches this error, prints out the input parameters to the user, and then exits by raising the exception.
test_P = P_conc_assert2(t=10.0, P0=1.0, k1prime=0.0, k2=0.0)
print("With negative k1prime and negative k2, P(t=10) =", test_P)
Also, the function has a generic except statement that will catch any other errors in the try block. Because the try block involves comparison of numbers, if we pass a string as a parameter, there will be an error when checking if that parameter is greater than 0. This exception will be caught by the generic except statement, and the code prints out an error message, and then quits:
test_P = P_conc_assert2(t=10.0, P0="hello world", k1prime=0.0, k2=0.0)
print("With negative k1prime and negative k2, P(t=10) =", test_P)