<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 A Crash Course in Python for Scientists
Rick Muller, Sandia National Laboratories
version 0.62, Updated Dec 15, 2016 by Ryan Smith, Cal State East Bay
version 0.63, Updated Oct 2017 by Ryan Smith, Cal State East Bay
Using Python 3.5.2 | Anaconda 4.1.1 
This work is licensed under a Creative Commons Attribution-ShareAlike 3.0 Unported License.
If skipping to other sections, it's a good idea to run this first
Step1:  Table of Contents
0. Preliminary
0.1 Why Python?
0.2 What You Need to Install
1. Python Overview
1.1 Using Python as a Calculator
1.2 Strings
1.3 Lists
1.4 Iteration, Indentation, and Blocks
1.5 Slicing
1.6 Booleans and Truth Testing
1.7 Code Example
Step2:  If you want to raise to an exponent, "^" is reserved for something else, so we use ** instead
Step3:  There used to be gotchas in division in python 2, like C or Fortran integer division, where division truncates the remainder and returns an integer. In version 3, Python returns a floating point number. If for some reason you are using Python 2, you can fix this by importing the module from the future features
Step4:  I prefer doing things with numpy though.. we already imported numpy as np in the first command box above, so we can just call it... I prefer things this way, so that code is always easily readable
Step5:  You can define variables using the equals (=) sign
Step6:  If you try to access a variable that you haven't yet defined, you get an error
Step7:  and you need to define it
Step8:  A few examples of common functions. Look up numpy library -- it has most mathematical functions you could want.
Step9:  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
Step10:  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.
1.2 Strings
Strings are lists of printable characters, and can be defined using either single quotes
Step11:  or double quotes
Step12:  But not both at the same time, unless you want one of the symbols to be part of the string.
Step13:  Just like the other two data objects we're familiar with (ints and floats), you can assign a string to a variable
Step14:  The print statement is often used for printing character strings
Step15:  But it can also print data types other than strings
Step16:  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
Step17:  Don't forget the space between the strings, if you want one there.
Step18:  You can use + to concatenate multiple strings in a single statement
Step19:  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.
1.3 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.
Step20:  You can access members of the list using the index of that item
Step21:  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
Step22:  You can add additional items to the list using the .append() command
Step23:  The range() command is a convenient way to make sequential lists of numbers
Step24:  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)
Step25:  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
Step26:  Lists do not have to hold the same data type. For example,
Step27:  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
Step28:  1.4 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
Step29:  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 ("
Step30:  The range() command is particularly useful with the for statement to execute loops of a specified length
Step31:  1.5 Slicing
Lists and strings have something in common that you might not suspect
Step32:  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
Step33:  If we want the list containing the first two elements of a list, we can do this via
Step34:  or simply
Step35:  If we want the last items of the list, we can do this with negative slicing
Step36:  which is somewhat logically consistent with negative indices accessing the last elements of the list.
You can do
Step37:  Since strings are sequences, you can also do this to them
Step38:  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)
Step39:  Note that in this example we omitted the first few arguments, so that the slice started at 6, went to the end of the list, and took every second element, to generate the list of even numbers lup to 20 (the last element in the original list).
1.6 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
Step40:  (Quick quiz
Step41:  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".
Try setting the day equal to "Sunday" and then running the above if/else statement. Did it work as you thought it would?
You can compare any data types in Python
Step42:  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
Step43:  Why is 1 not the same as 1.0? Different data type. You can check the data type
Step44:  We can do boolean tests on lists as well
Step45:  Finally, note that you can also string multiple comparisons together, which can result in very intuitive tests
Step46:  If statements can have elif parts ("else if"), in addition to if/else parts. For example
Step47:  Of course we can combine if statements with for loops, to make a snippet that is almost interesting
Step48:  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.
Step49:  1.7 Code Example
Step50:  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!
1.8 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
Step51:  We can now call fibonacci() for different sequence_lengths
Step52:  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
Step53:  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.
1.9 Recursion and Factorials
Functions can also call themselves, something that is often called recursion. We're going to experiment with recursion by computing the factorial function. The factorial is defined for a positive integer n as
$$ n! = n(n-1)(n-2)\cdots 1 $$
First, note that we don't need to write a function at all, since this is a function built into the standard math library. Let's use the help function to find out about it
Step54:  This is clearly what we want.
Step55:  However, if we did want to write a function ourselves, we could do recursively by noting that
$$ n! = n(n-1)!$$
The program then looks something like
Step56:  Recursion can be very elegant, and can lead to very simple programs.
1.10 Two More Data Structures
Step57:  Tuples are like lists, in that you can access the elements using indices
Step58:  However, tuples are immutable, you can't append to them or change the elements of them
Step59:  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
Step60:  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
Step61:  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
Step62:  Here we did two things with tuples you haven't seen before. First, we unpacked an object into a set of named variables using tuple assignment
Step63:  Dictionaries are an object called "mappings" or "associative arrays" in other languages. Whereas a list associates an integer index with a set of objects
Step64:  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 {}
Step65:  There's also a convenient way to create dictionaries without having to quote the keys.
Step66:  Notice in either case you are not choosing the ordering -- it is automagically grouped alphabetically.
The len() command works on both tuples and dictionaries
Step67:  1.11 Plotting with Matplotlib
We can generally understand trends in data by using a plotting program to chart it. Python has a wonderful plotting library called Matplotlib. The Jupyter notebook interface we are using for these notes has that functionality built in.
First off, it is important to import the library. We did this at the very beginning of this whole jupyter notebook, but here it is in case you've jumped straight here without running the first code line
Step68:  The %matplotlib inline command makes it so plots are within this notebook. To plot to a separate window, use instead
Step69:  Next lets generate the factorials.
Step70:  Now we use the Matplotlib function plot to compare the two.
Step71:  The factorial function grows much faster. In fact, you can't even see the Fibonacci sequence. It's not entirely surprising
Step72:  There are many more things you can do with Matplotlib. We'll be looking at some of them in the sections to come. In the meantime, if you want an idea of the different things you can do, look at the Matplotlib Gallery. Rob Johansson's IPython notebook Introduction to Matplotlib is also particularly good.
1.12 Conclusion of the Python Overview
There is, of course, much more to the language than we'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
Step73:  No matter how experienced a programmer you are, these are words to meditate on.
2. 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.)
First off, it is important to import the library. Again, we did this at the very beginning of this whole jupyter notebook, but here it is in case you've jumped straight here without running the first code line
Step74:  2.1 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
Step75:  You can pass in a second argument to array that gives the numeric type. There are a number of types listed here that your matrix can be. Some of these are aliased to single character codes. The most common ones are 'd' (double precision floating point number), 'D' (double precision complex number), and 'i' (int32). Thus,
Step76:  To build matrices, you can either use the array command with lists of lists
Step77:  You can also form empty (zero) matrices of arbitrary shape (including vectors, which Numpy treats as vectors with one row), using the zeros command
Step78:  The first argument is a tuple containing the shape of the matrix, and the second is the data type argument, which follows the same conventions as in the array command. Thus, you can make row vectors
Step79:  or column vectors
Step80:  There's also an identity command that behaves as you'd expect
Step81:  as well as a ones command.
2.2 Linspace, matrix functions, and plotting
The linspace command makes a linear array of points from a starting to an ending value.
Step82:  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.
Step83:  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,
Step84:  In conjunction with matplotlib, this is a nice way to plot things
Step85:  2.3 Matrix operations
Matrix objects act sensibly when multiplied by scalars
Step86:  as well as when you add two matrices together. (However, the matrices have to be the same shape.)
Step87:  Something that confuses Matlab users is that the times (*) operator give element-wise multiplication rather than matrix multiplication
Step88:  To get matrix multiplication, you need the dot command
Step89:  dot can also do dot products (duh!)
Step90:  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
Step91:  There's also a diag() function that takes a list or a vector and puts it along the diagonal of a square matrix.
Step92:  We'll find this useful later on.
2.4 Matrix Solvers
You can solve systems of linear equations using the solve command in the linear algebra toolbox of the numpy library
Step93:  There are a number of routines to compute eigenvalues and eigenvectors
eigvals returns the eigenvalues of a matrix
eigvalsh returns the eigenvalues of a Hermitian matrix
eig returns the eigenvalues and eigenvectors of a matrix
eigh returns the eigenvalues and eigenvectors of a Hermitian matrix.
Step94:  2.5 Example
Step95:  Let's see whether this works for our sin example from above
Step96:  Pretty close!
2.6 One-Dimensional Harmonic Oscillator using Finite Difference
Now that we've convinced ourselves that finite differences aren't a terrible approximation, let's see if we can use this to solve the one-dimensional harmonic oscillator.
We want to solve the time-independent Schrodinger equation
$$ -\frac{\hbar^2}{2m}\frac{\partial^2\psi(x)}{\partial x^2} + V(x)\psi(x) = E\psi(x)$$
for $\psi(x)$ when $V(x)=\frac{1}{2}m\omega^2x^2$ is the harmonic oscillator potential. We're going to use the standard trick to transform the differential equation into a matrix equation by multiplying both sides by $\psi^*(x)$ and integrating over $x$. This yields
$$ -\frac{\hbar}{2m}\int\psi(x)\frac{\partial^2}{\partial x^2}\psi(x)dx + \int\psi(x)V(x)\psi(x)dx = E$$
We will again use the finite difference approximation. The finite difference formula for the second derivative is
$$ y'' = \frac{y_{i+1}-2y_i+y_{i-1}}{x_{i+1}-x_{i-1}} $$
We can think of the first term in the Schrodinger equation as the overlap of the wave function $\psi(x)$ with the second derivative of the wave function $\frac{\partial^2}{\partial x^2}\psi(x)$. Given the above expression for the second derivative, we can see if we take the overlap of the states $y_1,\dots,y_n$ with the second derivative, we will only have three points where the overlap is nonzero, at $y_{i-1}$, $y_i$, and $y_{i+1}$. In matrix form, this leads to the tridiagonal Laplacian matrix, which has -2's along the diagonals, and 1's along the diagonals above and below the main diagonal.
The second term turns leads to a diagonal matrix with $V(x_i)$ on the diagonal elements. Putting all of these pieces together, we get
Step97:  We've made a couple of hacks here to get the orbitals the way we want them. First, I inserted a -1 factor before the wave functions, to fix the phase of the lowest state. The phase (sign) of a quantum wave function doesn't hold any information, only the square of the wave function does, so this doesn't really change anything. 
But the eigenfunctions as we generate them aren't properly normalized. The reason is that finite difference isn't a real basis in the quantum mechanical sense. It's a basis of Dirac δ functions at each point; we interpret the space betwen the points as being "filled" by the wave function, but the finite difference basis only has the solution being at the points themselves. We can fix this by dividing the eigenfunctions of our finite difference Hamiltonian by the square root of the spacing, and this gives properly normalized functions.
2.7 Special Functions
The solutions to the Harmonic Oscillator are supposed to be Hermite polynomials. The Wikipedia page has the HO states given by
$$\psi_n(x) = \frac{1}{\sqrt{2^n n!}}
\left(\frac{m\omega}{\pi\hbar}\right)^{1/4}
\exp\left(-\frac{m\omega x^2}{2\hbar}\right)
H_n\left(\sqrt{\frac{m\omega}{\hbar}}x\right)$$
Let's see whether they look like those. There are some special functions in the Numpy library, and some more in Scipy. Hermite Polynomials are in Numpy
Step98:  Let's compare the first function to our solution.
Step99:  The agreement is almost exact.
We can use the subplot command to put multiple comparisons in different panes on a single plot (run %matplotlib qt on a separate line first to plot in a separate window)
Step100:  Other than phase errors (which I've corrected with a little hack
Step102:  As well as Jacobi, Laguerre, Hermite polynomials, Hypergeometric functions, and many others. There's a full listing at the Scipy Special Functions Page.
2.8 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
Step103:  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.
Step104:  Since we expect the data to have an exponential decay, we can plot it using a semi-log plot.
Step105:  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) = ax + \log(A) $$
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)
Step106:  Let's see whether this curve fits the data.
Step108:  If we have more complicated functions, we may not be able to get away with fitting to a simple polynomial. Consider the following data
Step109:  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.
Step110:  Now fit to it using curve_fit
Step111:  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.
2.9 Monte Carlo, random numbers, and computing $\pi$
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 from the numpy library gives pseudorandom numbers uniformly distributed between 0 and 1
Step112:  Or, more elegantly
Step113:  np.random.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
Step114:  We can check the distribution by using the histogram feature, as shown on the help page for numpy.random.normal
Step115:  Here's an interesting use of random numbers
Step116:  The idea behind the program is that the ratio of the area of the unit circle to the square that inscribes it is $\pi/4$, so by counting the fraction of the random points in the square that are inside the circle, we get increasingly good estimates to $\pi$. 
The above code uses some higher level Numpy tricks to compute the radius of each point in a single line, to count how many radii are below one in a single line, and to filter the x,y points based on their radii. To be honest, I rarely write code like this
Step117:  If you're interested in another great method, check out Ramanujan's method. This converges so fast you really need arbitrary precision math to display enough decimal places. You can do this with the Python decimal module, if you're interested.
2.10 Numerical Integration
Integration can be hard, and sometimes it's easier to work out a definite integral using an approximation. For example, suppose we wanted to figure out the integral
Step118:  Scipy has a numerical integration routine quad (since sometimes numerical integration is called quadrature), that we can use for this
Step119:  The first number in the tuple is the result, the second number is an estimate of the absolute error in the result.
There are also 2d and 3d numerical integrators in Scipy. See the docs for more information.
2.11 Fast Fourier Transform and Signal Processing
Very often we want to use FFT techniques to help obtain the signal from noisy data. Scipy has several different options for this.
Step121:  There are additional signal processing routines in Scipy (e.g. splines, filtering) that you can read about here.
3. Intermediate Python
3.1 Parsing data output
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 CSV (comma separated values) format, a format that originally came from Microsoft Excel, and is increasingly used as a data interchange format in big data applications. How would we parse that?
Step122:  This is a giant string. If we use splitlines(), we see that a list is created where line gets separated into a string
Step123:  Splitting is a big concept in text processing. We used splitlines() here, and next we'll use the more general .split(",") function below to split each line into comma-delimited words.
We now want to do three things
Step124:  What does split() do?
Step125:  Since the data is now in a list of lines, we can iterate over it, splitting up data where we see a comma
Step126:  We need to add these results at each step to a list
Step127:  Let's examine what we just did
Step128:  Hartrees (what most quantum chemistry programs use by default) are really stupid units. We really want this in kcal/mol or eV or something we use. So let's quickly replot this in terms of eV above the minimum energy, which will give us a much more useful plot
Step129:  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.
3.2 Reading in data files
Let's take a look at a perhaps easier approach to a common problem -- you have a data file with some header info and comma-delimited values and you want the data so you can start doing stuff with it. Let's use numpy's genfromtxt()
Step130:  That was easy! Why didn't we only learn that? Because not every data set is "nice" like that. Better to have some tools for when things aren't working how you'd like them to be.  That being said, much data coming from scientific equipment and computational tools can be cast into a format that can be read in through genfromtxt(). For larger data sets, the library pandas might be helpful.
3.3 More Sophisticated String Formatting and Processing
Strings are a big deal in most modern languages, and hopefully the previous sections helped underscore how versatile Python's string processing techniques are. We will continue this topic in this section.
We can print out lines in Python using the print command.
Step131:  In IPython we don't even need the print command, since it will display the last expression not assigned to a variable.
Step132:  print even converts some arguments to strings for us
Step133:  As versatile as this is, you typically need more freedom over the data you print out. For example, what if we want to print a bunch of data to exactly 4 decimal places? We can do this using formatted strings.
Formatted strings share a syntax with the C printf statement. We make a string that has some funny format characters in it, and then pass a bunch of variables into the string that fill out those characters in different ways.
For example,
Step134:  We use a percent sign in two different ways here. First, the format character itself starts with a percent sign. %d or %i are for integers, %f is for floats, %e is for numbers in exponential formats. All of the numbers can take number immediately after the percent that specifies the total spaces used to print the number. Formats with a decimal can take an additional number after a dot . to specify the number of decimal places to print.
The other use of the percent sign is after the string, to pipe a set of variables in. You can pass in multiple variables (if your formatting string supports it) by putting a tuple after the percent. Thus,
Step136:  This is a simple formatting structure that will satisfy most of your string formatting needs. More information on different format symbols is available in the string formatting part of the standard docs.
It's worth noting that more complicated string formatting methods are in development, but I prefer this system due to its simplicity and its similarity to C formatting strings.
Recall we discussed multiline strings. We can put format characters in these as well, and fill them with the percent sign as before.
Step138:  The problem with a long block of text like this is that it's often hard to keep track of what all of the variables are supposed to stand for. There's an alternate format where you can pass a dictionary into the formatted string, and give a little bit more information to the formatted string itself. This method looks like
Step140:  By providing a little bit more information, you're less likely to make mistakes, like referring to your customer as "alien attack".
As a scientist, you're less likely to be sending bulk mailings to a bunch of customers. But these are great methods for generating and submitting lots of similar runs, say scanning a bunch of different structures to find the optimal configuration for something.
For example, you can use the following template for NWChem input files
Step142:  If you want to submit a sequence of runs to a computer somewhere, it's pretty easy to put together a little script, maybe even with some more string formatting in it
Step143:  This is a very bad geometry for a water molecule, and it would be silly to run so many geometry optimizations of structures that are guaranteed to converge to the same single geometry, but you get the idea of how you can run vast numbers of simulations with a technique like this.
We used the enumerate function to loop over both the indices and the items of a sequence, which is valuable when you want a clean way of getting both. enumerate is roughly equivalent to
Step144:  Although enumerate uses generators (see below) so that it doesn't have to create a big list, which makes it faster for really long sequenes.
3.4 Optional arguments of a function
You will recall that the linspace function can take either two arguments (for the starting and ending points)
Step145:  or it can take three arguments, for the starting point, the ending point, and the number of points
Step146:  You can also pass in keywords to exclude the endpoint
Step147:  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
Step148:  We can add an optional argument by specifying a default value in the argument list
Step149:  This gives exactly the same result if we don't specify anything
Step150:  But also let's us override the default value with a third argument
Step151:  We can add arbitrary keyword arguments to the function definition by putting a keyword argument **kwargs handle in
Step152:  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
Step153:  Note that we have defined a few new things you haven't seen before
Step154:  3.5 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
Step155:  You can also put some boolean testing into the construct
Step156:  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
Step157:  We can always turn an iterator into a list using the list command
Step158:  There's a special syntax called a generator expression that looks a lot like a list comprehension
Step159:  3.6 Factory Functions
A factory function is a function that returns a function. They have the fancy name lexical closure, which makes you sound really intelligent in front of your CS friends. But, despite the arcane names, factory functions can play a very practical role.
Suppose you want the Gaussian function centered at 0.5, with height 99 and width 1.0. You could write a general function.
Step160:  But what if you need a function with only one argument, like f(x) rather than f(x,y,z,...)? You can do this with Factory Functions
Step162:  Everything in Python is an object, including functions. This means that functions can be returned by other functions. (They can also be passed into other functions, which is also useful, but a topic for another discussion.) In the gauss_maker example, the g function that is output "remembers" the A, a, x0 values it was constructed with, since they're all stored in the local memory space (this is what the lexical closure really refers to) of that function.
Factories are one of the more important of the Software Design Patterns, which are a set of guidelines to follow to make high-quality, portable, readable, stable software. It's beyond the scope of the current work to go more into either factories or design patterns, but I thought I would mention them for people interested in software design.
3.7 Serialization
Step163:  Your data sits in something that looks like a Python dictionary, and in a single line of code, you can load it into a Python dictionary for use later.
In the same way, you can, with a single line of code, put a bunch of variables into a dictionary, and then output to a file using json
Step164:  3.8 Functional programming
Functional programming is a very broad subject. The idea is to have a series of functions, each of which generates a new data structure from an input, without changing the input structure at all. By not modifying the input structure (something that is called not having side effects), many guarantees can be made about how independent the processes are, which can help parallelization and guarantees of program accuracy. There is a Python Functional Programming HOWTO in the standard docs that goes into more details on functional programming. I just wanted to touch on a few of the most important ideas here.
There is an operator module that has function versions of most of the Python operators. For example
Step165:  These are useful building blocks for functional programming.
The lambda operator allows us to build anonymous functions, which are simply functions that aren't defined by a normal def statement with a name. For example, a function that doubles the input is
Step166:  We could also write this as
Step167:  And assign it to a function separately
Step168:  lambda is particularly convenient (as we'll see below) in passing simple functions as arguments to other functions.
map is a way to repeatedly apply a function to a list
Step169:  reduce is a way to repeatedly apply a function to the first two items of the list. There already is a sum function in Python that is a reduction
Step170:  We can use reduce to define an analogous prod function
Step171:  3.9 Object Oriented Programming
We've seen a lot of examples of objects in Python. We create a string object with quote marks
Step172:  and we have a bunch of methods we can use on the object
Step174:  Object oriented programming simply gives you the tools to define objects and methods for yourself. It's useful anytime you want to keep some data (like the characters in the string) tightly coupled to the functions that act on the data (length, split, startswith, etc.).
As an example, we're going to bundle the functions we did to make the 1d harmonic oscillator eigenfunctions with arbitrary potentials, so we can pass in a function defining that potential, some additional specifications, and get out something that can plot the orbitals, as well as do other things with them, if desired.
Step175:  The init() function specifies what operations go on when the object is created. The self argument is the object itself, and we don't pass it in. The only required argument is the function that defines the QM potential. We can also specify additional arguments that define the numerical grid that we're going to use for the calculation.
For example, to do an infinite square well potential, we have a function that is 0 everywhere. We don't have to specify the barriers, since we'll only define the potential in the well, which means that it can't be defined anywhere else.
Step176:  We can similarly redefine the Harmonic Oscillator potential.
Step177:  Let's define a finite well potential
Step178:  A triangular well
Step179:  Or we can combine the two, making something like a semiconductor quantum well with a top gate
Step180:  There's a lot of philosophy behind object oriented programming. Since I'm trying to focus on just the basics here, I won't go into them, but the internet is full of lots of resources on OO programming and theory. The best of this is contained in the Design Patterns book, which I highly recommend.
4. Speeding Python
Step181:  The little % sign that we have in front of the timeit call is an example of an IPython magic function, which we don't have time to go into here, but it's just some little extra mojo that IPython adds to the functions to make it run better in the IPython environment. You can read more about it in the IPython tutorial.
In any case, the timeit function runs 3 loops, and tells us that it took on the average of 583 ns to compute 20!. In contrast
Step182:  the factorial function we wrote is about a factor of 10 slower. This is because the built-in factorial function is written in C code and called from Python, and the version we wrote is written in plain old Python. A Python program has a lot of stuff in it that make it nice to interact with, but all that friendliness slows down the code. In contrast, the C code is less friendly but more efficient. If you want speed with as little effort as possible, write your code in an easy to program language like Python, but dump the slow parts into a faster language like C, and call it from Python. We'll go through some tricks to do this in this section.
4.2 Profiling
Profiling complements what timeit does by splitting the overall timing into the time spent in each function. It can give us a better understanding of what our program is really spending its time on.
Suppose we want to create a list of even numbers. Our first effort yields this
Step183:  Is this code fast enough? We find out by running the Python profiler on a longer run
Step184:  This looks okay, 0.05 seconds isn't a huge amount of time, but looking at the profiling shows that the append function is taking almost 20% of the time. Can we do better? Let's try a list comprehension.
Step185:  By removing a small part of the code using a list comprehension, we've doubled the overall speed of the code! 
It seems like range is taking a long time, still. Can we get rid of it? We can, using the xrange generator
Step187:  This is where profiling can be useful. Our code now runs 3x faster by making trivial changes. We wouldn't have thought to look in these places had we not had access to easy profiling. Imagine what you would find in more complicated programs.
4.3 Other Ways to Speed Python
When we compared the fact and factorial functions, above, we noted that C routines are often faster because they're more streamlined. Once we've determined that one routine is a bottleneck for the performance of a program, we can replace it with a faster version by writing it in C. This is called extending Python, and there's a good section in the standard documents. This can be a tedious process if you have many different routines to convert. Fortunately, there are several other options.
Swig (the simplified wrapper and interface generator) is a method to generate binding not only for Python but also for Matlab, Perl, Ruby, and other scripting languages. Swig can scan the header files of a C project and generate Python binding for it. Using Swig is substantially easier than writing the routines in C.
Cython is a C-extension language. You can start by compiling a Python routine into a shared object libraries that can be imported into faster versions of the routines. You can then add additional static typing and make other restrictions to further speed the code. Cython is generally easier than using Swig.
PyPy is the easiest way of obtaining fast code. PyPy compiles Python to a subset of the Python language called RPython that can be efficiently compiled and optimized. Over a wide range of tests, PyPy is roughly 6 times faster than the standard Python Distribution.
4.4 Fun
Step188:  You might think that Python is a bad choice for something like this, but, in terms of time, it really doesn't take long
<ASSISTANT_TASK:>
Python Code:
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
Explanation: A Crash Course in Python for Scientists
Rick Muller, Sandia National Laboratories
version 0.62, Updated Dec 15, 2016 by Ryan Smith, Cal State East Bay
version 0.63, Updated Oct 2017 by Ryan Smith, Cal State East Bay
Using Python 3.5.2 | Anaconda 4.1.1 
This work is licensed under a Creative Commons Attribution-ShareAlike 3.0 Unported License.
If skipping to other sections, it's a good idea to run this first:
(If you're using an IPython (jupyter) notebook, or otherwise using notebook file, you hit shift-Enter to evaluate a cell.)
End of explanation
2+2
(50-5*6)/4
Explanation: Table of Contents
0. Preliminary
0.1 Why Python?
0.2 What You Need to Install
1. Python Overview
1.1 Using Python as a Calculator
1.2 Strings
1.3 Lists
1.4 Iteration, Indentation, and Blocks
1.5 Slicing
1.6 Booleans and Truth Testing
1.7 Code Example: The Fibonacci Sequence
1.8 Functions
1.9 Recursion and Factorials
1.10 Two More Data Structures: Tuples and Dictionaries
1.11 Plotting with Matplotlib
1.12 Conclusion of the Python Overview
2. Numpy and Scipy
2.1 Making vectors and matrices
2.2 Linspace, matrix functions, and plotting
2.3 Matrix operations
2.4 Matrix Solvers
2.5 Example: Finite Differences
2.6 One-Dimensional Harmonic Oscillator using Finite Difference
2.7 Special Functions
2.8 Least squares fitting
2.9 Monte Carlo, random numbers, and computing $\pi$
2.10 Numerical Integration
2.11 Fast Fourier Transform and Signal Processing
3. Intermediate Python
3.1 Parsing data output
3.2 Reading in data files
3.3 More Sophisticated String Formatting and Processing
3.4 Optional arguments of a function
3.5 List Comprehensions and Generators
3.6 Factory Functions
3.7 Serialization: Save it for later
3.8 Functional programming
3.9 Object Oriented Programming
4. Speeding Python: Timeit, Profiling, Cython, SWIG, and PyPy
4.1 Timeit
4.2 Profiling
4.3 Other Ways to Speed Python
4.4 Fun: Finding Primes
5. References
6. Acknowledgements
0. Preliminary
0.1 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 your 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 IPython 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:
Rob Johansson's excellent notebooks, including Scientific Computing with Python and Computational Quantum Physics with QuTiP lectures;
XKCD style graphs in matplotlib;
A collection of Notebooks for using IPython effectively
A gallery of interesting IPython Notebooks
I find IPython 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. I find myself endlessly sweeping the IPython subreddit hoping someone will post a new notebook. 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, Matplotlib, and IPython development, as well as my own experience in using Python almost every day of this time.
IPython notebooks are now called Jupyter notebooks.
0.2 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. We're now (2016) past the halfway point, and people are moving to python 3.
These notes are written with Python 3 in mind.
If you are new to python, try installing Anaconda Python 3.5 (supported by Continuum) and you will automatically have all libraries installed with your distribution. These notes assume you have a Python distribution that includes:
Python version 3;
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.
Here are some other options for various ways to run python:
Continuum supports a bundled, multiplatform Python package called Anaconda 
Entought Python Distribution, also known as EPD. You can either purchase a license to use EPD, or there is also a free version that you can download and install.
Linux Most distributions have an installation manager. Redhat has yum, Ubuntu has apt-get. To my knowledge, all of these packages should be available through those installers.
Mac I use Macports, which has up-to-date versions of all of these packages.
Windows The PythonXY package has everything you need: install the package, then go to Start > PythonXY > Command Prompts > IPython notebook server.
Cloud This notebook is currently running on the IPython notebook viewer, which allows the notebook to be viewed but not interactively.
1. 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.
1.1 Using Python as a Calculator
Many of the things I used to use a calculator for, I now use Python for:
End of explanation
2**3
7/33
Explanation: If you want to raise to an exponent, "^" is reserved for something else, so we use ** instead
End of explanation
from math import sqrt
sqrt(81)
Explanation: There used to be gotchas in division in python 2, like C or Fortran integer division, where division truncates the remainder and returns an integer. In version 3, Python returns a floating point number. If for some reason you are using Python 2, you can fix this by importing the module from the future features:
from __future__ import division
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 numpy module containing many useful functions. To access, say, the square root function, you can either first import the sqrt function from the math library:
End of explanation
np.sqrt(81)
Explanation: I prefer doing things with numpy though.. we already imported numpy as np in the first command box above, so we can just call it... I prefer things this way, so that code is always easily readable: you will know from what library a function is being called immediately.
End of explanation
width = 20
length = 30
area = length*width
area
Explanation: You can define variables using the equals (=) sign:
End of explanation
volume
Explanation: If you try to access a variable that you haven't yet defined, you get an error:
End of explanation
depth = 10
volume = area*depth
volume
Explanation: and you need to define it:
End of explanation
np.exp(-1)
Explanation: A few examples of common functions. Look up numpy library -- it has most mathematical functions you could want.
End of explanation
return = 0
Explanation: 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:
End of explanation
'Hello, World!'
Explanation: 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.
1.2 Strings
Strings are lists of printable characters, and can be defined using either single quotes
End of explanation
"Hello, World!"
Explanation: or double quotes
End of explanation
"He's a Rebel"
'She asked, "How are you today?"'
Explanation: But not both at the same time, unless you want one of the symbols to be part of the string.
End of explanation
greeting = "Hello, World!"
Explanation: Just like the other two data objects we're familiar with (ints and floats), you can assign a string to a variable
End of explanation
print(greeting)
Explanation: The print statement is often used for printing character strings:
End of explanation
print("The area is ",area)
Explanation: But it can also print data types other than strings:
End of explanation
statement = "Hello," + "World!"
print(statement)
Explanation: 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:
End of explanation
statement = "Hello, " + "World!"
print(statement)
Explanation: Don't forget the space between the strings, if you want one there.
End of explanation
print( "This " + "is " + "a " + "longer " + "statement.")
Explanation: You can use + to concatenate multiple strings in a single statement:
End of explanation
days_of_the_week = ["Sunday","Monday","Tuesday","Wednesday","Thursday","Friday","Saturday"]
Explanation: 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.
1.3 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.
End of explanation
days_of_the_week[2]
Explanation: You can access members of the list using the index of that item:
End of explanation
days_of_the_week[-1]
Explanation: 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:
End of explanation
languages = ["Fortran","C","C++"]
languages.append("Python")
print(languages)
Explanation: You can add additional items to the list using the .append() command:
End of explanation
list(range(10))
Explanation: The range() command is a convenient way to make sequential lists of numbers:
End of explanation
list(range(2,8))
Explanation: 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)
End of explanation
evens = list(range(0,20,2))
evens
evens[3]
Explanation: 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:
End of explanation
["Today",7,99.3,""]
Explanation: Lists do not have to hold the same data type. For example,
End of explanation
help(len)
len(evens)
Explanation: 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:
End of explanation
for day in days_of_the_week:
    print(day)
Explanation: 1.4 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:
End of explanation
for day in days_of_the_week:
    statement = "Today is " + day
    print(statement)
Explanation: 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:
End of explanation
for i in range(20):
    print("The square of ",i," is ",i*i)
Explanation: The range() command is particularly useful with the for statement to execute loops of a specified length:
End of explanation
for letter in "Sunday":
    print(letter)
Explanation: 1.5 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:
End of explanation
days_of_the_week[0]
Explanation: 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:
End of explanation
days_of_the_week[0:2]
Explanation: If we want the list containing the first two elements of a list, we can do this via
End of explanation
days_of_the_week[:2]
Explanation: or simply
End of explanation
days_of_the_week[-2:]
Explanation: If we want the last items of the list, we can do this with negative slicing:
End of explanation
workdays = days_of_the_week[1:5]
print(workdays)
Explanation: which is somewhat logically consistent with negative indices accessing the last elements of the list.
You can do:
End of explanation
day = "Saturday"
abbreviation = day[:3]
print(abbreviation)
Explanation: Since strings are sequences, you can also do this to them:
End of explanation
numbers = list(range(0,21))
evens = numbers[6::2]
evens
Explanation: 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):
End of explanation
if day == "Sunday":
    print("Sleep in")
else:
    print("Go to work")
Explanation: Note that in this example we omitted the first few arguments, so that the slice started at 6, went to the end of the list, and took every second element, to generate the list of even numbers lup to 20 (the last element in the original list).
1.6 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:
End of explanation
day == "Sunday"
Explanation: (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
End of explanation
1 == 2
50 == 2*25
3 < 3.14159
1 == 1.0
1 != 0
1 <= 2
1 >= 1
Explanation: 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".
Try setting the day equal to "Sunday" and then running the above if/else statement. Did it work as you thought it would?
You can compare any data types in Python:
End of explanation
1 is 1.0
Explanation: 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:
End of explanation
type(1)
type(1.0)
Explanation: Why is 1 not the same as 1.0? Different data type. You can check the data type:
End of explanation
[1,2,3] == [1,2,4]
Explanation: We can do boolean tests on lists as well:
End of explanation
hours = 5
0 < hours < 24
Explanation: Finally, note that you can also string multiple comparisons together, which can result in very intuitive tests:
End of explanation
if day == "Sunday":
    print ("Sleep in")
elif day == "Saturday":
    print ("Do chores")
else:
    print ("Go to work")
Explanation: If statements can have elif parts ("else if"), in addition to if/else parts. For example:
End of explanation
for day in days_of_the_week:
    statement = "On " + day + ":" 
    print (statement)
    if day == "Sunday":
        print ("   Sleep in")
    elif day == "Saturday":
        print ("   Do chores")
    else:
        print ("   Go to work")
Explanation: Of course we can combine if statements with for loops, to make a snippet that is almost interesting:
End of explanation
bool(1)
bool(0)
bool(["This "," is "," a "," list"])
Explanation: 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.
End of explanation
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)
Explanation: 1.7 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.
End of explanation
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
Explanation: 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!
1.8 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:
End of explanation
fibonacci(2)
fibonacci(12)
Explanation: We can now call fibonacci() for different sequence_lengths:
End of explanation
help(fibonacci)
Explanation: 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:
End of explanation
from math import factorial
help(factorial)
Explanation: 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.
1.9 Recursion and Factorials
Functions can also call themselves, something that is often called recursion. We're going to experiment with recursion by computing the factorial function. The factorial is defined for a positive integer n as
$$ n! = n(n-1)(n-2)\cdots 1 $$
First, note that we don't need to write a function at all, since this is a function built into the standard math library. Let's use the help function to find out about it:
End of explanation
factorial(20)
Explanation: This is clearly what we want.
End of explanation
def fact(n):
    if n <= 0:
        return 1
    return n*fact(n-1)
fact(20)
Explanation: However, if we did want to write a function ourselves, we could do recursively by noting that
$$ n! = n(n-1)!$$
The program then looks something like:
End of explanation
t = (1,2,'hi',9.0)
t
Explanation: Recursion can be very elegant, and can lead to very simple programs.
1.10 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:
End of explanation
t[1]
Explanation: Tuples are like lists, in that you can access the elements using indices:
End of explanation
t.append(7)
t[1]=77
Explanation: However, tuples are immutable, you can't append to them or change the elements of them:
End of explanation
('Bob',0.0,21.0)
Explanation: 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:
End of explanation
positions = [
             ('Bob',0.0,21.0),
             ('Cat',2.5,13.1),
             ('Dog',33.0,1.2)
             ]
Explanation: 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:
End of explanation
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)
Explanation: 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:
End of explanation
x,y = 1,2
y,x = x,y
x,y
Explanation: Here we did two things with tuples you haven't seen before. First, we unpacked an object into a set of named variables using tuple assignment:
&gt;&gt;&gt; name,x,y = obj
We also returned multiple values (minx,miny), which were then assigned to two other variables (x,y), again by tuple assignment. This makes what would have been complicated code in C++ rather simple.
Tuple assignment is also a convenient way to swap variables:
End of explanation
mylist = [1,2,9,21]
Explanation: Dictionaries are an object called "mappings" or "associative arrays" in other languages. Whereas a list associates an integer index with a set of objects:
End of explanation
ages = {"Rick": 46, "Bob": 86, "Fred": 21}
print("Rick's age is ",ages["Rick"])
Explanation: 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 {}:
End of explanation
dict(Rick=46,Bob=86,Fred=20)
Explanation: There's also a convenient way to create dictionaries without having to quote the keys.
End of explanation
len(t)
len(ages)
Explanation: Notice in either case you are not choosing the ordering -- it is automagically grouped alphabetically.
The len() command works on both tuples and dictionaries:
End of explanation
import matplotlib.pyplot as plt
%matplotlib inline
Explanation: 1.11 Plotting with Matplotlib
We can generally understand trends in data by using a plotting program to chart it. Python has a wonderful plotting library called Matplotlib. The Jupyter notebook interface we are using for these notes has that functionality built in.
First off, it is important to import the library. We did this at the very beginning of this whole jupyter notebook, but here it is in case you've jumped straight here without running the first code line:
End of explanation
fibs = fibonacci(10)
Explanation: The %matplotlib inline command makes it so plots are within this notebook. To plot to a separate window, use instead:
%matplotlib qt
As an example of plotting, we have looked at two different functions, the Fibonacci function, and the factorial function, both of which grow faster than polynomially. Which one grows the fastest? Let's plot them. First, let's generate the Fibonacci sequence of length 20:
End of explanation
facts = []
for i in range(10):
    facts.append(factorial(i))
Explanation: Next lets generate the factorials.
End of explanation
plt.plot(facts,'-ob',label="factorial")
plt.plot(fibs,'-dg',label="Fibonacci")
plt.xlabel("n")
plt.legend()
Explanation: Now we use the Matplotlib function plot to compare the two.
End of explanation
plt.semilogy(facts,label="factorial")
plt.semilogy(fibs,label="Fibonacci")
plt.xlabel("n")
plt.legend()
Explanation: The factorial function grows much faster. In fact, you can't even see the Fibonacci sequence. It's not entirely surprising: a function where we multiply by n each iteration is bound to grow faster than one where we add (roughly) n each iteration.
Let's plot these on a semilog plot so we can see them both a little more clearly:
End of explanation
import this
Explanation: There are many more things you can do with Matplotlib. We'll be looking at some of them in the sections to come. In the meantime, if you want an idea of the different things you can do, look at the Matplotlib Gallery. Rob Johansson's IPython notebook Introduction to Matplotlib is also particularly good.
1.12 Conclusion of the Python Overview
There is, of course, much more to the language than we'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:
End of explanation
import numpy as np
Explanation: No matter how experienced a programmer you are, these are words to meditate on.
2. 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.)
First off, it is important to import the library. Again, we did this at the very beginning of this whole jupyter notebook, but here it is in case you've jumped straight here without running the first code line:
End of explanation
np.array([1,2,3,4,5,6])
Explanation: 2.1 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:
End of explanation
np.array([1,2,3,4,5,6],'d')
np.array([1,2,3,4,5,6],'D')
np.array([1,2,3,4,5,6],'i')
Explanation: You can pass in a second argument to array that gives the numeric type. There are a number of types listed here that your matrix can be. Some of these are aliased to single character codes. The most common ones are 'd' (double precision floating point number), 'D' (double precision complex number), and 'i' (int32). Thus,
End of explanation
np.array([[0,1],[1,0]],'d')
Explanation: To build matrices, you can either use the array command with lists of lists:
End of explanation
np.zeros((3,3),'d')
Explanation: You can also form empty (zero) matrices of arbitrary shape (including vectors, which Numpy treats as vectors with one row), using the zeros command:
End of explanation
np.zeros(3,'d')
np.zeros((1,3),'d')
Explanation: The first argument is a tuple containing the shape of the matrix, and the second is the data type argument, which follows the same conventions as in the array command. Thus, you can make row vectors:
End of explanation
np.zeros((3,1),'d')
Explanation: or column vectors:
End of explanation
np.identity(4,'d')
Explanation: There's also an identity command that behaves as you'd expect:
End of explanation
np.linspace(0,1)
Explanation: as well as a ones command.
2.2 Linspace, matrix functions, and plotting
The linspace command makes a linear array of points from a starting to an ending value.
End of explanation
np.linspace(0,1,11)
Explanation: 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.
End of explanation
x = np.linspace(0,2*np.pi)
np.sin(x)
Explanation: 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,
End of explanation
plt.plot(x,np.sin(x))
plt.show()
Explanation: In conjunction with matplotlib, this is a nice way to plot things:
End of explanation
0.125*np.identity(3,'d')
Explanation: 2.3 Matrix operations
Matrix objects act sensibly when multiplied by scalars:
End of explanation
np.identity(2,'d') + np.array([[1,1],[1,2]])
Explanation: as well as when you add two matrices together. (However, the matrices have to be the same shape.)
End of explanation
np.identity(2)*np.ones((2,2))
Explanation: Something that confuses Matlab users is that the times (*) operator give element-wise multiplication rather than matrix multiplication:
End of explanation
np.dot(np.identity(2),np.ones((2,2)))
Explanation: To get matrix multiplication, you need the dot command:
End of explanation
v = np.array([3,4],'d')
np.sqrt(np.dot(v,v))
Explanation: dot can also do dot products (duh!):
End of explanation
m = np.array([[1,2],[3,4]])
m.T
Explanation: 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:
End of explanation
np.diag([1,2,3,4,5])
Explanation: There's also a diag() function that takes a list or a vector and puts it along the diagonal of a square matrix.
End of explanation
A = np.array([[1,1,1],[0,2,5],[2,5,-1]])
b = np.array([6,-4,27])
np.linalg.solve(A,b)
Explanation: We'll find this useful later on.
2.4 Matrix Solvers
You can solve systems of linear equations using the solve command in the linear algebra toolbox of the numpy library:
End of explanation
A = np.array([[13,-4],[-4,7]],'d')
np.linalg.eigvalsh(A)
np.linalg.eigh(A)
Explanation: There are a number of routines to compute eigenvalues and eigenvectors
eigvals returns the eigenvalues of a matrix
eigvalsh returns the eigenvalues of a Hermitian matrix
eig returns the eigenvalues and eigenvectors of a matrix
eigh returns the eigenvalues and eigenvectors of a Hermitian matrix.
End of explanation
def nderiv(y,x):
    "Finite difference derivative of the function f"
    n = len(y)
    d = np.zeros(n,'d') # assume double
    # Use centered differences for the interior points, one-sided differences for the ends
    for i in range(1,n-1):
        d[i] = (y[i+1]-y[i])/(x[i+1]-x[i])
    d[0] = (y[1]-y[0])/(x[1]-x[0])
    d[n-1] = (y[n-1]-y[n-2])/(x[n-1]-x[n-2])
    return d
Explanation: 2.5 Example: Finite Differences
Now that we have these tools in our toolbox, we can start to do some cool stuff with it. Many of the equations we want to solve in Physics involve differential equations. We want to be able to compute the derivative of functions:
$$ y' = \frac{y(x+h)-y(x)}{h} $$
by discretizing the function $y(x)$ on an evenly spaced set of points $x_0, x_1, \dots, x_n$, yielding $y_0, y_1, \dots, y_n$. Using the discretization, we can approximate the derivative by
$$ y_i' \approx \frac{y_{i+1}-y_{i-1}}{x_{i+1}-x_{i-1}} $$
We can write a derivative function in Python via
End of explanation
x = np.linspace(0,2*np.pi)
dsin = nderiv(np.sin(x),x)
plt.plot(x,dsin,label='numerical')
plt.plot(x,np.cos(x),label='analytical')
plt.title("Comparison of numerical and analytical derivatives of sin(x)")
plt.legend()
Explanation: Let's see whether this works for our sin example from above:
End of explanation
def Laplacian(x):
    h = x[1]-x[0] # assume uniformly spaced points
    n = len(x)
    M = -2*np.identity(n,'d')
    for i in range(1,n):
        M[i,i-1] = M[i-1,i] = 1
    return M/h**2
x = np.linspace(-3,3)
m = 1.0
ohm = 1.0
T = (-0.5/m)*Laplacian(x)
V = 0.5*(ohm**2)*(x**2)
H =  T + np.diag(V)
E,U = np.linalg.eigh(H)
h = x[1]-x[0]
# Plot the Harmonic potential
plt.plot(x,V,color='k')
for i in range(4):
    # For each of the first few solutions, plot the energy level:
    plt.axhline(y=E[i],color='k',ls=":")
    # as well as the eigenfunction, displaced by the energy level so they don't
    # all pile up on each other:
    plt.plot(x,-U[:,i]/np.sqrt(h)+E[i])
plt.title("Eigenfunctions of the Quantum Harmonic Oscillator")
plt.xlabel("Displacement (bohr)")
plt.ylabel("Energy (hartree)")
Explanation: Pretty close!
2.6 One-Dimensional Harmonic Oscillator using Finite Difference
Now that we've convinced ourselves that finite differences aren't a terrible approximation, let's see if we can use this to solve the one-dimensional harmonic oscillator.
We want to solve the time-independent Schrodinger equation
$$ -\frac{\hbar^2}{2m}\frac{\partial^2\psi(x)}{\partial x^2} + V(x)\psi(x) = E\psi(x)$$
for $\psi(x)$ when $V(x)=\frac{1}{2}m\omega^2x^2$ is the harmonic oscillator potential. We're going to use the standard trick to transform the differential equation into a matrix equation by multiplying both sides by $\psi^*(x)$ and integrating over $x$. This yields
$$ -\frac{\hbar}{2m}\int\psi(x)\frac{\partial^2}{\partial x^2}\psi(x)dx + \int\psi(x)V(x)\psi(x)dx = E$$
We will again use the finite difference approximation. The finite difference formula for the second derivative is
$$ y'' = \frac{y_{i+1}-2y_i+y_{i-1}}{x_{i+1}-x_{i-1}} $$
We can think of the first term in the Schrodinger equation as the overlap of the wave function $\psi(x)$ with the second derivative of the wave function $\frac{\partial^2}{\partial x^2}\psi(x)$. Given the above expression for the second derivative, we can see if we take the overlap of the states $y_1,\dots,y_n$ with the second derivative, we will only have three points where the overlap is nonzero, at $y_{i-1}$, $y_i$, and $y_{i+1}$. In matrix form, this leads to the tridiagonal Laplacian matrix, which has -2's along the diagonals, and 1's along the diagonals above and below the main diagonal.
The second term turns leads to a diagonal matrix with $V(x_i)$ on the diagonal elements. Putting all of these pieces together, we get:
End of explanation
from numpy.polynomial.hermite import Hermite
def ho_evec(x,n,m,ohm):
    vec = [0]*9
    vec[n] = 1
    Hn = Hermite(vec)
    return (1/np.sqrt(2**n*factorial(n)))*pow(m*ohm/np.pi,0.25)*np.exp(-0.5*m*ohm*x**2)*Hn(x*np.sqrt(m*ohm))
Explanation: We've made a couple of hacks here to get the orbitals the way we want them. First, I inserted a -1 factor before the wave functions, to fix the phase of the lowest state. The phase (sign) of a quantum wave function doesn't hold any information, only the square of the wave function does, so this doesn't really change anything. 
But the eigenfunctions as we generate them aren't properly normalized. The reason is that finite difference isn't a real basis in the quantum mechanical sense. It's a basis of Dirac δ functions at each point; we interpret the space betwen the points as being "filled" by the wave function, but the finite difference basis only has the solution being at the points themselves. We can fix this by dividing the eigenfunctions of our finite difference Hamiltonian by the square root of the spacing, and this gives properly normalized functions.
2.7 Special Functions
The solutions to the Harmonic Oscillator are supposed to be Hermite polynomials. The Wikipedia page has the HO states given by
$$\psi_n(x) = \frac{1}{\sqrt{2^n n!}}
\left(\frac{m\omega}{\pi\hbar}\right)^{1/4}
\exp\left(-\frac{m\omega x^2}{2\hbar}\right)
H_n\left(\sqrt{\frac{m\omega}{\hbar}}x\right)$$
Let's see whether they look like those. There are some special functions in the Numpy library, and some more in Scipy. Hermite Polynomials are in Numpy:
End of explanation
plt.plot(x,ho_evec(x,0,1,1),label="Analytic")
plt.plot(x,-U[:,0]/np.sqrt(h),label="Numeric")
plt.xlabel('x (bohr)')
plt.ylabel(r'$\psi(x)$')
plt.title("Comparison of numeric and analytic solutions to the Harmonic Oscillator")
plt.legend()
Explanation: Let's compare the first function to our solution.
End of explanation
phase_correction = [-1,1,1,-1,-1,1]
for i in range(6):
    plt.subplot(2,3,i+1)
    plt.plot(x,ho_evec(x,i,1,1),label="Analytic")
    plt.plot(x,phase_correction[i]*U[:,i]/np.sqrt(h),label="Numeric")
Explanation: The agreement is almost exact.
We can use the subplot command to put multiple comparisons in different panes on a single plot (run %matplotlib qt on a separate line first to plot in a separate window):
End of explanation
from scipy.special import airy,jn,eval_chebyt,eval_legendre
plt.subplot(2,2,1)
x = np.linspace(-1,1)
Ai,Aip,Bi,Bip = airy(x)
plt.plot(x,Ai)
plt.plot(x,Aip)
plt.plot(x,Bi)
plt.plot(x,Bip)
plt.title("Airy functions")
plt.subplot(2,2,2)
x = np.linspace(0,10)
for i in range(4):
    plt.plot(x,jn(i,x))
plt.title("Bessel functions")
plt.subplot(2,2,3)
x = np.linspace(-1,1)
for i in range(6):
    plt.plot(x,eval_chebyt(i,x))
plt.title("Chebyshev polynomials of the first kind")
plt.subplot(2,2,4)
x = np.linspace(-1,1)
for i in range(6):
    plt.plot(x,eval_legendre(i,x))
plt.title("Legendre polynomials")
# plt.tight_layout()
plt.show()
Explanation: Other than phase errors (which I've corrected with a little hack: can you find it?), the agreement is pretty good, although it gets worse the higher in energy we get, in part because we used only 50 points.
The Scipy module has many more special functions:
End of explanation
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
Explanation: As well as Jacobi, Laguerre, Hermite polynomials, Hypergeometric functions, and many others. There's a full listing at the Scipy Special Functions Page.
2.8 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:
End of explanation
data = []
for line in raw_data.splitlines():
    words = line.split(',')
    data.append(list(map(float,words)))
data = np.array(data)
plt.title("Raw Data")
plt.xlabel("Distance")
plt.plot(data[:,0],data[:,1],'bo')
Explanation: 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.
End of explanation
plt.title("Raw Data")
plt.xlabel("Distance")
plt.semilogy(data[:,0],data[:,1],'bo')
Explanation: Since we expect the data to have an exponential decay, we can plot it using a semi-log plot.
End of explanation
params = np.polyfit(data[:,0],np.log(data[:,1]),1)
a = params[0] # the coefficient of x**1
logA = params[1] # the coefficient of x**0
# plot if curious:
# plt.plot(data[:,0],np.log(data[:,1]),'bo')
# plt.plot(data[:,0],data[:,0]*a+logA,'r')
Explanation: 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) = ax + \log(A) $$
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)
End of explanation
x = np.linspace(1,45)
plt.title("Raw Data")
plt.xlabel("Distance")
plt.semilogy(data[:,0],data[:,1],'bo',label='data')
plt.semilogy(x,np.exp(logA)*np.exp(a*x),'r-',label='fit')
plt.legend()
Explanation: Let's see whether this curve fits the data.
End of explanation
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
gdata = []
for line in gauss_data.splitlines():
    words = line.split(',')
    gdata.append(list(map(float,words)))
gdata = np.array(gdata)
plt.plot(gdata[:,0],gdata[:,1],'bo')
Explanation: If we have more complicated functions, we may not be able to get away with fitting to a simple polynomial. Consider the following data:
End of explanation
def gauss(x,A,a): return A*np.exp(a*x**2)
Explanation: 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.
End of explanation
from scipy.optimize import curve_fit
params,conv = curve_fit(gauss,gdata[:,0],gdata[:,1])
x = np.linspace(-1,1)
plt.plot(gdata[:,0],gdata[:,1],'bo')
A,a = params
plt.plot(x,gauss(x,A,a),'r-')
Explanation: Now fit to it using curve_fit:
End of explanation
# from random import random
rands = []
for i in range(100):
    rands.append(np.random.random())
Explanation: 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.
2.9 Monte Carlo, random numbers, and computing $\pi$
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 from the numpy library gives pseudorandom numbers uniformly distributed between 0 and 1:
End of explanation
rands = np.random.rand(100)
plt.plot(rands,'o')
Explanation: Or, more elegantly:
End of explanation
mu, sigma = 0, 0.1 # mean and standard deviation
s=np.random.normal(mu, sigma,1000)
Explanation: np.random.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:
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:
End of explanation
count, bins, ignored = plt.hist(s, 30, normed=True)
plt.plot(bins, 1/(sigma * np.sqrt(2 * np.pi)) *
             np.exp( - (bins - mu)**2 / (2 * sigma**2) ),
             linewidth=2, color='r')
Explanation: We can check the distribution by using the histogram feature, as shown on the help page for numpy.random.normal:
End of explanation
npts = 5000
xs = 2*np.random.rand(npts)-1
ys = 2*np.random.rand(npts)-1
r = xs**2+ys**2
ninside = (r<1).sum()
plt.figure(figsize=(6,6)) # make the figure square
plt.title("Approximation to pi = %f" % (4*ninside/float(npts)))
plt.plot(xs[r<1],ys[r<1],'b.')
plt.plot(xs[r>1],ys[r>1],'r.')
plt.figure(figsize=(8,6)) # change the figsize back to standard size for the rest of the notebook
Explanation: Here's an interesting use of random numbers: compute $\pi$ by taking random numbers as x and y coordinates, and counting how many of them were in the unit circle. For example:
End of explanation
n = 100
total = 0
for k in range(n):
    total += pow(-1,k)/(2*k+1.0)
print(4*total)
Explanation: The idea behind the program is that the ratio of the area of the unit circle to the square that inscribes it is $\pi/4$, so by counting the fraction of the random points in the square that are inside the circle, we get increasingly good estimates to $\pi$. 
The above code uses some higher level Numpy tricks to compute the radius of each point in a single line, to count how many radii are below one in a single line, and to filter the x,y points based on their radii. To be honest, I rarely write code like this: I find some of these Numpy tricks a little too cute to remember them, and I'm more likely to use a list comprehension (see below) to filter the points I want, since I can remember that.
As methods of computing $\pi$ go, this is among the worst. A much better method is to use Leibniz's expansion of arctan(1):
$$\frac{\pi}{4} = \sum_k \frac{(-1)^k}{2*k+1}$$
End of explanation
def f(x): return np.exp(-x)
x = np.linspace(0,10)
plt.plot(x,np.exp(-x))
Explanation: If you're interested in another great method, check out Ramanujan's method. This converges so fast you really need arbitrary precision math to display enough decimal places. You can do this with the Python decimal module, if you're interested.
2.10 Numerical Integration
Integration can be hard, and sometimes it's easier to work out a definite integral using an approximation. For example, suppose we wanted to figure out the integral:
$$\int_0^\infty\exp(-x)dx$$
(It turns out that this is equal to 1, as you can work out easily with a pencil :) )
End of explanation
from scipy.integrate import quad
quad(f,0,np.inf)
Explanation: Scipy has a numerical integration routine quad (since sometimes numerical integration is called quadrature), that we can use for this:
End of explanation
from scipy.fftpack import fft,fftfreq
npts = 4000
nplot = int(npts/10)
t = np.linspace(0,120,npts)
def sig(t): return 50*np.sin(2*np.pi*2.0*t) +  20*np.sin(2*np.pi*5.0*t) + 10*np.sin(2*np.pi*8.0*t) + 2*np.random.rand(npts)
Vsignal = sig(t)
FFT = abs(fft(Vsignal))
freqs = fftfreq(npts, t[1]-t[0])
FFT_plot = FFT[0:int(len(freqs)/2)]
freqs_plot = freqs[0:int(len(freqs)/2)]
plt.subplot(211)
plt.plot(t[:nplot], Vsignal[:nplot])
plt.xlabel ('time (s)')
plt.ylabel ('voltage\nmeasured (V)')
plt.subplot(212)
plt.semilogy(freqs_plot,FFT_plot**2,'-')
plt.xlabel ('frequency (Hz)')
plt.ylabel ('power\nspectrum (a.u.)')
plt.ylim([1e-1,np.max(FFT_plot**2)])
plt.tight_layout()
Explanation: The first number in the tuple is the result, the second number is an estimate of the absolute error in the result.
There are also 2d and 3d numerical integrators in Scipy. See the docs for more information.
2.11 Fast Fourier Transform and Signal Processing
Very often we want to use FFT techniques to help obtain the signal from noisy data. Scipy has several different options for this.
End of explanation
csv = \
1, -6095.12544083, 0.03686, 1391.5
2, -6095.25762870, 0.00732, 10468.0
3, -6095.26325979, 0.00233, 11963.5
4, -6095.26428124, 0.00109, 13331.9
5, -6095.26463203, 0.00057, 14710.8
6, -6095.26477615, 0.00043, 20211.1
7, -6095.26482624, 0.00015, 21726.1
8, -6095.26483584, 0.00021, 24890.5
9, -6095.26484405, 0.00005, 26448.7
10, -6095.26484599, 0.00003, 27258.1
11, -6095.26484676, 0.00003, 28155.3
12, -6095.26484693, 0.00002, 28981.7
13, -6095.26484693, 0.00002, 28981.7
csv
Explanation: There are additional signal processing routines in Scipy (e.g. splines, filtering) that you can read about here.
3. Intermediate Python
3.1 Parsing data output
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 CSV (comma separated values) format, a format that originally came from Microsoft Excel, and is increasingly used as a data interchange format in big data applications. How would we parse that?
End of explanation
lines = csv.splitlines()
lines
Explanation: This is a giant string. If we use splitlines(), we see that a list is created where line gets separated into a string:
End of explanation
lines[4].split(",")
Explanation: Splitting is a big concept in text processing. We used splitlines() here, and next we'll use the more general .split(",") function below to split each line into comma-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.
To break apart each line, we will use .split(","). Let's see what it does to one of the lines:
End of explanation
help("".split)
Explanation: What does split() do?
End of explanation
for line in lines:
    # do something with each line
    words = line.split(",")
Explanation: Since the data is now in a list of lines, we can iterate over it, splitting up data where we see a comma:
End of explanation
data = []
for line in csv.splitlines()[2:]:
    words = line.split(',')
    data.append(list(map(float,words)))
data = np.array(data)
data
Explanation: We need to add these results at each step to a list:
End of explanation
plt.plot(data[:,0],data[:,1],'-o')
plt.xlabel('step')
plt.ylabel('Energy (hartrees)')
plt.title('Convergence of NWChem geometry optimization for Si cluster\n')
Explanation: 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. Similarly, [:5] instead would take the first five lines.
We pass the comma string ",'"into the split function, so that it breaks to a new word every time it sees a comma. Next, to simplify things a bit, we're using the map() command to repeatedly apply a single function (float()) to a list, and to return the output as a list. Finally, we turn the list of lists into a numpy arrray structure.
End of explanation
energies = data[:,1]
minE = np.min(energies)
energies_eV = 27.211*(energies-minE)
plt.plot(data[:,0],energies_eV,'-o')
plt.xlabel('step')
plt.ylabel('Energy (eV)')
plt.title('Convergence of NWChem geometry optimization for Si cluster')
Explanation: Hartrees (what most quantum chemistry programs use by default) are really stupid units. We really want this in kcal/mol or eV or something we use. So let's quickly replot this in terms of eV above the minimum energy, which will give us a much more useful plot:
End of explanation
filename= 'DS0004.csv'
data = np.genfromtxt(filename,delimiter=',',skip_header=17 ) 
x_values = data[:,0]  
y_values = data[:,1]
plt.plot(x_values, y_values)
Explanation: 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.
3.2 Reading in data files
Let's take a look at a perhaps easier approach to a common problem -- you have a data file with some header info and comma-delimited values and you want the data so you can start doing stuff with it. Let's use numpy's genfromtxt()
End of explanation
print("I have 3 errands to run")
Explanation: That was easy! Why didn't we only learn that? Because not every data set is "nice" like that. Better to have some tools for when things aren't working how you'd like them to be.  That being said, much data coming from scientific equipment and computational tools can be cast into a format that can be read in through genfromtxt(). For larger data sets, the library pandas might be helpful.
3.3 More Sophisticated String Formatting and Processing
Strings are a big deal in most modern languages, and hopefully the previous sections helped underscore how versatile Python's string processing techniques are. We will continue this topic in this section.
We can print out lines in Python using the print command.
End of explanation
"I have 3 errands to run"
Explanation: In IPython we don't even need the print command, since it will display the last expression not assigned to a variable.
End of explanation
a,b,c = 1,2,3
print("The variables are ",1,2,3)
Explanation: print even converts some arguments to strings for us:
End of explanation
print("Pi as a decimal = %d" % np.pi)
print("Pi as a float = %f" % np.pi)
print("Pi with 4 decimal places = %.4f" % np.pi)
print("Pi with overall fixed length of 10 spaces, with 6 decimal places = %10.6f" % np.pi)
print("Pi as in exponential format = %e" % np.pi)
Explanation: As versatile as this is, you typically need more freedom over the data you print out. For example, what if we want to print a bunch of data to exactly 4 decimal places? We can do this using formatted strings.
Formatted strings share a syntax with the C printf statement. We make a string that has some funny format characters in it, and then pass a bunch of variables into the string that fill out those characters in different ways.
For example,
End of explanation
print("The variables specified earlier are %d, %d, and %d" % (a,b,c))
Explanation: We use a percent sign in two different ways here. First, the format character itself starts with a percent sign. %d or %i are for integers, %f is for floats, %e is for numbers in exponential formats. All of the numbers can take number immediately after the percent that specifies the total spaces used to print the number. Formats with a decimal can take an additional number after a dot . to specify the number of decimal places to print.
The other use of the percent sign is after the string, to pipe a set of variables in. You can pass in multiple variables (if your formatting string supports it) by putting a tuple after the percent. Thus,
End of explanation
form_letter = \
          %s
Dear %s,
We regret to inform you that your product did not
ship today due to %s.
We hope to remedy this as soon as possible.
          From,
          Your Supplier
print(form_letter % ("July 1, 2016","Valued Customer Bob","alien attack"))
Explanation: This is a simple formatting structure that will satisfy most of your string formatting needs. More information on different format symbols is available in the string formatting part of the standard docs.
It's worth noting that more complicated string formatting methods are in development, but I prefer this system due to its simplicity and its similarity to C formatting strings.
Recall we discussed multiline strings. We can put format characters in these as well, and fill them with the percent sign as before.
End of explanation
form_letter = \
          %(date)s
Dear %(customer)s,
We regret to inform you that your product did not
ship today due to %(lame_excuse)s.
We hope to remedy this as soon as possible.
          From,
          Your Supplier
print(form_letter % {"date" : "July 1, 2016","customer":"Valued Customer Bob","lame_excuse":"alien attack"})
Explanation: The problem with a long block of text like this is that it's often hard to keep track of what all of the variables are supposed to stand for. There's an alternate format where you can pass a dictionary into the formatted string, and give a little bit more information to the formatted string itself. This method looks like:
End of explanation
nwchem_format = 
start %(jobname)s
title "%(thetitle)s"
charge %(charge)d
geometry units angstroms print xyz autosym
%(geometry)s
end
basis
  * library 6-31G**
end
dft
  xc %(dft_functional)s
  mult %(multiplicity)d
end
task dft %(jobtype)s
Explanation: By providing a little bit more information, you're less likely to make mistakes, like referring to your customer as "alien attack".
As a scientist, you're less likely to be sending bulk mailings to a bunch of customers. But these are great methods for generating and submitting lots of similar runs, say scanning a bunch of different structures to find the optimal configuration for something.
For example, you can use the following template for NWChem input files:
End of explanation
oxygen_xy_coords = [(0,0),(0,0.1),(0.1,0),(0.1,0.1)]
charge = 0
multiplicity = 1
dft_functional = "b3lyp"
jobtype = "optimize"
geometry_template = \
  O    %f     %f      0.0
  H    0.0    1.0     0.0
  H    1.0    0.0     0.0
for i,xy in enumerate(oxygen_xy_coords):
    thetitle = "Water run #%d" % i
    jobname = "h2o-%d" % i
    geometry = geometry_template % xy
    print("---------")
    print(nwchem_format % dict(thetitle=thetitle,charge=charge,jobname=jobname,jobtype=jobtype,
                               geometry=geometry,dft_functional=dft_functional,multiplicity=multiplicity))
Explanation: If you want to submit a sequence of runs to a computer somewhere, it's pretty easy to put together a little script, maybe even with some more string formatting in it:
End of explanation
def my_enumerate(seq):
    l = []
    for i in range(len(seq)):
        l.append((i,seq[i]))
    return l
my_enumerate(oxygen_xy_coords)
Explanation: This is a very bad geometry for a water molecule, and it would be silly to run so many geometry optimizations of structures that are guaranteed to converge to the same single geometry, but you get the idea of how you can run vast numbers of simulations with a technique like this.
We used the enumerate function to loop over both the indices and the items of a sequence, which is valuable when you want a clean way of getting both. enumerate is roughly equivalent to:
End of explanation
np.linspace(0,1)
Explanation: Although enumerate uses generators (see below) so that it doesn't have to create a big list, which makes it faster for really long sequenes.
3.4 Optional arguments of a function
You will recall that the linspace function can take either two arguments (for the starting and ending points):
End of explanation
np.linspace(0,1,5)
Explanation: or it can take three arguments, for the starting point, the ending point, and the number of points:
End of explanation
np.linspace(0,1,5,endpoint=False)
Explanation: You can also pass in keywords to exclude the endpoint:
End of explanation
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)
Explanation: 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:
End of explanation
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
Explanation: We can add an optional argument by specifying a default value in the argument list:
End of explanation
my_linspace(0,1)
Explanation: This gives exactly the same result if we don't specify anything:
End of explanation
my_linspace(0,1,5)
Explanation: But also let's us override the default value with a third argument:
End of explanation
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)
Explanation: We can add arbitrary keyword arguments to the function definition by putting a keyword argument **kwargs handle in:
End of explanation
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
Explanation: 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?
End of explanation
my_range()
Explanation: 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:
End of explanation
evens1 = [2*i for i in range(10)]
print(evens1)
Explanation: 3.5 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:
End of explanation
odds = [i for i in range(20) if i%2==1]
odds
Explanation: You can also put some boolean testing into the construct:
End of explanation
def evens_below(n):
    for i in range(n):
        if i%2 == 0:
            yield i
    return
for i in evens_below(9):
    print(i)
Explanation: 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:
End of explanation
list(evens_below(9))
Explanation: We can always turn an iterator into a list using the list command:
End of explanation
evens_gen = (i for i in range(9) if i%2==0)
for i in evens_gen:
    print(i)
Explanation: There's a special syntax called a generator expression that looks a lot like a list comprehension:
End of explanation
def gauss(x,A,a,x0):
    return A*np.exp(-a*(x-x0)**2)
Explanation: 3.6 Factory Functions
A factory function is a function that returns a function. They have the fancy name lexical closure, which makes you sound really intelligent in front of your CS friends. But, despite the arcane names, factory functions can play a very practical role.
Suppose you want the Gaussian function centered at 0.5, with height 99 and width 1.0. You could write a general function.
End of explanation
def gauss_maker(A,a,x0):
    def f(x):
        return A*np.exp(-a*(x-x0)**2)
    return f
x = np.linspace(0,1)
g = gauss_maker(99.0,20,0.5)
plt.plot(x,g(x))
Explanation: But what if you need a function with only one argument, like f(x) rather than f(x,y,z,...)? You can do this with Factory Functions:
End of explanation
# Data in a json format:
json_data = \
{
    "a": [1,2,3],
    "b": [4,5,6],
    "greeting" : "Hello"
}
import json
loaded_json=json.loads(json_data)
loaded_json
Explanation: Everything in Python is an object, including functions. This means that functions can be returned by other functions. (They can also be passed into other functions, which is also useful, but a topic for another discussion.) In the gauss_maker example, the g function that is output "remembers" the A, a, x0 values it was constructed with, since they're all stored in the local memory space (this is what the lexical closure really refers to) of that function.
Factories are one of the more important of the Software Design Patterns, which are a set of guidelines to follow to make high-quality, portable, readable, stable software. It's beyond the scope of the current work to go more into either factories or design patterns, but I thought I would mention them for people interested in software design.
3.7 Serialization: Save it for later
Serialization refers to the process of outputting data (and occasionally functions) to a database or a regular file, for the purpose of using it later on. In the very early days of programming languages, this was normally done in regular text files. Python is excellent at text processing, and you probably already know enough to get started with this.
When accessing large amounts of data became important, people developed database software based around the Structured Query Language (SQL) standard. I'm not going to cover SQL here, but, if you're interested, I recommend using the sqlite3 module in the Python standard library.
As data interchange became important, the eXtensible Markup Language (XML) has emerged. XML makes data formats that are easy to write parsers for, greatly simplifying the ambiguity that sometimes arises in the process. Again, I'm not going to cover XML here, but if you're interested in learning more, look into Element Trees, now part of the Python standard library.
Python has a very general serialization format called pickle that can turn any Python object, even a function or a class, into a representation that can be written to a file and read in later. But, again, I'm not going to talk about this, since I rarely use it myself. Again, the standard library documentation for pickle is the place to go.
What I am going to talk about is a relatively recent format call JavaScript Object Notation (JSON) that has become very popular over the past few years. There's a module in the standard library for encoding and decoding JSON formats. The reason I like JSON so much is that it looks almost like Python, so that, unlike the other options, you can look at your data and edit it, use it in another program, etc.
Here's a little example:
End of explanation
json.dumps({"a":[1,2,3],"b":[9,10,11],"greeting":"Hola"})
Explanation: Your data sits in something that looks like a Python dictionary, and in a single line of code, you can load it into a Python dictionary for use later.
In the same way, you can, with a single line of code, put a bunch of variables into a dictionary, and then output to a file using json:
End of explanation
from operator import add, mul
add(1,2)
mul(3,4)
Explanation: 3.8 Functional programming
Functional programming is a very broad subject. The idea is to have a series of functions, each of which generates a new data structure from an input, without changing the input structure at all. By not modifying the input structure (something that is called not having side effects), many guarantees can be made about how independent the processes are, which can help parallelization and guarantees of program accuracy. There is a Python Functional Programming HOWTO in the standard docs that goes into more details on functional programming. I just wanted to touch on a few of the most important ideas here.
There is an operator module that has function versions of most of the Python operators. For example:
End of explanation
def doubler(x): return 2*x
doubler(17)
Explanation: These are useful building blocks for functional programming.
The lambda operator allows us to build anonymous functions, which are simply functions that aren't defined by a normal def statement with a name. For example, a function that doubles the input is:
End of explanation
lambda x: 2*x
Explanation: We could also write this as:
End of explanation
another_doubler = lambda x: 2*x
another_doubler(19)
Explanation: And assign it to a function separately:
End of explanation
list(map(float,'1 2 3 4 5'.split()))
Explanation: lambda is particularly convenient (as we'll see below) in passing simple functions as arguments to other functions.
map is a way to repeatedly apply a function to a list:
End of explanation
sum([1,2,3,4,5])
Explanation: reduce is a way to repeatedly apply a function to the first two items of the list. There already is a sum function in Python that is a reduction:
End of explanation
from functools import reduce
def prod(l): return reduce(mul,l)
prod([1,2,3,4,5])
Explanation: We can use reduce to define an analogous prod function:
End of explanation
mystring = "Hi there"
Explanation: 3.9 Object Oriented Programming
We've seen a lot of examples of objects in Python. We create a string object with quote marks:
End of explanation
mystring.split()
mystring.startswith('Hi')
len(mystring)
Explanation: and we have a bunch of methods we can use on the object:
End of explanation
class Schrod1d:
    \
    Schrod1d: Solver for the one-dimensional Schrodinger equation.
    
    def __init__(self,V,start=0,end=1,npts=50,**kwargs):
        m = kwargs.get('m',1.0)
        self.x = np.linspace(start,end,npts)
        self.Vx = V(self.x)
        self.H = (-0.5/m)*self.laplacian() + np.diag(self.Vx)
        return
    
    def plot(self,*args,**kwargs):
        titlestring = kwargs.get('titlestring',"Eigenfunctions of the 1d Potential")
        xstring = kwargs.get('xstring',"Displacement (bohr)")
        ystring = kwargs.get('ystring',"Energy (hartree)")
        if not args:
            args = [3]
        x = self.x
        E,U = np.linalg.eigh(self.H)
        h = x[1]-x[0]
        # Plot the Potential
        plt.plot(x,self.Vx,color='k')
        for i in range(*args):
            # For each of the first few solutions, plot the energy level:
            plt.axhline(y=E[i],color='k',ls=":")
            # as well as the eigenfunction, displaced by the energy level so they don't
            # all pile up on each other:
            plt.plot(x,U[:,i]/np.sqrt(h)+E[i])
        plt.title(titlestring)
        plt.xlabel(xstring)
        plt.ylabel(ystring) 
        return
        
    def laplacian(self):
        x = self.x
        h = x[1]-x[0] # assume uniformly spaced points
        n = len(x)
        M = -2*np.identity(n,'d')
        for i in range(1,n):
            M[i,i-1] = M[i-1,i] = 1
        return M/h**2
Explanation: Object oriented programming simply gives you the tools to define objects and methods for yourself. It's useful anytime you want to keep some data (like the characters in the string) tightly coupled to the functions that act on the data (length, split, startswith, etc.).
As an example, we're going to bundle the functions we did to make the 1d harmonic oscillator eigenfunctions with arbitrary potentials, so we can pass in a function defining that potential, some additional specifications, and get out something that can plot the orbitals, as well as do other things with them, if desired.
End of explanation
square_well = Schrod1d(lambda x: 0*x,m=10)
square_well.plot(4,titlestring="Square Well Potential")
Explanation: The init() function specifies what operations go on when the object is created. The self argument is the object itself, and we don't pass it in. The only required argument is the function that defines the QM potential. We can also specify additional arguments that define the numerical grid that we're going to use for the calculation.
For example, to do an infinite square well potential, we have a function that is 0 everywhere. We don't have to specify the barriers, since we'll only define the potential in the well, which means that it can't be defined anywhere else.
End of explanation
ho = Schrod1d(lambda x: x**2,start=-3,end=3)
ho.plot(6,titlestring="Harmonic Oscillator")
Explanation: We can similarly redefine the Harmonic Oscillator potential.
End of explanation
def finite_well(x,V_left=1,V_well=0,V_right=1,d_left=10,d_well=10,d_right=10):
    V = np.zeros(x.size,'d')
    for i in range(x.size):
        if x[i] < d_left: 
            V[i] = V_left
        elif x[i] > (d_left+d_well):
            V[i] = V_right
        else:
            V[i] = V_well
    return V
        
fw = Schrod1d(finite_well,start=0,end=30,npts=100)
fw.plot()
Explanation: Let's define a finite well potential:
End of explanation
def triangular(x,F=30): return F*x
tw = Schrod1d(triangular,m=10)
tw.plot()
Explanation: A triangular well:
End of explanation
def tri_finite(x): return finite_well(x)+triangular(x,F=0.025)
tfw = Schrod1d(tri_finite,start=0,end=30,npts=100)
tfw.plot()
Explanation: Or we can combine the two, making something like a semiconductor quantum well with a top gate:
End of explanation
%timeit factorial(20)
Explanation: There's a lot of philosophy behind object oriented programming. Since I'm trying to focus on just the basics here, I won't go into them, but the internet is full of lots of resources on OO programming and theory. The best of this is contained in the Design Patterns book, which I highly recommend.
4. Speeding Python: Timeit, Profiling, Cython, SWIG, and PyPy
The first rule of speeding up your code is not to do it at all. As Donald Knuth said:
"We should forget about small efficiencies, say about 97% of the time: premature optimization is the root of all evil."
The second rule of speeding up your code is to only do it if you really think you need to do it. Python has two tools to help with this process: a timing program called timeit, and a very good code profiler. We will discuss both of these tools in this section, as well as techniques to use to speed up your code once you know it's too slow.
4.1 Timeit
timeit helps determine which of two similar routines is faster. Recall that some time ago we wrote a factorial routine, but also pointed out that Python had its own routine built into the math module. Is there any difference in the speed of the two? timeit helps us determine this. For example, timeit tells how long each method takes:
End of explanation
%timeit fact(20)
Explanation: The little % sign that we have in front of the timeit call is an example of an IPython magic function, which we don't have time to go into here, but it's just some little extra mojo that IPython adds to the functions to make it run better in the IPython environment. You can read more about it in the IPython tutorial.
In any case, the timeit function runs 3 loops, and tells us that it took on the average of 583 ns to compute 20!. In contrast:
End of explanation
def evens(n):
    "Return a list of even numbers below n"
    l = []
    for x in range(n):
        if x % 2 == 0:
            l.append(x)
    return l
Explanation: the factorial function we wrote is about a factor of 10 slower. This is because the built-in factorial function is written in C code and called from Python, and the version we wrote is written in plain old Python. A Python program has a lot of stuff in it that make it nice to interact with, but all that friendliness slows down the code. In contrast, the C code is less friendly but more efficient. If you want speed with as little effort as possible, write your code in an easy to program language like Python, but dump the slow parts into a faster language like C, and call it from Python. We'll go through some tricks to do this in this section.
4.2 Profiling
Profiling complements what timeit does by splitting the overall timing into the time spent in each function. It can give us a better understanding of what our program is really spending its time on.
Suppose we want to create a list of even numbers. Our first effort yields this:
End of explanation
import cProfile
cProfile.run('evens(100000)')
Explanation: Is this code fast enough? We find out by running the Python profiler on a longer run:
End of explanation
def evens2(n):
    "Return a list of even numbers below n"
    return [x for x in range(n) if x % 2 == 0]
import cProfile
cProfile.run('evens2(100000)')
Explanation: This looks okay, 0.05 seconds isn't a huge amount of time, but looking at the profiling shows that the append function is taking almost 20% of the time. Can we do better? Let's try a list comprehension.
End of explanation
def evens3(n):
    "Return a list of even numbers below n"
    return [x for x in range(n) if x % 2 == 0]
import cProfile
cProfile.run('evens3(100000)')
Explanation: By removing a small part of the code using a list comprehension, we've doubled the overall speed of the code! 
It seems like range is taking a long time, still. Can we get rid of it? We can, using the xrange generator:
End of explanation
def primes(n):
    \
    From python cookbook, returns a list of prime numbers from 2 to < n
    >>> primes(2)
    [2]
    >>> primes(10)
    [2, 3, 5, 7]
    
    if n==2: return [2]
    elif n<2: return []
    s=list(range(3,n+2,2))
    mroot = n ** 0.5
    half=(n+1)/2-1
    i=0
    m=3
    while m <= mroot:
        if s[i]:
            j=int((m*m-3)/2)
            s[j]=0
            while j<half:
                s[j]=0
                j+=m
        i=i+1
        m=2*i+3
    return [2]+[x for x in s if x]
number_to_try = 1000000
list_of_primes = primes(number_to_try)
print(list_of_primes[10001])
Explanation: This is where profiling can be useful. Our code now runs 3x faster by making trivial changes. We wouldn't have thought to look in these places had we not had access to easy profiling. Imagine what you would find in more complicated programs.
4.3 Other Ways to Speed Python
When we compared the fact and factorial functions, above, we noted that C routines are often faster because they're more streamlined. Once we've determined that one routine is a bottleneck for the performance of a program, we can replace it with a faster version by writing it in C. This is called extending Python, and there's a good section in the standard documents. This can be a tedious process if you have many different routines to convert. Fortunately, there are several other options.
Swig (the simplified wrapper and interface generator) is a method to generate binding not only for Python but also for Matlab, Perl, Ruby, and other scripting languages. Swig can scan the header files of a C project and generate Python binding for it. Using Swig is substantially easier than writing the routines in C.
Cython is a C-extension language. You can start by compiling a Python routine into a shared object libraries that can be imported into faster versions of the routines. You can then add additional static typing and make other restrictions to further speed the code. Cython is generally easier than using Swig.
PyPy is the easiest way of obtaining fast code. PyPy compiles Python to a subset of the Python language called RPython that can be efficiently compiled and optimized. Over a wide range of tests, PyPy is roughly 6 times faster than the standard Python Distribution.
4.4 Fun: Finding Primes
Project Euler is a site where programming puzzles are posed that might have interested Euler. Problem 7 asks the question:
By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.
What is the 10,001st prime number?
To solve this we need a very long list of prime numbers. First we'll make a function that uses the Sieve of Erastothenes to generate all the primes less than n.
End of explanation
cProfile.run('primes(1000000)')
Explanation: You might think that Python is a bad choice for something like this, but, in terms of time, it really doesn't take long:
End of explanation
<END_TASK>