[MRG] ENH: Function to get spike rates

[jasmainak]

closes #154

not sure if this is the right way to do it, just opening a "work in progress thing" and we can discuss the details. How is the spike rate typically calculated?

TODO

• tests
• example
• what's new page

[jasmainak]

@ntolley I think you should take over this PR after the other one is merged. You've been working the spikes object recently, so this is right up your alley. Just add my fork as a remote on git and start pushing to the branch. Let me know if you need help.

[ntolley]

That sounds great to me! Perhaps not for this initial PR, but something to consider in the future is adding a convolution kernel to calculate spike rates. Here is a great example of how that sort function would work.

[jasmainak]

As long as you don't create a dependency on Neo, feel free to do it in this PR :-) Just copy the relevant code (they are BSD-3) and keep it simple. You can have an argument called 'method'

[ntolley]

Reading the issue more closely it actually seems like convolved spike rates gets at a different purpose.

Estimation of instantaneous rates is usually done at the single unit level, whereas @stephanie-r-jones proposes aggregating across cell types to describe the whole trial.

In my mind these are rather distinct and may be better split into two functions. mean_cell_rate() could return mean firing rate (of each trial separately?) across cells. Then spike_rate() could implement instantaneous rate estimation for each cell separately, returning an array of convolved rates indexed by gid. Still not obvious to me if data from each trial should be kept separate.

[jasmainak]

humm ... do you think the mean spiking rates (across trials and cells) could be added to this plot: https://jonescompneurolab.github.io/hnn-core/generated/hnn_core.viz.plot_spikes_raster.html#hnn_core.viz.plot_spikes_raster instead of a function since it's just 4 scalar values?

and we can also have another function for instantaneous rate, but I'm wondering if that should be reserved for another PR? What's not totally obvious to me in that case is:
a) what's the use case for this as opposed to the just mean spiking rates
b) what would the function signature look like and what would be the data structure of what the function returns? I'm assuming you'd need a parameter to control the width of a convolution kernel? And the returned data would be an array of shape n_cells x n_trials ?

[ntolley]

That could possibly work, but it doesn't lend itself easily to accessing the data. This may be a good opportunity to create an analysis.py module to give us more flexibility with feature requests in the future. If it seems important enough to create a member function on the Spikes class, we can just call analysis.py in the same way viz.py is implemented.

wrt plotting I think a box-plot makes the most sense for visualizing the rates across cells so perhaps we can append a subplot to the raster.

For the use case I personally intend to analyze single unit activity in my own projects. Generally it will be extremely useful to have available for highly mechanistic hypotheses. And yes the function would generally look like spike_rates(spikes, kernel='exponential', width=0.1, gids='all'). A matrix return could work, it could also return a dictionary indexed by gid. In any case yes I think this is sufficiently distinct and large enough for a separate PR (I'll open one soon).

[jasmainak]

Let's worry about analysis module later :) I tend to move things into a module once there are enough features to fit in there. If you plan for "future", it will lead to an overengineered solution that is hard to maintain.

So, if I understand correctly are you leaning towards:

1. function for computing instantaneous spike rate?
2. a box plot to go with it

If the user wanted average spike rates across cell types, how would they get it? would they average the return value from that function? or do you want another function for that? or should they just look at the plot which will have the values?

[ntolley]

Honestly I think keeping these separate is the best path. For this PR I'm leaning towards a function to return mean spike rate grouped by cell. Example return would be {'L5_basket: [trial1_mean, trial2_mean], ...}. Since the issue raised specifically refers to using this as an optimization heuristic it should be easy to access the data.

To be clear, for this PR there would be 2 additions to the Spikes class: mean_rate() and plot_mean_rate(). Perhaps eliminate "spike_rate" naming convention since it lives in the Spikes class?

After that I can open a separate PR to introduce instantaneous_rate().

[jasmainak]

okay, let's start coding :) We can iterate the API once we have something down. Can you push to this branch?

[ntolley]

Sorry it looks like this branch was behind your master so the rebase included the latest PR's merged. I'll work on cleaning it up tomorrow morning.

In either case all my work is on the last commit. It now computes the mean spike rates for every individual gid, you can then specify what type of average you want returned by setting mean_type='all', 'trial', or 'cell'.

[jasmainak]

Good luck with the rebase. Maybe the force be with you.

You can go ahead and add the tests as well.

[ntolley]

I've taken a chainsaw to the commit tree and it's looking beautiful again.

Regarding the tests, obvious ones are checking for correct mean_type calls. Since this relies on data members of the Spikes class I believe the existing tests will catch any problems that could arise in the rate calculations. Can you think of anything that should be caught?

[jasmainak]

For tests, you should create a Spikes object where you know what the spike rate is and then check that the function returns the same spike rate.

[jasmainak]

I've taken a chainsaw to the commit tree and it's looking beautiful again.

amazing!! :)

[ntolley]

Surprisingly beefy commit. I've implemented the tstart and tstop attributes. Let me know if you disagree, but it seemed like the most logical place to instantiate the attributes was _gather_trial_data() so all the data is stored from the same function.

Also there is now an issue with legacy spk_*.txt files and storing the new attributes; the write methods now include tstart and tstop. I've added additional logic to check for these additional columns.

TODO:

• Update examples
• Add test for precomputed rates

[jasmainak]

It's converging @ntolley slowly but surely :)

[codecov-commenter]

Codecov Report

Merging #155 into master will decrease coverage by 0.57%.
The diff coverage is 36.84%.

@@            Coverage Diff             @@
##           master     #155      +/-   ##
==========================================
- Coverage   67.80%   67.23%   -0.58%     
==========================================
  Files          19       19              
  Lines        2022     2060      +38     
==========================================
+ Hits         1371     1385      +14     
- Misses        651      675      +24     

Impacted Files	Coverage Δ
hnn_core/network.py	53.71% <4.00%> (-8.29%)	⬇️
hnn_core/tests/test_network.py	100.00% <100.00%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update a1cc05d...f24257d. Read the comment docs.

[ntolley]

Definitely procrastinated on rebasing for too long. In line with the previous discussion I removed all of the code related to storing tstart and tstop as attributes and writing to a file. Instead they are passed to mean_rates() as parameters along with gid_dict. This additionally fixes the case where gids that never spike across trials aren't included in the final calculation.

This heavily simplifies the code and cleans up the awkward logical checks.

[jasmainak]

@ntolley looks pretty clean to me and I see you have added tests too! Fantastic.

Would you mind addressing the last few comments and then we can merge. Don't forget to change the title from WIP to MRG once you are ready

[jasmainak]

Oh, one last thing is that can you update an example to show the use of this new function?

[ntolley]

That should be everything! Leaving a few conversations unresolved as I think some of the comments will be useful to reference in future PR's.

[jasmainak]

It's in, thanks @ntolley ! 🍻 🧠