Scientific Programming I
Berian James / [email protected]

1. Introduction to SciPy

Import NumPy and SciPy

In [ ]:
import numpy as np
In [ ]:
import scipy as sp

Differences in temperament between NumPy and SciPy; see, e.g., http://www.scipy.org/NegativeSquareRoot

In [ ]:
exp(pi*np.sqrt(-1)) + 1
In [ ]:
exp(pi*sp.sqrt(-1)) + 1

Getting data in and out of SciPy

Import Matlab data into Python

In [ ]:
import scipy.io as sio
In [ ]:
struct = sio.loadmat('testbox.mat') # Imports Matlab data structure
In [ ]:
struct
In [ ]:
box = struct['box'] # Extracts data object from structure

Exercise: What are the dimensions of the imported array?

Vector and array manipulation

Construct sequence of n (linearly-spaced) numbers, from a to b

In [ ]:
a = 50; b = 100-0.1; n = 57;
In [ ]:
sp.linspace(a,b,n)

Construct sequence of n (base-10 logarithmically-spaced) numbers between $10^a$ and $10^b$

In [ ]:
a = -1; b = 1; n = 20;
In [ ]:
sp.logspace(a,b,n)

Construct coordinate array

In [ ]:
x,y = np.mgrid[0:5,0:5]
In [ ]:
x
In [ ]:
y
In [ ]:
np.sqrt(x**2 + y**2)

Construct tiled array

In [ ]:
x = np.linspace(0,10,11);
In [ ]:
x
In [ ]:
np.c_[x,x**2]

Symbolic computation with SymPy

Import SymPy (CAUTION!)

In [ ]:
from sympy import *
In [ ]:
1/2 + 1/3
In [ ]:
Rational(1,2) + Rational(1,3)
In [ ]:
Rational(5,6).evalf(6)
In [ ]:
1./2. + 1./3.

Calculus with symbolic variables

In [ ]:
x = Symbol("x")
In [ ]:
limit( sin(x) / x, x, 0)
In [ ]:
diff(erf(x),x)
In [ ]:
integrate(1./sqrt(pi) * exp(-x**2),x)
In [ ]:
integrate(1./sqrt(pi) * exp(-x**2), (x,-oo,0) ).evalf(3)

Matrix compuation

In [ ]:
M = Matrix( ([1, x], [x, 1]) )
pprint(M)
In [ ]:
%load_ext sympyprinting
In [ ]:
M.inv()
In [ ]:
M.cholesky()

Interfacing with other languages

See files in f2py_example. Things to note:

i) differences between subroutines and functions in Fortran are not really important here

ii) passing arrays is fine, but you need to pass their dimensions as well

iii) In general, SciPy contains many useful functions that work as quickly as Fortran/C. But if you need to crunch something inside a loop, or nested loops, a considerable speed-up is possible

In [ ]:
more example.f90
In [ ]:
!f2py --fcompiler=gfortran -c example.f90 -m example
In [ ]:
import example
In [ ]:
a = 14; b = 78;
In [ ]:
example.arithmetic(a, b)
In [ ]:
c, d = example.arithmetic(a, b)

A weave example (courtesy of Nat Butler) is below

In [ ]:
from scipy import weave
In [ ]:
x = np.arange(1000)/1000.
y = np.zeros(1000,dtype=np.double)
n=len(x)
In [ ]:
code = """
    int i;
    for (i=0;i<n;i++) {
      if (*x<0.5) *y = exp(-(*x)*2);
      else *y = exp(-(*x));
      x++; y++;
    }
"""
In [ ]:
weave.inline(code,['n','x','y'])
In [ ]:
print y
In [ ]:
plot(x,y)

2. SciPy Packages

I. Integration and interpolation

The SciPy interpolation routines are available through the scipy.interpolate module.

In [ ]:
import scipy.interpolate
In [ ]:
scipy.interpolate.
In [ ]:
from scipy.interpolate import interp1d

Create (x,y) points following the function y = x**2

In [ ]:
x = linspace(0,5,10)
y = x**2

Interpolate routines return a function that accepts input at new abscissa (x) locations

In [ ]:
yfn = interp1d(x,y, kind='cubic')
#yfn = interp1d(x,y, kind='cubic', bounds_error = False, fill_value=0)

Plot the interpolated function along a denser line

In [ ]:
xs = linspace(0,5,100)
ys = yfn(xs)
plot(x,y,'k.',xs,ys,'k--')

The scipy.intergrate module provides numerical integration routines, i.e., for summing functions across a 1-, 2- or 3-dimensional domain.

In [ ]:
from scipy import integrate

The following lambda function returns the value $F(x_0) = \int_0^{x_0} y(x)\ dx$ for the function $y(x)$ we defined above with interpolate.

In [ ]:
inty = lambda x0: integrate.quad(yfn,0,x0)
#inty = lambda x0: integrate.quad(yfn,0,x0)[0]
In [ ]:
inty(3)[0]

Alternatively, we could have used trapezoidal integration over the original points, with no iterpolations. Would this make a difference?

In [ ]:
scipy.integrate.trapz?

II. Random sampling and FFTs

In [ ]:
from numpy.random import randn
In [ ]:
print randn()
In [ ]:
y = randn(25)
print y
In [ ]:
from scipy.fftpack import fft
In [ ]:
fy = fft(y)
print fy
In [ ]:
semilogy(range(len(y)),fy * conj(fy))
In [ ]:
sp.fftpack?

III. Optimization and linear algebra

SciPy provides numerical optimization methods that can tackle a broad range of fitting problems. The simplest generic optimization routine in SciPy 0.10 is fmin (which uses Simplex optimization).

In [ ]:
from scipy.optimize import fmin

Optimization routines operate on functions, finding the extreme (minimum) value location in parameter space by stepping through the space in a clever fashion, testing the value of the function at each point and using that information to progress to the next point.

In [ ]:
x = linspace(-5,5,100)
y = lambda x: x**3 - x**2 - x
plot(x,y(x))
In [ ]:
fmin(y, 5)
In [ ]:
fmin(y,-4) # Understand the concept of local minima

Linear algebra is handled in SciPy by the scipy.linalg module. The module wraps Fortran routines from LAPACK into Python.

In [ ]:
from scipy.io import loadmat
data = loadmat('C.mat')
C = data['C']
imshow(C)

Matrix inversion

In [ ]:
from scipy.linalg import inv
imshow(inv(C))

Matrix eigendecomposition

In [ ]:
from scipy.linalg import eigh
E,V = eigh(C)