580 likes | 594 Views
Python for Scientists. by Chris Seddon. Python for Scientists. 1. NumPy 2. Plotting with gnuplot 3. Plotting with matplotlib 4. SciPy. 1. Chapter 1. NumPy. multi-dimensional arrays creating arrays reshaping arrays slicing and iterating stacking arrays splitting arrays
E N D
Python for Scientists byChris Seddon
Python for Scientists • 1. NumPy • 2. Plotting with gnuplot • 3. Plotting with matplotlib • 4. SciPy
1 Chapter 1
NumPy • multi-dimensional arrays • creating arrays • reshaping arrays • slicing and iterating • stacking arrays • splitting arrays • shallow and deep copies • broadcasting • indexing • grids • matrices
01 creating arrays • several different ways to create an array a = empty( (3,5) ) # create empty 3 x 5 array a = ones( (3,5) ) # filled with 1's a = zeros( (3,5) ) # filled with 0's a = array( [2,3,5,7,11,13,17] ) # from a list # create 1D array using arange a = arange(4,17) # create equally spaced 1D array a = linspace(-50.0,50.0,5) # use a function a = fromfunction(lambda i,j,k: i + j + k, (4,2,3))
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)
03 array types • dtype defines the type of ndarray • bool_ uint8 uint16 uint32 uint64 int8 int16 int32 int64 float32 float64 float96 complex64 complex128 complex192 # 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="float64" ) # create array of complex a = array( [2,3,5,7,11,13,17], dtype="complex64" )
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)
05 array methods • Very few methods ... • sum() return sum of all elements • max() return largest element • min() return smallest element a = array( random.random(10) * 100, dtype="int" ) print "sum=", a.sum() print "max=", a.max() print "min=", a.min()
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)
07 slicing and iterating • Slicing similar to regular Python • Iteration only works on first dimension • but you can use flat to convert to a 1D 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,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,
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 1D array a = r_[10:20, 55, 66, 80:90]
09 splitting arrays • Splitting arrays ... • by each row (1st dimension) • by each column (lat dimension) • by unequal amounts (as specified by a tuple) a = array( random.random((3,5)) * 100 ); # split horizontally a0,a1,a2,a3,a4 = hsplit(a,5) # split vertically a0,a1,a2 = vsplit(a,3) # split unequally horizontally" a0,a1,a2 = hsplit(a,(1,2)) # 2 splits
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
11 deep copies • Use copy() to ... • create a new array • and copy the underlying data a = arange(12) b = a.copy() # changes to a do not affect b a[3] = 99
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;
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 3rd element
14 ix function • ix function ... • used to create multi-dimensional grids • useful in 3D plotting # set up arrays to be used for 3D 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])
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 2D grid of values m = xm * ym
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 # 5x + 3y = 31 # 2x - 7y = -45 # solution x=2, y=7 a = matrix( [[ 5, 3 ],[ 2, -7 ]] ) v = matrix( [[ 31 ],[ -45 ]] ) print linalg.solve(a,v)
2 Chapter 2
plotting data points using data files multiple plots vector fields hardcopies surface plots Plotting with gnuplot
01 data plot • Supply list of 2D or 3D points import Gnuplot gp = Gnuplot.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)
02 data files • Read data from a file import Gnuplot gp = Gnuplot.Gnuplot() f = open("data/f1.dat") # read in the data points data = [] for line in f: point = line.split() data.append(point) # plot them gp.plot(data)
03 multiple plots • Multiple plots on same graph import Gnuplot gp = Gnuplot.Gnuplot() gp('set data style lines') data1 = [[0, 0], [1, 1], [2, 4], [3, 9], [4, 16]] # The first data set (a quadratic) data2 = [[0, 0], [1, 1], [2, 2], [3, 3], [4, 4]] # The second data set (a straight line) # Make the plot items plot1 = Gnuplot.PlotItems.Data(data1, with_="lines", title="Quadratic") plot2 = Gnuplot.PlotItems.Data(data2, with_="points 3", title="Linear") gp.plot(plot1, plot2)
04 vector fields • Read data from a file import Gnuplot gp = Gnuplot.Gnuplot() # Generate the vector field data 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.PlotItems.Data(data, with_ = 'vectors') gp.plot(plot)
05 hardcopies • gnuplot can convert plots into hardcopy • but only postscript available import Gnuplot gp = Gnuplot.Gnuplot() ... gp.plot(plot) gp.hardcopy(filename="data/mygraph.ps",terminal="postscript")
06 delayed plotting • Plots can be cached in memory • displayed later with plot() from numpy import * import Gnuplot, math gp = Gnuplot.Gnuplot(debug=1) x = arange(10, dtype='float_') y = x**2 # save plots for displaying later plot1 = Gnuplot.Data(x, y, title='squares', with_='points 3 3') plot2 = Gnuplot.Func('x**2', title='calculated by gnuplot') # Plot gp.title('The function x**2') gp.xlabel('x'); gp.ylabel('x**2') gp.plot(plot1, plot2)
07 plotting styles • Many different plotting styles available mport Gnuplot import math gp = Gnuplot.Gnuplot(debug=1) x = arange(10, dtype='float_') y = x**2 # save plots for displaying later plot1 = Gnuplot.Func('x**2/100.', title='x**2', with_='boxes') plot2 = Gnuplot.Func('cos(x)', title='cos(x)', with_='lines') plot3 = Gnuplot.Func('norm(x)', title='norm(x)', with_='errorbars') # Plot gp.title('Various functions and styles') gp.xlabel('x'); gp.ylabel('y') gp.plot(plot1, plot2, plot3)
08 surface plots • 3D surface plots available • can rotate views with mouse import Gnuplot, math g = Gnuplot.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.GridData(m,x,y))
3 Chapter 3
2D plots only data points polynomials multiple plots subplots polar plots histograms Plotting with matplotlib
01 simple plot • Supply x and y data ... • as separate lists import matplotlib.pyplot as plt redCircles = "ro" plt.plot([1,2,3,4], [1,4,9,16], redCircles) plt.ylabel("squares") plt.show()
02 polynomials • Plotting functiond is easy! from matplotlib.mlab import normpdf import matplotlib.numerix as nx import pylab as p """ calculate polynomial: y = -2x^2 + 1 """ def poly(input): output = [(-2 * a * 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()
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 200ms intervals t = np.arange(0.0, 5.0, 0.200) redDashes = "r--" blueSquares = "bs" greenTriangles = "g^" plt.plot( t, t, redDashes, t, t**2, blueSquares, t, t**3, greenTriangles) plt.show()
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) t1 = np.arange(0.0, 5.0, 0.1) t2 = np.arange(0.0, 5.0, 0.02) plt.figure(1) plt.subplot(211) # 2 rows, 1 column, fig.1 plt.plot(t1, f(t1), 'bo', t2, f(t2), 'k') plt.subplot(212) # 2 rows, 1 column, fig.2 plt.plot(t2, np.cos(2*np.pi*t2), 'r--') plt.show()
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="#ee8d18", lw=3) ind = 800 thisr, thistheta = r[ind], theta[ind] ax.plot([thistheta], [thisr], "o") ... ... ax.annotate( "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()
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()
4 Chapter 4
SciPy • science constants • integration • linear equations • special functions • solving non-linear equations • ordinary differential equations
SciPy 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
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