Minimal Python for Scientific Computing

©M.D. Johnson, September 26, 2002; Revised August 19, 2003

This document is intended to help non-programmers get started using Python for scientific computing. I don’t assume any programming experience – but I do assume that Python and either SciPy or NumPy are already installed on your computer.

This tutorial is designed for participative learning: type the examples into a Python session, think about the output, and try examples as they occur to you. This is not even close to a complete explanation of the language. References to more comprehensive introductions and documentation are shown at the end.

Starting Python

Start an interactive Python session using one of:

Idle	- Windows or Unix
Pythonwin	- Windows
Python IDE	- MacPython under MacOS X

Each opens an interactive Python window which shows the prompt >>>.

Using Python like a calculator

To get a feel for how Python works, try the following. The hash symbol # marks the beginning of a comment in Python. I’ve used that here to explain some of the commands. Don’t type in the comments.

4./3.	# The default floating point type is double precision.
4/3	# Integer division truncates back to integers
4 / 3.	# Converts type.
5*3
5.**3
5*3-2
5*(3-2)
10.2 % 3	# Remainder.
x=3.	# Use of variables is very natural.
y=x**3	# After assignment (=), Python doesn’t show the value.
x,y
print x,y
print X	# Case sensitive.
Vinitial = 4.	# Use meaningful names (begin with a letter).
x == 1	# Logic tests return 1 (true) or 0 (false).
x == 3.
x != y	# Not equal.
x <= y
z=3+2J	# Complex type.
print z
w=2.-2j	# Either J or j works.
z*w
_*2	# The last computed value is in variable “_”.
z.real	# The first hint that Python is object-oriented: real and imag
z.imag	# are “attributes” of complex numbers.
z.conjugate()	# Objects can have associated functions, called “methods:”
	# conjugate() is a method of complex numbers.

Idle and Pythonwin have limited command-line editing abilities. Try moving to a higher line using the up arrow, hitting return, editing, and hitting return again. With IDE you have to cut and paste.

Data types

You’ve seen the simplest numerical types:

a=3	# Integer.
b=4.	# Floating point (double precision).
c=4+2j	# Complex (double precision).
type(a)	# A function that shows the type. Notice that function
	# arguments are enclosed within parentheses.

A string consists of characters within single or double quotation marks:

a = “An”
b = ‘ex’	# Single or double quotation marks.
c = “parrot”
print a,b,c	# Adds a space between each pair.
print a+b+c	# Concatenates strings.
a[0]	# You can index each character (beginning with zero).
a[1]	# Notice the square brackets.
len(a)	# Length.

Python has three other built-in container data types besides strings: lists, tuples, and dictionaries. Lists are written in square brackets, with the elements separated by commas:

x = [ 6,7,8,9,10 ]
type(x)
x[0]	# Lists can be indexed. Notice the square brackets.
x[1]	# The first index is 0, not 1.
len(x)	# Length.
x[-1]	# A negative index begins from the right.
range(5)	# A built-in function which produces a simple list.
range(3,9)

Objects which can be indexed (such as lists and strings) can also be accessed by slices, which give a range of indices:

x[1:3]	# The first index is included but the last is not!
x[:3]	# Dropping one of the indices means “everything in that direction.”
x[3:]
x[:]	# Dropping both indices means “everything.”

The way Python works with ranges takes some getting used to: x[first:last] means the elements x[first], x[first+1],…,x[last-1]. That is, the last element is not included. This means something like “first <= index < last.” This same behavior showed up above when you used range.

Lists are mutable. That means that elements of the list can be changed while leaving other elements alone.

x[0] = 0

print x

Lists can be made up of anything, not just numbers.

y = [ “Nobody” , “expects” ]
z = [ “the” , “Spanish” , “Inquisition” ]
y + z	# What does addition mean for lists?

The second built-in container data type is the tuple. These are similar to lists, but without the square brackets:

y = 1,2,3
type(y)
print y
y[0]	# Can be indexed …
y[0:2]	# … and sliced.
z = (1,2,3)	# The parentheses are optional.
print z

Unlike lists, tuples are immutable. That means that individual elements cannot be modified.

y[0] = 3

# Gives an error.

Tuples can be quickly packed and unpacked into other variables:

a,b,c = y	# Unpacks y (which can be a tuple or a list) into a,b,c.
print a,b,c
a = 2*a
print a
print y
y = a,b,c	# Packs a,b,c into y (a tuple).
print y

Dictionaries are like lists without ordering. They are indexed by “key” rather than position. They are denoted with braces:

review = { “Inquisition” : “Funny” , “Parrot” : “Hysterical”}
print review[0]	# Not indexed by number …
print review[‘Inquisition’]	# … but by key.
print review[“Parrot”]

Dictionaries are a very convenient way to organize large collections of data, process them, and keep track of everything along the way. You won’t usually use them directly in numerical work, but sometimes they are very efficient. (They are used in SciPy to implement sparse matrix operations). Since the entries can be any object (including functions), dictionaries can be remarkably powerful.

Scientific and Numerical Python

The Python core doesn’t have the basic data types and functions most needed for scientific programming. Many are available in an add-on module called Numerical Python (or NumPy). This in turn is part of a more comprehensive module called SciPy which includes NumPy plus much more. As of this writing SciPy is only available on Windows and Unix machines. I assume you have already installed one of these.

To get access to these features they need to be imported. Once per session type:

from scipy import *	# If SciPy is installed; or
from Numeric import *	# If NumPy (but not SciPy) is installed.

You will also need to do this in any module (see below) that needs these features.

NumPy defines a new container data type which is important for scientific programming: the multiarray, or, for short, the array. Arrays are similar to lists but have properties closer to matrices than to strings. Define an array using the function array() with (usually) one argument:

a = array( (1,2,3) )
type(a)
a
print a
b = array( [1,2,3] )	# Brackets or parentheses are both OK.
print b
c = array( 1,2,3 )	# But you need one or the other.
a / 2	# Oops … integer arithmetic
a = array( (1.,2,3) )	# Chooses the “bigger” type automatically.
a / 2
a = array((1,2,3),Float)	# Or set the type with an optional second argument.
print a
b = array( [4,5,6] )
a + b	# Add like matrices.
a=array( [[1,2],[3,4]] )	# 2 by 2 array.
print a

Many of the imported numerical functions are “universal functions,” or ufuncs, which means that they act element-by-element on arrays.

a = arange( 10 )	# “array range”: more efficient than range().
print a
a**2	# Squares each element.
a * (pi/9)	# pi is part of NumPy.
sin(a)	# Takes the sine of each element of a.
sin(2+3j)	# Works with any data type.

Elemental operations are very useful, but make sure that you don’t erroneously expect ordinary matrix multiplication.

a = array( (1,2) )
print a*2	# Does what you’d expect.
b = array( (4,5) )
print a*b	# Elemental multiplication
print dot(a,b)	# What you’d expect.
a = array( ((1,2),(3,4)) )
print a
print b
print matrixmultiply(a,b)	# What you’d expect
print a*b

This last operation needs explanation. When two arrays don’t have compatible shapes, Python tries to broadcast the smaller to the larger. Sometimes that’s easy to understand:

a = array( ((1,2),(3,4)) )

print a+2

Here the scalar 2 is added to each element of the matrix. Now try

b = array( (4,5) )

print a+b

You might be able to see what happened. The first element of b was added to the first column of a, and the second to the second. The same thing happens with a*b.

Programming

You won’t use the interactive window very much. Instead you will write programs in files and then run the programs. A collection of Python code saved in a file is called a module. The Python interfaces have special editors to make this easier. To begin a new module, click File and then New (or New Window). (Using Pythonwin, you will then need to choose “Python Script.”) In each case this starts an editor. Type in some very simple Python code, such as:

t=2

y=t**2

print y

You will of course need to save your modules. Be careful where you put them. Don’t use the default location (somewhere down inside the Python directory structure). Use a directory of your own. And use reasonable names (not “Script1”). The details of how to save and run depend on which interface you’re using:

· Idle: Save from the File menu. You must explicitly type a suffix .py when you name the module (e.g., “first.py”). Then run by using “Edit->Run Script.”

· Pythonwin: You can save either from the File menu or by clicking on the disk icon in the tool bar. Do not explicitly add the suffix .py (e.g., just give it the name “first”). Once it’s saved you can run it using “File->Run,” or by clicking the little running person icon in the tool bar. These bring up a dialogue box. Click OK and the code runs. (Holding down the shift button while clicking skips the dialogue box.)

· IDE: Save from the File menu. Certain things work better if you add the suffix .py to the name. Run your code by clicking “Run all” on the editor window, or “Python->Run window” from the menu bar.

Change your code and run it again. (For instance, make t=4.) With Idle you need to explicitly save before running again. With Pythonwin the act of running automatically saves the modified version. With IDLE saving is unnecessary (until you’re finished!).

Python will complain if you try to run a module containing a syntax error. For example, change the middle line above to y=t***2 and try to run.

The editors can help you find syntax errors. Using Idle, select “File->Check Module.” This will highlight the first error it finds. Using Pythonwin, select “File->Check” or click the check mark in the tool bar. The cursor moves to the point where a syntax error is found. With IDE, click the little triangle just below the upper right corner in the editor’s window, and then select “Modularize.” With any of these, once you fix the error, checking again will find the next error (if any). [Below you will see that editors also help catch syntax errors by adding automatic indentation and (maybe) colors. For example, if you notice that the indentation is wrong in some block, you probably omitted a colon.]

To get all the syntax bugs out of your program: make a change and then check or run. Do it again; and again; and again … until it’s perfect.

Control statements

Python has for and while loops, and if statements. The syntax is peculiar: Python uses indentation to determine how far a block extends instead of “end if” or “end for” statements. Make and run a module with the following:

x=5
if x<4 or x>27:		Notice the colon!
	print “Bring”
	print “us”	This line and the above have identical indentation.
print “finished”		Aligned with “if” so outside the block.

Change the first line to x=3 and run again.

The comparison operators most often used in if statements are ==, >, >=, <, <=, !=. Multiple tests use elif (meaning “else if”) and else constructs:

x=5
if x<4 or x>27:
	print “Bring”
elif x<7:		# Notice the colon
	print “us”
elif x<12:
	print “a”
else:
	print “shrubbery”
		# Can end a block with a blank line.

Try changing x to 5, 3, 8, 20, and 30. Make sure you understand the results.

All Python blocks have the same basic structure, a line ending with a colon followed by a block of indented code. The amount of indentation is arbitrary, but be careful to remain consistent or you’ll get an error.

The basic loop is the for loop:

for k in [ 1,2,3 ]:		# Notice the keyword “in” and the colon.
	print “k=”,k

In the above example the iteration is over the members of a list. In fact, a for loop can iterate over any sequence (a list, tuple, string, or array).

for k in [“I” , “feel” , “fine” ]:
	print k
for k in 1,2,3:		# A tuple
	print k
for k in “swallow”:		# A string
	print k
for k in arange(5):		# An array. To use this in a module you need
	print k	# to put this as the first line of the module:
		from Numeric import *
for k in range(5):		# A list
	print k

The last two examples are the most important for numerical programming. They show quick ways to iterate over a sequence of numbers. Both arange() and range() permit arbitrary starting points, ending points, and stride. Moreover, in arange() these don’t have to be integers:

print range(-10,10,2)
for x in arange(-pi/2,pi/2,pi/4):
	print x/pi

Notice again that the upper limit in the range is not included.

Another type of loop is provided by while:

a,b = 0,1
while b<10:		# Notice the colon.
	print a,b
	a,b = b,a+b	# Cute multiple assignment.

Sometimes within loops you have if tests and, depending on the outcome, want to break out of the loop or instead continue the next iteration. The