Modify error handling in function spike_contrast (spike_train_synchrony.py) #626
Conversation
Replace TypeError with a warning for cases where all input spike trains contain no more than 1 spike. This change ensures that empty spike trains return a complete data structure (trace object) for batch analysis, improving consistency and usability.
Remove the test of "All input spiketrains contain no more than 1 spike" as the ValueError was replaced by a warning.
|
Hello @ManuelCiba! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:
Comment last updated at 2024-04-17 12:54:48 UTC |
|
Hello @ManuelCiba I ran the following minimal example, with the code on branch ManuelCiba:modify_error_handling_in_spike-contrast: import neo
from quantities import ms
from elephant.spike_train_synchrony import spike_contrast
# Create list of spiektrains, where each spiketrain contains only one spike
spiketrains = [
neo.SpikeTrain([10] * ms, t_stop=1000 * ms),
neo.SpikeTrain([20] * ms, t_stop=1000 * ms),
]
for idx, spiketrain in enumerate(spiketrains):
print(f"Spike train no. {idx}: {spiketrain}")
# Run spike_contrast from elephant.spike_train_synchrony
synchrony, spike_contrast_trace = spike_contrast(
spiketrains, return_trace=True
)
print(f"Synchrony: {synchrony}")
print(f"Contrast: {spike_contrast_trace.contrast}")
print(f"Active Spiketrains: {spike_contrast_trace.active_spiketrains}")
print(f"Synchrony: {spike_contrast_trace.synchrony}")
print(f"Binsize: {spike_contrast_trace.bin_size}")which yields: Did I get the described scenario correctly with this minimal example? I assume batch processing refers to processing a number of lists of spike trains. i.e. Is this the desired output?, or should the value of The current implementation (Elephant v1.0.0) raises an error for the above script: I was expecting: PS: |
|
Hello Moritz, Thank you for the quick reply and for testing the code! You are right, I'm the author. Recently, I switched from Matlab to Python, so I'm finally using Elephant for my data analysis :) Regarding batch analysis: I have, for example, 100 recordings, each containing 60 parallel recorded spike trains. For each of the 100 recordings, I want to calculate the synchrony, having all 60 spike trains in one list. Some of the 100 recordings do not contain any spikes. However, I want to save the synchrony curve (trace.synchrony) also for these empty files, in order to generate a consistent data structure that can be used for machine learning, etc. Of course, it is also possible for the user to handle this case in their script. But I noticed that it is more convenient to receive a result from the function (since the number of bins depends on t_start and t_stop). Regarding your example: A synchrony value of 0.5 is mathematically correct considering how Spike-contrast is defined (since both spikes are within the first bin). This particular spike train could be considered as non-stationary since all the activity is just in the beginning. In the paper, I mentioned that non-stationary spike trains will lead to unreliable synchrony values. So, we could leave it to the responsibility of the user to ensure "good" spike trains(?). |
|
Hello @ManuelCiba , Thank you for your insights, now I understand more. I didn't consider the stationarity of spike trains, I was wondering about this line: elephant/elephant/spike_train_synchrony.py Line 193 in fe5053a For me the distinction between stationary and non-stationary spike trains is not obvious in this line, especially when dealing with cases where there are only one or two spikes per spiketrain. Given this ambiguity, I agree, to leave it to the user to ensure good spike trains, only the "obvious cases" are handled by returning NaN? So in conclusion, we could proceed with returning NaN for empty spike trains, as there is no sensible definition of synchrony in such cases? and raise a Warning? By adopting this approach, we maintain clarity in cases where synchrony cannot be computed due to empty spike trains while providing results for cases where spikes are present. This aligns with the preference to leave the responsibility of ensuring "nice" spike trains to the user. Currently I don't see an obvious way to handle non-stationary spike trains which lead to unreliable synchrony values. Maybe we could add a section to the documentation explaining/warning about this? Looking forward to hear your thoughts on this. |
|
For completion: Here is an example with empty spike trains using the code on branch ManuelCiba:modify_error_handling_in_spike-contrast: import neo
from quantities import ms
from elephant.spike_train_synchrony import spike_contrast
# Create list of spiektrains, where each spiketrain contains only one spike
spiketrains = [
neo.SpikeTrain([]*ms, t_stop=1000 * ms),
neo.SpikeTrain([]*ms, t_stop=1000 * ms),
]
for idx, spiketrain in enumerate(spiketrains):
print(f"Spike train no. {idx}: {spiketrain}")
# Run spike_contrast from elephant.spike_train_synchrony
synchrony, spike_contrast_trace = spike_contrast(
spiketrains, return_trace=True
)
print(f"Synchrony: {synchrony}")
print(f"Contrast: {spike_contrast_trace.contrast}")
print(f"Active Spiketrains: {spike_contrast_trace.active_spiketrains}")
print(f"Synchrony: {spike_contrast_trace.synchrony}")
print(f"Binsize: {spike_contrast_trace.bin_size}")Output: |
| try: | ||
| isi_min = min(np.diff(st).min() for st in spiketrains if len(st) > 1) | ||
| except TypeError: | ||
| raise ValueError("All input spiketrains contain no more than 1 spike.") | ||
| except: | ||
| # Do not raise an error but just a warning and continue the code which will return nan-values as synchrony. | ||
| # This is useful when running spike-contrast in a batch script, where it is necessary | ||
| # to get a complete data structure (trace for all bin-sizes) also for empty spike trains. | ||
| warnings.warn("All input spiketrains contain no more than 1 spike.") | ||
| isi_min = 0 # continue with isi_min = 0, causing to return spike-contrast with nan as synchrony values but | ||
| bin_min = max(isi_min / 2, min_bin) |
There was a problem hiding this comment.
| try: | |
| isi_min = min(np.diff(st).min() for st in spiketrains if len(st) > 1) | |
| except TypeError: | |
| raise ValueError("All input spiketrains contain no more than 1 spike.") | |
| except: | |
| # Do not raise an error but just a warning and continue the code which will return nan-values as synchrony. | |
| # This is useful when running spike-contrast in a batch script, where it is necessary | |
| # to get a complete data structure (trace for all bin-sizes) also for empty spike trains. | |
| warnings.warn("All input spiketrains contain no more than 1 spike.") | |
| isi_min = 0 # continue with isi_min = 0, causing to return spike-contrast with nan as synchrony values but | |
| bin_min = max(isi_min / 2, min_bin) | |
| if n_spikes_total == 0 : | |
| warnings.warn("All input spike trains contain no spikes.") | |
| isi_min = 0 | |
| else: | |
| isi_min = min(np.diff(st).min() for st in spiketrains if len(st) > 1) | |
| bin_min = max(isi_min / 2, min_bin) |
Depending on whether or not to handle the case with only 1 spike per spike train, we could drop the try except block, by handling empty spike trains this way?
There was a problem hiding this comment.
Hello @Moritz-Alexander-Kern
Thank you for your suggestions!
The result for the empty spike trains, looks good!
Regarding "stationarity" in this line:
elephant/elephant/spike_train_synchrony.py
Line 193 in fe5053a
You are right, (non-)stationarity is not related how we should handle the case of a single spike. It just explains, why we can get different synchrony values even if we just have one spike per spike train.
For example, these two spike trains (non-stationary)
spiketrains = [
neo.SpikeTrain([1.0]*ms, t_stop=1000 * ms),
neo.SpikeTrain([1.2]*ms, t_stop=1000 * ms),
]
will lead to another synchrony value compared to these two spike trains (stationary):
spiketrains = [
neo.SpikeTrain([250.0]*ms, t_stop=1000 * ms),
neo.SpikeTrain([750.0]*ms, t_stop=1000 * ms),
]
Since it is possible to obtain a synchrony value even with just a single spike per spike train, I would suggest, to handle this case the same way as you suggested for the empty spiketrain by setting isi_min = 0 plus the warning (If I'm not wrong, we could use my original suggestion with the try-except block?).
Just to explain, what happens if we set isi_min = 0:
Spike-contrast calculates a synchrony for many different bin sizes ("bin size sweep") and at the end selects the optimum as the final synchrony.
The smallest ISI from all spiketrains is used to define the lower limit for the bin-size sweep. If isi_min=0, the lower limit will be set to min_bin:
elephant/elephant/spike_train_synchrony.py
Line 196 in fe5053a
So, even in the case of just one spike per spike train, the algorithm can calculate a synchrony.
update copyright year
Description:
This pull request suggests an enhancement for the
spike_contrastfunction. Currently, the function raises aTypeErrorfor input spike trains containing no more than 1 spike. The proposed modification replaces theTypeErrorwith a warning to improve the user experience during batch analysis.More precisely, the function will return a complete data structure 'trace' with all bin_size values, even if the synchrony calculation was not possible. The synchrony values will be 'nan'. This will lead to consistent data structures among full spike trains and empty spike trains.
Changes:
TypeErrorwith a warning in thespike_train_synchrony.pyfunctionspike_contrast.test_spike_train_synchrony.pyfunctiontest_invalid_data.Testing:
spike_contrasthas been tested.test_spike_train_synchronywere made directly on GitHub, so no local tests have been run yet.