eFEL
eFEL copied to clipboard
BlueBrain
burst_mean_freq does not behave as expected
the feature burst_mean_freq does not behave as described in the efeature-documentation.pdf. The indices used are mixed, and the feature does not "Iterate over the burst ISI indices and yield the number of peaks divided by the length of the burst."
@wvangeit should I modify burst_mean_freq to make it stick to documentation? Or should I create a new feature (e.g. LibV5::burst_mean_frequency) to not break consistency with previous runs using burst_mean_freq?
Btw, here is the implementation that I would use (here in python for testing, I would translate it in C++ later):
def py_mean_freq(peak_time, burst_ISI_indices):
"""Python implementation of burst_mean_freq."""
# There is a discrepancy between peak_time indices and ISI_burst indices
# because 1st ISI value is ignored.
if burst_ISI_indices is None:
return None
# if no burst detected, do not consider all peaks in a single burst
elif len(burst_ISI_indices) == 0:
return []
burst_mean_freq = []
burst_index_tmp = burst_ISI_indices
burst_index = numpy.insert(
burst_index_tmp, burst_index_tmp.size, len(peak_time) - 1
)
burst_index = burst_index.astype(int)
# 1st burst
span = peak_time[burst_index[0]] - peak_time[0]
# + 1 because 1st ISI is ignored
N_peaks = burst_index[0] + 1
burst_mean_freq.append(N_peaks * 1000 / span)
for i, burst_idx in enumerate(burst_index[:-1]):
if burst_index[i + 1] - burst_idx != 1:
span = peak_time[burst_index[i + 1]] - peak_time[burst_idx + 1]
N_peaks = burst_index[i + 1] - burst_idx
burst_mean_freq.append(N_peaks * 1000 / span)
return burst_mean_freq
Do you consider it a bug in burst_mean_freq, or is it just behaving differently than the documentation? If not a bug, we should adapt the documentation. I'm also ok with a new feature with the requested implementation, but the name should be more different than your proposal though. I don't know if there is anything meaningful you could add after an underscore?
It follows the algorithm given in the documentation, but that algorithm is flawed. It does not use the right burst duration, nor the right number of peaks in the burst. So I guess I would consider it a bug.
Is it a problem of the detection of bursts, or is it only the mean_freq that is affected?
The previous implementation of the feature assumed that if a burst is detected, there was probably a 1st burst before that one. I kept that in the new implementation, but after thinking about it with @DrTaDa , we realized that it might not always be the case. We'll try to fix that asap.
@wvangeit
The problem: For the efeature burst_mean_freq, if a trace looks like spike1 -> spike2 -> spike3 -> burst1, the code will assume that there is a burst made of spike2 and spike3, which will mess up the computation of the feature. This is done (we think) because BurstIndex cannot detect a burst if it is at the very beginning of a trace because of its implementation, therefore, this is "fixed" in burst_mean_freq by assuming that there was a burst before.
Solution: one way to solve this would be to change the way bursts are detected in the first place (for example by keeping in memory the indices of which spikes belong to which burst and using a global median instead of a local one to tell the bursts apart). However, that could also have its lot of issues.
We are not sure if we should try to (1) patch-up the current implementation, (2) leave it as is and accept that it is flawed from time to time (3) rework everything linked to burst detection entirely using the approach described above.
yes, isn't the current definition of burst anyway that it splits the trace in blocks that start/stop when there is significant change in ISI length? Doesn't it in that context make sense that spike2 and spike3 are considered a burst themselves?
In the example @DrTaDa provided, I think the current implementation could actually create another problem: the long ISI between the 1st spikes could bias the median, and thus we might not be able to even detect the 1st burst occurring after the spikes.
Another clearer example of failure would be an isolated spike, followed by a first burst, and then a second burst, etc. The current implementation would detect the second burst, and assume that everything before was also one single burst, thus merging the isolated spike and the 1st burst.
I have a new strategy to detect the bursts. Tell me if it looks decent enough to you.
In the current implementation, we detect the beginning of a burst when the current ISI is larger than the burst factor multiplied by the median of previous ISIs, and when the next ISI is smaller than the current ISI divided by the burst factor.
if (ISIValues[i] > (BurstFactor * dMedian) && (ISIValues[i + 1] < ISIValues[i] / BurstFactor))
{BurstIndex.push_back(i + 1);}
The end of a burst is the peak before the beginning of the next burst.
This implementation has the following problem: individual spikes that should not belong to any bursts are automatically assigned to a burst.
I thought about decoupling the two parts of the previous formula to detect the beginning and the end of bursts:
- we define the end of a burst when the current ISI is larger than the burst factor multiplied by the median of previous ISIs
- and we define the beginning of a burst when the next ISI is smaller than the current ISI divided by the burst factor.
if (in_burst && ISI_values[i] > (burst_factor * dMedian)){
burst_end_indices.push_back(i);}
if (ISI_values[i] < ISI_values[i - 1] / burst_factor){
if (in_burst){
burst_begin_indices.back() = i;
} else {
burst_begin_indices.push_back(i);
}
}
With this implementation, we precisely know where each burst begins and ends, and that solves our isolated spike problem.
However, with the current default value of the burst_factor of 1.5, it sometimes exclude spikes from bursts they should belong to. Setting burst_factor to at least 2 solves this with the traces I tested this on.
Is there a deep reason we set the default value to 1.5 in the first place? What do you think of this implementation idea?
I agree with this implementation, it is similar to the custom one I made for Polina to compute the number of spike per burst.