compute_synchrony_metrics update #2605

chrishalcrow · 2024-03-20T10:57:24Z

(Another go of #2528, after the SortingAnalyzer refactor)

Planning to standardise, correctness test and optimise the quality metrics module. First one: the compute_synchrony_metrics function.

Two main things in this:

The generating function synthesize_random_firings and synthesize_poisson_spike_vectors were adjusted so that you can insert specific spikes, useful for testing. These are used in the new tests.
A new function compute_synchrony_counts is introduced, which computes the synchrony counts.

One thing I'm not sure about: The documentation said that compute_synchrony_metrics returns a dict, but I think we should export a namedtuple, which matches most other quality metrics.

I adjusted the algorithm in compute_synchrony_counts, which should be easier to understand and is significantly faster. Here are some benchmarking results, with times and create a log-log plot to check the algorithm's complexity. About 30x faster!

Benchmarked using the following code:

import numpy as np
from time import perf_counter

from spikeinterface.core import synthesize_random_firings
from spikeinterface.qualitymetrics import get_synchrony_counts

unit_ids = np.arange(100)

Ns = 10**np.array(np.arange(3,5.1,0.1))

times_by_size=[]
how_many = 100

for train_length in Ns:

    times_list = []

    for _ in range(0,how_many):

        spike_train = synthesize_random_firings(num_units=100, duration=train_length/1000, firing_rates=10.0)

        start = perf_counter()
        get_synchrony_counts(spike_train, np.array((2,4,8)), list(unit_ids))        
        end = perf_counter()

        times_list.append(end - start)

    times_by_size.append( np.median( times_list ) )

for more information, see https://pre-commit.ci

alejoe91 · 2024-03-20T14:04:29Z

Thanks @chrishalcrow

Can you heck why tests are failing?

chrishalcrow · 2024-03-20T15:45:46Z

Thanks @chrishalcrow

Can you heck why tests are failing?

Yes, it's to do with me originally changing the output from namedtuple to a dict (as documented). I'm changing back now...

for more information, see https://pre-commit.ci

src/spikeinterface/core/generate.py

h-mayorquin · 2024-03-20T16:38:19Z

src/spikeinterface/core/generate.py


-    return (times, labels)
+    return structured_insertions


Same coment here about the output. Input and output are your contracts with your users. We should be careful with changing them (and careful to define them in the first place)

Indeed, I'll add some documentation. I think this (and the add_insertions) functions don't have much use outside of testing. Should I "hide" them (i.e. make them _add_insertions etc)?

If you think that they will only be used for internal testing then the best is to have them in the corresponding testing module I feel.

I think hiding by default is a good option and then you make it public once you feel confident about the interface and purpose. That minimizes the risk if you want to change the function later. But this is my opinion, I don't think that anyone else would find any problems if you want to add some functions.

I'd prefer to keep them close to the synthesize_stuff functions, since they are closely related. I'm a fan of the "related code should be close together" mindset. Maybe they should all be in the testing module? But maybe for another time...

Great, I've hidden it and also added some documentation. Best of all worlds.

I think this is a good decision. If we ever use the function more and we are decided on the signature, then we can make it public.

src/spikeinterface/core/generate.py

h-mayorquin · 2024-03-20T17:15:35Z

src/spikeinterface/qualitymetrics/misc_metrics.py

+    ----------
+    spikes : memmap
+        The spike train
+    synchrony_sizes : numpy array


Is this given in seconds, samples, milliseconds?

I've made the description more accruate (since it's actually a structured numpy array). One of the fields is "sample_index" which clarifies the unit.

I was talking about the synchrony_sizes

Oh, sorry. These are the number of synchronous events you want to count. So if you want to see when two or four spikes fire at the same time you use synchrony_sizes = (2,4). So it's an integer, and I think it will be clear for anyone who knows enough about the metric that their using it.

Thanks. I guess my lack of familarity with the metric is evident then. I somehow imagined that it was the windows in which an event would be counted as syncronous. Thanks for explaining it.

h-mayorquin · 2024-03-20T17:15:47Z

Tip:
You have https://numpy.org/doc/stable/reference/generated/numpy.logspace.html to generate your samples.

Benchamark comments:

I think you should fix your seed in your benchmark explictly.
I find the spike_length variable hard to think about. I suggest using time and show us the units. Does it scale well to the relevant domain of hour long sortings?
Finally, I would be curios to see how this scales with the number of units as well.

chrishalcrow · 2024-03-21T14:25:10Z

Tip: You have https://numpy.org/doc/stable/reference/generated/numpy.logspace.html to generate your samples.

Benchamark comments:

I think you should fix your seed in your benchmark explictly.

I find the spike_length variable hard to think about. I suggest using time and show us the units. Does it scale well to the relevant domain of hour long sortings?

Finally, I would be curios to see how this scales with the number of units as well.

Oooh, thanks for the logspace function.

Maybe it would be helpful to say: the new algorithm is expected to be faster. I was able to replace the use of np.complexity, isin and count_nonzero (I think complexity is the heavy function of these) with some simple selections (like unique_spike_index[counts >= 2]).

Spike_length is just the number of spikes in the spike vector. So for generic data, this increases linearly with time.

I’ve been using a 30 minute 16 channel recording from my lab, which contains 70,000 spikes after sorting. Is this realistic? I’ve also been using this to make sure my algorithm is consistent with the old results. My tests go up to 100,000 spikes. What would a good number to go up to?

I'm testing with unit number now, going up to 1000 units. Is that high enough? The complexity function changes it’s method at some point, and I’m trying to figure out how and when and will report back soon! But the new code is still substantially faster at e.g. 1k units with 100k spikes.

h-mayorquin · 2024-03-21T20:34:31Z

Spike_length is just the number of spikes in the spike vector. So for generic data, this increases linearly with time.

On a second thought, this seems more amenable to control and more relevant than duration. It is just probably me that does not have a good intuitive grasp of this metric while duration I understand better. How are you controlling this though, given that you are passing this parameter in duration divided by 1000?

I'm testing with unit number now, going up to 1000 units. Is that high enough? The complexity function changes it’s method at some point, and I’m trying to figure out how and when and will report back soon! But the new code is still substantially faster at e.g. 1k units with 100k spikes.

1000 I have been told than 1000 units is like a very high number.

Bear in mind that the comments I made about the benchmark is just me being curios and suggestions. As long as Sam and Alessio are OK with the rest I would not bother if you don't care about it. The thing that I really feel strongly about is for not changing the signature and output of synthesize_poisson_spike_vector.

for more information, see https://pre-commit.ci

chrishalcrow · 2024-03-26T09:41:55Z

Spike_length is just the number of spikes in the spike vector. So for generic data, this increases linearly with time.

On a second thought, this seems more amenable to control and more relevant than duration. It is just probably me that does not have a good intuitive grasp of this metric while duration I understand better. How are you controlling this though, given that you are passing this parameter in duration divided by 1000?

I'm testing with unit number now, going up to 1000 units. Is that high enough? The complexity function changes it’s method at some point, and I’m trying to figure out how and when and will report back soon! But the new code is still substantially faster at e.g. 1k units with 100k spikes.

1000 I have been told than 1000 units is like a very high number.

Bear in mind that the comments I made about the benchmark is just me being curios and suggestions. As long as Sam and Alessio are OK with the rest I would not bother if you don't care about it. The thing that I really feel strongly about is for not changing the signature and output of synthesize_poisson_spike_vector.

Hello,

I learned a few things from your comments and trying to recreate my previous benchmarking. Thanks! Before: I refactored the code, creating a new function compute_synchrony_count. Then I did a speed test on the new function using the old algorithm against an optimised version of it. This is hard for others to reproduce since you can’t easily run compute_synchrony_count. So now I’m benchmarking against the full compute_synchrony_metrics function. This should be closer to “real world” use than before and easier for anyone to check. And I’ve hopefully made the wording/details of the plots more obvious.

The old benchmarking results are basically unchanged.

And I’ve now done more benchmarking, with changing units. First, here’s a plot of timings with increasing units N but fixed spike train length (of about 50k: in our lab, this is the number of spikes in a typical 20 minute recording). The new algorithm grows linearly with the number of spikes. The old algorithm, which uses more complex functions, changes tactic at ~600 units (I think it’s np.unique that changes its tactic…). Both its tactics look like O(N) algorithms. But even when it changes, it is still about 10x slower than the new one.

I also tried keeping the duration and firing rates fixed. So here, as the number of units N increases the length of the spike train increases at the same rate. But the density of spikes per time unit stays constant. Both algorithms grow as N^2, as expected (since units and train length are both increasing with N). Again, the old algorithm changes tactic but the new code is still significantly faster.

The code for both is below:

times_by_size=[]
sizes = []

how_many = 100

duration = 10**1.7
# duration = 10**1.5 for the first plot
nums_units = np.arange(100,1201,100)

for num_units in nums_units:

    recording, sorting = generate_ground_truth_recording(
            durations=[
               # duration # for the first plot
               # duration/num_units*nums_units[0] # for the second plot
            ],
            sampling_frequency=30_000.0,
            num_channels=6,
            num_units=num_units,
            generate_sorting_kwargs=dict(firing_rates=10.0, refractory_period_ms=4.0),
            noise_kwargs=dict(noise_level=5.0, strategy="tile_pregenerated"),
            seed=2205,
        )

    sorting_analyzer = create_sorting_analyzer(sorting, recording, format="memory", sparse=False)

    times_list = []
    for _ in range(0,how_many):

        start = perf_counter()
        compute_synchrony_metrics(sorting_analyzer)        
        end = perf_counter()

        times_list.append(end - start)

    times_by_size.append( np.median( times_list ) )
    sizes.append( len(sorting_analyzer.sorting.to_spike_vector()) )

h-mayorquin

Hey, thanks for removing the change of the output and doing benchmark with units.

I sitll would like to keep the function for generating spikes as simple as possible so I request a removal or the arguments. As I mentioned before:

The thing that I really feel strongly about is for not changing the signature and output of synthesize_poisson_spike_vector.

Unfortunately, I can't review the synchrony metrics changes as I am not really well acquainted with the method. Maybe @DradeAW or @zm711 can take a look?

src/spikeinterface/core/generate.py

h-mayorquin · 2024-03-26T16:31:49Z

src/spikeinterface/qualitymetrics/misc_metrics.py

+    ----------
+    spikes : memmap
+        The spike train
+    synchrony_sizes : numpy array


I was talking about the synchrony_sizes

src/spikeinterface/core/generate.py

zm711 · 2024-03-26T17:01:43Z

Unfortunately, I can't review the synchrony metrics changes as I am not really well acquainted with the method. Maybe @DradeAW or @zm711 can take a look?

@alejoe91 actually ported this from Elephant so he would be the best to review unless @DradeAW knows this super well. I never actually use this metric :)

chrishalcrow · 2024-03-27T15:15:21Z

Unfortunately, I can't review the synchrony metrics changes as I am not really well acquainted with the method. Maybe @DradeAW or @zm711 can take a look?

@alejoe91 actually ported this from Elephant so he would be the best to review unless @DradeAW knows this super well. I never actually use this metric :)

An important part of the update are the four test_synchrony_counts_blah in test_metrics_functions.py . These are correctness tests for the main part of the algorithm. I make a small spike_train with certain spikes (e.g. two spikes which fire at the same time) and make sure the metric gives the right answer. If anyone can think of any good edge cases to test here, please comment!

for more information, see https://pre-commit.ci

h-mayorquin

Thanks for bearing with me and for adding thorough testing.

checking if this removes the changes requested tag

src/spikeinterface/qualitymetrics/misc_metrics.py

src/spikeinterface/core/generate.py

for more information, see https://pre-commit.ci

from int64 to int16 Co-authored-by: Alessio Buccino <alejoe9187@gmail.com>

alejoe91 · 2024-04-15T08:36:04Z

src/spikeinterface/core/generate.py

+def _add_spikes_to_spiketrain(
+    spike_indices,
+    spike_labels,
+    segment_indices=[],


Suggested change

segment_indices=[],

segment_indices=None,

See #2345

To quote Joe: "Thanks I did not know about that, I am learning so much working on the SI codebase 😍!"

alejoe91 · 2024-04-15T08:36:13Z

src/spikeinterface/core/generate.py

+    segment_indices=[],
+    added_spikes_indices=None,
+    added_spikes_labels=None,
+    added_segment_indices=[],


Suggested change

added_segment_indices=[],

added_segment_indices=None,

for more information, see https://pre-commit.ci

chrishalcrow · 2024-04-15T12:49:28Z

Hello, realised that the _add_spikes_to_spiketrain function was basically just checking a bunch of array lengths then concatenating some arrays together, that were then used for testing. So it wasn't really adding much, and it was actually much simpler to get rid of it and concatenate in the tests directly. Got rid of 200 lines of code, and simplified the pull request a lot!

alejoe91 · 2024-04-15T15:27:30Z

src/spikeinterface/qualitymetrics/misc_metrics.py

+
+    synchrony_metrics_dict = {}
+    for sync_idx, synchrony_size in enumerate(synchrony_sizes_np):
+        synch_id_metrics_dict = {}


Suggested change

synch_id_metrics_dict = {}

sync_id_metrics_dict = {}

Can we stick to the same naming? :)

Sorry for being overly annoying @chrishalcrow !!!

You can name it anything you want if I can stop working on this ;) Just you wait: my next quality_metrics update is gonna be so careful!

samuelgarcia · 2024-04-22T18:01:43Z

Hi.
When we will have time I would be happy to discuss this function.
The synchrony here if I understand correctly is the exact same timestamp but we could have synchrony with very small lags. A lag can be due to small noise and small time shift.
So my main question is how robut is that mertics regarding theses lags ?

samuelgarcia · 2024-04-22T18:02:47Z

src/spikeinterface/qualitymetrics/misc_metrics.py

@@ -496,7 +496,51 @@ def compute_sliding_rp_violations(
 )


-def compute_synchrony_metrics(sorting_analyzer, synchrony_sizes=(2, 4, 8), unit_ids=None, **kwargs):
+def get_synchrony_counts(spikes, synchrony_sizes, all_unit_ids):


I would excpet a delta in the signature no ?

samuelgarcia · 2024-04-22T18:03:11Z

src/spikeinterface/qualitymetrics/misc_metrics.py

+    This code was adapted from `Elephant - Electrophysiology Analysis Toolkit <https://github.com/NeuralEnsemble/elephant/blob/master/elephant/spike_train_synchrony.py#L245>`_
+    """
+
+    synchrony_counts = np.zeros((np.size(synchrony_sizes), len(all_unit_ids)), dtype=np.int64)


In this the synhcony is higly dependant on the sampling rate no ?

chrishalcrow · 2024-04-23T07:06:49Z

Hi. When we will have time I would be happy to discuss this function. The synchrony here if I understand correctly is the exact same timestamp but we could have synchrony with very small lags. A lag can be due to small noise and small time shift. So my main question is how robut is that mertics regarding theses lags ?

Hello, indeed the metric is not robust to lags at all. This PR was just optimizing what already existed. I'm happy to add a delta, which would check if there were spikes within delta samples of the first one?

alejoe91 · 2024-04-23T07:10:46Z

Let's keep it like this for now, since this was the original design. :)

chrishalcrow and others added 4 commits March 20, 2024 10:50

Update and test compute_synchrony_metrics

e8e4ca4

[pre-commit.ci] auto fixes from pre-commit.com hooks

220add7

for more information, see https://pre-commit.ci

Make namedtuple and check div by zero

cb45927

[pre-commit.ci] auto fixes from pre-commit.com hooks

b91f5f2

for more information, see https://pre-commit.ci

alejoe91 added the qualitymetrics Related to qualitymetrics module label Mar 20, 2024

alejoe91 added refactor Refactor of code, with no change to functionality performance Performance issues/improvements labels Mar 20, 2024

chrishalcrow and others added 2 commits March 20, 2024 16:18

Redo tests to use namedtuple

bfbc5e3

[pre-commit.ci] auto fixes from pre-commit.com hooks

f92e2a9

for more information, see https://pre-commit.ci

h-mayorquin previously requested changes Mar 20, 2024

View reviewed changes

h-mayorquin reviewed Mar 20, 2024

View reviewed changes

chrishalcrow and others added 3 commits March 26, 2024 09:27

Merge branch 'SpikeInterface:main' into sync_counts_update

4f5b119

Revert generate_stuff outputs, add tests + docstrings

1bb06b6

[pre-commit.ci] auto fixes from pre-commit.com hooks

0706b4d

for more information, see https://pre-commit.ci

chrishalcrow added 2 commits March 26, 2024 09:46

delete mentions of segments in in synthesize_things doc

635d330

Change "spikes" description

f22698b

chrishalcrow requested a review from h-mayorquin March 26, 2024 14:59

h-mayorquin reviewed Mar 26, 2024

View reviewed changes

type bins for numba

56f0559

chrishalcrow and others added 3 commits March 27, 2024 15:23

revert to old synthesize_stuff functions (before cleanup)

3de7c27

Merge branch 'main' into sync_counts_update

c53fedd

[pre-commit.ci] auto fixes from pre-commit.com hooks

c9244a2

for more information, see https://pre-commit.ci

chrishalcrow added 2 commits March 27, 2024 15:34

oops, remove isi.py edits

7cfd0ea

Merge branch 'main' into sync_counts_update

f3217f4

h-mayorquin reviewed Mar 27, 2024

View reviewed changes

Merge branch 'main' into sync_counts_update

51f0039

alejoe91 reviewed Apr 12, 2024

View reviewed changes

src/spikeinterface/qualitymetrics/misc_metrics.py Outdated Show resolved Hide resolved

alejoe91 reviewed Apr 12, 2024

View reviewed changes

src/spikeinterface/core/generate.py Outdated Show resolved Hide resolved

alejoe91 reviewed Apr 12, 2024

View reviewed changes

src/spikeinterface/core/generate.py Outdated Show resolved Hide resolved

chrishalcrow and others added 4 commits April 12, 2024 12:25

rename insertions and unpack spike train into indices, labels, segments

f5d20fa

[pre-commit.ci] auto fixes from pre-commit.com hooks

23c3355

for more information, see https://pre-commit.ci

Merge branch 'main' into sync_counts_update

7719ee1

Update src/spikeinterface/qualitymetrics/misc_metrics.py

66077e2

from int64 to int16 Co-authored-by: Alessio Buccino <alejoe9187@gmail.com>

alejoe91 reviewed Apr 15, 2024

View reviewed changes

chrishalcrow and others added 6 commits April 15, 2024 10:58

Change to None, add some ifs

11d5d39

[pre-commit.ci] auto fixes from pre-commit.com hooks

5de7b66

for more information, see https://pre-commit.ci

Revert namedtuple name

378976f

Merge branch 'main' into sync_counts_update

ecc6a9a

Delete _add_spikes_to_spiketrain function, for simplicity

f6652f0

[pre-commit.ci] auto fixes from pre-commit.com hooks

0b9a58a

for more information, see https://pre-commit.ci

alejoe91 reviewed Apr 15, 2024

View reviewed changes

synch -> sync

16b344c

alejoe91 approved these changes Apr 16, 2024

View reviewed changes

alejoe91 merged commit fa57fee into SpikeInterface:main Apr 16, 2024
11 checks passed

samuelgarcia reviewed Apr 22, 2024

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

compute_synchrony_metrics update #2605

compute_synchrony_metrics update #2605

chrishalcrow commented Mar 20, 2024 •

edited

Loading

alejoe91 commented Mar 20, 2024

chrishalcrow commented Mar 20, 2024

h-mayorquin Mar 20, 2024

chrishalcrow Mar 21, 2024

h-mayorquin Mar 21, 2024

chrishalcrow Mar 26, 2024

h-mayorquin Mar 27, 2024

h-mayorquin Mar 20, 2024

chrishalcrow Mar 26, 2024

h-mayorquin Mar 26, 2024

chrishalcrow Mar 27, 2024

h-mayorquin Mar 27, 2024 •

edited

Loading

h-mayorquin commented Mar 20, 2024

chrishalcrow commented Mar 21, 2024

h-mayorquin commented Mar 21, 2024

chrishalcrow commented Mar 26, 2024

h-mayorquin left a comment •

edited

Loading

h-mayorquin Mar 26, 2024

zm711 commented Mar 26, 2024

chrishalcrow commented Mar 27, 2024

h-mayorquin left a comment

alejoe91 Apr 15, 2024

chrishalcrow Apr 15, 2024

alejoe91 Apr 15, 2024

chrishalcrow commented Apr 15, 2024

alejoe91 Apr 15, 2024

chrishalcrow Apr 15, 2024

samuelgarcia commented Apr 22, 2024

samuelgarcia Apr 22, 2024

samuelgarcia Apr 22, 2024

chrishalcrow commented Apr 23, 2024

alejoe91 commented Apr 23, 2024

compute_synchrony_metrics update #2605

compute_synchrony_metrics update #2605

Conversation

chrishalcrow commented Mar 20, 2024 • edited Loading

alejoe91 commented Mar 20, 2024

chrishalcrow commented Mar 20, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

h-mayorquin Mar 27, 2024 • edited Loading

Choose a reason for hiding this comment

h-mayorquin commented Mar 20, 2024

chrishalcrow commented Mar 21, 2024

h-mayorquin commented Mar 21, 2024

chrishalcrow commented Mar 26, 2024

h-mayorquin left a comment • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

zm711 commented Mar 26, 2024

chrishalcrow commented Mar 27, 2024

h-mayorquin left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

chrishalcrow commented Apr 15, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

samuelgarcia commented Apr 22, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

chrishalcrow commented Apr 23, 2024

alejoe91 commented Apr 23, 2024

chrishalcrow commented Mar 20, 2024 •

edited

Loading

h-mayorquin Mar 27, 2024 •

edited

Loading

h-mayorquin left a comment •

edited

Loading