(c) 2017 Justin Bois. This work is licensed under a Creative Commons Attribution License CC-BY 4.0. All code contained herein is licensed under an MIT license.
This tutorial was generated from an Jupyter notebook. You can download the notebook here.
In Tutorial 0, we installed the Anaconda Python distribution. Anaconda contains most of what we need to do scientific computing with Python. At the most basic level, it has Python 3.6, which we will use as the core language for our analysis. It contains other modules we will make heavy use of, the three most important ones being NumPy, Pandas, and Bokeh. We will also make heavy use of SciPy and scikit-image throughout the course. We will use the conda package manager to keep our software updated and to install new packages we may need. For example, you already installed HoloViews using conda in Tutorial 0. For some packages, we will use pip. We will be doing Markov chain Monte Carlo later in the course, mostly using PyMC3, which you also installed in Tutorial 0.
Assuming you have completed Tutorial 0 and all necessary software is installed, we will proceed to learning how to use some of the useful packages for scientific computing using Python. 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 Bokeh 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, in your first code cell, do the following.
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.
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 attribute (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 function!
As an example, we will write a function that finds the roots of the quadratic equation
\begin{align} ax^2 + bx + c = 0. \end{align}
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.
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.**.Now, let's test our new module out!
# 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 = 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 = 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 np.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.
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 roots() function as follows.
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 RuntimeError('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.)
# Pass in parameters that will give imaginary roots
a, b, c = 1.0, -2.0, 2.0
root_1, root_2 = roots(a, b, c)
This threw the appropriate exception.
Congrats! You wrote a working function. 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 np.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:
np.roots?
While not shown on a static version of this document, if you executed that cell in a live Jupyter notebook, a window pops up with documentation on how np.roots() works. 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.
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 100 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 Bokeh!
Bokeh is a data visualization tool. We will use it, and a high-level package called HoloViews, to make our plots. Important, Bokeh's plots display in browsers and are interactive. After we import the submodules we need, we execute bokeh.io.output_notebook() to tell Bokeh that the plot is to be displayed in the notebook.
# Import necessary Bokeh submodules
import bokeh.plotting
import bokeh.io
bokeh.io.output_notebook()
Now that we have Bokeh read to go, we can make our plot. We first make a figure on which to populate the glyphs of our plot. We specify the dimensions of the figure, its axis labels, and any other useful details.
# Create a figure
p = bokeh.plotting.figure(plot_height=350,
plot_width=450,
x_axis_label='x',
y_axis_label='exp(sin(x))')
Notice that no plot has yet appeared. It will not appear until we call
bokeh.io.show(p)
But, we need to finish constructing the figure first, by populating it with glyphs. Let's say we wanted to plot our curve as a line. In this case, we add a line to the plot using the p.line() method.
# Put a line on the figure
p.line(x, y, line_width=2)
Now notice that this function returned a GlyphRenderer, which we still have to show if we want to look at it. Notice also that I used keyword arguments (called "kwargs" for short) when I populated the glyphs. Keyword arguments do not need to be specified in a 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. (The default for line_width is 1.)
Ok, now we are ready to show our plot!
bokeh.io.show(p)
Notice the toolbar to the right of the plot. This enables you to interact with the plot by panning, box zooming, wheel zooming, and resetting. You can also save the image as a PNG by clicking on the floppy disk.
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.
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):
"""Compute derivative by forward differencing."""
# Use np.empty to make an empty array 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 integers 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.
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)
# Create figure
p = bokeh.plotting.figure(plot_height=350,
plot_width=450,
x_axis_label='x',
y_axis_label='dy/dx')
# Plot approximate derivative with gray dots; deriv is approximately between
# grid points in x
p.circle((x[1:] + x[:-1]) / 2.0,
deriv,
color='gray',
legend='finite difference')
# Plot analytical derivative
p.line(x, deriv_exact, color='black', legend='exact')
# Put legend at top center
p.legend.location = 'top_center'
# Show the plot
bokeh.io.show(p)
Notice that we added the legend kwarg to specify how the glyphs were labeled and to create a legend. Then, we specified the legend's location by setting the p.legend.location attribute.
With the interactive plot, you can use the zoom feature to see that the finite difference derivative doesn't quiiiite match the exact analytical one.
I also want to draw your attention to the call to p.circle for the derivative:
p.circle((x[1:] + x[:-1]) / 2.0,
deriv,
color='gray',
legend='finite difference')
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.
The style of the Bokeh plot might not quite be to your favorite specifications. For example, you might want a larger font size for the axis labels. Just as you set the p.legend.location parameter, so too can you set other parameters of the plot. Let's re-style the plot by setting various attributes.
# Set the axis labels to be 14 pt font
p.xaxis.axis_label_text_font_size = '14pt'
p.yaxis.axis_label_text_font_size = '14pt'
# Change the background color to be light gray
p.background_fill_color = '#EEEEEE'
# Make the grid white
p.grid.grid_line_color = '#FFFFFF'
# Make the legend a little more transparent
p.legend.background_fill_alpha = 0.5
bokeh.io.show(p)
I don't particularly like this styling, but I show it to give you a feel of the flexibility and how to set attributes of your plot. As you can see, the attributes are intuitively and descriptively named. You can explore what is available by looking at the Bokeh documentation or by tab-completion.
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()
I contend that the scientific paper of the future is interactive; the reader will be avble to interact with the plots and data. Bokeh is an excellent tool in this case. If you want to export your Bokeh plot to HTML, you can do that as follows.
bokeh.io.save(p, filename='my_pretty_plot.html', title='dy/dx')
We get a warning that we did not specify an output file a priori, but the file is generated nonetheless and we don't really have to worry about it. A string is returned with the location of the HTML file containing the plot, in this case /Users/Justin/git/bebi103_course/2017/tutorials/my_pretty_plot.html. This HTML file has full interactivity, just like the plot in the notebook.
Unfortunately, we are not yet in the future, and many journals want static graphics. In this case, you will want to export your plot as a scalable vector graphic, or SVG. You can then edit this in Inkscape, Illustrator, or whatever vector graphics editor your prefer. Alternatively, you can also convert it directly to PDF using CairoSVG, which you installed in Tutorial 0 while practicing using pip. I show these steps below.
Importantly, you need to have some extra tools installed to do the conversion.
pip install selenium
pip install cairosvg
conda install phantomjs pillow
With those installed, you can proceed to output the SVG.
import cairosvg
# Specify that p's output is SVG
p.output_backend = 'svg'
# Export to SVG
bokeh.io.export_svgs(p, 'my_pretty_plot.svg')
# Convert the SVG to PDF
cairosvg.svg2pdf(url='my_pretty_plot.svg', write_to='my_pretty_plot.pdf')
# Switch p's output back to HTML canvas, which is more performant for interactivity
# in the notebook
p.output_backend = 'canvas'
I got a warning, but nonetheless, the appropriate plot was generated.
This concludes our introductory tour of Python with some NumPy, SciPy, and Bokeh 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, Bokeh, pandas, and the other packages we will use throughout the term, are just tools we are using to analyze data. They may not be 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.