Python for Scientists by Chris Seddon Copyright 2000

Python for Scientists 1. 2. 3. 4. Num. Py Plotting with gnuplot Plotting with

Chapter 1 1 Copyright © 2000 -9 CRS Enterprises Ltd 5

Num. Py multi-dimensional arrays creating arrays reshaping arrays slicing and iterating stacking arrays splitting arrays shallow and deep copies broadcasting indexing grids matrices Copyright © 2000 -9 CRS Enterprises Ltd 6

01 creating arrays several different ways to create an array a = empty( (3, 5) ) a = ones( (3, 5) ) a = zeros( (3, 5) ) a = array( [2, 3, 5, 7, 11, 13, 17] ) # create empty 3 x 5 array # filled with 1's # filled with 0's # from a list # create 1 D array using arange a = arange(4, 17) # create equally spaced 1 D array a = linspace(-50. 0, 5) # use a function a = fromfunction(lambda i, j, k: i + j + k, (4, 2, 3)) Copyright © 2000 -9 CRS Enterprises Ltd 7

02 reshaping arrays numpy. reshape(a, newshape, order='C') gives a new shape to an array without changing its data C or Fortran style arrays available Fortran stores arrays in row order with the first subscript increasing most rapidly C/C++ in column order with the last subscript increasing most rapidly # create array a = ones( (24, ) ) # reshape it b = a. reshape(2, 3, 4) Copyright © 2000 -9 CRS Enterprises Ltd 8

03 array types dtype defines the type of ndarray bool_ uint 8 uint 16 uint 32 uint 64 int 8 int 16 int 32 int 64 float 32 float 64 float 96 complex 64 complex 128 complex 192 # create array of int a = array( [2, 3, 5, 7, 11, 13, 17] ) # create array of float a = array( [2, 3, 5, 7, 11, 13, 17], dtype="float 64" ) # create array of complex a = array( [2, 3, 5, 7, 11, 13, 17], dtype="complex 64" ) Copyright © 2000 -9 CRS Enterprises Ltd 9

04 operating on arrays operations are performed on each and every element arrays are mutable # operations are performed on each element c=a*b+2 # dot and cross product c = dot(a, b) c = cross(a, b) Copyright © 2000 -9 CRS Enterprises Ltd 10

05 array methods Very few methods. . . sum() max() min() return sum of all elements return largest element return smallest element a = array( random(10) * 100, dtype="int" ) print "sum=", a. sum() print "max=", a. max() print "min=", a. min() Copyright © 2000 -9 CRS Enterprises Ltd 11

06 basic functions Although arrays are mutable. . . these methods return new arrays # perform operations along specified axis result = average(a, axis=0) # create new arrays b = msort(a) b = insert(a, (1, 5, 9), (-1, -2, -3) ). reshape(3, 5) b = append(a, (-1, -2, -3, -4) ). reshape(4, 4) b = round_(a ** 1. 1, 2) Copyright © 2000 -9 CRS Enterprises Ltd 12

07 slicing and iterating Slicing similar to regular Python Iteration only works on first dimension but you can use flat to convert to a 1 D array # one dimensional arrays a = arange(20); print a[7: 14] # multi-dimensional arrays a = arange(24). reshape(4, 3, 2) print a[0: 2, 0: 2] # iterate (works on first dimension) for row in a: print row # flatten array (not a function) to iterate over every element for element in a. flat: print element, Copyright © 2000 -9 CRS Enterprises Ltd 13

08 stacking arrays Arrays can be stacked. . . horizontally or vertically a = ones( (3, 5) ) b = zeros( (3, 5) ) # stacking h = hstack( (a, b) ) v = vstack( (a, b) ) # create a 1 D array a = r_[10: 20, 55, 66, 80: 90] Copyright © 2000 -9 CRS Enterprises Ltd 14

09 splitting arrays Splitting arrays. . . by each row (1 st dimension) by each column (lat dimension) by unequal amounts (as specified by a tuple) a = array( random((3, 5)) * 100 ); # split horizontally a 0, a 1, a 2, a 3, a 4 = hsplit(a, 5) # split vertically a 0, a 1, a 2 = vsplit(a, 3) # split unequally horizontally" a 0, a 1, a 2 = hsplit(a, (1, 2)) # 2 splits Copyright © 2000 -9 CRS Enterprises Ltd 15

10 shallow copies Views are shallow copies the new array is a different from the old array, but the underlying data is shared a = arange(12) c = a. view() # a and c point to different arrays that share data # reshape one array - data unchanged c. shape = 3, 4 # slice creates a view - data unchanged d = a[2: 5] # d sees the change made by a a[3] = 88 Copyright © 2000 -9 CRS Enterprises Ltd 16

11 deep copies Use copy() to. . . create a new array and copy

12 broadcasting Allows operations to be performed on different sized arrays Extra dimensions are added to smaller arrays each dimension has size 1 The new dimensions are then adjusted to the size of the larger array by broadcasting values duplicated across new columns a = array([100, 200, 300, 400]) b = arange(16). reshape(4, 4) # a has shape = 4 # converted to 1, 4 # then broadcasted to 4, 4 c = a + b; Copyright © 2000 -9 CRS Enterprises Ltd 18

13 indexing Indexing is an advanced form of slicing used to select items from an array by index or by expressions # set up an array to be used in indexing a = arange(10)**2 # setup index array index = array( [2, 3, 5, 9] ) # apply index to a (to extract items 2, 3, 5 and 9 b = a[index] # set up a boolean filter for a a = arange(24). reshape(6, 4) filter = a % 3 == 0 a[filter] = 99 # only change every 3 rd element Copyright © 2000 -9 CRS Enterprises Ltd 19

14 ix function. . . used to create multi-dimensional grids useful in 3 D plotting # set up arrays to be used for 3 D grid a = array([2, 5, 7]) b = array([4, 8, 1]) c = array([3, 5, 2]) ax, bx, cx = ix_(a, b, c) result = ax * (bx + cx) # these are equivalent print result[1, 2, 1] print a[1] * (b[2] + c[1]) Copyright © 2000 -9 CRS Enterprises Ltd 20

15 grids Alternative to ix_function x = arange(1, 11) y = arange(1, 11) # Make a 2 -d array containing a function of x and y. First create # xm and ym which contain the x and y values in a matrix form that # can be `broadcast' into a matrix of the appropriate shape: xm = x[ : , newaxis ] ym = y[ newaxis , : ] # generate 2 D grid of values m = xm * ym Copyright © 2000 -9 CRS Enterprises Ltd 21

16 matrices Subclass of Array # matrix multiplication (not element wise multiplication) a = matrix( [[3, 4, 5], [2, 3, 8], [4, 1, 7]] ) b = matrix( [[2, 3, 4], [1, 2, 7], [3, 0, 6]] ) c=a*b a = matrix( [[3, 5], [4, 1]] ) print "Transpose", a. T print "Inverse", a. I # linear algebra # 5 x + 3 y = 31 # 2 x - 7 y = -45 # solution x=2, y=7 a = matrix( [[ 5, 3 ], [ 2, -7 ]] ) v = matrix( [[ 31 ], [ -45 ]] ) print linalg. solve(a, v) Copyright © 2000 -9 CRS Enterprises Ltd 22

Chapter 2 2 Copyright © 2000 -9 CRS Enterprises Ltd 25

Plotting with gnuplotting data points using data files multiple plots vector fields hardcopies surface

01 data plot Supply list of 2 D or 3 D points import Gnuplot gp = Gnuplot() gp('set data style lines') # Pairs of x and y values data = [[0, 0], [1, 1], [2, 4], [3, 9], [4, 16]] gp. plot(data) # Triplets of x, y, z values data = [[0, 0, 0], [1, 1, 1], [2, 4, 8], [3, 9, 27], [4, 16, 64]] gp. splot(data) Copyright © 2000 -9 CRS Enterprises Ltd 27

02 data files Read data from a file import Gnuplot gp = Gnuplot() f = open("data/f 1. dat") # read in the data points data = [] for line in f: point = line. split() data. append(point) # plot them gp. plot(data) Copyright © 2000 -9 CRS Enterprises Ltd 28

03 multiple plots Multiple plots on same graph import Gnuplot gp = Gnuplot() gp('set data style lines') data 1 = [[0, 0], [1, 1], [2, 4], [3, 9], [4, 16]] # The first data set (a quadratic) data 2 = [[0, 0], [1, 1], [2, 2], [3, 3], [4, 4]] # The second data set (a straight line) # Make the plot items plot 1 = Gnuplot. Plot. Items. Data(data 1, with_="lines", title="Quadratic") plot 2 = Gnuplot. Plot. Items. Data(data 2, with_="points 3", title="Linear") gp. plot(plot 1, plot 2) Copyright © 2000 -9 CRS Enterprises Ltd 29

04 vector fields Read data from a file import Gnuplot gp = Gnuplot() # Generate the vector field data = [] for i in range(-10, 11): for j in range(-10, 11): data. append([i/10. , j/10. , -i/100. , -j/100. ]) # Each data entry has the form [x, y, u, v] # Plot plot = Gnuplot. Plot. Items. Data(data, with_ = 'vectors') gp. plot(plot) Copyright © 2000 -9 CRS Enterprises Ltd 30

05 hardcopies gnuplot can convert plots into hardcopy but only postscript available import Gnuplot gp = Gnuplot(). . . gp. plot(plot) gp. hardcopy(filename="data/mygraph. ps", terminal="postscript") Copyright © 2000 -9 CRS Enterprises Ltd 31

06 delayed plotting Plots can be cached in memory displayed later with plot() from numpy import * import Gnuplot, math gp = Gnuplot(debug=1) x = arange(10, dtype='float_') y = x**2 # save plots for displaying later plot 1 = Gnuplot. Data(x, y, title='squares', with_='points 3 3') plot 2 = Gnuplot. Func('x**2', title='calculated by gnuplot') # Plot gp. title('The function x**2') gp. xlabel('x'); gp. ylabel('x**2') gp. plot(plot 1, plot 2) Copyright © 2000 -9 CRS Enterprises Ltd 32

07 plotting styles Many different plotting styles available mport Gnuplot import math gp = Gnuplot(debug=1) x = arange(10, dtype='float_') y = x**2 # save plots for displaying later plot 1 = Gnuplot. Func('x**2/100. ', title='x**2', with_='boxes') plot 2 = Gnuplot. Func('cos(x)', title='cos(x)', with_='lines') plot 3 = Gnuplot. Func('norm(x)', title='norm(x)', with_='errorbars') # Plot gp. title('Various functions and styles') gp. xlabel('x'); gp. ylabel('y') gp. plot(plot 1, plot 2, plot 3) Copyright © 2000 -9 CRS Enterprises Ltd 33

08 surface plots 3 D surface plots available can rotate views with mouse import Gnuplot, math g = Gnuplot(debug=1) # Demonstrate a 3 -d plot: x = arange(25) y = arange(25) # Make a 2 -d array containing a function of x and y xm = x[: , newaxis] ym = y[newaxis, : ] m = xm * ym**2 g('set parametric'); g('set data style lines'); g('set hidden') g('set contour base'); g. title('An example of a surface plot') g. xlabel('x'); g. ylabel('y') g. splot(Gnuplot. Grid. Data(m, x, y)) Copyright © 2000 -9 CRS Enterprises Ltd 34

Chapter 3 3 Copyright © 2000 -9 CRS Enterprises Ltd 37

Plotting with matplotlib 2 D plots only data points polynomials multiple plots subplots polar

01 simple plot Supply x and y data. . . as separate lists import matplotlib. pyplot as plt red. Circles = "ro" plt. plot([1, 2, 3, 4], [1, 4, 9, 16], red. Circles) plt. ylabel("squares") plt. show() Copyright © 2000 -9 CRS Enterprises Ltd 39

02 polynomials Plotting functiond is easy! from matplotlib. mlab import normpdf import matplotlib. numerix as nx import pylab as p """ calculate polynomial: y = -2 x^2 + 1 """ def poly(input): output = [(-2 * a + 1) for a in input] return output x = nx. arange(-4, 4, 0. 01) y = poly(x) p. plot(x, y, color='red', lw=2) p. show() Copyright © 2000 -9 CRS Enterprises Ltd 40

03 multiple plots Put several plot on one graph with different formats import numpy as np import matplotlib. pyplot as plt # evenly sampled time at 200 ms intervals t = np. arange(0. 0, 5. 0, 0. 200) red. Dashes = "r--" blue. Squares = "bs" green. Triangles = "g^" plt. plot( t, t, red. Dashes, t, t**2, blue. Squares, t, t**3, green. Triangles) plt. show() Copyright © 2000 -9 CRS Enterprises Ltd 41

04 subplots Combine different graphs on same figure import numpy as np import matplotlib. pyplot as plt def f(t): return np. exp(-t) * np. cos(2*np. pi*t) t 1 = np. arange(0. 0, 5. 0, 0. 1) t 2 = np. arange(0. 0, 5. 0, 0. 02) plt. figure(1) plt. subplot(211) # 2 rows, 1 column, fig. 1 plt. plot(t 1, f(t 1), 'bo', t 2, f(t 2), 'k') plt. subplot(212) # 2 rows, 1 column, fig. 2 plt. plot(t 2, np. cos(2*np. pi*t 2), 'r--') plt. show() Copyright © 2000 -9 CRS Enterprises Ltd 42

05 polar plots import numpy as np import matplotlib. pyplot as plt fig = plt. figure() ax = fig. add_subplot(111, polar=True) r = np. arange(0, 1, 0. 001) theta = 2*2*np. pi*r line, = ax. plot(theta, r, color="#ee 8 d 18", lw=3) ind = 800. . . thisr, thistheta = r[ind], theta[ind] ax. annotate( ax. plot([thistheta], [thisr], "o") "a polar annotation", . . . xy=(thistheta, thisr), # theta, radius xytext=(0. 05, 0. 05), # fraction, fraction textcoords="figure fraction", arrowprops=dict(facecolor="black", shrink=0. 05), horizontalalignment="left", verticalalignment="bottom", ) plt. show() Copyright © 2000 -9 CRS Enterprises Ltd 43

06 histograms Use the hist() function import numpy as np import matplotlib. pyplot as plt mu, sigma = 100, 15 x = mu + sigma * np. random. randn(10000) # the histogram of the data n, bins, patches = plt. hist(x, 50, normed=1, facecolor='g', alpha=0. 75) plt. xlabel('Smarts') plt. ylabel('Probability') plt. title('Histogram of IQ') plt. text(60, . 025, r'$mu=100, sigma=15$') plt. axis([40, 160, 0, 0. 03]) plt. grid(True) plt. show() Copyright © 2000 -9 CRS Enterprises Ltd 44

Chapter 4 4 Copyright © 2000 -9 CRS Enterprises Ltd 47

Sci. Py science constants integration linear equations special functions solving non-linear equations ordinary differential

Sci. Py Packages Clustering package scipy. cluster Constants scipy. constants Fourier transforms scipy. fftpack Integration and ODEs scipy. integrate Interpolation scipy. interpolate Input and output scipy. io Linear algebra scipy. linalg Miscellaneous routines scipy. misc Multi-dimensional image processing scipy. ndimage Optimization and root finding scipy. optimize Signal processing scipy. signal Sparse matrices scipy. sparse Sparse linear algebra scipy. sparse. linalg Special functions scipy. special Statistical functions scipy. stats Copyright © 2000 -9 CRS Enterprises Ltd 49

01 constants Many mathematical and physical constants defined #SI prefixes print c. giga, c. mega, c. kilo # physical constants print "speed of light", c. c print "Plank's constant", c. h print "Boltzmann constant", c. k #weight in kg print "gram", c. gram print "lb", c. lb print "electron mass", c. m_e #length in meter print "inch", c. inch print "yard", c. yard print "mile", c. mile Copyright © 2000 -9 CRS Enterprises Ltd 50

02 miscellaneous Miscellaneous functions group in misc module from numpy import * from scipy import misc # factorial print "factorial(23) =", misc. factorial(23) # differentiation print "derivative of x**5 at x=5", print misc. derivative(lambda x: x**5, 5, dx=0. 0001) # combinations print "2 from 5 = ", misc. comb(5, 2) print "13 from 52 = ", misc. comb(52, 13) Copyright © 2000 -9 CRS Enterprises Ltd 51

03 integration Integration using quadrature from scipy import * from scipy. integrate import quad def f(x): return exp(-x) # integrate exp(-x) between 0 and 5 result = quad(f, 0, 5) print "result", result[0] print "error estimate", result[1] # more tests print quad(lambda x: x, 0, 10) print quad(lambda x: x**2, 0, 1) Copyright © 2000 -9 CRS Enterprises Ltd 52

04 linear equations from numpy import * from scipy. linalg import solve # x = 5, y = 3, z = 1 # 5 x - 2 y + 3 z = 22 # 6 x + y + z = 34 # 2 x - 5 y - 2 z = -7 a = array( [ [5, -2, 3], [6, 1, 1], [2, -5, -2] ]) b = array( [22, 34, -7] ) result = solve(a, b) print result Copyright © 2000 -9 CRS Enterprises Ltd 53

05 special functions Many special functions available from numpy import * from scipy import special # Bessel function of real order v at complex z print special. jv(6, 2+3 j) # Gamma function print special. gamma(5) print special. gamma(5. 1) # Zeta function print special. zeta(10, 2) Copyright © 2000 -9 CRS Enterprises Ltd 54

06 solving non linear equations Must supply an initial guess may be multiple roots finds closest root from scipy import * from scipy. optimize import fsolve def f(x): return (exp(x)-1)-cos(x) # use fsolve to solve exp(x)-1 - cos(x) = 0 guess = 2. 0; print fsolve(f, guess) guess = -1. 0; print fsolve(f, guess) guess = -5. 0; print fsolve(f, guess) guess = 0. 0; print fsolve(f, guess) guess = -10. 0; print fsolve(f, guess) Copyright © 2000 -9 CRS Enterprises Ltd 55

07 ordinary differential equations scipy. integrate. odeint(func, y 0, t, . . . ) dy/dt = func(y, t 0, . . . ) func : computes the derivative of y at t 0. y 0 : initial condition on y (can be a vector). t : sequence of time points for which to solve for y from numpy import * from scipy. integrate import odeint from math import * # evaluate dy/dt = y for initial value of y = y 0 # and a time range from t 0. . . t 1 y 0 = 7; t 0 = 2; t 1 = 6 # solution: y = e**(t + k) # where k = log(y 0) - t 0 print "calculated result", odeint(lambda y, t: y, y 0, [t 0, t 1]) k = log(y 0) - t 0 print "analytical result", [ y 0, e**(t 1 + k)] Copyright © 2000 -9 CRS Enterprises Ltd 56