Does spike vector need a full lexsort ?

[samuelgarcia]

After trying to fasten the kilosort/phy reading and subselecting units with sevral way:

• Implement get_unit_spike_trains and performance improvements #4502
• Implement cache spike vector for phy sorting reader #4579

It appears that the main slow down is to force the spike vector to be fully lexsorted:

spikes = spikes[np.lexsort((spikes["unit_index"], spikes["sample_index"], spikes["segment_index"]))]

It use to be partially lexsorted :
spikes = spikes[np.lexsort((spikes["sample_index"], spikes["segment_index"]))]

This commit recent commit 456265c has introduced this "need of fully lexsorted".

In short does 2 spikes with same sample_index need to be ordered by unit_index ?

Yes : we need full lexsort

Pro :

• It ensure always an order of the spiek vector because spike amplitude, amplitude scaling, spike location is a vector related to this spike vector.

Cons :

• it is a big slow down when selecting units
• constructing the spike vector (with many spikes) can be very slow

No : we need partial lexsort

Pro :

• We do not care 2 spikes of same sample_index to be ordered by units for any computation (CCG)
• We only need to ensure that the spike vector is always constructed/read/cached in the same consistent order (for instance by reading the phy/kilosort spike_clusters.npy)

Cons :

• some future part of code could rely on the fact this vector is fully lexsorted

My personnal taste is second case.

The main problem : the spikeinterface codebase use to be No : we need partial lexsort and in february 2025 we switched to Yes : we need full lexsort without too much backward compatibility in mind. Zarr is now opening the internal spike vector and force lexsort in the __init__() which is a big big slow down but also can break the order of spike amplitude, amplitide scaling, ... it this extension was computed before this change (v 0.104.0) and the analyzer was saved in the zarr format before this version.

I would strongly vote for removing the behavior (force full lexsort) and go back to previous situation. But we need to find a better way for backward compatibility to avoid the same mistake.

@grahamfindlay @alejoe91 @chrishalcrow @yger @tayheau
Any opinion on this ?

PS : the first reader able to find more than 10 typos in this post will have a glass of wine.

[alejoe91]

PS : the first reader able to find more than 10 typos in this post will have a glass of wine.

Ahahahaha too easy!

[tayheau]

ill go with the second one if most of the cases relieves on the partial one, but silly question, dont we have an helper method to query the spike vector ? If so maybe we can also add a full_lexsort boolean flag to it with default to False ?

[chrishalcrow]

I'm ok with it being partially lexsorted, I think. We should double check that none of the optimized metrics computations rely on it being fully lexsorted.

@alejoe91 or @yger did you originally introduce the extra layer of lexsorting? What was the motivation? I can understand wanting to have a unique reproducible sorting for any spike vector. Anything else?

[tayheau]

@alejoe91 or @yger did you originally introduce the extra layer of lexsorting? What was the motivation? I can understand wanting to have a unique reproducible sorting for any spike vector. Anything else?

I introduced it in the cited commit with goal in mind to make check_sortings_equal test use the spikevector. That said test needed a full lexsort and we assumed that it could be the case for future methods so we kept it.

[alejoe91]

Partial is ok for me. I think we introduced #4315, and probvably the goal was to make things consistent. But as long as we are clear about it, I think that partial lexsort is ok.

For zarr, we can simply add an _handle_backward_compatibility function for the versions it's affected (should only be 0.104.*).

[grahamfindlay]

It sounds like the real, for-sure requirements for the spike vector are:

1. Deterministic, stable construction, for check_sortings_equal.
2. Partial ordering by ascending segment_index, then by ascending sample_index.

I guess we should check every metric to see if it requires fully lexsorted data (or ask an agent to check), but maybe it is simpler to ask: why would we need full lexsorting? If we need a deterministic "winner" unit in the case of co-temporal spikes, that can be satisfied by any deterministic construction -- not just an ascending sort.

Does the deterministic unit ordering of co-temporal spikes need to be the same for different sorters? For example, would someone want to call check_sortings_equal to compare e.g. a Phy/Kilosort sorting and a DARTsort? If so, then there is an argument for having some canonicalized form of the spike vector. But also, I don't see why it would need to be computed up-front, every time, on load. It could just be done on-demand. We could keep the optimizations around for this.

My 2 cents: get rid of the full lexsorting, pending an exhaustive review of current spike vector consumers.

[grahamfindlay]

I asked Mr. Claude xHigh to review the code and check whether any consumer of the spike vector relies on a unit_index tie-break and, if so, whether the reliance is on an ascending sort vs any old deterministic ordering. Here's what it came back with:

Verdict: none rely on the unit_index tie-break.

 Computation: Synchrony _get_synchrony_counts
 File:line: metrics/quality/misc_metrics.py:1946,1965
 Depends on: np.unique(sample_index) + contiguous slice — partial contiguity only; collects the set of units, tie order irrelevant
 Needs tie-break?: No
 ────────────────────────────────────────
 Computation: CCG / ACG (numba + numpy)
 File:line: postprocessing/correlograms.py:704,776,531,973
 Depends on: pairwise / per-unit; symmetric counting
 Needs tie-break?: No
 ────────────────────────────────────────
 Computation: ISI hist (numba + numpy)
 File:line: postprocessing/isi.py:208,132
 Depends on: per-unit filter spike_clusters==i, ascending samples (partial)
 Needs tie-break?: No
 ────────────────────────────────────────
 Computation: sliding_rp_violation
 File:line: metrics/quality/misc_metrics.py:607
 Depends on: explicit to_reordered_spike_vector(["sample_index","segment_index","unit_index"])
 Needs tie-break?: No
 ────────────────────────────────────────
 Computation: amplitude_cv, drift
 File:line: misc_metrics.py:917,1291
 Depends on: per-unit get_unit_spike_train + searchsorted
 Needs tie-break?: No
 ────────────────────────────────────────
 Computation: amplitude_scalings collision
 File:line: postprocessing/amplitude_scalings.py:416
 Depends on: searchsorted window on sample_index (partial); self found by exact match
 Needs tie-break?: No
 ────────────────────────────────────────
 Computation: PCA
 File:line: postprocessing/principal_component.py:700
 Depends on: positional alignment to random_spikes order
 Needs tie-break?: No
 ────────────────────────────────────────
 Computation: node pipeline / SpikeRetriever
 File:line: core/node_pipeline.py:213,226,275
 Depends on: searchsorted on segment_index/sample_index (partial) + positional copy
 Needs tie-break?: No
 ────────────────────────────────────────
 Computation: check_sortings_equal
 File:line: core/testing.py:115
 Depends on: assert_array_equal on whole vector
 Needs tie-break?: Yes

I checked on the sliding_rp_violation. What it really needs is a unit-major layout where each unit's spike train is contiguous and sample-ascending, so it can slice out spikes[u0:u1] per unit. I actually think my optimization for to_reordered_spike_vector() is going to work just fine for this even if the input isn't fully lexsorted. The stable counting sort over (unit, segment) buckets should (I'd test) give exactly the layout needed, because it only needs the input to be partially sorted (segment, then sample). Each bucket holds exactly one (unit, segment) pair, so there are no cross-unit ties inside a bucket.

So, seems safe to drop the full lexsort guarantee. And some of the optimizations are even still useful 😃 .

[samuelgarcia]

My plan was to make the lextsort before assert_array_equal() only in check_sortings_equal().
In my mind 2 sorting with differents lexsort are the same no ?

[grahamfindlay]

If check_sortings_equal() is only used to compare sortings made with the same sorter, and the spike vector construction is deterministic (but not fully lexsorted), then there is no need to make any lexsort, partial or full, before assert_array_equal(), right?

[samuelgarcia]

If check_sortings_equal() is only used to compare sortings made with the same sorter, and the spike vector construction is deterministic (but not fully lexsorted), then there is no need to make any lexsort, partial or full, before assert_array_equal(), right?

Not sure. Because the spike_vector can be taken from different formats, that can have diffrent spike vector order.

For me the deterministic of ordering is for the same format and file source.
If we save a kilosort to any other format then the spike vector will not be garanty to have the order with this partial lexsort.

[grahamfindlay]

Ah, then yes, I agree that full lexsorting to canonicalize the vector is needed for check_sortings_equal(). 👍