Exercises: Scientific Programming with Python

Exercises for the second day of the FacetsPythonCourse2008.

Author: Eilif
Tutors: Eilif, Jens, Jochen

Ex 1. Renewal processes

  1. Generate Poisson process event times with a rate of 10 Hz and length of 10 seconds using numpy.random.exponential

  2. plot an inter-spike interval histogram using both numpy.histogram and pylab.hist

    Hints:

    • Histograms are better plotted with: linestyle, ls = 'steps'

  3. Generate gamma renewal process event times with a rate of 10 Hz, a length of 10 seconds and Coefficient of Variation (CV) of 0.5,1.0,1.5 using scipy.stats.gamma

    Hints:

  4. Plot the peri-stimulus time histogram and ISI distribution of a sinusoidal modulated poisson process

    Hint:

    • discretize time, the poisson process has a mean number of events of dt*rate(t) and is poisson distributed

    • at each time step draw a poisson random number using the command numpy.random.poisson(dt*rate(t)):

      • if it is 0, then no spike occurred at time t
      • if it is 1, then a spike occured.
      • if it is 2, then dt is too large!

    • This approach is somewhat slow.

    • More efficient approaches, such as time warping, or thinning, are given in chapter 6 of this book: http://cg.scs.carleton.ca/~luc/rnbookindex.html

Ex 2. scipy.weave to write an OU process generator

  1. draw gaussian white noise using numpy.random.normal
  2. write C code to solve the discrete form of the OU stochastic differential equation using euler's method
  3. verify mean and std
  4. Compute the auto-correlation of the OU process
  5. Write an OU process generator using scipy.signal.convolve and compare performance using time.time() Hint: it's an example in the matplotlib user quide.

TakeItToTheNextLevel(tm):

  1. modulation by arbitrary functions
  2. doubly stochastic process
    • record your voice.
    • load it using scipy.io.wavfile
    • use voice as modulation of OU process mean

Ex 3. 3D plotting

Plot this function

def f(x, y):
    from numpy import sqrt,cos
    return scipy.cos(sqrt(x**2+y**2)*2)

in 3D with MayaVi!

Hint:

load python as follows:

$ ipython -wthread

for wxWindows treads (the GUI toolkit used by MayaVI)

Resources:

http://neuralensemble.org/cookbook/wiki/InstallMayaVI2

http://scipy.org/Cookbook/MayaVi/mlab

http://scipy.org/Cookbook/MayaVi

http://code.enthought.com/downloads/

easy_install -f http://code.enthought.com/enstaller/eggs/source \
  --prefix=/home/emuller/opt/facets ets

or choose an egg for your system

Ex 4. Curve fitting

Following the example from the lecture

from scipy.optimize import leastsq

# generate data
from numpy import random,histogram,arange,sqrt,exp,nonzero

n = 1000; isi = random.exponential(0.1,size=n)

db = 0.01; bins = arange(0,1.0,db)
h = histogram(isi,bins)[0]

p = h.astype(float)/n/db

# function to be fit
# x - independent variable
# p - tuple of parameters
fitfunc = lambda p, x: exp(-x/p[0])/p[0]
# standard form, here err is absolute error
errfunc = lambda p, x, y, err: (y - fitfunc(p, x)) / err


# initial values for fit parameters
pinit = array([0.2])

# hist count less than 4 has poor estimate of the weight
# don't use in the fitting process
idx = nonzero(h>4)

out = leastsq(errfunc, pinit,\
              args=(bins[idx]+0.01/2, p[idx], \
              p[idx]/sqrt(h[idx])),
              full_output = 1)


l1 = 'data'

errorbar(bins[idx],p[idx],\
         yerr=p[idx]/sqrt(h[idx]),fmt='ko',label=l1)

l2 = r"""$\frac{1}{\lambda} exp(-x/\lambda)$,
$\lambda = %.3f \pm %.3f$""" % (out[0], sqrt(out[1]))

plot(bins,fitfunc((out[0],),bins),'r--',lw=2,label=l2)

legend()
  1. Generate and fit gamma distributed ISI events. Be careful to properly weight the data points to yield std err on fit parameters
  2. Write a loop to repeat the fit 1000 times for different realizations of the gamma ISIs. Keep the resulting fit parameters and check that they are indeed distributed with a variance approximately predicted by the leastsq optimizer output.

Ex 5. Interpolation

# Image you have the following sparse data
x = arange(0,10)
y = sin(x)
# and you want intermediate points ... interpolation

import scipy.interpolate as interp

interp ?
interp.interp1d ?
i1d = interp.interp1d

f_l = i1d(x,y)
f_c = i1d(x,y,'cubic')

figure()
plot(x,y,'k+',markersize=16.0)
x_refined = arange(0,9,0.1)
plot(x_refined,f_l(x_refined),'b-',linewidth=2)
plot(x_refined,f_c(x_refined),'r--',lw=2)
legend(('data','linear-interp','cubic-interp'))
xlabel(r'$\tau_{\theta}$ [$\mu$S]',fontsize=20)
xticks(size=16)
yticks(size=16)

Ex 6. Spikes from voltage traces

  1. Get neuron data from: http://lcn.epfl.ch ... Modeling Challenge A
  2. Write a weave program to find the spike times of the voltage data.

Answer:

def v_to_spikes(v,threshold=0.0):
  """ Get events which exceed spike_threshold in v"""

  spikes = []

  code = """

  const double th = threshold;
  int i;

  for(i=0;i<Nv[0];) {
    if (v(i)>th) {
      spikes.append(i);
      // proceed to where spike stops
      // and continue search there
      i++;
      while ((v(i)>th) && (i<Nv[0])) i++;
    }
    else i++;

  }

  """

  scipy.weave.inline(code,['v','threshold','spikes'],type_converters=scipy.weave.converters.blitz)

  return numpy.array(spikes)

Ex 7. ODE solver

  1. Solve and plot the system:

    dx/dt = x - y
dy/dt = -2*y*(y-1)

    for various initial conditions.

  2. plot the vector field