$$ \newcommand{\tp}{\thinspace .} $$

Study guide: Software engineering with exponential decay models

Hans Petter Langtangen [1, 2]

[1] Center for Biomedical Computing, Simula Research Laboratory
[2] Department of Informatics, University of Oslo

Feb 12, 2016










Table of contents

From flat program to module with functions
      Mathematical model problem
      Many will make a rough, flat program first
      There are major issues with this solution
      New flat program
      Such flat programs are ideal for IPython notebooks!
      But: Further development of such flat programs require many scattered edits - easy to make mistakes!
      The DRY principle: Don't repeat yourself!
      Make sure any program file is a valid Python module
      The requirements of a module are so simple
      The module file decay.py for our example
      The module file decay.py for our example w/prefix
      How do we add code for comparing schemes visually?
      We just add a new function with the tailored plotting
      The comparison
      Prefixing imported functions by the module name
      Documenting functions and modules
      Example on NumPy-style doc string
      Example on an autogenerated nice HTML manual
      Logging intermediate results
      Introductory example on using a logger
      A message is tied to a level, and one can specify one many levels that get printed
      Using a logger in our solver function
      Monitoring messages
User interfaces
      Accessing command-line arguments
      Reading a sequence of command-line arguments
      Implementation
      Working with an argument parser
      A graphical user interface
      The Parampool package
      Making a compute function
      The compute function must return HTML code
      Generating the user interface
      Running the web application
      More advanced use
Tests for verifying implementations
      Doctests
      Running doctests
      Unit testing with nose
      Basic use of nose and pytest
      Example on a test function in the source code
      Example on test functions in a separate file
      Test function for solver
      Can test that potential integer division is avoided too
Packaging the software for other users
      Python convention: install your software with setup.py
      split.py for several modules in a package
      The __init__.py file can be empty
      Always develop software and write reports with Git
      More pro use with Git
Classes for problem and solution method
      Collect physical problem and parameters in class Problem
      Collect numerical parameters and methods in class Solver
      Get input from the command line; class Problem
      Get input from the command line; class Problem
      How to combine class Problem and class Solver
Performning scientific experiments
      Model problem and numerical solution method
      Plan for the experiments
      Available software
      Required new results
      Reproducible science is key!
      What actions are needed in the script?
      Run a program from a program with subprocess
      Interpreting the output from an operating system command
      Combining plot files: PNG and PDF solutions
      Making a report
      Publishing a complete project









From flat program to module with functions









Mathematical model problem

$$ \begin{align*} u'(t) &= -au(t), \quad t \in (0,T]\\ u(0) &= I \end{align*} $$

Solution by \( \theta \)-scheme: $$ \begin{equation*} u^{n+1} = \frac{1 - (1-\theta) a\Delta t}{1 + \theta a\Delta t}u^n \end{equation*} $$

\( \theta =0 \): Forward Euler, \( \theta =1 \): Backward Euler, \( \theta =1/2 \): Crank-Nicolson (midpoint method)









Many will make a rough, flat program first 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
from numpy import *
from matplotlib.pyplot import *

A = 1
a = 2
T = 4
dt = 0.2
N = int(round(T/dt))
y = zeros(N+1)
t = linspace(0, T, N+1)
theta = 1
y[0] = A
for n in range(0, N):
    y[n+1] = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)*y[n]

y_e = A*exp(-a*t) - y
error = y_e - y
E = sqrt(dt*sum(error**2))
print 'Norm of the error: %.3E' % E
plot(t, y, 'r--o')
t_e = linspace(0, T, 1001)
y_e = A*exp(-a*t_e)
plot(t_e, y_e, 'b-')
legend(['numerical, theta=%g' % theta, 'exact'])
xlabel('t')
ylabel('y')
show()









There are major issues with this solution

  1. The notation in the program does not correspond exactly to the notation in the mathematical problem: the solution is called y and corresponds to \( u \) in the mathematical description, the variable A corresponds to the mathematical parameter \( I \), N in the program is called \( N_t \) in the mathematics.

  2. There are no comments in the program.








New flat program 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
from numpy import *
from matplotlib.pyplot import *

I = 1
a = 2
T = 4
dt = 0.2
Nt = int(round(T/dt))     # no of time intervals
u = zeros(Nt+1)           # array of u[n] values
t = linspace(0, T, Nt+1)  # time mesh
theta = 1                 # Backward Euler method

u[0] = I                  # assign initial condition
for n in range(0, Nt):    # n=0,1,...,Nt-1
    u[n+1] = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)*u[n]

# Compute norm of the error
u_e = I*exp(-a*t) - u     # exact u at the mesh points
error = u_e - u
E = sqrt(dt*sum(error**2))
print 'Norm of the error: %.3E' % E

# Compare numerical (u) and exact solution (u_e) in a plot
plot(t, u, 'r--o')
t_e = linspace(0, T, 1001)       # very fine mesh for u_e
u_e = I*exp(-a*t_e)
plot(t_e, u_e, 'b-')
legend(['numerical, theta=%g' % theta, 'exact'])
xlabel('t')
ylabel('u')
show()









Such flat programs are ideal for IPython notebooks!









But: Further development of such flat programs require many scattered edits - easy to make mistakes!

The solution formula for \( u^{n+1} \) is completely general and should be available as a Python function with all input data as function arguments and all output data returned to the calling code 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
def solver(I, a, T, dt, theta):
    """Solve u'=-a*u, u(0)=I, for t in (0,T] with steps of dt."""
    dt = float(dt)               # avoid integer division
    Nt = int(round(T/dt))        # no of time intervals
    T = Nt*dt                    # adjust T to fit time step dt
    u = np.zeros(Nt+1)           # array of u[n] values
    t = np.linspace(0, T, Nt+1)  # time mesh

    u[0] = I                  # assign initial condition
    for n in range(0, Nt):    # n=0,1,...,Nt-1
        u[n+1] = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)*u[n]
    return u, t

Call: 

1
u, t = solver(I=1, a=2, T=4, dt=0.2, theta=0.5)









The DRY principle: Don't repeat yourself!

DRY:

When implementing a particular functionality in a computer program, make sure this functionality and its variations are implemented in just one piece of code. That is, if you need to revise the implementation, there should be one and only one place to edit. It follows that you should never duplicate code (don't repeat yourself!), and code snippets that are similar should be factored into one piece (function) and parameterized (by function arguments).









Make sure any program file is a valid Python module 

1
2
3
from decay import solver
# Solve a decay problem
u, t = solver(I=1, a=2, T=4, dt=0.2, theta=0.5)

or prefix function names by the module name: 

1
2
3
import decay
# Solve a decay problem
u, t = decay.solver(I=1, a=2, T=4, dt=0.2, theta=0.5)









The requirements of a module are so simple

  1. The filename without .py must be a valid Python variable name.

  2. The main program must be executed (through statements or a function call) in the test block.
The test block is normally placed at the end of a module file: 

1
2
if __name__ == '__main__':
    # Statements

If the file is imported, the if test fails and no main program is run, otherwise, the file works as a program









The module file decay.py for our example 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
from numpy import *
from matplotlib.pyplot import *

def solver(I, a, T, dt, theta):
    ...

def u_exact(t, I, a):
    return I*exp(-a*t)

def experiment_compare_numerical_and_exact():
    I = 1;  a = 2;  T = 4;  dt = 0.4;  theta = 1
    u, t = solver(I, a, T, dt, theta)

    t_e = linspace(0, T, 1001)       # very fine mesh for u_e
    u_e = u_exact(t_e, I, a)

    plot(t,   u,   'r--o')           # dashed red line with circles
    plot(t_e, u_e, 'b-')             # blue line for u_e
    legend(['numerical, theta=%g' % theta, 'exact'])
    xlabel('t')
    ylabel('u')
    plotfile = 'tmp'
    savefig(plotfile + '.png');  savefig(plotfile + '.pdf')

    error = u_exact(t, I, a) - u
    E = sqrt(dt*sum(error**2))
    print 'Error norm:', E

if __name__ == '__main__':
    experiment_compare_numerical_and_exact()

Complete file: decay.py









The module file decay.py for our example w/prefix 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
import numpy as np
import matplotlib.pyplot as plt

def solver(I, a, T, dt, theta):
    ...

def u_exact(t, I, a):
    return I*np.exp(-a*t)

def experiment_compare_numerical_and_exact():
    I = 1;  a = 2;  T = 4;  dt = 0.4;  theta = 1
    u, t = solver(I, a, T, dt, theta)

    t_e = np.linspace(0, T, 1001)       # very fine mesh for u_e
    u_e = u_exact(t_e, I, a)

    plt.plot(t,   u,   'r--o')       # dashed red line with circles
    plt.plot(t_e, u_e, 'b-')         # blue line for u_e
    plt.legend(['numerical, theta=%g' % theta, 'exact'])
    plt.xlabel('t')
    plt.ylabel('u')
    plotfile = 'tmp'
    plt.savefig(plotfile + '.png');  plt.savefig(plotfile + '.pdf')

    error = u_exact(t, I, a) - u
    E = np.sqrt(dt*np.sum(error**2))
    print 'Error norm:', E

if __name__ == '__main__':
    experiment_compare_numerical_and_exact()









How do we add code for comparing schemes visually?

Think of edits in the flat program that are required to produce this plot (!)









We just add a new function with the tailored plotting 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
def experiment_compare_schemes():
    """Compare theta=0,1,0.5 in the same plot."""
    I = 1;  a = 2;  T = 4;  dt = 0.4
    legends = []
    for theta in [0, 1, 0.5]:
        u, t = solver(I, a, T, dt, theta)
        plt.plot(t, u, '--o')
        legends.append('theta=%g' % theta)
    t_e = np.linspace(0, T, 1001)        # very fine mesh for u_e
    u_e = u_exact(t_e, I, a)
    plt.plot(t_e, u_e, 'b-')
    legends.append('exact')
    plt.legend(legends, loc='upper right')
    plotfile = 'tmp'
    plt.savefig(plotfile + '.png');  plt.savefig(plotfile + '.pdf')

import logging
# Define a default logger that does nothing
logging.getLogger('decay').addHandler(logging.NullHandler())

def solver_with_logging(I, a, T, dt, theta):
    """Solve u'=-a*u, u(0)=I, for t in (0,T] with steps of dt."""
    dt = float(dt)               # avoid integer division
    Nt = int(round(T/dt))        # no of time intervals
    T = Nt*dt                    # adjust T to fit time step dt
    u = np.zeros(Nt+1)           # array of u[n] values
    t = np.linspace(0, T, Nt+1)  # time mesh
    logging.debug('solver: dt=%g, Nt=%g, T=%g' % (dt, Nt, T))

    u[0] = I                  # assign initial condition
    for n in range(0, Nt):    # n=0,1,...,Nt-1
        u[n+1] = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)*u[n]

        logging.info('u[%d]=%g' % (n, u[n]))
        logging.debug('1 - (1-theta)*a*dt: %g, %s' %
                      (1-(1-theta)*a*dt,
                       str(type(1-(1-theta)*a*dt))[7:-2]))
        logging.debug('1 + theta*dt*a: %g, %s' %
                      (1 + theta*dt*a,
                       str(type(1 + theta*dt*a))[7:-2]))
    return u, t

def configure_basic_logger():
    logging.basicConfig(
        filename='decay.log', filemode='w', level=logging.DEBUG,
        format='%(asctime)s - %(levelname)s - %(message)s',
        datefmt='%Y.%m.%d %I:%M:%S %p')









The comparison









Prefixing imported functions by the module name

MATLAB-style names (linspace, plot): 

1
2
from numpy import *
from matplotlib.pyplot import *

Python community convention is to prefix with module name (np.linspace, plt.plot): 

1
2
import numpy as np
import matplotlib.pyplot as plt









Documenting functions and modules









Example on NumPy-style doc string 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
def solver(I, a, T, dt, theta):
    """
    Solve :math:`u'=-au` with :math:`u(0)=I` for :math:`t \in (0,T]`
    with steps of `dt` and the method implied by `theta`.

    Parameters
    ----------
    I: float
        Initial condition.
    a: float
        Parameter in the differential equation.
    T: float
        Total simulation time.
    theta: float, int
        Parameter in the numerical scheme. 0 gives
        Forward Euler, 1 Backward Euler, and 0.5
        the centered Crank-Nicolson scheme.

    Returns
    -------
    `u`: array
        Solution array.
    `t`: array
        Array with time points corresponding to `u`.

    Examples
    --------
    Solve :math:`u' = -\\frac{1}{2}u, u(0)=1.5`
    with the Crank-Nicolson method:

    >>> u, t = solver(I=1.5, a=0.5, T=9, theta=0.5)
    >>> import matplotlib.pyplot as plt
    >>> plt.plot(t, u)
    >>> plt.show()
    """









Example on an autogenerated nice HTML manual









Logging intermediate results

Most programming languages has a logging object for this purpose:







Introductory example on using a logger 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
import logging
import logging
logging.basicConfig(
    filename='myprog.log', filemode='w', level=logging.WARNING,
    format='%(asctime)s - %(levelname)s - %(message)s',
    datefmt='%m/%d/%Y %I:%M:%S %p')
logging.info('Here is some general info.')
logging.warning('Here is a warning.')
logging.debug('Here is some debugging info.')
logging.critical('Dividing by zero!')
logging.error('Encountered an error.')

Output in myprog.log: 

1
2
3
4
09/26/2015 09:25:10 AM - INFO - Here is some general info.
09/26/2015 09:25:10 AM - WARNING - Here is a warning.
09/26/2015 09:25:10 AM - CRITICAL - Dividing by zero!
09/26/2015 09:25:10 AM - ERROR - Encountered an error.









A message is tied to a level, and one can specify one many levels that get printed

Levels: critical, error, warning, info, debug

  1. level=logging.CRITICAL: print critical messages

  2. level=logging.ERROR: print critical and error messages

  3. level=logging.WARNING: print critical, error, and warning messages

  4. level=logging.INFO: print critical, error, warning, and info messages

  5. level=logging.DEBUG: print critical, error, warning, info, and debug messages








Using a logger in our solver function 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
import logging
# Define a default logger that does nothing
logging.getLogger('decay').addHandler(logging.NullHandler())

def solver_with_logging(I, a, T, dt, theta):
    """Solve u'=-a*u, u(0)=I, for t in (0,T] with steps of dt."""
    dt = float(dt)               # avoid integer division
    Nt = int(round(T/dt))        # no of time intervals
    T = Nt*dt                    # adjust T to fit time step dt
    u = np.zeros(Nt+1)           # array of u[n] values
    t = np.linspace(0, T, Nt+1)  # time mesh
    logging.debug('solver: dt=%g, Nt=%g, T=%g' % (dt, Nt, T))

    u[0] = I                  # assign initial condition
    for n in range(0, Nt):    # n=0,1,...,Nt-1
        u[n+1] = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)*u[n]

        logging.info('u[%d]=%g' % (n, u[n]))
        logging.debug('1 - (1-theta)*a*dt: %g, %s' %
                      (1-(1-theta)*a*dt,
                       str(type(1-(1-theta)*a*dt))[7:-2]))
        logging.debug('1 + theta*dt*a: %g, %s' %
                      (1 + theta*dt*a,
                       str(type(1 + theta*dt*a))[7:-2]))
    return u, t

def configure_basic_logger():
    logging.basicConfig(
        filename='decay.log', filemode='w', level=logging.DEBUG,
        format='%(asctime)s - %(levelname)s - %(message)s',
        datefmt='%Y.%m.%d %I:%M:%S %p')









Monitoring messages

One terminal window (1M steps!): 

1
2
3
>>> import decay
>>> u, t = decay.solver_with_logging(I=1, a=0.5, T=10, \ 
           dt=0.5, theta=0.5)

Another terminal window: 

1
2
3
4
5
6
7
Terminal> tail -f decay.log
2015.09.26 05:37:41 AM - INFO - u[0]=1
2015.09.26 05:37:41 AM - INFO - u[1]=0.777778
2015.09.26 05:37:41 AM - INFO - u[2]=0.604938
2015.09.26 05:37:41 AM - INFO - u[3]=0.470508
2015.09.26 05:37:41 AM - INFO - u[4]=0.36595
2015.09.26 05:37:41 AM - INFO - u[5]=0.284628

Or if level=logging.DEBUG: 

1
2
3
4
5
6
7
8
Terminal> tail -f decay.log
2015.09.26 05:40:01 AM - DEBUG - solver: dt=0.5, Nt=20, T=10
2015.09.26 05:40:01 AM - INFO - u[0]=1
2015.09.26 05:40:01 AM - DEBUG - 1 - (1-theta)*a*dt: 0.875, float
2015.09.26 05:40:01 AM - DEBUG - 1 + theta*dt*a: 1.125, float
2015.09.26 05:40:01 AM - INFO - u[1]=0.777778
2015.09.26 05:40:01 AM - DEBUG - 1 - (1-theta)*a*dt: 0.875, float
2015.09.26 05:40:01 AM - DEBUG - 1 + theta*dt*a: 1.125, float









User interfaces









Accessing command-line arguments 

1
2
Terminal> python myprog.py 1.5 2 0.5 0.8 0.4
Terminal> python myprog.py --I 1.5 --a 2 --dt 0.8 0.4









Reading a sequence of command-line arguments

Required input:

Give these on the command line in correct sequence 

1
Terminal> python decay_cml.py 1.5 0.5 4 CN 0.1 0.2 0.05









Implementation 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
def define_command_line_options():
    import argparse
    parser = argparse.ArgumentParser()
    parser.add_argument(
        '--I', '--initial_condition', type=float,
        default=1.0, help='initial condition, u(0)',
        metavar='I')
    parser.add_argument(
        '--a', type=float, default=1.0,
        help='coefficient in ODE', metavar='a')
    parser.add_argument(
        '--T', '--stop_time', type=float,
        default=1.0, help='end time of simulation',
        metavar='T')
    parser.add_argument(
        '--scheme', type=str, default='CN',
        help='FE, BE, or CN')
    parser.add_argument(
        '--dt', '--time_step_values', type=float,
        default=[1.0], help='time step values',
        metavar='dt', nargs='+', dest='dt_values')
    return parser

Note:









Working with an argument parser

Set option-value pairs on the command line if the default value is not suitable: 

1
Terminal> python decay_argparse.py --I 1.5 --a 2 --dt 0.8 0.4

Code: 

1
2
3
4
5
6
7
def read_command_line_argparse():
    parser = define_command_line_options()
    args = parser.parse_args()
    scheme2theta = {'BE': 1, 'CN': 0.5, 'FE': 0}
    data = (args.I, args.a, args.T, scheme2theta[args.scheme],
            args.dt_values)
    return data

(metavar is the symbol used in help output)









A graphical user interface

Normally very much programming required - and much competence on graphical user interfaces.

Here: use a tool to automatically create it in a few minutes (!)









The Parampool package

Remark.

The forthcoming material aims at those with particular interest in equipping their programs with a GUI - others can safely skip it.









Making a compute function









The compute function must return HTML code 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
def main_GUI(I=1.0, a=.2, T=4.0,
             dt_values=[1.25, 0.75, 0.5, 0.1],
             theta_values=[0, 0.5, 1]):
    # Build HTML code for web page. Arrange plots in columns
    # corresponding to the theta values, with dt down the rows
    theta2name = {0: 'FE', 1: 'BE', 0.5: 'CN'}
    html_text = '<table>\n'
    for dt in dt_values:
        html_text += '<tr>\n'
        for theta in theta_values:
            E, html = compute4web(I, a, T, dt, theta)
            html_text += """
<td>
<center><b>%s, dt=%g, error: %.3E</b></center><br>
%s
</td>
""" % (theta2name[theta], dt, E, html)
        html_text += '</tr>\n'
    html_text += '</table>\n'
    return html_text









Generating the user interface

Make a file decay_GUI_generate.py: 

1
2
3
4
5
6
from parampool.generator.flask import generate
from decay import main_GUI
generate(main_GUI,
         filename_controller='decay_GUI_controller.py',
         filename_template='decay_GUI_view.py',
         filename_model='decay_GUI_model.py')

Running decay_GUI_generate.py results in

  1. decay_GUI_model.py defines HTML widgets to be used to set input data in the web interface,

  2. templates/decay_GUI_views.py defines the layout of the web page,

  3. decay_GUI_controller.py runs the web application.
Good news: we only need to run decay_GUI_controller.py and there is no need to look into any of these files!









Running the web application

Start the GUI 

1
Terminal> python decay_GUI_controller.py

Open a web browser at 127.0.0.1:5000









More advanced use









Tests for verifying implementations









Doctests

Doc strings can be equipped with interactive Python sessions for demonstrating usage and automatic testing of functions. 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
def solver(I, a, T, dt, theta):
    """
    Solve u'=-a*u, u(0)=I, for t in (0,T] with steps of dt.


    >>> u, t = solver(I=0.8, a=1.2, T=2, dt=0.5, theta=0.5)
    >>> for t_n, u_n in zip(t, u):
    ...     print 't=%.1f, u=%.14f' % (t_n, u_n)
    t=0.0, u=0.80000000000000
    t=0.5, u=0.43076923076923
    t=1.0, u=0.23195266272189
    t=1.5, u=0.12489758761948
    t=2.0, u=0.06725254717972
    """
    ...









Running doctests

Automatic check that the code reproduces the doctest output: 

1
Terminal> python -m doctest decay.py

Floats are difficult to compare.

Limit the number of digits in the output in doctests! Otherwise, round-off errors on a different machine may ruin the test.









Unit testing with nose









Basic use of nose and pytest

  1. Implement tests in test functions with names starting with test_.

  2. Test functions cannot have arguments.

  3. Test functions perform assertions on computed results using assert functions from the nose.tools module.

  4. Test functions can be in the source code files or be collected in separate files test*.py.








Example on a test function in the source code

Very simple module mymod (in file mymod.py): 

1
2
def double(n):
    return 2*n

Write test function in mymod.py: 

1
2
3
4
5
6
7
8
def double(n):
    return 2*n

def test_double():
    n = 4
    expected = 2*4
    computed = double(n)
    assert expected == computed

Running one of 

1
2
Terminal> nosetests -s -v mymod
Terminal> py.test   -s -v mymod

makes the framework run all test_*() functions in mymod.py.









Example on test functions in a separate file

Write the test in a separate file, say test_mymod.py: 

1
2
3
4
5
6
7
import mymod

def test_double():
    n = 4
    expected = 2*4
    computed = double(n)
    assert expected == computed

Running one of 

1
2
Terminal> nosetests -s -v
Terminal> py.test   -s -v

makes the frameworks run all test_*() functions in all files test*.py in the current directory and in all subdirectories (pytest) or just those with names tests or *_tests (nose)

Tip.

Start with test functions in the source code file. When the file contains many tests, or when you have many source code files, move tests to separate files.









Test function for solver

Use exact discrete solution of the \( \theta \) scheme as test: $$ u^n = I\left( \frac{1 - (1-\theta) a\Delta t}{1 + \theta a \Delta t} \right)^n$$ 

1
2
3
4
5
def u_discrete_exact(n, I, a, theta, dt):
    """Return exact discrete solution of the numerical schemes."""
    dt = float(dt)  # avoid integer division
    A = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)
    return I*A**n


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
def test_u_discrete_exact():
    """Check that solver reproduces the exact discr. sol."""
    theta = 0.8; a = 2; I = 0.1; dt = 0.8
    Nt = int(8/dt)  # no of steps
    u, t = solver(I=I, a=a, T=Nt*dt, dt=dt, theta=theta)

    # Evaluate exact discrete solution on the mesh
    u_de = np.array([u_discrete_exact(n, I, a, theta, dt)
                     for n in range(Nt+1)])

    # Find largest deviation
    diff = np.abs(u_de - u).max()
    tol = 1E-14
    success = diff < tol
    assert success









Can test that potential integer division is avoided too

Warning.

If \( a \), \( \Delta t \), and \( \theta \) are integers, the formula for \( u^{n+1} \) in the solver function may lead to 0 because of unintended integer division. 

1
2
3
4
5
6
7
8
9
def test_potential_integer_division():
    """Choose variables that can trigger integer division."""
    theta = 1; a = 1; I = 1; dt = 2
    Nt = 4
    u, t = solver(I=I, a=a, T=Nt*dt, dt=dt, theta=theta)
    u_de = np.array([u_discrete_exact(n, I, a, theta, dt)
                     for n in range(Nt+1)])
    diff = np.abs(u_de - u).max()
    assert diff < 1E-14









Packaging the software for other users

Python convention: install your software with setup.py

Installation of a single module file decay.py: 

1
2
3
4
5
6
from distutils.core import setup
setup(name='decay',
      version='0.1',
      py_modules=['decay'],
      scripts=['decay.py'],
      )

Installation: 

1
Terminal> sudo python setup.py install

(Many variants!)









split.py for several modules in a package

setup.py: 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
from distutils.core import setup
import os

setup(name='decay',
      version='0.1',
      author='Hans Petter Langtangen',
      author_email='hpl@simula.no',
      url='https://github.com/hplgit/decay-package/',
      packages=['decay'],
      scripts=[os.path.join('decay', 'decay.py')]
     )









The __init__.py file can be empty

Empty __init__.py: 

1
2
import decay
u, t = decay.decay.solver(...)

Do this in __init__.py to avoid decay.decay.solver: 

1
from decay import *

Can now write 

1
2
3
4
5
6
import decay
u, t = decay.solver(...)

# or
from decay import solver
u, t = solver(...)









Always develop software and write reports with Git

The Git work cycle: 

1
2
3
4
5
git pull                # before starting a new session
# edit files
git add mynewfile       # remember to add new files!
git commit -am 'Short description of what I did'
git push origin master  # before end of day or a break









More pro use with Git

See what others have done in the project: 

1
2
3
git fetch origin         # instead of git pull
git diff origin/master   # what are the changes?
git merge origin/master  # update my files

Develop new features in a separate branch: 

1
2
3
4
5
git branch newstuff
git checkout newstuff
# edit files
git commit -am 'Changed ...'
git push origin newstuff

When newstuff is tested and matured, merge back in master: 

1
2
git checkout master
git merge newstuff

















Classes for problem and solution method









Collect physical problem and parameters in class Problem 

1
2
3
4
5
6
7
8
9
from numpy import exp

class Problem(object):
    def __init__(self, I=1, a=1, T=10):
        self.T, self.I, self.a = I, float(a), T

    def u_exact(self, t):
        I, a = self.I, self.a
        return I*exp(-a*t)









Collect numerical parameters and methods in class Solver 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
class Solver(object):
    def __init__(self, problem, dt=0.1, theta=0.5):
        self.problem = problem
        self.dt, self.theta = float(dt), theta

    def solve(self):
        self.u, self.t = solver(
            self.problem.I, self.problem.a, self.problem.T,
            self.dt, self.theta)

    def error(self):
        """Return norm of error at the mesh points."""
        u_e = self.problem.u_exact(self.t)
        e = u_e - self.u
        E = np.sqrt(self.dt*np.sum(e**2))
        return E









Get input from the command line; class Problem 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
class Problem(object):
    def __init__(self, I=1, a=1, T=10):
        self.T, self.I, self.a = I, float(a), T

    def define_command_line_options(self, parser=None):
        """Return updated (parser) or new ArgumentParser object."""
        if parser is None:
            import argparse
            parser = argparse.ArgumentParser()

        parser.add_argument(
            '--I', '--initial_condition', type=float,
            default=1.0, help='initial condition, u(0)',
            metavar='I')
        parser.add_argument(
            '--a', type=float, default=1.0,
            help='coefficient in ODE', metavar='a')
        parser.add_argument(
            '--T', '--stop_time', type=float,
            default=1.0, help='end time of simulation',
            metavar='T')
        return parser

    def init_from_command_line(self, args):
        """Load attributes from ArgumentParser into instance."""
        self.I, self.a, self.T = args.I, args.a, args.T









Get input from the command line; class Problem 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
class Solver(object):
    def __init__(self, problem, dt=0.1, theta=0.5):
        self.problem = problem
        self.dt, self.theta = float(dt), theta

    def define_command_line_options(self, parser):
        """Return updated (parser) or new ArgumentParser object."""
        parser.add_argument(
            '--scheme', type=str, default='CN',
            help='FE, BE, or CN')
        parser.add_argument(
            '--dt', '--time_step_values', type=float,
            default=[1.0], help='time step values',
            metavar='dt', nargs='+', dest='dt_values')
        return parser

    def init_from_command_line(self, args):
        """Load attributes from ArgumentParser into instance."""
        self.dt, self.theta = args.dt, args.theta









How to combine class Problem and class Solver 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
def experiment_classes():
    problem = Problem()
    solver = Solver(problem)

    # Read input from the command line
    parser = problem.define_command_line_options()
    parser = solver. define_command_line_options(parser)
    args = parser.parse_args()
    problem.init_from_command_line(args)
    solver. init_from_command_line(args)

    # Solve and plot
    solver.solve()
    import matplotlib.pyplot as plt
    t_e = np.linspace(0, T, 1001)    # very fine mesh for u_e
    u_e = problem.u_exact(t_e)

    plt.plot(t,   u,   'r--o')       # dashed red line with circles
    plt.plot(t_e, u_e, 'b-')         # blue line for u_e
    plt.legend(['numerical, theta=%g' % theta, 'exact'])
    plt.xlabel('t')
    plt.ylabel('u')
    plt.show()









Performning scientific experiments

Goals:

  1. Explore the behavior of a numerical method for an ODE

  2. Show how a program can set up, execute, and report scientific investigations

  3. Demonstrate how to write a scientific report

  4. Demonstrate various technologies for reports: HTML w/MathJax, LaTeX, Sphinx, IPython notebooks, ...








Model problem and numerical solution method

Problem: $$ \begin{equation} u'(t) = -au(t),\quad u(0)=I,\ 0 < t \leq T, \label{decay:experiments:model} \end{equation} $$

Solution method (\( \theta \)-rule): $$ u^{n+1} = \frac{1 - (1-\theta) a\Delta t}{1 + \theta a\Delta t}u^n, \quad u^0=I\tp $$









Plan for the experiments

For fixed \( I \), \( a \), and \( T \), we run the three schemes for various values of \( \Delta t \), and present in a report the following results:

  1. visual comparison of the numerical and exact solution in a plot for each \( \Delta t \) and \( \theta=0,1,\frac{1}{2} \),

  2. a table and a plot of the norm of the numerical error versus \( \Delta t \) for \( \theta=0,1,\frac{1}{2} \).








Available software

model.py: 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
Terminal> python model.py --I 1.5 --a 0.25 --T 6 --dt 1.25 0.75 0.5
0.0   1.25:    5.998E-01
0.0   0.75:    1.926E-01
0.0   0.50:    1.123E-01
0.0   0.10:    1.558E-02
0.5   1.25:    6.231E-02
0.5   0.75:    1.543E-02
0.5   0.50:    7.237E-03
0.5   0.10:    2.469E-04
1.0   1.25:    1.766E-01
1.0   0.75:    8.579E-02
1.0   0.50:    6.884E-02
1.0   0.10:    1.411E-02

+ a set of plot files of numerial vs exact solution









Required new results

Must write a script exper1.py to automate running model.py and generating these results 

1
Terminal> python exper1.py 0.5 0.25 0.1 0.05

(\( \Delta t \) values on the comand line)









Reproducible science is key!

Let your scientific investigations be automated by scripts!









What actions are needed in the script?

  1. Run model.py program with appropriate input

  2. Interpret the output and make table and plot of numerical errors

  3. Combine plot files to new figures
Complete script: exper1.py









Run a program from a program with subprocess

Command to be run: 

1
python model.py --I 1.2 --a 0.2 --T 8 -dt 1.25 0.75 0.5 0.1

Constructed in Python: 

1
2
3
4
# Given I, a, T, and a list dt_values
cmd = 'python model.py --I %g --a %g --T %g' % (I, a, T)
dt_values_str = ' '.join([str(v) for v in dt_values])
cmd += ' --dt %s' % dt_values_str

Run under the operating system: 

1
2
3
4
5
6
7
from subprocess import Popen, PIPE, STDOUT
p = Popen(cmd, shell=True, stdout=PIPE, stderr=STDOUT)
output, dummy = p.communicate()

failure = p.returncode
if failure:
    print 'Command failed:', cmd; sys.exit(1)









Interpreting the output from an operating system command

The output if the previous command run by subprocess is in a string output: 

1
2
3
4
5
6
7
8
errors = {'dt': dt_values, 1: [], 0: [], 0.5: []}
for line in output.splitlines():
    words = line.split()
    if words[0] in ('0.0', '0.5', '1.0'):  # line with E?
        # typical line: 0.0   1.25:    7.463E+00
        theta = float(words[0])
        E = float(words[2])
        errors[theta].append(E)









Combining plot files: PNG and PDF solutions

PNG: 

1
2
3
4
Terminal> montage -background white -geometry 100% -tile 2x \ 
          f1.png f2.png f3.png f4.png f.png
Terminal> convert -trim f.png f.png
Terminal> convert f.png -transparent white f.png

PDF: 

1
2
3
4
Terminal> pdftk f1.pdf f2.pdf f3.pdf f4.pdf output tmp.pdf
Terminal> pdfnup --nup 2x2 --outfile tmp.pdf tmp.pdf
Terminal> pdfcrop tmp.pdf f.pdf
Terminal> rm -f tmp.pdf

Easy to build these commands in Python and execute them with subprocess or os.system: os.system(cmd)









Making a report









Publishing a complete project