[Fix] spike time tiling coefficient for unsorted spiketrains, added validation test

[Moritz-Alexander-Kern]

This is related to Issue #563 .

The example from the Issue now returns:

• STTC for E7 and E3: 1.0 (Spiketrains identical)
• STTC for E8 and B3: -0.02 (Spiketrains different)

Previously, this function encountered a bug when using unsorted spiketrains as input, leading to the generation of values exceeding the desired range of -1 to 1.

This fix involves several improvements:

• the spiketrains used as input are now sorted
• additional input checks were added
• added a regression test, verifying that the sttc function returns values within the desired range of -1 to 1, with the example given in the Issue

Example Code:

import neo

import quantities as pq

from elephant.spike_train_correlation import spike_time_tiling_coefficient

spiketrain_E7 = neo.SpikeTrain([ 1678. , 23786.3, 34641.8, 71520.7, 73606.9, 78383.3,
97387.9, 144313.4, 4607.6, 19275.1, 152894.2, 44240.1], units='ms', t_stop=160000*pq.ms)

spiketrain_E3 = neo.SpikeTrain([ 1678. , 23786.3, 34641.8, 71520.7, 73606.9, 78383.3,
97387.9, 144313.4, 4607.6, 19275.1, 152894.2, 44240.1], units='ms', t_stop=160000*pq.ms)

print(spike_time_tiling_coefficient(spiketrain_E7, spiketrain_E3, dt=0.10 * pq.s))

spiketrain_E8 = neo.SpikeTrain([ 20646.8, 25875.1, 26154.4, 35121. , 55909.7, 79164.8,
110849.8, 117484.1, 3731.5, 4213.9, 119995.1, 123748.1,
171016.8, 172989. , 185145.2, 12043.5, 185995.9, 186740.1,
12629.8, 23394.3, 34993.2], units='ms', t_stop=210000*pq.ms)

spiketrain_B3 = neo.SpikeTrain([ 10600.7, 19699.6, 22803. , 40769.3, 121385.7, 127402.9,
130829.2, 134363.8, 1193.5, 8012.7, 142037.3, 146628.2,
165925.3, 168489.3, 175194.3, 10339.8, 178676.4, 180807.2,
201431.3, 22231.1, 38113.4], units='ms', t_stop=210000*pq.ms)

print(spike_time_tiling_coefficient(spiketrain_E8, spiketrain_B3, dt=0.10 * pq.s))

todo:

• added validation tests against c-implementation from https://github.com/CCutts/Detecting_pairwise_correlations_in_spike_trains/blob/master/spike_time_tiling_coefficient.c
• merge elephant-data PR: https://gin.g-node.org/NeuralEnsemble/elephant-data/pulls/4

A description of the sttc algorithm can be found here: https://doi.org/10.1523/jneurosci.2767-14.2014 .

[coveralls]

Coverage: 86.768% (+0.4%) from 86.418% when pulling 64c8ad9 on INM-6:fix/sttc into 471bf0f on NeuralEnsemble:master.

[pep8speaks]

Hello @Moritz-Alexander-Kern! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

There are currently no PEP 8 issues detected in this Pull Request. Cheers! 🍻

Comment last updated at 2023-10-30 17:08:20 UTC

[mdenker]

After pair programming review, the following was identified:

• new implementation was checked against pen and paper solution, based on the spike trains generated in the test setup. The validation test against the c implementation should check against this example as well.
• Adding a third spike train to the validation unit test would be good. For the existing validation unit test based on the regression test, it should be noted that this test has no overlaps between spikes -- therefore it tests for PA=PB=0.
• Additional comments for the new improved implementation and for the validation code calling the original c code would be nice.

[Moritz-Alexander-Kern]

* new implementation was checked against pen and paper solution, based on the spike trains generated in the test setup. The validation test against the c implementation should check against this example as well.

Added this example to generate_results.py , see PR on GIN: https://gin.g-node.org/NeuralEnsemble/elephant-data/pulls/4

* Adding a third spike train to the validation unit test would be good. For the existing validation unit test based on the regression test, it should be noted that this test has no overlaps between spikes -- therefore it tests for PA=PB=0.

Generated two spiketrains with:

# Generate two spiketrains with different firing rates
np.random.seed(0)
spiketrain_5_hz = StationaryPoissonProcess(rate=5*pq.Hz, t_start=0*pq.ms, t_stop=5000*pq.ms, refractory_period=3*pq.ms).generate_spiketrain()
spiketrain_30_hz = StationaryPoissonProcess(rate=30*pq.Hz, t_start=0*pq.ms, t_stop=5000*pq.ms, refractory_period=3*pq.ms).generate_spiketrain()

The STTC for this example using dt=100ms is: 0.5445476264987906

* Additional comments for the new improved implementation and for the validation code calling the original c code would be nice.

Added comments for run_t and run_p functions, added comments to the Cython code in wrapper_spike_time_tiling_coefficient.pyx, see PR on GIN.

[Moritz-Alexander-Kern]

build wheels