Fix Llobet expected-violations formula in sliding RP metric

[oliche]

Fix Llobet expected-violations formula in the sliding RP metric

_compute_violations computed expected violations as Ve = 2·τ/D · Nc·N (i.e. Nc·(Nb + Nc)), which overestimates Ve and inflates the contamination confidence. Corrected to the Llobet et al. (2022) form used by the reference implementation:

Ve = 2·τ/D · Nc·(Nb + (Nc − 1)/2)     Nc = C·N,  Nb = (1 − C)·N,  D = N / firing_rate

Also extends the default contamination grid to include 35 % (np.arange(0.5, 35.5, 0.5), previously stopped at 34.5 %) and adds a unit test for _compute_violations against the analytical formula.

Reference: https://github.com/SteinmetzLab/slidingRefractory (metrics.computeViol).

Known non-equivalences (out of scope)

The metric definition now matches the reference, but results are not bit-identical because SI (by existing design) uses a coarser ACG bin size (bin_size_ms=0.25 vs sample-resolution) and a different correlogram engine (_compute_correlograms_* vs the reference histdiff-equivalent). Full numerical parity would be a separate change.

[oliche]

The original repo fix is here for reference: SteinmetzLab/slidingRefractory#8

[zm711]

I think @DradeAW wrote the initial implementation based on working with Llobet right Aurélien? If you have a moment to comment on this that would be great.

[chrishalcrow]

Hi @oliche, thanks for pushing your changes to si! I thought I'd try to verify the change in the Llobet formula but I'm getting confused when reading the original paper/calculation. Some more specific referencing would help me: where in the paper is the discussion about the expected violations? And doesn't this depend on the modeling assumptions: e.g. whether the contamination is due to a single contaminating neuron or totally random noise? Thanks!

[DradeAW]

Hi all!
I see you referencing Llobet et al. 2022, but we have a (very) updated version of the pre-print : https://www.biorxiv.org/content/10.1101/2022.02.08.479192v2
Everything is explained in the methods.
In short: the formula is 100% correct while the censored period is 0. But with a censored period, the correction actually depends on the ISI without contamination (which is unknown!). Hence, we computed an approximation which isn’t perfect as perfection cannot be achieved (to our knowledge). Probably a better approximation can be achieved.
I’ll look at your code later this afternoon :)

[DradeAW]

So I looked over some code:
The Llobet formula is something I implemented in SpikeInterface with the functions compute_refrac_period_violations, _compute_rp_contamination_one_unit and _compute_nb_violations_numba.
I looked over the code, and it seems correct.

You are proposing a change in the function _compute_violations, which is used in slidingRP_violations.
I never wrote this function and never used it, so I cannot help you guys on this subject :/

However, if there are questions about the Llobet formula, I would be happy to answer you :)

Best to you guys!