Running MATLAB(R) scripts in pure Python

The demonstration described in this article shows a partial pure-Python implementation of MATLAB(R)'s basic numerical objects that are sufficient to run the "Simple Model of Spiking Neurons" (http://vesicle.nsi.edu/users/izhikevich/publications/spikes.htm) of Izhikevich(2003).

The results of running the MATLAB(R) and the OMPC version of the Izhikevich model are on the following figure:
 
MATLAB(R)
 
OMPC (plot with Matplotlib)
You can download the code along with other OMPC examples in a zip file. No installation is required, all but one example in the archive does not require any additional modules and will work on a pure Python installation.


Requirements

First you will need a Python interpreter. The development and most of the testing is done with Python 2.5 (get it here). 
For the more adventurous readers, you could try any other Python as well, I did. Particularly the script described in this article was tested with the still fresh and warm Python 2.6. Surprisingly to me Jython2.5.a1 that I happened to have installed for some time worked as well. I didn't have so much luck with IronPython 2.0, but I gave it just one try. Python 3.0 starts running but hangs somewhere in the places where generators are used heavily.

The second requirement is of course the source file. The source m-file is also available in the examples zip file.
In the next step you need to translate the m-file using OMPC. There is actually a number of ways. Choose one of them:
  1. use the following command from the root of the unzipped Examples archive: ompc/ompcply.py mfiles/izhikevich.m > izhikevich.pym
  2. use the online compiler at http://ompclib.appspot.com/ 
  3. use the translate.py script provided with the examples, this scripts demonstrate how you can use the online compiler from Python without having the OMPC installed.
After a successful translation you have have a script similar to this:
# Created by Eugene M. Izhikevich, February 25, 2003
# Excitatory neurons   Inhibitory neurons
Ne = 800
Ni = 200

re = rand(Ne, 1)
ri = rand(Ni, 1)

a = mcat([0.02 * ones(Ne, 1), OMPCSEMI, 0.02 + 0.08 * ri])
b = mcat([0.2 * ones(Ne, 1), OMPCSEMI, 0.25 - 0.05 * ri])
c = mcat([-65 + 15 * re **elpow** 2, OMPCSEMI, -65 * ones(Ni, 1)])
d = mcat([8 - 6 * re **elpow** 2, OMPCSEMI, 2 * ones(Ni, 1)])
S = mcat([0.5 * rand(Ne + Ni, Ne), -rand(Ne + Ni, Ni)])

v = -65 * ones(Ne + Ni, 1)# Initial values of v
u = b *elmul* v# Initial values of u
firings = mcat([])# spike timings

for t in mslice[1:1000]:# simulation of 1000 ms
    I = mcat([5 * randn(Ne, 1), OMPCSEMI, 2 * randn(Ni, 1)])# thalamic input
    fired = find(v >= 30)# indices of spikes
    if not isempty(fired):
        firings = mcat([firings, OMPCSEMI, t + 0 * fired, fired])
        v(fired).lvalue = c(fired)
        u(fired).lvalue = u(fired) + d(fired)
        I = I + sum(S(mslice[:], fired), 2)
    end
    v = v + 0.5 * (0.04 * v **elpow** 2 + 5 * v + 140 - u + I)
    v = v + 0.5 * (0.04 * v **elpow** 2 + 5 * v + 140 - u + I)
    u = u + a *elmul* (b *elmul* v - u)
end
plot(firings(mslice[:], 1), firings(mslice[:], 2), mstring('.'))

Running it

At this point it is better to move into the folder izhikevich/. This folder contains all the files necessary to tun this demo. I suggest putting a print t statement at the beginning of the for loop, otherwise you might think that the interpreter is dead. Once there run this command
    python -m ompc_pure izhikevich.pym
To look at a version that requires numpy and matplotlib:
    python -m ompc_numpy izhikevich.pym

The numpy version runs slower that a translation by hand, but there are many ways to improve the script. Take it as the first prototype.


It is slow

The pure Python version is slow because of the many assignment and read operations on the arrays. The script uses the very slow _ndilin and ndilin1 function. This is the next thing on my schedule to correct. To read more about the bottle neck read the section called Optimizing element access at the end of this article.
The numpy version shows that things can be made reasonably fast. It is not as fast as the numpy version translated by hand (izhikevich/izhikevich.py) but it is a good start.
To show how a fast the script can get in Python I have include a translation by hand into Python+numpy. This file is called izhikevich/izhikevich.py. On my system this version runs in 2 seconds. Compared to MATLAB(R)'s incredible 0.5s (but with the interpreter initialized,  which took 20 seconds). I think given the overhead of __?etitem1__ and elementwise operations it should be possible to get somewhere between 2 and 3 seconds, but let's see.

Why is it working

To execute this script somebody needs to implement a number of functions and a numerical array object that behaves similar to the MATLAB(R)'s mxArray. I did it. I call the numerical object marray. A number of additional objects are ready to be used in the file izhikevich/ompclib_demo.py, namely mslice, mstring and their base object mvar, a helper for element-wise operations elmul and elpow, and a special object end for slicing.
The marray object support assignment through its overloaded __call__ method (actually mvar.__call__). This function accepts a variable number of arguments and passes them to objects __setitem1__. This function has the same interface __setitem__ but assumes that indices passed to it are based on 1 (the 1st element of the array is 1 as in MATLAB(R)). Python doesn't allow slices in argument list, that is why OMPC generated the last line with firings(mslice[:],1).  For the retrieval of element funcntion __getitem1__ is called. 
Actually if you look at a code a(1,1).lvalue = 1; you might realize that there is no way to tell that __call__ operator is followed by an assignment. Actually for each assignment this happens:
  1. a.__getitem1__(1,1) is called
  2. __getitem1__ returns a new object, let's call it temp, that keeps reference to the original object it was retrieved from. This is neccessary to get code like x = a(1,1) working as well.
  3. temp.lvalue = 1 causes a call to the mvar._lvalue_set function which in turn calls the a.__setitem1__((1,1),1) of the parent object.
  4. The result of the __getitem1__ operation is discarded by the garbage collector.

The mslice object is a generator that can turn into a real marray when necessary. For example when you ask for the ctypes buffer, data have to be generated. For slicing of the array however the object can be passed as a Python generator. In the case when the end object is used, for example mslice[1:end-1], this object is evaluated at the time it is used for slicing. This is because this is the only lace where the shape of the slices object is know to the mslice.

The element-wise operations are mediated by a special object. The elmul object, for example, in the expression a*elmul*b, will execute as follows:
  1. a*elmul is called, elmul remembers a -> elmul.left = a
  2. elmul*b is called, elmul realizes that is has a left operand ready a tries calling a.__elmul__(b), if it fails, calls b.__elmul__(a).
The function mcat concatenates the contents of its single list argument in the way MATLAB(R) concatenates arrays in pair of brackets. The OMPCSEMI is a special symbol the emulates the function of ';' in MATLAB(R).

Functions rand, find, ones and sum have a straightforward implementations (in the file ompclib_demo.py). They all return marrays to provide a close system for the script, otherwise all arguments passed to any functions have to be checked for their type and converted into marrays.

Optimizing element access

After running this example you will realize that it is (painfully) slow. Actually all this is cause by long loops and use of generators. Generators are roughly 2x slower than real loops in Python 2.5, but they are somehow much easier and cleaner to work with.
After running a profiler you will realize that the bottleneck it the function _ndilin and _ndilin1. These function take N slices and a shape of an N dimensional array and calculate the linear index into a FORTRAN ordered array that is at the core of the marray object.
For example the statement S = mcat([0.5 * rand(Ne + Ni, Ne), -rand(Ne + Ni, Ni)]) requires an operation S[1:1000,1:800] = 0.5 * rand(Ne + Ni, Ne), that copies 800,000 elements. On my system with Python 2.5 _ndilin takes good 9s to generate all the indices. I guess the generators are a thing of the future I tried a similar function in other Pythons, IronPython 2.0 won with 2.0s and Python 3.0 did it just a little bit slower. The jython2.5a1 that I happened  to have was crunching for almost 12s. 
A C++ implementation of these 2 functions is ready in the file ompc_fast.cpp and needs to be incorporated. For orientation, on my system the 800,000 indices take 0.45s.I know that there is a much simple algorithm, but I simply love the idea of generators.
But more ... another day in another article.

References

Izhikevich, E. M. (2003). "Simple Model of Spiking Neurons", IEEE Transactions on Neural Networks, Vol. 14, No. 6, 1569-1572


Comments