Spike train metrics in utils

[arvoelke]

Looking for some early feedback on this one.

I needed the Victor-Purpura distance metric for a paper I'm working on, in order to find some neurons that "match" spike recordings (Bensmaia (2010) used this metric for the same model).

I wasn't sure how our code licensing works in these cases, so I made sure to not look at any code (even pseudo-code) to implement this algorithm. I just used the description from the above link and convinced myself that it's doing the right thing. If someone is curious enough to cross-check it that might be nice.

Are there any other metrics that we think would be good to include here?

[jgosmann]

As it happens, I implemented a number of spike train metrics in Python two years ago as part of spykeutils. The code is BSD licensed in case we want to include it in Nengo, but I am not sure whether we really want to maintain spike train metrics in our code base.

There is also another package, Elephant, in the making which is supposed to include those spike train metrics with implementations based on the aforementioned implementation.

In case you decide to look at other code, the commonly used efficient algorithm for the Victor-Purpura distance is described in
Aronov, Dmitriy. “Fast Algorithm for the Metric-Space Analysis of Simultaneous Responses of Multiple Single Neurons.” Journal of Neuroscience Methods 124, no. 2 (April 2003): 175–79. doi:10.1016/S0165-0270(03)00006-2.

It is worth noting that I spend a significant amount of time optimizing the implementation of the Victor-Purpura distance so that it uses mostly vectorized Numpy operations (keeping only one for loop). In theory this also gives a reduction of memory usage from quadratic to linear (compared to the paper mentioned), but in practice it is faster to precompute the quadratic matrix.

[hunse]

Haha, spykeutils is a great name.

I am not sure whether we really want to maintain spike train metrics in our code base.

In general, I would say we don't. Not only would this add extra weight and re-invent the wheel, but it would be nice for people to be able to use the same analysis code whether spikes come from a nengo model, a petri dish, or a live animal. I could see us outputting to some standard spike storage format to facilitate this.

That said, I would be okay with us having a handful of basic analysis functions if there are some that are pretty standard, can be implemented fairly easily, and would make it nicer to do quick analyses in nengo. I know nothing about spike train analysis, so I can't suggest what those might be.

That said, my only comment on this PR is that I would call the file spike_metrics.py instead of metrics.py, since there are many more metrics than those for spike train analysis.

[jgosmann]

I could see us outputting to some standard spike storage format to facilitate this.

Neo would be a good candidate for that.

I would be okay with us having a handful of basic analysis functions if there are some that are pretty standard

Spike train metrics don't seem very basic nor very standard to me.

[arvoelke]

I've been using spykeutils with neo and it's a lot better than what I have here, mainly because my code is too inefficient for sparse trains.

I do still think that stuff like this would be a valuable addition to Nengo, since it provides a more well-rounded toolkit with little maintenance cost and no interdependencies. But this exercise has shown that implementing the algorithm well is not quite so simple.

For now we should direct people to spykeutils with neo if they need to do this sort of analysis in Nengo. Here's some example usage for future reference:

import quantities as pq
from neo.core import SpikeTrain
from spykeutils.spike_train_metrics import victor_purpura_dist

...
p_spikes = nengo.Probe(model.x.neurons, synapse=None)
...

simt = sim.trange()
t_stop = simt[-1] * quantities.s
trains = [SpikeTrain(simt[train > 0], units=pq.s, t_stop=t_stop)
          for train in sim.data[p_spikes].T]
dists = victor_purpura_dist(trains, q=pq.Hz)

[jgosmann]

It's nice to see that my efforts a few years back were worthwhile. 😃

I could definitely see some Nengo functions to make it easy to get the appropiate neo objects. But I'm not yet convinced about the actual analysis functions.

[tbekolay]

Probing to NeoHDF5 files is definitely on the roadmap! 🌈 🐙