This tutorial was generated from an IPython notebook. You can download the notebook here.
In Tutorial 0, we installed Enthought Canopy. Canopy contains most of what we need to do scientific computing with Python. At the most basic level, it has Python 2.7. It contains other modules we will make heavy use of, the three most important ones being NumPy, matplotlib, and IPython. We will also make heavy use of pandas, SciPy, and scikit-image throughout the course.
In this tutorial, we will first learn some of the basics of the Python programming language at the same time exploring the properties of NumPy's very useful (and as we'll see, ubiquitous) ndarray data structure. Finally, we'll load some data and use matplotlib to generate plots.
Hello, world.¶Normally we would start with a Hello, world program and then move systematically forward, performing a series of simple tasks in Python. This tutorial will plow through that and go directly toward using NumPy. You will pick up Python's syntax as we go along.
Even though this is the plan, we have to start with Hello, world. So, fire up an IPython session using Canopy or whatever else you like and enter the following at the prompt.
print('Hello, world.')
We see the syntax for function calls in Python. Function arguments are enclosed in parentheses. Ok, now to modules and NumPy.
It is common that scientific software packages such as Matlab and Mathematica are optimized for a flavor of scientific computing (such as matrix computation in the case of Matlab) and are rather full-featured. On the other hand, Python is a programming language. It was not specifically designed to do scientific computing. So, plain old Python is very limited in scientific computing capability.
However, Python is very flexible and allows use of modules. A module contains classes, functions, attributes, data types, etc., beyond what is built in to Python. In order to use a module, you need to import it to make it available for use. So, as we begin working on data analysis, we need to import modules we will use.
The first things we will import come from the __future__ module. This is a special module that enables use of Python 3 standards while running Python 2. Having these in place in your code will help you in the future when you eventually migrate to Python 3. Before we do that, we'll do an experiment in Python 2 division.
# The print function is smart: it will display numerical results.
print(4 / 3)
print(4.0 / 3.0)
Notice that division of two integers in Python 2 results in the integer floor of the result. This can be a problem if you are not careful! In Python 3, this is not the case. So, we have to be sure to import division from the __future__ module. We also import the other Python 3 compatibility modules.
from __future__ import division, absolute_import, \
print_function, unicode_literals
The construction from <module> import <attributes> puts <attributes> (the things you imported) in the namespace. This construction allows you to pick and choose what attributes you want to import from a given module. Let's try division again.
# Python 3 divides integers to return a float
print(4 / 3)
# If you want to return an integer (floor of result), use //
print(4 // 3)
That's better!
It is important to note that until we imported the __future__ module, its capacities were not available. Keep that in mind: Plain old Python won't do much until you import a module.
Let's now import one of the major workhorses of our class, NumPy!
# Importing is done with the import statement
import numpy as np
# We now have access to some handy things
print('circumference / diameter = ', np.pi)
print('cos(pi) = ', np.cos(np.pi))
Notice that we used the import ... as construction. This enabled us to abbreviate the name of the module so we do not have to type numpy each time.
Also, notice that to access the (approximate) value of $\pi$ in the numpy module, we prefaced the name of the attiribute (pi) with the module name followed by a dot (np.). This is generally how you access attributes in modules.
We also learned an important bit of syntax: strings are enclosed in single (or double) quotes.
We're already getting dangerous with Python. So dangerous, in fact, that we'll write our own module!
Modules are stored in files ending in .py. As an example, we will create a module that finds the roots of the quadratic equation
\begin{align} ax^2 + bx + c = 0. \end{align}
Using the Canopy editing window (or your favorite text editor), we'll create a file called quadratic.py containing the code below. The file should be saved in a directory that is part of your PYTHONPATH environment variable (which usually contains the present working directory) so that the interpreter can find it.
"""
*** This should be stored in a file quadratic.py. ***
Quadratic formula module
"""
from __future__ import division, print_function
import numpy as np
# ############
def discriminant(a, b, c):
"""
Returns the discriminant of a quadratic polynomial
a * x**2 + b * x + c = 0.
"""
return b**2 - 4.0 * a * c
# ############
def roots(a, b, c):
"""
Returns the roots of the quadratic equation
a * x**2 + b * x + c = 0.
"""
delta = discriminant(a, b, c)
root_1 = (-b + np.sqrt(delta)) / (2.0 * a)
root_2 = (-b - np.sqrt(delta)) / (2.0 * a)
return root_1, root_2
There is a whole bunch of syntax in there to point out.
.py file.def statement. It has the function prototype, followed by a colon.def statement is part of the function. Once the indentation goes back to the level of the def statement, you are no longer in the function..py file).**.Now, let's test our new module out!
import quadratic as qd
# Python has nifty syntax for making multiple definitions on the same line
a, b, c = 3.0, -7.0, -6.0
# Call the function and print the result
root_1, root_2 = qd.roots(a, b, c)
print('roots:', root_1, root_2)
Very nice! Now, let's try another example. This one might have a problem....
# Specify a, b, and c that will give imaginary roots
a, b, c = 1.0, -2.0, 2.0
# Call the function and print the result
root_1, root_2 = qd.roots(a, b, c)
print('roots:', root_1, root_2)
Oh no! It gave us nan, which means "not a number," as our roots. It also gave some warning that it encountered invalid (negative) arguments for the sqrt function. The roots should be $1 \pm i$, where $i = \sqrt{-1}$. We will use this opportunity to introduce Python's control flow, starting with an if statement.
if statement¶We will decree that our quadratic equation solver only handles real roots, so it will raise an exception if an imaginary root is encountered. So, we modify the contents of the file quadratic.py as follows.
"""
*** This should be stored in a file quadratic.py. ***
Quadratic formula module
"""
from __future__ import division, print_function
import numpy as np
# ############
def discriminant(a, b, c):
"""
Returns the discriminant of a quadratic polynomial
a * x**2 + b * x + c = 0.
"""
return b**2 - 4.0 * a * c
# ############
def roots(a, b, c):
"""
Returns the roots of the quadratic equation
a * x**2 + b * x + c = 0.
"""
delta = discriminant(a, b, c)
if delta < 0.0:
raise ValueError('Imaginary roots! We only do real roots!')
else:
root_1 = (-b + np.sqrt(delta)) / (2.0 * a)
root_2 = (-b - np.sqrt(delta)) / (2.0 * a)
return root_1, root_2
We have now exposed the syntax for a Python if statement. The conditional expression ends with a colon, just like the def statement. Note the indentation of blocks of code after the conditionals. (We actually did not need the else statement, because the program would just continue without the exception, but I left it there for illustrative purposes. It is actually preferred not to have the else statement.)
Now if we re-import the module (we can use the Python function reload for this), the if statement will catch our imaginary roots and raise an exception. Note that you must reload (or start over again and import) the module before your changes take effect.
# Reload the quadratic module using its abbeviated name we already defined
reload(qd)
# Pass in parameters that will give imaginary roots
a, b, c = 1.0, -2.0, 2.0
root_1, root_2 = qd.roots(a, b, c)
This threw the appropriate exception.
Congrats! You wrote a functioning module. But now is an important lesson....
If you are trying to do a task that you think might be common, it's probably part of NumPy or some other package. Look, or ask Google, first. In this case, NumPy has a function called roots that computes the roots of a polynomial. To figure out how to use it, we can either look at the doc string, or look in the NumPy and SciPy documentation online (the documentation for np.roots is available here). To look at the doc string, you can enter the following at an IPython prompt:
np.roots?
We see that we need to pass the coefficients of the polynomial we would like the roots of using an "array_like" object. We will discuss what this means in a moment, but for now, we will just use a list to specify our coefficients and call the np.roots function.
# Define the coefficients in a list (using square brackets)
coeffs = [3.0, -7.0, -6.0]
# Call np.roots. It returns an np.ndarray with the roots
roots = np.roots(coeffs)
print('Roots for (a, b, c) = (3, -7, -6):', roots)
# It even handles complex roots!
roots = np.roots([1.0, -2.0, 2.0])
print('Roots for (a, b, c) = (1, -2, 2): ', roots)
array_like data types¶In the previous example, we used a list as an array_like data type. Python has several native data types. We have already mentioned ints and floats. We just were not very explicit about it. Python's native array_like data types are lists and tuples. There are collections of items (in this class, usually numbers) separated by commas. They can be nested to generate lists of lists, lists of tuples, tuples of tuples, etc. The primary difference between a list and a tuple is that a list is mutable, meaning that it can be changed in-place, where a tuple is immutable.
# This is a list
my_list = [1, 2, 3, 4]
# This is a tuple
my_tuple = (1, 2, 3, 4)
# We can change a list in place
my_list[1] = 1.3
print(my_list)
# We cannot change a tuple in place. Below will raise an exception
# my_tuple[1] = 1.3
Notice how we indexed the list. The index uses square brackets. Importantly, indexing in Python starts with zero.
np.ndarray: maybe your new best friend¶Lists and tuples can be useful, but for many many applications in data analysis, the np.ndarray, which we will colloquially call a "NumPy array," is most often used. They are created using the np.array function with a list or tuple as an argument. Once created, we can do all sorts of things with them. Let's play!
a = np.array([1.0, 2.0, 3.0])
b = np.array([4.0, 5.0, 6.0])
# Arithmetic operations are done elementwise
print('a: ', a)
print('a + b: ', a + b)
print('a * b: ', a * b)
print('1.0 + a:', 1.0 + a)
print('a**2: ', a**2)
print('b**a: ', b**a)
Let's generate some longer arrays to play with. We will use the function np.linspace. Check out its doc string to see how it works.
# Make 20 evenly spaced points from 0 to 2*pi
x = np.linspace(0.0, 2.0 * np.pi, 100)
# NumPy functions also work elementwise. Let's make a function exp(sin(x))
y = np.exp(np.sin(x))
# Let's look at them
print('x: ', x)
print('y: ', y)
That's not very useful. It would be nice to plot the function. So, now is a very nice time to introduce matplotlib!
matplotlib: our primary plotting tool¶matplotlib is a tool for plotting the data stored in NumPy arrays. We will mostly use the interactive plotting module, matplotlib.pyplot, which we will abberviate as plt. Its syntax is quite simple, and best learned by example. We will plot our function and look at the graphics.
# Import Matplotlib
import matplotlib.pyplot as plt
# For the purposes of generating this tutorial, I use inline matplotlib.
# You do not need to do this, but this line of code will be in all tutorials.
%matplotlib inline
# Make the plot!
plt.plot(x, y)
# Label the axes. The r'' construction makes the string a literal, which
# enables use of $, \, etc. matplotlib then interprets TeX!
plt.xlabel(r'$x$')
plt.ylabel(r'$\mathrm{e}^{\sin x}$')
# We would like the limits in the x axis to be 0 to 2*pi
plt.xlim((0.0, 2.0 * np.pi)) # Limits passes as a tuple argument
This is all very exciting! Next, we will numerically compute the derivative of the function we just plotted. We know the analytical expression.
\begin{align} \frac{\mathrm{d}y}{\mathrm{d}x} = \mathrm{e}^{\sin x}\,\cos x = y \cos x, \end{align}
but we will also compute it numerically. We will do it using forward differencing.
\begin{align} \frac{\mathrm{d}y(x_i)}{\mathrm{d}x} \approx \frac{y(x_{i+1}) - y(x_i)}{x_{i+1} - x_i}. \end{align}
We will use this opportunity to look at another aspect of Python's control flow, the for loop.
for loop¶Use of a for loop is best seen through example. We will define a function to do finite difference differentiation.
# Define forward differencing function
def forward_diff(y, x):
# Use np.empty to make an empty ndarray to put our derivatives in
deriv = np.empty(len(y) - 1)
# Use a for loop to go through each point and compute forw. diff. deriv.
for i in range(len(y)-1):
deriv[i] = (y[i+1] - y[i]) / (x[i+1] - x[i])
# Return the derivative (remember, it is a NumPy array)
return deriv
# Call the function to perform finite differencing
deriv = forward_diff(y, x)
Let's go over some of the functions we just used there. The len function returns the length of a one-dimensional array. In our case, len(y) is 100, since y has 100 elements.
The range function is built in to Python. It generates intgers that we can use to iterate for loops. In this case, it generates all integers from 0 to 98. That's right, even though len(y) - 1 is 99, the range function does not generate the last number. Don't get burned by this!
We used the for loop to successively compute the elements of the array. We took care to note that there are only 99 difference of points and 100 total points.
Now, we can plot our derivatives. Let's plot them right on top of the analytical expression.
# Compute exact derivative
deriv_exact = y * np.cos(x)
# Plot exact derivative as a black line (given by 'k-').
plt.plot(x, deriv_exact, 'k-')
# Plot approximate derivative with red x's; derive is approximate between
# grid points in x
plt.plot((x[1:] + x[:-1]) / 2.0, deriv, 'rx')
# Label axes, set axis limits
plt.xlabel(r'$x$')
plt.ylabel(r'$\mathrm{d}y/\mathrm{d}x$')
plt.xlim((0.0, 2.0 * np.pi))
# Put in a legend
plt.legend(('exact', 'fin. diff'), loc='upper center', numpoints=1)
In the code above, notice how we included a legend in the plot. The first argument is a tuple containing the names of the labels for the two curves. The following arguments are keyword arguments. Keyword arguments do not need to be specified in the function call. If left unspecified, default values are used for them. They are specified with the syntax above, as kwarg=value within the function call.
I also want to draw your attention to the call to plt.plot for the derivative:
plt.plot((x[1:] + x[:-1]) / 2.0, deriv, 'rx')
Remember that the derivative is only defined for 99 of the 100 points we sample along the $x$-axis. We therefore sliced the x NumPy array. The colon implies filling in all indices. Python allows negative indices, and -1 means the last index of the array. Like the range function, slicing does not include the last index. So, x[:-1] means that we take all values of the array except the last one. Similarly, x[1:] means that we take all values of the array except the first one.
You might think finite differencing is something that one would commonly want to do, right? So, of course, NumPy has a function for it!
# Use np.diff function to compute forward differences
np_deriv = np.diff(y) / np.diff(x)
# Verify that the results is the same we got
print('Did NumPy give what we got?:', np.allclose(np_deriv, deriv))
That code is much cleaner than our clunky for loop. It is also much faster. We can use IPython %timeit magic function to test.
%timeit np_deriv = np.diff(y) / np.diff(x)
%timeit deriv = forward_diff(y, x)
We will now integrate our function $y(x)$. From the integral formulas for modified Bessel functions, we know
\begin{align} \int_0^{2\pi}\mathrm{d} x\, \mathrm{e}^{\sin x} = 2\pi \,I_0(1), \end{align}
where $I_0$ is the modified Bessel function of the first kind. Fortunately, SciPy has a module that allows for calculation of special functions, this particular Bessel function included.
import scipy.special
# Call scipy.special.iv to compute the value of the integral
exact_integral = 2.0 * np.pi * scipy.special.iv(0, 1.0)
print('Exact integral:', exact_integral)
Now, let's approximately compute the integral by the trapezoidal rule. We could write a for loop to do this. An easier way to do this integral is to note that the integral is over the periodic domain of the function. If we were to write the function as a Fourier series, it consists of a constant term, plus a bunch of terms with periods of multiples of $2\pi$. So, only the constant term, which is the mean of the function over the interval, remains. All of the added periodic terms integrate to zero. The integral is then just like computing the area of a rectangle: $2\pi \langle y(x)\rangle_x$.
\begin{align} \int_0^{2\pi}\mathrm{d} x\, y(x) \approx \frac{2\pi}{n}\sum_{i=0}^{n} y(x_i), \end{align}
where we have sampled $y$ at $n$ evenly spaced points over the interval of length $2\pi$. We can then use the mean method of NumPy arrays to quickly carry out the sum.
# Compute integral numerically (do not use last point because it would
# then appear twice in the mean)
approx_integral = 2.0 * np.pi * y[:-1].mean()
print('Approx integral:', approx_integral)
print('exact - approx: ', exact_integral - approx_integral)
It is exact to machine precision! This is a feature known as spectral accuracy, which we will not cover in this class, but is fascinating nonetheless.
Note how we used the construction y[:-1].mean(). y is a NumPy array, so it has associated with it many methods. A method is a function that can operate on an instance of a class. y is an instance of a NumPy array, so y.mean() computes the mean of the NumPy array y. It is advisable to use these native methods on NumPy arrays. To prove this point, we will compare use of Python's built-in sum function to compute a sum, compared to a computation using NumPy. We will also compare it to summing using a for loop.
# Make a sum function that uses a for loop
def sum_with_for_loop(y):
my_sum = 0.0
for i in range(len(y)):
my_sum += y[i] # += adds the right hand side to the left hand side
return my_sum
# Time the different ways of evaluating the sum the entries of a NumPy array
%timeit sum_with_for_loop(y)
%timeit sum(y)
%timeit np.sum(y)
%timeit y.sum()
PEPs (Python Enhancement Proposals) document suggestions and guidelines for the Python language, its development, etc. My personal favorite is PEP8, the Style Guide for Python Code. This was largely written by Guido von Rossum, the inventor of Python and its benevolent dictator for life. It details how you should style your code. As Guido says, code is much more often read than written. I strongly urge you to follow PEP8 the best you can. It's a lot to read, so I will highlight important points here. If you follow these guidelines, your code will be much more readable than otherwise. This is particularly useful because you are working in groups and the TAs need to grade your work.
\).**. E.g., a = b + c * d**2.my_fun(a, b, c=True).a[i], not a [ i ].This concludes our introductory tour of Python with some NumPy, SciPy, and matplotlib thrown in for good measure. There is still much to learn, but you will pick up more and more tricks and syntax as we go along.
As we work through the course, it is important to remember that Python with NumPy, SciPy, matplotlib, pandas, and the other packages we will use throughout the term, are just tools we are using to analyze data. They not the best tool for each type of data set. However, the Python-based tools constitute a nice Swiss Army knife. They can almost always get the job done.
Most importantly, we are using Python to learn data analysis concepts and methods. This is not a class about Python; it is a class about data analysis in the biological sciences.