# 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)

In [ ]:
inty(3)


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)