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.
The
spike_vectorwas 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 literallyto_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 againstspike_positions.npy, though phy already knows the count fromspike_times.npy) and a quality metric (only forfiring_rate = n_spikes / duration).count_num_spikes_per_unit: per-unit counts fall back toto_spike_vector()plusnp.unique, yet they are per-unit sums a unit-indexed backend already knows without building anything.get_last_spike_framematerializes the vector only to take the maxsample_indexof one segment, andget_end_time(get duration) calls it, so computing a sorting's end time reads every spike.get_unit_spike_trainfor 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().sizebut 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 isSortingGenerator, a lazy sorting that never stores aspike_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.spike_vectorcount_total_num_spikescount_num_spikes_per_unitget_last_spike_frameget_unit_spike_train(one unit)Every current path materializes the same ~11 GB vector to return a scalar or a single unit;
get_unit_spike_trainalso 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_vectorrouting assumes the vector is already materialized in memory. That assumption holds for the in-memoryNumpySorting, 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 aspike_vectornatively, and I am thinking about how to improve our sorting representation.