nwbwidgets icon indicating copy to clipboard operation
nwbwidgets copied to clipboard

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

Open felixp8 opened this issue 1 year ago • 2 comments

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 avatar Oct 02 '23 14:10 felixp8

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

felixp8 avatar Oct 02 '23 14:10 felixp8

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

felixp8 avatar Oct 03 '23 11:10 felixp8