Add time-ordered and unit-grouped spike access to BaseSorting #4851

I started #4817 so the caches inside BaseSorting say which order of spikes they hold. That makes the order explicit internally but the more I think about it I think we should make it explicit in the public API too.

Right now the public method names on BaseSorting do not make the spike order explicit. Most things go through the time-ordered spike vector of to_spike_vector() or through a reordering of it with to_reordered_spike_vector(lexsort=...).
However, the methods in spikeinterface differ in which order works better for them. Some want the spikes grouped by unit, for example:

• get_unit_spike_train
• the per-unit metrics

And others want them ordered by time, for example:

• waveform extraction
• correlograms

This has two costs:

• Memory: requesting unit-grouped spikes through to_reordered_spike_vector() first builds the time-ordered vector, even when the format already stores spikes grouped by unit (#4702). Even counting the spikes of a sorting builds it:
  

  spikeinterface/src/spikeinterface/core/basesorting.py

  Lines 718 to 730 in 7a4fc1d

  	def count_total_num_spikes(self) -> int:
  	"""
  	Get total number of spikes in the sorting.
  	This is the sum of all spikes in all segments across all units.
  	Returns
  	-------
  	total_num_spikes : int
  	The total number of spike
  	"""
  	return self.to_spike_vector().size

• Confusion: to_reordered_spike_vector does not say which order it returns either, and #4683 is a whole thread on which orders the lexsort parameter should accept and what each one should return.

The proposal is two methods on BaseSorting, time_ordered_spikes(segment_index: int | None) and unit_grouped_spikes(segment_index: int | None), and a computation declares the order it needs by calling one of them. I first suggested them in #4683:
#4683 (comment)

A per-unit computation like the ISI violations wants each unit's spikes together:

spikes = sorting.unit_grouped_spikes()
for unit_index, unit_id in enumerate(sorting.unit_ids):
    start, stop = np.searchsorted(spikes["unit_index"], [unit_index, unit_index + 1])
    train = spikes[start:stop]
    violations[unit_id] = np.sum(np.diff(train["sample_index"]) < refractory_samples)

And population activity in a time window wants every unit's spikes in time order:

spikes = sorting.time_ordered_spikes()
start, stop = np.searchsorted(spikes["sample_index"], [t0, t1])
window = spikes[start:stop]

The order is the guarantee (in the unit-grouped one each unit's spikes are sorted by sample_index). Beyond the order, they support what the spike vector supports: columns (at least sample_index and unit_index) and slicing.

These are new methods and they can start clean without breaking the old ones. As each returns the spikes of a single segment the returned spikes don't need a segment_index field. This will simplify moving forward with #4816. And as they don't promise a numpy array a format can return a lazy view.

The main advantage is that once the order is in the method name, each format can implement each method its own way. The formats we own (zarr and the numpy binary folder) can store both orders internally and serve both methods lazily from disk, instead of reading the whole sorting into memory as we do right now. NWB stores the unit-grouped order today and I am working on storing the time-ordered one in ndx-spikesorting:
catalystneuro/ndx-spikesorting#25

The other formats (neo, phy and the rest) just declare which order they have and the other one is materialized on demand, only when a method asks for it.

to_spike_vector() stays exactly as it is and keeps returning the cached array without a copy. to_reordered_spike_vector() is deprecated toward unit_grouped_spikes() which is way more clear.

@samuelgarcia @alejoe91 @grahamfindlay @chrishalcrow @JoeZiminski what do you think?