from IPython.core.display import HTML

def _set_css_style(css_file_path):
   """
   Read the custom CSS file and load it into Jupyter.
   Pass the file path to the CSS file.
   """
   styles = open(css_file_path, "r").read()
   s = '<style>%s</style>' % styles
   return HTML(s)

_set_css_style('rise.css')

Principles of code design and functions

print view
notebook

• Best practices for writing functions
• Code documentation
• assert and other checks
• Errors and exceptions

Modeling biochemical reactions with Euler's method

For an ODE $\,\dfrac{d\mathbf{x}}{dt} = \mathbf{f}(\mathbf{x}, t)\,$ and small step $\Delta t$:

$$ \mathbf{x}(t+\Delta t) \;\approx\; \mathbf{x}(t) + \Delta t \cdot \mathbf{f}(\mathbf{x}(t), t)\,. $$

Last time, we used this approach to model $A \xrightarrow{k} \varnothing$, mathematically,

$$ \frac{dA}{dt} = -k\,A, \qquad A(0)=A_0.$$

Today, let's write a function to do this

Scaffolding

def euler_decay_step():
    pass

• What arguments should we include?
• What should we call them?
• What should the function do? Return a number, operate on input?

General principles

• Where possible, prefer pure functions with clear input/output and without global variables
• Clear names for functions and variables
• Clear defaults where appropriate
• Sufficient comments and documentation
• Make it easy to test (more on this soon)

Laying out the function in words

def euler_decay_step(A, k_per_s, dt_s):  # Units included for parameters

    # Define the time derivative of concentration
    def dAdt(A, k):
        pass

    # Get the change in concentration in the time step dt_s

    # Return the new value of the concentration

def euler_decay_step(A, k_per_s, dt_s):  # Units included for parameters

    # Define the time derivative of concentration
    def dAdt(A, k):
        return -k * A

    # Get the change in concentration in the time step dt_s
    dA = dAdt(A, k_per_s) * dt_s
        
    # Return the new value of the concentration
    return A + dA

Testing our function

import numpy as np

# Parameters
A0 = 1.0        # initial concentration [arbitrary units]
k_per_s = 0.8   # rate constant [1/s]
t_end_s = 6.0   # total simulation time [s]
dt_s = 0.1      # simulation time step size [s]

t_s = np.arange(0, t_end_s, dt_s)
A_euler = [A0]

# Iterate through time, updating the last concentration and adding it to the list
for t in t_s[1:]:
    A_euler.append(euler_decay_step(A_euler[-1], k_per_s, dt_s))

import matplotlib.pyplot as plt

plt.plot(t_s, A_euler)
plt.xlabel('time [s]')
plt.ylabel('A [a.u.]');

Documentation with docstrings

def euler_decay_step(A, k_per_s, dt_s):  # Units included for parameters
    '''
    Integrates a differential equation for irreversible decay (A -> 0) with Euler's method.

    Implements the update for the model:
        dA/dt = -k_per_s * A
    
    Args:
        A (float):       concentration, units of [µM]
        k_per_s (float): decay rate, units of [1/s]
        dt_s (float):    time step, units of [s]
        
    Returns:
        A + dA (float): new concentration after one time step
    '''
    
    # Define the time derivative of concentration
    def dAdt(A, k):
        return -k * A

    # Get the change in concentration in the time step dt_s
    dA = dAdt(A, k_per_s) * dt_s
        
    # Return the new value of the concentration
    return A + dA

help(euler_decay_step)

Checking for and catching errors

• It's helpful for code to fail loudly when something goes wrong and guide the user to fix the problem.
• Good error messages save hours of guessing—include what went wrong and with which value.

1 / 0 

When something fails, read from the bottom line up:

• Exception type: ZeroDivisionError
• Short message: division by zero
• File and line number that failed

Handle an error with try / except

try:
    # code that might fail
except SomeError as e:
    # what to do if that error happens

Let’s catch the divide-by-zero and print a useful message instead of crashing.

try:
    x = 1 / 0
    print("result:", x)
except ZeroDivisionError as e:
    print("Oh no! I tried to divide by zero. Details:", e)

Catch “any” error (carefully)

You can catch any exception:

except Exception as e:
    ...

Wise to use this sparingly. Prefer specific exceptions when possible so real bugs don’t get hidden.

def risky(divisor):
    return 10 / divisor

for d in [5, 0, 'zero']:
    try:
        print(f"10 / {d} =", risky(d))
    except ZeroDivisionError:
        print('Cannot divide by zero. Please choose a nonzero divisor.')
    except TypeError:
        print('Divisor must be a number (got a different type).')

Raise your own errors

When inputs are invalid, stop early and explain the problem. Two common choices:

• ValueError — wrong value (e.g., negative number when positive expected).
• TypeError — wrong type (e.g., list instead of float).

def validate_string(test_string):
    if not (isinstance(test_string, str)):
        raise TypeError(f'test_string must be a string (got type {type(test_string).__name__}).')

for s in ['zero', 0]:
    try:
        validate_string(s)
        print(f'passed: s={s} is OK.')
    except Exception as e:
        print(f'failed: s={s} ->', type(e).__name__ + ":", e)

try / except / else / finally

• try: code that might fail
• except: how to respond to a specific failure
• else: runs only when no error happened
• finally: runs always (cleanup, closing files, etc.)

def safe_divide(a, b):
    try:
        out = a / b
    except ZeroDivisionError:
        return float("inf")  # or raise a helpful error
    else:
        return out
    finally:
        pass  # nothing to clean up here

print("safe_divide(6, 3) ->", safe_divide(6, 3))
print("safe_divide(6, 0) ->", safe_divide(6, 0))

assert shorthand

Quick way of checking conditions, but not guaranteed to execute.

assert expression

is equivalent to

if __debug__:
    if not expression: raise AssertionError

__debug__ is True by default, but can sometimes be set to False (e.g., when using code optimization)

x = 0
assert x>0, f'Oh no! Got x={x}'

How to use it

Check to make sure the input is formatted the way that you expect and that output has the right properties.

A few checks/tests can catch errors early and make them easier to fix.

Activity

Part 1. Team up in small groups and write a function that will perform one Euler integration step for an arbitrary ODE.
Input:

• x: the variable that is being integrated
• dxdt: a function that computes the derivative of x with respect to time
• theta: a set of parameters that can be passed to dxdt to help compute the derivative (e.g., k_per_s in the decay equation)
• dt: the time step

Output:

• x + dx: the value of the variable after one Euler step

Make sure to document your function and catch potential errors. Save your function in a python file euler.py.

Part 2. Swap files with a neighboring team (Slack DMs are useful for this!). Place the neighboring team's Python file in your current directory and import their integration function with

from euler import <their function name>

Test the function! Can you understand the documentation? Does it catch errors if something is wrong with the input?

Part 3. Discuss what you found with your neighbors.