What follows below is a short and very basic practical introduction to python/numpy. This is an informal introduction – you can find a large amount of information on the web. I have used information from the following websites:
Python is a high-level, dynamically typed multiparadigm programming language. Python code is almost like pseudocode. It allows you to express very powerful ideas in very few lines of code while being very readable. It is a great general-purpose programming language on its own, but with the help of a few popular libraries (numpy, scipy, matplotlib) it becomes a powerful environment for scientific computing. Some of you may already have some experience with Python and numpy; for the rest of you, this page will serve as a quick crash course on the use of Python for scientific computing.
For detailed documentation, please visit numpy, scipy and matplotlib.
Basics
Python (even without numpy) is already capable of performing numerical calculations in an easy and straightforward way. It is text-based so you type in commands, hit return, and this creates output. For example, if you type 2+2 and hit return, you get 4.
>>> 2+2
4
Operators
The usual mathematical operations of interest are +, -, * (multiply), / (divide), and ** (exponentiation). So,
>>> 2**3
8
>>> (3.0+4.0)/14.0
0.5
These can be made to act on numpy arrays as well (see below).
Variables
Python spits out an the answer in the two examples above in the next line. The variable in which it holds the answer to the calculation is _. If you type the name of a variable by itself and hit return, python prints its value (e.g. type _ and see its value). One can create and assign variables. For example, typing x=2**3 will result in 8. It is not printed out on the screen – but, the variable x is defined and currently has value 8. If you type x by itself, python will print its value.
>>> x=2**3
>>> x
8
Complex numbers and complex conjugation
To create complex numbers use j or J as $\sqrt{-1}$. Complex conjugation can be using the function conjugate:
>>> z = 2+3j
>>> z.conjugate()
2-3j
Numpy
Numpy is the core library for scientific computing in Python. It provides a high-performance multidimensional array object, and tools for working with these arrays. To use numpy, one has to first import it.
import numpy as np
This allows one to use all the functions of numpy (but instead of typing numpy, you can just type np)
Arrays
A numpy array is a grid of values, all of the same type, and is indexed by a tuple of non-negative integers. The number of dimensions is the rank of the array; the shape of an array is a tuple of integers giving the size of the array along each dimension.
import numpy as np
a = np.array([1, 2, 3]) # Create a rank 1 array
print type(a) # Prints "<type 'numpy.ndarray'>"
print a.shape # Prints "(3,)"
print a[0], a[1], a[2] # Prints "1 2 3"
a[0] = 5 # Change an element of the array
print a # Prints "[5, 2, 3]"
b = np.array([[1,2,3],[4,5,6]]) # Create a rank 2 array
print b.shape # Prints "(2, 3)"
print b[0, 0], b[0, 1], b[1, 0] # Prints "1 2 4"
Transposition and Hermitian conjugation
Transpose of a numpy array is easily done in numpy using transpose function:
import numpy as np
a = np.array([1, 2, 3]) # Create a rank 1 array
print a.shape # Prints "(3,)"
b = a.transpose() # This does not change anything as you cannot transpose a 1D array!
print b.shape # Prints "(3,)"
a = array([[0, 1], [2, 3]]) # Create a rank 2 array
print a # Prints "[[0, 1]
# [2, 3]]"
b = a.transpose() # generates the transpose
c = a.T # This does the same thing
print b # Prints "[[0, 2]
# [1, 3]]"
For Hermitian conjugate – just use a.conj().T
Creating an array of equally spaced numbers
Creating a numpy array of numbers can be achieved by using loops. But the simplest way is to use numpy.arange(start, stop, step) which returns equally spaced values in the open half interval [start, stop) with step size – step
import numpy as np
np.arange(1, 5) # Prints "array([1, 2, 3, 4])"
np.arange(1, 6, 2) # Prints "array([1, 3, 5])"
Creating an array of zeros or ones or any constant
import numpy as np
a = np.zeros((2,2)) # Create an array of all zeros
print a # Prints "[[ 0. 0.]
# [ 0. 0.]]"
b = np.ones((1,2)) # Create an array of all ones
print b # Prints "[[ 1. 1.]]"
c = np.full((2,2), 3) # Create a constant array
print c # Prints "[[ 3. 3.]
# [ 3. 3.]]"
Creating an identity matrix or one full of random numbers
import numpy as np
d = np.eye(2) # Create a 2x2 identity matrix
print d # Prints "[[ 1. 0.]
# [ 0. 1.]]"
e = np.random.random((2,2)) # Create an array filled with random values
print e # Might print "[[ 0.73168019, 0.6214697 ]
# [ 0.64961548, 0.54018523]]"
Array indexing
Numpy offers several ways to index into arrays. Numpy arrays can be sliced. Since arrays may be multidimensional, you must specify a slice for each dimension of the array:
import numpy as np
# Create the following rank 2 array with shape (3, 4)
# [[ 1 2 3 4]
# [ 5 6 7 8]
# [ 9 10 11 12]]
a = np.array([[1,2,3,4], [5,6,7,8], [9,10,11,12]])
# Use slicing to pull out the subarray consisting of the first 2 rows
# and columns 1 and 2; b is the following array of shape (2, 2):
# [[2 3]
# [6 7]]
b = a[:2, 1:3]
# A slice of an array is a view into the same data, so modifying it
# will modify the original array.
print a[0, 1] # Prints "2"
b[0, 0] = 77 # b[0, 0] is the same piece of data as a[0, 1]
print a[0, 1] # Prints "77"
Creating a diagonal array
Numpy offers two functions to create and extract diagonals of arrays nump.diag and numpy.diagonal. Here we will discuss the numpy.diag.
import numpy as np
np.diag([1, 2, 3])
# Prints "array([[1, 0, 0],
# [0, 2, 0],
# [0, 0, 3]])"
np.diag(np.full(3), 2.0)
# Prints "array([[2.0, 0, 0],
# [0, 2.0, 0],
# [0, 0, 2.0]])"
Creating offdiagonal/triangular arrays
Often, one wants to fill the non-diagonal entries in an array in a rather regular way. For example, one may want to fill the diagonal above or below the main diagonal of a matrix with some numbers. The numpy.diag function performs this when we give it an extra argument. For example, let us take a 3x3 diagonal matrix with the numbers [10.0 11.0 12.0] on its diagonal and make its upper off diagonal have the numbers [-3.5 -4.5]:
import numpy as np
M = np.diag([10.0, 11.0, 12.0]) + np.diag([-3.5, -4.5],1)
# Prints "array([[ 10. , -3.5, 0. ],
# [ 0. , 11. , -4.5],
# [ 0. , 0. , 12. ]])
For the lower diagonal use M = np.diag([10.0, 11.0, 12.0]) + np.diag([-3.5, -4.5],-1)
Array Math
Basic mathematical functions operate elementwise on arrays, and are available both as operator overloads and as functions in the numpy module:
import numpy as np
x = np.array([[1,2],[3,4]], dtype=np.float64)
y = np.array([[5,6],[7,8]], dtype=np.float64)
# Elementwise sum; both produce the array
# [[ 6.0 8.0]
# [10.0 12.0]]
print x + y
print np.add(x, y)
# Elementwise difference; both produce the array
# [[-4.0 -4.0]
# [-4.0 -4.0]]
print x - y
print np.subtract(x, y)
# Elementwise product; both produce the array
# [[ 5.0 12.0]
# [21.0 32.0]]
print x * y
print np.multiply(x, y)
# Elementwise division; both produce the array
# [[ 0.2 0.33333333]
# [ 0.42857143 0.5 ]]
print x / y
print np.divide(x, y)
# Elementwise square root; produces the array
# [[ 1. 1.41421356]
# [ 1.73205081 2. ]]
print np.sqrt(x)
For matrix multiplication, one has to use the numpy.dot function. This can be used to compute inner products of vectors, to multiply a vector by a matrix, and to multiply matrices. dot is available both as a function in the numpy module and as an instance method of array objects:
import numpy as np
x = np.array([[1,2],[3,4]])
y = np.array([[5,6],[7,8]])
v = np.array([9,10])
w = np.array([11, 12])
# Inner product of vectors; both produce 219
print v.dot(w)
print np.dot(v, w)
# Matrix / vector product; both produce the rank 1 array [29 67]
print x.dot(v)
print np.dot(x, v)
# Matrix / matrix product; both produce the rank 2 array
# [[19 22]
# [43 50]]
print x.dot(y)
print np.dot(x, y)
Inverting or diagonalizing an array
The linear algebra package offers methods to invert an array and find its eigenvalues.
import numpy as np
A = np.diag([10.0, 11.0, 12.0]) + np.diag([3, 4], 1) + np.diag([3, 4], -1)
np.linalg.inv(A) # Calculate the inverse
# Prints "array([[ 0.11026616, -0.03422053, 0.01140684],
# [-0.03422053, 0.11406844, -0.03802281],
# [ 0.01140684, -0.03802281, 0.0960076 ]])"
np.linalg.eig(A) # Calculate the eigenvalues and eigenvectors
# Prints "(array([ 6.041338 , 10.7300123 , 16.22864971]), array([[ 0.5325576 , 0.77929762, 0.330269 ],
# [-0.70273852, 0.18963228, 0.68570998],
# [ 0.47174249, -0.59727281, 0.64863257]]))"
# where the first array is eigenvalues and the next array is eigenvectors.
Testing the orthogonality of eigenvectors
import numpy as np
A = np.diag([10.0, 11.0, 12.0]) + np.diag([3, 4], 1) + np.diag([3, 4], -1)
[e,V] = np.linalg.eig(A) # Now V contains the eigenvectors and e the eigenvalues
np.dot(V.T, V)
# Prints "array([[ 1.00000000e+00, 0.00000000e+00, -5.55111512e-17],
# [ 0.00000000e+00, 1.00000000e+00, 1.66533454e-16],
# [ -5.55111512e-17, 1.66533454e-16, 1.00000000e+00]])"
Creating an outer product
import numpy as np
A = np.diag([10.0, 11.0, 12.0]) + np.diag([3, 4], 1) + np.diag([3, 4], -1)
[e,V] = np.linalg.eig(A) # Now V contains the eigenvectors and e the eigenvalues
np.dot(np.dot(h,np.diag(e)),h.T) # Outer product including the diagonal -- V e V^T
# This produces the matrix A
# Prints "array([[ 1.00000000e+01, 3.00000000e+00, 8.88178420e-16],
# [ 3.00000000e+00, 1.10000000e+01, 4.00000000e+00],
# [ 4.44089210e-16, 4.00000000e+00, 1.20000000e+01]])"
Loops
While in principle a lot of what has been written above can be achieved with simple loops, it is best to avoid them when using numpy. Loops in numpy lead to slow code. In some cases, loops are unavoidable. If one has to loop the simplest way is
for i in range(0,2):
print i
# Prints "0
# 1"
File handling
File handling in python is similar to C
fp = open('workfile', 'r') # This opens a file 'workfile' for reading
line = fp.readline() # This reads the first line of the file
lines = fp.readlines() # This reads all the rest of the lines in the file
for line in lines: # Iterates over the 'lines' (read from the file)
... do something ...
fp.close() # Close the file
General and file IO
The formatting for input and output is given very well at Python IO formatting
Functions
Python functions are defined using the def keyword. For example:
def sign(x):
if x > 0:
return 'positive'
elif x < 0:
return 'negative'
else:
return 'zero'
for x in [-1, 0, 1]:
print sign(x)
# Prints "negative", "zero", "positive"
Functions can get optional arguments as shown below
def hello(name, loud=False):
if loud:
print 'HELLO, %s!' % name.upper()
else:
print 'Hello, %s' % name
hello('Bob') # Prints "Hello, Bob"
hello('Fred', loud=True) # Prints "HELLO, FRED!"
More numpy exercises
To learn more about numpy, a good option is to go over the problem set at 100 Numpy problems.