Biological Computing in Python

Biological Computing in Python#

Learning goals

By the end of this chapter you should be able to:

Describe a simple algorithm in plain language and check it against a worked example
Write and execute Python scripts and understand basic syntax
Work with Python data structures (lists, tuples, dictionaries, sets)
Implement control flow (conditionals, loops) in Python
Write, document, and test functions
Handle files and directories in Python
Use list comprehensions and understand functional programming basics
Debug Python code and handle errors/exceptions
Use regular expressions for pattern matching
Work with key Python packages for scientific computing
Understand Python’s object-oriented features (see OOP chapter)

Prerequisites

Before starting this chapter, you should:

Have experience using UNIX/Linux command line (covered in UNIX and Linux)
Understand basic file operations
Have Python installed

No prior Python knowledge is assumed. Variables, functions and loops are introduced in this chapter. The shell-scripting chapter introduces algorithms, but the introduction below also provides a self-contained example.

Use your existing bootcamp project, or choose a learning project if you are studying Python independently. Follow Organising your learning project: Python scripts can share code/ with Shell and R scripts, and examples normally run from the project root. You do not need another language-specific project or prior Git study to begin; explicitly temporary exercises remain separate.

Introduction#

In the previous chapter, you used shell commands and scripts to manipulate files, combine programs and automate workflows. In this chapter, we build on those skills using Python.

Shell scripting is particularly effective for connecting existing programs and managing files. Python complements it by providing richer tools for representing data, expressing more complex logic, organizing reusable code and developing complete analyses or simulations.

From a task to an algorithm#

In shell scripting, we described a task as a sequence of precise steps before translating it into commands. That sequence is an algorithm. The same principle applies in Python: first determine what the program must do, then express those steps in Python syntax.

If you are starting here, consider four plant heights: 8, 12, 10 and 15 cm. How many exceed 10 cm? Work out the answer, then describe the steps that someone else could follow on a thousand measurements:

Start a counter at zero.
Examine each height in turn.
If the height is greater than 10 cm, increase the counter by one.
After the last height, report the counter.

These steps form an algorithm: a precise, ordered procedure for solving a specified problem. Its input is a collection of recorded heights in centimetres; its output is the number greater than 10 cm. Here, the expected answer is 2. The height of exactly 10 cm is not counted because the question says “exceed”, not “equal or exceed”.

The algorithm does not depend on a particular programming language. A program implements its steps as instructions a computer can run. In a shell script, we might implement the procedure by combining commands; in Python, we can express the data, repetition and decision directly using a collection, a loop and a conditional. Why choose Python for this work?

Why Python?#

Python complements the tools introduced in the Shell scripting chapter. Its readable syntax makes algorithms relatively easy to recognize in the resulting code, while its scientific packages support data analysis, visualization, statistics and numerical modelling. Python is therefore useful both for small standalone programs and for components of larger shell-based workflows.

Python also offers a useful balance between development time and computational performance, particularly when numerical work is handled by optimized scientific packages (see the figure below).

Language speeds — Fig. 11 **Python’s numerical computing performance compared to some others**. Smaller numbers are better. Note that the y-axis is in $\log_{10}$ scale.
(Source: http://julialang.org/)#

What you will learn#

Turning plans into reliable programs is the focus of this chapter. You will learn to use Python to implement, test and refine algorithms for biological questions, from small calculations to complete scientific workflows. Specifically, you will learn to:

use basic Python syntax and data structures;
work with Python objects and their methods;
write and run code interactively and as programs;
implement decisions and repetition with control-flow tools;
write, use, test and debug functions;
use regular expressions and scientific Python packages;
improve the efficiency of numerical code; and
combine Python with other programs in reproducible analysis and simulation workflows.

Some terminology#

What does “float” mean? You will inevitably run into some such jargon in this chapter. The main ones you need to know are (you will learn more about these along the way):

Term	Meaning
Workspace	The “environment” of your current Python session, including all variables, functions, objects, etc.
Variable	A named number, text string, boolean (`True` or `False`), or data structure that can change (more on variable and data types later)
Function	A computer procedure or routine that performs operations and returns some value(s), and which can be used again and again
Module	Variables and functions packaged into a single set of programs that can be invoked as a re-useable command (potentially with sub-commands)
Class	A way of grouping Variables and functions into a single object with specific properties that are inherited when you create its copy. Unlike modules, you can create (“spawn”) many copies of a class within a Python session or program
Object	A particular instance of a class (every object belongs to a class) that is created in a session and eventually destroyed; everything in your workspace is an object in Python!

This Module vs. Class vs. Object business is confusing. These constructs are created to make an (object-oriented) programming language like Python more flexible and user friendly (though it might not seem so to you currently!). In practice, at least for your current purposes, you will not build you own Python classes much (but will use the inbuilt Python classes).You will however write your own modules. More on all this later in this Chapter.

Note

Data “structures” vs. “objects”: You will often see the terms “object” and “data structure” in this and other chapters. A data structure is a way of organising data, such as a list or dictionary. In Python, values, functions and fitted models can all be objects. Each kind of object has its own rules for the operations that Python can perform on it.

Getting started with Python#

All examples in this chapter use Python 3. Run them with python3 unless your environment already maps python to Python 3.

$\star$ In your bash terminal (or by opening a new one with ctrl+alt+t), type:

python3

You will get a new command prompt that looks like this:

>>>

Baby steps#

Now, try some simple operations:

2 + 2 # Summation; note that comments still start with #

2 * 2 # Multiplication

2 / 2 # division

1.0

To specify an integer division, use //:

2//2

2 > 3 # logical operation

False

2 >= 2 # another one

True

The Zen of Python and good programming practices#

Type this in your Python terminal:

import this

The Zen of Python, by Tim Peters

Beautiful is better than ugly.
Explicit is better than implicit.
Simple is better than complex.
Complex is better than complicated.
Flat is better than nested.
Sparse is better than dense.
Readability counts.
Special cases aren't special enough to break the rules.
Although practicality beats purity.
Errors should never pass silently.
Unless explicitly silenced.
In the face of ambiguity, refuse the temptation to guess.
There should be one-- and preferably only one --obvious way to do it.
Although that way may not be obvious at first unless you're Dutch.
Now is better than never.
Although never is often better than *right* now.
If the implementation is hard to explain, it's a bad idea.
If the implementation is easy to explain, it may be a good idea.
Namespaces are one honking great idea -- let's do more of those!

They are Python’s programming principles. The Dutch reference is to Guido van Rossum, Python’s creator.

Use these aphorisms as prompts for discussion, not as an exhaustive checklist:

Explicit is better than implicit: state the inputs, assumptions and steps before coding. This builds on the shell-scripting treatment of a program as an explicit recipe.
Readability counts: choose clear names and layout so another person can follow the steps. See Python variables and writing readable Python code, and consult the official PEP 8 style guide.
Check before trusting a result: predict what a short expression will do, run it, and compare the output with your prediction. Try Baby steps and testing/running blocks of code.
Keep a useful work history: use Git to record changes, with clear commit messages. This supports the workflow developed in Shell scripting: inspect, correct, rerun and share.

In Functions, you will learn to break a larger task into smaller reusable steps; Unit testing then introduces formal checks for those steps.

ipython#

Let’s switch to the interactive Python shell, ipython that you installed above.

$\star$ Type ctrl+D in the terminal at the Python prompt: this will exit you from the Python shell and you will see the bash prompt again.

Now type

ipython3

After some text, you should now see the ipython prompt:

In [1]:

The ipython shell has many advantages over the bare-bones, non-interactive Python shell (with the >>> prompt). For example, as in the bash shell, TAB leads to auto-completion of a command or file name (try it).

Magic commands#

Ipython has “magic commands” (which start with % ; e.g., %run). Here are some useful magic commands:

Command	What it does
`%who`	Shows current namespace (all variables, modules and functions)
`%whos`	Also display the type of each variable; typing `%whos function` only displays functions etc.
`%pwd`	Print working directory
`%history`	Print recent commands

You can try any or all of these now. For example:

%who

this

That is, there are no objects in your workspace yet. Let’s create one:

a = 1

%who

a	 this

%whos

Variable   Type      Data/Info
------------------------------
a          int       1
this       module    <module 'this' from '/usr<...>/lib/python3.14/this.py'>

Determining an object’s type#

Another useful Ipython feature is the question mark, which can be used to find what a particular Python object is, including variables you created. For example, try: ?a

This will give you detailed information about this variable (which is an object, belonging to a particular class, because this is Python!).

You can also check a variable’s type:

type(a)

int

Tip

You can configure ipython’s environment and behavior by editing the ipython_config.py file, which is located in the .ipython directory of your home (on Linux/Ubuntu). This file does not initially exist, but you can create it by running ipython profile create [profilename] in a bash terminal. Then, edit it. For example, on Ubuntu you can

nano ~/.ipython/profile_default/ipython_config.py &

And then make the changes you want to the default ipython configuration. For example, If you don’t like the blue ipython prompt, you can type %colors linux (once inside the shell). If you want to make this color the default, then edit the ipython_config.py — search for “Set the color scheme” option in the file.

Python variables#

Now, let’s continue our Python intro. We will first learn about the Python variable types that were mentioned above. The types are:

a = 2 #integer
type(a)

int

a = 2. #Float
type(a)

float

a = "Two" #String
type(a)

str

a = True #Boolean
type(a)

bool

Thus, Python has integer, float (real numbers, with different precision levels), string and Boolean variables.

Also try ?a after defining a as a Boolean variable, and note this output in particular:

Tip

The builtins True and False are the only two instances of the class bool. The class bool is a subclass of the class int, but with only two possible values.

The idea of what a class is should be clearer now.

Note

In Python, a variable’s type is determined while the program runs (dynamic typing), as in R and unlike C or Fortran. Dynamic typing is convenient, but it can add runtime overhead. More on efficient computing later.

Python operators#

Here are common arithmetic, comparison, and Boolean operators in Python:

Operator	Meaning
`+`	Addition
`-`	Subtraction
`*`	Multiplication
`/`	True division
`**`	Exponentiation
`%`	Remainder (modulo)
`//`	Floor division
`==`	Equal to
`!=`	Not equal to
`>`	Greater than
`>=`	Greater than or equal to
`and`	Logical AND
`or`	Logical OR
`not`	Logical NOT

Use and, or, and not for Boolean logic. & is bitwise AND, while ! is not a Python operator.

At the Python prompt:

>>> 7 / 2
3.5
>>> 7 // 2
3
>>> -7 // 2
-4
>>> 7 % 2
1

/ keeps the fractional part; // rounds down to the next whole number, so -7 // 2 is -4, not -3. % gives the remainder, so 7 % 2 is 1.

Assigning and manipulating variables#

Try the following:

2 == 2

True

2 != 2

False

3 / 2

1.5

3 // 2

What happened here? This is an integer division, so the decimal part is lost.

'hola, ' + 'me llamo Samraat' #why not learn two human languages at the same time?! 

'hola, me llamo Samraat'

x = 5

x + 3

y = 2

x + y

x = 'My string'

x + ' now has more stuff'

'My string now has more stuff'

x + y

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[28], line 1
----> 1 x + y

TypeError: can only concatenate str (not "int") to str

Doesn’t work. No problem, we can convert from one type to another:

x + str(y)

'My string2'

z = '88'

x + z

'My string88'

y + int(z)

The modulo operator % is used to find the remainder after a division.

9 % 4

Writing readable Python code#

Before building complete .py programs, let’s establish some guidelines for readable Python. These practices apply to code entered in IPython, notebooks and saved scripts alike (also see the official Python style guide, PEP 8):

Wrap lines to be <80 characters long. You can use parentheses () or signal that the line continues using a backslash \
Use either 4 spaces for indentation or tabs, but not both. (Spaces are the preferred indentation method according to pep8)
Separate functions using a blank line
When possible, write comments on separate lines

Make sure you have chosen a particular indent type (space or tab) in whatever code IDE/editor you are using — indentation is all-important in python.

Tip

IDEs / code editors, by default, will typically impose consistency of which indentation (tab or 4 spaces) is used in and across python scripts. For example, 4 spaces is usually the default, and if you use a tab (easier, quicker) to indent code while writing, the editor will automatically convert it to 4 spaces. If your code editor does not do this automatically, you should be able to configure it to do so.

Furthermore,

Use “docstrings” to document how to use the code, and comments to explain why and how the code works (we will learn about docstrings soon, below)
Follow naming conventions, especially:
- a_variable (this is “snake case”; use underscores, not spaces in variable names)
- _internal_variable (allowed, but use as a module-specific variable only)
- SOME_CONSTANT
- a_function (and don’t use capital letters in function names)
Never call a variable l or O or o (why not? – you are likely to confuse it with 1 or 0!)
Use spaces around operators and after commas:

a = func(x, y) + other(3, 4)

Testing/running blocks of code#

Now that you have seen how all-important indentation of Python code is. You can test a block of code, indentation and all, by pasting it directly into the ipython terminal. Let’s try it.

Type the following code in a temporary text file:

for i in range(x):
    if i > 3: #4 spaces or 2 tabs in this case
        print(i)

Now, assign some integer value to a variable x:

x = 11

Then, paste this code at the ipython prompt (ctrl+shift+v), and hit enter:

for i in range(x):
    if i > 3: #4 spaces or 2 tabs in this case
        print(i)

Of course, this code is simple, so directly pasting works. For more complex code, you may need to use the ipython %cpaste magic function.

Python data structures#

Python variables can be stored and manipulated in:


List:	most versatile, can contain mixed data, are “mutable” (elements can can be modified), and can contain duplicates; created by enclosing in square brackets, `[ ]`
Tuple:	like a list, but “immutable” — like a read only list; created by enclosing in round brackets, `( )`
Dictionary:	a kind of “hash table” of key-value pairs where the key can be number or string and values can be any Python object, but there can be only unique key-value pairs; created by enclosing in curly brackets, `{ }`
numpy arrays:	Fast, compact, convenient for numerical computing — more on this later!

Note

What about Pandas arrays as a data structure? We will learn about them later. These are inherently slower (less computationally efficient) than numpy arrays, and are best used in the right context (e.g., they make data exploration or visualizations more convenient, or apply functions designed to run optiomally on them).

Lists#

These are the most versatile, and can contain compound data. They are “mutable”, as will be illustrated below. Try this:

MyList = [3,2.44,'green',True]

MyList[1]

2.44

MyList[0]

Note that Python “indexing” starts at 0, not 1!

MyList[4]

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[39], line 1
----> 1 MyList[4]

IndexError: list index out of range

As expected!

MyList[2] = 'blue'

MyList

[3, 2.44, 'blue', True]

MyList.append('a new item')

Note .append. This is an operation (a “method”) that can be applied to any “object” with the “class” list. You can check the type of any object:

%whos

Variable   Type      Data/Info
------------------------------
MyList     list      n=5
a          bool      True
i          int       10
this       module    <module 'this' from '/usr<...>/lib/python3.14/this.py'>
x          int       11
y          int       2
z          str       88

type(MyList)

list

print(type(MyList))

<class 'list'>

Note

In Python, every class is an object, usually created by the metaclass type. Calling type(value) reports the class of value. Classes and types are closely related, but the terms are not interchangeable.

MyList

[3, 2.44, 'blue', True, 'a new item']

del MyList[2]

MyList

[3, 2.44, True, 'a new item']

Tip

Note that in ipython you can suffix an . to a particular object (e.g., MyList.), and then hit tab to see the methods that can be applied to that object.

Tuples#

Tuples are like a list, but “immutable”, that is, a particular pair or sequence of strings or numbers cannot be modified after it is created. So a tuple is like a read-only list.

MyTuple = ("a", "b", "c")

print(MyTuple)

('a', 'b', 'c')

type(MyTuple)

tuple

MyTuple[0]

'a'

len(MyTuple)

Tuples are particularly useful when you have sets (pairs, triplets, etc) of items associated with each other that you want to store such that those associations cannot be modified.

Try this:

FoodWeb=[('a','b'),('a','c'),('b','c'),('c','c')]
FoodWeb

[('a', 'b'), ('a', 'c'), ('b', 'c'), ('c', 'c')]

So we created a list of tuples (note the use of square brackets [])

FoodWeb[0]

('a', 'b')

FoodWeb[0][0]

'a'

FoodWeb[0][0] = "bbb"

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[57], line 1
----> 1 FoodWeb[0][0] = "bbb"

TypeError: 'tuple' object does not support item assignment

Thus, tuples are “immutable”!

However, you can change a whole pairing:

FoodWeb[0] = ("bbb","ccc") 

FoodWeb[0]

('bbb', 'ccc')

In the above example, why assign these food web data to a list of tuples and not a list of lists? — because we want to maintain the species associations, no matter what — they are sacrosanct.

Thus, you cannot:

add elements to a tuple. Tuples have no inbuilt append or extend method (but you can append to them with a bit of extra work - see below).
remove elements from a tuple. Tuples have no remove method.

But you can:

find elements in a tuple, since this doesn’t change the tuple.
use the in operator to check if an element exists in the tuple or iterate over them (explained in the control flow tools sections below).

The key point is that tuples are faster than lists, as you might expect for an immutable object (it has a fixed memory space, which makes it more efficient to retrieve). using tuples also makes your code safer as it effectively write-protects data (as long as you don’t plan to modify those particular data).

Tuples may be immutable, but you can append to them by first creating an “empty space” for the new item:

a = (1, 2, []) 
a

(1, 2, [])

a[2].append(1000)
a

(1, 2, [1000])

a[2].append(1000)
a

(1, 2, [1000, 1000])

a[2].append((100,10))
a

(1, 2, [1000, 1000, (100, 10)])

You can also concatenate, slice and dice them as long as they contain a single sequence or set of items:

a = (1, 2, 3)

b = a + (4, 5, 6)
b

(1, 2, 3, 4, 5, 6)

c = b[1:]
c

(2, 3, 4, 5, 6)

b = b[1:]
b

(2, 3, 4, 5, 6)

They can be a heterogeneous set as well.

a = ("1", 2, True)
a

('1', 2, True)

Sets#

You can convert a list to an mutable “set” — an unordered collection with no duplicate elements. Once you create a set you can perform set operations on it:

a = [5,6,7,7,7,8,9,9]

b = set(a)

{5, 6, 7, 8, 9}

c = set([3,4,5,6])

b & c # intersection

{5, 6}

b | c # union

{3, 4, 5, 6, 7, 8, 9}

The key set operations in Python are:

Operation	Command
`a - b`	a.difference(b)
`a <= b`	a.issubset(b)
`a >= b`	b.issubset(a)
`a & b`	a.intersection(b)
`a` $\vert$ `b`	a.union(b)

You can also convert tuples to sets:

# Define a tuple
my_tuple = (1, 2, 3, 4, 5)

# Convert tuple to set
my_set = set(my_tuple)

# Print the set
print(my_set)

{1, 2, 3, 4, 5}

Dictionaries#

A dictionary is a set of values (any Python object) indexed by keys (string or number). So they are a bit like R lists.

genome_size = {'Homo sapiens': 3200.0, 'Escherichia coli': 4.6, 'Arabidopsis thaliana': 157.0}

 genome_size

{'Homo sapiens': 3200.0,
 'Escherichia coli': 4.6,
 'Arabidopsis thaliana': 157.0}

genome_size['Arabidopsis thaliana']

157.0

genome_size['Saccharomyces cerevisiae'] = 12.1

genome_size

{'Homo sapiens': 3200.0,
 'Escherichia coli': 4.6,
 'Arabidopsis thaliana': 157.0,
 'Saccharomyces cerevisiae': 12.1}

genome_size['Escherichia coli'] = 4.6 

genome_size

{'Homo sapiens': 3200.0,
 'Escherichia coli': 4.6,
 'Arabidopsis thaliana': 157.0,
 'Saccharomyces cerevisiae': 12.1}

Because ‘Escherichia coli’ is already in the dictionary, it is not repeated.

genome_size['Homo sapiens'] = 3201.1

genome_size

{'Homo sapiens': 3201.1,
 'Escherichia coli': 4.6,
 'Arabidopsis thaliana': 157.0,
 'Saccharomyces cerevisiae': 12.1}

Note

Tuples that contain immutable values like strings, numbers, and other tuples can be used as dictionary keys. Lists can never be used as dictionary keys, because they are mutable.

Also, every dictionary key must be unique. If you try to add a new key-value pair to a dictionary using a key that already exists, the new value will overwrite the previous value associated with that key.

Try this:

# Define a dictionary with duplicate keys
my_dict = {'a': 1, 'b': 2, 'a': 3}

# Print the dictionary
print(my_dict)

{'a': 3, 'b': 2}

In summary, the guidelines for choosing a Python data structure are:

If your elements/data are unordered and indexed by numbers use a list
If you’re defining a constant set of values (or ordered sequences) and all you’re ever going to do with them is iterate through them, use a tuple.
If you want to perform set operations on data, use a set
If they are unordered and indexed by keys (e.g., names), use a dictionary

But why not use dictionaries for everything? – because it can slow down your code!

Copying mutable objects#

Copying mutable objects can be tricky because by default, when you create a new variable based on an existing one, Python only creates a reference to the original (that is it does not create a new, duplicate variable in memory as such). To understand this, let’s see an example.

First, try this:

a = [1, 2, 3]
b = a

Here, you have not really copied, but merely created a new “tag” (like a label) for a, called b.

a.append(4)

print(a)
print(b)

[1, 2, 3, 4]
[1, 2, 3, 4]

So b changed as well! This is because b is just a “pointer” or “reference” to a, not an actual copy in memory.

Now, try:

a = [1, 2, 3]
b = a[:]  # This is a "shallow" copy; one level deep

a.append(4)
print(a)
print(b)

[1, 2, 3, 4]
[1, 2, 3]

That worked! But what about more complex lists? Try this nested list:

a = [[1, 2], [3, 4]]
b = a[:]
print(a)
print(b)

[[1, 2], [3, 4]]
[[1, 2], [3, 4]]

Now, modify a, and then inspect both a and b:

a[0][1] = 22 # Note how I accessed this 2D list
print(a)
print(b)

[[1, 22], [3, 4]]
[[1, 22], [3, 4]]

So b still got modified!

This is because shallow copy is not recursive, that is, it does not copy beyond the first level of the list, leaving the values in the nested list still linked in memory to the original object a.

The solution is to do a “deep” copy:

import copy

a = [[1, 2], [3, 4]]
b = copy.deepcopy(a)
a[0][1] = 22
print(a)
print(b)

[[1, 22], [3, 4]]
[[1, 2], [3, 4]]

So, you need to employ deepcopy to really copy an existing object or variable and assign a new name to the copy. So, in summary, shallow copying an object won’t create objects that are independent clones, i.e., the copy is not fully independent of the original. A deep copy of an object will recursively clone “child” objects (like nested parts of a list). The clone is fully independent of the original, but creating a deep copy is slower, as it involves assigning new memory space. Keep in mind that this shallow vs. deep copy business does not just apply to lists. You can copy arbitrary objects (including custom classes) with the copy module.

Note

Why python “shallow” copies objects: This is a bit of a technical detail, but important to keep in mind: python does shallow copying of mutable objects for (computing) performance considerations. By not copying the underlying object when you re-assign a mutable object to a new (“variable”) name, python avoids unnecessary memory usage. This is known as “passing by reference” (in contrast to passing by “value”, where a new variable would be actually created in memory). That does not change the fact that shallow vs. deep copying can be confusing, of course!

Python with strings#

One of the things that makes Python so useful and versatile, is that it has a powerful set of inbuilt commands to perform string manipulations. For example, try these:

s = " this is a string "
len(s) # length of s -> 18

s.replace(" ","-") # Substitute spaces " " with dashes

'-this-is-a-string-'

s.find("s") # First occurrence of s (remember, indexing starts at 0)

s.count("s")# Count the number of "s"

t = s.split() # Split the string using spaces and make a list 
t

['this', 'is', 'a', 'string']

t = s.split(" is ") # Split the string using " is " and make a list out of it
t

[' this', 'a string ']

t = s.strip() # remove trailing spaces
t

'this is a string'

s.upper()

' THIS IS A STRING '

s.upper().strip() # can perform sequential operations

'THIS IS A STRING'

'WORD'.lower() # can perform operations directy on a literal string 

'word'

Getting help#

You can do this:

?s.upper

Also try help() at the Python/ipython prompt.

Exercises: Python data structures#

Complete tuple.py and dictionary.py from the MQB code/ directory. For each task:

Predict the result before running the script.
Complete the unfinished code without changing the supplied input data.
Run the script with python3 and compare its output with your prediction.
Use clear variable names and briefly comment any decision that is not obvious from the code.

You are finished when both scripts run without errors and you can explain why a tuple or dictionary is appropriate for each task.

Loops#

What exactly is going on in the piece of code above? What is i? What does range(x) do?

Basically, this piece of code runs a loop (“loops”) over the full range of x numbers, printing each one of them.

First, let’s understand the range() function. This function generates, as the name suggests, a range of integers depending on the input to it. So, for example, range(10) generates 10 numbers, starting at 0:

for i in range(10):
    print(i)

The start point is 0 because this is Python (it will start at 1 in R, for example). Note that if you try to runrange() by itself, it will not actually produce a range of numbers. For example:

a = range(10)
a

range(0, 10)

Printing a shows a compact description of the range, not a list of all its values. range(10) represents the integers from 0 up to, but not including, 10.

A range is a reusable iterable, not an iterator or a generator. It stores a compact description of the sequence rather than all the integers, so a large range uses little memory. A loop obtains an iterator from the range to retrieve its values. Use list(a) if you actually need a list of those values; that list stores them all.

You can also use range() to generate numbers (and loop over) from a specific range of integers. For example, to generate a range from 1 to 5, do:

for i in range(1, 6):
    print(i)

Yes, it is slightly counter-intuitive that you have to use range(1, 6) to generate numbers from 1 to 5, but that’s inevitable (and something to get used to) because of the fact that Python’s indexing starts at 0!

You can also generate a set of indices that skips values using range() like so:

for i in range(2, 10, 2): # skip odd numbers
    print(i)

Play around with range a bit, and also check out its documentation. This is a very important function that you will use again and again!

Note

In Python 2, range() creates a list, while xrange() provides a lazy sequence. In Python 3, range() is lazy, like Python 2’s xrange(), but it is not simply a renamed version. Use range() in Python 3.

The variable i in our loop is the loop variable: Python assigns the next value from the iterable to it on each iteration. In for i in range(10), its first value is 0. That value happens to match the first position in this range, but a loop variable is not necessarily an index and is not the iterator. Choose a meaningful name for the values you are working with, such as height_cm when looping over plant heights.

Iterator vs Iterable#

An iterable is an object from which Python can obtain an iterator using the built-in function iter(). Lists, tuples, dictionaries, strings and ranges are examples. A list is iterable, but it is not itself an iterator.

An iterator remembers progress through values. The built-in function next() retrieves its next value; when there are no more values, it raises StopIteration. Advancing an iterator consumes values from that iterator, not from the original list. An iterator is also iterable: iter() on it returns the same iterator, rather than restarting it.

A for loop obtains an iterator from its iterable, repeatedly retrieves the next value and assigns it to the loop variable, then ends when the iterator raises StopIteration. You do not normally need to call iter() or next() yourself when writing a for loop.

To see the difference, try the list and iterator examples below:

my_iterable = [1,2,3]

type(my_iterable)

list

my_iterator = iter(my_iterable)

type(my_iterator)

list_iterator

next(my_iterator) # same as my_iterator.__next__()

next(my_iterator)

next(my_iterator)

next(my_iterator)

---------------------------------------------------------------------------
StopIteration                             Traceback (most recent call last)
Cell In[113], line 1
----> 1 next(my_iterator)

StopIteration: 

Once, when you iterated all items in an iterator, and no more data are available, and a StopIteration exception is raised.

Note

Generators are a kind of iterator. A generator expression, such as (height_cm * 2 for height_cm in heights), produces values as they are requested. A generator function uses yield to supply a value and pause until the next request. Every generator is an iterator, and every iterator is iterable; the reverse relationships are not guaranteed. iter([8, 12, 10]) is an iterator but not a generator, and range(3) is iterable but neither an iterator nor a generator.

Here is a comparison using plant heights. Predict each output before running it:

heights = [8, 12, 10]
height_iterator = iter(heights)
print(next(height_iterator))
print(list(height_iterator))
print(list(height_iterator))
print(heights)

height_generator = (height_cm * 2 for height_cm in heights)
print(next(height_generator))
print(list(height_generator))
print(list(height_generator))

height_range = range(3)
print(list(height_range))
print(list(height_range))

The output is:

8
[12, 10]
[]
[8, 12, 10]
16
[24, 20]
[]
[0, 1, 2]
[0, 1, 2]

next() consumes one value; list() consumes the remaining values from that iterator and stores them in a new list. Here, consumption means advancing the iterator, not deleting values from the original heights list. An exhausted iterator has no remaining values to return; consumption does not itself imply a memory saving. The generator also remembers its progress and cannot restart after exhaustion; make a new generator to repeat the calculation. The range, like the list, can supply a fresh iterator for another pass. Being iterable does not itself guarantee reuse: iterators, including generators, are iterable too.

Generators can save memory by producing values on demand instead of storing a complete result list, but that is distinct from consumption. For example, a loop can process each generated value without collecting every result. Calling list(height_generator), as above, instead creates and stores a list of the remaining results. The original input list heights remains in memory in this example. Not every iterator generates new values or saves memory: the list iterator above still refers to an existing list.

See Python’s official documentation on iterator types and generator expressions.

Some loops examples#

Write the following, and save them to loops.py:

# FOR loops
for i in range(5):
    print(i)

my_list = [0, 2, "geronimo!", 3.0, True, False]
for k in my_list:
    print(k)

total = 0
summands = [0, 1, 11, 111, 1111]
for s in summands:
    total = total + s
    print(total)

# WHILE loop
z = 0
while z < 100:
    z = z + 1
    print(z)

Functions#

In Python, you delineate a function (recall what a function means from the table above) by using indentation. For example:

def foo(x):
    x *= x # same as x = x*x
    print (x)
    return x

Now you will have a function object called foo in your workspace. You can check this using the %whos magic command, which lists and describes all the objects in your workspace:

%whos

Variable      Type             Data/Info
----------------------------------------
FoodWeb       list             n=4
MyList        list             n=4
MyTuple       tuple            n=3
a             range            range(0, 10)
b             list             n=2
c             set              {3, 4, 5, 6}
copy          module           <module 'copy' from '/usr<...>/lib/python3.14/copy.py'>
foo           function         <function foo at 0x7433f91e5220>
genome_size   dict             n=4
i             int              8
my_dict       dict             n=2
my_iterable   list             n=3
my_iterator   list_iterator    <list_iterator object at 0x74340d5e0910>
my_set        set              {1, 2, 3, 4, 5}
my_tuple      tuple            n=5
s             str               this is a string 
t             str              this is a string
this          module           <module 'this' from '/usr<...>/lib/python3.14/this.py'>
x             int              11
y             int              2
z             str              88

So, foo is a function stored in memory (at address given by the value 0x... in the Data/Info column), and ready to serve you!

Now “call it”:

foo(2)

Note that the first, print command only outputs the value of x to the terminal, whereas, the second return command actually outputs it so that you can “capture” and store it.

To see this distinction, let’s try the following.

def foo(x):
    x *= x # same as x = x*x
    print (x)
    return x

y = foo(2)

type(y)

int

Thus, the output of foo was stored as a new variable y.

def foo(x):
    x *= x # same as x = x*x
    print (x)
    # return x

y = foo(2)

type(y)

NoneType

So, if we don’t explicitly return the value of x, the output of foo cannot be stored.

Hints about inputs and outputs#

As your functions become more substantial, it helps readers if the definition gives a clue about the values going in and the value coming back. For example:

def double(value: float) -> float:
    return value * 2

Here, value: float says that value is expected to be a floating-point number, while -> float says that the function is expected to return one. These type hints are notes for the reader (and for checking tools); they do not convert or check a value for you.

You will meet combinations of the same notation in larger programs:

Hint	Read it as
`dict[str, str]`	A dictionary with strings for its keys and values, such as a row read from a CSV file.
`list[dict[str, str]]`	A list of those row dictionaries.
`set[str]`	A set of strings.
`tuple[A, B]`	Two returned items: one of type A and one of type B.
`dict[str, str] \| None`	Either a row dictionary or `None`.
`Iterable[dict[str, str]]`	Any object you can loop over to get row dictionaries; it need not be a list.
`-> None`	The function does not return a useful value; it may, for example, write a file.

Iterable comes from Python’s typing module. You may also see from __future__ import annotations near the top of a script; this lets Python store the hints without evaluating them immediately. Neither line changes what the function actually does.

Control flow tools#

OK, let’s get deeper into Python code. A computer script or program’s control flow is the order in which the code executes. Upto now, you have written scripts with simple control flows, with the code executing statements from the top to bottom. But very often, you want more flexible flows of commands and statements, for example, where you can switch between alternative commands depending on some condition. This is possible using control flow tools. Let’s learn Python’s control flow tools hands-on.

Conditionals#

Now that we know how to define functions in Python, let’s look at conditionals that allow you fine-grained control over the function’s operations.

Let’s start with the if statement, which allows you to control when to execute a statement. This is the basic syntax:

if <expr>:
    <statement>

<expr> is the (boolean) condition (True / False) that needs to be met for the <statement> to be executed.

Let’s look at a basic example.

x = 0; y = 2

if x < y: # "Truthy"
    print('yes')

yes

if x:
    print('yes')

No output appears because x is zero. Zero is an integer, but it is falsy, so the body of if x: does not run.

if x==0:
    print('yes')

yes

if y:
    print('yes')

yes

if y == 2:
    print('yes')

yes

x = True # make x boolean

if x: # now it works
    print('yes')

yes

if x == True:
    print('yes')

yes

We can now try combining functions and the if conditional.

$\star$ run the following functions one by one, by pasting the block in the ipython command line. First, type all them all in a script and save it as cfexercises1.py. Then you can send them block by block easily to the command line assuming you have set your code editor to allow selections of code to be sent to terminal directly using a key binding (typically , ctrl+enter).

def foo_1(x):
    return x ** 0.5

def foo_2(x, y):
    if x > y:
        return x
    return y

def foo_3(x, y, z):
    if x > y:
        x, y = y, x
    if x > z:
        x, z = z, x
    if y > z:
        y, z = z, y
    return [x, y, z]

def foo_4(x):
    result = 1
    for i in range(1, x + 1):
        result = result * i
    return result

def foo_5(x): # a recursive function that calculates the factorial of x
    if x == 1:
        return 1
    return x * foo_5(x - 1)
     
def foo_6(x): # Calculate the factorial of x in a different way; no if statement involved
    facto = 1
    while x >= 1:
        facto = facto * x
        x = x - 1
    return facto

Think about what each of the foo_x function does before running it. Note that foo_5 is a recursive function, meaning that the function calls itself. Note that the factorial of 0 is also 1, i.e., $0! = 1$, but foo_5 and foo_6 cannot handle this. You may wish to modify the functions to fix this!

More examples of loops and conditionals combined#

$\star$ Write the following functions and save them to cfexercises2.py:

########################
def hello_1(x):
    for j in range(x):
        if j % 3 == 0:
            print('hello')
    print(' ')

hello_1(12)

########################
def hello_2(x):
    for j in range(x):
        if j % 5 == 3:
            print('hello')
        elif j % 4 == 3:
            print('hello')
    print(' ')

hello_2(12)

########################
def hello_3(x, y):
    for i in range(x, y):
        print('hello')
    print(' ')

hello_3(3, 17)

########################
def hello_4(x):
    while x != 15:
        print('hello')
        x = x + 3
    print(' ')

hello_4(0)

########################
def hello_5(x):
    while x < 100:
        if x == 31:
            for k in range(7):
                print('hello')
        elif x == 18:
            print('hello')
        x = x + 1
    print(' ')

hello_5(12)

# WHILE loop with BREAK
def hello_6(x, y):
    while x: # while x is True
        print("hello! " + str(y))
        y += 1 # increment y by 1 
        if y == 6:
            break
    print(' ')

hello_6 (True, 0)

Try to predict how many times “hello” will be printed before testing each of these functions.

Note

Note how, in the last function above, the break directive exits the loop when the condition is met. If you did not have this, you would get an infinite loop! (and would need to use Ctrl+c to stop it). Note also that break only breaks out of the current loop. It does not stop the execution of the rest of the code that may be in that program or script.

From an algorithm to Python#

Return to the plant-height question from the introduction: how many of the heights 8, 12, 10 and 15 cm exceed 10 cm? First trace this complete program by hand, then run it. Write down the value of above_threshold after each height.

heights_cm = [8, 12, 10, 15]
above_threshold = 0
for height_cm in heights_cm:
    if height_cm > 10:
        above_threshold += 1
print(above_threshold)

The counter starts at zero. The loop examines each height, and the conditional decides whether to increase the counter:

Height (cm)	Is `height_cm > 10` true?	Counter after this height
8	No	0
12	Yes	1
10	No	1
15	Yes	2

The program prints 2, matching the answer predicted before coding. Now change > to >= to answer a different question: how many heights are at least 10 cm? Predict the answer before rerunning it, and explain which row changes.

Comprehensions#

Python offers a way to combine loops and logical tests / conditionals in a single line of code to transform any iterable object (list, set, or dictionary, over which you can iterate) into another object, after performing some operations on the elements in the original object. That is, they are a compact way to create a new list, dictionary or object from an existing one. As you might expect, there are three types of comprehensions, each corresponding to what the target object is (list, set, dictionary).

Let’s look at how list comprehensions work:

x = [i for i in range(10)]
print(x)

[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

This is the same as writing the following loop:

x = []
for i in range(10):
    x.append(i)
print(x)

[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

Here’s another example:

x = [i.lower() for i in ["LIST","COMPREHENSIONS","ARE","COOL"]]
print(x)

['list', 'comprehensions', 'are', 'cool']

Which is same as the loop:

x = ["LIST","COMPREHENSIONS","ARE","COOL"]
for i in range(len(x)): # explicit loop
    x[i] = x[i].lower()
print(x)

['list', 'comprehensions', 'are', 'cool']

Or this loop:

x = ["LIST","COMPREHENSIONS","ARE","COOL"]
x_new = []
for i in x: # implicit loop
    x_new.append(i.lower())
print(x_new)

['list', 'comprehensions', 'are', 'cool']

How about a nested loop? Let’s try an example where we “flatten” a 3 × 3 matrix into a one-dimensional list (all numbers across rows in a single sequence):

matrix = [[1,2,3],[4,5,6],[7,8,9]]
flattened_matrix = []
for row in matrix:
    for n in row:
        flattened_matrix.append(n)
print(flattened_matrix)

[1, 2, 3, 4, 5, 6, 7, 8, 9]

A list comprehension to do the same:

matrix = [[1,2,3],[4,5,6],[7,8,9]]
flattened_matrix = [n for row in matrix for n in row]
print(flattened_matrix)

[1, 2, 3, 4, 5, 6, 7, 8, 9]

This list comprehension uses a nested loop, so its syntax is a bit harder to understand than the list comprehension using a single loop. So let’s dissect it:

We need to loop over each row in the given 2D list and append it to the previous row in the new flattened_matrix object (list). To understand how the list comprehension achieves this, let’s divide it into its three parts:

flatten_matrix = [n
                  for row in matrix
                  for n in row]

The first line n states what we want to append in the new flattened list.
The second line is the outer loop (same as for row in matrix: in the nested loop above); it returns the rows (sub-lists) inside the matrix one by one ([1, 2, 3], [4, 5, 6], [7, 8, 9]
The third line is the inner loop (same as for n in row: in the nested loop above); it returns all the values inside the row (the sub-list) and appends them to the growing flattened list flattened_matrix. That is, if row = [1, 2, 3], for n in row returns 1, 2, 3 one by one, and appends it to flattened_matrix in their original sequence.

Set and Dictionary comprehensions work in an analogous way. For example, create a set of all the first letters in a sequence of words using a loop:

words = ["These", "are", "some", "words"]

first_letters = set()
for w in words:
    first_letters.add(w[0])

print(first_letters)

type(first_letters)

{'T', 'a', 's', 'w'}

set

Note that sets are unordered (the first letters don’t appear in the order you might expect).

Now, the same as a set comprehension:

words = ["These", "are", "some", "words"]

first_letters = {w[0] for w in words} # note the curly brackets

print(first_letters)

type(first_letters)

{'T', 'a', 's', 'w'}

set

Now, type the following in a script file called oaks.py and test it:

## Finds just those taxa that are oak trees from a list of species

taxa = [ 'Quercus robur',
         'Fraxinus excelsior',
         'Pinus sylvestris',
         'Quercus cerris',
         'Quercus petraea',
       ]

def is_an_oak(name):
    return name.lower().startswith('quercus ')

##Using for loops
oaks_loops = set()
for species in taxa:
    if is_an_oak(species):
        oaks_loops.add(species)
print(oaks_loops)

##Using list comprehensions   
oaks_lc = set([species for species in taxa if is_an_oak(species)])
print(oaks_lc)

##Get names in UPPER CASE using for loops
oaks_loops = set()
for species in taxa:
    if is_an_oak(species):
        oaks_loops.add(species.upper())
print(oaks_loops)

##Get names in UPPER CASE using list comprehensions
oaks_lc = set([species.upper() for species in taxa if is_an_oak(species)])
print(oaks_lc)

Carefully compare the looping vs list comprehension way for the two tasks (find oak tree species names and get names in upper case) to make sure you understand what’s going on.

Note

Don’t go mad with list comprehensions — code readability is more important than squeezing lots into a single line! They can also make your code run more slowly or use more memory in some cases (we will learn about this more in the second python Chapter).

Exercises: loops and comprehensions#

Complete lc1.py and lc2.py from the MQB code/ directory. Each file asks you to solve the same transformation using a conventional loop and a comprehension.

For each result, first work out a small example by hand, then run the script with python3. Confirm that the loop and comprehension produce equal results. Keep both implementations so their readability can be compared.

You are finished when both scripts run without errors, the paired results agree, and you can explain the order of the clauses in each comprehension.

Variable scope#

One important thing to note about functions, in any programming language, is that variables created inside functions are invisible outside of it, nor do they persist once the function has run unless they are explicitly returned. These are called “local” variables, and are only accessible inside their function.

Here is an example. First, type and run this block of code:

i = 1
x = 0
for i in range(10):
    x += 1
print(i)
print(x)

9
10

Let’s break this down to understand the outputs and what’s going on in terms of variable scope:

i = 1
x = 0

i, x # or print(i, x)

(1, 0)

So we created two new variables, which are now in memory (your current “workspace”).

for i in range(10):
    x += 1

i, x

(9, 10)

So i, which was originally stored with a value of 1 in your workspace has now been updated to 9 (because the loop ran over all the numbers from 0 to 9), and x is has been updated to 10 due to the successive summation that the loop performs.

That is, the operations on i and x inside the loop were in fact on the variables in the main workspace (they were changed everywhere).

Now, let’s encapsulate the above code in a function:

i = 1
x = 0
def a_function(y):
    x = 0
    for i in range(y):
        x += 1
    return x
a_function(10)
print(i)
print(x)

1
0

So i and x did not get updated in your workspace despite the fact that they were within the function! This is because the scope of the variables inside a function is retricted to within that function. That is, within the function,i and x were created and modified in a separate location in your compiter’s memory, leaving the original variables with the same names unchanged.

Again, let’s break this code down to understand whats’s going on.

First run this in your ipython commandline:

i = 1
x = 0
def a_function(y):
    x = 0
    for i in range(y):
        x += 1
    return x

Now try whos in your ipython3 commandline, which will show you that a_function has indeed been saved as a Python function object (and so is raring to go!).

Now,

a_function(10)

This returns 10 because you asked the function to, resturn the final value of x with return x

But x in your main workspace is still 0 because even if you asked the function to return x, that made it pront x that was modified within the function to the screen, but this new, modified x is not actually saved to oyur main workspace (i.e., the original x in your main workspace remains untouched).

To replace the old x with the new one resulting from running our function, you have to explicitly “catch” it by reassigning the variable name x to the output of the function:

x = a_function(10)

Thus, the things to note from this second example:

Both x and i are variables localised to the function (their scope is limited to within it)
x was only updated in the main workspace, outside the function, when it was explicitly returned from the function and you explicitly reassigned the variable name x to the function’s output.
i remained unchanged outside the function because it was not returned.

Try returning both x and i (with a return x,y instead of return x in the function).

Global variables#

In contrast, you can designate certain variables to be “global” so that they visible both inside and outside of functions in Python, like any other programming language.

To understand this, let’s look at an example.

First try this:

_a_global = 10 # a global variable

if _a_global >= 5:
    _b_global = _a_global + 5 # also a global variable
    
print("Before calling a_function, outside the function, the value of _a_global is", _a_global)
print("Before calling a_function, outside the function, the value of _b_global is", _b_global)

def a_function():
    _a_global = 4 # a local variable
    
    if _a_global >= 4:
        _b_global = _a_global + 5 # also a local variable
    
    _a_local = 3
    
    print("Inside the function, the value of _a_global is", _a_global)
    print("Inside the function, the value of _b_global is", _b_global)
    print("Inside the function, the value of _a_local is", _a_local)
    
a_function()

print("After calling a_function, outside the function, the value of _a_global is (still)", _a_global)
print("After calling a_function, outside the function, the value of _b_global is (still)", _b_global)
print("After calling a_function, outside the function, the value of _a_local is ", _a_local)

Before calling a_function, outside the function, the value of _a_global is 10
Before calling a_function, outside the function, the value of _b_global is 15
Inside the function, the value of _a_global is 4
Inside the function, the value of _b_global is 9
Inside the function, the value of _a_local is 3
After calling a_function, outside the function, the value of _a_global is (still) 10
After calling a_function, outside the function, the value of _b_global is (still) 15

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[152], line 25
     21 a_function()
     22 
     23 print("After calling a_function, outside the function, the value of _a_global is (still)", _a_global)
     24 print("After calling a_function, outside the function, the value of _b_global is (still)", _b_global)
---> 25 print("After calling a_function, outside the function, the value of _a_local is ", _a_local)

NameError: name '_a_local' is not defined

The things to note from this example:

Although _a_global was overwritten inside the function, what happened inside the function remained inside the function (What happens in Vegas…)
The variable _a_local does not persist outside the function (therefore you get the NameError at the end)
Also note that _a_global is just a naming convention – nothing special about this variable as such.

Of course, if you assign a variable outside a function, it will be available inside it even if you don’t assign it inside that function:

_a_global = 10

def a_function():
    _a_local = 4
    
    print("Inside the function, the value _a_local is", _a_local)
    print("Inside the function, the value of _a_global is", _a_global)
    
a_function()

print("Outside the function, the value of _a_global is", _a_global)

Inside the function, the value _a_local is 4
Inside the function, the value of _a_global is 10
Outside the function, the value of _a_global is 10

So _a_global was available to the function, and you were able to use it in the print command.

If you really want to modify or assign a global variable from inside a function (that is, and make it available outside the function), you can use the global keyword:

_a_global = 10

print("Before calling a_function, outside the function, the value of _a_global is", _a_global)

def a_function():
    global _a_global
    _a_global = 5
    _a_local = 4
    
    print("Inside the function, the value of _a_global is", _a_global)
    print("Inside the function, the value _a_local is", _a_local)
    
a_function()

print("After calling a_function, outside the function, the value of _a_global now is", _a_global)

Before calling a_function, outside the function, the value of _a_global is 10
Inside the function, the value of _a_global is 5
Inside the function, the value _a_local is 4
After calling a_function, outside the function, the value of _a_global now is 5

So, using the global specification converted _a_global to a truly global variable that became available outside that function (overwriting the original _a_global).

The global keyword also works from inside nested functions, but it can be slightly confusing:

def a_function():
    _a_global = 10

    def _a_function2():
        global _a_global
        _a_global = 20
    
    print("Before calling a_function2, value of _a_global is", _a_global)

    _a_function2()
    
    print("After calling a_function2, value of _a_global is", _a_global)
    
a_function()

print("The value of a_global in main workspace / namespace now is", _a_global)

Before calling a_function2, value of _a_global is 10
After calling a_function2, value of _a_global is 10
The value of a_global in main workspace / namespace now is 20

That is, using the global keyword inside the inner function _a_function2 resulted in changing the value of _a_global in the main workspace / namespace to 20, but within the scope of _a_function, its value remained 10!

Compare the above with this:

_a_global = 10

def a_function():

    def _a_function2():
        global _a_global
        _a_global = 20
    
    print("Before calling a_function2, value of _a_global is", _a_global)

    _a_function2()
    
    print("After calling a_function2, value of _a_global is", _a_global)

a_function()

print("The value of a_global in main workspace / namespace is", _a_global)

Before calling a_function2, value of _a_global is 10
After calling a_function2, value of _a_global is 20
The value of a_global in main workspace / namespace is 20

Now, because _a_global was defined in advance (outside the first function), when a_function was run,

This value was “inherited” within a_function from the main workspace / namespace,
It was then given a global designation in the inner function _a_function2,
And then, in the inner function _a_function2, when it was changed to a different value, it was modified everywhere (both within the a_function’s scope/namespace and main workspace / namespace) .

Warning

In general, avoid assigning globals because you run the risk of “exposing” unwanted variables to all functions within your workspace / namespace. Furthermore, avoid assigning globals within functions or sub-functions, as we did in the last two examples above!

Tip

But in some cases you may find it useful to assign one or more global variables that are shared across multiple modules/functions. You can do this by assigning those variables as global at the start of the script/program, but a better, safer option is to create a separate module (say, called config.py) to hold the global variables and then import it where needed.

$\star$ Collect all blocks of code above illustrating variable scope into one script called scope.py and test it (run and check for errors).

Importance of the `return` directive#

In the context of scope of variables, it is also important to keep in mind that in Python, arguments are passed to a function by assignment. This is a bit of a technical detail that we don’t need to go into here, but basically, in practice, this means that for mutable objects such as lists, unless you do something special, if a function modifies the (mutable) variable inside it, the original variable outside the function remains unchanged.

Let’s look at an example to understand this:

def modify_list_1(some_list):
    print('got', some_list)
    some_list = [1, 2, 3, 4]
    print('set to', some_list)

my_list = [1, 2, 3]

print('before, my_list =', my_list)

before, my_list = [1, 2, 3]

modify_list_1(my_list)

got [1, 2, 3]
set to [1, 2, 3, 4]

print('after, my_list =', my_list)

after, my_list = [1, 2, 3]

The original list remains the same even though it is changed inside the function, as you would expect (what happens in Vegas…)

This is where the return directive becomes important. Now modify the function to return the value of the input list:

def modify_list_2(some_list):
    print('got', some_list)
    some_list = [1, 2, 3, 4]
    print('set to', some_list)
    return some_list

my_list = modify_list_2(my_list)

got [1, 2, 3]
set to [1, 2, 3, 4]

print('after, my_list =', my_list)

after, my_list = [1, 2, 3, 4]

So now the original my_list is changed because you explicitly replaced it. This reinforces the fact that explicit return statements are important.

And if we do want to modify the original list in place, use append:

def modify_list_3(some_list):
    print('got', some_list)
    some_list.append(4) # an actual modification of the list
    print('changed to', some_list)

my_list = [1, 2, 3]

print('before, my_list =', my_list)

before, my_list = [1, 2, 3]

modify_list_3(my_list)

got [1, 2, 3]
changed to [1, 2, 3, 4]

print('after, my_list =', my_list)

after, my_list = [1, 2, 3, 4]

That did it. So append will actually change the original list object. However, the fact still remains that you should use a return statement at the end of the function to be safe and be able to capture the output (and use it to replace an existing variable if needed).

Note

returning a None: Even if you do not add a return directive at the end of a function, python does in fact return something: a None value, which stands for a NULL value (no value at all; so that something is a nothing!). You can use an return None in a python function if you want to be explicit, or use just return to completely end the execution of the code (like exit in a shell script). This is different from using the break directive, which you were introduced above under control flow tools.

Exercises: functions and control flow#

Return to cfexercises1.py, which you created in the control-flow section. Make every foo_* function accept arguments and return a result without calling input().

A Python file can be executed directly or imported by another file. Python sets its special __name__ value to "__main__" only during direct execution. The guard if __name__ == "__main__": checks that value, so calls indented under it run during direct execution but not during import. The next section develops this structure in detail.

Add representative calls under that guard. Include boundary checks where they matter—for example, confirm that your factorial functions handle zero—and compare functions that are intended to calculate the same result.

From the coursework project root, run the file with:

python3 code/cfexercises1.py

You are finished when the calls produce the results you predicted, equivalent implementations agree, and importing the file does not run the example calls.

Writing and running Python programs#

This section connects the functions you have written to complete, reproducible Python programs. A program is code written to perform a task; a script is a program saved in a file and commonly run from a shell. A .py file is also a module: Python can import it so another script, test or notebook can reuse its definitions. These terms describe related roles rather than mutually exclusive kinds of file. A function defined with def is contained in a module; it is not itself a module.

What to study first

Work through Running Python scripts, Structuring Python programs, the explanation of the main guard under Components of the Python program, and Modules and control_flow.py. Unless an exercise says otherwise, run commands from the coursework project root and include code/ in the script name, as in python3 code/example.py. This gives every exercise the same starting point: files in data/, code/ and results/ are all named from the project root.

Passing inputs from the command line is an extension, not a prerequisite for writing reusable functions. The shebang and IPython execution examples are also further reference.

Supplied-file links in this chapter are bundled with the book build, so they match the listings in the version you are reading. Type the short introductory examples when instructed; download longer scripts and datasets directly into the stated code/ or data/ directory.

Modules and code compartmentalization#

A module can contain reusable functions, constants and other definitions. Importing it lets another script or notebook use those definitions without copying them. Keeping calculations in functions and startup activity in a small guarded block makes code easier to test and reuse.

There are several common import forms:

import my_module makes the module available; call a function as my_module.my_function().
from my_module import my_function imports one named function; call it as my_function().
import my_module as mm gives the module a shorter local name; call its function as mm.my_function().

Avoid from my_module import *: it can replace existing names and hides where definitions came from. The later control_flow.py example uses explicit imports.

Running Python scripts#

Instead of pasting code at an interactive Python prompt, save related instructions in a script so you can run the same program again. Open the coursework project folder in VS Code and type the following as code/MyExampleScript.py. Because this is a short first example, typing it yourself is part of the exercise. If you need to recover or compare an exact copy, use the supplied MyExampleScript.py:

"""A first example of saving and running Python code."""


def square(x):
    """Return the square of x."""
    return x * x


result = square(2)
print(f"The square of 2 is {result}")

Using the UNIX shell#

In a shell, use pwd and ls to confirm that you are at the coursework project root. You should see the code/, data/ and results/ directories. Then run the script with:

python3 code/MyExampleScript.py

The shell starts a fresh Python process, Python executes the file from top to bottom, and the process then ends. The function calculates and returns a value; the separate print call displays it. This first script still performs that call at the top level. The next example introduces the main guard so importing reusable functions does not also run the program’s startup work. If the command fails, read the traceback before changing directory or trying a different execution method. Do not type this shell command at Python’s >>> prompt.

Let the script find the project#

Running programs from the project root keeps commands and file locations predictable. A longer program can also work out where its own project lives. Python’s pathlib module makes this quite readable:

from pathlib import Path

project_root = Path(__file__).resolve().parents[1]
input_file = project_root / "data" / "example.csv"

In a saved script, __file__ is the script’s filename. If that script is in code/, .parents[1] takes us up to the project root. The / then joins the pieces of the path; it does not open the file. This is more portable than writing a personal absolute path such as /home/yourname/.... (__file__ is normally unavailable in an interactive prompt or notebook cell.)

When a program needs an output directory, output_file.parent.mkdir(parents=True, exist_ok=True) creates it if necessary and is also happy if it already exists. Other filesystem problems will still be reported.

Optional: running the script in IPython

Inside IPython started from the project root, %run code/MyExampleScript.py executes the file using IPython’s special %run command. This can be useful for interactive exploration, but it does not show how normal Python imports work. Use the shell command above when following this chapter and when checking that the program runs from a clean start.

Structuring Python programs#

A short script can contain calculations and output in one sequence. As a program grows, separate reusable functions from the small amount of startup code that coordinates them. This structure allows a test or another module to import the functions without automatically running the complete program.

The example below is a boilerplate, or reusable structural template, similar to the template introduced in the shell scripting chapter. For now, focus on three parts: definitions, a main function that coordinates the program, and the final if __name__ == "__main__": guard. The shebang, metadata, command-line arguments and exit status are explained later for anyone who needs them.

Type the displayed example as code/boilerplate.py in the coursework project. If you need to recover or compare an exact copy, use the supplied boilerplate.py.

#!/usr/bin/env python3

"""A minimal structure for an import-safe Python program."""

def main():
    """Run the program's main task."""
    print("This is a boilerplate")

if __name__ == "__main__":
    main()

From a shell at the project root, run:

python3 code/boilerplate.py

You should get:

This is a boilerplate

This is the usual way to run the script in this chapter. It starts a fresh process and works whether or not you have previously opened IPython.

The next two cells prepare an optional IPython session started at the project root. The first runs the file; the second tells this interactive session to look in code/ when you later practise importing modules. Normal shell commands such as python3 code/test_control_flow.py do this automatically because Python looks beside the script being run. You do not need either cell for the later testing or observations exercises.

import sys
sys.path.insert(0, "code")

Components of the Python program#

Now examine the roles of the boilerplate’s main components.

The `main` function and main guard#

main() is an ordinary function used by convention to coordinate a program. Defining it does not run it. When Python executes a file directly, Python sets that module’s special __name__ value to "__main__"; when another file imports the module, Python sets __name__ to the module’s import name. The guard checks that value:

if __name__ == "__main__":
    main()

Code indented under the guard therefore runs during direct execution but is skipped during import. The guard does not set __name__, and it does not stop other unguarded top-level statements from running during import. Keep reusable function definitions above the guard and keep startup actions—such as loading a particular file or printing a report—inside main() or the guarded block.

Module and function docstrings#

A docstring is a string placed as the first statement in a module, function or class. Python stores it in __doc__, so tools such as help() can display it. A module docstring states the file’s purpose; a function docstring describes its inputs, returned value and important behaviour. A # comment instead explains a source-code decision that is not obvious from the names and structure. A docstring is documentation, not a special kind of comment.

After importing the example in the optional IPython session prepared above, help(boilerplate) displays the module docstring. IPython also supports the convenience syntax boilerplate?.

Further reference: the shebang#

The first line #!/usr/bin/env python3 is a shebang. On Unix-like systems it tells the operating system which interpreter to use only when an executable file is invoked directly, for example as ./code/boilerplate.py. The env command searches the current PATH for python3; it does not specify a fixed installation location. The shebang is ignored when you run python3 code/boilerplate.py, so you do not need executable permissions or direct invocation for this chapter. A shebang does not compile Python into a standalone application.

import boilerplate

help(boilerplate)

Help on module boilerplate:

NAME
    boilerplate - A minimal structure for an import-safe Python program.

FUNCTIONS
    main()
        Run the program's main task.

FILE
    /home/mhasoba/Documents/Projects/Software/MulQuaBio/MQB/content/code/boilerplate.py

Note

Docstrings and comments have different jobs. Use docstrings to document how a module or function is used. Use comments sparingly to explain a non-obvious decision. Neither replaces clear names, small functions or tests.

Special names and ordinary names#

Python defines some special dunder names with two leading and two trailing underscores, including __name__ and __doc__. Learn the meaning of a special name before using it, and do not invent new dunder names for ordinary variables. This is different from a single leading underscore, which is a naming convention for an internal implementation detail rather than a reserved Python variable.

Functions and modules are different#

def begins a function definition, and its body is indented using four spaces. A function is one reusable operation. A module is the complete .py file that can contain several function definitions, constants and other statements. Import a module to reuse its functions; a function is not a module or submodule.

In the optional IPython session prepared above, import the example module with:

import boilerplate

Typing the module name at an interactive prompt displays its representation; it does not run the program again. Use help(boilerplate) for its documentation or call a named function explicitly.

boilerplate

<module 'boilerplate' from '/home/mhasoba/Documents/Projects/Software/MulQuaBio/MQB/content/code/boilerplate.py'>

When python3 code/boilerplate.py executes the file directly, Python has already set __name__ to "__main__", so the condition is true and the guarded call to main runs. If another file executes import boilerplate, Python sets the imported module’s __name__ to "boilerplate"; the condition is false and the guarded call is skipped. Function definitions and any other unguarded top-level statements are still processed during import.

The following small diagnostic makes those two __name__ values visible. Type it as code/using_name.py, or use the supplied using_name.py to recover or compare an exact copy:

#!/usr/bin/env python3

"""Demonstrate the value Python assigns to __name__."""

def module_name():
    """Return the name Python assigned to this module."""
    return __name__

def main():
    """Report the module name during direct execution."""
    print(f"This module's name is: {module_name()}")

if __name__ == "__main__":
    main()

From the project root, python3 code/using_name.py reports __main__. The optional IPython example above does the same. Both commands execute the file directly, so the guarded call to main() runs.

%run code/using_name.py

This module's name is: __main__

Now import the same file in that IPython session, then call its function. Importing alone produces no report; the explicit function call returns the module name using_name:

import using_name
using_name.module_name()

'using_name'

Tip

The official Python tutorial on modules provides a longer reference on imports, module search paths and packages.

A program-with-control-flows example#

control_flow.py brings earlier functions into one module and uses a small startup section to display example results. Read it in two layers: the functions calculate and return values; main selects example inputs and presents results. This separation is the pattern later used for data loading, summaries and tests.

Download the supplied control_flow.py into your coursework project’s code/ directory. The complete listing is shown below so that you can inspect and annotate it without switching to the MQB repository. The indented lines beginning with >>> are documentation examples that can also be checked automatically; the testing section explains them shortly.

#!/usr/bin/env python3

"""Small, reusable functions that demonstrate control flow."""


def even_or_odd(x):
    """Return a description of whether integer x is even or odd.

    >>> even_or_odd(10)
    '10 is Even!'
    >>> even_or_odd(5)
    '5 is Odd!'
    >>> even_or_odd(-2)
    '-2 is Even!'
    """
    if not isinstance(x, int):
        raise TypeError("x must be an integer")
    if x % 2 == 0:
        return f"{x} is Even!"
    return f"{x} is Odd!"


def largest_divisor_five(x):
    """Return the largest divisor of x among 2, 3, 4 and 5, or None."""
    if not isinstance(x, int):
        raise TypeError("x must be an integer")
    for divisor in (5, 4, 3, 2):
        if x % divisor == 0:
            return divisor
    return None


def is_prime(x):
    """Return True when integer x is prime, and False otherwise."""
    if not isinstance(x, int):
        raise TypeError("x must be an integer")
    if x < 2:
        return False
    for divisor in range(2, x):
        if x % divisor == 0:
            return False
    return True


def find_all_primes(limit):
    """Return a list of all prime integers from 2 through limit."""
    if not isinstance(limit, int):
        raise TypeError("limit must be an integer")
    return [number for number in range(2, limit + 1) if is_prime(number)]


def count_above(heights_cm, threshold_cm):
    """Return the number of heights strictly above a threshold."""
    count = 0
    for height_cm in heights_cm:
        if height_cm > threshold_cm:
            count += 1
    return count


def main():
    """Display representative results from the functions above."""
    print(even_or_odd(22))
    print(even_or_odd(33))
    print(f"Largest divisor among 2, 3, 4 and 5: {largest_divisor_five(120)}")
    print(f"Is 59 prime? {is_prime(59)}")
    print(f"Primes up to 20: {find_all_primes(20)}")


if __name__ == "__main__":
    main()

Now run the program from the project root:

python3 code/control_flow.py

The next cell shows how to do the same thing with IPython’s optional %run command.

%run code/control_flow.py

22 is Even!
33 is Odd!
Largest divisor among 2, 3, 4 and 5: 5
Is 59 prime? True
Primes up to 20: [2, 3, 5, 7, 11, 13, 17, 19]

After the optional IPython %run code/control_flow.py cell above, IPython places the file’s definitions in the interactive namespace, so this call works in that same session. This is a %run convenience, not what happens after the core shell command python3 code/control_flow.py:

even_or_odd(11)

'11 is Odd!'

Ordinary shell execution starts a separate process; when that process ends, its functions do not become available in a different Python session. To reuse a function from another .py file, import its module or import the named function explicitly, as shown below.

Tip

An f-string inserts the value of an expression written inside braces. If x is 5, then f"The value of x is {x}." produces "The value of x is 5.". Use f-strings for clear labelled output; keep calculations in functions rather than embedding complex logic in a display string.

Modules and `control_flow.py`#

Importing control_flow executes the module’s definitions and makes them available through the module name. Because its demonstration calls belong under the main guard, importing it should not run the complete demonstration or print its report.

The core reusable pattern appears in the later separate test script: from the project root, python3 code/test_control_flow.py imports its sibling control_flow.py. The following interactive example is optional and assumes that you have already used %run code/control_flow.py in the same IPython session:

import control_flow
control_flow.even_or_odd(5)
control_flow.largest_divisor_five(100)

To import one named function directly, use an explicit from import. Prefer either this form or the qualified control_flow.even_or_odd(...) form over from control_flow import *, which hides where names originated:

from control_flow import even_or_odd
even_or_odd(5)

'5 is Odd!'

Passing inputs from the command line#

Importing a function lets you call it with different inputs from Python. But what if you want to run a program from the shell with a different input each time, without editing its code? You can supply values after the script’s filename. These are command-line arguments: they let the same program work with different filenames, species names or numerical settings.

Inspecting `sys.argv`#

After import sys, sys.argv is a list of strings describing how the current script was invoked. sys.argv[0] identifies the invoked script; elements from index 1 onwards are the arguments supplied after the script name. These are text values, even when an argument looks like a number. argv is traditionally short for argument vector; it is not a collection of Python variable names.

Type the displayed example as code/sysargv.py, or use the supplied sysargv.py to recover or compare an exact copy:

#!/usr/bin/env python3

"""Display the command-line arguments supplied to this program."""

import sys


def main():
    """Report the script name and its command-line arguments."""
    print(f"Script name: {sys.argv[0]}")
    print(f"Number of arguments: {len(sys.argv) - 1}")
    print(f"Arguments: {sys.argv[1:]}")


if __name__ == "__main__":
    main()

From a shell at the project root, compare these direct executions:

python3 code/sysargv.py
python3 code/sysargv.py oak 12

The second run has two user-supplied arguments, the strings "oak" and "12". Python does not pass sys.argv into main() automatically. A program may read sys.argv as above, or explicitly call a function such as main(sys.argv); in the latter form, argv is simply an ordinary parameter name.

Tip

Direct indexing is suitable for this small demonstration. For a real command-line interface, use Python’s standard-library argparse module to define arguments, convert values, validate input and generate --help text. Do not add a command-line interface when a task does not require one.

Reporting success or failure#

What should happen if someone runs the program without a required argument? Printing a helpful message tells the person what went wrong. But the program might also be running as one step in a shell script, where nobody is watching the output. We therefore need a way to tell the shell whether it completed successfully.

A process reports an integer exit status to the shell: zero conventionally means success and a non-zero value means failure. In the following example, exactly one argument is required. Save it as code/require_one_argument.py:

import sys


def main():
    if len(sys.argv) != 2:
        print("Usage: python3 code/require_one_argument.py VALUE", file=sys.stderr)
        return 1
    print(f"Received: {sys.argv[1]}")
    return 0


if __name__ == "__main__":
    sys.exit(main())

Run both paths from a shell at the project root; echo $? displays the status of the immediately preceding command:

python3 code/require_one_argument.py
echo $?
python3 code/require_one_argument.py oak
echo $?

The missing-argument run prints usage information to standard error and reports status 1; the valid run prints Received: oak and reports status 0. Notice that returning a value from a function does not automatically set the program’s exit status. Here, main() returns an integer, and sys.exit(main()) uses that integer as the exit status when the Python process ends. This lets a shell script check whether the program succeeded before deciding what to do next.

Note

Checking exit status in IPython? The commands above belong at the shell prompt, not the IPython prompt. With %run code/require_one_argument.py, the script runs within the IPython process rather than in a separate Python process. Without an argument, you will see SystemExit: 1. Don’t panic! This is the expected failure signal: sys.exit(1) raises a SystemExit exception, which IPython catches so that your interactive session stays open.

If you then type !echo $?, you may see 0. Has the script somehow succeeded? No: ! starts a separate shell, so its $? does not refer to the preceding %run.

Use the shell commands above to check the exit status. Alternatively, from IPython at the project root, run the script and check its status in the same shell invocation:

!python3 code/require_one_argument.py; echo $?
!python3 code/require_one_argument.py oak; echo $?

The first line prints the usage message followed by 1; the second prints Received: oak followed by 0. Keep each pair on one line: two separate ! commands do not share a shell. This distinction applies on both macOS and Linux.

Python input and output#

Let’s learn to import and export data in Python (and write code to do it).

Make a text file called test.txt in your coursework project’s root-level sandbox/ directory with the following content (including the empty lines):

First Line
Second Line

Third Line

Fourth Line

Then, type the following in code/basic_io1.py:

#############################
# FILE INPUT
#############################
# Open a file for reading
f = open('sandbox/test.txt', 'r')
# use "implicit" for loop:
# if the object is a file, python will cycle over lines
for line in f:
    print(line)

# close the file
f.close()

# Same example, skip blank lines
f = open('sandbox/test.txt', 'r')
for line in f:
    if len(line.strip()) > 0:
        print(line)

f.close()

You can try the two blocks separately in IPython before running the complete script. Start IPython at the project root so that sandbox/test.txt names the same file in both cases. Then check the saved script from a shell at the project root:

python3 code/basic_io1.py

The output demonstrates file iteration and blank-line handling; it does not change test.txt.

Note the following:

The for line in f is an implicit loop — implicit because stating the range of things in f to loop over in this way allows Python to handle any kind of objects to loop through.
- For example, if f was an array of numbers 1 to 10, it would loop through them.
- Another example: if f is a file, as in the case of the script above, it will loop through the lines in the file.
if len(line.strip()) > 0 checks if the line is empty. Try ? to see what .strip() does.
There are indentations in the code that determine what is and is not inside the for and if statements. If you get errors or unexpected outputs, it will very likely be because of wrong or missing indentations.

Next, type the following code in code/basic_io2.py. It creates a temporary output in sandbox/.

#############################
# FILE OUTPUT
#############################
# Save the elements of a list to a file
list_to_save = range(100)

f = open('sandbox/testout.txt', 'w')
for i in list_to_save:
    f.write(str(i) + '\n') ## Add a new line at the end

f.close()

Run it from the project root with python3 code/basic_io2.py.

Finally, type the following code in code/basic_io3.py. It saves a Python object in sandbox/, then loads it again.

#############################
# STORING OBJECTS
#############################
# To save an object (even complex) for later use
my_dictionary = {"a key": 10, "another key": 11}

import pickle

f = open('sandbox/testp.p', 'wb')  # b means binary mode
pickle.dump(my_dictionary, f)
f.close()

## Load the data again
f = open('sandbox/testp.p', 'rb')
another_dictionary = pickle.load(f)
f.close()

print(another_dictionary)

Run it from the project root with python3 code/basic_io3.py.

Note

The b flag for reading the file above stands for “binary”. Basically, binary files are machine readable, but not human readable. For example, try opening testp.p in a text reader (e.g., your code editor) and reading it - you will see considerable gibberish (compare with testout.txt).

Safely opening files using `with open()`#

Whilst open() and close() are useful to remember, it can be very problematic if f.close() is missed.

Luckily Python has your back here. By using the with command, you can make sure that no matter what, the file is closed after you have finished working with it.

with is typically used in the following manner:

with open("data/file.txt", "r") as myfile:
    # do things to myfile
    ...
    

Note that running myfile.close() here is not necessary as the file is closed once you drop out of the with block.

Here’s an example of basic_io1.py rewritten using the with statement.

#############################
# FILE INPUT
#############################
# Open a file for reading
with open('sandbox/test.txt', 'r') as f:
    # use "implicit" for loop:
    # if the object is a file, python will cycle over lines
    for line in f:
        print(line)

# Once you drop out of the with, the file is automatically closed

# Same example, skip blank lines
with open('sandbox/test.txt', 'r') as f:
    for line in f:
        if len(line.strip()) > 0:
            print(line)

Note

The rest of this session will use the with method of opening files

Handling CSV files#

The csv package makes it easy to manipulate CSV files. Let’s try it.

$\star$ Download testcsv.csv into your coursework project’s data/ directory. Then type the following as code/basic_csv.py:

import csv

# read a file containing:
# 'Species','Infraorder','Family','Distribution','Body mass male (Kg)'
with open('data/testcsv.csv', 'r', newline='') as f:

    csvread = csv.reader(f)
    temp = []
    for row in csvread:
        temp.append(tuple(row))
        print(row)
        print("The species is", row[0])

# write a file containing only species name and Body mass
with open('data/testcsv.csv', 'r', newline='') as f:
    with open('results/bodymass.csv', 'w', newline='') as g:

        csvread = csv.reader(f)
        csvwrite = csv.writer(g)
        for row in csvread:
            print(row)
            csvwrite.writerow([row[0], row[4]])

Run the script from the project root with python3 code/basic_csv.py. It reads the original table from data/ and creates results/bodymass.csv; it does not add a generated file to the source-data directory. Open the new CSV and check that each row contains a species name and its body mass.

Using column names instead of positions#

For a small table, row[0] and row[4] may be manageable. With a wider table, names such as row["species"] are much easier to follow. csv.DictReader uses the CSV header to turn each row into a dictionary:

with open("data/observations.csv", newline="", encoding="utf-8") as source:
    reader = csv.DictReader(source)
    records = list(reader)

reader.fieldnames contains the header names in their original order. The values themselves are still text, so a number such as 3.5 must be converted before you calculate with it. A missing or extra field can also produce an unusual row, so reading a file successfully does not guarantee that every record is sensible.

The matching csv.DictWriter writes dictionaries back out in a chosen column order:

columns = ["site", "species", "count"]
with open("results/summary.csv", "w", newline="", encoding="utf-8") as destination:
    writer = csv.DictWriter(destination, fieldnames=columns)
    writer.writeheader()
    writer.writerows(records)

Opening the file with "w" replaces its previous contents. writeheader() writes the column names once, and writerows(...) writes the records beneath them.

Practical: programs, modules, and file input/output#

Turn cfexercises1.py into a reusable module and a small reproducible program:

Keep the foo_* function definitions importable, with docstrings describing their inputs and returned values. Do not call input() inside these functions.
Put direct example calls under an if __name__ == "__main__": guard.
Create run_cfexercises.py. Import the functions from cfexercises1, call them with representative inputs, and write labelled results to results/cfexercises_results.txt using with open(...).
Run both files from the coursework project root:

python3 code/cfexercises1.py
python3 code/run_cfexercises.py

You are finished when both commands succeed, importing cfexercises1 does not run its examples, and the results file contains the labelled values you predicted. Keep the two scripts; if you use Git, track them and ignore generated result files, which can be recreated.

Testing, debugging, and errors#

What do you want from your code? rank the following by importance:

it is very fast
it gives me the right answer
it is easy to read
it uses lots of ‘clever’ programming techniques
it uses cool features of the language

Then, think about this:

If you are very lucky, your program will crash when you run it
If you are lucky, you will get an answer that is obviously wrong
If you are unlucky, you won’t notice until after publication
If you are very unlucky, someone else will notice it after publication

Ultimately, most of your time could well be spent finding errors and fixing them (“debugging”), not writing code. You can debug when errors appear, but why not just “nip” as many as you can in the “bud”? For this, you would use unit testing.

Unit testing#

Unit testing prevents the most common mistakes and helps you write reliable code. Indeed, there are many reasons for testing:

Can you prove (to yourself) that your code does what you think it does?
Did you think about the things that might go wrong?
Can you prove to other people that your code works?
Does it still all work if you fix a bug?
Does it still all work if you add a feature?
Does it work with that new dataset?
Does it work on the latest version of the language (e.g., Python 3.x vs. 2.7.x)?
Does it work on Mac? on Linux? on Windows?
Does it work on 64 bit and 32 bit?
Does it work on an old version of a Mac?
Does it work on any Linux server, or Imperial College’s unix cluster?

The idea is to write independent tests for the smallest units of code. Why the smallest units? — to increase the likelihood of retaining the tests upon code modification.

Testing vocabulary for novices

A test case supplies an input and states the expected behaviour. An assertion is an executable check that compares the observed behaviour with that expectation. A unit test checks one small unit of code, usually a function, separately from the larger workflow.

A normal case uses ordinary valid input and checks the main intended behaviour.
Boundary testing, also called Boundary Value Analysis (BVA), checks values at the edges of a rule and, where meaningful, immediately on either side. For a threshold of 10, useful cases may include 9, 10 and 11. Empty input, an all-missing group or the first and last permitted values can also be boundaries. The function’s documented rules determine which boundaries matter.
A failure case uses invalid or unsupported input and checks that the code fails clearly in the expected way, rather than silently returning a misleading result.

Use several kinds of case: one passing normal case does not establish that every boundary or invalid input is handled correctly. See Key computing terms: testing and checking results for a short glossary.

Unit testing with `doctest`#

Let’s start with doctest, a simple testing tool in Python. Its examples sit in a function’s docstring and serve as both documentation and executable checks. Keep the implementation in code/control_flow.py; do not copy it into a test file. Edit even_or_odd so that its accepted inputs, examples and validation are clear:

def even_or_odd(x):
    """Return a description of whether integer x is even or odd.

    >>> even_or_odd(10)
    '10 is Even!'
    >>> even_or_odd(5)
    '5 is Odd!'
    >>> even_or_odd(-2)
    '-2 is Even!'
    """
    if not isinstance(x, int):
        raise TypeError("x must be an integer")
    if x % 2 == 0:
        return f"{x} is Even!"
    return f"{x} is Odd!"

Run the doctests from the project root:

python3 -m doctest -v code/control_flow.py

doctest imports the module before checking its docstrings. The existing if __name__ == "__main__": guard prevents the program’s demonstration output from running during that import. The -v option shows each example, its expected result and the final pass/fail summary.

From doctest to a separate test script#

Doctests work well for small examples that also document a function. As a test suite grows, keep the implementation in one module and import its functions into a separate test script. Download the supplied test_control_flow.py into code/. Its complete listing is shown below; notice that it imports functions rather than copying their definitions:

#!/usr/bin/env python3

"""Checks for the small functions in control_flow.py."""

from control_flow import (
    count_above,
    even_or_odd,
    find_all_primes,
    is_prime,
    largest_divisor_five,
)


# Normal cases: representative even and odd integers.
assert even_or_odd(10) == "10 is Even!"
assert even_or_odd(5) == "5 is Odd!"

# Edge cases: zero and a negative integer remain valid inputs.
assert even_or_odd(0) == "0 is Even!"
assert even_or_odd(-2) == "-2 is Even!"

# Failure case: the function accepts integers, not text.
try:
    even_or_odd("10")
except TypeError:
    pass
else:
    raise AssertionError("A non-integer input was accepted")

# The divisor function returns an integer or None, not display text.
assert largest_divisor_five(120) == 5
assert largest_divisor_five(8) == 4
assert largest_divisor_five(7) is None

# Integers below 2 are not prime; 2 is the first prime.
assert not is_prime(-1)
assert not is_prime(0)
assert not is_prime(1)
assert is_prime(2)
assert is_prime(3)
assert not is_prime(4)
assert find_all_primes(1) == []
assert find_all_primes(10) == [2, 3, 5, 7]

# A normal biological example.
assert count_above([8, 12, 10, 15], 10) == 2

# Values immediately below, at and above the threshold expose > versus >=.
assert count_above([9, 10, 11], 10) == 1

# An empty collection contains no values above the threshold.
assert count_above([], 10) == 0

The first assertions express the same expectations as the doctests, while the remaining checks cover the other functions’ normal, boundary and failure behaviour. In particular, is_prime is checked below, at and above 2. The implementation and checks now have separate responsibilities. try calls even_or_odd with invalid input, except TypeError accepts the promised failure, and else raises AssertionError if the invalid input is wrongly accepted. Other exception types remain visible.

Run the separate test script from the project root:

python3 code/test_control_flow.py

Passing assertions produce no output; the return of the shell prompt means that every check passed.

Examine a meaningful boundary test#

Parity gives useful normal and failure cases, but it has no threshold at which the intended rule changes. The count_above function in control_flow.py provides a simple biological example for genuine boundary testing:

def count_above(heights_cm, threshold_cm):
    """Return the number of heights strictly above a threshold."""
    count = 0
    for height_cm in heights_cm:
        if height_cm > threshold_cm:
            count += 1
    return count

The complete test_control_flow.py listing above imports it and includes these checks:

from control_flow import count_above


# Normal case: some values are above the threshold.
assert count_above([8, 12, 10, 15], 10) == 2

# Boundary case: values immediately below, at and above the threshold.
assert count_above([9, 10, 11], 10) == 1

# Edge case: there are no observations to count.
assert count_above([], 10) == 0

Explain why the value exactly at the threshold is not counted. A test containing only 8 and 12 would not distinguish > from the plausible but incorrect >=; the 10 in the boundary test does. From the project root, rerun python3 code/test_control_flow.py after examining the function and checks.

At this point the same test script checks normal behaviour, an important boundary, empty input and an expected failure. You now have the reusable program pattern used later: an implementation module defines functions, a guarded main() handles the report, and a separate test script imports the functions without triggering that report. The integrated observations practical retains these responsibilities while introducing CSV records, missing values and duplicates.

Other unit testing approaches#

For more complex testing, see the official doctest documentation.

Python’s standard library also provides unittest, while pytest is a widely used external testing framework. The plain assertions here keep the first workflow visible; the same test cases can later be expressed in either framework.

Start testing early, and keep functions small enough that each test has a clear target. Choose cases from the function’s documented behaviour rather than trying arbitrary inputs.

Debugging#

A failing test tells you that observed behaviour differs from the expected behaviour, but it may not reveal why. A debugger lets you pause execution and inspect the values that led to the failure, without scattering temporary print calls throughout the program.

Python’s built-in debugger is pdb. Type the following in a file named debugme.py in your code/ directory:

def buggyfunc(x):
    y = x
    for i in range(x):
        y = y-1
        z = x/y
    return z

buggyfunc(20)

Start IPython from the project root, then run it:

%run code/debugme.py

---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
File ~/Documents/Projects/Software/MulQuaBio/MQB/content/code/debugme.py:8
      5         z = x/y
      6     return z
----> 8 buggyfunc(20)

File ~/Documents/Projects/Software/MulQuaBio/MQB/content/code/debugme.py:5, in buggyfunc(x)
      3 for i in range(x):
      4     y = y-1
----> 5     z = x/y
      6 return z

ZeroDivisionError: division by zero

The traceback points to code/debugme.py and the line where execution failed. Now turn on automatic debugging with %pdb:

%pdb

Now run the script again:

%run code/debugme.py

which gives

---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
code/debugme.py in <module>
      6     return z
      7 
----> 8 buggyfunc(20)

code/debugme.py in buggyfunc(x)
      3     for i in range(x):
      4         y = y-1
----> 5         z = x/y
      6     return z
      7 

ZeroDivisionError: division by zero
> code/debugme.py(5)buggyfunc()
      3     for i in range(x):
      4         y = y-1
----> 5         z = x/y
      6     return z
      7 

Now you will be in the debugger shell, which has a different command prompt: ipdb>.

Here, you can use the following commands to navigate and test the code line by line or block by block:

Command	What it does
`n`	continue execution until the next line
`ENTER`	repeat the previous command
`s`	“step” into next function or procedure (i.e., continue the debugging inside the function, as opposed to simply execute it till it returns)
`p x`	print variable x
`pp locals()`	pretty print all variables and objects in current workspace scope
`c`	continue until next break-point
`l`	print the code surrounding the current position (you can specify how many)
`R`	continue until the end of the function
`q`	quit

So let’s continue our debugging:

ipdb> p x                                                                                                                     
20

ipdb> p y                                                                                                                     
0

ipdb> p z                                                                                                                     
20.0

ipdb> p x/y                                                                                                                   
*** ZeroDivisionError: division by zero

ipdb> l                                                                                                                       
      1 def buggyfunc(x):
      2     y = x
      3     for i in range(x):
      4         y = y-1
----> 5         z = x/y
      6     return z
      7 
      8 buggyfunc(20)
ipdb> q

In []: %pdb
Automatic pdb calling has been turned OFF

Paranoid programming: pause before the failure#

The %pdb example starts the debugger after an exception. You can also stop just before a suspected failure and inspect the values that lead to it. In code/debugme.py, temporarily put a breakpoint between the subtraction and division:

        y = y - 1
        if y == 0:
            breakpoint()
        z = x / y

Run %run code/debugme.py again from IPython. The condition lets earlier loop iterations run; the built-in breakpoint() pauses when y first becomes zero. At the debugger prompt, use p i, p x and p y. What value will the next line try to divide by? Enter q to stop without executing that division. Remove the temporary breakpoint when you finish investigating.

If you have the separate ipdb package installed in the Python environment running the script, you can use its IPython-style debugger at the same point. Replace breakpoint() in the conditional block with:

            import ipdb
            ipdb.set_trace()

The same inspection commands work there. ipdb is an alternative for this exercise, not a package you need to install to use breakpoint() or IPython’s %pdb.

You can inspect the same problem in VS Code: set a breakpoint on the division line and examine the variables before continuing. See the VS Code Python debugging guide. A breakpoint helps explain a failure; a repeatable test checks whether your eventual change gives the intended result.

Handling expected errors#

In many scenarios, you do not want to, or cannot, debug away an error in your program. The two most common such scenarios are:

An error that arises when a particular condition is not met in your program.
You want the code to continue running despite an error.

to “handle” such errors, you can use the try - except keywords.

Let’s expand our previous debugme.py example to understand how this works.

def buggyfunc(x):
    y = x
    for i in range(x):
        try: 
            y = y-1
            z = x/y
        except:
            print(f"This didn't work;{x = }; {y = }")
    return z

buggyfunc(20)

This didn't work;x = 20; y = 0

20.0

$\star$ run this block of code yourself in the ipython console using %cpaste.

So, try does what the name suggests - tries to do something (the division), and if it does not work (any error happens), transfers control to the except block and whatever you ask Python to do in that block gets executed. This is a good strategy if you want your program to give feedback, or that particular error is not common.

You can also “catch” specific types of errors. For example, modify the code further:

def buggyfunc(x):
    y = x
    for i in range(x):
        try: 
            y = y-1
            z = x/y
        except ZeroDivisionError:
            print(f"The result of dividing a number by zero is undefined")
        except:
            print(f"This didn't work;{x = }; {y = }")
        else:
            print(f"OK; {x = }; {y = }, {z = };")
    return z

buggyfunc(20)

OK; x = 20; y = 19, z = 1.0526315789473684;
OK; x = 20; y = 18, z = 1.1111111111111112;
OK; x = 20; y = 17, z = 1.1764705882352942;
OK; x = 20; y = 16, z = 1.25;
OK; x = 20; y = 15, z = 1.3333333333333333;
OK; x = 20; y = 14, z = 1.4285714285714286;
OK; x = 20; y = 13, z = 1.5384615384615385;
OK; x = 20; y = 12, z = 1.6666666666666667;
OK; x = 20; y = 11, z = 1.8181818181818181;
OK; x = 20; y = 10, z = 2.0;
OK; x = 20; y = 9, z = 2.2222222222222223;
OK; x = 20; y = 8, z = 2.5;
OK; x = 20; y = 7, z = 2.857142857142857;
OK; x = 20; y = 6, z = 3.3333333333333335;
OK; x = 20; y = 5, z = 4.0;
OK; x = 20; y = 4, z = 5.0;
OK; x = 20; y = 3, z = 6.666666666666667;
OK; x = 20; y = 2, z = 10.0;
OK; x = 20; y = 1, z = 20.0;
The result of dividing a number by zero is undefined

20.0

So here we anticipated a ZeroDivisionError and gave feedback on that. As such, ZeroDivisionError is a pre-defined error type in python. Another example is TypeError, which you will get if you try to run the above script with something non-numeric, such as a string.

Note that we also used else (just like in the case of if-else) here to give feedback on every successful calculation. There is also a finally keyword that enables you to execute sections of code that should always run, with or without any previously encountered exceptions.

A number that converts may still be unusable#

When numbers come from a data file, successful conversion is only the first check. float("unknown") raises ValueError, which you can catch as above. However, Python also accepts special values such as infinity and ‘not a number’. The math.isfinite function helps us spot them:

import math

print(math.isfinite(float("3.5")))  # True
print(math.isfinite(float("inf")))  # False
print(math.isfinite(float("nan")))  # False

A finite value can still be biologically impossible or outside the range allowed for a particular dataset. That is a separate decision, made using the question and the accompanying metadata.

Tip

Notice again in the above script we create an f string but this time we save ourselves even more time by incorporating the variable and an = sign:

f"Here is a string, {x = }"

This is the same (but even more concise) as writing:

f"Here is a string, x = {x}"

This works for expressions as well! Try something like the following to see:

from math import cos, radians

theta = 30

print(f”{theta = } {cos(radians(theta)) = :.3f}”)

Integrated practical: summarising observations#

This practical brings together functions, control flow, file handling, exceptions and testing. Plan the calculation, implement it in stages, and check the results against small inputs worked out by hand.

Designing an algorithm for imperfect data#

The height example assumed that all measurements were present and in centimetres. Real biological data require additional decisions before we write code.

Example: summarize observations by site#

Row	Site	Date	Count
1	A	2026-10-01	2
2	B	2026-10-01	0
3	A	2026-10-02	missing
4	A	2026-10-01	2

What is the total known count at each site, and how many observations have missing counts? For this example, keep one copy of each exact repeated (site, date, count) record. Count missing observations separately and sum only known counts. These are teaching assumptions, not universal ecological rules: an identical row could be a genuine replicate in another study. A known-count total is a subtotal of recorded counts, not an estimate of what the missing observations would have contributed.

Trace the calculation for site A. Start with a running sum of 0, a missing-count counter of 0 and no known count encountered:

Row	Action for site A	Running sum	Missing count	Any known count encountered?
1	Add the known count 2	2	0	Yes
2	Valid row, but for another site	2	0	Yes
3	Count one missing observation	2	1	Yes
4	Exact duplicate of row 1; skip it	2	1	Yes

Site A therefore has a known-count total of 2 and one missing count. Repeat the trace for B: its observed zero is a known count, so its total is 0, with no missing counts. If every count for a site is missing, report an unavailable total, not zero. We must remember whether any known count was encountered; the running sum alone cannot tell us this.

We can record the overall plan as pseudocode: an informal sketch of an algorithm, not executable Python. There is no required syntax; make the steps and decisions explicit.

INPUT: observations with site, date and count
Validate every record; stop and report an error for invalid input
Identify exact duplicate records and keep one copy of each
For each site:
    count the observations with missing counts
    if at least one count is known:
        add the known counts and report their total
    otherwise:
        report that no known total is available
OUTPUT: known-count total (or unavailable) and missing-count number per site

This plan describes a report for all sites. In the exercise below, first write a reusable function that calculates the result for one requested site. Only after testing it will you write the loop that calls it for each site. Keep the hand-worked results as checks while developing the program.

Minimum requirements for pseudocode

There is no universal pseudocode syntax. For MQB, use these minimum conventions so someone else can follow and check your algorithm:

Inputs and outputs: name what the algorithm receives and what it reports.
Assumptions: state relevant units and rules, such as how to treat duplicates and missing values.
Precise, ordered steps: say what each step does; use consistent names and indentation to show which steps belong together.
Decisions and repetition: give the condition for each choice and explain which items are repeated over, or when repetition stops.
Relevant edge cases: say what happens for invalid or empty input and missing values; distinguish an unavailable result from an observed zero.
A worked check: trace a small input by hand and compare the result with the expected output.

Plain language is enough; Python syntax and capitalised keywords are not required. These conventions make the plan easier to read and check.

Independent exercise: summarise observations by site#

Start with your own plan and attempt, without a coding assistant. If you later use a coding assistant for a hint or explanation, test its suggestions yourself. Use the plan, implement, test, explain and rerun routine.

Your goal: write and test a function that returns the known-count total and missing-count number for one site. Work through the checkpoints in order. Code labelled complete example can be run unchanged; code labelled unfinished template must be completed before its checks can pass. The CSV loader is provided; you write the count converter and site summary.

Reuse the module-and-test-script pattern#

This exercise extends the structure used in From doctest to a separate test script:

Testing example	Observations exercise	Responsibility
`control_flow.py`	`observations_practice.py`	Define the functions.
`test_control_flow.py`	`test_observations_practice.py`	Import the functions and check them.
Normal, boundary and failure cases	Normal, boundary and failure cases	Compare what the function does with its documented behaviour.
From the project root: `python3 code/test_control_flow.py`	From the project root: `python3 code/test_observations_practice.py`	Start a fresh Python process from the same familiar place.

The testing method has not changed. The observations problem adds CSV records, missing values and duplicates. Keep implementation code out of the test script and test code out of the implementation functions.

Inputs and project layout#

Download the ordinary observations, boundary observations and invalid observations. Put them in a project’s data/ directory. Use its code/ directory for the two scripts below; create these directories if needed. Use the filenames shown here.

coursework-project/
    code/
        observations_practice.py
        test_observations_practice.py
    data/
        bootcamp_observations.csv
        bootcamp_observations_boundary.csv
        bootcamp_observations_invalid.csv

Open the coursework project folder in VS Code. As in the earlier exercises, run terminal commands from the project root, where code/ and data/ are siblings. Use pwd before you begin if you are unsure. The examples below use paths such as data/bootcamp_observations.csv, named from the project root.

Field	Teaching input rule
`site`	Non-empty text identifying a site.
`date`	Non-empty text identifying the sampling date. Supplied dates use `YYYY-MM-DD`; calendar validation is outside this exercise.
`count`	Empty text `""` means missing; otherwise text that `int` can convert to a non-negative integer. Zero is a known observation.

Complete example: copy this loader into code/observations_practice.py. It uses Python’s standard-library csv module; no installation or hand-written CSV parser is needed.

import csv


def load_observations(path):
    with open(path, newline="", encoding="utf-8") as stream:
        return list(csv.DictReader(stream))

with open(...) opens the file and closes it when the block ends. csv.DictReader uses the column headings as dictionary keys, and list(...) collects the rows. A blank count stays as the empty string ""; the loader does not convert counts to numbers or validate them.

To inspect the data, add this complete block at the bottom of the loader script:

if __name__ == "__main__":
    rows = load_observations("data/bootcamp_observations.csv")
    print(rows[0])
    print(rows[2])

Run it in a shell terminal, not at Python’s >>> or IPython’s In [...] prompt:

python3 code/observations_practice.py

You should see:

{'site': 'A', 'date': '2026-10-01', 'count': '2'}
{'site': 'A', 'date': '2026-10-02', 'count': ''}

For example, rows[0]["count"] is the text "2", not the integer 2. We will convert it to 2; the empty text will become None, Python’s value for an unavailable count. None, "" and 0 are different values.

The if __name__ == "__main__": guard makes this inspection run only when you execute the script directly. Importing its functions from your checks will not run the inspection. Add new function definitions above this guard, not inside it.

Exercise: convert count text#

Your task: complete parse_count. Add the unfinished template below the loader function, above the startup guard. Keep the function definition and the existing missing-value check. Delete the line beginning raise NotImplementedError and replace it with your own conversion, validation and return statements, indented inside the function. Do not add input(); the argument count_text receives its value when the function is called.

def parse_count(count_text):
    if count_text == "":
        return None
    raise NotImplementedError("Replace this line with conversion and validation")

The placeholder deliberately stops execution for every non-empty input. It marks unfinished code, not invalid data. Your finished function must have these behaviours:

Call	Required result
`parse_count("2")`	Return the integer `2`.
`parse_count("0")`	Return the integer `0`.
`parse_count("")`	Return `None`.
`parse_count("-1")`	Raise `ValueError`.
`parse_count("many")`	Raise `ValueError`.
`parse_count(2)`	Raise `ValueError`: the input must be text.

Use isinstance(count_text, str) to check whether the input is text and int(count_text) to convert non-empty text. int("many") already raises ValueError; int("-1") succeeds, so you must reject a negative result. Return the converted integer when it is non-negative.

This is the same failure-testing pattern used for the invalid text passed to even_or_odd. Here the promised exception is ValueError: raise stops the current call, except ValueError accepts that expected failure, and the else block makes the check fail if invalid data are accepted. Do not catch every exception; a different error should remain visible.

Create code/test_observations_practice.py and put these checks in it. This is test code for your completed function, not the missing implementation:

from observations_practice import load_observations, parse_count

# Normal case: ordinary count text becomes an integer.
assert parse_count("2") == 2

# Boundary case: zero is a valid known count.
assert parse_count("0") == 0

# Special valid case: blank text represents a missing count.
assert parse_count("") is None

# Failure cases: each invalid value must raise ValueError.
for bad_count in ["-1", "many", 2]:
    try:
        parse_count(bad_count)
    except ValueError:
        pass
    else:
        raise AssertionError(f"Invalid count accepted: {bad_count!r}")

pass means no further action is needed for the expected exception; !r shows the offending value with its representation, including quotes for text. Run the following shell command from the project root after completing the template:

python3 code/test_observations_practice.py

Checkpoint: all count-conversion checks pass, no placeholder remains, and you can explain why missing and zero have different results. Only then continue.

Now implement the site summary#

In the same implementation script, above the startup guard, write summarise_site(rows, site). rows is a list of dictionaries returned by the loader; site is the name of one requested site. Return exactly the two-item tuple (known_total, missing_count), without changing the input rows. For A in the ordinary file this is (2, 1). Do not write the all-site reporting loop yet.

First write pseudocode for this one-site function and trace the A example by hand. Develop one function in four stages, extending it and rerunning the earlier checks each time. The early stages are temporary versions, not the finished specification.

Stage 1: known counts in valid, unique rows. Select rows whose row["site"] matches the requested site, call parse_count(row["count"]), and add their known counts. Return the sum and a missing count of 0. At this stage, use only valid rows with non-missing counts and no duplicates.

Add summarise_site to the import in your test script and append this check:

from observations_practice import summarise_site

unique_rows = [
    {"site": "A", "date": "2026-10-01", "count": "2"},
    {"site": "A", "date": "2026-10-02", "count": "3"},
    {"site": "B", "date": "2026-10-01", "count": "9"},
]
assert summarise_site(unique_rows, "A") == (5, 0)

Stage 2: missing counts and observed zeros. Count missing observations separately. Track whether any known count has been encountered: a Boolean such as has_known_count, initially False, can record this. If none was encountered, return None for the total. Use count is None to identify missing values, not if count: zero is known even though it is false in a Boolean test. An empty list or absent site must return (None, 0).

missing_rows = [
    {"site": "C", "date": "2026-10-01", "count": ""},
    {"site": "C", "date": "2026-10-02", "count": ""},
]
zero_rows = [{"site": "D", "date": "2026-10-01", "count": "0"}]
assert summarise_site(missing_rows, "C") == (None, 2)
assert summarise_site(zero_rows, "D") == (0, 0)
assert summarise_site([], "A") == (None, 0)
assert summarise_site(unique_rows, "Z") == (None, 0)

Stage 3: exact duplicate records. Keep one copy of each exact (site, date, count) record, including missing-count records. Different dates are distinct. A set of record tuples can remember which records have already been encountered; do not remove or modify rows in the input list. Test a repeated non-zero count and a repeated missing count: repeating zero alone cannot expose double-counting.

assert summarise_site(unique_rows + [unique_rows[0]], "A") == (5, 0)
assert summarise_site(missing_rows + [missing_rows[0]], "C") == (None, 2)

Stage 4: validation before filtering or deduplication. The finished function must require a list of row dictionaries and a non-empty text site argument. Validate every row, including rows for other sites: exactly the fields site, date and count, all with text values; non-empty site/date; and a count accepted by parse_count. Raise ValueError for invalid input instead of silently skipping it. No date parsing or whitespace normalisation is required.

For field validation, set(row) gives a dictionary’s keys, which you can compare with {"site", "date", "count"}. Check that a row is a dictionary before reading its fields, and that the fields are present before accessing them. Reuse parse_count rather than writing the conversion again.

Check, explain and rerun#

After all four stages, append the file-based checks to your test script. These call your implementation; do not copy its body into the test file.

rows = load_observations("data/bootcamp_observations.csv")
assert summarise_site(rows, "A") == (2, 1)
assert summarise_site(rows, "B") == (0, 0)

boundary_rows = load_observations("data/bootcamp_observations_boundary.csv")
assert summarise_site(boundary_rows, "C") == (None, 2)
assert summarise_site(boundary_rows, "D") == (0, 0)
assert summarise_site(rows, "Z") == (None, 0)

A checks duplicate removal and missing counts; B and D check observed zeros; C checks an unavailable total. Retain the earlier non-zero and missing-duplicate checks, since D’s repeated zero alone cannot prove duplicate removal.

Complete example: add this helper and failure checks to the test script. Each call must raise ValueError; any other exception remains visible. The invalid file has two rows, which must be checked separately. Requesting A must still reject the invalid B row.

def expect_summary_error(rows, site="A"):
    try:
        summarise_site(rows, site)
    except ValueError:
        pass
    else:
        raise AssertionError("Expected ValueError, but the call succeeded")

invalid_rows = load_observations("data/bootcamp_observations_invalid.csv")
for invalid_row in invalid_rows:
    expect_summary_error([invalid_row])

expect_summary_error([{"site": "A", "count": "2"}])
expect_summary_error([{"site": "A", "date": "2026-10-01", "count": "2", "extra": "x"}])
expect_summary_error([{"site": "", "date": "2026-10-01", "count": "2"}])
expect_summary_error([{"site": "A", "date": "", "count": "2"}])
expect_summary_error([{"site": "A", "date": "2026-10-01", "count": 2}])
expect_summary_error(["not a row dictionary"])
expect_summary_error(None)
expect_summary_error(rows, "")
expect_summary_error(rows, 2)

Also confirm the input rows are unchanged. Copy the dictionaries, not just the outer list, so a change to a field is detectable:

original_rows = [row.copy() for row in rows]
summarise_site(rows, "A")
assert rows == original_rows

Run all checks from the project root in a fresh Python process:

python3 code/test_observations_practice.py

No output and a return to the shell prompt means all these checks passed. AssertionError means a check failed; NotImplementedError means a template remains unfinished. ModuleNotFoundError or FileNotFoundError usually means you should check filenames, file locations and the working directory. Do not type shell commands at a Python prompt.

Finish the algorithm exercise: make a new valid input with a different site and dates, predict its result by hand, then add a check. Explain why (None, 2) differs from (0, 0), how you recognise duplicates, and why an invalid row for another site still raises an error. You are ready for the prepared-suggestion critique when your normal, boundary, failure and unchanged-input checks pass and you can explain these decisions.

Follow-on: report and reproduce#

After the function checks pass, replace the inspection block under if __name__ == "__main__": in your implementation script with a report that loads the ordinary file and calls your tested function once for each distinct site. For example, sorted({row["site"] for row in rows}) gives the distinct site names in a consistent order. The expected summaries are A: (2, 1) and B: (0, 0); label each site’s output.

Run python3 code/observations_practice.py from the project root. Compare a peer’s fresh run with your prediction, explain one assumption and one failure check, and record a meaningful Git checkpoint containing the two scripts and three input files. To check reproducibility, use a separate clean checkout containing those files and run both the test command and report command from its root. Do not reset or discard your working copy.

Can your tests detect a plausible wrong answer?#

Complete the independent observations exercise, including its checks and your new predicted input, before this critique.

Before continuing, confirm that your test script now checks ordinary valid input, an important boundary, invalid input that must raise an exception, and empty, missing and duplicated data. You should be able to explain the expected result of every check independently of the implementation. The next task asks whether those checks can distinguish a correct result from a plausible but incorrect one; do not replace or weaken your own implementation.

The function below was written by the teaching team and deliberately contains a fault. It looks reasonable and passes the ordinary checks, but a boundary test can expose its incorrect assumption. The purpose of this exercise is to practise checking unfamiliar code, regardless of where it came from. It is an exercise, not a reference solution.

The prepared function returns only the known-count total for one site. Use valid rows and focus on the missing-value rule, not on adding validation or the separate missing-count output. Add this function above the startup guard in your implementation script so you can compare its behaviour with your own function.

Prepared faulty example: do not copy this function into your solution as if it were correct. First predict its behaviour, then use tests to identify the case it handles incorrectly.

def suggested_total(rows, site):
    total = 0
    seen = set()
    for row in rows:
        record = (row["site"], row["date"], row["count"])
        if record in seen:
            continue
        seen.add(record)
        if row["site"] == site and row["count"] != "":
            total += int(row["count"])
    return total

In your test script, import suggested_total alongside your other functions. With the ordinary rows loaded, these checks pass:

from observations_practice import suggested_total

assert suggested_total(rows, "A") == 2
assert suggested_total(rows, "B") == 0

These checks pass because the ordinary data do not distinguish an unavailable total from a genuine total of zero. Passing them therefore does not establish that the function implements the full data rule. Before running the prepared function on boundary_rows, predict its result for C and compare it with the rule: no known count means an unavailable total (None), not zero.

Write a check that passes for the known-total part of your own correct summary, summarise_site(boundary_rows, "C")[0], but fails for suggested_total(boundary_rows, "C"). The [0] selects the total from your two-item tuple.
Explain why the ordinary input did not expose the fault. Also check D: a genuinely observed zero must remain zero.
Record the deliberate failure separately, then remove or comment out that failing assertion before rerunning your own normal, boundary and failure checks in a fresh process. Do not weaken the expected result to make the suggestion pass.

The important question is not where unfamiliar code came from. Check whether its behaviour matches the stated problem and whether your tests cover the important cases. Retain, revise or reject it on that evidence, not because it looks convincing or because a few tests pass. The same process applies to code suggested by another person, found online or produced by a coding assistant. Explain the missing-versus-zero distinction without assistance.

Capstone practical: testing biological programs#

This capstone contains two main practicals: sequence alignment and tree-survey validation. The two practicals can be worked through independently. The alignment extensions develop the same techniques further.

For each practical, work out the expected result for a small example before running the code. Keep reusable functions in the implementation script and import them into a separate test script. This makes it easier to organise several inputs and compare the program’s results with expectations you have worked out for yourself.

Align DNA sequences#

Download align_seqs.py into your project’s code/ directory. This script aligns two DNA sequences using the simple offset-and-score algorithm described below. Run it and make sure you understand what every line does. If you want to inspect the program while it runs, insert Python’s built-in breakpoint() just after the start of the for loop in calculate_score; use n, p, l, c and q as in the debugging section. This is optional, and no additional debugger package is required.

The aligning algorithm is simple. Start by positioning the beginning of the shorter sequence at all positions (bases) of the longer one (the start position), and count the number of bases matched downstream. Then, for each start position, count the “score” as total of number of bases matched. The alignment with the highest score wins. Ties are possible, in which case, you just take the arbitrary alignment (e.g., first or last) with the highest score.

Your tasks:

Convert align_seqs.py into a non-interactive program that reads two DNA sequences from one input file and saves the best alignment and score as labelled plain text in results/. Do not use input(); running the program should require only this command from the project root:

python3 code/align_seqs.py

A two-row CSV containing the example sequences is sufficient. Validate that exactly two non-empty sequences were read. Put program startup under an if __name__ == "__main__": guard so that test_align_seqs.py can import functions without running the whole program. In the test script, include a small normal case with a hand-calculated score, a boundary or edge case such as an alignment reaching the end of a sequence, and failure cases for missing, empty or surplus input sequences. From the project root, run python3 code/test_align_seqs.py. You are finished when those checks pass and the saved alignment and score match the original example’s result.

Extension: align FASTA sequences#

Align all the .fasta sequences from the unix chapter. Call the new script align_seqs_fasta.py. Unlike align_seqs.py, this script should take any two fasta sequences (in separate files) to be aligned as input. So this script would typically run by using explicit inputs, for example, when called with something like:

python3 code/align_seqs_fasta.py data/seq1.fasta data/seq2.fasta

However, it should still run if no inputs were given, using two fasta sequences from the data directory as defaults.

Extension: retain tied best alignments#

The current script/program runs through all possible starting points on the main sequence and then just takes the first of the alignments with the highest score. This should be apparent if you closely examine this part of the script:

for i in range(l1):
    z = calculate_score(s1, s2, l1, l2, i)
    if z > my_best_score:
        my_best_align = "." * i + s2
        my_best_score = z

The if z > my_best_score condition means that when a new alignment has the same score as the current one, the current one will be kept and the new one discarded. As a result, multiple alignments with the same score will be lost even if they are all equally good alignments (highly likely in longer sequences), with only the first one retained.

Modify the script so that all equally best alignments are recorded. Save them with their shared score in a human-readable CSV or TSV file in results/; do not use pickle for portable results. Call the new script align_seqs_better.py. Add a small input with a known tie and verify that every tied alignment is retained.

Practical: validating tree-survey records#

This practical brings together functions, data structures, control flow, CSV handling and testing in one tree-record cleaning workflow. The survey is made up to illustrate a realistic mixture of spelling mistakes, missing entries, invalid measurements and duplicate identifiers that you will likely encounter in real data recorded by real people!

Use the linked metadata whenever you need the exact validation rules or output format.

Download the files into the root-level data/ and code/ directories of the project you use for MQB practicals:

Tree-survey data to data/.
Data metadata and validation rules to data/.
Starter Python program to code/.

Complete the cleaner#

Edit your downloaded copy of code/clean_silwood_trees.py and keep that filename. It is an unfinished starter. Complete its three unfinished functions: validate_record, process_records and summarise_oaks. Keep your completed cleaner alongside a separate test script.

Keep the original survey data, data/silwood_trees_raw.csv, unchanged. The instruction not to edit the data does not apply to the Python starter! Use the metadata rather than guessing how a doubtful record should be corrected.

Keep reusable function definitions outside main() so that your test script can import them. The supplied main() calls the functions that read the survey and write the results; recall that the if __name__ == "__main__": guard calls main() only when you run the script directly. Reading and writing helpers can also be defined outside main()—defining a function does not run it! Leave the supplied startup code in place: importing the module should not read the survey or create any result files.

See how the pieces fit together#

Before writing any code, read through the function names and their docstrings. You do not need to run the program or use a debugger at this point. The supplied code takes care of reading records and writing results; your three unfinished functions make the decisions in between.

load_records reads the CSV header and returns one dictionary per row. Its default is the supplied survey, but tests can give it a different file. Each row then reaches validate_record, together with its original row number and the identifiers accepted so far. This function returns either a cleaned record or None, plus any reasons for rejecting the row.

process_records will repeat that step for the whole survey, gathering the accepted records in one list and the rejected rows in another. summarise_oaks will work only with the records that made it through, producing the species counts for the third output. Finally, main joins these stages together and passes each result to the supplied write_csv helper. The main guard means this sequence begins when you run the script, but not when your test script imports its functions.

The three NotImplementedError lines show you exactly where your work goes. Until you replace them, running the starter will stop at the first unfinished function; that is expected. Leave the loader, writer and startup code in place. Before you begin, make sure you can say what each of your three functions receives and returns, then sketch the same flow as pseudocode in your README.

The ideas used in the starter have already appeared in the chapter: type hints, project paths, CSV rows as dictionaries, checking converted numbers, and the main guard. Return to any of these as you need them rather than trying to absorb the whole file at once.

Write and run your tests#

Create a separate file, code/test_silwood_trees.py, alongside the cleaner. Import the cleaner’s functions into this file and check them using small records made for the tests, rather than changing the survey itself. The three unfinished functions belong in the cleaner, not in this test script.

Include ordinary text cleaning, missing and unrecognised values, valid and invalid measurements, measurements outside the allowed range, duplicate identifiers and rows with more than one problem. For the oak summary, check both an exact oak name and a near match, as well as the case where no oaks are accepted. Finally, check that importing the module does not start the program and that running it twice gives the same results. You have already met the testing needed here; the debugging and expected-error sections are there if you want to revisit them.

From the project root, run python3 code/test_silwood_trees.py to run your checks. This tests your implementation; running the cleaner itself produces the three survey output files described below.

What your program should produce#

From the project root, run python3 code/clean_silwood_trees.py. It should create (or replace) three files: results/silwood_trees_clean.csv for accepted records, results/silwood_tree_issues.csv for rejected rows and their reasons, and results/oak_species_counts.csv for the oak summary. Your program should create results/ if it is missing. Check the metadata for the precise columns, order and wording. These files can always be regenerated, so you would not normally commit them.

Keep a record of your approach#

Before you implement the functions, add a short plan to the root README.md. Name the input and outputs, then use Markdown pseudocode to show how rows become accepted records, issue reports and the oak summary. Designing an algorithm for imperfect data has a worked example. There is no need to copy the detailed rules from the metadata.

Also record the commands for running and testing the program, what each output contains, one useful test and the fault it would reveal, any outside help you used, and any limitation you know about. Keep this concise: it should help you think through the program and help somebody else run it.

Use the familiar MQB project layout: README.md and .gitignore at the root, scripts in code/, and the unchanged input and metadata in data/. Keep your modified, completed code/clean_silwood_trees.py and your own code/test_silwood_trees.py, together with the unchanged input and metadata and your root README.md. If you use Git, track those source and input files and use .gitignore for generated files. Keep quick experiments in sandbox/, and leave generated results, caches and temporary files out of Git. When you think you are finished, run the tests from the project root, run the cleaner twice, inspect all three outputs, check that the input has not changed, and try the commands exactly as you wrote them in the README.

Before starting, revisit From a task to an algorithm, Structuring Python programs, Handling CSV files, and the separate test-script example.

Before you finish#

Test the normal cases and the relevant edge or invalid cases named in the practical. Work out small expected results before running the code.
Keep reusable functions separate from code that starts the program, so they can be imported and tested safely.
Add short docstrings that explain what each module and function does. Comment decisions that are not obvious from the code itself.
From a new terminal at the project root, rerun the documented commands, inspect the outputs and check that running the program again gives the intended result.
If you use Git, check git status, track source code and intentional input files, and leave generated results, caches and temporary experiments untracked.
If possible, ask someone else to try the documented commands from a clean checkout. This is a useful final check, not a condition for completing the practical.

XKCD on Python — Fig. 12 **Is Python the most common answer to your daily programming needs?** Possibly!
(Source: XKCD)#

Numerical computing in Python#

Python is a good choice for numerical computing (recall this comparison). Using the right packages, you can do some serious number crunching in python.

We will use the numpy and scipy packages. The latter offers a data structure called numpy array that is suitable for efficient computing, along with basic operations on these arrays.

On the other hand, scipy is proper numerical computing package (which uses numpy arrays) that can do a lot of things, including:

Linear algebra (matrix and vector operations) using scipy.linalg
Dealing with sparse matrix problems using scipy.sparse
Numerical integration (including solving of Ordinary Differential Equations (ODEs)) using scipy.integrate
random number generation and using statistical functions and transformations using scipy.stats
Optimization using scipy.optimize
Fourier transforms using scipy.fft
Signal Processing using scipy.signal

We will learn more about scipy further below. First, let’s start with numpy.

Numpy#

Numpy provides the numpy array data structure, along with methods for for data creation, manipulations and basic numerical calculations. Numpy arrays are similar in some respects to Python lists, but are homogeneous in type (the default is float), allow efficient (fast) manipulations, and are more naturally multidimensional (e.g., you can store multiple matrices in one array).

Note

numpy arrays are analogous to the r matrix data structure.

Let’s start by importing numpy:

import numpy as np

a = np.array(range(5)) # a one-dimensional array
a

array([0, 1, 2, 3, 4])

print(type(a))

<class 'numpy.ndarray'>

print(type(a[0]))

<class 'numpy.int64'>

The last two outputs show that a is a numpy.ndarray and that a[0] is a NumPy integer scalar (often numpy.int64 on a 64-bit system), not a built-in Python int. The array’s dtype describes the type used for its elements. range() yields Python integers, but NumPy converts them to an array dtype when it creates a.

Python numpy array — Fig. 13 **A graphical depiction of NumPy arrays.** These can have multiple dimensions (even greater than 3). Source: Cornell Computational Methods#

You can also specify the data type of the array:

a = np.array(range(5), float)
a

array([0., 1., 2., 3., 4.])

a.dtype # Check type 

dtype('float64')

You can also get a 1-D arrays as follows:

x = np.arange(5)
x

array([0, 1, 2, 3, 4])

x = np.arange(5.) #directly specify float using decimal
x

array([0., 1., 2., 3., 4.])

As with other Python variables (e.g., created as a list or a dictionary), you can apply methods to variables created as numpy arrays. For example, type x. and hit TAB to see all methods you can apply tox. To see dimensions of x:

x.shape

(5,)

Tip

remember, you can type :?x.methodname to get info on a particular method (command). For example, try ?x.shape.

You can also convert to and from Python lists (recall list comprehensions from the Python I chapter):

b = np.array([i for i in range(10) if i % 2 == 1]) #odd numbers between 1 and 10 
b

array([1, 3, 5, 7, 9])

c = b.tolist() #convert back to list
c

[1, 3, 5, 7, 9]

To make a matrix, you need a 2-D numpy array:

mat = np.array([[0, 1], [2, 3]])
mat

array([[0, 1],
       [2, 3]])

mat.shape

(2, 2)

Indexing and accessing arrays#

As with other Python data objects such as lists, numpy array elements can be accessed using square brackets ([ ]) with the usual [row,column] reference. Indexing of numpy arrays works like that for other data structures, with index values starting at 0. So, you can obtain all the elements of a particular row as:

mat[1] # accessing whole 2nd row, remember indexing starts at  0

array([2, 3])

mat[:,1] #accessing whole second column  

array([1, 3])

And accessing particular elements:

mat[0,0] # 1st row, 1st column element

np.int64(0)

mat[1,0] # 2nd row, 1st column element

np.int64(2)

Note that (like all other programming languages) row index always comes before column index. That is, mat[1] is always going to mean “whole second row”, and mat[1,1] means 1st row and 1st column element. Therefore, to access the whole second column, you need:

mat[:,0] #accessing whole first column  

array([0, 2])

Python indexing also accepts negative values for going back to the start from the end of an array:

mat[0,1]

np.int64(1)

mat[0,-1]

np.int64(1)

mat[-1,0]

np.int64(2)

mat[0,-2]

np.int64(0)

Very interesting, but rather useless for this simple matrix!

Manipulating arrays#

Manipulating numpy arrays is pretty straightforward.

Note

Why NumPy arrays are computationally efficient: The data associated with a NumPy array object — its metadata (number of dimensions, shape, data type, etc.) and its values — are stored in a homogeneous, contiguous block of Random Access Memory (RAM) called a data buffer. This can make NumPy arrays more efficient than Python data structures such as lists, whose contents need not be stored contiguously.

Replacing, adding or deleting elements#

Let’s look at how you can replace, add, or delete an array element (a single entry, or whole row(s) or whole column(s)):

mat[0,0] = -1 #replace a single element
mat

array([[-1,  1],
       [ 2,  3]])

mat[:,0] = [12,12] #replace whole column
mat

array([[12,  1],
       [12,  3]])

np.append(mat, [[12,12]], axis = 0) #append row, note axis specification

array([[12,  1],
       [12,  3],
       [12, 12]])

np.append(mat, [[12],[12]], axis = 1) #append column

array([[12,  1, 12],
       [12,  3, 12]])

newrow = [[12,12]] #create new row

mat = np.append(mat, newrow, axis = 0) #append that existing row
mat

array([[12,  1],
       [12,  3],
       [12, 12]])

np.delete(mat, 2, 0) #Delete 3rd row

array([[12,  1],
       [12,  3]])

And concatenation:

mat = np.array([[0, 1], [2, 3]])
mat0 = np.array([[0, 10], [-1, 3]])
np.concatenate((mat, mat0), axis = 0)

array([[ 0,  1],
       [ 2,  3],
       [ 0, 10],
       [-1,  3]])

Flattening or reshaping arrays#

You can also “flatten” or “melt” arrays, that is, change array dimensions (e.g., from a matrix to a vector):

mat.ravel()

array([0, 1, 2, 3])

mat.reshape((4,1))

array([[0],
       [1],
       [2],
       [3]])

This is different from ravel: check the documentation by using ?np.reshape.

mat.reshape((1,4))

array([[0, 1, 2, 3]])

Note

You might have noticed that flattening and reshaping is “row-priority”: elements of the vector are allocated to a matrix row-wise, and vice versa(e.g., with ravel unraveling also happens row by row).

mat.reshape((3, 1))

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[215], line 1
----> 1 mat.reshape((3, 1))

ValueError: cannot reshape array of size 4 into shape (3,1)

This gives an error because total elements must remain the same!

This is a bit different than how R behaves, where you won’t get an error (r “recycles” data) (which can be dangerous!)

Pre-allocating arrays#

As in other computer languages, it is usually more efficient to preallocate an array than to append, insert or concatenate additional elements, rows or columns. Resizing may require copying the array to a larger contiguous block of RAM, which can be slow for a large array. Preallocation requests the required memory once.

For example, if you know the size of your matrix or array, you can initialize it with ones or zeros:

np.ones((4,2))

array([[1., 1.],
       [1., 1.],
       [1., 1.],
       [1., 1.]])

Note that (4,2) are the (row,col) array dimensions.

np.zeros((4,2)) # or zeros

array([[0., 0.],
       [0., 0.],
       [0., 0.],
       [0., 0.]])

m = np.identity(4) #create an identity matrix
m

array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]])

m.fill(16) #fill the matrix with 16
m

array([[16., 16., 16., 16.],
       [16., 16., 16., 16.],
       [16., 16., 16., 16.],
       [16., 16., 16., 16.]])

`numpy` matrices#

Scipy/Numpy also has a matrix data structure class. Numpy matrices are strictly 2-dimensional, while numpy arrays are N-dimensional. Matrix objects are a subclass of numpy arrays, so they inherit all the attributes and methods of numpy arrays (ndarrays).

The main advantage of scipy matrices is that they provide a convenient notation for matrix multiplication: for example, if a and b are matrices, then a * b is their matrix product.

Matrix-vector operations#

Now let’s perform some common matrix-vector operations on arrays (you can also try the same using matrices instead of arrays):

mm = np.arange(16)
mm = mm.reshape(4,4) #Convert to matrix
mm

array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11],
       [12, 13, 14, 15]])

mm.transpose()

array([[ 0,  4,  8, 12],
       [ 1,  5,  9, 13],
       [ 2,  6, 10, 14],
       [ 3,  7, 11, 15]])

mm + mm.transpose()

array([[ 0,  5, 10, 15],
       [ 5, 10, 15, 20],
       [10, 15, 20, 25],
       [15, 20, 25, 30]])

mm - mm.transpose()

array([[ 0, -3, -6, -9],
       [ 3,  0, -3, -6],
       [ 6,  3,  0, -3],
       [ 9,  6,  3,  0]])

mm * mm.transpose() # Note that this is element-wise multiplication

array([[  0,   4,  16,  36],
       [  4,  25,  54,  91],
       [ 16,  54, 100, 154],
       [ 36,  91, 154, 225]])

mm // mm.transpose()

/tmp/ipykernel_317113/1630225829.py:1: RuntimeWarning: divide by zero encountered in floor_divide
  mm // mm.transpose()

array([[0, 0, 0, 0],
       [4, 1, 0, 0],
       [4, 1, 1, 0],
       [4, 1, 1, 1]])

Note that we used integer division //. Note also the warning you get (because of zero division). So let’s avoid the divide by zero:

mm // (mm + 1).transpose()

array([[0, 0, 0, 0],
       [2, 0, 0, 0],
       [2, 1, 0, 0],
       [3, 1, 1, 0]])

mm * np.pi

array([[ 0.        ,  3.14159265,  6.28318531,  9.42477796],
       [12.56637061, 15.70796327, 18.84955592, 21.99114858],
       [25.13274123, 28.27433388, 31.41592654, 34.55751919],
       [37.69911184, 40.8407045 , 43.98229715, 47.1238898 ]])

mm.dot(mm) # No this is matric multiplication, or the dot product

array([[ 56,  62,  68,  74],
       [152, 174, 196, 218],
       [248, 286, 324, 362],
       [344, 398, 452, 506]])

There is also a numpy matrix class:

mm = np.matrix(mm) # convert to scipy/numpy matrix class
mm

matrix([[ 0,  1,  2,  3],
        [ 4,  5,  6,  7],
        [ 8,  9, 10, 11],
        [12, 13, 14, 15]])

print(type(mm))

<class 'numpy.matrix'>

This data structure makes matrix multiplication syntactically easier:

mm * mm # instead of mm.dot(mm)

matrix([[ 56,  62,  68,  74],
        [152, 174, 196, 218],
        [248, 286, 324, 362],
        [344, 398, 452, 506]])

Warning

However, it is not recommended that you use the numpy matrix class because it may be removed in the future.

Tip

You can do a lot more with matrices and vectors by importing the linalg sub-package from scipy: scipy.linalg.

The `scipy` package#

Now let’s move on to scipy.

Note

scipy vs numpy: It’s a bit confusing. Please have a look at this. Basically, there is some overlap between what these two packages can do, and this redundancy will eventually be phased out completely. The recommended approach is to use numpy for creating and manipulating data and scipy for more complex numerical operations.

We will look at two particularly useful scipy sub-packages here: scipy.stats and scipy.integrate.

Scipy stats#

Let’s take a quick spin in scipy.stats.

Why not use R for stats? — because often you might just need to calculate some summary stats of your simulation results within Python, or you simulations may require the generation of random numbers.

First, import scipy:

import scipy as sc

Or you can use from scipy import stats.

Let’s generate 10 samples from the normal distribution ($\mathcal{N}(\mu,\,\sigma^{2})$):

sc.stats.norm.rvs(size = 10)

array([ 0.03379392, -1.10285228,  1.33339903,  1.29748676, -1.35256438,
       -0.52885516,  1.15613648, -0.25302346,  0.16701658, -1.49492858])

By default, as in R, these are numbers from the standard normal distribution ($\mathcal{N}(0,\,1)$).

Tip

Continuing on the numpy vs scipy theme: Both scipy and numpy can generate random numbers (e.g., np.random.normal(size=10) would work equally well above). In fact, scipy uses the numpy.random package under the hood for drawing random numbers. You may choose to use either, but for sanity, its probably a good idea to just stick with scipy.stats for all your stats number crunching.

Also, as you learned in R, you can “seed” random numbers to get the same sequence every time (important for reproducibility – when you need to know what specific random numbers were input into a particular program routine or algorithm).

np.random.seed(1234)
sc.stats.norm.rvs(size = 10)

array([ 0.47143516, -1.19097569,  1.43270697, -0.3126519 , -0.72058873,
        0.88716294,  0.85958841, -0.6365235 ,  0.01569637, -2.24268495])

But setting a global random number state is not always needed or in some many cases, recommended, because using something like np.random.seed(1234) will set the seed for all following instances where a random number is generated. In many scenarios a more robust way is to use the random_state argument for each specific generation of a set of random numbers:

sc.stats.norm.rvs(size=5, random_state=1234)

array([ 0.47143516, -1.19097569,  1.43270697, -0.3126519 , -0.72058873])

We will not move on from scipy.stats. Before we do so, here’s an example of generating random integers between 0 and 10:

sc.stats.randint.rvs(0, 10, size = 7)

array([5, 0, 9, 6, 2, 0, 5])

And again with a random seed:

sc.stats.randint.rvs(0, 10, size = 7, random_state=1234)

array([3, 6, 5, 4, 8, 9, 1])

sc.stats.randint.rvs(0, 10, size = 7, random_state=3445) # a different seed

array([8, 1, 9, 6, 8, 6, 0])

Numerical integration using `scipy`#

OK, on to to and scipy.integrate.

Numerical integration is the approximate computation of an integral using numerical techniques. You need numerical integration whenever you have a complicated function that cannot be integrated analytically using anti-derivatives. For example, calculating the area under a curve is a particularly useful application. Another one is solving ordinary differential equations (ODEs), commonly used for modelling biological systems.

import scipy.integrate as integrate

Area under a curve#

Let’s calculate the area under an arbitrary curve.

y = np.array([5, 20, 18, 19, 18, 7, 4]) # The y values; can also use a python list here

Let’s visualize the curve. We can use the matplotlib package for this:

import matplotlib.pylab as p

p.plot(y)

[<matplotlib.lines.Line2D at 0x7433aefaa270>]

../_images/2f28a941039d4ed271b884880c555c5430e6a0a480d0859021e8190f0133d4dd.png

Now compute the area using the composite trapezoidal rule:

area = integrate.trapezoid(y, dx = 2)
print("area =", area)

area = 173.0

The argument dx defines the spacing between points of the curve (the x-axis values). The default is 1 (don’t forget to check out the documentation: ?integrate.trapezoid). Changing this will change the area, of course:

area = integrate.trapezoid(y, dx = 1)
print("area =", area)

area = 86.5

area = integrate.trapezoid(y, dx = 3)
print("area =", area)

area = 259.5

Now, the same, using Simpson’s rule:

area = integrate.simpson(y, dx = 2)
print("area =", area)

area = 176.66666666666666

area = integrate.simpson(y, dx = 1)
print("area =", area)

area = 88.33333333333333

area = integrate.simpson(y, dx = 3)
print("area =", area)

area = 265.0

Note

There are alternatives to matplotlib for visualizations in python. Check out bokeh, plotly and altair. Also, there is seaborn, which is built directly on Matplotlib. Examples using seaborn can be found in some of the appendices of these notes.

The Lotka-Volterra model#

Now let’s try numerical integration in Python for solving a classical model in biology — the Lotka-Volterra (LV) model for a predator-prey system in two-dimensional space (e.g., on land). The LV model is:

(2)#\[\begin{align} \frac{dR}{dt} &= r R - a C R \nonumber \\ \frac{dC}{dt} &= - z C + e a C R \end{align}\]

where,

$C$ and $R$ are the consumer (e.g., predator) and resource (e.g., prey) population abundances, measured here as number $\times$ area$^{-1}$.
$r$ is the intrinsic (per-capita) growth rate of the resource population (time$^{-1}$).
$a$ is per-capita “search rate” for the resource ($\text{area}\times \text{time}^{-1}$) multiplied by its attack success probability, which determines the encounter and consumption rate of the consumer on the resource.
$z$ is mortality rate ($\text{time}^{-1}$) and $e$ is the consumer’s efficiency (a fraction) in converting resource to consumer biomass.

We have already imported scipy above (import scipy as sc) so we can proceed to solve the LV model using numerical integration.

First, import scipy’s integrate submodule:

Now define a function that returns the rates of change of the resource and consumer populations at any given time.

def dRC_dt(pops, t=0):

    R = pops[0]
    C = pops[1]
    dRdt = r * R - a * R * C 
    dCdt = -z * C + e * a * R * C
    
    return np.array([dRdt, dCdt])

type(dRC_dt)

function

So dRC_dt has been stored as a function object in the current Python session, all ready to go.

Now assign some parameter values:

r = 1.
a = 0.1 
z = 1.5
e = 0.75

Define the time vector; let’s integrate from time point 0 to 15, using 1000 sub-divisions of time:

t = np.linspace(0, 15, 1000)

Note that the units of time are arbitrary here.

Set the initial conditions for the two populations (10 resources and 5 consumers per unit area), and put them into an array because dRC_dt takes an array as input.

R0 = 10
C0 = 5 
RC0 = np.array([R0, C0])

Now numerically integrate this system forward from those starting conditions:

pops, infodict = integrate.odeint(dRC_dt, RC0, t, full_output=True)

pops

array([[10.        ,  5.        ],
       [10.07578091,  4.94421976],
       [10.1529783 ,  4.88948321],
       ...,
       [ 9.99869712, 17.56204194],
       [ 9.8872779 , 17.3642589 ],
       [ 9.78000354, 17.16658946]], shape=(1000, 2))

So pops contains the result (the population trajectories). Also check what’s in infodict (it’s a dictionary with additional information)

type(infodict)

dict

infodict.keys()

dict_keys(['hu', 'tcur', 'tolsf', 'tsw', 'nst', 'nfe', 'nje', 'nqu', 'imxer', 'lenrw', 'leniw', 'mused', 'message'])

Check what the infodict output is by reading the help documentation with ?scipy.integrate.odeint. For example, you can return a message to screen about whether the integration was successful:

infodict['message']

'Integration successful.'

So it worked, great!

But we would like to visualize the results. Let’s do it using the matplotlib package.

Again, to visualize the results of your numerical simulations in Python (or for data exploration/analyses), you can use matplotlib, which uses Matlab like plotting syntax.

First import the package:

import matplotlib.pylab as p

Now open an empty figure object (analogous to an R graphics object).

f1 = p.figure()

<Figure size 640x480 with 0 Axes>

p.plot(t, pops[:,0], 'g-', label='resource density') # Plot
p.plot(t, pops[:,1]  , 'b-', label='consumer density')
p.grid()
p.legend(loc='best')
p.xlabel('Time')
p.ylabel('Population density')
p.title('consumer-resource population dynamics')
p.show()# To display the figure

../_images/90a7c567acf4b023b91ba3eaf97293a98d9646c6b1d7e0154e8641223b9edab9.png

Finally, save the figure as a pdf:

f1.savefig('results/LV_model.pdf')  # Save the figure

Tip

You can use many other graphics output formats in matplotlib; check the documentation of p.savefig.

Scientific-computing practical: Lotka–Volterra dynamics#

Create a self-contained script from the example above and save it as LV1.py in code/. It must generate the time-series figure and the following consumer–resource phase plot:

Lotka-Volterra phase plot — Fig. 14 **Generate this figure as part of the `LV1.py` script.**#

Save both figures as clearly named PDF files in results/ without displaying them on screen. From the project root, run python3 code/LV1.py. You are finished when both files are created, are non-empty, have labelled axes and legends where appropriate, and a fresh run replaces them without manual intervention.

The need for speed: profiling code#

Donald Knuth says: Premature optimization is the root of all evil.

Indeed, computational speed may not be your initial concern. Also, you should focus on developing clean, reliable, reusable code rather than worrying first about how fast your code runs. However, speed will become an issue when and if your analysis or modeling becomes complex enough (e.g., food web or large network simulations). In that case, knowing which parts of your code take the most time is useful – optimizing those parts may save you lots of time.

Profiling#

To find out what is slowing down your code you need to “profile” your code: locate the sections of your code where speed bottlenecks exist.

Profiling is easy in ipython – simply use the command:

%run -p code/your_script.py

Let’s write an illustrative program (name it profileme.py) and run it:

def my_squares(iters):
    out = []
    for i in range(iters):
        out.append(i ** 2)
    return out

def my_join(iters, string):
    out = ''
    for i in range(iters):
        out += string.join(", ")
    return out

def run_my_funcs(x,y):
    print(x,y)
    my_squares(x)
    my_join(x,y)
    return 0

run_my_funcs(10000000,"My string")

10000000 My string

Look carefully at what each of these functions does.

From IPython started at the project root, run the script with %run -p code/profileme.py. You should see something like:

        20000063 function calls (20000062 primitive calls) in 9.026 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    3.335    3.335    3.732    3.732 profileme.py:1(my_squares)
        1    2.698    2.698    5.200    5.200 profileme.py:7(my_join)
 10000001    2.502    0.000    2.502    0.000 {method 'join' of 'str' objects}
 10000008    0.397    0.000    0.397    0.000 {method 'append' of 'list' objects}
        1    0.093    0.093    9.025    9.025 profileme.py:13(run_my_funcs)
[more output]

The column headings in the profiling results are:

Column name	Meaning
`filename:lineno(function)`	The filename, line number, and function name
`ncalls`	Number of times the function was called
`tottime`	Total time spent in the function (excluding time spent in calls to sub-functions of that function)
`percall`	`tottime` divided by `ncalls`
`cumtime`	Cumulative time spent in the function and all its sub-functions
`percall`	`cumtime` divided by `ncalls`

(Please also see the official Python documentation)

The difference between tottime and cumtime columns often holds the most important information about performance. For example, in the above output, cumtime is much higher than tottime for my_join than for my_squares. This is because the sub-function (method) join is taking an additional 2.502 seconds (due to 10000001 calls), giving a cumtime of 2.698 + 2.502 = 5.2 seconds. In comparison, the sub-function (method) append is taking an additional 0.397 seconds (due to 10000008 calls), giving a cumtime of 3.335 + 0.397 = 3.732 seconds. Note also that the tottimes for join and append include the “overhead” of running their respective loop.

Thus, we can infer that the my_join function is hogging most of the time, followed by my_squares, and furthermore, that its the string method join that is slowing my_join down, and list method append that is slowing my_squares down. In other words, .joining the string, and .appending values to a list repeatedly are both not particularly fast, though the latter is less costly.

Can we do better? Yes!

Let’s try this alternative approach to writing the program (save it as profileme2.py, and again, run it):

def my_squares(iters):
    out = [i ** 2 for i in range(iters)]
    return out

def my_join(iters, string):
    out = ''
    for i in range(iters):
        out += ", " + string
    return out

def run_my_funcs(x,y):
    print(x,y)
    my_squares(x)
    my_join(x,y)
    return 0

run_my_funcs(10000000,"My string")

10000000 My string

We did two things: converted the loop to a list comprehension, and replaced the .join with an explicit string concatenation.

Now profile this program with %run -p code/profileme2.py. You should get something like:

        64 function calls (63 primitive calls) in 4.585 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    2.497    2.497    2.497    2.497 profileme2.py:2(<listcomp>)
        1    1.993    1.993    1.993    1.993 profileme2.py:5(my_join)
        1    0.094    0.094    4.584    4.584 profileme2.py:11(run_my_funcs)
[more output]

Woo hoo! So we about halved the time! Not quite enough to grab a pint, but ah well…

$\star$ Another approach would be to preallocate a numpy array instead of using a list for my_squares - Try it

You can also modify how the profiling results are displayed, and more, by using flags. For example, -s allows sorting the report by a particular column, -l limits the number of lines displayed or filters the results by function name, and -T saves the report in a text file.

For example, try

%run -p -s cumtime code/profileme2.py

This will do the profiling and display the results sorted by cumtime.

Profiling without ipython#

You may want to execute the code that you need to profile from outside Ipython (e.g., using a bash script). In this case, we you can use the module cProfile (ipython actually uses this when you do %run -p).

here is an example (run from bash):

python3 -m cProfile code/profileme2.py

Or,

python3 -m cProfile -o results/profile_results code/myscript.py

Here the profiling data are stored as results/profile_results, which can be read using the pstats module. Read the documentation for cProfile and pstats; they support several ways to examine a saved profiling report.

Quick profiling with `timeit`#

Additionally, you can use the timeit module if you want to figure out what the best way to do something specific as part of a larger program (say a particular command or a loop) might be.

Type the following code in code/timeitme.py, then run it from IPython at the project root with %run code/timeitme.py:

##############################################################################
# loops vs. list comprehensions: which is faster?
##############################################################################

iters = 1000000

import timeit

from profileme import my_squares as my_squares_loops

from profileme2 import my_squares as my_squares_lc

##############################################################################
# loops vs. the join method for strings: which is faster?
##############################################################################

mystring = "my string"

from profileme import my_join as my_join_join

from profileme2 import my_join as my_join

Note how we imported the functions using from profileme import my_squares as my_squares_loops, etc., which highlights the convenience of Python’s elegant object-oriented approach.

Now run the two sets of comparisons using timeit() in ipython and make sure every line makes sense:

%timeit my_squares_loops(iters)
%timeit my_squares_lc(iters)
%timeit (my_join_join(iters, mystring))
%timeit (my_join(iters, mystring))

Note that embedding the %timeit() commands within the script won’t work calling magic commands from inside a script is not a good idea!

Of course, a simple approach would have been to time the functions like this:

import time
start = time.time()
my_squares_loops(iters)
print(f"my_squares_loops takes {time.time() - start} s to run.")

start = time.time()
my_squares_lc(iters)
print(f"my_squares_lc takes {time.time() - start} s to run.")

But you’ll notice that if you run it multiple times, the time taken changes each time. So timeit takes a sample of runs and returns the average, which is better.

But remember, don’t go crazy with profiling for the sake of shaving a couple of milliseconds, tempting as that may be!

Vectorization#

We have now had fairly extensive practice in iteratively creating solutions to problems using for loops. Thus far all our problems have been mathematically quite straightforward, and not very computationally intensive. As you begin to move on from your taught modules into project work, you may find yourselves solving larger and more complex problems, at which point you will start to discover that for-loops have a fundamental weakness - speed!

In a nutshell, there are two issues keeping loops slow:

Producing solutions to a large problem one loop iteration at a time means that our scripts and functions spend a lot of time doing stuff other than actually solving the problem we want them to solve - e.g. shuffling data around between variables in each loop iteration, or storing the result of the current loop’s calculations in a (sometimes temporary) variable.
loops (at least in Python) are not generally able to take advantage of the parallel computing capabilities of modern-day computers, and thus must literally compute a solution one loop iteration at a time. vectorized functions tend to be optimised to spread their computations over multiple processors/threads, which can be much faster!

These issues become especially pronounced in the case of nested loops - which often appear in more spatially-explicit problems or time-dependent ones.

The following two examples will showcase the difference in runtime between a loop method and a vectorized method using numpy. The first is a relatively simple (if artificial) problem, intended to demonstrate basically at-a-glace the difference between the two approaches. The second is taken from current research on metabolic models of bacterial communities.

An example#

Let us imagine we have two simple 1D arrays $a = (a_1, a_2, ... , a_N)$ and $b = (b_1, b_2, ... ,b_N)$, each of length N, and that we want to calculate a new array $c$ in which each entry is just the product of the two corresponding entries in $a$ and $b$:

\[c = (a_1 \times b_1, a_2 \times b_2, ... , a_N \times b_N)\]

This operation is called the entrywise product of $a$ and $b$.

Below are a loop-based function and a vectorized function to calculate the entrywise product of two 1D arrays of the same length. We will test them both on larger and larger 1D arrays to see how the vectorized approach is faster.

import numpy as np

def loop_product(a, b):
    N = len(a)
    c = np.zeros(N)
    for i in range(N):
        c[i] = a[i] * b[i]   
    return c


def vect_product(a, b):
    return np.multiply(a, b)

The multiply function from numpy is a vectorized implementation of the elementwise product that we have explicitly written in the function loop_product above it. In general, numpy is an excellent choice for vectorized implementations of functions involving matrix maths (or maths using higher-dimensional analogues of matricies).

Let’s try comparing the runtimes of loop_product and vect_product on increasingly large randomly-generated 1D arrays:

import timeit

array_lengths = [1, 100, 10000, 1000000, 10000000]
t_loop = []
t_vect = []

for N in array_lengths:
    print(f"\nSet {N=}")
    #randomly generate our 1D arrays of length N
    a = np.random.rand(N)
    b = np.random.rand(N)
    
    # time loop_product 3 times and save the mean execution time.
    timer = timeit.repeat('loop_product(a, b)', globals=globals().copy(), number=3)
    t_loop.append(1000 * np.mean(timer))
    print(f"Loop method took {t_loop[-1]} ms on average.")
    
    # time vect_product 3 times and save the mean execution time.
    timer = timeit.repeat('vect_product(a, b)', globals=globals().copy(), number=3)
    t_vect.append(1000 * np.mean(timer))
    print(f"vectorized method took {t_vect[-1]} ms on average.")

Set N=1
Loop method took 0.004096003249287605 ms on average.
vectorized method took 0.0025806017220020294 ms on average.

Set N=100
Loop method took 0.0344651925843209 ms on average.
vectorized method took 0.0017308047972619534 ms on average.

Set N=10000
Loop method took 3.228893002960831 ms on average.
vectorized method took 0.01196939847432077 ms on average.

Set N=1000000

Loop method took 333.40105700772256 ms on average.
vectorized method took 2.8119451948441565 ms on average.

Set N=10000000

Loop method took 3509.266687592026 ms on average.

vectorized method took 50.03122999914922 ms on average.

Phew! That last one just exploded in terms of the time it took!

Now let’s compare the timings on a plot:

p.figure()
p.plot(array_lengths, t_loop, label="loop method")
p.plot(array_lengths, t_vect, label="vect method")
p.xlabel("Array length")
p.ylabel("Execution time (ms)")
p.legend()
p.show()

../_images/626fce959971a133c885978155154e74b903e1aa972fb8dd7cf3b1836a775b7a.png

When to vectorize?#

Thus vectorizing your code can have it running in a fraction of the time it otherwise would. Why not always vectorize then?

Generally, you should follow the same principles as with any code profiling: don’t spend time speeding up code that isn’t slow in the first place, or code which you will probably not need to run more than a small number of times.

“No free lunch!”#

There are trade-offs to vectorizing, most notably memory usage. One downside of calculating many steps simultaneously is that your computer needs to hold much more in memory in order to do it. If you try to vectorize a problem thats too large, you will probably run into memory errors. One easy example is to re-run the above example, but make it even bigger:

N = 1000000000

a = np.random.rand(N)
b = np.random.rand(N)
c = vect_product(a, b)

# if no error, remove a, b, c from memory.
del a
del b
del c

This will almost certainly return a memory error: your computer ran out of RAM while trying to store the very large arrays $a$, $b$ and $c$. If you did not get an error, try again with an extra zero in $N$.

Again, this is a rather contrived example - you probably won’t be taking element-wise products of arrays with a billion entries in your research - but more complex problems can easily become too big for memory while still remaining biologically reasonable!

Check out the CMEE module on High Performance Computing to learn about using Imperial College’s supercomputing cluster to run extremely large problems

Tip

You can use the multiprocessing package for parallelizing your code on your own computer.

Profiling and vectorization practicals#

The Lotka–Volterra model revisited#

In the original Lotka–Volterra model, the resource population grows without bound when consumers are absent. Add logistic density dependence to make resource growth slow as the resource density $R$ approaches its carrying capacity $K$:

(3)#\[\begin{align} \frac{dR}{dt} &= rR\left(1 - \frac{R}{K}\right) - aCR\\ \frac{dC}{dt} &= -zC + eaCR \end{align}\]

Here $r$ is the resource’s intrinsic growth rate, $a$ is the consumption coefficient, $z$ is consumer mortality, $e$ is conversion efficiency and $K$ is carrying capacity. If consumers disappear, the resource approaches $K$, not zero. A positive coexistence equilibrium requires $K > z/(ea)$; below that threshold, the consumer cannot persist and the resource approaches $K$.

Copy LV1.py to LV2.py and:

Accept the five parameters $r$, $a$, $z$, $e$ and $K$ as command-line arguments, converting them to numbers and rejecting missing, non-numeric or non-positive values with a clear message.
Integrate the density-dependent model and save a labelled PDF plot in results/.
Include the parameter values in the plot title or legend and print the final resource and consumer densities.
Choose one parameter set satisfying $K > z/(ea)$, predict whether both populations persist, and compare that prediction with the numerical result. Then try a value of $K$ below the threshold and explain the different outcome.

For example, from the project root:

python3 code/LV2.py 1 0.1 0.5 0.75 50

Profile the scripts with Python’s standard-library profiler rather than IPython magic:

python3 -m cProfile -s cumulative code/LV1.py
python3 -m cProfile -s cumulative code/LV2.py 1 0.1 0.5 0.75 50

Put these commands in code/run_LV.sh, make the script stop if a command fails, and run it from the project root with bash code/run_LV.sh. Compare the functions with the largest cumulative times. You are finished when both model scripts run non-interactively, produce their expected files, and the profiling commands complete successfully.

Group practical: compare R and Python vectorization#

Implement Python versions of Vectorize1.R and Vectorize2.R from the R chapter, calling them Vectorize1.py and Vectorize2.py.

Use the same input sizes and, where randomness is involved, documented seeds in both languages. Time only the equivalent calculation—not interpreter startup, imports or plotting. Run at least five timed repetitions after one untimed warm-up and report the median elapsed time. Make each R and Python script print one CSV-formatted result row containing the language, task, input size, repetitions and median seconds. Write a shell script that runs all four scripts and combines those rows into a labelled comparison table. Explain one limitation of the comparison.

Extension: discrete-time LV model#

Write a discrete-time version called LV3.py:

(4)#\[\begin{align} R_{t+1} &= R_t \left(1 + r\left(1 - \frac{R_t}{K}\right) - aC_t\right)\\ C_{t+1} &= C_t \left(1 - z + eaR_t\right) \end{align}\]

Choose and state how the program handles a step that would produce a negative population. Add the script to run_LV.sh and profile it with cProfile.

Extension: stochastic discrete-time LV model#

Copy LV3.py to LV4.py and add a normally distributed fluctuation to the resource growth rate at each step. Use numpy.random.default_rng(seed) so a supplied seed makes the run reproducible. If you add noise to both populations, draw separate fluctuations unless perfectly correlated environmental noise is an explicit modelling assumption. Save a plot, record the seed and noise standard deviation, and add the script to run_LV.sh.

Test the scripts and record a meaningful Git checkpoint. Generated plots and profiler output can be recreated and should remain untracked.

Regular expressions in Python#

Let’s shift gears now, and look at a very important tool set that you should learn, or at least be aware of — regular expressions.

regular expressions (regex) are a tool to find patterns (not just a particular sequence of characters) in strings. For example, your@email.com is a specific sequence of characters, but, in fact, all email addresses have such a pattern: alphanumeric characters, a “@”, alphanumeric characters, a “.”, alphanumeric characters. Using regex, you can search for all email addresses in a text file by searching for this pattern.

There are many uses of regex, such as:

Parsing (reading) text files and finding and replacing or deleting specific patterns
Finding DNA motifs in sequence data
Navigating through files in a directory
Extracting information from html and xml files

Thus, if you are interested in data mining, need to clean or process data in any other way, or convert a bunch of information into usable data, knowing regex is absolutely necessary.

XKCD on regex — Fig. 15 **regular expressions can really improve your quality of life.**
(Source: XKCD)#

regex packages are available for most programming languages (recall grep in UNIX; that is how regex first became popular).

Metacharacters vs. regular characters#

A regex may consist of a combination of special “metacharacters” (modifiers) and “regular” or literal characters. There are 14 metacharacters:

Metacharacter	Description
`[` `]`	Used to enclose a specific character “class” — the set of characters that you wish to match. For example, `[12]` means match target to “1” and if that does not match then match target to “2”; `[0-9]` means match to any character in range “0” to “9”
`\`	Inhibits the “specialness” of a (meta)character so that it can be interpreted literally. So, for example, use `\.` to match an actual period, and `\\` to match an actual back slash.
`.`	Match any character except line break (newline); e.g., `he..o` will match hello as well as he12o
`^`	Indicates that the string to be matched is at the start of a longer string; e.g., `^hello` will match “hello” in “hello fellow!”, but not in “fellow, hello!”
`$`	Match the end of a string; for example, `world$` will match “world” in “Hello world”, but not in “Hello world!”
`*`	Match zero or more occurrences of the character or pattern that precedes it.
`+`	Match 1 or more occurrences of the character or pattern that precedes it.
`?`	Match the preceding pattern element zero or one times
`{` `}`	Match exactly the specified number of occurrences; e.g., `.{2}` finds the first two instances of any character (except newline)
`	`
`(` `)`	Capture and group; examples of this appear below

Everything other than these metacharacters is interpreted literally (e.g., a is matched by entering a in the regex) – a regular character.

Note

The difference between *, +, and ?: * matches zero or more times, so whatever’s being repeated may not be present at all, while + requires at least one occurrence. So, for example, ra+t will match “rat” (one ‘a’ is present) and “raaat” (three “a”s), but won’t match “rt”. On the other hand, ? matches a pattern either once or zero times, so it makes the pattern matching more flexible. For example, home-?brew matches either “homebrew” or “home-brew”.

Regex special sequences#

Along with inhibiting the “specialness” of a metacharacter so that it can be interpreted literally (see examples in table above), the backslash (\) can be followed by various standard character types to denote various special sequences.

Below is a list of commonly encountered special sequences in Unicode string patterns. For a complete list look here.

Sequence	Description
`\d`	Matches any numeric (integer); this is equivalent to the regex class [0-9]
`\D`	Matches any non-digit character not covered by `\d` (i.e., match a non-digit); this is equivalent to the class [^0-9]
`\n`	Matches a newline
`\t`	Matches a tab space
`\s`	Matches any whitespace character; this is equivalent to the class [ \t\n\r\f\v]
`\S`	Matches any non-whitespace character; this is equivalent to the class [^ \t\n\r\f\v]
`\w`	Matches any “alphanumeric” character (including underscore); this is equivalent to the class [a-zA-Z0-9_]
`\W`	Matches any non-alphanumeric character not covered by `\w`, i.e., any non-alphanumeric character excluding underscore, such as `?`, `!`, `+`, `<`, etc. ; this is equivalent to the class [^a-zA-Z0-9_]

In the above table,

\r stands for a “carriage-return”, which is usually (but not always) the same as as a newline (\n);
\f stands fior “form feed” (or a page break character)
\v stands for “vertical whitespace”, which includes all characters treated as line breaks in the Unicode standard.

These are rarely encountered, but can exist in certain text files.

Note

The reason why we have specified underscore as belonging to \w (i.e., the regex set [A-Za-z0-9_]) is because this is the specification of this class in python regex in particular. This class definition is shared by most, but not all regex “flavors” (i.e., regex in other languages, such as Perl or Bash (recall grep)). The goal is to not to worry about it, but to keep in mind that \w will also match any pattern that includes one or more _’s!

Some regex examples#

So combining metacharacters, regular characters and special sequences allows you to find pretty much any pattern. Here are some examples:

regex	Description
`aX9`	match the character string aX9 exactly (case sensitively)
`aX9\n`	match the character string aX9 (case sensitively) followed by a newline
`8`	match the number 8
`\d8`	match the number 8 preceded by any decimal number
`[atgc]`	match any character listed: `a`, `t`, `g`, `c`
`at	gc`
`[^atgc]`	match any character not listed: any character except `a`, `t`, `g`, `c`
`[^a-z]`	match everything except lower case a to z
`[^Ff]`	match anything except upper or lower case f
`\w{n}`	match the preceding pattern element (any alphanumeric character) exactly `n` times
`\w{n,}`	match the preceding pattern element (any alphanumeric character) at least `n` times
`\w{n,m}`	match the preceding pattern element (any alphanumeric character) at least `n` but not more than `m` times

Tip

Metacharacters are not active inside classes. For example, [a-z$] will match any of the characters a to z, but also $, because inside a character class it loses its special metacharacter status.

Regex in Python#

regex functions in Python are in the module re.

Let’s import it:

import re

The simplest python regex function is re.search, which searches the string for match to a given pattern — returns a match object if a match is found and None if not. Thus, the command match = re.search(pat, str) finds matches of the pattern pat in the given string str and stores the search result in a variable named match.

Tip

Always put r in front of your regex — it tells python to read the regex in its “raw” (literal) form. Without raw string notation (r"text"), every backslash (\) in a regular expression would have to be prefixed with another one to escape it. read more about this here.

OK, let’s try some regexes (type all that follows in regexs.py):

my_string = "a given string"

Find a space in the string:

match = re.search(r'\s', my_string)
print(match)

<re.Match object; span=(1, 2), match=' '>

That’s only telling you that a match was found (the object was created successfully).

To see the match, use:

match.group()

' '

Now let’s try another pattern:

match = re.search(r'\d', my_string)

print(match)

None

No surprise, because there are no numeric characters in our string!

To know whether a pattern was matched, we can use an if:

MyStr = 'an example'

match = re.search(r'\w*\s', MyStr) # what pattern is this?

if match:                      
    print('found a match:', match.group()) 
else:
    print('did not find a match')    

found a match: an

Here are some more regexes (add all that follows to regexs.py):

match = re.search(r'2' , "it takes 2 to tango")
match.group()

'2'

match = re.search(r'\d' , "it takes 2 to tango")
match.group()

'2'

match = re.search(r'\d.*' , "it takes 2 to tango")
match.group()

'2 to tango'

match = re.search(r'\s\w{1,3}\s', 'once upon a time')
match.group()

' a '

match = re.search(r'\s\w*$', 'once upon a time')
match.group()

' time'

Let’s switch to a more compact syntax by directly returning the matched group (by directly appending .group() to the result).

re.search(r'\w*\s\d.*\d', 'take 2 grams of H2O').group()

'take 2 grams of H2'

re.search(r'^\w*.*\s', 'once upon a time').group() # 'once upon a '

'once upon a '

Note that *, +, and { } are all “greedy”: They repeat the previous regex token as many times as possible.

As a result, they may match more text than you want. To make it non-greedy and terminate at the first found instance of a pattern, use ?:

re.search(r'^\w*.*?\s', 'once upon a time').group()

'once '

To further illustrate greediness in regexes, let’s try matching an HTML tag:

re.search(r'<.+>', 'This is a <EM>first</EM> test').group()

'<EM>first</EM>'

But we wanted just <EM>!

It’s because + is greedy. Instead, we can make + “lazy”:

re.search(r'<.+?>', 'This is a <EM>first</EM> test').group()

'<EM>'

OK, moving on from greed and laziness…

re.search(r'\d*\.?\d*','1432.75+60.22i').group()

'1432.75'

Note \ before the ., to be able to find a literal .

Otherwise, re.search will consider it to be a regex element (. means “match any character except newline”).

A couple more examples:

re.search(r'[AGTC]+', 'the sequence ATTCGT').group()

'ATTCGT'

re.search(r'\s+[A-Z]\w+\s*\w+', "The bird-shit frog's name is Theloderma asper.").group()

' Theloderma asper'

Bird-shit Frog — Fig. 16 **In case you were wondering what *Theloderma asper*, the “bird-shit frog”, looks like.** Samraat snapped this one in a North-East Indian rainforest ages ago.#

How about looking for email addresses in a string? For example, let’s try matching a string consisting of an academic’s name, email address and research area or interest (no need to type this into any Python file):

MyStr = 'Samraat Pawar, s.pawar@imperial.ac.uk, Systems biology and ecological theory'
match = re.search(r"[\w\s]+,\s[\w\.@]+,\s[\w\s]+",MyStr)
match.group()

'Samraat Pawar, s.pawar@imperial.ac.uk, Systems biology and ecological theory'

Note the use of [ ]’s: for example, [\w\s] ensures that any combination of word characters and spaces is found.

Let’s see if this regex works on a different pattern of email addresses:

MyStr = 'Samraat Pawar, s-pawar@imperial.ac.uk, Systems biology and ecological theory'

match = re.search(r"[\w\s]+,\s[\w\.@]+,\s[\w\s]+",MyStr)
match.group()

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[290], line 2
      1 match = re.search(r"[\w\s]+,\s[\w\.@]+,\s[\w\s]+",MyStr)
----> 2 match.group()

AttributeError: 'NoneType' object has no attribute 'group'

Nope! So let’s make the email address part of the regex more robust:

match = re.search(r"[\w\s]+,\s[\w\.-]+@[\w\.-]+,\s[\w\s]+",MyStr)
match.group()

'Samraat Pawar, s-pawar@imperial.ac.uk, Systems biology and ecological theory'

Formative exercises: regular expressions#

Use these exercises to practise and record any examples or failures you find.

Try the name pattern r'[\w\s]+' on names containing punctuation such as ?, +, a hyphen or an apostrophe. Record which cases match and revise the pattern for the name rules you choose.
Translate these regular expressions into plain English:
- r'^abc[ab]+\s\t\d'
- r'^\d{1,2}\/\d{1,2}\/\d{4}$'
- r'\s*[a-zA-Z,\s]+\s*'
Write a regex for dates in YYYYMMDD format with a year from 1901 to 2099, a month from 01 to 12 and a day from 01 to 31. Test valid boundary examples and invalid values such as month 00, month 13, day 00 and day 32. This checks the numeric format only; rejecting impossible calendar dates such as 20230231 requires date parsing rather than regex alone.

For every exercise, save at least one example that should match and one that should not. You are finished when your predicted and observed results agree and you can explain each part of your patterns.

Grouping regex patterns#

You can group regex patterns into meaningful blocks using parentheses. Let’s look again at the example of finding email addresses.

MyStr = 'Samraat Pawar, s.pawar@imperial.ac.uk, Systems biology and ecological theory'
match = re.search(r"[\w\s]+,\s[\w\.-]+@[\w\.-]+,\s[\w\s]+",MyStr)
match.group()

'Samraat Pawar, s.pawar@imperial.ac.uk, Systems biology and ecological theory'

Without grouping the regex:

match.group(0)

'Samraat Pawar, s.pawar@imperial.ac.uk, Systems biology and ecological theory'

Now create groups using ( ):

match = re.search(r"([\w\s]+),\s([\w\.-]+@[\w\.-]+),\s([\w\s&]+)",MyStr)
if match:
    print(match.group(0))
    print(match.group(1))
    print(match.group(2))
    print(match.group(3))

Samraat Pawar, s.pawar@imperial.ac.uk, Systems biology and ecological theory
Samraat Pawar
s.pawar@imperial.ac.uk
Systems biology and ecological theory

Nice! This is very handy for extracting specific patterns from text data. Note that we excluded the ,’s and the \s’s from the grouping parentheses because we don’t want them to be returned in the match group list.

Download and examine re4.py for another example of parsing email addresses using regular expressions.

Useful `re` commands#

Here are some important functions in the re module:

Command	What it does
`re.search(reg, text)`	Scans the string and finds the first match of the pattern, returning a `match` object if successful and `None` otherwise.
`re.match(reg, text)`	Like `re.search`, but only matches the beginning of the string.
`re.compile(reg)`	Compiles (stores) a regular expression for repeated use, improving efficiency.
`re.split(ref, text)`	Splits the text by the occurrence of the pattern described by the regular expression.
`re.findall(ref, text)`	Like `re.search`, but returns a list of all matches. If groups are present, returns a list of groups.
`re.finditer(ref, text)`	Like `re.findall`, but returns an iterator containing the match objects over which you can iterate. Useful for “crawling” efficiently through text till you find all necessary number of matches.
`re.sub(ref, repl, text)`	Substitutes each non-overlapping occurrence of the match with the text in `repl`.

Many of these commands also work on whole contents of files. We will look at an example of this below. Let us try some particularly useful applications of some of these commands.

Finding all matches#

Above we used re.search() to find the first match for a pattern. In many scenarios, you will need to find all the matches of a pattern. The function re.findall() does precisely this and returns all matches as a list of strings, with each string representing one match.

Let’s try this on an extension of the email example above for some data with multiple addresses:

MyStr = "Samraat Pawar, s.pawar@imperial.ac.uk, Systems biology and ecological theory; Another academic, a-academic@imperial.ac.uk, Some other stuff thats equally boring; Yet another academic, y.a_academic@imperial.ac.uk, Some other stuff thats even more boring"

Now re.findall() returns a list of all the emails found:

emails = re.findall(r'[\w\.-]+@[\w\.-]+', MyStr) 
for email in emails:
    print(email)

s.pawar@imperial.ac.uk
a-academic@imperial.ac.uk
y.a_academic@imperial.ac.uk

Nice!

Finding in files#

You will generally be wanting to apply regex searches to whole files. You might be tempted to write a loop to iterate over the lines of the file, calling re.findall() on each line. However, re.findall() can return a list of all the matches in a single step.

Let’s try finding all species names that correspond to Oaks in a data file:

with open('data/test_oaks_data.csv', 'r') as f:
    oak_data = f.read()

found_oaks = re.findall(r"Q[\w\s].*\s", oak_data)

found_oaks

['Quercus, robur\n', 'Quercus, cerris\n', 'Quercus, petraea\n']

This works because f.read() returns the complete contents as one string. The with block closes the file before the regular-expression search begins.

Groups within multiple matches#

Grouping pattern matches using ( ) as you learned above, can be combined with re.findall(). If the pattern includes two or more groups, then instead of returning a list of strings, re.findall() returns a list of tuples. Each tuple represents one match of the pattern, and inside the tuple is group(1), group(2), etc.

Let’s try it:

MyStr = "Samraat Pawar, s.pawar@imperial.ac.uk, Systems biology and ecological theory; Another academic, a.academic@imperial.ac.uk, Some other stuff thats equally boring; Yet another academic, y.a.academic@imperial.ac.uk, Some other stuff thats even more boring"

found_matches = re.findall(r"([\w\s]+),\s([\w\.-]+@[\w\.-]+)", MyStr)
found_matches

[('Samraat Pawar', 's.pawar@imperial.ac.uk'),
 (' Another academic', 'a.academic@imperial.ac.uk'),
 (' Yet another academic', 'y.a.academic@imperial.ac.uk')]

for item in found_matches:
    print(item)

('Samraat Pawar', 's.pawar@imperial.ac.uk')
(' Another academic', 'a.academic@imperial.ac.uk')
(' Yet another academic', 'y.a.academic@imperial.ac.uk')

Extracting text from webpages#

OK, let’s step up the ante here. How about extracting text from a web page to create your own data? Let’s try extracting data from this page.

You will need a new package urllib3. Install it, and import it (also import re if needed).

import urllib3

conn = urllib3.PoolManager() # open a connection
r = conn.request('GET', 'https://www.imperial.ac.uk/silwood-park/academic-staff/') 
webpage_html = r.data #read in the webpage's contents

This is returned as bytes (not strings).

type(webpage_html)

bytes

So decode it (remember, the default decoding that this method applies is utf-8):

My_Data  = webpage_html.decode()
#print(My_Data)

That’s a lot of potentially useful information! Let’s extract all the names of academics:

pattern = r"Dr\s+\w+\s+\w+"
regex = re.compile(pattern) # example use of re.compile(); you can also ignore case  with re.IGNORECASE 
for match in regex.finditer(My_Data): # example use of re.finditer()
    print(match.group())

Again, nice! However, its’ not perfect.

You can improve this by:

Extracting Prof names as well
Eliminating the repeated matches
Grouping to separate title from first and second names
Extracting names that have unexpected characters, such as in hyphenated names (a “-” in the name)

Try making these improvements.

Of course, you can match and extract other types of patterns as well, such as urls and email addresses (though this example web page does not have email addresses).

Replacing text#

Using the same web page data, let’s try using the re.sub command on the same web page data (My_Data) to replace text:

New_Data = re.sub(r'\t'," ", My_Data) # replace all tabs with a space
# print(New_Data)

Regex capstone: Blackbirds#

Download blackbirds.py into code/ and blackbirds.txt into data/. Complete the script by extracting the kingdom, phylum and two-word scientific name for every species record.

Required result:

four records are found;
each result contains exactly three text fields: kingdom, phylum and scientific name;
the scientific name contains the genus and species epithet;
the script prints the records in a consistent, readable format without unrelated matches.

Keep at least one positive and one negative check for the core pattern. From the project root, run python3 code/blackbirds.py and compare the four results with the examples in the starter file.

Using Python to build workflows#

You can use Python to build an automated data analysis or simulation workflow that involves multiple languages, especially the ones you have already learnt: r, $\LaTeX$, and UNIX bash. For example, you could, in theory, write a single Python script to generate and update your masters dissertation, tables, plots, and all. Python is ideal for building such workflows because it has packages for practically every purpose.

Thus this topic may be useful for your Miniproject, which will involve building a reproducible computational workflow.

Using `subprocess`#

For building a workflow in Python the subprocess module is key. With this module you can run non-Python commands and scripts, obtain their outputs, and also crawl through and manipulate directories.

First, import the module (this is part of the Python standard library, so you won’t need to install it):

import subprocess

Running processes#

There are two main ways to run commands through subprocess: run (available in Python 3.5 onwards) for basic usage, and Popen (Process open) for more advanced usage. We will work directly with popen because run() is a wrapper around Popen. Using Popen directly gives more control over how the command is run, and how its input and output are processed.

Let’s try running some commands in the UNIX bash.

$\star$ In a terminal, first cd to your code directory, launch ipython3, then and type:

p = subprocess.Popen(["echo", "I'm talkin' to you, bash!"], stdout=subprocess.PIPE, stderr=subprocess.PIPE)

This creates an object p, from which you can extract the output and other information of the command you ran. Before we do anything more, let’s look at our subprocess.popen call carefully.

The command line arguments were passed as a list of strings, which avoids the need for escaping quotes or other special characters that might be interpreted by the shell (for example, in this case, there are apostrophes in the string that is being echoed in bash).
stdout is the output from the process “spawned” by your command. This is a sequence of bytes (which you will need to decode - more on this below).
stderr is the process’s error-output stream; it is not an error code. After the process finishes, p.returncode contains its exit status (zero conventionally means success). subprocess.PIPE creates a pipe through which Python can capture a child process’s output.

stdout, stderr = p.communicate()

stderr

b''

b'' is an empty bytes object: echo wrote nothing to its error stream. The b prefix denotes bytes, not bits. By default, stdout and stderr returned by p.communicate() are bytes. After communicate() returns, inspect p.returncode to determine whether the command succeeded.

Now check what’s in stdout:

stdout

b"I'm talkin' to you, bash!\n"

Let’s decode the bytes to text and print them.

print(stdout.decode())

I'm talkin' to you, bash!

You can instead pass text=True (the modern spelling of universal_newlines=True) so that these outputs are returned as decoded strings, with line endings converted to \n. You can specify the character encoding explicitly with encoding="utf-8". For more information, see the subprocess documentation.

Let’s try something else:

p = subprocess.Popen(["ls", "-l"], stdout=subprocess.PIPE)
stdout, stderr = p.communicate()

Now use print(stdout.decode()) to see the output.

recall that the ls -l command lists all files in a long listing format.

You can also call Python itself from bash (!):

p = subprocess.Popen(["python3", "boilerplate.py"], stdout=subprocess.PIPE, stderr=subprocess.PIPE) # A bit silly! 
stdout, stderr = p.communicate()

print(stdout.decode())

This should give you the output of your boilerplate.py script (remember, your current path must be the code/ directory containing the script; see Writing and running Python programs).

Similarly, to compile a $\LaTeX$ document (using pdflatex in this case), you can do something like:

subprocess.run(["pdflatex", "yourlatexdoc.tex"], check=True)

Running `R`#

R is likely to be an important part of your project’s analysis, particularly for statistical analyses and plotting. You can run R from Python easily:

$\star$ Create an R script file called test_r.R in your code directory with the following content:

print("Hello, this is R!")

Then run it from your bash terminal (not from the Python terminal), starting at the coursework project root: Rscript code/test_r.R. This should return in your terminal:

[1] "Hello, this is R!"

Then try Rscript --verbose code/test_r.R. This should return something like:

  '/usr/lib/r/bin/r --no-echo --no-restore --file=code/test_r.R'

[1] "Hello, this is R!"

Thus the --verbose flag to Rscript prints additional information about the R script’s execution.

Note

The location of Rscript can differ between Unix/Linux systems. To locate it, run command -v Rscript in the shell (not in Python). You can then use the returned path if necessary.

Now try:

Rscript --verbose code/test_r.R > results/test_r.rout 2> results/test_r_errFile.rout

Then check results/test_r.rout and results/test_r_errFile.rout. The former should contain the output of the R script, and the latter the additional information produced by the --verbose flag.

Now try running this using subprocess instead:

with open("results/test_r.rout", "w") as output_file, \
     open("results/test_r_errFile.rout", "w") as error_file:
    result = subprocess.run(
        ["Rscript", "--verbose", "code/test_r.R"],
        stdout=output_file,
        stderr=error_file,
        check=False,
    )
print(result.returncode)

Here subprocess.run() waits for Rscript to finish. The command and its arguments are passed as a list, so no shell is needed. Standard output and standard error are sent directly to separate files, and result.returncode records the exit status. A return code of zero conventionally means success; a non-zero value indicates a failure. Passing check=True instead would raise subprocess.CalledProcessError automatically for a non-zero return code.

Now, create a script test_r.py with the following content :

import subprocess

with open("results/test_r.rout", "w") as output_file, \
     open("results/test_r_errFile.rout", "w") as error_file:
    result = subprocess.run(
        ["Rscript", "--verbose", "code/test_r.R"],
        stdout=output_file,
        stderr=error_file,
        check=False,
    )

print("Return code:", result.returncode)
if result.returncode != 0:
    raise SystemExit(result.returncode)

Return code: 0

From the project root, run python3 code/test_r.py (or %run code/test_r.py in IPython) and check results/test_r.rout and results/test_r_errFile.rout.

Also check what happens if you run (type directly in ipython or python console):

result = subprocess.run(
    ["Rscript", "--verbose", "code/NonExistScript.R"],
    capture_output=True,
    text=True,
)
print("Return code:", result.returncode)
print(result.stderr)

Return code: 2
running
  '/usr/lib/R/bin/R --no-echo --no-restore --file=code/NonExistScript.R'

What do you see on the screen?

The non-zero return code signals failure, while result.stderr contains R’s diagnostic message.

Handling directory and file paths#

Use the standard-library os module directly for operating-system-independent path handling. For example:

import os
os.path.join('directory', 'subdirectory', 'file')

'directory/subdirectory/file'

The result would be appropriately different on Windows (with backslashes instead of forward slashes).

Note that in all cases you can “catch” the output of subprocess so that you can then use the output within your Python script. A simple example, where the output is a platform-dependent directory path, is:

MyPath = os.path.join('directory', 'subdirectory', 'file')
MyPath

'directory/subdirectory/file'

Tip

Explore available names in IPython by pressing the Tab key after typing subprocess., os. or os.path.. These are separate standard-library modules with different responsibilities.

Workflow practical: R and filesystem automation#

Run an R analysis from Python#

Download fmr.R into code/ and nagy_et_al_1999.csv into data/. Read the R script, then write run_fmr_r.py that runs Rscript code/fmr.R from the coursework project root using subprocess.run() without shell=True.

Capture standard output and standard error as text, print both streams with labels, and print whether the command succeeded based on its return code. On success, verify that results/fmr_plot.pdf exists and is non-empty; on failure, exit with a non-zero status. Your wrapper is complete when a normal run reports success and deliberately using a nonexistent R script reports failure without a Python traceback.

Traverse directories with `os`#

Download using_os.py into code/ and complete its three tasks using the os module directly, including os.walk()—do not access os through subprocess. For a small temporary test directory, predict and check the returned files and directories before scanning a larger tree. Keep files and directories as separate results and make the case-sensitive and case-insensitive name rules explicit.

Readings & Resources#

Python packages#

A Python package is simply a directory of Python modules (quite like an R package). A number of packages, such as the following particularly useful ones, are always available as standard libraries (just require import from within Python or ipython):

io: file input-output with .csv, .txt, etc.
subprocess: to run other programs, including multiple ones at the same time, including operating system-dependent functionality
math: for mathematical functions
sqlite3: for manipulating and querying sqlite databases

Scores of other packages are accessible by explicitly installing them using sudo apt install python-packagename (as you did previously) or by using pip. Some particularly mentionable ones are:

sciPy contains a wide array of numerical tools for scientific computing, including NumPy for efficient data crunching
matplotlib for plotting (very matlab-like, requires scipy) (all packaged in pylab)
pandas for manipulating data, and comes with a DataFrame object similar to the R data frame.
scikit-learn for applying different machine learning algorithms to data
ipython is an enhanced Python terminal (which you are currently using)
jupyter is an interactive notebook environment for data analysis, visualization, and creation of documents that can be shared. This course is written entirely as Jupyter notebooks.
scrapy allows you to write efficient programs that crawl through web sites and extract data from them
beautifulsoup4 can be used for parsing HTML and XML to extract data (can do a lot of what scrapy does)
biopython for bioinformatics. Check out the worked examples in particular.

We will use many of these in the Advanced Python Chapter. You can look up https://docs.python.org/3/library/index.html – read about the packages you think will be important to you.

General#

M. Fowler, refactoring: improving the design of existing code (2000), Addison-Wesley Professional
R. C. Martin, Clean Code: A Handbook of Agile Software Craftsmanship (2009), Prentice Hall
Browse the Python tutorial
Shaw, Zed. 2017. Learn Python 3 the Hard Way: A Very Simple Introduction to the Terrifyingly Beautiful World of Computers and Code. Boston : Addison-Wesley (likely available in your University library)
Python style guide from the Hitchhiker’s Guide to Python
A very detailed style guide: the Google Python Style Guide

Ipython#

The Ipython documentation
Cookbooks can be very useful

Python data structures#

read more about native Python data structures here

Elements of the Python program#

read more about the Python shebang here
See the official docstring conventions

Errors and Debugging#

https://docs.python.org/3/tutorial/errors.html

Functions, Modules, and Classes#

For functions vs. modules vs. classes in Python, have a look at this

Scientific computing#

In general, scores of good module/package-specific cookbooks are out there — google “cookbook” along with the name of the package you are interested in (e.g., “scipy cookbook”).
For SciPy, read the official documentation; in particular, read about the scipy modules you think will be important to you.
The “ecosystem” for Scientific computing in Python: https://scipy-lectures.org/
Many great examples of applications in the scipy cookbook
Scipy stats: https://docs.scipy.org/doc/scipy/tutorial/stats.html
A Primer on Scientific Programming with Python http://www.springer.com/us/book/9783642549595; Multiple copies of this book are available from the central library and can be requested to Silwood from the IC library website. You can also find a pdf - seach online.

Regular expressions#

Python regex documentation: https://docs.python.org/3.6/howto/regex.html
Python regex examples in the standard library docs: https://docs.python.org/3/library/re.html#regular-expression-examples
Python’s regular expression HOWTO
Good intro to regex, tips and a great array of canned solutions: http://www.regular-expressions.info
Use and abuse of regex: https://blog.codinghorror.com/regex-use-vs-regex-abuse/

Other stuff#

The matplotlib website
Some of you might find the Python package biopython particularly useful — check out http://biopython.org/, and especially, the cookbook
Alternatives to matplotlib for plotting in Python: Seaborn and Plotly

Biological Computing in Python

Contents

Biological Computing in Python#

Introduction#

From a task to an algorithm#

Why Python?#

What you will learn#

Some terminology#

Getting started with Python#

Baby steps#

The Zen of Python and good programming practices#

ipython#

Magic commands#

Determining an object’s type#

Python variables#

Python operators#

Assigning and manipulating variables#

Writing readable Python code#

Testing/running blocks of code#

Python data structures#

Lists#

Tuples#

Sets#

Dictionaries#

Copying mutable objects#

Python with strings#

Getting help#

Exercises: Python data structures#

Loops#

Iterator vs Iterable#

Some loops examples#

Functions#

Hints about inputs and outputs#

Control flow tools#

Conditionals#

More examples of loops and conditionals combined#

From an algorithm to Python#

Comprehensions#

Exercises: loops and comprehensions#

Variable scope#

Global variables#

Importance of the return directive#

Exercises: functions and control flow#

Writing and running Python programs#

Modules and code compartmentalization#

Running Python scripts#

Using the UNIX shell#

Let the script find the project#

Structuring Python programs#

Components of the Python program#

The main function and main guard#

Module and function docstrings#

Further reference: the shebang#

Special names and ordinary names#

Functions and modules are different#

A program-with-control-flows example#

Modules and control_flow.py#

Passing inputs from the command line#

Inspecting sys.argv#

Reporting success or failure#

Python input and output#

Safely opening files using with open()#

Handling CSV files#

Using column names instead of positions#

Practical: programs, modules, and file input/output#

Testing, debugging, and errors#

Unit testing#

Unit testing with doctest#

From doctest to a separate test script#

Examine a meaningful boundary test#

Other unit testing approaches#

Debugging#

Paranoid programming: pause before the failure#

Handling expected errors#

A number that converts may still be unusable#

Integrated practical: summarising observations#

Designing an algorithm for imperfect data#

Example: summarize observations by site#

Independent exercise: summarise observations by site#

Reuse the module-and-test-script pattern#

Importance of the `return` directive#

The `main` function and main guard#

Modules and `control_flow.py`#

Inspecting `sys.argv`#

Safely opening files using `with open()`#

Unit testing with `doctest`#

`numpy` matrices#

The `scipy` package#

Numerical integration using `scipy`#

Quick profiling with `timeit`#

Useful `re` commands#

Using `subprocess`#

Running `R`#

Traverse directories with `os`#