NeurodataWithoutBorders / nwbwidgets

Explore the hierarchical structure of NWB 2.0 files and visualize data with Jupyter widgets.
https://nwb-widgets.readthedocs.io/en/latest/
Other
47 stars 22 forks source link

[Bug]: Smoothed PSTHs ignore multiple spikes occurring in same bin #313

Open felixp8 opened 11 months ago

felixp8 commented 11 months ago

What happened?

nwbwidgets.analysis.spikes.compute_smoothed_firing_rate currently uses np.searchsorted to find bin indices for spike times, and then indexes into a zero array using those bin indices and increments that by one. When there a bin index is repeated, that bin is only incremented once.

A simple solution is to use np.add.at. From the docs:

ufunc.at(a, indices, b=None, /) Performs unbuffered in place operation on operand ‘a’ for elements specified by ‘indices’. For addition ufunc, this method is equivalent to a[indices] += b, except that results are accumulated for elements that are indexed more than once. For example, a[[0,0]] += 1 will only increment the first element once because of buffering, whereas add.at(a, [0,0], 1) will increment the first element twice.

As a side note, np.searchsorted returns the index at which an element should be inserted to preserve order, e.g. if the timestamps are [0.0, 1.0, 2.0], a spike time of 0.5 would return 1 from np.searchsorted and be placed in bin 1, although conventional histogram/binning practices would place it in bin 0. If this is unintentional, to fix this we could subtract 1 from all of these indices and set side='right' in np.searchsorted instead of the default side='left'.

Steps to Reproduce

from nwbwidgets.analysis.spikes import compute_smoothed_firing_rate

spike_times = np.array([0.5, 2.0, 2.5, 3.0])
tt = np.arange(5)
sigma_in_secs = 0.01 # effectively no smoothing
expected_binned_spikes = np.array([1., 0., 2., 1., 0.])

smoothed_fr = compute_smoothed_firing_rate(spike_times, tt, sigma_in_secs)
assert np.testing.assert_equal(smoothed_fr, expected_binned_spikes)

Traceback

AssertionError: 
Arrays are not equal

Mismatched elements: 3 / 5 (60%)
Max absolute difference: 1.
Max relative difference: 1.
 x: array([0., 1., 1., 1., 0.])
 y: array([1., 0., 2., 1., 0.])

Operating System

Linux

Python Version

3.10

Package Versions

nwbwidgets==0.11.3

Code of Conduct

felixp8 commented 11 months ago

happy to fix it myself, but would like confirmation whether the general binning behavior should be fixed as well

felixp8 commented 11 months ago

also is there a reason why zeros are returned even if there is one spike? based on https://github.com/NeurodataWithoutBorders/nwbwidgets/blob/master/nwbwidgets/analysis/spikes.py#L17