nwbwidgets
nwbwidgets copied to clipboard
[Bug]: Smoothed PSTHs ignore multiple spikes occurring in same bin
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, whereasadd.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
- [X] I agree to follow this project's Code of Conduct
- [X] Have you ensured this bug was not already reported?
happy to fix it myself, but would like confirmation whether the general binning behavior should be fixed as well
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