-
Notifications
You must be signed in to change notification settings - Fork 23
Description
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 toa[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
- I agree to follow this project's Code of Conduct
- Have you ensured this bug was not already reported?