From 53733b024994defa47ce2ef50925fad6d2c9a1ed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Aur=C3=A9lien=20Jaquier?= <72930209+AurelienJaquier@users.noreply.github.com> Date: Mon, 29 Apr 2024 16:54:02 +0200 Subject: [PATCH] new features: postburst_slow_ahp_values, time_to_postburst_slow_ahp (#289) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * new features: postburst_slow_ahp_values, time_to_postburst_slow_ahp * remove std::cout in LibV1 * fix edge case when index is stim_end_index or end_index * new feature: postburst_fast_ahp_values * new feature: postburst_adp_peak_values * add new features: interburst_XX_percent_values for 20, 40, 60% * adding test, docs for postburst_fast_ahp_values, postburst_adp_peak_values * fix time_to_postburst_adp_peak for edge case when there are less adp peaks than bursts * add test and doc for feature time_to_postburst_fast_ahp * add test and doc for feature time_to_postburst_adp_peak * add test and docs to interburst 20%, 40%, 60% features * turn 40%, 60% into 25%, 30% for interburst voltage * add interburst 10% and 15% voltage values features * interburst features: removed 10%, restored 40% and 60% * new feature: interburst_duration * update docs and tests with new features * update memoization of new features * fix docs, units and tests for latest new features * create interburst_percent_indices function to reduce code duplication * remove deprecated bind2nd * refactor interburst_XXpercent_indices * merging all interburst_XXpercent_values in the docs --------- Co-authored-by: Jaquier Aurélien Tristan --- docs/source/eFeatures.rst | 207 +++++- efel/DependencyV5.txt | 24 +- efel/cppcore/FillFptrTable.cpp | 70 ++ efel/cppcore/LibV5.cpp | 478 ++++++++++++- efel/cppcore/LibV5.h | 34 + efel/cppcore/cfeature.cpp | 22 + efel/units/units.json | 13 + tests/featurenames.json | 23 + tests/test_basic.py | 666 +++++++++++++++++- .../testdata/allfeatures/expectedresults.json | 130 +++- 10 files changed, 1643 insertions(+), 24 deletions(-) diff --git a/docs/source/eFeatures.rst b/docs/source/eFeatures.rst index 893ad81c..67bac2fc 100644 --- a/docs/source/eFeatures.rst +++ b/docs/source/eFeatures.rst @@ -534,7 +534,7 @@ The burst detection can be fine-tuned by changing the setting strict_burst_facto - **Units**: mV - **Pseudocode**: :: - interburst_min = [ + postburst_min = [ numpy.min( v[peak_indices[i]:peak_indices[i + 1]] ) for i in burst_end_indices if i + 1 < len(peak_indices) @@ -551,6 +551,36 @@ The burst detection can be fine-tuned by changing the setting strict_burst_facto v[peak_indices[burst_end_indices[-1]]:] )) +`LibV5`_ : postburst_slow_ahp_values +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The slow AHP voltage after the end of a burst. + +The number of ms to skip after the spike to skip fast AHP and look for slow AHP can be set with sahp_start. +Default is 5. + +This implementation does not assume that every spike belongs to a burst. + +The first spike is ignored by default. This can be changed by setting ignore_first_ISI to 0. + +The burst detection can be fine-tuned by changing the setting strict_burst_factor. Defalt value is 2.0. + +- **Required features**: peak_indices, burst_end_indices +- **Units**: mV +- **Pseudocode**: :: + + postburst_slow_ahp = [] + for i in burst_end_indices: + i_start = numpy.where(t >= t[peak_indices[i]] + sahp_start)[0][0] + if i + 1 < len(peak_indices): + postburst_slow_ahp.append(numpy.min(v[i_start:peak_indices[i + 1]])) + else: + if t[burst_end_indices[-1]] < stim_end: + end_idx = numpy.where(t >= stim_end)[0][0] + postburst_slow_ahp.append(numpy.min(v[i_start:end_idx])) + else: + postburst_slow_ahp.append(numpy.min(v[i_start:])) + `LibV5`_ : time_to_interburst_min ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -573,6 +603,181 @@ The burst detection can be fine-tuned by changing the setting strict_burst_facto for i in burst_end_indices if i + 1 < len(peak_indices) ] +`LibV5`_ : time_to_postburst_slow_ahp +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The time between the last spike of a burst and the slow ahp afterwards. + +The number of ms to skip after the spike to skip fast AHP and look for slow AHP can be set with sahp_start. +Default is 5. + +This implementation does not assume that every spike belongs to a burst. + +The first spike is ignored by default. This can be changed by setting ignore_first_ISI to 0. + +The burst detection can be fine-tuned by changing the setting strict_burst_factor. Defalt value is 2.0. + +- **Required features**: postburst_slow_ahp_indices, burst_end_indices, peak_time +- **Units**: ms +- **Pseudocode**: :: + + time_to_postburst_slow_ahp_py = t[postburst_slow_ahp_indices] - peak_time[burst_end_indices] + +`LibV5`_ : postburst_fast_ahp_values +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The fast AHP voltage after the end of a burst. + +This implementation does not assume that every spike belongs to a burst. + +The first spike is ignored by default. This can be changed by setting ignore_first_ISI to 0. + +The burst detection can be fine-tuned by changing the setting strict_burst_factor. Defalt value is 2.0. + +- **Required features**: peak_indices, burst_end_indices +- **Units**: mV +- **Pseudocode**: :: + + postburst_fahp = [] + for i in burst_end_indices: + if i + 1 < len(peak_indices): + stop_i = peak_indices[i + 1] + elif i + 1 < stim_end_index: + stop_i = stim_end_index + else: + stop_i = len(v) - 1 + + v_crop = v[peak_indices[i]:stop_i] + # get where the voltage is going up + crop_args = numpy.argwhere(numpy.diff(v_crop) >= 0)[:,0] + # the voltage should go up for at least two consecutive points + crop_arg_arg = numpy.argwhere(numpy.diff(crop_args) == 1)[0][0] + crop_arg = crop_args[crop_arg_arg] + end_i = peak_indices[i] + crop_arg + 1 + # the fast ahp is between last peak of burst and the point where voltage is going back up + postburst_fahp.append(numpy.min(v[peak_indices[i]:end_i])) + + return postburst_fahp + +`LibV5`_ : postburst_adp_peak_values +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The small ADP peak after the fast AHP after the end of a burst. + +This implementation does not assume that every spike belongs to a burst. + +The first spike is ignored by default. This can be changed by setting ignore_first_ISI to 0. + +The burst detection can be fine-tuned by changing the setting strict_burst_factor. Defalt value is 2.0. + +- **Required features**: postburst_fast_ahp_indices, postburst_slow_ahp_indices +- **Units**: mV +- **Pseudocode**: :: + + adp_peak_values = [] + for i, sahpi in enumerate(postburst_sahpi): + if sahpi < postburst_fahpi[i]: + continue + adppeaki = numpy.argmax(v[postburst_fahpi[i]:sahpi]) + postburst_fahpi[i] + if adppeaki != sahpi - 1: + adp_peak_values.append(v[adppeaki]) + + if len(adp_peak_values) == 0: + return None + return adp_peak_values + +`LibV5`_ : time_to_postburst_fast_ahp +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Time to the fast AHP after the end of a burst. + +This implementation does not assume that every spike belongs to a burst. + +The first spike is ignored by default. This can be changed by setting ignore_first_ISI to 0. + +The burst detection can be fine-tuned by changing the setting strict_burst_factor. Defalt value is 2.0. + +- **Required features**: postburst_fast_ahp_indices, burst_end_indices, peak_time +- **Units**: ms +- **Pseudocode**: :: + + [t[fahpi] - peak_time[burst_endi[i]] for i, fahpi in enumerate(postburst_fahpi)] + +`LibV5`_ : time_to_postburst_adp_peak +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Time to the small ADP peak after the fast AHP after the end of a burst. + +This implementation does not assume that every spike belongs to a burst. + +The first spike is ignored by default. This can be changed by setting ignore_first_ISI to 0. + +The burst detection can be fine-tuned by changing the setting strict_burst_factor. Defalt value is 2.0. + +- **Required features**: postburst_adp_peak_indices, burst_end_indices, peak_time +- **Units**: ms +- **Pseudocode**: :: + + time_to_postburst_adp_peaks = [] + n_peaks = len(peak_time) + for i, adppeaki in enumerate(postburst_adppeaki): + # there are not always an adp peak after each burst + # so make sure that the burst and adp peak indices are consistent + k = 0 + while ( + burst_endi[i] + k + 1 < n_peaks and peak_time[burst_endi[i] + k + 1] < t[adppeaki] + ): + k += 1 + + time_to_postburst_adp_peaks.append(t[adppeaki] - peak_time[burst_endi[i] + k]) + + return time_to_postburst_adp_peaks + + +`LibV5`_ : interburst_15percent_values, interburst_20percent_values, interburst_25percent_values, interburst_30percent_values, interburst_40percent_values, interburst_60percent_values +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Voltage value after a given percentage (15%, 20%, 25%, 30%, 40% or 60%) of the interburst duration after the fast AHP. + +This implementation does not assume that every spike belongs to a burst. + +The first spike is ignored by default. This can be changed by setting ignore_first_ISI to 0. + +The burst detection can be fine-tuned by changing the setting strict_burst_factor. Defalt value is 2.0. + +- **Required features**: postburst_fast_ahp_indices, burst_end_indices, peak_indices +- **Units**: mV +- **Pseudocode**: :: + + interburst_XXpercent_values = [] + for i, postburst_fahp_i in enumerate(postburst_fahpi): + if i < len(burst_endi) and burst_endi[i] + 1 < len(peaki): + time_interval = t[peaki[burst_endi[i] + 1]] - t[postburst_fahp_i] + time_at_XXpercent = t[postburst_fahp_i] + time_interval * percentage / 100. + index_at_XXpercent = numpy.argwhere(t >= time_at_XXpercent)[0][0] + interburst_XXpercent_values.append(v[index_at_XXpercent]) + +`LibV5`_ : interburst_duration +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Duration between the last spike of each burst and the next spike. + +This implementation does not assume that every spike belongs to a burst. + +The first spike is ignored by default. This can be changed by setting ignore_first_ISI to 0. + +The burst detection can be fine-tuned by changing the setting strict_burst_factor. Defalt value is 2.0. + +- **Required features**: burst_end_indices, peak_time +- **Units**: ms +- **Pseudocode**: :: + + interburst_duration = [ + peak_time[idx + 1] - peak_time[idx] + for idx in burst_end_indices + if idx + 1 < len(peak_time) + ] + `Python efeature`_ : single_burst_ratio ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/efel/DependencyV5.txt b/efel/DependencyV5.txt index 12a4c4cb..fe43da29 100644 --- a/efel/DependencyV5.txt +++ b/efel/DependencyV5.txt @@ -114,5 +114,27 @@ LibV5:ADP_peak_amplitude #LibV5:ADP_peak_values #LibV5:min_AHP_values #LibV1:in LibV5:interburst_min_indices #LibV5:peak_indices #LibV5:burst_end_indices #LibV1:interpolate LibV5:interburst_min_values #LibV5:interburst_min_indices #LibV1:interpolate LibV5:postburst_min_indices #LibV5:peak_indices #LibV5:burst_end_indices #LibV1:interpolate -LibV5:postburst_min_values #LibV5:postburst_min_indices #LibV1:interpolate +LibV5:postburst_min_values #LibV5:postburst_min_indices #LibV1: interpolate +LibV5:postburst_slow_ahp_indices #LibV5:peak_indices #LibV5:burst_end_indices #LibV1:interpolate +LibV5:postburst_slow_ahp_values #LibV5:postburst_slow_ahp_indices #LibV1:interpolate LibV5:time_to_interburst_min #LibV5:interburst_min_indices #LibV1:peak_time #LibV5:burst_end_indices #LibV1:interpolate +LibV5:time_to_postburst_slow_ahp #LibV5:postburst_slow_ahp_indices #LibV1:peak_time #LibV5:burst_end_indices #LibV1:interpolate +LibV5:postburst_fast_ahp_indices #LibV5:peak_indices #LibV5:burst_end_indices #LibV1:interpolate +LibV5:postburst_fast_ahp_values #LibV5:postburst_fast_ahp_indices #LibV1:interpolate +LibV5:postburst_adp_peak_indices #LibV5:postburst_fast_ahp_indices #LibV5:postburst_slow_ahp_indices #LibV1:interpolate +LibV5:postburst_adp_peak_values #LibV5:postburst_adp_peak_indices #LibV1:interpolate +LibV5:time_to_postburst_fast_ahp #LibV5:postburst_fast_ahp_indices #LibV1:peak_time #LibV5:burst_end_indices #LibV1:interpolate +LibV5:time_to_postburst_adp_peak #LibV5:postburst_adp_peak_indices #LibV1:peak_time #LibV5:burst_end_indices #LibV1:interpolate +LibV5:interburst_15percent_indices #LibV5:peak_indices #LibV5:burst_end_indices #LibV5:postburst_fast_ahp_indices #LibV1:interpolate +LibV5:interburst_15percent_values #LibV5:interburst_15percent_indices #LibV1:interpolate +LibV5:interburst_20percent_indices #LibV5:peak_indices #LibV5:burst_end_indices #LibV5:postburst_fast_ahp_indices #LibV1:interpolate +LibV5:interburst_20percent_values #LibV5:interburst_20percent_indices #LibV1:interpolate +LibV5:interburst_25percent_indices #LibV5:peak_indices #LibV5:burst_end_indices #LibV5:postburst_fast_ahp_indices #LibV1:interpolate +LibV5:interburst_25percent_values #LibV5:interburst_25percent_indices #LibV1:interpolate +LibV5:interburst_30percent_indices #LibV5:peak_indices #LibV5:burst_end_indices #LibV5:postburst_fast_ahp_indices #LibV1:interpolate +LibV5:interburst_30percent_values #LibV5:interburst_30percent_indices #LibV1:interpolate +LibV5:interburst_40percent_indices #LibV5:peak_indices #LibV5:burst_end_indices #LibV5:postburst_fast_ahp_indices #LibV1:interpolate +LibV5:interburst_40percent_values #LibV5:interburst_40percent_indices #LibV1:interpolate +LibV5:interburst_60percent_indices #LibV5:peak_indices #LibV5:burst_end_indices #LibV5:postburst_fast_ahp_indices #LibV1:interpolate +LibV5:interburst_60percent_values #LibV5:interburst_60percent_indices #LibV1:interpolate +LibV5:interburst_duration #LibV1:peak_time #LibV5:burst_end_indices #LibV1:interpolate diff --git a/efel/cppcore/FillFptrTable.cpp b/efel/cppcore/FillFptrTable.cpp index 6c60463a..333b4713 100644 --- a/efel/cppcore/FillFptrTable.cpp +++ b/efel/cppcore/FillFptrTable.cpp @@ -181,7 +181,77 @@ int FillFptrTable() { FptrTableV5["interburst_min_values"] = &LibV5::interburst_min_values; FptrTableV5["postburst_min_indices"] = &LibV5::postburst_min_indices; FptrTableV5["postburst_min_values"] = &LibV5::postburst_min_values; + FptrTableV5["postburst_slow_ahp_indices"] = &LibV5::postburst_slow_ahp_indices; + FptrTableV5["postburst_slow_ahp_values"] = &LibV5::postburst_slow_ahp_values; FptrTableV5["time_to_interburst_min"] = &LibV5::time_to_interburst_min; + FptrTableV5["time_to_postburst_slow_ahp"] = &LibV5::time_to_postburst_slow_ahp; + FptrTableV5["postburst_fast_ahp_indices"] = &LibV5::postburst_fast_ahp_indices; + FptrTableV5["postburst_fast_ahp_values"] = &LibV5::postburst_fast_ahp_values; + FptrTableV5["postburst_adp_peak_indices"] = &LibV5::postburst_adp_peak_indices; + FptrTableV5["postburst_adp_peak_values"] = &LibV5::postburst_adp_peak_values; + FptrTableV5["time_to_postburst_fast_ahp"] = &LibV5::time_to_postburst_fast_ahp; + FptrTableV5["time_to_postburst_adp_peak"] = &LibV5::time_to_postburst_adp_peak; + FptrTableV5["interburst_15percent_indices"] = [](mapStr2intVec& intData, + mapStr2doubleVec& doubleData, + mapStr2Str& strData) { + return LibV5::interburst_XXpercent_indices(intData, doubleData, strData, 15); + }; + FptrTableV5["interburst_15percent_values"] = [](mapStr2intVec& intData, + mapStr2doubleVec& doubleData, + mapStr2Str& strData) { + return LibV5::interburst_XXpercent_indices(intData, doubleData, strData, 15); + }; + FptrTableV5["interburst_20percent_indices"] = [](mapStr2intVec& intData, + mapStr2doubleVec& doubleData, + mapStr2Str& strData) { + return LibV5::interburst_XXpercent_indices(intData, doubleData, strData, 20); + }; + FptrTableV5["interburst_20percent_values"] = [](mapStr2intVec& intData, + mapStr2doubleVec& doubleData, + mapStr2Str& strData) { + return LibV5::interburst_XXpercent_indices(intData, doubleData, strData, 20); + }; + FptrTableV5["interburst_25percent_indices"] = [](mapStr2intVec& intData, + mapStr2doubleVec& doubleData, + mapStr2Str& strData) { + return LibV5::interburst_XXpercent_indices(intData, doubleData, strData, 25); + }; + FptrTableV5["interburst_25percent_values"] = [](mapStr2intVec& intData, + mapStr2doubleVec& doubleData, + mapStr2Str& strData) { + return LibV5::interburst_XXpercent_indices(intData, doubleData, strData, 25); + }; + FptrTableV5["interburst_30percent_indices"] = [](mapStr2intVec& intData, + mapStr2doubleVec& doubleData, + mapStr2Str& strData) { + return LibV5::interburst_XXpercent_indices(intData, doubleData, strData, 30); + }; + FptrTableV5["interburst_30percent_values"] = [](mapStr2intVec& intData, + mapStr2doubleVec& doubleData, + mapStr2Str& strData) { + return LibV5::interburst_XXpercent_indices(intData, doubleData, strData, 30); + }; + FptrTableV5["interburst_40percent_indices"] = [](mapStr2intVec& intData, + mapStr2doubleVec& doubleData, + mapStr2Str& strData) { + return LibV5::interburst_XXpercent_indices(intData, doubleData, strData, 40); + }; + FptrTableV5["interburst_40percent_values"] = [](mapStr2intVec& intData, + mapStr2doubleVec& doubleData, + mapStr2Str& strData) { + return LibV5::interburst_XXpercent_indices(intData, doubleData, strData, 40); + }; + FptrTableV5["interburst_60percent_indices"] = [](mapStr2intVec& intData, + mapStr2doubleVec& doubleData, + mapStr2Str& strData) { + return LibV5::interburst_XXpercent_indices(intData, doubleData, strData, 60); + }; + FptrTableV5["interburst_60percent_values"] = [](mapStr2intVec& intData, + mapStr2doubleVec& doubleData, + mapStr2Str& strData) { + return LibV5::interburst_XXpercent_indices(intData, doubleData, strData, 60); + }; + FptrTableV5["interburst_duration"] = &LibV5::interburst_duration; //****************** end of FptrTableV5 ***************************** diff --git a/efel/cppcore/LibV5.cpp b/efel/cppcore/LibV5.cpp index c6b4908a..55a80a51 100644 --- a/efel/cppcore/LibV5.cpp +++ b/efel/cppcore/LibV5.cpp @@ -26,6 +26,7 @@ #include #include #include +#include #include "EfelExceptions.h" @@ -2264,21 +2265,27 @@ static int __postburst_min_indices(const vector& t, [stim_end](double x) { return x >= stim_end; })); end_index = distance(t.begin(), t.end()); for (size_t i = 0; i < burst_end_indices.size(); i++) { - if (burst_end_indices[i] + 1 < peak_indices.size()) { - postburst_min_index = - min_element(v.begin() + peak_indices[burst_end_indices[i]], - v.begin() + peak_indices[burst_end_indices[i] + 1]) - - v.begin(); - } else if (peak_indices[burst_end_indices[i]] < stim_end_index) { - postburst_min_index = - min_element(v.begin() + peak_indices[burst_end_indices[i]], - v.begin() + stim_end_index) - - v.begin(); + if (burst_end_indices[i] + 1 < peak_indices.size()){ + postburst_min_index = min_element( + v.begin() + peak_indices[burst_end_indices[i]], + v.begin() + peak_indices[burst_end_indices[i] + 1] + ) - v.begin(); + } else if (peak_indices[burst_end_indices[i]] < stim_end_index){ + postburst_min_index = min_element( + v.begin() + peak_indices[burst_end_indices[i]], + v.begin() + stim_end_index + ) - v.begin(); + if (postburst_min_index == stim_end_index){ + continue; + } } else { - postburst_min_index = - min_element(v.begin() + peak_indices[burst_end_indices[i]], - v.begin() + end_index) - - v.begin(); + postburst_min_index = min_element( + v.begin() + peak_indices[burst_end_indices[i]], + v.begin() + end_index + ) - v.begin(); + if (postburst_min_index == end_index){ + continue; + } } postburst_min_indices.push_back(postburst_min_index); @@ -2317,6 +2324,106 @@ int LibV5::postburst_min_values(mapStr2intVec& IntFeatureData, return 1; } +static int __postburst_slow_ahp_indices(const vector& t, + const vector& v, + const vector& peak_indices, + const vector& burst_end_indices, + vector& postburst_slow_ahp_indices, + vector& postburst_slow_ahp_values, + const double stim_end, + const double sahp_start) { + unsigned postburst_slow_ahp_index, stim_end_index, end_index, t_start_index; + stim_end_index = + distance(t.begin(), + find_if(t.begin(), t.end(), + [stim_end](double x) { return x >= stim_end; })); + end_index = distance(t.begin(), t.end()); + for (size_t i = 0; i < burst_end_indices.size(); i++) { + double t_start = t[peak_indices[burst_end_indices[i]]] + sahp_start; + + if (burst_end_indices[i] + 1 < peak_indices.size()){ + t_start_index = find_if( + t.begin() + peak_indices[burst_end_indices[i]], + t.begin() + peak_indices[burst_end_indices[i] + 1], + [t_start](double x) { return x >= t_start; } + ) - t.begin(); + postburst_slow_ahp_index = min_element( + v.begin() + t_start_index, + v.begin() + peak_indices[burst_end_indices[i] + 1] + ) - v.begin(); + } else if (peak_indices[burst_end_indices[i]] < stim_end_index){ + t_start_index = find_if( + t.begin() + peak_indices[burst_end_indices[i]], + t.begin() + stim_end_index, + [t_start](double x) { return x >= t_start; } + ) - t.begin(); + if (t_start_index < stim_end_index){ + postburst_slow_ahp_index = min_element( + v.begin() + t_start_index, v.begin() + stim_end_index + ) - v.begin(); + } else { + // edge case: stim_end_index is 1 index after stim_end + continue; + } + } else { + t_start_index = find_if( + t.begin() + peak_indices[burst_end_indices[i]], + t.begin() + end_index, + [t_start](double x) { return x >= t_start; } + ) - t.begin(); + if (t_start_index < end_index){ + postburst_slow_ahp_index = min_element( + v.begin() + t_start_index, v.begin() + end_index + ) - v.begin(); + } else{ + // edge case: end_index is 1 index after end + continue; + } + } + + postburst_slow_ahp_indices.push_back(postburst_slow_ahp_index); + postburst_slow_ahp_values.push_back(v[postburst_slow_ahp_index]); + } + + return postburst_slow_ahp_indices.size(); +} + +int LibV5::postburst_slow_ahp_indices(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData) { + const auto& doubleFeatures = + getFeatures(DoubleFeatureData, {"T", "V", "stim_end", "sahp_start"}); + const auto& intFeatures = + getFeatures(IntFeatureData, {"peak_indices", "burst_end_indices"}); + + vector postburst_slow_ahp_values; + vector postburst_slow_ahp_indices; + int retVal = __postburst_slow_ahp_indices( + doubleFeatures.at("T"), + doubleFeatures.at("V"), + intFeatures.at("peak_indices"), + intFeatures.at("burst_end_indices"), + postburst_slow_ahp_indices, + postburst_slow_ahp_values, + doubleFeatures.at("stim_end").front(), + // time after the spike in ms after which to start searching for minimum + doubleFeatures.at("sahp_start").empty() ? 5.0 : doubleFeatures.at("sahp_start").front() + ); + if (retVal >= 0) { + setVec(IntFeatureData, StringData, "postburst_slow_ahp_indices", + postburst_slow_ahp_indices); + setVec(DoubleFeatureData, StringData, "postburst_slow_ahp_values", + postburst_slow_ahp_values); + } + return retVal; +} + +int LibV5::postburst_slow_ahp_values(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData) { + return 1; +} + int LibV5::time_to_interburst_min(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData) { @@ -2346,3 +2453,346 @@ int LibV5::time_to_interburst_min(mapStr2intVec& IntFeatureData, time_to_interburst_min); return time_to_interburst_min.size(); } + +int LibV5::time_to_postburst_slow_ahp(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData) { + const auto& doubleFeatures = + getFeatures(DoubleFeatureData, {"T", "peak_time"}); + const auto& intFeatures = + getFeatures(IntFeatureData, {"postburst_slow_ahp_indices", "burst_end_indices"}); + + vector time_to_postburst_slow_ahp; + const vector& time = doubleFeatures.at("T"); + const vector& peak_time = doubleFeatures.at("peak_time"); + const vector& burst_end_indices = intFeatures.at("burst_end_indices"); + const vector& postburst_slow_ahp_indices = intFeatures.at("postburst_slow_ahp_indices"); + + if (burst_end_indices.size() < postburst_slow_ahp_indices.size()){ + GErrorStr += + "\nburst_end_indices should not have less elements than postburst_slow_ahp_indices\n"; + return -1; + } + + for (size_t i = 0; i < postburst_slow_ahp_indices.size(); i++) { + time_to_postburst_slow_ahp.push_back(time[postburst_slow_ahp_indices[i]] - + peak_time[burst_end_indices[i]]); + } + setVec(DoubleFeatureData, StringData, "time_to_postburst_slow_ahp", + time_to_postburst_slow_ahp); + return (time_to_postburst_slow_ahp.size()); +} + +static int __postburst_fast_ahp_indices(const vector& t, const vector& v, + const vector& peak_indices, + const vector& burst_end_indices, + const double stim_end, + vector& postburst_fast_ahp_indices, + vector& postburst_fast_ahp_values) { + vector start_indices, end_indices; + for (size_t i = 0; i < burst_end_indices.size(); i++) { + start_indices.push_back(peak_indices[burst_end_indices[i]]); + if (burst_end_indices[i] + 1 < peak_indices.size()){ + end_indices.push_back(peak_indices[burst_end_indices[i] + 1]); + } + } + + unsigned end_index = 0; + if (t[start_indices.back()] < stim_end) { + end_index = + distance(t.begin(), + find_if(t.begin(), t.end(), + [stim_end](double x) { return x >= stim_end; })); + } else { + end_index = distance(t.begin(), t.end()); + } + + if (end_indices.size() < start_indices.size()){ + end_indices.push_back(end_index); + } + + size_t fahpindex = 0; + for (size_t i = 0; i < start_indices.size(); i++) { + // can use first_min_element because dv/dt is very steep before fash ahp + // and noise is very unlikely to make voltage go up before reaching fast ahp + fahpindex = distance( + v.begin(), first_min_element(v.begin() + start_indices[i], + v.begin() + end_indices[i])); + + if (fahpindex != end_index - 1) { + postburst_fast_ahp_indices.push_back(fahpindex); + + EFEL_ASSERT(fahpindex < v.size(), + "fast AHP index falls outside of voltage array"); + postburst_fast_ahp_values.push_back(v[fahpindex]); + } + } + + return postburst_fast_ahp_indices.size(); +} + +int LibV5::postburst_fast_ahp_indices(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData) { + const auto& doubleFeatures = + getFeatures(DoubleFeatureData, {"T", "V", "stim_end"}); + const auto& intFeatures = + getFeatures(IntFeatureData, {"peak_indices", "burst_end_indices"}); + + vector postburst_fast_ahp_indices; + vector postburst_fast_ahp_values; + int retVal = + __postburst_fast_ahp_indices( + doubleFeatures.at("T"), + doubleFeatures.at("V"), + intFeatures.at("peak_indices"), + intFeatures.at("burst_end_indices"), + doubleFeatures.at("stim_end").front(), + postburst_fast_ahp_indices, + postburst_fast_ahp_values + ); + + if (retVal > 0) { + setVec(IntFeatureData, StringData, "postburst_fast_ahp_indices", postburst_fast_ahp_indices); + setVec(DoubleFeatureData, StringData, "postburst_fast_ahp_values", + postburst_fast_ahp_values); + return postburst_fast_ahp_indices.size(); + } + return -1; +} + +int LibV5::postburst_fast_ahp_values(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData) { + return 1; +} + +static int __postburst_adp_peak_indices(const vector& t, const vector& v, + const vector& postburst_fast_ahp_indices, + const vector& postburst_slow_ahp_indices, + vector& postburst_adp_peak_indices, + vector& postburst_adp_peak_values) { + + if (postburst_slow_ahp_indices.size() > postburst_fast_ahp_indices.size()){ + GErrorStr += + "\n postburst_slow_ahp should not have more elements than " + "postburst_fast_ahp for postburst_adp_peak_indices calculation.\n"; + return -1; + } + size_t adppeakindex = 0; + for (size_t i = 0; i < postburst_slow_ahp_indices.size(); i++) { + if (postburst_slow_ahp_indices[i] < postburst_fast_ahp_indices[i]){ + continue; + } + adppeakindex = distance( + v.begin(), max_element(v.begin() + postburst_fast_ahp_indices[i], + v.begin() + postburst_slow_ahp_indices[i])); + + if (adppeakindex < postburst_slow_ahp_indices[i] - 1) { + postburst_adp_peak_indices.push_back(adppeakindex); + + EFEL_ASSERT(adppeakindex < v.size(), + "ADP peak index falls outside of voltage array"); + postburst_adp_peak_values.push_back(v[adppeakindex]); + } + } + + return postburst_adp_peak_indices.size(); +} + +int LibV5::postburst_adp_peak_indices(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData) { + const auto& doubleFeatures = + getFeatures(DoubleFeatureData, {"T", "V"}); + const auto& intFeatures = + getFeatures(IntFeatureData, {"postburst_fast_ahp_indices", "postburst_slow_ahp_indices"}); + vector postburst_adp_peak_indices; + vector postburst_adp_peak_values; + int retVal = + __postburst_adp_peak_indices( + doubleFeatures.at("T"), + doubleFeatures.at("V"), + intFeatures.at("postburst_fast_ahp_indices"), + intFeatures.at("postburst_slow_ahp_indices"), + postburst_adp_peak_indices, + postburst_adp_peak_values + ); + + if (retVal > 0) { + setVec(IntFeatureData, StringData, "postburst_adp_peak_indices", postburst_adp_peak_indices); + setVec(DoubleFeatureData, StringData, "postburst_adp_peak_values", + postburst_adp_peak_values); + return retVal; + } + return -1; +} + +int LibV5::postburst_adp_peak_values(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData) { + return 1; +} + +int LibV5::time_to_postburst_fast_ahp(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData) { + const auto& doubleFeatures = + getFeatures(DoubleFeatureData, {"T", "peak_time"}); + const auto& intFeatures = + getFeatures(IntFeatureData, {"postburst_fast_ahp_indices", "burst_end_indices"}); + + vector time_to_postburst_fast_ahp; + const vector& time = doubleFeatures.at("T"); + const vector& peak_time = doubleFeatures.at("peak_time"); + const vector& burst_end_indices = intFeatures.at("burst_end_indices"); + const vector& postburst_fast_ahp_indices = intFeatures.at("postburst_fast_ahp_indices"); + + if (burst_end_indices.size() < postburst_fast_ahp_indices.size()){ + GErrorStr += + "\nburst_end_indices should not have less elements than postburst_fast_ahp_indices\n"; + return -1; + } + + for (size_t i = 0; i < postburst_fast_ahp_indices.size(); i++) { + time_to_postburst_fast_ahp.push_back(time[postburst_fast_ahp_indices[i]] - + peak_time[burst_end_indices[i]]); + } + setVec(DoubleFeatureData, StringData, "time_to_postburst_fast_ahp", + time_to_postburst_fast_ahp); + return (time_to_postburst_fast_ahp.size()); +} + +int LibV5::time_to_postburst_adp_peak(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData) { + const auto& doubleFeatures = + getFeatures(DoubleFeatureData, {"T", "peak_time"}); + const auto& intFeatures = + getFeatures(IntFeatureData, {"postburst_adp_peak_indices", "burst_end_indices"}); + + vector time_to_postburst_adp_peak; + const vector& time = doubleFeatures.at("T"); + const vector& peak_time = doubleFeatures.at("peak_time"); + const vector& burst_end_indices = intFeatures.at("burst_end_indices"); + const vector& postburst_adp_peak_indices = intFeatures.at("postburst_adp_peak_indices"); + + if (burst_end_indices.size() < postburst_adp_peak_indices.size()){ + GErrorStr += + "\nburst_end_indices should not have less elements than postburst_adp_peak_indices\n"; + return -1; + } + + for (size_t i = 0; i < postburst_adp_peak_indices.size(); i++) { + // there are not always an adp peak after each burst + // so make sure that the burst and adp peak indices are consistent + size_t k = 0; + while (burst_end_indices[i] + k + 1 < peak_time.size() && + peak_time[burst_end_indices[i] + k + 1] < time[postburst_adp_peak_indices[i]]){ + k++; + } + time_to_postburst_adp_peak.push_back(time[postburst_adp_peak_indices[i]] - + peak_time[burst_end_indices[i] + k]); + } + setVec(DoubleFeatureData, StringData, "time_to_postburst_adp_peak", + time_to_postburst_adp_peak); + return (time_to_postburst_adp_peak.size()); +} + +// index and voltage value at a given percentage of the duration of the interburst after fast AHP +int __interburst_percent_indices(const vector& t, const vector& v, + const vector& postburst_fast_ahp_indices, + const vector& peak_indices, + const vector& burst_end_indices, + vector& interburst_percent_indices, + vector& interburst_percent_values, + // percentage should be a value between 0 and 1 + double fraction) { + + double time_interval, time_at_fraction; + size_t index_at_fraction; + for (size_t i = 0; i < postburst_fast_ahp_indices.size(); i++) { + if (i < burst_end_indices.size()){ + if (burst_end_indices[i] + 1 < peak_indices.size()){ + time_interval = t[peak_indices[burst_end_indices[i] + 1]] - t[postburst_fast_ahp_indices[i]]; + time_at_fraction = t[postburst_fast_ahp_indices[i]] + time_interval * fraction; + index_at_fraction = + distance(t.begin(), + find_if(t.begin(), t.end(), + [time_at_fraction](double x){ return x >= time_at_fraction; })); + interburst_percent_indices.push_back(index_at_fraction); + interburst_percent_values.push_back(v[index_at_fraction]); + } + } + } + return interburst_percent_indices.size(); +} + +int LibV5::interburst_XXpercent_indices(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData, int percent) { + const auto& doubleFeatures = + getFeatures(DoubleFeatureData, {"T", "V"}); + const auto& intFeatures = + getFeatures(IntFeatureData, {"peak_indices", "burst_end_indices", "postburst_fast_ahp_indices"}); + + vector interburst_XXpercent_indices; + vector interburst_XXpercent_values; + char featureNameIndices[30], featureNameValues[30]; + sprintf(featureNameIndices, "interburst_%dpercent_indices", percent); + sprintf(featureNameValues, "interburst_%dpercent_values", percent); + + int retVal = + __interburst_percent_indices( + doubleFeatures.at("T"), + doubleFeatures.at("V"), + intFeatures.at("postburst_fast_ahp_indices"), + intFeatures.at("peak_indices"), + intFeatures.at("burst_end_indices"), + interburst_XXpercent_indices, + interburst_XXpercent_values, + percent / 100. + ); + + if (retVal > 0) { + setVec(IntFeatureData, StringData, featureNameIndices, interburst_XXpercent_indices); + setVec(DoubleFeatureData, StringData, featureNameValues, + interburst_XXpercent_values); + return interburst_XXpercent_indices.size(); + } + return -1; +} + +static int __interburst_duration(const vector& peak_time, + const vector& burst_end_indices, + vector& interburst_duration) { + + double duration; + for (size_t i = 0; i < burst_end_indices.size(); i++) { + if (burst_end_indices[i] + 1 < peak_time.size()){ + duration = peak_time[burst_end_indices[i] + 1] - peak_time[burst_end_indices[i]]; + interburst_duration.push_back(duration); + } + } + return interburst_duration.size(); +} + +int LibV5::interburst_duration(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData) { + const auto& doubleFeatures = + getFeatures(DoubleFeatureData, {"peak_time"}); + const auto& intFeatures = + getFeatures(IntFeatureData, {"burst_end_indices"}); + + vector interburst_duration; + int retVal = + __interburst_duration( + doubleFeatures.at("peak_time"), intFeatures.at("burst_end_indices"), interburst_duration); + + if (retVal > 0) { + setVec(DoubleFeatureData, StringData, "interburst_duration", interburst_duration); + return interburst_duration.size(); + } + return -1; +} diff --git a/efel/cppcore/LibV5.h b/efel/cppcore/LibV5.h index 95c4d019..f4868045 100644 --- a/efel/cppcore/LibV5.h +++ b/efel/cppcore/LibV5.h @@ -259,8 +259,42 @@ int postburst_min_indices(mapStr2intVec& IntFeatureData, int postburst_min_values(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); +int postburst_slow_ahp_indices(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData); +int postburst_slow_ahp_values(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData); int time_to_interburst_min(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); +int time_to_postburst_slow_ahp(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData); +int postburst_fast_ahp_indices(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData); +int postburst_fast_ahp_values(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData); +int postburst_adp_peak_indices(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData); +int postburst_adp_peak_values(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData); +int time_to_postburst_fast_ahp(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData); +int time_to_postburst_adp_peak(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData); +int interburst_XXpercent_indices(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData, + int percent); +int interburst_duration(mapStr2intVec& IntFeatureData, + mapStr2doubleVec& DoubleFeatureData, + mapStr2Str& StringData); } #endif diff --git a/efel/cppcore/cfeature.cpp b/efel/cppcore/cfeature.cpp index f2b2c956..e2d9b125 100644 --- a/efel/cppcore/cfeature.cpp +++ b/efel/cppcore/cfeature.cpp @@ -204,7 +204,29 @@ void cFeature::fillfeaturetypes() { featuretypes["interburst_min_values"] = "double"; featuretypes["postburst_min_indices"] = "int"; featuretypes["postburst_min_values"] = "double"; + featuretypes["postburst_slow_ahp_indices"] = "int"; + featuretypes["postburst_slow_ahp_values"] = "double"; featuretypes["time_to_interburst_min"] = "double"; + featuretypes["time_to_postburst_slow_ahp"] = "double"; + featuretypes["postburst_fast_ahp_indices"] = "int"; + featuretypes["postburst_fast_ahp_values"] = "double"; + featuretypes["postburst_adp_peak_indices"] = "int"; + featuretypes["postburst_adp_peak_values"] = "double"; + featuretypes["time_to_postburst_fast_ahp"] = "double"; + featuretypes["time_to_postburst_adp_peak"] = "double"; + featuretypes["interburst_15percent_indices"] = "int"; + featuretypes["interburst_15percent_values"] = "double"; + featuretypes["interburst_20percent_indices"] = "int"; + featuretypes["interburst_20percent_values"] = "double"; + featuretypes["interburst_25percent_indices"] = "int"; + featuretypes["interburst_25percent_values"] = "double"; + featuretypes["interburst_30percent_indices"] = "int"; + featuretypes["interburst_30percent_values"] = "double"; + featuretypes["interburst_40percent_indices"] = "int"; + featuretypes["interburst_40percent_values"] = "double"; + featuretypes["interburst_60percent_indices"] = "int"; + featuretypes["interburst_60percent_values"] = "double"; + featuretypes["interburst_duration"] = "double"; // end of feature types } diff --git a/efel/units/units.json b/efel/units/units.json index c80ca33d..bf8d74b7 100644 --- a/efel/units/units.json +++ b/efel/units/units.json @@ -138,5 +138,18 @@ "initburst_sahp_vb": "mV", "depol_block_bool": "constant", "impedance": "Hz", + "postburst_slow_ahp_values": "mV", + "time_to_postburst_slow_ahp": "ms", + "postburst_fast_ahp_values": "mV", + "postburst_adp_peak_values": "mV", + "time_to_postburst_fast_ahp": "ms", + "time_to_postburst_adp_peak": "ms", + "interburst_15percent_values": "mV", + "interburst_20percent_values": "mV", + "interburst_25percent_values": "mV", + "interburst_30percent_values": "mV", + "interburst_40percent_values": "mV", + "interburst_60percent_values": "mV", + "interburst_duration": "ms", "phaseslope_max": "V/s" } diff --git a/tests/featurenames.json b/tests/featurenames.json index 47706097..f97a0514 100644 --- a/tests/featurenames.json +++ b/tests/featurenames.json @@ -147,7 +147,30 @@ "interburst_min_values", "postburst_min_indices", "postburst_min_values", + "postburst_slow_ahp_indices", + "postburst_slow_ahp_values", "time_to_interburst_min", + "time_to_postburst_slow_ahp", + "AHP_depth_slow", + "postburst_fast_ahp_indices", + "postburst_fast_ahp_values", + "postburst_adp_peak_indices", + "postburst_adp_peak_values", + "time_to_postburst_fast_ahp", + "time_to_postburst_adp_peak", + "interburst_15percent_indices", + "interburst_15percent_values", + "interburst_20percent_indices", + "interburst_20percent_values", + "interburst_25percent_indices", + "interburst_25percent_values", + "interburst_30percent_indices", + "interburst_30percent_values", + "interburst_40percent_indices", + "interburst_40percent_values", + "interburst_60percent_indices", + "interburst_60percent_values", + "interburst_duration", "AHP_depth_slow", "phaseslope_max" ] \ No newline at end of file diff --git a/tests/test_basic.py b/tests/test_basic.py index 7339c2dc..2948c0b6 100644 --- a/tests/test_basic.py +++ b/tests/test_basic.py @@ -65,6 +65,8 @@ ) spontaneous_url = testdata_dir / 'basic' / 'spontaneous.txt' +testdata_url = testdata_dir / 'allfeatures' / 'testdata.txt' + def load_data(data_name, interp=False, interp_dt=0.1): """Load data file""" @@ -3554,7 +3556,7 @@ def py_postburst_min_values(t, v, peak_indices, burst_end_indices, stim_end): def test_postburst_min_values(): - """basic: Test interburst_min_values""" + """basic: Test postburst_min_values""" urls = [burst1_url, burst2_url, burst3_url] for i, url in enumerate(urls): import efel @@ -3692,6 +3694,668 @@ def test_spikes_in_burst1_burstlast_diff(): assert list(spikes_in_burst1_burstlast_diff) == [1] +def test_postburst_slow_ahp_values(): + """basic: Test postburst_slow_ahp_values when no fast AHP is present""" + urls = [burst1_url, burst2_url, burst3_url] + for i, url in enumerate(urls): + import efel + efel.reset() + # use this to have all spikes in burst for burst3_url case + efel.setDoubleSetting('strict_burst_factor', 4.0) + + time, voltage = load_ascii_input(url) + + interp_time, interp_voltage = interpolate(time, voltage, 0.1) + + trace = {} + + trace['T'] = time + trace['V'] = voltage + if i in [0, 1]: + trace['stim_start'] = [250] + trace['stim_end'] = [1600] + elif i == 2: + trace['stim_start'] = [800] + trace['stim_end'] = [2150] + + features = [ + "postburst_slow_ahp_values", + "postburst_min_values", + ] + + feature_values = efel.getFeatureValues( + [trace], + features, + raise_warnings=False + ) + + postburst_slow_ahp_values = feature_values[0][ + "postburst_slow_ahp_values" + ] + postburst_min_values = feature_values[0]["postburst_min_values"] + + # convert to float so that None edge case get converted to nan + # and can pass in assert_allclose + postburst_slow_ahp_values = numpy.array( + postburst_slow_ahp_values, dtype=numpy.float64 + ) + postburst_min_values = numpy.array( + postburst_min_values, dtype=numpy.float64 + ) + + # are equal when there is no fast AHP + numpy.testing.assert_allclose( + postburst_slow_ahp_values, postburst_min_values + ) + + +def test_time_to_postburst_slow_ahp(): + """basic: Test time_to_postburst_slow_ahp""" + urls = [burst1_url, burst2_url, burst3_url] + for i, url in enumerate(urls): + import efel + efel.reset() + + time, voltage = load_ascii_input(url) + + interp_time, interp_voltage = interpolate(time, voltage, 0.1) + + trace = {} + + trace['T'] = time + trace['V'] = voltage + if i in [0, 1]: + trace['stim_start'] = [250] + trace['stim_end'] = [1600] + elif i == 2: + trace['stim_start'] = [800] + trace['stim_end'] = [2150] + + features = [ + "burst_end_indices", + "peak_time", + "time_to_postburst_slow_ahp", + "postburst_slow_ahp_indices", + ] + + feature_values = efel.getFeatureValues( + [trace], + features, + raise_warnings=False + ) + + peak_time = feature_values[0]["peak_time"] + burst_end_indices = feature_values[0]["burst_end_indices"] + postburst_slow_ahp_indices = feature_values[0][ + "postburst_slow_ahp_indices" + ] + time_to_postburst_slow_ahp = feature_values[0][ + "time_to_postburst_slow_ahp" + ] + + if ( + postburst_slow_ahp_indices is None or len(postburst_slow_ahp_indices) == 0 + or burst_end_indices is None or len(burst_end_indices) == 0 + ): + time_to_postburst_slow_ahp_py = None + else: + time_to_postburst_slow_ahp_py = interp_time[ + postburst_slow_ahp_indices + ] - peak_time[burst_end_indices] + + # convert to float so that None edge case get converted to nan + # and can pass in assert_allclose + time_to_postburst_slow_ahp_py = numpy.array( + time_to_postburst_slow_ahp_py, dtype=numpy.float64 + ) + time_to_postburst_slow_ahp = numpy.array( + time_to_postburst_slow_ahp, dtype=numpy.float64 + ) + + numpy.testing.assert_allclose( + time_to_postburst_slow_ahp_py, time_to_postburst_slow_ahp + ) + + +def py_postburst_fahp(v, peak_indices, burst_end_indices, stim_end_index): + """python implementation of post burst fast ahp.""" + if burst_end_indices is None: + return None, None + + postburst_fahp = [] + postburst_fahp_idxs = [] + for i in burst_end_indices: + if i + 1 < len(peak_indices): + stop_i = peak_indices[i + 1] + elif i + 1 < stim_end_index: + stop_i = stim_end_index + else: + stop_i = len(v) - 1 + + v_crop = v[peak_indices[i]:stop_i] + # get where the voltage is going up + crop_args = numpy.argwhere(numpy.diff(v_crop) >= 0)[:, 0] + # the voltage should go up for at least two consecutive points + crop_arg_arg = numpy.argwhere(numpy.diff(crop_args) == 1)[0][0] + crop_arg = crop_args[crop_arg_arg] + end_i = peak_indices[i] + crop_arg + 1 + # the fast ahp is between the last peak of the burst and + # the point where the voltage is going back up + postburst_fahp.append(numpy.min(v[peak_indices[i]:end_i])) + postburst_fahp_idxs.append( + numpy.argmin(v[peak_indices[i]:end_i]) + peak_indices[i] + ) + + return postburst_fahp_idxs, postburst_fahp + + +def test_postburst_fast_ahp_values(): + """basic: Test postburst_fast_ahp_values""" + urls = [burst1_url, burst2_url, burst3_url, testdata_url] + for i, url in enumerate(urls): + import efel + efel.reset() + + time, voltage = load_ascii_input(url) + + interp_time, interp_voltage = interpolate(time, voltage, 0.1) + + trace = {} + + trace['T'] = time + trace['V'] = voltage + if i in [0, 1]: + trace['stim_start'] = [250] + trace['stim_end'] = [1600] + elif i == 2: + trace['stim_start'] = [800] + trace['stim_end'] = [2150] + elif i == 3: + trace['stim_start'] = [700] + trace['stim_end'] = [2700] + + features = [ + "burst_end_indices", + "peak_indices", + "postburst_fast_ahp_indices", + "postburst_fast_ahp_values", + ] + + feature_values = efel.getFeatureValues( + [trace], + features, + raise_warnings=False + ) + + peak_indices = feature_values[0]["peak_indices"] + burst_end_indices = feature_values[0]["burst_end_indices"] + postburst_fahpi = feature_values[0][ + "postburst_fast_ahp_indices" + ] + postburst_fahp = feature_values[0][ + "postburst_fast_ahp_values" + ] + + stim_end_index = numpy.argwhere(trace['stim_end'] > interp_time)[0][0] + postburst_fahpi_py, postburst_fahp_py = py_postburst_fahp( + interp_voltage, peak_indices, burst_end_indices, stim_end_index + ) + + postburst_fahp = numpy.array( + postburst_fahp, dtype=numpy.float64 + ) + postburst_fahp_py = numpy.array( + postburst_fahp_py, dtype=numpy.float64 + ) + numpy.testing.assert_allclose( + postburst_fahp, postburst_fahp_py + ) + + postburst_fahpi = numpy.array( + postburst_fahpi, dtype=numpy.float64 + ) + postburst_fahpi_py = numpy.array( + postburst_fahpi_py, dtype=numpy.float64 + ) + numpy.testing.assert_allclose( + postburst_fahpi, postburst_fahpi_py + ) + + +def py_postburst_adppeak( + v, postburst_sahpi, postburst_fahpi, +): + """python implementation of adp peak after burst.""" + if postburst_sahpi is None or postburst_fahpi is None: + return None, None + + adp_peak_indices = [] + adp_peak_values = [] + for i, sahpi in enumerate(postburst_sahpi): + if sahpi < postburst_fahpi[i]: + continue + adppeaki = numpy.argmax( + v[postburst_fahpi[i]:sahpi] + ) + postburst_fahpi[i] + if adppeaki != sahpi - 1: + adp_peak_indices.append(adppeaki) + adp_peak_values.append(v[adppeaki]) + + if len(adp_peak_indices) == 0: + return None, None + return adp_peak_indices, adp_peak_values + + +def test_postburst_adp_peak_values(): + """basic: Test postburst_adp_peak_values""" + urls = [burst1_url, burst2_url, burst3_url, testdata_url] + for i, url in enumerate(urls): + import efel + efel.reset() + + time, voltage = load_ascii_input(url) + + interp_time, interp_voltage = interpolate(time, voltage, 0.1) + + trace = {} + + trace['T'] = time + trace['V'] = voltage + if i in [0, 1]: + trace['stim_start'] = [250] + trace['stim_end'] = [1600] + elif i == 2: + trace['stim_start'] = [800] + trace['stim_end'] = [2150] + elif i == 3: + trace['stim_start'] = [700] + trace['stim_end'] = [2700] + + features = [ + "postburst_fast_ahp_indices", + "postburst_slow_ahp_indices", + "postburst_adp_peak_indices", + "postburst_adp_peak_values", + ] + + feature_values = efel.getFeatureValues( + [trace], + features, + raise_warnings=False + ) + + postburst_fahpi = feature_values[0][ + "postburst_fast_ahp_indices" + ] + postburst_sahpi = feature_values[0][ + "postburst_slow_ahp_indices" + ] + postburst_adppeaki = feature_values[0][ + "postburst_adp_peak_indices" + ] + postburst_adppeak = feature_values[0][ + "postburst_adp_peak_values" + ] + + postburst_adppeaki_py, postburst_adppeak_py = py_postburst_adppeak( + interp_voltage, postburst_sahpi, postburst_fahpi, + ) + + postburst_adppeak = numpy.array( + postburst_adppeak, dtype=numpy.float64 + ) + postburst_adppeak_py = numpy.array( + postburst_adppeak_py, dtype=numpy.float64 + ) + numpy.testing.assert_allclose( + postburst_adppeak, postburst_adppeak_py + ) + + postburst_adppeaki = numpy.array( + postburst_adppeaki, dtype=numpy.float64 + ) + postburst_adppeaki_py = numpy.array( + postburst_adppeaki_py, dtype=numpy.float64 + ) + numpy.testing.assert_allclose( + postburst_adppeaki, postburst_adppeaki_py + ) + + +def py_time_to_postburst_fast_ahp(t, postburst_fahpi, burst_endi, peak_time): + """python implementation of time to post burst fast ahp.""" + if postburst_fahpi is None or burst_endi is None: + return None + return [t[fahpi] - peak_time[burst_endi[i]] for (i, fahpi) in + enumerate(postburst_fahpi)] + + +def test_time_to_postburst_fast_ahp(): + """basic: Test time_to_postburst_fast_ahp""" + urls = [burst1_url, burst2_url, burst3_url, testdata_url] + for i, url in enumerate(urls): + import efel + efel.reset() + + time, voltage = load_ascii_input(url) + + interp_time, interp_voltage = interpolate(time, voltage, 0.1) + + trace = {} + + trace['T'] = time + trace['V'] = voltage + if i in [0, 1]: + trace['stim_start'] = [250] + trace['stim_end'] = [1600] + elif i == 2: + trace['stim_start'] = [800] + trace['stim_end'] = [2150] + elif i == 3: + trace['stim_start'] = [700] + trace['stim_end'] = [2700] + + features = [ + "time_to_postburst_fast_ahp", + "postburst_fast_ahp_indices", + "burst_end_indices", + "peak_time", + ] + + feature_values = efel.getFeatureValues( + [trace], + features, + raise_warnings=False + ) + + time_to_postburst_fast_ahp = feature_values[0][ + "time_to_postburst_fast_ahp" + ] + postburst_fahpi = feature_values[0][ + "postburst_fast_ahp_indices" + ] + burst_endi = feature_values[0][ + "burst_end_indices" + ] + peak_time = feature_values[0][ + "peak_time" + ] + + time_to_postburst_fast_ahp_py = py_time_to_postburst_fast_ahp( + interp_time, postburst_fahpi, burst_endi, peak_time + ) + + time_to_postburst_fast_ahp = numpy.array( + time_to_postburst_fast_ahp, dtype=numpy.float64 + ) + time_to_postburst_fast_ahp_py = numpy.array( + time_to_postburst_fast_ahp_py, dtype=numpy.float64 + ) + numpy.testing.assert_allclose( + time_to_postburst_fast_ahp, time_to_postburst_fast_ahp_py + ) + + +def py_time_to_postburst_adp_peak( + t, postburst_adppeaki, burst_endi, peak_time +): + """python implementation of time to post burst ADP peak.""" + if postburst_adppeaki is None or burst_endi is None: + return None + + time_to_postburst_adp_peaks = [] + n_peaks = len(peak_time) + for i, adppeaki in enumerate(postburst_adppeaki): + # there are not always an adp peak after each burst + # so make sure that the burst and adp peak indices are consistent + k = 0 + while (burst_endi[i] + k + 1 < n_peaks and + peak_time[burst_endi[i] + k + 1] < t[adppeaki]): + k += 1 + + time_to_postburst_adp_peaks.append( + t[adppeaki] - peak_time[burst_endi[i] + k] + ) + + return time_to_postburst_adp_peaks + + +def test_time_to_postburst_adp_peak(): + """basic: Test time_to_postburst_adp_peak""" + urls = [burst1_url, burst2_url, burst3_url, testdata_url] + for i, url in enumerate(urls): + import efel + efel.reset() + + time, voltage = load_ascii_input(url) + + interp_time, interp_voltage = interpolate(time, voltage, 0.1) + + trace = {} + + trace['T'] = time + trace['V'] = voltage + if i in [0, 1]: + trace['stim_start'] = [250] + trace['stim_end'] = [1600] + elif i == 2: + trace['stim_start'] = [800] + trace['stim_end'] = [2150] + elif i == 3: + trace['stim_start'] = [700] + trace['stim_end'] = [2700] + + features = [ + "time_to_postburst_adp_peak", + "postburst_adp_peak_indices", + "burst_end_indices", + "peak_time", + ] + + feature_values = efel.getFeatureValues( + [trace], + features, + raise_warnings=False + ) + + time_to_postburst_adp_peak = feature_values[0][ + "time_to_postburst_adp_peak" + ] + postburst_adppeaki = feature_values[0][ + "postburst_adp_peak_indices" + ] + burst_endi = feature_values[0][ + "burst_end_indices" + ] + peak_time = feature_values[0][ + "peak_time" + ] + + time_to_postburst_adp_peak_py = py_time_to_postburst_adp_peak( + interp_time, postburst_adppeaki, burst_endi, peak_time + ) + + time_to_postburst_adp_peak = numpy.array( + time_to_postburst_adp_peak, dtype=numpy.float64 + ) + time_to_postburst_adp_peak_py = numpy.array( + time_to_postburst_adp_peak_py, dtype=numpy.float64 + ) + numpy.testing.assert_allclose( + time_to_postburst_adp_peak, time_to_postburst_adp_peak_py + ) + + +def py_XXpercent_interburst_values( + t, v, postburst_fahpi, burst_endi, peaki, percent=20. +): + """python implementation of 20% / 25% / 30% interburst values.""" + if postburst_fahpi is None or burst_endi is None: + return None, None + + XXpercent_interburst = [] + XXpercent_interburst_i = [] + for i, postburst_fahp_i in enumerate(postburst_fahpi): + if i < len(burst_endi) and burst_endi[i] + 1 < len(peaki): + time_interval = t[peaki[burst_endi[i] + 1]] - t[postburst_fahp_i] + time_at_XXpercent = (t[postburst_fahp_i] + + time_interval * percent / 100.) + index_at_XXpercent = numpy.argwhere(t >= time_at_XXpercent)[0][0] + XXpercent_interburst_i.append(index_at_XXpercent) + XXpercent_interburst.append(v[index_at_XXpercent]) + + return XXpercent_interburst_i, XXpercent_interburst + + +def interburst_XXpercent_values_compare_with_py( + feature_values, interp_time, interp_voltage, percent +): + """compare python implementation with efel output for interburst values""" + postburst_fahpi = feature_values[0]["postburst_fast_ahp_indices"] + burst_endi = feature_values[0]["burst_end_indices"] + peaki = feature_values[0]["peak_indices"] + interbursti = feature_values[0][f"interburst_{percent}percent_indices"] + interburst = feature_values[0][f"interburst_{percent}percent_values"] + + interbursti_py, interburst_py = py_XXpercent_interburst_values( + interp_time, + interp_voltage, + postburst_fahpi, + burst_endi, + peaki, + percent, + ) + + interburst = numpy.array(interburst, dtype=numpy.float64) + interburst_py = numpy.array(interburst_py, dtype=numpy.float64) + numpy.testing.assert_allclose(interburst, interburst_py) + + interbursti = numpy.array(interbursti, dtype=numpy.float64) + interbursti_py = numpy.array(interbursti_py, dtype=numpy.float64) + numpy.testing.assert_allclose(interbursti, interbursti_py) + + +def test_interburst_XXpercent_values(): + """basic: Test interburst_15percent_values, interburst_20percent_values, + interburst_25percent_values, interburst_30percent_values, + interburst_40percent_values, interburst_60percent_values + """ + urls = [burst1_url, burst2_url, burst3_url, testdata_url] + for i, url in enumerate(urls): + import efel + efel.reset() + + time, voltage = load_ascii_input(url) + + interp_time, interp_voltage = interpolate(time, voltage, 0.1) + + trace = {} + + trace['T'] = time + trace['V'] = voltage + if i in [0, 1]: + trace['stim_start'] = [250] + trace['stim_end'] = [1600] + elif i == 2: + trace['stim_start'] = [800] + trace['stim_end'] = [2150] + elif i == 3: + trace['stim_start'] = [700] + trace['stim_end'] = [2700] + + features = [ + "postburst_fast_ahp_indices", + "burst_end_indices", + "peak_indices", + "interburst_15percent_indices", + "interburst_15percent_values", + "interburst_20percent_indices", + "interburst_20percent_values", + "interburst_25percent_indices", + "interburst_25percent_values", + "interburst_30percent_indices", + "interburst_30percent_values", + "interburst_40percent_indices", + "interburst_40percent_values", + "interburst_60percent_indices", + "interburst_60percent_values", + ] + + feature_values = efel.getFeatureValues( + [trace], + features, + raise_warnings=False + ) + + interburst_XXpercent_values_compare_with_py( + feature_values, interp_time, interp_voltage, 15 + ) + interburst_XXpercent_values_compare_with_py( + feature_values, interp_time, interp_voltage, 20 + ) + interburst_XXpercent_values_compare_with_py( + feature_values, interp_time, interp_voltage, 25 + ) + interburst_XXpercent_values_compare_with_py( + feature_values, interp_time, interp_voltage, 30 + ) + interburst_XXpercent_values_compare_with_py( + feature_values, interp_time, interp_voltage, 40 + ) + interburst_XXpercent_values_compare_with_py( + feature_values, interp_time, interp_voltage, 60 + ) + + +def test_interburst_duration(): + """basic: Test interburst_duration""" + urls = [burst1_url, burst2_url, burst3_url, testdata_url] + for i, url in enumerate(urls): + import efel + efel.reset() + + time, voltage = load_ascii_input(url) + + trace = {} + trace['T'] = time + trace['V'] = voltage + if i in [0, 1]: + trace['stim_start'] = [250] + trace['stim_end'] = [1600] + elif i == 2: + trace['stim_start'] = [800] + trace['stim_end'] = [2150] + elif i == 3: + trace['stim_start'] = [700] + trace['stim_end'] = [2700] + + features = ["burst_end_indices", "peak_time", "interburst_duration"] + + feature_values = efel.getFeatureValues( + [trace], + features, + raise_warnings=False + ) + + burst_endi = feature_values[0]["burst_end_indices"] + peak_time = feature_values[0]["peak_time"] + inter_dur = feature_values[0]["interburst_duration"] + + inter_dur_py = None + if burst_endi is not None: + print(peak_time) + print(burst_endi) + inter_dur_py = [ + peak_time[idx + 1] - peak_time[idx] + for idx in burst_endi + if idx + 1 < len(peak_time) + ] + + inter_dur = numpy.array(inter_dur, dtype=numpy.float64) + inter_dur_py = numpy.array(inter_dur_py, dtype=numpy.float64) + numpy.testing.assert_allclose(inter_dur, inter_dur_py) + + def test_spontaneous_firing(): """basic: Test features with spontaneous firing""" import efel diff --git a/tests/testdata/allfeatures/expectedresults.json b/tests/testdata/allfeatures/expectedresults.json index 06db917d..27463db9 100644 --- a/tests/testdata/allfeatures/expectedresults.json +++ b/tests/testdata/allfeatures/expectedresults.json @@ -1,4 +1,28 @@ { + "ADP_peak_amplitude": [ + 0.0, + 2.6873666528663236, + 0.48122323161050673, + 0.06877493979540361, + 0.0, + 2.99957475153942 + ], + "ADP_peak_indices": [ + 7115, + 9317, + 14162, + 17198, + 23943, + 27000 + ], + "ADP_peak_values": [ + -47.716407064147894, + -41.60422947192351, + -42.066708105903004, + -41.97919073015028, + -41.12299346923828, + -38.12341871769886 + ], "AHP1_depth_from_peak": [ 66.46542716365163 ], @@ -42,6 +66,12 @@ 44.685313698492095, 45.59806337795396 ], + "AHP_depth_slow": [ + 28.81047723194662, + 32.02907786238346, + 32.65404519917785, + 33.435288866699345 + ], "AHP_slow_time": [ 0.14554275318387552, 0.07352941176470588, @@ -290,6 +320,14 @@ 3.2999999999969987, 3.2999999999969987 ], + "AP_width_between_threshold": [ + 1.7000000000003865, + 2.700000000000614, + 2.8999999999973625, + 2.9999999999972715, + 3.2999999999969987, + 3.2999999999969987 + ], "APlast_amp": [ 37.848929573728334 ], @@ -345,6 +383,12 @@ "burst_ISI_indices": [ 3 ], + "burst_begin_indices": [ + 4 + ], + "burst_end_indices": [ + 5 + ], "burst_mean_freq": [ 3.9840637450221186, 7.990411506199836 @@ -355,7 +399,6 @@ "decay_time_constant_after_stim": [ 13.45188514197874 ], - "multiple_decay_time_constant_after_stim": null, "depolarized_base": [ -43.027576049717304, -41.06055206163736, @@ -379,6 +422,21 @@ -2.3296693271049533, -2.3944223175083863 ], + "interburst_15percent_indices": null, + "interburst_15percent_values": null, + "interburst_20percent_indices": null, + "interburst_20percent_values": null, + "interburst_25percent_indices": null, + "interburst_25percent_values": null, + "interburst_30percent_indices": null, + "interburst_30percent_values": null, + "interburst_40percent_indices": null, + "interburst_40percent_values": null, + "interburst_60percent_indices": null, + "interburst_60percent_values": null, + "interburst_duration": null, + "interburst_min_indices": [], + "interburst_min_values": [], "interburst_voltage": [ -38.820613862535545 ], @@ -447,6 +505,22 @@ -41.12299346923828, -41.12299346923828 ], + "min_between_peaks_indices": [ + 7115, + 9833, + 14285, + 17215, + 23960, + 28300 + ], + "min_between_peaks_values": [ + -47.716407064147894, + -45.90401077270508, + -42.68541014226824, + -42.06044280547385, + -41.279199137952354, + -80.43351009404397 + ], "min_voltage_between_spikes": [ -47.716407064147894, -45.90401077270508, @@ -457,6 +531,7 @@ "minimum_voltage": [ -75.87104945779376 ], + "multiple_decay_time_constant_after_stim": null, "number_initial_spikes": [ 1 ], @@ -486,6 +561,30 @@ 3.562320229253816, 4.475069908715679 ], + "postburst_adp_peak_indices": [ + 26438 + ], + "postburst_adp_peak_values": [ + -41.12299346923828 + ], + "postburst_fast_ahp_indices": [ + 26438 + ], + "postburst_fast_ahp_values": [ + -41.12299346923828 + ], + "postburst_min_indices": [ + 26465 + ], + "postburst_min_values": [ + -41.52919675004766 + ], + "postburst_slow_ahp_indices": [ + 26465 + ], + "postburst_slow_ahp_values": [ + -41.52919675004766 + ], "sag_amplitude": null, "sag_ratio1": null, "sag_ratio2": [ @@ -510,12 +609,6 @@ 2.6602184713542556, 2.8238770245407068 ], - "spikes_per_burst": [ - 2 - ], - "spikes_per_burst_diff": [], - "spikes_in_burst1_burst2_diff": null, - "spikes_in_burst1_burstlast_diff": null, "steady_state_hyper": [ -39.7063975216648 ], @@ -525,6 +618,13 @@ "steady_state_voltage_stimend": [ -38.28601014909481 ], + "strict_burst_mean_freq": [ + 7.990411506199836 + ], + "strict_burst_number": [ + 1 + ], + "strict_interburst_voltage": [], "time": [ 0.0, 0.1, @@ -30532,9 +30632,19 @@ "time_to_first_spike": [ 8.000000000092427 ], + "time_to_interburst_min": [], "time_to_last_spike": [ 1937.7999999986964 ], + "time_to_postburst_adp_peak": [ + 5.999999999994543 + ], + "time_to_postburst_fast_ahp": [ + 5.999999999994543 + ], + "time_to_postburst_slow_ahp": [ + 8.699999999992087 + ], "time_to_second_spike": [ 211.30000000013865 ], @@ -60567,6 +60677,12 @@ "initburst_sahp": null, "initburst_sahp_vb": null, "initburst_sahp_ssse": null, + "spikes_per_burst": [ + 2 + ], + "spikes_per_burst_diff": [], + "spikes_in_burst1_burst2_diff": null, + "spikes_in_burst1_burstlast_diff": null, "depol_block": null, "depol_block_bool": [ 1