SpikeInterface / SpikeInterface/spikeinterface

Sorting simple queries should not force full data reading

Open
#4,702 0 comments 1 reaction 0 assignees View on GitHub

Nobody has claimed this yet.

Dominant language
Python
Stars
847
Forks
280
Avg merge
3d 9h
Merged PRs (30d)
29

Description

The spike_vector was designed for one job: a fully-in-memory, time-ordered structured array to iterate for time-sweeping computations (peeling, waveform extraction). It is a good representation for that job. The problem is that it has quietly become the mandatory path for computations that do not need it. A large number of methods materialize the entire merged vector (and sometimes a second unit-grouped copy) even when the answer is a narrow property that a unit-indexed backend already knows without building anything. That is wasteful of memory, and it is actively hostile to streaming, where the whole point is to avoid pulling every spike into RAM.

Concrete chokepoints where a narrow operation forces full materialization today:

  • count_total_num_spikes: the total count is literally to_spike_vector().size, building the whole array only to read its length. The same "materialize everything just to count" appears in the phy exporter (a length-check against spike_positions.npy, though phy already knows the count from spike_times.npy) and a quality metric (only for firing_rate = n_spikes / duration).
  • count_num_spikes_per_unit: per-unit counts fall back to to_spike_vector() plus np.unique, yet they are per-unit sums a unit-indexed backend already knows without building anything.
  • get_last_spike_frame materializes the vector only to take the max sample_index of one segment, and get_end_time (get duration) calls it, so computing a sorting's end time reads every spike.
  • get_unit_spike_train for a single unit warms the full lexsorted cache (the merged vector plus a second unit-grouped copy) just to return one unit's train.

The distinction is whether a count (or narrow property) is the sole need: sites like the principal-components projection call the same to_spike_vector().size but then iterate every spike anyway, so there it is a free by-product, not a chokepoint. For a backend that stores spike trains per unit (NWB's units table) the genuine chokepoints are answerable directly, the count from row lengths, one unit's train as one row read, and even a columnar store like zarr can answer a total count from a dataset shape without reading data.

Benchmark: four narrow methods, each computed two ways and compared on wall time and peak RAM, the current implementation (which builds the merged spike_vector) versus the lazy per-unit path (get_unit_spike_train(..., use_cache=False), materializes nothing). The instrument is SortingGenerator, a lazy sorting that never stores a spike_vector, scaled to a full day of one IBL-rate probe (400 units, 24 h, 4.37 Hz, ~151M spikes; the rate is grounded in DANDI 000409: 867 units, 16.4M spikes, 72 min). Each measurement runs in its own process; the import + construct baseline is 0.05 GB.

method current, builds spike_vector lazy, per-unit
count_total_num_spikes 18.3 s, 11.4 GB 3.3 s, 0.06 GB
count_num_spikes_per_unit 19.0 s, 11.4 GB 3.4 s, 0.06 GB
get_last_spike_frame 17.2 s, 11.4 GB 2.7 s, 0.06 GB
get_unit_spike_train (one unit) 52.6 s, 11.4 GB 0.01 s, 0.06 GB

Every current path materializes the same ~11 GB vector to return a scalar or a single unit; get_unit_spike_train also builds the unit-grouped reorder on top, so one unit's train costs 52 s. The cost is linear in duration, so a week-long recording would need on the order of 55 GB just to return an integer. Long chronic recordings are the motivating case, a workflow I know @chrishalcrow and @alejoe91 care about. Benchmark script: gist.

I'm filing this as the problem rather than a finished proposal, since I'm still working out the fix. My current read is that the root cause is narrow: the spike_vector routing assumes the vector is already materialized in memory. That assumption holds for the in-memory NumpySorting, but for large and streaming sortings it means a scalar query pays a full O(N) materialization for an O(1) answer. I also think it is problematic for sortings that do not hold a spike_vector natively, and I am thinking about how to improve our sorting representation.

Contributor guide

No contributing guide indexed for this repository

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

Research direction

Start with the linked methods in src/spikeinterface/core/basesorting.py: get_unit_spike_train, count_num_spikes_per_unit, count_total_num_spikes, and get_last_spike_frame. Review the phy exporter and quality-metric call sites, then run the linked SortingGenerator benchmark to compare wall time and peak RAM. Done means narrow queries no longer materialize the merged spike_vector when it is their sole need, while retaining their existing results.

Written by the indexing model from the issue text.

Assessment

Tech stack
python
Domain
data, performance
Issue type
Refactor
Difficulty
5/5
Estimated time
Over a week
Activity status
Quiet
Clarity
Needs clarification
Newbie friendliness
35/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.