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:
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. RequirementsFirst 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:
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 itAt 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 slowThe 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 workingTo 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:
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:
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 accessAfter 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.
ReferencesIzhikevich, E. M. (2003). "Simple Model of Spiking Neurons", IEEE Transactions on Neural Networks, Vol. 14, No. 6, 1569-1572 |
|
|
