We are studying inflammation in patients who have been given a new treatment for arthritis, and need to analyze the first dozen data sets. The data sets are stored in comma-separated values (CSV) format. Each row holds data for a single patient. The columns represent arthritis inflammations per day on successive days. The first few rows of our first file look like this:

0,0,1,3,1,2,4,7,8,3,3,3,10,5,7,4,7,7,12,18,6,13,11,11,7,7,4,6,8,8,4,4,5,7,3,4,2,3,0,0
0,1,2,1,2,1,3,2,2,6,10,11,5,9,4,4,7,16,8,6,18,4,12,5,12,7,11,5,11,3,3,5,4,4,5,5,1,1,0,1
0,1,1,3,3,2,6,2,5,9,5,7,4,5,4,15,5,11,9,10,19,14,12,17,7,12,11,7,4,2,10,5,4,2,2,3,2,2,1,1
0,0,2,0,4,2,2,1,6,7,10,7,9,13,8,8,15,10,10,7,17,4,4,7,6,15,6,4,9,11,3,5,6,3,3,4,2,3,2,1
0,1,1,3,3,1,3,5,2,4,4,7,6,5,3,10,8,10,6,17,9,14,9,7,13,9,12,6,7,7,9,6,3,2,2,4,2,0,1,1

The master copy of the data for this study is in a git repository. You can browse that repository here: https://github.com/SWC-OSG-Workshop/ExampleData. Let's take a look. Inside the Python folder are several CSV files. We'll use them in the upcoming examples.

N.B. You can also get a copy of this data in your OSG Connect account by cloning it from github:

 
$ git clone https://github.com/SWC-OSG-Workshop/ExampleData.git

We want to:

• load the CSV data into memory,
• calculate the average inflammation per day across all patients, and
• plot the result.

To do all that, we'll learn a little bit about programming in Python.

The essence of practical programming — as with mathematics — is learning to use fundamental tools and techniques to accomplish more complex tasks, without feeling overwhelmed by the magnitude of the greater task.

Many of powerful tools are built into the core of languages like Python, but as with a volume of mathematical proofs, even more power resides in the libraries that developers build. Different tasks that these libraries implement can be combined sequentially or conditionally to produce programs that do a lot of work, but with relatively simple rules.

In Python, a library is loaded using the import statement. Let's do a simple demonstration by looking for the files in the current folder. (Recall from the section on shell programming that . represents the current directory.) Type the following:

import os
os.listdir('.')

os is a library that implements some useful procedures for interacting with the operating system. listdir is a function, residing in the os library (or "module") which finds all the files in a directory. It returns these as a list that you can assign to a variable, or just display on the screen.

In order to work with our inflammation data, we will import a library called NumPy that knows how to operate on matrices:

import numpy

numpy.loadtxt(fname='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.]])

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

weight_kg = 55

print weight_kg

55

print 'weight in pounds:', 2.2 * weight_kg

weight in pounds: 121.0

weight_kg = 57.5
print 'weight in kilograms is now:', weight_kg

weight in kilograms is now: 57.5

In Python, as in most other common programming languages, assignment refers to the linking of a name (also called a variable, symbol, or label) to a value. It's a simple one-way relation. Assigning a value to one variable does not change the values of other variables.

For example, let's now store the subject's weight in pounds in a variable:

weight_lb = 2.2 * weight_kg
print 'weight in kilograms:', weight_kg, 'and in pounds:', weight_lb

weight in kilograms: 57.5 and in pounds: 126.5

weight_kg = 100.0
print 'weight in kilograms is now:', weight_kg, 'and weight in pounds is still:', weight_lb

weight in kilograms is now: 100.0 and weight in pounds is still: 126.5

Even though we initially defined weight_lb in terms of weight_kg, that relation does not persist, and weight_lb doesn't change when we redefine weight_kg. The conversion of weight_kg to weight_lb happens at the moment that it's assigned, and never again after. The assignment itself is binding. This may feel unnatural, especially if you're used to working with spreadsheets. In a spreadsheet the only way to define one cell in terms of another is through a function. Python has functions too, but they don't look quite the same as assignments. We'll go over these later.

Now that we know how to assign things to variables, let's re-run numpy.loadtxt and save its result:

This statement doesn't produce any output because assignment doesn't return anything — the value previously returned and displayed is now bound up in the variable data.. If we want to check that our data has been loaded, we can print the variable's value:

Now that our data is in memory, we can start doing things with it. First, let's ask what type of thing data refers to:

The output tells us that data currently refers to an N-dimensional array created by the NumPy library. We can see what its shape is like this:

This tells us that data has 60 rows and 40 columns. data.shape is a member of data, i.e., a value that is stored as part of a larger value. We use the same dotted notation for the members of values that we use for the functions in libraries because they have the same part-and-whole relationship.

If we want to get a single value from the matrix, we must provide an index in square brackets, just as we do in math:

The expression data[30, 20] may not surprise you, but data[0, 0] might. Programming languages like Fortran and MATLAB start counting at 1, because that's what humans are most accustomed to. Languages in the C family (including C++, Java, Perl, and Python) count from 0 because it's simpler for computers to do, which in turn makes it easier for programmers to sequence multiple steps that deal with the same arrays. (There's a lot less adding or subtracting 1 overall.) As a result, if we have an M×N array in Python, its indices go from 0 to M-1 on the first axis and 0 to N-1 on the second. It takes a bit of getting used to, but one way to remember the rule is that the index is how many steps we have to take from the start to get the item we want.

In the Corner

What may also surprise you is that when Python displays an array, it shows the element with index [0, 0] in the upper left corner rather than the lower left. This is consistent with the way mathematicians draw matrices, but different from the Cartesian coordinates. The indices are (row, column) instead of (column, row) for the same reason.

An index like [30, 20] selects a single element of an array, but we can select whole sections as well, using an extended form of range notation. For example, we can select the first ten days (columns) of values for the first four (rows) patients like this:

The slice 0:4 means, "Start at index 0 and go up to, but not including, index 4." Again, the up-to-but-not-including takes a bit of getting used to, but the rule is that the difference between the upper and lower bounds is the number of values in the slice.

We don't have to start slices at 0:

and we don't have to take all the values in the slice---if we provide a stride, Python takes values spaced that far apart:

Here, we have taken rows 0, 3, 6, and 9, and columns 0, 2, 4, 6, and 8. (Again, we always include the lower bound, but stop when we reach or cross the upper bound.)

We also don't have to include the upper and lower bound on the slice. If we don't include the lower bound, Python uses 0 by default; if we don't include the upper, the slice runs to the end of the axis, and if we don't include either (i.e., if we just use ':' on its own), the slice includes everything:

Arrays also know how to perform common mathematical operations on their values. If we want to find the average inflammation for all patients on all days, for example, we can just ask the array for its mean value.

mean is a method of the array, i.e., a function that belongs to it in the same way that the member shape does. If variables are nouns, methods are verbs: they are what the thing in question knows how to do. This is why data.shape doesn't need to be called (it's just a thing) but data.mean() does (it's an action). It is also why we need empty parentheses for data.mean(): even when we're not passing in any parameters, parentheses are how we tell Python to go and do something for us.

NumPy arrays have lots of useful methods:

When analyzing data, though, we often want to look at partial statistics, such as the maximum value per patient or the average value per day. One way to do this is to select the data we want to create a new temporary array, then ask it to do the calculation:

We don't actually need to store the row in a variable of its own. Instead, we can combine the selection and the method call:

Note that this is effectively the same mechanism as shell pipelines. You can call methods on top of methods so long as those methods return valid Python objects, without needing to assign these intermediary objects to variables.

What if we need the maximum inflammation for all patients, or the average for each day? As the diagram below shows, we want to perform the operation across an axis:

To support this, most array methods allow us to specify the axis we want to work on. If we ask for the average across axis 0, we get:

As a quick check, we can ask this array what its shape is:

The expression (40,) tells us we have an N×1 vector, so this is the average inflammation per day for all patients. If we average across axis 1, we get:

which is the average inflammation per patient across all days.

The mathematician Richard Hamming once said, "The purpose of computing is insight, not numbers," and the best way to develop insight is often to visualize data. Visualization deserves an entire lecture (or course) of its own, but we can explore a few features of Python's matplotlib here. First, let's tell the IPython Notebook that we want our plots displayed inline, rather than in a separate viewing window:

The % at the start of the line signals that this is a command for the notebook, rather than a statement in Python. Next, we will import the pyplot module from matplotlib and use two of its functions to create and display a heat map of our data:

Blue regions in this heat map are low values, while red shows high values. As we can see, inflammation rises and falls over a 40-day period. Let's take a look at the average inflammation over time:

Here, we have put the average per day across all patients in the variable ave_inflammation, then asked pyplot to create and display a line graph of those values. The result is roughly a linear rise and fall, which is suspicious: based on other studies, we expect a sharper rise and slower fall. Let's have a look at two other statistics:

The maximum value rises and falls perfectly smoothly, while the minimum seems to be a step function. Neither result seems particularly likely, so either there's a mistake in our calculations or something is wrong with our data.

It's very common to create an alias for a library when importing it in order to reduce the amount of typing we have to do. Here are our three plots side by side using aliases for numpy and pyplot:

The first two lines re-load our libraries as np and plt, which are the aliases most Python programmers use. The call to loadtxt reads our data, and the rest of the program tells the plotting library how large we want the figure to be, that we're creating three sub-plots, what to draw for each one, and that we want a tight layout. (Perversely, if we leave out that call to plt.tight_layout(), the graphs will actually be squeezed together more closely.)

Next Steps

Our work so far has convinced us that something's wrong with our first data file. We would like to check the other 11 the same way, but typing in the same commands repeatedly is tedious and error-prone. Since computers don't get bored (that we know of), we should create a way to do a complete analysis with a single command, and then figure out how to repeat that step once for each file. These operations are the subjects of the next two lessons.