Scientific Python

This tutorial introduces

  • ipython, a powerful CLI for python providing TAB completion, history, helper functions
  • numpy, the main python package to deal with multidimensional arrays, masked arrays, random numbers, linear algebra and fast fourier transform
  • scipy, a collection of mathematical functions and operations on top of numpy
  • matplotlib, a plotting package which provide high-quality output on the screen and for publications. It handles 1d, 2d and 3d data and is accessible through a functional (matlab like) and object oriented interface.

IPython

Start ipython in a plotting event loop, this allows to interactively plot, without python being blocked. 

ipython -pylab
  • Run shell commands using ! 
!firefox
  • Use <TAB> for completion of objects, variables, filenames 
In [3]: x.s<TAB>
x.split       x.splitlines  x.startswith  x.strip  x.swapcase
# get help
x.swapcase?
  • history through up/down arrows <Ctrl-r> (search the history)

The history is written to ~/.ipython/history on exit.

  • macros 
# print the history with numbers
%history
# define a macro from lines 2 to 5 in the history
%macro mymacro 2-5

Numpy - Arrays

from numpy import *
### Numpy arrays
a = array(range(16)).reshape(4,4)
print a
b = array(range(16,32)).reshape(4,4)
print b
msk = logical_and(a>2, a<7)
print msk
amasked = ma.MaskedArray(a, mask=msk)
print amasked
c = amasked+b
print c.data
print c.mask

# auto-ranges
rng = linspace(-1.0,1.0,10)
print rng, len(rng)

Matplotlib - functional interface

Attach:testdata

Read in some data via numpy.loadtxt modify and plot it. 

from scipy.signal import medfilt

# numpy utility function to read column data from textfile
data = loadtxt("testdata.txt", usecols=(0,1))
x = data[:,0] # get first column using ":" to get all elements in that axis
y = data[:,1] # gte second column

plot(x, y)
# x/y label as latex
xlbl = r'$v_{\rm{LSR}}\ (\rm{kms}^{-1})$'
xlabel(xlbl)
ylbl=r"$T_{a}^{*}\ (\rm{Jy})$"
ylabel(ylbl)
# set the viewport
xlim(min(x),max(x))

yinterf = y.copy() # make a copy x=y references y and any change to x changes y

# introduce a spike
yinterf[200] = 300.0

# median filter the data using a kernel of width 5
yfilt = medfilt(yinterf, 5)
cla() # clear plot
# plot several curves plus legend using labels
plot(x, yinterf, label='interference')
plot(x, yfilt, label='filtered', color="red")
xlabel(xlbl)
ylabel(ylbl)
xlim(min(x),max(x))
legend()

Use FFTs 

cla()
plot(abs(yfft))
xlim(1,len(yfft)/2) # set viewport

# flag the 'bad frequencies'
flag0 = 26 # xvalue
dflag = 3 # xoffset
yfft[flag0:flag0+dflag] = 0+0j # use python/numpy's slicing capability
yfft[-(flag0+dflag):-flag0] = 0+0j # and from the end using "-"
xlim(0,100)

# plot thr flagged regions
axvspan(flag0,flag0+dflag, alpha=0.2, facecolor="green")
xfft = arange(len(yfft))
axvspan(xfft[-(flag0+dflag)],xfft[-flag0], alpha=0.2, facecolor="green")

yifft= fft.ifft(yfft).real # only interested in teh real part, use attribute
cla()
plot(y, label='signal')
plot(yifft, label="fft'ed")
plot(y2, label="signal+sin")
plot(y2-yifft, label="residual")
legend()

Fitting 

# Fit a gaussian to profile
import scipy
from scipy.optimize import leastsq

# make copies just in case
xdat = x.copy()
ydat = y.copy()

# the gaussian to fit and its associated residual to minimize
# this uses the lambda keyword an "inline" function
gaussfunc = lambda p, x: p[0]*scipy.exp(-(x-p[1])**2/(2.0*p[2]**2))
residfunc = lambda p, x, y: gaussfunc(p,x)-y

# initial estimates required
# take the peak of the spectrum
peak0 = max(ydat)
# and it's associated velocity
v0 = xdat[scipy.where(ydat==max(ydat))[0]]
# use 50 channel(pixel) width as sigma
sigma0 = xdat[-1]-xdat[0]/len(xdat)*50

#The inital estimates
initial = [peak0, v0, sigma0]

# least-square fit the residuals
fitparams, success = leastsq(residfunc, initial, args=(xdat,ydat))

# plot data, fit and residual
yfit = gaussfunc(fitparams, xdat)
yresid = residfunc(fitparams, xdat, ydat)

cla()
plot(xdat, ydat, label='data')
plot(xdat, yfit, label='fit')
plot(xdat, yresid, label='residual')
legend()
xlim(-10,20)

Summary

scipy contains a lot of packages and functions. Explore it! help(scipy) is your friend.

Example

Attach:ai2009_scipy.py

Presentation

Recommended reading