Please do use the course Etherpad:

• https://pad.carpentries.org/2019-05-02-standrews

• To learn basic concepts of programming (in Python)

• we need something ;)
• free, well-documented, and cross-platform
• large academic userbase
• many libraries for specialist work

• No. ;)
  • The lesson is general, it’s just taught in Python
  • The principles are the same in nearly all languages
  • If your colleagues/field settled on another language(s), maybe learn that (instead or as well)
  • (language wars are unproductive… ;) )

Before we begin…

• make a neat working environment using the terminal
• obtain data

cd ~/Desktop        # or any convenient location
mkdir pni           # (python novice inflammation)
cd pni

Before we begin…

• make a neat working environment
• obtain data

cp 2019-05-02-standrews/lessons/python/files/python-novice-inflammation-data.zip ./
unzip python-novice-inflammation-data.zip
cp 2019-05-02-standrews/lessons/python/files/python-novice-inflammation-code.zip ./
unzip python-novice-inflammation-code.zip

(you can download files via Etherpad: http://pad.software-carpentry.org/2019-05-02-standrews)

$ python
Python 3.6.3 |Anaconda custom (64-bit)| (default, Oct  6 2017, 12:04:38) 
[GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> 

• Python’s console is a read-evaluate-print-loop, just like the shell

>>> 3 + 5
8
>>> 12 / 7
1.7142857142857142
>>> 2 ** 16
65536

• To do interesting things, we want persistent values

>>> name = "Samia"
>>> name
'Samia'
>>> print(name)
Samia

weight_kg = 55
print(weight_kg)
2.2 * weight_kg
print("weight in pounds", 2.2 * weight_kg)

mass = 47.5
age = 122
mass = mass * 2.0
age = age - 20

1. mass == 47.5, age == 122
2. mass == 95.0, age == 102
3. mass == 47.5, age == 102
4. mass == 95.0, age == 122

first, second = 'Grace', 'Hopper'
third, fourth = second, first
print(third, fourth)

1. Hopper Grace
2. Grace Hopper
3. "Grace Hopper"
4. "Hopper Grace"

• Python contains many powerful, general tools
• Specialised tools are contained in libraries or packages
• We call on libraries/packages, when needed
• Packages are loaded with import
• Packages are shared via repositories, e.g. PyPI and conda

>>> import numpy

• numpy provides a function loadtxt() to load tabular data:

numpy.loadtxt(fname='data/inflammation-01.csv', delimiter=',')

>>> numpy.loadtxt(fname='data/inflammation-01.csv', delimiter=',')
array([[ 0.,  0.,  1., ...,  3.,  0.,  0.],
       [ 0.,  1.,  2., ...,  1.,  0.,  1.],
       [ 0.,  1.,  1., ...,  2.,  1.,  1.],
       ..., 
       [ 0.,  1.,  1., ...,  1.,  1.,  1.],
       [ 0.,  0.,  0., ...,  0.,  2.,  0.],
       [ 0.,  0.,  1., ...,  1.,  1.,  0.]])

• The matrix is truncated to fit the screen
• ... indicate missing rows or columns
• If there are no significant digits, they are not shown (1 == 1. == 1.0)

>>> type(data)
>>> data.dtype

• Creating the array created information, too
• Info stored in attributes (or members) that belong to data
• data.<attribute> e.g. data.shape

>>> print(data.dtype)
>>> print(data.shape)

• We often work with subsets of data
  • individual rows (patients)
  • individual columns (days)

>>> print('first value in data:', data[0, 0])
first value in data: 0.0
>>> print('middle value in data:', data[30, 20])
middle value in data: 13.0

• To get a range of data from the array, index with [] and specify start and end indices

>>> print(data[0:4, 0:10])

>>> small = data[:3, 36:]
>>> print(small)

We can take slices of any series, not just arrays.

>>> element = 'oxygen'

What is the value of element[:4]?

1. oxyg
2. gen
3. oxy
4. en

• arrays know how to perform operations on their values
• +, -, *, /, etc. are elementwise

>>> doubledata = data * 2.0
>>> print(data[:3, 36:])
 [[ 2.  3.  0.  0.]
 [ 1.  1.  0.  1.]
 [ 2.  2.  1.  1.]]
>>> print(doubledata[:3, 36:])
 [[ 4.  6.  0.  0.]
 [ 2.  2.  0.  2.]
 [ 4.  4.  2.  2.]]

• numpy provides functions to operate on arrays

>>> print(numpy.mean(data))
6.14875
>>> maxval = numpy.max(data)
>>> print('maximum inflammation:', maxval)
maximum inflammation: 20.0
>>> minval = data.min()
>>> print('minimum inflammation:', minval)
minimum inflammation: 0.0

• We want to summarise inflammation by patient

>>> patient_0 = data[0, :] # temporary variable
>>> print('maximum inflammation for patient 0:', numpy.max(data[0, :]))
maximum inflammation for patient 0: 18.0

• What if we need data for each patient or average on each day?

• numpy functions take an axis= parameter: 0 (columns) or 1 (rows)

>>> print(numpy.max(data, axis=1))    # max by patient
>>> print(data.mean(axis=0))          # mean by day

“The purpose of computing is insight, not numbers” - Richard Hamming

• Visualisation deserves a course to itself
• Here’s one I prepared earlier (for the Software Sustainability Institute):

• Research Data Visualisation Workshop

>>> import matplotlib.pyplot

>>> image = matplotlib.pyplot.imshow(data)
>>> matplotlib.pyplot.show()

• small values are dark blue, large values are yellow
• inflammation rises and falls over a 40-day period

• matplotlib.pyplot.plot() renders a line graph

>>> ave_inflammation = numpy.mean(data, axis=0)
>>> ave_plot = matplotlib.pyplot.plot(ave_inflammation)
>>> matplotlib.pyplot.show()

• The plot of .mean() looks artificial
• Look at other statistics to gain insight

>>> max_plot = matplotlib.pyplot.plot(numpy.max(data, axis=0))
>>> matplotlib.pyplot.show()
>>> min_plot = matplotlib.pyplot.plot(numpy.min(data, axis=0))
>>> matplotlib.pyplot.show()

• create a figure (fig) with fig = matplotlib.pyplot.figure()
• add subplots to fig with ax = fig.add_subplot()
• set labels on a subplot with ax.set_ylabel()
• plot data to a subplot with ax.plot()

• We wrote some code that plots values of interest from a single dataset

• BUT we’re soon going to receive dozens of datasets to plot
• So we need to make the computer iterate over the data

• Suppose we wanted to spell a word out with one letter on each line

word = "lead"
print(word[0])
print(word[1])
print(word[2])
print(word[3])

• for loops perform actions for every item in a collection

>>> word = "lead"
>>> for char in word:
...     print(char)
... 
l
e
a
d

for variable in collection:
    <do things using variable>

• The for loop statement ends in a colon, :
• The code block is indented with a tab (\t)

>>> length = 0
>>> for vowel in 'aeiou':
...     length = length + 1
... 
>>> print("There are", length, "vowels") 

>>> letter = "z"
>>> print(letter)
z
>>> for letter in "abc":
...     print(letter)
... 
>>> print("after the loop, letter is:", letter)

• It returns a range type that can be iterated over in a loop

>>> seq = range(3)
>>> print("Range is:", seq)
>>> for val in seq:
...     print(val)
>>> seq = range(2, 5)
>>> seq = range(3, 10, 3)
>>> seq = range(10, 0, -1)

• lists are a built in Python datatype for storing multiple values
• Denoted by square brackets, comma-separated
  • iterable lists of values
  • indexed and sliced like arrays

>>> odds = [1, 3, 5, 7]
>>> print("odds are:", odds)
>>> print('first and last:', odds[0], odds[-1])
>>> for number in odds:
...     print(number)

• lists, like strings, are sequences

>>> names = ["Curie", "Darwing", "Turing"] # typo in Darwin's name
>>> print("names is originally:", names)
>>> names[1] = 'Darwin'    # correct the name
>>> print('final value of names:', names)
>>> name = "darwin"
>>> name[0] = "D"

>>> my_list = [1, 2, 3, 4]
>>> your_list = my_list
>>> print("my list:", my_list)
>>> my_list[1] = 0
>>> print("your list:", your_list)

>>> my_list = [1, 2, 3, 4]           # original list
>>> your_list = my_list[:]           # copy 1
>>> your_other_list = list(my_list)  # copy 2
>>> print("my_list:", my_list)
>>> my_list[1] = 0                   # change element
>>> print("my_list:", my_list)
>>> print("your_list:", your_list)
>>> print("your_other_list:", your_list)

• lists are Python objects and have useful functions (methods)

>>> print(odds)
[1, 3, 5, 7]
>>> odds.append(9)
>>> print("odds after adding a value:", odds)
>>> odds.reverse()
>>> print("odds after reversing the list:", odds)
>>> odds.pop()
>>> print("odds after popping:", odds)

• For numbers, + means add
• For lists, + means concatenate

>>> vowels = ['a', 'e', 'i', 'o', 'u']
>>> vowels_welsh = ['a', 'e', 'i', 'o', 'u', 'w', 'y']
>>> print(vowels + vowels_welsh)
>>> counts = [2, 4, 6, 8, 10]
>>> repeats = counts * 2
>>> print(repeats)

• We often want to do <something> if some condition is true
• To do this, we can use an if statement:

if <condition>:
    <executed if condition is True>

>>> num = 37
>>> if num > 100:
...     print('greater')
... 
>>> num = 149
>>> if num > 100:
...     print('greater')
... 

• An if statement executes code if the condition evaluates as true
  • But what if the condition evaluates as false?

if <condition>:
    <executed if condition is True>
else:
    <executed if condition is not True>

>>> num = 37
>>> if num > 100:
...     print('greater')
... else:
...     print('not greater')
... 

• We can chain tests together using elif (else if)

if <condition1>:
    <executed if condition1 is True>
elif <condition2>:
    <executed if condition2 is True and condition1 is not True>
else:
    <executed if no conditions True>

>>> num = -3
>>> if num > 0:
...     print(num, "is positive")
... elif num == 0:
...     print(num, "is zero")
... else:
...     print(num, "is negative")
... 

• Operators include and, or and not

>>> if (4 > 0) or (2 > 0):
...     print('at least one part is true')
... else:
...     print('both parts are false')
... 

>>> if 4 > 5:
...     print('A')
... elif 4 == 5:
...     print('B')
... elif 4 < 5:
...     print('C')

1. A
2. B
3. C
4. B and C

>>> print(1 == 1)
>>> print(1 == 2)
>>> print('a' in 'toast')
>>> print('b' in 'toast')
>>> print(1 in [1, 2, 3])
>>> print(1 in range(3))

• We have several files of inflammation study data
  • We want to visualise/analyse each of them
  • We know how to load, visualise, loop over, and make decisions on the data

• The list can be filtered with a for loop or list comprehension
• Our data is in the data directory

import os

# Get a list of inflammation data files
files = []
for fname in os.listdir('data'):
  if 'inflammation' in fname:
    files.append(fname)
print("Inflammation data files:", files)

$ nano analyse_files.py

import matplotlib.pyplot
import numpy as np
import os
...

$ python analyse_files.py

• The os.listdir() function only returns filenames, not the path (relative or absolute)

files = []
for fname in os.listdir('data'):
  if 'inflammation' in fname:
    files.append(os.path.join('data', fname))
print("Inflammation data files:", files)

$ python analyse_files.py
Inflammation data files: ['data/inflammation-05.csv', …]

• list of paths to the data files with os
• load data from a file with np.loadtxt()
• calculate summary statistics with np.mean(), np.max(), etc.
• create figures with matplotlib
• create arrays of figures with .add_subplot()

$ nano analyse_files.py
$ python analyse_files.py 
Analysing data/inflammation-05.csv
Writing image to data/inflammation-05.png
Analysing data/inflammation-11.csv
Writing image to data/inflammation-11.png
[…]

1. The maximum values rose and fell as straight lines
2. The minimum values are consistently zero

$ nano analyse_files.py

    # Test for suspicious maxima
    if numpy.max(data, axis=0)[0] == 0 and numpy.max(data, axis=0)[20] == 20:
        print("Suspicious-looking maxima!")

$ python analyse_files.py

$ nano analyse_files.py

    # Test for suspicious maxima
    if numpy.max(data, axis=0)[0] == 0 and numpy.max(data, axis=0)[20] == 20:
        print("Suspicious-looking maxima!")
    elif numpy.sum(numpy.min(data, axis=0)) == 0:
        print("Minima sum to zero!")

$ python analyse_files.py

$ nano analyse_files.py

    # Test for suspicious maxima
    if numpy.max(data, axis=0)[0] == 0 and numpy.max(data, axis=0)[20] == 20:
        print("Suspicious-looking maxima!")
    elif numpy.sum(numpy.min(data, axis=0)) == 0:
        print("Minima sum to zero!")
    else:
        print("Seems OK!")

$ python analyse_files.py

At the command-line, start Jupyter notebook:

jupyter notebook

• Give your notebook a name (functions)

• A Jupyter cell can have one of several types

• We have code to plot values of interest from multiple datasets

• It’s not flexible enough to deal with thousands of files
• We can’t modify it easily

\[y = f(x)\]

• \(f()\) is the function
• \(x\) is an input (or inputs)

• \(y\) is the returned value, or output(s)

• The output \(y\) depends in some way on the value of \(x\) - defined by \(f()\).

\[f(x) = ((x - 32) \times \frac{5}{9}) + 273.15\]

• Calling fahr_to_kelvin() in the notebook is the same as calling any other function

print('freezing point of water:', fahr_to_kelvin(32))
print('boiling point of water:', fahr_to_kelvin(212))

def kelvin_to_celsius(temp):
  return temp - 273.15

print('freezing point of water', kelvin_to_celsius(273.15))

• We could convert F (temp_f) to C (temp_c) by executing the code:

temp_c = kelvin_to_celsius(fahr_to_kelvin(temp_f))

def fahr_to_celsius(temp_f):
    return kelvin_to_celsius(fahr_to_kelvin(temp_f))
print('freezing point of water in Celsius:', fahr_to_celsius(32.0))

print(outer("helium"))
hm

• This is called function scope

a = "Hello"

def my_fn(a):
  a = "Goodbye"

my_fn(a)  
print(a)

What would be printed if you ran the code below?

a, b = 3, 7

def swap(a, b):
    temp = a
    a = b
    b = temp

swap(a, b)
print(b, a)

1. 7 3
2. 3 7
3. 3 3
4. 7 7

• Open the analyse_files.py notebook from the first lesson

1. check the data for problems
2. plot the data

• We noticed that some data was questionable
• This function spots problems with the data
  • Call the function after loading, before plotting

def detect_problems(data):
    if numpy.max(data, axis=0)[0] == 0 and numpy.max(data, axis=0)[20] == 20:
        print("Suspicious-looking maxima!")
    elif numpy.sum(numpy.min(data, axis=0)) == 0:
        print("Minima sum to zero!")
    else:
        print("Seems OK!")

def plot_data(data, fname):
    # create figure and three axes
    fig = plt.figure(figsize=(10.0, 3.0))
    [...]

• Loop over the files, load data, detect_problems() and plot_data()

# Analyse each file in turn
for fname in files:
    print("Analysing", fname)

    # load data
    data = np.loadtxt(fname=fname, delimiter=',')

    # identify problems in the data
    detect_problems(data)

    # plot image in file
    imgname = fname[:-4] + '.png'
    plot_data(data, imgname)        

• After 6 months, the referee’s report arrives and you need to rerun experiments
• Another student is continuing the project
• Some random person reads your article and asks for the code
• Helps spot errors quickly
• Clarifies structure in your mind as well as in the code
• Saves you time in the long run! (“Future You” will back this up)

• Use the values of command-line arguments in a program.
• Handle flags and files separately in a command-line program.
• Read data from standard input in a program so that it can be used in a pipeline (with pipes: |)

$ nano sys_version.py

import sys
print('version is', sys.version)

$ python sys_version.py 

$ nano sys_argv.py

import sys
print('sys.argv is', sys.argv)

$ python sys_argv.py 
$ python sys_argv.py item1 item2 somefile.txt

$ python readings.py mydata.csv

• We will make it take options --min, --max, --mean
  • The script will report one of these

$ python readings.py --min mydata.csv

• We will make it handle multiple files

$ python readings.py --min mydata.csv myotherdata.csv

• We will make it take STDIN so we can use it with pipes

$ cat mydata.csv | readings.py --min

• We’ll build features in one-by-one

$ nano readings.py

import sys
import numpy

def main():
    script = sys.argv[0]
    filename = sys.argv[1]
    data = numpy.loadtxt(filename, delimiter=',')
    for m in numpy.mean(data, axis=1):
        print(m)

• If we do this, we can use the same file as:
  • a module (import readings)
  • a script ($ python readings.py)
• The Python code has __name__ == '__main__' only when run as a script

if __name__ == '__main__':
   main()

$ ls data/small-*
data/small-01.csv  data/small-02.csv  data/small-03.csv
$ python sys_argv.py data/small-*

• All arguments from index 1 onwards are filenames

def main():
    script = sys.argv[0]
    for filename in sys.argv[1:]:
        print(filename)
        data = numpy.loadtxt(filename, delimiter=',')
        for m in numpy.mean(data, axis=1):
            print(m)

$ python readings.py --max myfile.csv

• We should check that flags are valid

def main():
    script = sys.argv[0]
    action = sys.argv[1]
    filenames = sys.argv[2:]
    if action not in ['--min', '--mean', '--max']:
        print('Action is not one of --min, --mean, or --max: ' + action)
        sys.exit(1)
    for f in filenames:
        process(f, action)

• The process() function returns the summarised data

def process(filename, action):
    data = numpy.loadtxt(filename, delimiter=',')

    if action == '--min':
        values = numpy.min(data, axis=1)
    elif action == '--mean':
        values = numpy.mean(data, axis=1)
    elif action == '--max':
        values = numpy.max(data, axis=1)

    for m in values:
        print(m)

• sys.stdin catches STDIN from the operating system

    if len(filenames) == 0:
        process(sys.stdin, action)
    else:
        for f in filenames:
            process(f, action)

$ python readings.py --max < data/small-01.csv

• Once written, functions are reused

• Example: centring a numerical array

• Add the function

import numpy as np

def centre(data, value):
    return (data - np.mean(data)) + value

• We could try centre() on real data
  • but we don’t know the answer!

z = np.zeros((2, 2))
print(centre(z, 3.0))

data = np.loadtxt(fname='data/inflammation-01.csv', delimiter=',')
print(centre(data, 0))

• But how do we know it worked?

• We can check properties of the original and centred data
  • mean, min, max, std

centred = centre(data, 0)
print('original min, mean, and max are:', 
      np.min(data), np.mean(data), np.max(data))
print('min, mean, and max of centered data are:', 
      np.min(centred), np.mean(centred), np.max(centred))
print('std dev before and after:', 
      np.std(data), np.std(centred))      

• Writing comments in the code (using the hash #) is a good thing

def centre(data, value):
    """Returns the array in data, recentered around value."""
    return (data - np.mean(data)) + value

help(centre)

def centre(data, value=0.0):
    """Returns the array in data, recentered around the value"""
    return (data - np.mean(data)) + value

centre(data, 0.0)
centre(data, desired=0.0)
centre(data)

• All programmers make errors
• Identifying, fixing, and coping with errors is a valuable skill

def favourite_ice_cream():
    ice_creams = [
        "chocolate",
        "vanilla",
        "strawberry"
    ]
    print(ice_creams[3])

favourite_ice_cream()

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-2-713948c7cbba> in <module>
----> 1 favourite_ice_cream()

<ipython-input-1-d5749f56e901> in favourite_ice_cream()
      1 def favourite_ice_cream():
      2     ice_creams = ["chocolate", "vanilla", "strawberry"]
----> 3     print(ice_creams[3])

IndexError: list index out of range

• (mostly, you can just look at the last couple of levels)

def some_function()
    msg = "hello, world!"
    print(msg)
     return msg

  File "<ipython-input-3-dbf32ad5d3e8>", line 1
    def some_function()
                       ^
SyntaxError: invalid syntax

def some_function():
    msg = "hello, world!"
    print(msg)
     return msg

  File "<ipython-input-4-e169556d667b>", line 4
    return msg
    ^
IndentationError: unexpected indent

print(a)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-5-c5a4f3535135> in <module>()
----> 1 print(a)

NameError: name 'a' is not defined

letters = ['a', 'b']
[...]
print("Letter #3 is", letters[2])

Letter #1 is a
Letter #2 is b
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-9-62bced7460d2> in <module>()
      2 print("Letter #1 is", letters[0])
      3 print("Letter #2 is", letters[1])
----> 4 print("Letter #3 is", letters[2])

IndexError: list index out of range

• Can you read the code below, and (without running it) identify what the errors are?
• Can you fix all the errors so the code prints abbabbabba?

for number in range(10):
    # use a if the number is a multiple of 3, otherwise use b
    if (Number % 3) = 0:
        message = message + a
    else:
        message = message + "b"
print(message)

def s(p):
    a = 0
    for v in p:
        a += v
    m = a / len(p)
    d = 0
    for v in p:
        d += (v - m) * (v - m)
    return numpy.sqrt(d / (len(p) - 1))

def std_dev(sample):
    sample_sum = 0
    for value in sample:
        sample_sum += value

    sample_mean = sample_sum / len(sample)

    sum_squared_devs = 0
    for value in sample:
        sum_squared_devs += (value - sample_mean) * (value - sample_mean)

    return numpy.sqrt(sum_squared_devs / (len(sample) - 1))

• We’ve not focused on whether the code is ‘correct’
• Defensive programming is expecting your code to have mistakes, and guarding against them

• This is good practice
  • speeds up software development
  • helps ensure that your code is doing what you intend

• We assert that a condition is True
  • If it’s True, the code may be correct
  • If it’s False, the code is not correct

assert <condition>, "Some text describing the problem"

numbers = [1.5, 2.3, 0.7, -0.001, 4.4]
total = 0.0
for n in numbers:
    assert n > 0.0, 'Data should only contain positive values'
    total += n
print('total is:', total)

def normalise_rectangle(rect):
    """Normalises a rectangle to the origin, longest axis 1.0 units."""
    x0, y0, x1, y1 = rect
    dx = x1 - x0
    dy = y1 - y0
    
    if dx > dy:
        scaled = float(dy) / dx
        upper_x, upper_y = 1.0, scaled
    else:
        scaled = float(dx) / dy
        upper_x, upper_y = scaled, 1.0
        
    return (0, 0, upper_x, upper_y)

• Here, we want to ensure that rect has four values

def normalise_rectangle(rect):
    """Normalises a rectangle to the origin, longest axis 1.0 units."""
    assert len(rect) == 4, "Rectangle must have four co-ordinates"
    x0, y0, x1, y1 = rect
    [...]

• Here, we want to assert that the upper x and y values are in the range [0, 1]

def normalise_rectangle(rect):
    """Normalises a rectangle to the origin, longest axis 1.0 units."""
    [...]

    assert 0 < upper_x <= 1.0, "Calculated upper x-coordinate invalid"
    assert 0 < upper_y <= 1.0, "Calculated upper y-coordinate invalid"    
        
    return (0, 0, upper_x, upper_y)

• assertions declare what the program should be doing
• assertions help the person reading the program match their understanding of the code to what the code expects

• Turn bugs into assertions or tests: if you’ve made the mistake once, you might make it again

• variables
• data types: numpy.arrays, lists, strings, numbers
• file IO: loading data, listing files, manipulating filenames
• computing statistics
• plotting data: plots and subplots
• program flow: loops and conditionals
• automating multiple analyses
• Python scripts: edit-save-execute

Please do use the course Etherpad:

• http://pad.software-carpentry.org/2019-05-02-standrews

• To learn basic concepts of programming (in Python)