Introduction to Python for Data Analysis

by Alejandro Correa Bahnsen

version 0.1, Feb 2016

Part of the class Practical Machine Learning

This notebook is licensed under a Creative Commons Attribution-ShareAlike 3.0 Unported License. Special thanks goes to Rick Muller, Sandia National Laboratories

Why Python?

Python is the programming language of choice for many scientists to a large degree because it offers a great deal of power to analyze and model scientific data with relatively little overhead in terms of learning, installation or development time. It is a language you can pick up in a weekend, and use for the rest of one's life.

The Python Tutorial is a great place to start getting a feel for the language. To complement this material, I taught a Python Short Course years ago to a group of computational chemists during a time that I was worried the field was moving too much in the direction of using canned software rather than developing one's own methods. I wanted to focus on what working scientists needed to be more productive: parsing output of other programs, building simple models, experimenting with object oriented programming, extending the language with C, and simple GUIs.

I'm trying to do something very similar here, to cut to the chase and focus on what scientists need. In the last year or so, the Jupyter Project has put together a notebook interface that I have found incredibly valuable. A large number of people have released very good IPython Notebooks that I have taken a huge amount of pleasure reading through. Some ones that I particularly like include:

I find Jupyter notebooks an easy way both to get important work done in my everyday job, as well as to communicate what I've done, how I've done it, and why it matters to my coworkers. In the interest of putting more notebooks out into the wild for other people to use and enjoy, I thought I would try to recreate some of what I was trying to get across in the original Python Short Course, updated by 15 years of Python, Numpy, Scipy, Pandas, Matplotlib, and IPython development, as well as my own experience in using Python almost every day of this time.

Why Python for Data Analysis?

  • Python is great for scripting and applications.
  • The pandas library offers imporved library support.
  • Scraping, web APIs
  • Strong High Performance Computation support
    • Load balanceing tasks
    • MPI, GPU
    • MapReduce
  • Strong support for abstraction
    • Intel MKL
    • HDF5
  • Environment

What You Need to Install

There are two branches of current releases in Python: the older-syntax Python 2, and the newer-syntax Python 3. This schizophrenia is largely intentional: when it became clear that some non-backwards-compatible changes to the language were necessary, the Python dev-team decided to go through a five-year (or so) transition, during which the new language features would be introduced and the old language was still actively maintained, to make such a transition as easy as possible.

Nonetheless, I'm going to write these notes with Python 3 in mind, since this is the version of the language that I use in my day-to-day job, and am most comfortable with.

With this in mind, these notes assume you have a Python distribution that includes:

  • Python version 3.4;
  • Numpy, the core numerical extensions for linear algebra and multidimensional arrays;
  • Scipy, additional libraries for scientific programming;
  • Matplotlib, excellent plotting and graphing libraries;
  • IPython, with the additional libraries required for the notebook interface.
  • Pandas, Python version of R dataframe
  • Seaborn, used mainly for plot styling
  • scikit-learn, Machine learning library!

A good, easy to install option that supports Mac, Windows, and Linux, and that has all of these packages (and much more) is the Anaconda.

Checking your installation

You can run the following code to check the versions of the packages on your system:

(in IPython notebook, press shift and return together to execute the contents of a cell)


In [1]:
import sys

print('Python version:', sys.version)

import IPython
print('IPython:', IPython.__version__)

import numpy
print('numpy:', numpy.__version__)

import scipy
print('scipy:', scipy.__version__)

import matplotlib
print('matplotlib:', matplotlib.__version__)

import pandas
print('pandas:', pandas.__version__)

import sklearn
print('scikit-learn:', sklearn.__version__)

import seaborn
print('seaborn', seaborn.__version__)


Python version: 3.5.1 |Anaconda 2.4.0 (64-bit)| (default, Dec  7 2015, 11:16:01) 
[GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
IPython: 4.0.0
numpy: 1.10.2
scipy: 0.16.0
matplotlib: 1.5.1
pandas: 0.17.0
scikit-learn: 0.17
seaborn 0.6.0
/home/al/anaconda3/lib/python3.5/site-packages/matplotlib/__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.
  warnings.warn(self.msg_depr % (key, alt_key))

I. Python Overview

This is a quick introduction to Python. There are lots of other places to learn the language more thoroughly. I have collected a list of useful links, including ones to other learning resources, at the end of this notebook. If you want a little more depth, Python Tutorial is a great place to start, as is Zed Shaw's Learn Python the Hard Way.

The lessons that follow make use of the IPython notebooks. There's a good introduction to notebooks in the IPython notebook documentation that even has a nice video on how to use the notebooks. You should probably also flip through the IPython tutorial in your copious free time.

Briefly, notebooks have code cells (that are generally followed by result cells) and text cells. The text cells are the stuff that you're reading now. The code cells start with "In []:" with some number generally in the brackets. If you put your cursor in the code cell and hit Shift-Enter, the code will run in the Python interpreter and the result will print out in the output cell. You can then change things around and see whether you understand what's going on. If you need to know more, see the IPython notebook documentation or the IPython tutorial.

Using Python as a Calculator

Many of the things I used to use a calculator for, I now use Python for:


In [2]:
2+2


Out[2]:
4

In [3]:
(50-5*6)/4


Out[3]:
5.0

(If you're typing this into an IPython notebook, or otherwise using notebook file, you hit shift-Enter to evaluate a cell.)

In the last few lines, we have sped by a lot of things that we should stop for a moment and explore a little more fully. We've seen, however briefly, two different data types: integers, also known as whole numbers to the non-programming world, and floating point numbers, also known (incorrectly) as decimal numbers to the rest of the world.

We've also seen the first instance of an import statement. Python has a huge number of libraries included with the distribution. To keep things simple, most of these variables and functions are not accessible from a normal Python interactive session. Instead, you have to import the name. For example, there is a math module containing many useful functions. To access, say, the square root function, you can either first

from math import sqrt

and then


In [4]:
sqrt(81)


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-4-e1bf934cd6c7> in <module>()
----> 1 sqrt(81)

NameError: name 'sqrt' is not defined

In [5]:
from math import sqrt
sqrt(81)


Out[5]:
9.0

or you can simply import the math library itself


In [6]:
import math
math.sqrt(81)


Out[6]:
9.0

You can define variables using the equals (=) sign:


In [7]:
radius = 20
pi = math.pi
area = pi * radius ** 2 
area


Out[7]:
1256.6370614359173

If you try to access a variable that you haven't yet defined, you get an error:


In [8]:
volume


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-8-0c7fc58f9268> in <module>()
----> 1 volume

NameError: name 'volume' is not defined

and you need to define it:


In [9]:
volume = 4/3*pi*radius**3
volume


Out[9]:
33510.32163829113

You can name a variable almost anything you want. It needs to start with an alphabetical character or "_", can contain alphanumeric charcters plus underscores ("_"). Certain words, however, are reserved for the language:

and, as, assert, break, class, continue, def, del, elif, else, except, 
exec, finally, for, from, global, if, import, in, is, lambda, not, or,
pass, print, raise, return, try, while, with, yield

Trying to define a variable using one of these will result in a syntax error:


In [10]:
return = 0


  File "<ipython-input-10-2b99136d4ec6>", line 1
    return = 0
           ^
SyntaxError: invalid syntax

The Python Tutorial has more on using Python as an interactive shell. The IPython tutorial makes a nice complement to this, since IPython has a much more sophisticated iteractive shell.

Strings

Strings are lists of printable characters, and can be defined using either single quotes


In [11]:
'Hello, World!'


Out[11]:
'Hello, World!'

or double quotes


In [12]:
"Hello, World!"


Out[12]:
'Hello, World!'

But not both at the same time, unless you want one of the symbols to be part of the string.


In [13]:
"He's a Rebel"


Out[13]:
"He's a Rebel"

In [14]:
'She asked, "How are you today?"'


Out[14]:
'She asked, "How are you today?"'

Just like the other two data objects we're familiar with (ints and floats), you can assign a string to a variable


In [15]:
greeting = "Hello, World!"

The print statement is often used for printing character strings:


In [16]:
print(greeting)


Hello, World!

But it can also print data types other than strings:


In [17]:
print("The area is " + area)


---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-17-c74a322ab524> in <module>()
----> 1 print("The area is " + area)

TypeError: Can't convert 'float' object to str implicitly

In [18]:
print("The area is " + str(area))


The area is 1256.6370614359173

In the above snipped, the number 600 (stored in the variable "area") is converted into a string before being printed out.

You can use the + operator to concatenate strings together:


In [19]:
statement = "Hello," + "World!"
print(statement)


Hello,World!

Don't forget the space between the strings, if you want one there.


In [20]:
statement = "Hello, " + "World!"
print(statement)


Hello, World!

You can use + to concatenate multiple strings in a single statement:


In [21]:
print("This " + "is " + "a " + "longer " + "statement.")


This is a longer statement.

If you have a lot of words to concatenate together, there are other, more efficient ways to do this. But this is fine for linking a few strings together.

Lists

Very often in a programming language, one wants to keep a group of similar items together. Python does this using a data type called lists.


In [22]:
days_of_the_week = ["Sunday","Monday","Tuesday","Wednesday","Thursday","Friday","Saturday"]

You can access members of the list using the index of that item:


In [23]:
days_of_the_week[2]


Out[23]:
'Tuesday'

Python lists, like C, but unlike Fortran, use 0 as the index of the first element of a list. Thus, in this example, the 0 element is "Sunday", 1 is "Monday", and so on. If you need to access the nth element from the end of the list, you can use a negative index. For example, the -1 element of a list is the last element:


In [24]:
days_of_the_week[-1]


Out[24]:
'Saturday'

You can add additional items to the list using the .append() command:


In [25]:
languages = ["Fortran","C","C++"]
languages.append("Python")
print(languages)


['Fortran', 'C', 'C++', 'Python']

The range() command is a convenient way to make sequential lists of numbers:


In [26]:
list(range(10))


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

Note that range(n) starts at 0 and gives the sequential list of integers less than n. If you want to start at a different number, use range(start,stop)


In [27]:
list(range(2,8))


Out[27]:
[2, 3, 4, 5, 6, 7]

The lists created above with range have a step of 1 between elements. You can also give a fixed step size via a third command:


In [28]:
evens = list(range(0,20,2))
evens


Out[28]:
[0, 2, 4, 6, 8, 10, 12, 14, 16, 18]

In [29]:
evens[3]


Out[29]:
6

Lists do not have to hold the same data type. For example,


In [30]:
["Today",7,99.3,""]


Out[30]:
['Today', 7, 99.3, '']

However, it's good (but not essential) to use lists for similar objects that are somehow logically connected. If you want to group different data types together into a composite data object, it's best to use tuples, which we will learn about below.

You can find out how long a list is using the len() command:


In [31]:
help(len)


Help on built-in function len in module builtins:

len(obj, /)
    Return the number of items in a container.


In [32]:
len(evens)


Out[32]:
10

Iteration, Indentation, and Blocks

One of the most useful things you can do with lists is to iterate through them, i.e. to go through each element one at a time. To do this in Python, we use the for statement:


In [33]:
for day in days_of_the_week:
    print(day)


Sunday
Monday
Tuesday
Wednesday
Thursday
Friday
Saturday

This code snippet goes through each element of the list called days_of_the_week and assigns it to the variable day. It then executes everything in the indented block (in this case only one line of code, the print statement) using those variable assignments. When the program has gone through every element of the list, it exists the block.

(Almost) every programming language defines blocks of code in some way. In Fortran, one uses END statements (ENDDO, ENDIF, etc.) to define code blocks. In C, C++, and Perl, one uses curly braces {} to define these blocks.

Python uses a colon (":"), followed by indentation level to define code blocks. Everything at a higher level of indentation is taken to be in the same block. In the above example the block was only a single line, but we could have had longer blocks as well:


In [34]:
for day in days_of_the_week:
    statement = "Today is " + day
    print(statement)


Today is Sunday
Today is Monday
Today is Tuesday
Today is Wednesday
Today is Thursday
Today is Friday
Today is Saturday

The range() command is particularly useful with the for statement to execute loops of a specified length:


In [35]:
for i in range(20):
    print("The square of ",i," is ",i*i)


The square of  0  is  0
The square of  1  is  1
The square of  2  is  4
The square of  3  is  9
The square of  4  is  16
The square of  5  is  25
The square of  6  is  36
The square of  7  is  49
The square of  8  is  64
The square of  9  is  81
The square of  10  is  100
The square of  11  is  121
The square of  12  is  144
The square of  13  is  169
The square of  14  is  196
The square of  15  is  225
The square of  16  is  256
The square of  17  is  289
The square of  18  is  324
The square of  19  is  361

Slicing

Lists and strings have something in common that you might not suspect: they can both be treated as sequences. You already know that you can iterate through the elements of a list. You can also iterate through the letters in a string:


In [36]:
for letter in "Sunday":
    print(letter)


S
u
n
d
a
y

This is only occasionally useful. Slightly more useful is the slicing operation, which you can also use on any sequence. We already know that we can use indexing to get the first element of a list:


In [37]:
days_of_the_week[0]


Out[37]:
'Sunday'

If we want the list containing the first two elements of a list, we can do this via


In [38]:
days_of_the_week[0:2]


Out[38]:
['Sunday', 'Monday']

or simply


In [39]:
days_of_the_week[:2]


Out[39]:
['Sunday', 'Monday']

If we want the last items of the list, we can do this with negative slicing:


In [40]:
days_of_the_week[-2:]


Out[40]:
['Friday', 'Saturday']

which is somewhat logically consistent with negative indices accessing the last elements of the list.

You can do:


In [41]:
workdays = days_of_the_week[1:6]
print(workdays)


['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday']

Since strings are sequences, you can also do this to them:


In [42]:
day = "Sunday"
abbreviation = day[:3]
print(abbreviation)


Sun

If we really want to get fancy, we can pass a third element into the slice, which specifies a step length (just like a third argument to the range() function specifies the step):


In [43]:
numbers = list(range(0,40))
evens = numbers[2::2]
evens


Out[43]:
[2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38]

Note that in this example I was even able to omit the second argument, so that the slice started at 2, went to the end of the list, and took every second element, to generate the list of even numbers less that 40.

Booleans and Truth Testing

We have now learned a few data types. We have integers and floating point numbers, strings, and lists to contain them. We have also learned about lists, a container that can hold any data type. We have learned to print things out, and to iterate over items in lists. We will now learn about boolean variables that can be either True or False.

We invariably need some concept of conditions in programming to control branching behavior, to allow a program to react differently to different situations. If it's Monday, I'll go to work, but if it's Sunday, I'll sleep in. To do this in Python, we use a combination of boolean variables, which evaluate to either True or False, and if statements, that control branching based on boolean values.

For example:


In [44]:
if day == "Sunday":
    print("Sleep in")
else:
    print("Go to work")


Sleep in

(Quick quiz: why did the snippet print "Go to work" here? What is the variable "day" set to?)

Let's take the snippet apart to see what happened. First, note the statement


In [45]:
day == "Sunday"


Out[45]:
True

If we evaluate it by itself, as we just did, we see that it returns a boolean value, False. The "==" operator performs equality testing. If the two items are equal, it returns True, otherwise it returns False. In this case, it is comparing two variables, the string "Sunday", and whatever is stored in the variable "day", which, in this case, is the other string "Saturday". Since the two strings are not equal to each other, the truth test has the false value.

The if statement that contains the truth test is followed by a code block (a colon followed by an indented block of code). If the boolean is true, it executes the code in that block. Since it is false in the above example, we don't see that code executed.

The first block of code is followed by an else statement, which is executed if nothing else in the above if statement is true. Since the value was false, this code is executed, which is why we see "Go to work".

You can compare any data types in Python:


In [46]:
1 == 2


Out[46]:
False

In [47]:
50 == 2*25


Out[47]:
True

In [48]:
3 < 3.14159


Out[48]:
True

In [49]:
1 == 1.0


Out[49]:
True

In [50]:
1 != 0


Out[50]:
True

In [51]:
1 <= 2


Out[51]:
True

In [52]:
1 >= 1


Out[52]:
True

We see a few other boolean operators here, all of which which should be self-explanatory. Less than, equality, non-equality, and so on.

Particularly interesting is the 1 == 1.0 test, which is true, since even though the two objects are different data types (integer and floating point number), they have the same value. There is another boolean operator is, that tests whether two objects are the same object:


In [53]:
1 is 1.0


Out[53]:
False

We can do boolean tests on lists as well:


In [54]:
[1,2,3] == [1,2,4]


Out[54]:
False

In [55]:
[1,2,3] < [1,2,4]


Out[55]:
True

Finally, note that you can also string multiple comparisons together, which can result in very intuitive tests:


In [56]:
hours = 5
0 < hours < 24


Out[56]:
True

If statements can have elif parts ("else if"), in addition to if/else parts. For example:


In [57]:
if day == "Sunday":
    print("Sleep in")
elif day == "Saturday":
    print("Do chores")
else:
    print("Go to work")


Sleep in

Of course we can combine if statements with for loops, to make a snippet that is almost interesting:


In [58]:
for day in days_of_the_week:
    statement = "Today is " + day
    print(statement)
    if day == "Sunday":
        print("   Sleep in")
    elif day == "Saturday":
        print("   Do chores")
    else:
        print("   Go to work")


Today is Sunday
   Sleep in
Today is Monday
   Go to work
Today is Tuesday
   Go to work
Today is Wednesday
   Go to work
Today is Thursday
   Go to work
Today is Friday
   Go to work
Today is Saturday
   Do chores

This is something of an advanced topic, but ordinary data types have boolean values associated with them, and, indeed, in early versions of Python there was not a separate boolean object. Essentially, anything that was a 0 value (the integer or floating point 0, an empty string "", or an empty list []) was False, and everything else was true. You can see the boolean value of any data object using the bool() function.


In [59]:
bool(1)


Out[59]:
True

In [60]:
bool(0)


Out[60]:
False

In [61]:
bool(["This "," is "," a "," list"])


Out[61]:
True

Code Example: The Fibonacci Sequence

The Fibonacci sequence is a sequence in math that starts with 0 and 1, and then each successive entry is the sum of the previous two. Thus, the sequence goes 0,1,1,2,3,5,8,13,21,34,55,89,...

A very common exercise in programming books is to compute the Fibonacci sequence up to some number n. First I'll show the code, then I'll discuss what it is doing.


In [62]:
n = 10
sequence = [0,1]
for i in range(2,n): # This is going to be a problem if we ever set n <= 2!
    sequence.append(sequence[i-1]+sequence[i-2])
print(sequence)


[0, 1, 1, 2, 3, 5, 8, 13, 21, 34]

Let's go through this line by line. First, we define the variable n, and set it to the integer 20. n is the length of the sequence we're going to form, and should probably have a better variable name. We then create a variable called sequence, and initialize it to the list with the integers 0 and 1 in it, the first two elements of the Fibonacci sequence. We have to create these elements "by hand", since the iterative part of the sequence requires two previous elements.

We then have a for loop over the list of integers from 2 (the next element of the list) to n (the length of the sequence). After the colon, we see a hash tag "#", and then a comment that if we had set n to some number less than 2 we would have a problem. Comments in Python start with #, and are good ways to make notes to yourself or to a user of your code explaining why you did what you did. Better than the comment here would be to test to make sure the value of n is valid, and to complain if it isn't; we'll try this later.

In the body of the loop, we append to the list an integer equal to the sum of the two previous elements of the list.

After exiting the loop (ending the indentation) we then print out the whole list. That's it!

Functions

We might want to use the Fibonacci snippet with different sequence lengths. We could cut an paste the code into another cell, changing the value of n, but it's easier and more useful to make a function out of the code. We do this with the def statement in Python:


In [63]:
def fibonacci(sequence_length):
    "Return the Fibonacci sequence of length *sequence_length*"
    sequence = [0,1]
    if sequence_length < 1:
        print("Fibonacci sequence only defined for length 1 or greater")
        return
    if 0 < sequence_length < 3:
        return sequence[:sequence_length]
    for i in range(2,sequence_length): 
        sequence.append(sequence[i-1]+sequence[i-2])
    return sequence

We can now call fibonacci() for different sequence_lengths:


In [64]:
fibonacci(2)


Out[64]:
[0, 1]

In [65]:
fibonacci(12)


Out[65]:
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89]

We've introduced a several new features here. First, note that the function itself is defined as a code block (a colon followed by an indented block). This is the standard way that Python delimits things. Next, note that the first line of the function is a single string. This is called a docstring, and is a special kind of comment that is often available to people using the function through the python command line:


In [66]:
help(fibonacci)


Help on function fibonacci in module __main__:

fibonacci(sequence_length)
    Return the Fibonacci sequence of length *sequence_length*

If you define a docstring for all of your functions, it makes it easier for other people to use them, since they can get help on the arguments and return values of the function.

Next, note that rather than putting a comment in about what input values lead to errors, we have some testing of these values, followed by a warning if the value is invalid, and some conditional code to handle special cases.

Two More Data Structures: Tuples and Dictionaries

Before we end the Python overview, I wanted to touch on two more data structures that are very useful (and thus very common) in Python programs.

A tuple is a sequence object like a list or a string. It's constructed by grouping a sequence of objects together with commas, either without brackets, or with parentheses:


In [67]:
t = (1,2,'hi',9.0)
t


Out[67]:
(1, 2, 'hi', 9.0)

Tuples are like lists, in that you can access the elements using indices:


In [68]:
t[1]


Out[68]:
2

However, tuples are immutable, you can't append to them or change the elements of them:


In [69]:
t.append(7)


---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-69-50c7062b1d5f> in <module>()
----> 1 t.append(7)

AttributeError: 'tuple' object has no attribute 'append'

In [70]:
t[1]=77


---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-70-03cc8ba9c07d> in <module>()
----> 1 t[1]=77

TypeError: 'tuple' object does not support item assignment

Tuples are useful anytime you want to group different pieces of data together in an object, but don't want to create a full-fledged class (see below) for them. For example, let's say you want the Cartesian coordinates of some objects in your program. Tuples are a good way to do this:


In [71]:
('Bob',0.0,21.0)


Out[71]:
('Bob', 0.0, 21.0)

Again, it's not a necessary distinction, but one way to distinguish tuples and lists is that tuples are a collection of different things, here a name, and x and y coordinates, whereas a list is a collection of similar things, like if we wanted a list of those coordinates:


In [72]:
positions = [
             ('Bob',0.0,21.0),
             ('Cat',2.5,13.1),
             ('Dog',33.0,1.2)
             ]

Tuples can be used when functions return more than one value. Say we wanted to compute the smallest x- and y-coordinates of the above list of objects. We could write:


In [73]:
def minmax(objects):
    minx = 1e20 # These are set to really big numbers
    miny = 1e20
    for obj in objects:
        name,x,y = obj
        if x < minx: 
            minx = x
        if y < miny:
            miny = y
    return minx,miny

x,y = minmax(positions)
print(x,y)


0.0 1.2

Dictionaries are an object called "mappings" or "associative arrays" in other languages. Whereas a list associates an integer index with a set of objects:


In [74]:
mylist = [1,2,9,21]

The index in a dictionary is called the key, and the corresponding dictionary entry is the value. A dictionary can use (almost) anything as the key. Whereas lists are formed with square brackets [], dictionaries use curly brackets {}:


In [75]:
ages = {"Rick": 46, "Bob": 86, "Fred": 21}
print("Rick's age is ",ages["Rick"])


Rick's age is  46

There's also a convenient way to create dictionaries without having to quote the keys.


In [76]:
dict(Rick=46,Bob=86,Fred=20)


Out[76]:
{'Bob': 86, 'Fred': 20, 'Rick': 46}

The len() command works on both tuples and dictionaries:


In [77]:
len(t)


Out[77]:
4

In [78]:
len(ages)


Out[78]:
3

Conclusion of the Python Overview

There is, of course, much more to the language than I've covered here. I've tried to keep this brief enough so that you can jump in and start using Python to simplify your life and work. My own experience in learning new things is that the information doesn't "stick" unless you try and use it for something in real life.

You will no doubt need to learn more as you go. I've listed several other good references, including the Python Tutorial and Learn Python the Hard Way. Additionally, now is a good time to start familiarizing yourself with the Python Documentation, and, in particular, the Python Language Reference.

Tim Peters, one of the earliest and most prolific Python contributors, wrote the "Zen of Python", which can be accessed via the "import this" command:


In [79]:
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!

No matter how experienced a programmer you are, these are words to meditate on.

II. Numpy and Scipy

Numpy contains core routines for doing fast vector, matrix, and linear algebra-type operations in Python. Scipy contains additional routines for optimization, special functions, and so on. Both contain modules written in C and Fortran so that they're as fast as possible. Together, they give Python roughly the same capability that the Matlab program offers. (In fact, if you're an experienced Matlab user, there a guide to Numpy for Matlab users just for you.)

Making vectors and matrices

Fundamental to both Numpy and Scipy is the ability to work with vectors and matrices. You can create vectors from lists using the array command:


In [80]:
import numpy as np
import scipy as sp

In [140]:
array = np.array([1,2,3,4,5,6])
array


Out[140]:
array([1, 2, 3, 4, 5, 6])

size of the array


In [141]:
array.shape


Out[141]:
(6,)

To build matrices, you can either use the array command with lists of lists:


In [135]:
mat = np.array([[0,1],[1,0]])
mat


Out[135]:
array([[0, 1],
       [1, 0]])

Add a column of ones to mat


In [142]:
mat2 = np.c_[mat, np.ones(2)]
mat2


Out[142]:
array([[ 0.,  1.,  1.],
       [ 1.,  0.,  1.]])

size of a matrix


In [143]:
mat2.shape


Out[143]:
(2, 3)

You can also form empty (zero) matrices of arbitrary shape (including vectors, which Numpy treats as vectors with one row), using the zeros command:


In [83]:
np.zeros((3,3))


Out[83]:
array([[ 0.,  0.,  0.],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.]])

There's also an identity command that behaves as you'd expect:


In [84]:
np.identity(4)


Out[84]:
array([[ 1.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.],
       [ 0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  1.]])

as well as a ones command.

Linspace, matrix functions, and plotting

The linspace command makes a linear array of points from a starting to an ending value.


In [85]:
np.linspace(0,1)


Out[85]:
array([ 0.        ,  0.02040816,  0.04081633,  0.06122449,  0.08163265,
        0.10204082,  0.12244898,  0.14285714,  0.16326531,  0.18367347,
        0.20408163,  0.2244898 ,  0.24489796,  0.26530612,  0.28571429,
        0.30612245,  0.32653061,  0.34693878,  0.36734694,  0.3877551 ,
        0.40816327,  0.42857143,  0.44897959,  0.46938776,  0.48979592,
        0.51020408,  0.53061224,  0.55102041,  0.57142857,  0.59183673,
        0.6122449 ,  0.63265306,  0.65306122,  0.67346939,  0.69387755,
        0.71428571,  0.73469388,  0.75510204,  0.7755102 ,  0.79591837,
        0.81632653,  0.83673469,  0.85714286,  0.87755102,  0.89795918,
        0.91836735,  0.93877551,  0.95918367,  0.97959184,  1.        ])

If you provide a third argument, it takes that as the number of points in the space. If you don't provide the argument, it gives a length 50 linear space.


In [86]:
np.linspace(0,1,11)


Out[86]:
array([ 0. ,  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9,  1. ])

linspace is an easy way to make coordinates for plotting. Functions in the numpy library (all of which are imported into IPython notebook) can act on an entire vector (or even a matrix) of points at once. Thus,


In [87]:
x = np.linspace(0,2*np.pi)
np.sin(x)


Out[87]:
array([  0.00000000e+00,   1.27877162e-01,   2.53654584e-01,
         3.75267005e-01,   4.90717552e-01,   5.98110530e-01,
         6.95682551e-01,   7.81831482e-01,   8.55142763e-01,
         9.14412623e-01,   9.58667853e-01,   9.87181783e-01,
         9.99486216e-01,   9.95379113e-01,   9.74927912e-01,
         9.38468422e-01,   8.86599306e-01,   8.20172255e-01,
         7.40277997e-01,   6.48228395e-01,   5.45534901e-01,
         4.33883739e-01,   3.15108218e-01,   1.91158629e-01,
         6.40702200e-02,  -6.40702200e-02,  -1.91158629e-01,
        -3.15108218e-01,  -4.33883739e-01,  -5.45534901e-01,
        -6.48228395e-01,  -7.40277997e-01,  -8.20172255e-01,
        -8.86599306e-01,  -9.38468422e-01,  -9.74927912e-01,
        -9.95379113e-01,  -9.99486216e-01,  -9.87181783e-01,
        -9.58667853e-01,  -9.14412623e-01,  -8.55142763e-01,
        -7.81831482e-01,  -6.95682551e-01,  -5.98110530e-01,
        -4.90717552e-01,  -3.75267005e-01,  -2.53654584e-01,
        -1.27877162e-01,  -2.44929360e-16])

In conjunction with matplotlib, this is a nice way to plot things:


In [88]:
%matplotlib inline
import matplotlib.pyplot as plt

In [89]:
plt.plot(x,np.sin(x))


Out[89]:
[<matplotlib.lines.Line2D at 0x7f2c958344e0>]

Matrix operations

Matrix objects act sensibly when multiplied by scalars:


In [90]:
0.125*np.identity(3)


Out[90]:
array([[ 0.125,  0.   ,  0.   ],
       [ 0.   ,  0.125,  0.   ],
       [ 0.   ,  0.   ,  0.125]])

as well as when you add two matrices together. (However, the matrices have to be the same shape.)


In [91]:
np.identity(2) + np.array([[1,1],[1,2]])


Out[91]:
array([[ 2.,  1.],
       [ 1.,  3.]])

Something that confuses Matlab users is that the times (*) operator give element-wise multiplication rather than matrix multiplication:


In [92]:
np.identity(2)*np.ones((2,2))


Out[92]:
array([[ 1.,  0.],
       [ 0.,  1.]])

To get matrix multiplication, you need the dot command:


In [93]:
np.dot(np.identity(2),np.ones((2,2)))


Out[93]:
array([[ 1.,  1.],
       [ 1.,  1.]])

dot can also do dot products (duh!):


In [94]:
v = np.array([3,4])
np.sqrt(np.dot(v,v))


Out[94]:
5.0

as well as matrix-vector products.

There are determinant, inverse, and transpose functions that act as you would suppose. Transpose can be abbreviated with ".T" at the end of a matrix object:


In [95]:
m = np.array([[1,2],[3,4]])
m.T


Out[95]:
array([[1, 3],
       [2, 4]])

In [160]:
np.linalg.inv(m)


Out[160]:
array([[-2. ,  1. ],
       [ 1.5, -0.5]])

There's also a diag() function that takes a list or a vector and puts it along the diagonal of a square matrix.


In [96]:
np.diag([1,2,3,4,5])


Out[96]:
array([[1, 0, 0, 0, 0],
       [0, 2, 0, 0, 0],
       [0, 0, 3, 0, 0],
       [0, 0, 0, 4, 0],
       [0, 0, 0, 0, 5]])

We'll find this useful later on.

Least squares fitting

Very often we deal with some data that we want to fit to some sort of expected behavior. Say we have the following:


In [97]:
raw_data = """\
3.1905781584582433,0.028208609537968457
4.346895074946466,0.007160804747670053
5.374732334047101,0.0046962988461934805
8.201284796573875,0.0004614473299618756
10.899357601713055,0.00005038370219939726
16.295503211991434,4.377451812785309e-7
21.82012847965739,3.0799922117601088e-9
32.48394004282656,1.524776208284536e-13
43.53319057815846,5.5012073588707224e-18"""

There's a section below on parsing CSV data. We'll steal the parser from that. For an explanation, skip ahead to that section. Otherwise, just assume that this is a way to parse that text into a numpy array that we can plot and do other analyses with.


In [98]:
data = []
for line in raw_data.splitlines():
    words = line.split(',')
    data.append(words)
data = np.array(data, dtype=np.float)

In [99]:
data


Out[99]:
array([[  3.19057816e+00,   2.82086095e-02],
       [  4.34689507e+00,   7.16080475e-03],
       [  5.37473233e+00,   4.69629885e-03],
       [  8.20128480e+00,   4.61447330e-04],
       [  1.08993576e+01,   5.03837022e-05],
       [  1.62955032e+01,   4.37745181e-07],
       [  2.18201285e+01,   3.07999221e-09],
       [  3.24839400e+01,   1.52477621e-13],
       [  4.35331906e+01,   5.50120736e-18]])

In [100]:
data[:, 0]


Out[100]:
array([  3.19057816,   4.34689507,   5.37473233,   8.2012848 ,
        10.8993576 ,  16.29550321,  21.82012848,  32.48394004,  43.53319058])

In [101]:
plt.title("Raw Data")
plt.xlabel("Distance")
plt.plot(data[:,0],data[:,1],'bo')


Out[101]:
[<matplotlib.lines.Line2D at 0x7f2c95865be0>]

Since we expect the data to have an exponential decay, we can plot it using a semi-log plot.


In [102]:
plt.title("Raw Data")
plt.xlabel("Distance")
plt.semilogy(data[:,0],data[:,1],'bo')


Out[102]:
[<matplotlib.lines.Line2D at 0x7f2c8ff81470>]

For a pure exponential decay like this, we can fit the log of the data to a straight line. The above plot suggests this is a good approximation. Given a function $$ y = Ae^{-ax} $$ $$ \log(y) = \log(A) - ax$$ Thus, if we fit the log of the data versus x, we should get a straight line with slope $a$, and an intercept that gives the constant $A$.

There's a numpy function called polyfit that will fit data to a polynomial form. We'll use this to fit to a straight line (a polynomial of order 1)


In [103]:
params = sp.polyfit(data[:,0],np.log(data[:,1]),1)
a = params[0]
A = np.exp(params[1])

Let's see whether this curve fits the data.


In [104]:
x = np.linspace(1,45)
plt.title("Raw Data")
plt.xlabel("Distance")
plt.semilogy(data[:,0],data[:,1],'bo')
plt.semilogy(x,A*np.exp(a*x),'b-')


Out[104]:
[<matplotlib.lines.Line2D at 0x7f2c8fca0438>]

If we have more complicated functions, we may not be able to get away with fitting to a simple polynomial. Consider the following data:


In [150]:
gauss_data = """\
-0.9902286902286903,1.4065274110372852e-19
-0.7566104566104566,2.2504438576596563e-18
-0.5117810117810118,1.9459459459459454
-0.31887271887271884,10.621621621621626
-0.250997150997151,15.891891891891893
-0.1463309463309464,23.756756756756754
-0.07267267267267263,28.135135135135133
-0.04426734426734419,29.02702702702703
-0.0015939015939017698,29.675675675675677
0.04689304689304685,29.10810810810811
0.0840994840994842,27.324324324324326
0.1700546700546699,22.216216216216214
0.370878570878571,7.540540540540545
0.5338338338338338,1.621621621621618
0.722014322014322,0.08108108108108068
0.9926849926849926,-0.08108108108108646"""

data = []
for line in gauss_data.splitlines():
    words = line.split(',')
    data.append(words)
data = np.array(data, dtype=np.float)

plt.plot(data[:,0],data[:,1],'bo')


Out[150]:
[<matplotlib.lines.Line2D at 0x7f2c8fbd7240>]

This data looks more Gaussian than exponential. If we wanted to, we could use polyfit for this as well, but let's use the curve_fit function from Scipy, which can fit to arbitrary functions. You can learn more using help(curve_fit).

First define a general Gaussian function to fit to.


In [151]:
def gauss(x,A,a): 
    return A*np.exp(a*x**2)

Now fit to it using curve_fit:


In [152]:
from scipy.optimize import curve_fit

params,conv = curve_fit(gauss,data[:,0],data[:,1])
x = np.linspace(-1,1)
plt.plot(data[:,0],data[:,1],'bo')
A,a = params
plt.plot(x,gauss(x,A,a),'b-')


Out[152]:
[<matplotlib.lines.Line2D at 0x7f2c8fac29e8>]

The curve_fit routine we just used is built on top of a very good general minimization capability in Scipy. You can learn more at the scipy documentation pages.

Monte Carlo and random numbers

Many methods in scientific computing rely on Monte Carlo integration, where a sequence of (pseudo) random numbers are used to approximate the integral of a function. Python has good random number generators in the standard library. The random() function gives pseudorandom numbers uniformly distributed between 0 and 1:


In [153]:
from random import random
rands = []
for i in range(100):
    rands.append(random())
plt.plot(rands)


Out[153]:
[<matplotlib.lines.Line2D at 0x7f2c8fa3ee80>]

random() uses the Mersenne Twister algorithm, which is a highly regarded pseudorandom number generator. There are also functions to generate random integers, to randomly shuffle a list, and functions to pick random numbers from a particular distribution, like the normal distribution:


In [154]:
from random import gauss
grands = []
for i in range(100):
    grands.append(gauss(0,1))
plt.plot(grands)


Out[154]:
[<matplotlib.lines.Line2D at 0x7f2c8f9a5b70>]

It is generally more efficient to generate a list of random numbers all at once, particularly if you're drawing from a non-uniform distribution. Numpy has functions to generate vectors and matrices of particular types of random distributions.


In [155]:
plt.plot(np.random.rand(100))


Out[155]:
[<matplotlib.lines.Line2D at 0x7f2c8f985208>]

Slicing numpy arrays and matrices


In [156]:
data.shape


Out[156]:
(16, 2)

Select second column


In [157]:
data[:, 1]


Out[157]:
array([  1.40652741e-19,   2.25044386e-18,   1.94594595e+00,
         1.06216216e+01,   1.58918919e+01,   2.37567568e+01,
         2.81351351e+01,   2.90270270e+01,   2.96756757e+01,
         2.91081081e+01,   2.73243243e+01,   2.22162162e+01,
         7.54054054e+00,   1.62162162e+00,   8.10810811e-02,
        -8.10810811e-02])

Select the first 5 rows


In [158]:
data[:5, :]


Out[158]:
array([[ -9.90228690e-01,   1.40652741e-19],
       [ -7.56610457e-01,   2.25044386e-18],
       [ -5.11781012e-01,   1.94594595e+00],
       [ -3.18872719e-01,   1.06216216e+01],
       [ -2.50997151e-01,   1.58918919e+01]])

Select the second row and the last column


In [159]:
data[1, -1]


Out[159]:
2.2504438576596563e-18

III. Intermediate Python

Output Parsing

As more and more of our day-to-day work is being done on and through computers, we increasingly have output that one program writes, often in a text file, that we need to analyze in one way or another, and potentially feed that output into another file.

Suppose we have the following output:


In [111]:
myoutput = """\
@ Step       Energy      Delta E   Gmax     Grms     Xrms     Xmax   Walltime
@ ---- ---------------- -------- -------- -------- -------- -------- --------
@    0   -6095.12544083  0.0D+00  0.03686  0.00936  0.00000  0.00000   1391.5
@    1   -6095.25762870 -1.3D-01  0.00732  0.00168  0.32456  0.84140  10468.0
@    2   -6095.26325979 -5.6D-03  0.00233  0.00056  0.06294  0.14009  11963.5
@    3   -6095.26428124 -1.0D-03  0.00109  0.00024  0.03245  0.10269  13331.9
@    4   -6095.26463203 -3.5D-04  0.00057  0.00013  0.02737  0.09112  14710.8
@    5   -6095.26477615 -1.4D-04  0.00043  0.00009  0.02259  0.08615  20211.1
@    6   -6095.26482624 -5.0D-05  0.00015  0.00002  0.00831  0.03147  21726.1
@    7   -6095.26483584 -9.6D-06  0.00021  0.00004  0.01473  0.05265  24890.5
@    8   -6095.26484405 -8.2D-06  0.00005  0.00001  0.00555  0.01929  26448.7
@    9   -6095.26484599 -1.9D-06  0.00003  0.00001  0.00164  0.00564  27258.1
@   10   -6095.26484676 -7.7D-07  0.00003  0.00001  0.00161  0.00553  28155.3
@   11   -6095.26484693 -1.8D-07  0.00002  0.00000  0.00054  0.00151  28981.7
@   11   -6095.26484693 -1.8D-07  0.00002  0.00000  0.00054  0.00151  28981.7"""

This output actually came from a geometry optimization of a Silicon cluster using the NWChem quantum chemistry suite. At every step the program computes the energy of the molecular geometry, and then changes the geometry to minimize the computed forces, until the energy converges. I obtained this output via the unix command

% grep @ nwchem.out

since NWChem is nice enough to precede the lines that you need to monitor job progress with the '@' symbol.

We could do the entire analysis in Python; I'll show how to do this later on, but first let's focus on turning this code into a usable Python object that we can plot.

First, note that the data is entered into a multi-line string. When Python sees three quote marks """ or ''' it treats everything following as part of a single string, including newlines, tabs, and anything else, until it sees the same three quote marks (""" has to be followed by another """, and ''' has to be followed by another ''') again. This is a convenient way to quickly dump data into Python, and it also reinforces the important idea that you don't have to open a file and deal with it one line at a time. You can read everything in, and deal with it as one big chunk.

The first thing we'll do, though, is to split the big string into a list of strings, since each line corresponds to a separate piece of data. We will use the splitlines() function on the big myout string to break it into a new element every time it sees a newline (\n) character:


In [112]:
lines = myoutput.splitlines()
lines


Out[112]:
['@ Step       Energy      Delta E   Gmax     Grms     Xrms     Xmax   Walltime',
 '@ ---- ---------------- -------- -------- -------- -------- -------- --------',
 '@    0   -6095.12544083  0.0D+00  0.03686  0.00936  0.00000  0.00000   1391.5',
 '@    1   -6095.25762870 -1.3D-01  0.00732  0.00168  0.32456  0.84140  10468.0',
 '@    2   -6095.26325979 -5.6D-03  0.00233  0.00056  0.06294  0.14009  11963.5',
 '@    3   -6095.26428124 -1.0D-03  0.00109  0.00024  0.03245  0.10269  13331.9',
 '@    4   -6095.26463203 -3.5D-04  0.00057  0.00013  0.02737  0.09112  14710.8',
 '@    5   -6095.26477615 -1.4D-04  0.00043  0.00009  0.02259  0.08615  20211.1',
 '@    6   -6095.26482624 -5.0D-05  0.00015  0.00002  0.00831  0.03147  21726.1',
 '@    7   -6095.26483584 -9.6D-06  0.00021  0.00004  0.01473  0.05265  24890.5',
 '@    8   -6095.26484405 -8.2D-06  0.00005  0.00001  0.00555  0.01929  26448.7',
 '@    9   -6095.26484599 -1.9D-06  0.00003  0.00001  0.00164  0.00564  27258.1',
 '@   10   -6095.26484676 -7.7D-07  0.00003  0.00001  0.00161  0.00553  28155.3',
 '@   11   -6095.26484693 -1.8D-07  0.00002  0.00000  0.00054  0.00151  28981.7',
 '@   11   -6095.26484693 -1.8D-07  0.00002  0.00000  0.00054  0.00151  28981.7']

Splitting is a big concept in text processing. We used splitlines() here, and we will use the more general split() function below to split each line into whitespace-delimited words.

We now want to do three things:

  • Skip over the lines that don't carry any information
  • Break apart each line that does carry information and grab the pieces we want
  • Turn the resulting data into something that we can plot.

For this data, we really only want the Energy column, the Gmax column (which contains the maximum gradient at each step), and perhaps the Walltime column.

Since the data is now in a list of lines, we can iterate over it:


In [113]:
for line in lines[2:]:
    # do something with each line
    words = line.split()

Let's examine what we just did: first, we used a for loop to iterate over each line. However, we skipped the first two (the lines[2:] only takes the lines starting from index 2), since lines[0] contained the title information, and lines[1] contained underscores.

We then split each line into chunks (which we're calling "words", even though in most cases they're numbers) using the string split() command. Here's what split does:


In [114]:
lines[2].split()


Out[114]:
['@',
 '0',
 '-6095.12544083',
 '0.0D+00',
 '0.03686',
 '0.00936',
 '0.00000',
 '0.00000',
 '1391.5']

This is almost exactly what we want. We just have to now pick the fields we want:


In [115]:
for line in lines[2:]:
    # do something with each line
    words = line.split()
    energy = words[2]
    gmax = words[4]
    time = words[8]
    print(energy,gmax,time)


-6095.12544083 0.03686 1391.5
-6095.25762870 0.00732 10468.0
-6095.26325979 0.00233 11963.5
-6095.26428124 0.00109 13331.9
-6095.26463203 0.00057 14710.8
-6095.26477615 0.00043 20211.1
-6095.26482624 0.00015 21726.1
-6095.26483584 0.00021 24890.5
-6095.26484405 0.00005 26448.7
-6095.26484599 0.00003 27258.1
-6095.26484676 0.00003 28155.3
-6095.26484693 0.00002 28981.7
-6095.26484693 0.00002 28981.7

This is fine for printing things out, but if we want to do something with the data, either make a calculation with it or pass it into a plotting, we need to convert the strings into regular floating point numbers. We can use the float() command for this. We also need to save it in some form. I'll do this as follows:


In [116]:
data = []
for line in lines[2:]:
    # do something with each line
    words = line.split()
    energy = float(words[2])
    gmax = float(words[4])
    time = float(words[8])
    data.append((energy,gmax,time))
data = np.array(data)

We now have our data in a numpy array, so we can choose columns to print:


In [117]:
plt.plot(data[:,0])
plt.xlabel('step')
plt.ylabel('Energy (hartrees)')
plt.title('Convergence of NWChem geometry optimization for Si cluster')


Out[117]:
<matplotlib.text.Text at 0x7f2c8fc8b550>

In [118]:
energies = data[:,0]
minE = min(energies)
energies_eV = 27.211*(energies-minE)
plt.plot(energies_eV)
plt.xlabel('step')
plt.ylabel('Energy (eV)')
plt.title('Convergence of NWChem geometry optimization for Si cluster')


Out[118]:
<matplotlib.text.Text at 0x7f2c8fdfcc50>

This gives us the output in a form that we can think about: 4 eV is a fairly substantial energy change (chemical bonds are roughly this magnitude of energy), and most of the energy decrease was obtained in the first geometry iteration.

We mentioned earlier that we don't have to rely on grep to pull out the relevant lines for us. The string module has a lot of useful functions we can use for this. Among them is the startswith function. For example:


In [119]:
lines = """\
                 ----------------------------------------
                 |  WALL  |       0.45   |     443.61   |
                 ----------------------------------------

@ Step       Energy      Delta E   Gmax     Grms     Xrms     Xmax   Walltime
@ ---- ---------------- -------- -------- -------- -------- -------- --------
@    0   -6095.12544083  0.0D+00  0.03686  0.00936  0.00000  0.00000   1391.5
                                                       ok       ok



                                Z-matrix (autoz)
                                --------
""".splitlines()

for line in lines:
    if line.startswith('@'):
        print(line)


@ Step       Energy      Delta E   Gmax     Grms     Xrms     Xmax   Walltime
@ ---- ---------------- -------- -------- -------- -------- -------- --------
@    0   -6095.12544083  0.0D+00  0.03686  0.00936  0.00000  0.00000   1391.5

and we've successfully grabbed all of the lines that begin with the @ symbol.

The real value in a language like Python is that it makes it easy to take additional steps to analyze data in this fashion, which means you are thinking more about your data, and are more likely to see important patterns.

Optional arguments

You will recall that the linspace function can take either two arguments (for the starting and ending points):


In [120]:
np.linspace(0,1)


Out[120]:
array([ 0.        ,  0.02040816,  0.04081633,  0.06122449,  0.08163265,
        0.10204082,  0.12244898,  0.14285714,  0.16326531,  0.18367347,
        0.20408163,  0.2244898 ,  0.24489796,  0.26530612,  0.28571429,
        0.30612245,  0.32653061,  0.34693878,  0.36734694,  0.3877551 ,
        0.40816327,  0.42857143,  0.44897959,  0.46938776,  0.48979592,
        0.51020408,  0.53061224,  0.55102041,  0.57142857,  0.59183673,
        0.6122449 ,  0.63265306,  0.65306122,  0.67346939,  0.69387755,
        0.71428571,  0.73469388,  0.75510204,  0.7755102 ,  0.79591837,
        0.81632653,  0.83673469,  0.85714286,  0.87755102,  0.89795918,
        0.91836735,  0.93877551,  0.95918367,  0.97959184,  1.        ])

or it can take three arguments, for the starting point, the ending point, and the number of points:


In [121]:
np.linspace(0,1,5)


Out[121]:
array([ 0.  ,  0.25,  0.5 ,  0.75,  1.  ])

You can also pass in keywords to exclude the endpoint:


In [122]:
np.linspace(0,1,5,endpoint=False)


Out[122]:
array([ 0. ,  0.2,  0.4,  0.6,  0.8])

Right now, we only know how to specify functions that have a fixed number of arguments. We'll learn how to do the more general cases here.

If we're defining a simple version of linspace, we would start with:


In [123]:
def my_linspace(start,end):
    npoints = 50
    v = []
    d = (end-start)/float(npoints-1)
    for i in range(npoints):
        v.append(start + i*d)
    return v
my_linspace(0,1)


Out[123]:
[0.0,
 0.02040816326530612,
 0.04081632653061224,
 0.061224489795918366,
 0.08163265306122448,
 0.1020408163265306,
 0.12244897959183673,
 0.14285714285714285,
 0.16326530612244897,
 0.18367346938775508,
 0.2040816326530612,
 0.22448979591836732,
 0.24489795918367346,
 0.26530612244897955,
 0.2857142857142857,
 0.3061224489795918,
 0.32653061224489793,
 0.3469387755102041,
 0.36734693877551017,
 0.3877551020408163,
 0.4081632653061224,
 0.42857142857142855,
 0.44897959183673464,
 0.4693877551020408,
 0.4897959183673469,
 0.5102040816326531,
 0.5306122448979591,
 0.5510204081632653,
 0.5714285714285714,
 0.5918367346938775,
 0.6122448979591836,
 0.6326530612244897,
 0.6530612244897959,
 0.673469387755102,
 0.6938775510204082,
 0.7142857142857142,
 0.7346938775510203,
 0.7551020408163265,
 0.7755102040816326,
 0.7959183673469387,
 0.8163265306122448,
 0.836734693877551,
 0.8571428571428571,
 0.8775510204081632,
 0.8979591836734693,
 0.9183673469387754,
 0.9387755102040816,
 0.9591836734693877,
 0.9795918367346939,
 0.9999999999999999]

We can add an optional argument by specifying a default value in the argument list:


In [124]:
def my_linspace(start,end,npoints = 50):
    v = []
    d = (end-start)/float(npoints-1)
    for i in range(npoints):
        v.append(start + i*d)
    return v

This gives exactly the same result if we don't specify anything:


In [125]:
my_linspace(0,1)


Out[125]:
[0.0,
 0.02040816326530612,
 0.04081632653061224,
 0.061224489795918366,
 0.08163265306122448,
 0.1020408163265306,
 0.12244897959183673,
 0.14285714285714285,
 0.16326530612244897,
 0.18367346938775508,
 0.2040816326530612,
 0.22448979591836732,
 0.24489795918367346,
 0.26530612244897955,
 0.2857142857142857,
 0.3061224489795918,
 0.32653061224489793,
 0.3469387755102041,
 0.36734693877551017,
 0.3877551020408163,
 0.4081632653061224,
 0.42857142857142855,
 0.44897959183673464,
 0.4693877551020408,
 0.4897959183673469,
 0.5102040816326531,
 0.5306122448979591,
 0.5510204081632653,
 0.5714285714285714,
 0.5918367346938775,
 0.6122448979591836,
 0.6326530612244897,
 0.6530612244897959,
 0.673469387755102,
 0.6938775510204082,
 0.7142857142857142,
 0.7346938775510203,
 0.7551020408163265,
 0.7755102040816326,
 0.7959183673469387,
 0.8163265306122448,
 0.836734693877551,
 0.8571428571428571,
 0.8775510204081632,
 0.8979591836734693,
 0.9183673469387754,
 0.9387755102040816,
 0.9591836734693877,
 0.9795918367346939,
 0.9999999999999999]

But also let's us override the default value with a third argument:


In [126]:
my_linspace(0,1,5)


Out[126]:
[0.0, 0.25, 0.5, 0.75, 1.0]

We can add arbitrary keyword arguments to the function definition by putting a keyword argument **kwargs handle in:


In [127]:
def my_linspace(start,end,npoints=50,**kwargs):
    endpoint = kwargs.get('endpoint',True)
    v = []
    if endpoint:
        d = (end-start)/float(npoints-1)
    else:
        d = (end-start)/float(npoints)
    for i in range(npoints):
        v.append(start + i*d)
    return v
my_linspace(0,1,5,endpoint=False)


Out[127]:
[0.0, 0.2, 0.4, 0.6000000000000001, 0.8]

What the keyword argument construction does is to take any additional keyword arguments (i.e. arguments specified by name, like "endpoint=False"), and stick them into a dictionary called "kwargs" (you can call it anything you like, but it has to be preceded by two stars). You can then grab items out of the dictionary using the get command, which also lets you specify a default value. I realize it takes a little getting used to, but it is a common construction in Python code, and you should be able to recognize it.

There's an analogous *args that dumps any additional arguments into a list called "args". Think about the range function: it can take one (the endpoint), two (starting and ending points), or three (starting, ending, and step) arguments. How would we define this?


In [128]:
def my_range(*args):
    start = 0
    step = 1
    if len(args) == 1:
        end = args[0]
    elif len(args) == 2:
        start,end = args
    elif len(args) == 3:
        start,end,step = args
    else:
        raise Exception("Unable to parse arguments")
    v = []
    value = start
    while True:
        v.append(value)
        value += step
        if value > end: break
    return v

Note that we have defined a few new things you haven't seen before: a break statement, that allows us to exit a for loop if some conditions are met, and an exception statement, that causes the interpreter to exit with an error message. For example:


In [129]:
my_range()


---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
<ipython-input-129-0e8004dab150> in <module>()
----> 1 my_range()

<ipython-input-128-c34e09da2551> in my_range(*args)
      9         start,end,step = args
     10     else:
---> 11         raise Exception("Unable to parse arguments")
     12     v = []
     13     value = start

Exception: Unable to parse arguments

List Comprehensions and Generators

List comprehensions are a streamlined way to make lists. They look something like a list definition, with some logic thrown in. For example:


In [130]:
evens1 = [2*i for i in range(10)]
print(evens1)


[0, 2, 4, 6, 8, 10, 12, 14, 16, 18]

You can also put some boolean testing into the construct:


In [131]:
odds = [i for i in range(20) if i%2==1]
odds


Out[131]:
[1, 3, 5, 7, 9, 11, 13, 15, 17, 19]

Here i%2 is the remainder when i is divided by 2, so that i%2==1 is true if the number is odd. Even though this is a relative new addition to the language, it is now fairly common since it's so convenient.

iterators are a way of making virtual sequence objects. Consider if we had the nested loop structure:

for i in range(1000000):
    for j in range(1000000):

Inside the main loop, we make a list of 1,000,000 integers, just to loop over them one at a time. We don't need any of the additional things that a lists gives us, like slicing or random access, we just need to go through the numbers one at a time. And we're making 1,000,000 of them.

iterators are a way around this. For example, the xrange function is the iterator version of range. This simply makes a counter that is looped through in sequence, so that the analogous loop structure would look like:

for i in xrange(1000000):
    for j in xrange(1000000):

Even though we've only added two characters, we've dramatically sped up the code, because we're not making 1,000,000 big lists.

We can define our own iterators using the yield statement:


In [132]:
def evens_below(n):
    for i in range(n):
        if i%2 == 0:
            yield i
    return

for i in evens_below(9):
    print(i)


0
2
4
6
8

We can always turn an iterator into a list using the list command:


In [133]:
list(evens_below(9))


Out[133]:
[0, 2, 4, 6, 8]

There's a special syntax called a generator expression that looks a lot like a list comprehension:


In [134]:
evens_gen = (i for i in range(9) if i%2==0)
for i in evens_gen:
    print(i)


0
2
4
6
8