diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 5f058375..dbe0e96f 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -10,7 +10,7 @@ jobs: runs-on: ubuntu-22.04 strategy: matrix: - python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"] + python-version: ["3.9", "3.10", "3.11", "3.12"] steps: - uses: actions/checkout@v4 diff --git a/.gitignore b/.gitignore index 016e9905..2244d766 100644 --- a/.gitignore +++ b/.gitignore @@ -22,3 +22,4 @@ bin lib fllog.txt *.DS_Store +pyvenv.cfg diff --git a/CHANGELOG.rst b/CHANGELOG.rst index f298230c..876a23e2 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -5,6 +5,15 @@ All notable changes to this project will be documented in this file. The format is based on `Keep a Changelog `_, and this project adheres to `Semantic Versioning `_. +5.6.0 - 2024-02 +---------------- + +- Reduce 3 alternative implementations to get ISIs into 1. +- "all_ISI_values" is recommended, "ISI_values" and "ISIs" are deprecated. +- The features depending on "ISI_values" are moved to Python and now they depend on "all_ISI_values". +- BUGFIX: single_burst_ratio, irregularity_index, burst_mean_freq, interburst_voltage features were ignoring the first two ISIs when the ignore_first_ISI was set. +- Added new feature: inv_ISI_values that computes and returns all of the inverse isi values. + 5.5.5 - 2024-01 ---------------- - Type annotate api.py's functions. diff --git a/README.md b/README.md index c19bfc1f..66dd2ec9 100644 --- a/README.md +++ b/README.md @@ -105,7 +105,7 @@ When you use this eFEL software for your research, we ask you to cite the follow Requirements ============ -* [Python 3.8+](https://www.python.org/downloads/) +* [Python 3.9+](https://www.python.org/downloads/) * [Pip](https://pip.pypa.io) (installed by default in newer versions of Python) * C++ compiler that can be used by pip * [Numpy](http://www.numpy.org) (will be installed automatically by pip) diff --git a/docs/source/api.rst b/docs/source/api.rst index 399d2bef..4e6b6b9c 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -15,6 +15,8 @@ Submodules api io + pyfeatures.cppfeature_access + pyfeatures.isi pyfeatures.multitrace pyfeatures.pyfeatures pyfeatures.validation diff --git a/docs/source/eFeatures.rst b/docs/source/eFeatures.rst index 0e71329d..6ed4057a 100644 --- a/docs/source/eFeatures.rst +++ b/docs/source/eFeatures.rst @@ -57,11 +57,10 @@ Time from the start of the stimulus to the maximum of the second peak inv_time_to_first_spike = 0 -`LibV1`_ : ISI_values -~~~~~~~~~~~~~~~~~~~~~ +`Python efeature`_ : ISI_values +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The interspike intervals (i.e. time intervals) between adjacent peaks. -If ignore_first_ISI is True, the 1st spike will not be taken into account, because some cells spike right after the stimulus onset and then stay silent for a while. - **Required features**: peak_time (ms) - **Units**: ms @@ -94,8 +93,8 @@ The interspike intervals, i.e., the time intervals between adjacent peaks. all_isi_values_vec = numpy.diff(peak_time) -`LibV5`_ : inv_first_ISI, inv_second_ISI, inv_third_ISI, inv_fourth_ISI, inv_fifth_ISI, inv_last_ISI -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +`Python efeature`_ : inv_first_ISI, inv_second_ISI, inv_third_ISI, inv_fourth_ISI, inv_fifth_ISI, inv_last_ISI +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1.0 over first/second/third/fourth/fith/last ISI; returns 0 when no ISI @@ -135,6 +134,17 @@ The interspike intervals, i.e., the time intervals between adjacent peaks. else: inv_last_ISI = 0 +`Python efeature`: inv_ISI_values +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Computes all inverse spike interval values. + +- **Required features**: peak_time (ms) +- **Units**: Hz +- **Pseudocode**: :: + + all_isi_values_vec = numpy.diff(peak_time) + inv_isi_values = 1000.0 / all_isi_values_vec `LibV5`_ : time_to_last_spike ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -214,7 +224,7 @@ The mean frequency of the firing rate The slope of a linear fit to a semilog plot of the ISI values. Attention: the 1st ISI is not taken into account unless ignore_first_ISI is set to 0. -See LibV1: ISI_values feature for more details. +See Python efeature: ISIs feature for more details. - **Required features**: t, V, stim_start, stim_end, ISI_values - **Units**: ms @@ -226,13 +236,13 @@ See LibV1: ISI_values feature for more details. ISI_semilog_slope = slope -`LibV5`_ : ISI_log_slope -~~~~~~~~~~~~~~~~~~~~~~~~ +`Python efeature`_ : ISI_log_slope +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The slope of a linear fit to a loglog plot of the ISI values. Attention: the 1st ISI is not taken into account unless ignore_first_ISI is set to 0. -See LibV1: ISI_values feature for more details. +See Python efeature: ISIs feature for more details. - **Required features**: t, V, stim_start, stim_end, ISI_values - **Units**: ms @@ -244,8 +254,8 @@ See LibV1: ISI_values feature for more details. ISI_log_slope = slope -`LibV5`_ : ISI_log_slope_skip -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +`Python efeature`_ : ISI_log_slope_skip +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The slope of a linear fit to a loglog plot of the ISI values, but not taking into account the first ISI values. @@ -265,27 +275,25 @@ However, if this number of ISI values to skip is higher than max_spike_skip, the ISI_log_slope = slope -`LibV1`_ : ISI_CV -~~~~~~~~~~~~~~~~~ +`Python efeature`_ : ISI_CV +~~~~~~~~~~~~~~~~~~~~~~~~~~~ The coefficient of variation of the ISIs. Attention: the 1st ISI is not taken into account unless ignore_first_ISI is set to 0. -See LibV1: ISI_values feature for more details. +See Python efeature: ISIs feature for more details. -- **Required features**: ISI_values +- **Required features**: ISIs - **Units**: constant - **Pseudocode**: :: ISI_mean = numpy.mean(ISI_values) - ISI_variance = numpy.sum(numpy.square(ISI_values-ISI_mean)) / (len(ISI_values)-1) - ISI_std = math.sqrt(ISI_variance) - ISI_CV = ISI_std / ISI_mean + ISI_CV = np.std(isi_values, ddof=1) / ISI_mean -`LibV5`_ : irregularity_index -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +`Python efeature`_ : irregularity_index +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -Mean of the absolute difference of all ISIs, except the first one (see LibV1: ISI_values feature for more details.) +Mean of the absolute difference of all ISIs, except the first one (see Python efeature: ISIs feature for more details.) The first ISI can be taken into account if ignore_first_ISI is set to 0. @@ -349,8 +357,8 @@ The adaptation index is zero for a constant firing rate and bigger than zero for ISI_sub = ISI_values[1:] - ISI_values[:-1] adaptation_index = numpy.mean(ISI_sum / ISI_sub) -`LibV1`_ : burst_mean_freq -~~~~~~~~~~~~~~~~~~~~~~~~~~ +`Python efeature`_ : burst_mean_freq +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The mean frequency during a burst for each burst @@ -436,8 +444,8 @@ The burst detection can be fine-tuned by changing the setting strict_burst_facto burst_number = len(strict_burst_mean_freq) -`LibV1`_ : interburst_voltage -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +`Python efeature`_ : interburst_voltage +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The voltage average in between two bursts @@ -565,8 +573,8 @@ 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) ] -`LibV1`_ : single_burst_ratio -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +`Python efeature`_ : single_burst_ratio +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Length of the second isi over the median of the rest of the isis. The first isi is not taken into account, because it could bias the feature. diff --git a/efel/DependencyV5.txt b/efel/DependencyV5.txt index 4e8320b7..12a4c4cb 100644 --- a/efel/DependencyV5.txt +++ b/efel/DependencyV5.txt @@ -1,9 +1,7 @@ LibV1:interpolate LibV5:peak_indices #LibV1:interpolate -LibV1:ISI_values #LibV1:peak_time #LibV1:interpolate LibV1:doublet_ISI #LibV1:peak_time #LibV1:interpolate LibV1:peak_voltage #LibV5:peak_indices #LibV1:interpolate -LibV1:burst_ISI_indices #LibV5:peak_indices #LibV1:ISI_values #LibV1:interpolate LibV1:mean_frequency #LibV1:peak_time #LibV1:interpolate LibV1:peak_time #LibV5:peak_indices #LibV1:interpolate LibV1:time_to_first_spike #LibV1:peak_time #LibV1:interpolate @@ -11,8 +9,6 @@ LibV1:adaptation_index #LibV1:peak_time #LibV1:interpolate LibV1:adaptation_index2 #LibV1:peak_time #LibV1:interpolate LibV1:spike_width2 #LibV5:min_AHP_indices #LibV1:interpolate LibV1:AP_width #LibV5:peak_indices #LibV5:min_AHP_indices #LibV1:interpolate -LibV1:burst_mean_freq #LibV1:burst_ISI_indices #LibV1:peak_time #LibV1:interpolate -LibV1:interburst_voltage #LibV1:burst_ISI_indices #LibV1:interpolate LibV1:AP_height #LibV1:peak_voltage #LibV1:interpolate LibV1:AP_amplitude #LibV5:AP_begin_indices #LibV1:peak_voltage #LibV1:peak_time #LibV1:interpolate LibV1:AHP_depth_abs_slow #LibV5:peak_indices #LibV1:interpolate @@ -26,7 +22,6 @@ LibV1:maximum_voltage #LibV1:interpolate LibV1:minimum_voltage #LibV1:interpolate LibV1:steady_state_voltage #LibV1:interpolate LibV3:depolarized_base #LibV5:AP_end_indices #LibV5:AP_begin_indices #LibV1:interpolate -LibV1:ISI_CV #LibV1:ISI_values #LibV1:interpolate LibV1:AHP_depth #LibV5:voltage_base #LibV5:min_AHP_values #LibV1:interpolate LibV1:AHP_depth_slow #LibV5:voltage_base #LibV1:AHP_depth_abs_slow #LibV1:interpolate LibV2:AP_rise_indices #LibV5:peak_indices #LibV5:AP_begin_indices #LibV1:interpolate @@ -45,21 +40,16 @@ LibV2:AP_rise_rate_change #LibV2:AP_rise_rate #LibV1:interpolate LibV2:AP_fall_rate_change #LibV2:AP_fall_rate #LibV1:interpolate LibV2:fast_AHP_change #LibV2:fast_AHP #LibV1:interpolate LibV2:AP_duration_half_width_change #LibV2:AP_duration_half_width #LibV1:interpolate -LibV1:single_burst_ratio #LibV1:ISI_values #LibV1:interpolate LibV2:steady_state_hyper #LibV1:interpolate LibV2:amp_drop_first_second #LibV1:peak_voltage #LibV1:interpolate LibV2:amp_drop_first_last #LibV1:peak_voltage #LibV1:interpolate LibV2:amp_drop_second_last #LibV1:peak_voltage #LibV1:interpolate LibV2:max_amp_difference #LibV1:peak_voltage #LibV1:interpolate LibV1:AP_amplitude_diff #LibV1:AP_amplitude #LibV1:interpolate -LibV5:ISI_log_slope #LibV1:ISI_values #LibV1:interpolate -LibV5:ISI_semilog_slope #LibV1:ISI_values #LibV1:interpolate -LibV5:ISI_log_slope_skip #LibV1:ISI_values #LibV1:interpolate LibV1:AHP_depth_diff #LibV1:AHP_depth #LibV1:interpolate LibV5:min_AHP_indices #LibV5:peak_indices #LibV1:interpolate LibV5:min_AHP_values #LibV5:min_AHP_indices #LibV1:interpolate LibV5:number_initial_spikes #LibV1:peak_time #LibV1:interpolate -LibV5:irregularity_index #LibV1:ISI_values #LibV1:interpolate LibV5:AP1_amp #LibV1:AP_amplitude #LibV1:interpolate LibV5:APlast_amp #LibV1:AP_amplitude #LibV1:interpolate LibV5:AP2_amp #LibV1:AP_amplitude #LibV1:interpolate @@ -76,12 +66,6 @@ LibV5:AHP1_depth_from_peak #LibV5:AHP_depth_from_peak #LibV1:interpolate LibV5:AHP2_depth_from_peak #LibV5:AHP_depth_from_peak #LibV1:interpolate LibV5:time_to_second_spike #LibV1:peak_time #LibV1:interpolate LibV5:time_to_last_spike #LibV1:peak_time #LibV1:interpolate -LibV5:inv_first_ISI #LibV5:all_ISI_values #LibV1:interpolate -LibV5:inv_second_ISI #LibV5:all_ISI_values #LibV1:interpolate -LibV5:inv_third_ISI #LibV5:all_ISI_values #LibV1:interpolate -LibV5:inv_fourth_ISI #LibV5:all_ISI_values #LibV1:interpolate -LibV5:inv_fifth_ISI #LibV5:all_ISI_values #LibV1:interpolate -LibV5:inv_last_ISI #LibV5:all_ISI_values #LibV1:interpolate LibV5:inv_time_to_first_spike #LibV1:time_to_first_spike #LibV1:interpolate LibV5:spike_half_width #LibV5:min_AHP_indices #LibV5:peak_indices #LibV1:interpolate LibV5:AP_begin_indices #LibV5:min_AHP_indices #LibV5:peak_indices #LibV1:interpolate diff --git a/efel/cppcore/FillFptrTable.cpp b/efel/cppcore/FillFptrTable.cpp index 5cb1a6a0..6c60463a 100644 --- a/efel/cppcore/FillFptrTable.cpp +++ b/efel/cppcore/FillFptrTable.cpp @@ -20,16 +20,12 @@ int FillFptrTable() { //****************** for FptrTableV1 ***************************** FptrTableV1["interpolate"] = &LibV1::interpolate; - FptrTableV1["ISI_values"] = &LibV1::ISI_values; FptrTableV1["peak_voltage"] = &LibV1::peak_voltage; FptrTableV1["mean_frequency"] = &LibV1::firing_rate; FptrTableV1["peak_time"] = &LibV1::peak_time; FptrTableV1["time_to_first_spike"] = &LibV1::first_spike_time; - FptrTableV1["burst_ISI_indices"] = &LibV1::burst_ISI_indices; FptrTableV1["adaptation_index"] = &LibV1::adaptation_index; FptrTableV1["spike_width2"] = &LibV1::spike_width2; - FptrTableV1["burst_mean_freq"] = &LibV1::burst_mean_freq; - FptrTableV1["interburst_voltage"] = &LibV1::interburst_voltage; // passive properties FptrTableV1["time_constant"] = &LibV1::time_constant; FptrTableV1["voltage_deflection"] = &LibV1::voltage_deflection; @@ -40,11 +36,9 @@ int FillFptrTable() { FptrTableV1["AP_height"] = &LibV1::AP_height; FptrTableV1["AP_amplitude"] = &LibV1::AP_amplitude; - FptrTableV1["single_burst_ratio"] = &LibV1::single_burst_ratio; FptrTableV1["AP_width"] = &LibV1::AP_width; FptrTableV1["doublet_ISI"] = &LibV1::doublet_ISI; FptrTableV1["adaptation_index2"] = &LibV1::adaptation_index2; - FptrTableV1["ISI_CV"] = &LibV1::ISI_CV; FptrTableV1["AHP_depth_abs_slow"] = &LibV1::AHP_depth_abs_slow; FptrTableV1["AHP_slow_time"] = &LibV1::AHP_slow_time; FptrTableV1["AHP_depth"] = &LibV1::AHP_depth; @@ -94,38 +88,8 @@ int FillFptrTable() { //****************** end of FptrTableV3 ***************************** //****************** FptrTableV5 ***************************** - - FptrTableV5["ISI_log_slope"] = &LibV5::ISI_log_slope; - FptrTableV5["ISI_semilog_slope"] = &LibV5::ISI_semilog_slope; - FptrTableV5["ISI_log_slope_skip"] = &LibV5::ISI_log_slope_skip; FptrTableV5["time_to_second_spike"] = &LibV5::time_to_second_spike; FptrTableV5["time_to_last_spike"] = &LibV5::time_to_last_spike; - FptrTableV5["inv_first_ISI"] = [](mapStr2intVec& intData, - mapStr2doubleVec& doubleData, - mapStr2Str& strData) { - return LibV5::inv_ISI_generic(intData, doubleData, strData, 0); - }; - FptrTableV5["inv_second_ISI"] = [](mapStr2intVec& intData, - mapStr2doubleVec& doubleData, - mapStr2Str& strData) { - return LibV5::inv_ISI_generic(intData, doubleData, strData, 1); - }; - FptrTableV5["inv_third_ISI"] = [](mapStr2intVec& intData, - mapStr2doubleVec& doubleData, - mapStr2Str& strData) { - return LibV5::inv_ISI_generic(intData, doubleData, strData, 2); - }; - FptrTableV5["inv_fourth_ISI"] = [](mapStr2intVec& intData, - mapStr2doubleVec& doubleData, - mapStr2Str& strData) { - return LibV5::inv_ISI_generic(intData, doubleData, strData, 3); - }; - FptrTableV5["inv_fifth_ISI"] = [](mapStr2intVec& intData, - mapStr2doubleVec& doubleData, - mapStr2Str& strData) { - return LibV5::inv_ISI_generic(intData, doubleData, strData, 4); - }; - FptrTableV5["inv_last_ISI"] = &LibV5::inv_last_ISI; FptrTableV5["inv_time_to_first_spike"] = &LibV5::inv_time_to_first_spike; FptrTableV5["min_AHP_indices"] = &LibV5::min_AHP_indices; FptrTableV5["min_AHP_values"] = &LibV5::min_AHP_values; @@ -133,7 +97,6 @@ int FillFptrTable() { FptrTableV5["spike_half_width"] = &LibV5::spike_width1; FptrTableV5["AP_begin_indices"] = &LibV5::AP_begin_indices; FptrTableV5["AP_end_indices"] = &LibV5::AP_end_indices; - FptrTableV5["irregularity_index"] = &LibV5::irregularity_index; FptrTableV5["number_initial_spikes"] = &LibV5::number_initial_spikes; FptrTableV5["AP1_amp"] = &LibV5::AP1_amp; FptrTableV5["APlast_amp"] = &LibV5::APlast_amp; diff --git a/efel/cppcore/LibV1.cpp b/efel/cppcore/LibV1.cpp index b73f9cf2..579b1fa4 100644 --- a/efel/cppcore/LibV1.cpp +++ b/efel/cppcore/LibV1.cpp @@ -76,61 +76,6 @@ int LibV1::interpolate(mapStr2intVec& IntFeatureData, return 1; } -int LibV1::ISI_values(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - const auto& doubleFeatures = getFeatures(DoubleFeatureData, {"peak_time"}); - if (doubleFeatures.at("peak_time").size() < 3) - throw FeatureComputationError("Three spikes required for calculation of ISI_values."); - - const auto& intFeatures = getFeatures(IntFeatureData, {"ignore_first_ISI"}); - int IgnoreFirstISI = 1; - if (intFeatures.at("ignore_first_ISI").size() > 0 && - intFeatures.at("ignore_first_ISI")[0] == 0) { - IgnoreFirstISI = 0; - } - - vector VecISI; - for (size_t i = IgnoreFirstISI + 1; i < doubleFeatures.at("peak_time").size(); - i++) - VecISI.push_back(doubleFeatures.at("peak_time")[i] - - doubleFeatures.at("peak_time")[i - 1]); - setVec(DoubleFeatureData, StringData, "ISI_values", VecISI); - return VecISI.size(); -} - -// *** ISI_CV *** -// the coefficient of variation of the ISI -static int __ISI_CV(const vector& isivalues, vector& isicv) { - // mean - double isi_mean = 0.; - for (size_t i = 0; i < isivalues.size(); i++) { - isi_mean += isivalues[i]; - } - isi_mean /= isivalues.size(); - - // sigma^2 - double variance = 0.; - for (size_t i = 0; i < isivalues.size(); i++) { - double dev = isivalues[i] - isi_mean; - variance += dev * dev; - } - // variation coefficient cv = sigma / mean - isicv.push_back(sqrt(variance / (isivalues.size() - 1)) / isi_mean); - return isicv.size(); -} -int LibV1::ISI_CV(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData) { - const auto& doubleFeatures = getFeatures(DoubleFeatureData, {"ISI_values"}); - if (doubleFeatures.at("ISI_values").size() < 2) return -1; - - vector isicv; - int retval = __ISI_CV(doubleFeatures.at("ISI_values"), isicv); - if (retval >= 0) { - setVec(DoubleFeatureData, StringData, "ISI_CV", isicv); - } - return retval; -} int LibV1::peak_voltage(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, @@ -302,168 +247,6 @@ int LibV1::AHP_slow_time(mapStr2intVec& IntFeatureData, return -1; } -static int __burst_ISI_indices(double BurstFactor, const vector& PeakIndex, - const vector& ISIValues, - vector& BurstIndex) { - vector ISIpcopy; - int n, count = -1; - double dMedian; - - for (size_t i = 1; i < (ISIValues.size() - 1); i++) { - ISIpcopy.clear(); - for (size_t j = count + 1; j < i; j++) ISIpcopy.push_back(ISIValues[j]); - sort(ISIpcopy.begin(), ISIpcopy.end()); - n = ISIpcopy.size(); - if ((n % 2) == 0) { - dMedian = - (ISIpcopy[int((n - 1) / 2)] + ISIpcopy[int((n - 1) / 2) + 1]) / 2; - } else { - dMedian = ISIpcopy[int(n / 2)]; - } - if (ISIValues[i] > (BurstFactor * dMedian) && - (ISIValues[i + 1] < ISIValues[i] / BurstFactor)) { - BurstIndex.push_back(i + 1); - count = i - 1; - } - } - return BurstIndex.size(); -} - -int LibV1::burst_ISI_indices(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - const auto& doubleFeatures = - getFeatures(DoubleFeatureData, {"ISI_values", "burst_factor"}); - const auto& intFeatures = getFeatures(IntFeatureData, {"peak_indices"}); - if (intFeatures.at("peak_indices").size() < 5) - throw FeatureComputationError("More than 5 spikes are needed for burst calculation."); - double BurstFactor = (doubleFeatures.at("burst_factor").empty()) - ? 2 - : doubleFeatures.at("burst_factor")[0]; - vector BurstIndex; - int retVal = __burst_ISI_indices(BurstFactor, intFeatures.at("peak_indices"), - doubleFeatures.at("ISI_values"), BurstIndex); - if (retVal >= 0) { - setVec(IntFeatureData, StringData, "burst_ISI_indices", BurstIndex); - } - return retVal; -} - -// discrepancy betwen PVTime indices and BurstIndex indices because -// first ISI value is ignored -static int __burst_mean_freq(const vector& PVTime, - const vector& BurstIndex, int IgnoreFirstISI, - vector& BurstMeanFreq) { - // if no burst detected, do not consider all peaks in a single burst - if (BurstIndex.size() == 0) return BurstMeanFreq.size(); - double span; - size_t i; - - // 1st burst - span = PVTime[BurstIndex[0] - 1 + IgnoreFirstISI] - PVTime[0]; - BurstMeanFreq.push_back((BurstIndex[0] + IgnoreFirstISI) * 1000 / span); - - for (i = 0; i < BurstIndex.size() - 1; i++) { - if (BurstIndex[i + 1] - BurstIndex[i] > 1) { - span = PVTime[BurstIndex[i + 1] - 1 + IgnoreFirstISI] - - PVTime[BurstIndex[i] + IgnoreFirstISI]; - BurstMeanFreq.push_back((BurstIndex[i + 1] - BurstIndex[i]) * 1000 / - span); - } - } - - // last burst - if (PVTime.size() - IgnoreFirstISI - BurstIndex[i] > 1) { - span = PVTime[PVTime.size() - 1] - PVTime[BurstIndex[i] + IgnoreFirstISI]; - BurstMeanFreq.push_back((PVTime.size() - IgnoreFirstISI - BurstIndex[i]) * - 1000 / span); - } - - return BurstMeanFreq.size(); -} - -int LibV1::burst_mean_freq(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - const auto& doubleFeatures = getFeatures(DoubleFeatureData, {"peak_time"}); - const auto& intFeatures = - getFeatures(IntFeatureData, {"burst_ISI_indices"}); - int IgnoreFirstISI = 1; - int retValIgnore; - vector retIgnore; - retValIgnore = getParam(IntFeatureData, "ignore_first_ISI", retIgnore); - if ((retValIgnore == 1) && (retIgnore.size() > 0) && (retIgnore[0] == 0)) { - IgnoreFirstISI = 0; - } - vector BurstMeanFreq; - int retVal = __burst_mean_freq(doubleFeatures.at("peak_time"), - intFeatures.at("burst_ISI_indices"), - IgnoreFirstISI, BurstMeanFreq); - if (retVal >= 0) { - setVec(DoubleFeatureData, StringData, "burst_mean_freq", BurstMeanFreq); - } - return retVal; -} - -// reminder: first ISI value is ignored in burst_ISI_indices if IgnoreFirstISI=1 -static int __interburst_voltage(const vector& BurstIndex, - const vector& PeakIndex, - const vector& T, - const vector& V, int IgnoreFirstISI, - vector& IBV) { - if (BurstIndex.size() < 1) return 0; - int j, pIndex, tsIndex, teIndex, cnt; - double tStart, tEnd, vTotal = 0; - for (size_t i = 0; i < BurstIndex.size(); i++) { - pIndex = BurstIndex[i] + IgnoreFirstISI - 1; - tsIndex = PeakIndex[pIndex]; - tStart = T[tsIndex] + 5; // 5Millisecond after - pIndex = BurstIndex[i] + IgnoreFirstISI; - teIndex = PeakIndex[pIndex]; - tEnd = T[teIndex] - 5; // 5Millisecond before - - for (j = tsIndex; j < teIndex; j++) { - if (T[j] > tStart) break; - } - tsIndex = --j; - - for (j = teIndex; j > tsIndex; j--) { - if (T[j] < tEnd) break; - } - teIndex = ++j; - vTotal = 0; - for (j = tsIndex, cnt = 1; j <= teIndex; j++, cnt++) vTotal = vTotal + V[j]; - if (cnt == 0) continue; - IBV.push_back(vTotal / (cnt - 1)); - } - return IBV.size(); -} - -int LibV1::interburst_voltage(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - const auto& doubleFeatures = getFeatures(DoubleFeatureData, {"T", "V"}); - const auto& intFeatures = - getFeatures(IntFeatureData, {"peak_indices", "burst_ISI_indices"}); - int IgnoreFirstISI; - int retValIgnore; - vector retIgnore; - retValIgnore = getParam(IntFeatureData, "ignore_first_ISI", retIgnore); - if ((retValIgnore == 1) && (retIgnore.size() > 0) && (retIgnore[0] == 0)) { - IgnoreFirstISI = 0; - } else { - IgnoreFirstISI = 1; - } - vector IBV; - int retVal = __interburst_voltage( - intFeatures.at("burst_ISI_indices"), intFeatures.at("peak_indices"), - doubleFeatures.at("T"), doubleFeatures.at("V"), IgnoreFirstISI, IBV); - if (retVal >= 0) { - setVec(DoubleFeatureData, StringData, "interburst_voltage", IBV); - } - return retVal; -} - static int __adaptation_index(double spikeSkipf, int maxnSpike, double StimStart, double StimEnd, double Offset, const vector& peakVTime, @@ -1057,37 +840,6 @@ int LibV1::steady_state_voltage(mapStr2intVec& IntFeatureData, return retVal; } -// *** single_burst_ratio *** -// according to Shaul: measures the length of the first isi over the median of -// the rest of the isis -static int __single_burst_ratio(const vector& isivalues, - vector& singleburstratio) { - if (isivalues.size() < 2) { - return 0; - } - // calculate the average instead of the median - double average = 0.; - for (size_t i = 1; i < isivalues.size(); i++) { - average += isivalues[i]; - } - average /= isivalues.size() - 1; - singleburstratio.push_back(isivalues[0] / average); - return singleburstratio.size(); -} - -int LibV1::single_burst_ratio(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - const auto& doubleFeatures = getFeatures(DoubleFeatureData, {"ISI_values"}); - vector singleburstratio; - int retval = - __single_burst_ratio(doubleFeatures.at("ISI_values"), singleburstratio); - if (retval > 0) { - setVec(DoubleFeatureData, StringData, "single_burst_ratio", - singleburstratio); - } - return retval; -} // *** AP_width *** // diff --git a/efel/cppcore/LibV1.h b/efel/cppcore/LibV1.h index ed2781ae..0431d1e4 100644 --- a/efel/cppcore/LibV1.h +++ b/efel/cppcore/LibV1.h @@ -30,9 +30,6 @@ namespace LibV1 { int interpolate(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); -int ISI_values(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); - int peak_voltage(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); @@ -49,18 +46,6 @@ int first_spike_time(mapStr2intVec& IntFeatureData, int spike_width2(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); -int burst_ISI_indices(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData); - -int burst_mean_freq(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData); - -int interburst_voltage(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData); - int adaptation_index(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); @@ -94,11 +79,6 @@ int AP_amplitude(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); int AHP_depth(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); - -int single_burst_ratio(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData); - int AP_width(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); int doublet_ISI(mapStr2intVec& IntFeatureData, @@ -106,8 +86,6 @@ int doublet_ISI(mapStr2intVec& IntFeatureData, int adaptation_index2(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); -int ISI_CV(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData); int AHP_depth_abs_slow(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); diff --git a/efel/cppcore/LibV5.cpp b/efel/cppcore/LibV5.cpp index 15c457fd..28a7f01c 100644 --- a/efel/cppcore/LibV5.cpp +++ b/efel/cppcore/LibV5.cpp @@ -32,105 +32,6 @@ using std::distance; using std::find_if; -// slope of loglog of ISI curve -static int __ISI_log_slope(const vector& isiValues, - vector& slope, bool skip, double spikeSkipf, - size_t maxnSpike, bool semilog) { - std::deque skippedISIValues; - - vector log_isivalues; - vector x; - - for (size_t i = 0; i < isiValues.size(); i++) { - skippedISIValues.push_back(isiValues[i]); - } - - if (skip) { - // Remove n spikes given by spike_skipf or max_spike_skip - size_t isisToRemove = (size_t)((isiValues.size() + 1) * spikeSkipf + .5); - - isisToRemove = std::min(maxnSpike, isisToRemove); - - // Remove spikeToRemove spike from SpikeTime list - for (size_t i = 0; i < isisToRemove; i++) { - skippedISIValues.pop_front(); - } - } - - for (size_t i = 0; i < skippedISIValues.size(); i++) { - log_isivalues.push_back(log(skippedISIValues[i])); - if (semilog) { - x.push_back((double)i + 1); - } else { - x.push_back(log((double)i + 1)); - } - } - - if (x.size() == 0 || log_isivalues.size() == 0) return -1; - - linear_fit_result fit; - fit = slope_straight_line_fit(x, log_isivalues); - - if (fit.slope == 0. || is_nan(fit.slope)) return -1; - - slope.push_back(fit.slope); - - return slope.size(); -} - -int LibV5::ISI_log_slope(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - const auto& isivalues = getFeature(DoubleFeatureData, {"ISI_values"}); - vector slope; - bool semilog = false; - int retval = __ISI_log_slope(isivalues, slope, false, 0.0, 0, semilog); - if (retval >= 0) { - setVec(DoubleFeatureData, StringData, "ISI_log_slope", slope); - return static_cast(slope.size()); - } else { - return retval; - } -} - -int LibV5::ISI_semilog_slope(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - const auto& isivalues = getFeature(DoubleFeatureData, {"ISI_values"}); - vector slope; - bool semilog = true; - int retval = __ISI_log_slope(isivalues, slope, false, 0.0, 0, semilog); - if (retval >= 0) { - setVec(DoubleFeatureData, StringData, "ISI_semilog_slope", slope); - return static_cast(slope.size()); - } else { - return retval; - } -} - -int LibV5::ISI_log_slope_skip(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - const auto& isivalues = getFeature(DoubleFeatureData, {"ISI_values"}); - const auto spikeSkipf = - getFeature(DoubleFeatureData, {"spike_skipf"}).front(); - const auto maxnSpike = getFeature(IntFeatureData, {"max_spike_skip"}).front(); - - // Check the validity of spikeSkipf value - if (spikeSkipf < 0 || spikeSkipf >= 1) - throw FeatureComputationError("spike_skipf should lie between [0 1)."); - - vector slope; - bool semilog = false; - int retVal = - __ISI_log_slope(isivalues, slope, true, spikeSkipf, maxnSpike, semilog); - if (retVal >= 0) { - setVec(DoubleFeatureData, StringData, "ISI_log_slope_skip", slope); - return static_cast(slope.size()); - } else { - return retVal; - } -} // time from stimulus start to second threshold crossing int LibV5::time_to_second_spike(mapStr2intVec& IntFeatureData, @@ -163,66 +64,6 @@ int LibV5::time_to_last_spike(mapStr2intVec& IntFeatureData, return 1; } -double calculateInvISI(const std::vector& all_isi_values_vec, - size_t index) { - if (index < all_isi_values_vec.size()) { - return 1000.0 / all_isi_values_vec[index]; - } - throw FeatureComputationError("inverse ISI index out of range"); -} - -int LibV5::inv_ISI_generic(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData, size_t index) { - const auto& all_isi_values_vec = - getFeature(DoubleFeatureData, {"all_ISI_values"}); - double inv_ISI = calculateInvISI(all_isi_values_vec, index); - std::string featureName; - - switch (index) { - case 0: - featureName = "inv_first_ISI"; - break; - case 1: - featureName = "inv_second_ISI"; - break; - case 2: - featureName = "inv_third_ISI"; - break; - case 3: - featureName = "inv_fourth_ISI"; - break; - case 4: - featureName = "inv_fifth_ISI"; - break; - default: - if (index == all_isi_values_vec.size() - 1) { - featureName = "inv_last_ISI"; - } else { - featureName = "inv_" + std::to_string(index + 1) + "th_ISI"; - } - break; - } - - vector inv_ISI_vec = {inv_ISI}; - setVec(DoubleFeatureData, StringData, featureName, inv_ISI_vec); - return 1; -} - -// 1.0 over last ISI (in Hz); returns 0 when no ISI -int LibV5::inv_last_ISI(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - const auto& all_isi_values_vec = - getFeature(DoubleFeatureData, "all_ISI_values"); - - double inv_last_ISI = calculateInvISI( - all_isi_values_vec, all_isi_values_vec.size() - 1); // Last ISI - vector inv_last_ISI_vec = {inv_last_ISI}; - setVec(DoubleFeatureData, StringData, "inv_last_ISI", inv_last_ISI_vec); - return 1; -} - // 1.0 over time to first spike (in Hz); returns 0 when no spike int LibV5::inv_time_to_first_spike(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, @@ -595,35 +436,6 @@ int LibV5::AP_end_indices(mapStr2intVec& IntFeatureData, return retVal; } -static int __irregularity_index(const vector& isiValues, - vector& irregularity_index) { - double ISISub, iRI; - iRI = ISISub = 0; - if (isiValues.size() == 0) return -1; - - for (size_t i = 1; i < isiValues.size(); i++) { - ISISub = std::abs(isiValues[i] - isiValues[i - 1]); - iRI = iRI + (ISISub); - } - iRI = iRI / isiValues.size(); - irregularity_index.clear(); - irregularity_index.push_back(iRI); - return 1; -} - -int LibV5::irregularity_index(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - const vector& isiValues = getFeature(DoubleFeatureData, "ISI_values"); - vector irregularity_index; - int retVal = __irregularity_index(isiValues, irregularity_index); - if (retVal >= 0) { - setVec(DoubleFeatureData, StringData, "irregularity_index", - irregularity_index); - } - return retVal; -} - static int __number_initial_spikes(const vector& peak_times, double stimstart, double stimend, double initial_perc, diff --git a/efel/cppcore/LibV5.h b/efel/cppcore/LibV5.h index 3556a9c3..95c4d019 100644 --- a/efel/cppcore/LibV5.h +++ b/efel/cppcore/LibV5.h @@ -28,14 +28,6 @@ using std::vector; namespace LibV5 { -int ISI_log_slope(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); -int ISI_semilog_slope(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); -int ISI_log_slope_skip(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData); - int time_to_second_spike(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); @@ -72,10 +64,6 @@ int AP_end_indices(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); -int irregularity_index(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData); - int number_initial_spikes(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); diff --git a/efel/cppcore/Utils.cpp b/efel/cppcore/Utils.cpp index 397497e0..03bb8d58 100644 --- a/efel/cppcore/Utils.cpp +++ b/efel/cppcore/Utils.cpp @@ -256,5 +256,13 @@ double vec_median(vector v) { return (double)(v[(n - 1) / 2] + v[n / 2]) / 2.0; // even } +void removeFirstISI(vector& vec) { + if (vec.empty()) { + throw std::runtime_error("Cannot remove from an empty vector."); + } + + vec.erase(vec.begin()); // Remove the first element from the vector. +} + template double vec_mean(const vector& v); template double vec_median(vector v); diff --git a/efel/cppcore/Utils.h b/efel/cppcore/Utils.h index 463d36db..a49ed8f3 100644 --- a/efel/cppcore/Utils.h +++ b/efel/cppcore/Utils.h @@ -57,6 +57,8 @@ double vec_median(vector v); template double vec_mean(const vector &v); +void removeFirstISI(vector& vec); + std::pair get_time_index(const std::vector &t, double startTime, double endTime, double precisionThreshold); diff --git a/efel/cppcore/cfeature.cpp b/efel/cppcore/cfeature.cpp index 57a56f72..f2b2c956 100644 --- a/efel/cppcore/cfeature.cpp +++ b/efel/cppcore/cfeature.cpp @@ -75,9 +75,7 @@ const vector cFeature::getmapDoubleData(string strName) { void cFeature::fillfeaturetypes() { // initialize feature type information featuretypes["peak_indices"] = "int"; - featuretypes["ISI_values"] = "double"; featuretypes["peak_voltage"] = "double"; - featuretypes["burst_ISI_indices"] = "int"; featuretypes["mean_frequency"] = "double"; featuretypes["peak_time"] = "double"; featuretypes["time_to_first_spike"] = "double"; @@ -86,8 +84,6 @@ void cFeature::fillfeaturetypes() { featuretypes["adaptation_index"] = "double"; featuretypes["spike_width2"] = "double"; featuretypes["spike_half_width"] = "double"; - featuretypes["burst_mean_freq"] = "double"; - featuretypes["interburst_voltage"] = "double"; featuretypes["voltage_base"] = "double"; featuretypes["AHP_depth_abs"] = "double"; featuretypes["time_constant"] = "double"; @@ -119,12 +115,10 @@ void cFeature::fillfeaturetypes() { featuretypes["AP_duration_half_width_change"] = "double"; featuretypes["AP_height"] = "double"; featuretypes["AP_amplitude"] = "double"; - featuretypes["single_burst_ratio"] = "double"; featuretypes["steady_state_hyper"] = "double"; featuretypes["AP_width"] = "double"; featuretypes["doublet_ISI"] = "double"; featuretypes["adaptation_index2"] = "double"; - featuretypes["ISI_CV"] = "double"; featuretypes["AHP_depth_abs_slow"] = "double"; featuretypes["AHP_slow_time"] = "double"; featuretypes["AHP_depth"] = "double"; @@ -136,19 +130,9 @@ void cFeature::fillfeaturetypes() { featuretypes["depolarized_base"] = "double"; // LibV5 - featuretypes["ISI_log_slope"] = "double"; - featuretypes["ISI_semilog_slope"] = "double"; - featuretypes["ISI_log_slope_skip"] = "double"; featuretypes["time_to_second_spike"] = "double"; featuretypes["time_to_last_spike"] = "double"; - featuretypes["inv_first_ISI"] = "double"; - featuretypes["inv_second_ISI"] = "double"; - featuretypes["inv_third_ISI"] = "double"; - featuretypes["inv_fourth_ISI"] = "double"; - featuretypes["inv_fifth_ISI"] = "double"; - featuretypes["inv_last_ISI"] = "double"; featuretypes["inv_time_to_first_spike"] = "double"; - featuretypes["irregularity_index"] = "double"; featuretypes["number_initial_spikes"] = "int"; featuretypes["AP1_amp"] = "double"; featuretypes["APlast_amp"] = "double"; diff --git a/efel/pyfeatures/cppfeature_access.py b/efel/pyfeatures/cppfeature_access.py new file mode 100644 index 00000000..c7f97509 --- /dev/null +++ b/efel/pyfeatures/cppfeature_access.py @@ -0,0 +1,26 @@ +"""Module containing access functions to C++ features for Python features.""" +from __future__ import annotations +import warnings +import numpy as np +from efel import cppcore + + +def get_cpp_feature(feature_name: str, raise_warnings=False) -> np.ndarray | None: + """Return value of feature implemented in cpp.""" + cppcoreFeatureValues: list[int | float] = list() + exitCode = cppcore.getFeature(feature_name, cppcoreFeatureValues) + if exitCode < 0: + if raise_warnings: + warnings.warn( + f"Error while calculating {feature_name}, {cppcore.getgError()}", + RuntimeWarning) + return None + return np.array(cppcoreFeatureValues) + + +def _get_cpp_data(data_name: str) -> float | int: + """Get cpp data value.""" + try: + return cppcore.getMapDoubleData(data_name)[0] + except Exception: + return cppcore.getMapIntData(data_name)[0] diff --git a/efel/pyfeatures/isi.py b/efel/pyfeatures/isi.py new file mode 100644 index 00000000..f90f0b42 --- /dev/null +++ b/efel/pyfeatures/isi.py @@ -0,0 +1,434 @@ +"""Features that are depending on the inter-spike intervals.""" +from __future__ import annotations +import logging +from typing_extensions import deprecated +import warnings +import numpy as np +from efel.pyfeatures.cppfeature_access import ( + _get_cpp_data, get_cpp_feature +) + + +logger = logging.getLogger(__name__) + + +@deprecated("Use all_ISI_values instead") +def ISIs() -> np.ndarray | None: + """Get all ISIs, inter-spike intervals.""" + return get_cpp_feature("all_ISI_values") + + +@deprecated("Use ISIs instead.") +def ISI_values() -> np.ndarray | None: + """Get all ISIs, inter-spike intervals.""" + isi_values = get_cpp_feature("all_ISI_values") + if isi_values is None: + return None + + # Check "ignore_first_ISI" flag + ignore_first_ISI = _get_cpp_data("ignore_first_ISI") + if ignore_first_ISI: + isi_values = isi_values[1:] + return isi_values + + +def __ISI_CV(isi_values) -> float | None: + if len(isi_values) < 2: + return None + + # Calculate mean + isi_mean = np.mean(isi_values) + # Calculate coefficient of variation + cv = np.std(isi_values, ddof=1) / isi_mean # ddof 1 to replicate C++ impl + return cv + + +def ISI_CV() -> np.ndarray | None: + """Coefficient of variation of ISIs. + + If the ignore_first_ISI flag is set, the first ISI will be ignored. + """ + isi_values = get_cpp_feature("all_ISI_values") + if isi_values is None: + return None + + # Check "ignore_first_ISI" flag + ignore_first_ISI = _get_cpp_data("ignore_first_ISI") + if ignore_first_ISI: + isi_values = isi_values[1:] + + result = __ISI_CV(isi_values) + return np.array([result]) + + +def single_burst_ratio() -> np.ndarray | None: + """Calculates the single burst ratio. + + The ratio is the length of the first ISI over the average of the rest. + If the ignore_first_ISI flag is set, the first ISI will be ignored. + """ + isi_values = get_cpp_feature("all_ISI_values") + if isi_values is None: + return None + + # Check "ignore_first_ISI" flag + ignore_first_ISI = _get_cpp_data("ignore_first_ISI") + if ignore_first_ISI: + isi_values = isi_values[1:] + + if len(isi_values) < 2: + return None + + single_burst_ratio_value = isi_values[0] / np.mean(isi_values) + return np.array([single_burst_ratio_value]) + + +def irregularity_index() -> np.ndarray | None: + """Calculate the irregularity index of ISI values. + + If the ignore_first_ISI flag is set, the first ISI will be ignored. + """ + isi_values = get_cpp_feature("all_ISI_values") + if isi_values is None: + return None + + # Check "ignore_first_ISI" flag + ignore_first_ISI = _get_cpp_data("ignore_first_ISI") + if ignore_first_ISI: + isi_values = isi_values[1:] + + # Calculate the absolute differences between consecutive ISI values + isi_differences = np.abs(np.diff(isi_values)) + result = np.mean(isi_differences) + + return np.array([result]) + + +def _isi_log_slope_core( + isi_values, skip=False, spike_skipf=0.0, max_spike_skip=0, semilog=False +) -> np.ndarray | None: + if isi_values is None or len(isi_values) == 0: + return None + + if skip: + isisToRemove = min( + max_spike_skip, int((len(isi_values) + 1) * spike_skipf + 0.5) + ) + isi_values = isi_values[isisToRemove:] + + log_isi_values = np.log(isi_values) + x = np.arange(1, len(log_isi_values) + 1) + if not semilog: + x = np.log(x) + + try: + slope, _ = np.polyfit(x, log_isi_values, 1) + except np.linalg.LinAlgError as e: + warnings.warn(f"Error in polyfit: {e}") + return None + + return np.array([slope]) + + +def ISI_log_slope() -> np.ndarray | None: + """The slope of a linear fit to a loglog plot of the ISI values. + + If the ignore_first_ISI flag is set, the first ISI will be ignored. + """ + isi_values = get_cpp_feature("all_ISI_values") + if isi_values is None: + return None + + # Check "ignore_first_ISI" flag + ignore_first_ISI = _get_cpp_data("ignore_first_ISI") + if ignore_first_ISI: + isi_values = isi_values[1:] + + return _isi_log_slope_core(isi_values, False, 0.0, 0, False) + + +def ISI_semilog_slope() -> np.ndarray | None: + """The slope of a linear fit to a semilog plot of the ISI values. + + If the ignore_first_ISI flag is set, the first ISI will be ignored. + """ + isi_values = get_cpp_feature("all_ISI_values") + if isi_values is None: + return None + + # Check "ignore_first_ISI" flag + ignore_first_ISI = _get_cpp_data("ignore_first_ISI") + if ignore_first_ISI: + isi_values = isi_values[1:] + + return _isi_log_slope_core(isi_values, False, 0.0, 0, True) + + +def ISI_log_slope_skip() -> np.ndarray | None: + """The slope of a linear fit to a loglog plot of the ISI values, + but not taking into account the first ISI values. + + Uses the spike_skipf and max_spike_skip settings to determine how many + ISIs to skip. + .""" + isi_values = get_cpp_feature("all_ISI_values") + if isi_values is None: + return None + # Check "ignore_first_ISI" flag + ignore_first_ISI = _get_cpp_data("ignore_first_ISI") + if ignore_first_ISI: + isi_values = isi_values[1:] + + spike_skipf = _get_cpp_data("spike_skipf") + if spike_skipf < 0 or spike_skipf >= 1: + raise ValueError("spike_skipf should lie between [0, 1).") + max_spike_skip = _get_cpp_data("max_spike_skip") + return _isi_log_slope_core(isi_values, True, spike_skipf, max_spike_skip, False) + + +def burst_ISI_indices() -> np.ndarray | None: + """Calculate burst ISI indices based on burst factor and ISI values.""" + # Fetching necessary data + isi_values = ISI_values() + if isi_values is None: + return None + + burst_factor = _get_cpp_data("burst_factor") + + if len(isi_values) < 4: + logger.warning("4 or more spikes are needed for burst calculation.") + return None + + burst_indices = [] + count = -1 + + for i in range(1, len(isi_values) - 1): + isi_p_copy = isi_values[count + 1: i] + n = len(isi_p_copy) + + if n == 0: + continue + + # Compute median + d_median = np.median(isi_p_copy) + + # Check burst condition + if isi_values[i] > (burst_factor * d_median) and isi_values[i + 1] < ( + isi_values[i] / burst_factor + ): + burst_indices.append(i + 1) + count = i - 1 + + if burst_indices == []: + return None + return np.array(burst_indices) + + +def burst_mean_freq() -> np.ndarray | None: + """Calculate the mean frequency of bursts.""" + # Fetching required features + peak_time = get_cpp_feature("peak_time") + burst_isi_idx = burst_ISI_indices() + if burst_isi_idx is None or peak_time is None: + return None + + burst_mean_freq = [] + burst_index_tmp = burst_isi_idx + burst_index = np.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 np.array(burst_mean_freq) + + +def interburst_voltage() -> np.ndarray | None: + """The voltage average in between two bursts.""" + peak_idx = get_cpp_feature("peak_indices") + v = get_cpp_feature("voltage") + t = get_cpp_feature("time") + burst_isi_idx = burst_ISI_indices() + if peak_idx is None or v is None or t is None or burst_isi_idx is None: + return None + + interburst_voltage = [] + for idx in burst_isi_idx: + ts_idx = peak_idx[idx] + t_start = t[ts_idx] + 5 + start_idx = np.argwhere(t < t_start)[-1][0] + + te_idx = peak_idx[idx + 1] + t_end = t[te_idx] - 5 + end_idx = np.argwhere(t > t_end)[0][0] + + interburst_voltage.append(np.mean(v[start_idx:end_idx + 1])) + + return np.array(interburst_voltage) + + +def initburst_sahp() -> np.ndarray | None: + """SlowAHP voltage after initial burst.""" + # Required cpp features + voltage = get_cpp_feature("voltage") + time = get_cpp_feature("time") + if voltage is None or time is None: + return None + time = time[:len(voltage)] + peak_times = get_cpp_feature("peak_time") + if peak_times is None: + return None + + # Required python features + all_isis = get_cpp_feature("all_ISI_values") + + # Required trace data + stim_end = _get_cpp_data("stim_end") + + # Required settings + initburst_freq_thresh = _get_cpp_data("initburst_freq_threshold") + initburst_sahp_start = _get_cpp_data("initburst_sahp_start") + initburst_sahp_end = _get_cpp_data("initburst_sahp_end") + + last_isi = None + + # Loop over ISIs until frequency higher than initburst_freq_threshold + if all_isis is None: + return None + for isi_counter, isi in enumerate(all_isis): + # Convert to Hz + freq = 1000.0 / isi + if freq < initburst_freq_thresh: + # Threshold reached + break + else: + # Add isi to initburst + last_isi = isi_counter + + if last_isi is None: + # No initburst found + return None + else: + # Get index of second peak of last ISI + last_peak = last_isi + 1 + + # Get time of last peak + last_peak_time = peak_times[last_peak] + + # Determine start of sahp interval + sahp_interval_start = min( + last_peak_time + + initburst_sahp_start, + stim_end) + + # Get next peak, we wont search beyond that + next_peak = last_peak + 1 + + # Determine end of sahp interval + # Add initburst_slow_ahp_max to last peak time + # If next peak or stim_end is earlier, use these + # If no next peak, use stim end + if next_peak < len(peak_times): + next_peak_time = peak_times[next_peak] + + sahp_interval_end = min( + last_peak_time + initburst_sahp_end, next_peak_time, stim_end) + else: + sahp_interval_end = min( + last_peak_time + initburst_sahp_end, stim_end) + + if sahp_interval_end <= sahp_interval_start: + return None + else: + sahp_interval = voltage[np.where( + (time <= sahp_interval_end) & + (time >= sahp_interval_start))] + + if len(sahp_interval) > 0: + min_volt_index = np.argmin(sahp_interval) + else: + return None + + slow_ahp = sahp_interval[min_volt_index] + + return np.array([slow_ahp]) + + +def strict_burst_number() -> np.ndarray: + """Calculate the strict burst number. + + 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. Default value is 2.0.""" + burst_mean_freq = get_cpp_feature("strict_burst_mean_freq") + if burst_mean_freq is None: + return np.array([0]) + return np.array([burst_mean_freq.size]) + + +def inv_ISI_values() -> np.ndarray | None: + """Calculate the inverse of ISI values.""" + all_isi_values_vec = get_cpp_feature("all_ISI_values") + return None if all_isi_values_vec is None else 1000.0 / all_isi_values_vec + + +def inv_first_ISI() -> np.ndarray | None: + """Calculate the inverse of the first ISI.""" + inv_isi_values = inv_ISI_values() + if inv_isi_values is None or len(inv_isi_values) == 0: + return None + return np.array([inv_isi_values[0]]) + + +def inv_second_ISI() -> np.ndarray | None: + """Calculate the inverse of the second ISI.""" + inv_isi_values = inv_ISI_values() + if inv_isi_values is None or len(inv_isi_values) < 2: + return None + return np.array([inv_isi_values[1]]) + + +def inv_third_ISI() -> np.ndarray | None: + """Calculate the inverse of the third ISI.""" + inv_isi_values = inv_ISI_values() + if inv_isi_values is None or len(inv_isi_values) < 3: + return None + return np.array([inv_isi_values[2]]) + + +def inv_fourth_ISI() -> np.ndarray | None: + """Calculate the inverse of the fourth ISI.""" + inv_isi_values = inv_ISI_values() + if inv_isi_values is None or len(inv_isi_values) < 4: + return None + return np.array([inv_isi_values[3]]) + + +def inv_fifth_ISI() -> np.ndarray | None: + """Calculate the inverse of the fifth ISI.""" + inv_isi_values = inv_ISI_values() + if inv_isi_values is None or len(inv_isi_values) < 5: + return None + return np.array([inv_isi_values[4]]) + + +def inv_last_ISI() -> np.ndarray | None: + """Calculate the inverse of the last ISI.""" + inv_isi_values = inv_ISI_values() + if inv_isi_values is None or len(inv_isi_values) == 0: + return None + return np.array([inv_isi_values[-1]]) diff --git a/efel/pyfeatures/pyfeatures.py b/efel/pyfeatures/pyfeatures.py index acdddee5..097f1a8b 100644 --- a/efel/pyfeatures/pyfeatures.py +++ b/efel/pyfeatures/pyfeatures.py @@ -25,12 +25,11 @@ (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. """ - +from efel.pyfeatures.cppfeature_access import _get_cpp_data, get_cpp_feature +from efel.pyfeatures.isi import * from typing_extensions import deprecated -import warnings import numpy as np -from efel import cppcore from numpy.fft import * @@ -39,6 +38,16 @@ 'time', 'current', 'ISIs', + 'ISI_values', + 'ISI_CV', + 'single_burst_ratio', + 'irregularity_index', + 'burst_ISI_indices', + 'burst_mean_freq', + 'interburst_voltage', + 'ISI_log_slope', + 'ISI_semilog_slope', + 'ISI_log_slope_skip', 'initburst_sahp', 'initburst_sahp_vb', 'initburst_sahp_ssse', @@ -57,16 +66,23 @@ 'strict_burst_number', 'trace_check', 'phaseslope_max', + 'inv_ISI_values', + 'inv_first_ISI', + 'inv_second_ISI', + 'inv_third_ISI', + 'inv_fourth_ISI', + 'inv_fifth_ISI', + 'inv_last_ISI' ] -def voltage(): - """Get voltage trace""" +def voltage() -> np.ndarray | None: + """Get voltage trace.""" return get_cpp_feature("voltage") -def time(): - """Get time trace""" +def time() -> np.ndarray | None: + """Get time trace.""" return get_cpp_feature("time") @@ -119,25 +135,8 @@ def trace_check() -> np.ndarray | None: def burst_number() -> np.ndarray: """The number of bursts.""" - burst_mean_freq = get_cpp_feature("burst_mean_freq") - if burst_mean_freq is None: - return np.array([0]) - return np.array([burst_mean_freq.size]) - - -def strict_burst_number() -> np.ndarray: - """Calculate the strict burst number. - - 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. Default value is 2.0.""" - burst_mean_freq = get_cpp_feature("strict_burst_mean_freq") - if burst_mean_freq is None: - return np.array([0]) - return np.array([burst_mean_freq.size]) + mean_freq = burst_mean_freq() + return np.array([0]) if mean_freq is None else np.array([mean_freq.size]) def impedance(): @@ -145,7 +144,7 @@ def impedance(): dt = _get_cpp_data("interp_step") Z_max_freq = _get_cpp_data("impedance_max_freq") - voltage_trace = voltage() + voltage_trace = get_cpp_feature("voltage") holding_voltage = get_cpp_feature("voltage_base") normalized_voltage = voltage_trace - holding_voltage current_trace = current() @@ -179,15 +178,6 @@ def current(): return get_cpp_feature("current") -def ISIs(): - """Get all ISIs.""" - peak_times = get_cpp_feature("peak_time") - if peak_times is None: - return None - else: - return np.diff(peak_times) - - def initburst_sahp_vb(): """SlowAHP voltage from voltage base after initial burst""" @@ -216,90 +206,6 @@ def initburst_sahp_ssse(): return np.array([initburst_sahp_value[0] - ssse[0]]) -def initburst_sahp(): - """SlowAHP voltage after initial burst""" - - # Required cpp features - voltage = get_cpp_feature("voltage") - time = get_cpp_feature("time") - time = time[:len(voltage)] - peak_times = get_cpp_feature("peak_time") - - # Required python features - all_isis = ISIs() - - # Required trace data - stim_end = _get_cpp_data("stim_end") - - # Required settings - initburst_freq_thresh = _get_cpp_data("initburst_freq_threshold") - initburst_sahp_start = _get_cpp_data("initburst_sahp_start") - initburst_sahp_end = _get_cpp_data("initburst_sahp_end") - - last_isi = None - - # Loop over ISIs until frequency higher than initburst_freq_threshold - if all_isis is None: - return None - for isi_counter, isi in enumerate(all_isis): - # Convert to Hz - freq = 1000.0 / isi - if freq < initburst_freq_thresh: - # Threshold reached - break - else: - # Add isi to initburst - last_isi = isi_counter - - if last_isi is None: - # No initburst found - return None - else: - # Get index of second peak of last ISI - last_peak = last_isi + 1 - - # Get time of last peak - last_peak_time = peak_times[last_peak] - - # Determine start of sahp interval - sahp_interval_start = min( - last_peak_time + - initburst_sahp_start, - stim_end) - - # Get next peak, we wont search beyond that - next_peak = last_peak + 1 - - # Determine end of sahp interval - # Add initburst_slow_ahp_max to last peak time - # If next peak or stim_end is earlier, use these - # If no next peak, use stim end - if next_peak < len(peak_times): - next_peak_time = peak_times[next_peak] - - sahp_interval_end = min( - last_peak_time + initburst_sahp_end, next_peak_time, stim_end) - else: - sahp_interval_end = min( - last_peak_time + initburst_sahp_end, stim_end) - - if sahp_interval_end <= sahp_interval_start: - return None - else: - sahp_interval = voltage[np.where( - (time <= sahp_interval_end) & - (time >= sahp_interval_start))] - - if len(sahp_interval) > 0: - min_volt_index = np.argmin(sahp_interval) - else: - return None - - slow_ahp = sahp_interval[min_volt_index] - - return np.array([slow_ahp]) - - def depol_block(): """Check for a depolarization block""" @@ -441,21 +347,3 @@ def phaseslope_max() -> np.ndarray | None: return np.array([np.max(phaseslope)]) except ValueError: return None - - -def get_cpp_feature(feature_name: str, raise_warnings=False) -> np.ndarray | None: - """Return value of feature implemented in cpp.""" - cppcoreFeatureValues: list[int | float] = list() - exitCode = cppcore.getFeature(feature_name, cppcoreFeatureValues) - if exitCode < 0: - if raise_warnings: - warnings.warn( - f"Error while calculating {feature_name}, {cppcore.getgError()}", - RuntimeWarning) - return None - return np.array(cppcoreFeatureValues) - - -def _get_cpp_data(data_name: str) -> float: - """Get cpp data value.""" - return cppcore.getMapDoubleData(data_name)[0] diff --git a/efel/units/units.json b/efel/units/units.json index df18f17c..c80ca33d 100644 --- a/efel/units/units.json +++ b/efel/units/units.json @@ -75,6 +75,7 @@ "fast_AHP_change": "constant", "interburst_min_values": "mV", "interburst_voltage": "mV", + "inv_ISI_values": "Hz", "inv_fifth_ISI": "Hz", "inv_first_ISI": "Hz", "inv_fourth_ISI": "Hz", diff --git a/pyproject.toml b/pyproject.toml index 5d204893..fa5343b7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -8,3 +8,6 @@ exclude_lines = [ "pragma: no cover", "raise NotImplementedError" ] + +[tool.pytest.ini_options] +filterwarnings = ["ignore::DeprecationWarning"] \ No newline at end of file diff --git a/setup.py b/setup.py index 96b978c7..33de406b 100644 --- a/setup.py +++ b/setup.py @@ -78,8 +78,17 @@ name="efel", version=versioneer.get_version(), cmdclass=versioneer.get_cmdclass(), - install_requires=['numpy>=1.6', 'neo>=0.5.2', 'typing-extensions>=4.8.0'], - packages=['efel', 'efel.pyfeatures', 'efel.units'], + install_requires=[ + 'numpy>=1.6', + 'neo>=0.5.2', + 'typing-extensions>=4.8.0', + 'scipy>=1.12.0,<2.0.0', + ], + packages=[ + 'efel', + 'efel.pyfeatures', + 'efel.units' + ], author="BlueBrain Project, EPFL", maintainer="Werner Van Geit", maintainer_email="werner.vangeit@epfl.ch", diff --git a/tests/DependencyV5_LibV5peakindices.txt b/tests/DependencyV5_LibV5peakindices.txt index fe49caac..d632aa2b 100644 --- a/tests/DependencyV5_LibV5peakindices.txt +++ b/tests/DependencyV5_LibV5peakindices.txt @@ -1,8 +1,6 @@ LibV5:peak_indices #LibV1:interpolate -LibV1:ISI_values #LibV1:peak_time LibV1:doublet_ISI #LibV1:peak_time LibV1:peak_voltage #LibV5:peak_indices -LibV1:burst_ISI_indices #LibV5:peak_indices #LibV1:ISI_values LibV1:mean_frequency #LibV1:peak_time LibV1:peak_time #LibV5:peak_indices LibV1:time_to_first_spike #LibV1:peak_time @@ -10,8 +8,6 @@ LibV1:adaptation_index #LibV1:peak_time LibV1:adaptation_index2 #LibV1:peak_time LibV1:spike_width2 #LibV5:min_AHP_indices LibV1:AP_width #LibV5:peak_indices #LibV5:min_AHP_indices -LibV1:burst_mean_freq #LibV1:burst_ISI_indices #LibV1:peak_time -LibV1:interburst_voltage #LibV1:burst_ISI_indices LibV1:AP_height #LibV1:peak_voltage LibV1:AP_amplitude #LibV5:AP_begin_indices #LibV1:peak_voltage LibV1:AHP_depth_abs_slow #LibV5:peak_indices @@ -24,7 +20,6 @@ LibV1:minimum_voltage LibV1:interpolate LibV1:steady_state_voltage LibV3:depolarized_base -LibV1:ISI_CV #LibV1:ISI_values LibV1:AHP_depth #LibV5:voltage_base #LibV5:min_AHP_values LibV2:AP_rise_indices #LibV5:peak_indices #LibV5:AP_begin_indices LibV2:AP_end_indices #LibV5:peak_indices @@ -42,20 +37,16 @@ LibV2:AP_rise_rate_change #LibV2:AP_rise_rate LibV2:AP_fall_rate_change #LibV2:AP_fall_rate LibV2:fast_AHP_change #LibV2:fast_AHP LibV2:AP_duration_half_width_change #LibV2:AP_duration_half_width -LibV1:single_burst_ratio #LibV1:ISI_values LibV2:steady_state_hyper LibV2:amp_drop_first_second #LibV1:peak_voltage LibV2:amp_drop_first_last #LibV1:peak_voltage LibV2:amp_drop_second_last #LibV1:peak_voltage LibV2:max_amp_difference #LibV1:peak_voltage LibV1:AP_amplitude_diff #LibV1:AP_amplitude -LibV5:ISI_log_slope #LibV1:ISI_values -LibV5:ISI_log_slope_skip #LibV1:ISI_values LibV1:AHP_depth_diff #LibV1:AHP_depth LibV5:min_AHP_indices #LibV5:peak_indices LibV5:min_AHP_values #LibV5:min_AHP_indices LibV5:number_initial_spikes #LibV1:peak_time -LibV5:irregularity_index #LibV1:ISI_values LibV5:AP1_amp #LibV1:AP_amplitude LibV5:APlast_amp #LibV1:AP_amplitude LibV5:AP2_amp #LibV1:AP_amplitude @@ -71,12 +62,6 @@ LibV5:AHP1_depth_from_peak #LibV5:AHP_depth_from_peak LibV5:AHP2_depth_from_peak #LibV5:AHP_depth_from_peak LibV5:time_to_second_spike #LibV1:peak_time LibV5:time_to_last_spike #LibV1:peak_time -LibV5:inv_first_ISI #LibV5:all_ISI_values -LibV5:inv_second_ISI #LibV5:all_ISI_values -LibV5:inv_third_ISI #LibV5:all_ISI_values -LibV5:inv_fourth_ISI #LibV5:all_ISI_values -LibV5:inv_fifth_ISI #LibV5:all_ISI_values -LibV5:inv_last_ISI #LibV5:all_ISI_values LibV5:inv_time_to_first_spike #LibV1:time_to_first_spike LibV5:spike_half_width #LibV5:min_AHP_indices #LibV5:peak_indices LibV5:AP_begin_indices #LibV5:min_AHP_indices diff --git a/tests/featurenames.json b/tests/featurenames.json index 0e4425e9..47706097 100644 --- a/tests/featurenames.json +++ b/tests/featurenames.json @@ -81,6 +81,7 @@ "initburst_sahp_ssse", "initburst_sahp_vb", "interburst_voltage", + "inv_ISI_values", "inv_fifth_ISI", "inv_first_ISI", "inv_fourth_ISI", diff --git a/tests/test_basic.py b/tests/test_basic.py index c8889e4b..fb190b38 100644 --- a/tests/test_basic.py +++ b/tests/test_basic.py @@ -2410,37 +2410,6 @@ def test_AP_width_between_threshold_strict(): assert AP_width is None -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 - - def test_burst_mean_freq(): """basic: Test burst_mean_freq""" urls = [burst1_url, burst2_url, burst3_url] @@ -2471,18 +2440,12 @@ def test_burst_mean_freq(): features, raise_warnings=False) burst_mean_freq = feature_values[0]['burst_mean_freq'] - burst_ISI_indices = feature_values[0]['burst_ISI_indices'] - peak_time = feature_values[0]['peak_time'] - - burst_mean_freq_py = py_mean_freq(peak_time, burst_ISI_indices) if i == 1: - assert burst_mean_freq == burst_mean_freq_py assert burst_mean_freq is None elif i == 2: assert burst_mean_freq is None else: - numpy.testing.assert_allclose(burst_mean_freq, burst_mean_freq_py) numpy.testing.assert_allclose(burst_mean_freq, expected_value) @@ -2514,23 +2477,6 @@ def test_segfault_in_AP_begin_width(): feature_values[0]['AP_begin_width'], expected_values) -def py_interburst_voltage(burst_ISI_idxs, peak_idxs, t, v): - """Python implementation of interburst_voltage""" - interburst_voltage = [] - for idx in burst_ISI_idxs: - ts_idx = peak_idxs[idx] - t_start = t[ts_idx] + 5 - start_idx = numpy.argwhere(t < t_start)[-1][0] - - te_idx = peak_idxs[idx + 1] - t_end = t[te_idx] - 5 - end_idx = numpy.argwhere(t > t_end)[0][0] - - interburst_voltage.append(numpy.mean(v[start_idx:end_idx + 1])) - - return numpy.array(interburst_voltage) - - def test_interburst_voltage(): """basic: Test interburst_voltage""" import efel @@ -2553,14 +2499,7 @@ def test_interburst_voltage(): features, raise_warnings=False) interburst_voltage = feature_values[0]['interburst_voltage'] - burst_ISI_indices = feature_values[0]['burst_ISI_indices'] - peak_indices = feature_values[0]['peak_indices'] - - interburst_voltage_py = py_interburst_voltage( - burst_ISI_indices, peak_indices, time, voltage - ) - numpy.testing.assert_allclose(interburst_voltage, interburst_voltage_py) numpy.testing.assert_allclose(interburst_voltage, -63.234682) diff --git a/tests/test_isi.py b/tests/test_isi.py new file mode 100644 index 00000000..6a8c1253 --- /dev/null +++ b/tests/test_isi.py @@ -0,0 +1,156 @@ +from __future__ import annotations +from unittest.mock import patch +import numpy as np +import pytest +import efel +from efel import get_feature_names, get_feature_values + + +def generate_spike_data( + duration=1000, spike_interval=100, spike_duration=1, baseline=-80, spike_value=0 +) -> tuple[np.ndarray, np.ndarray]: + """Generate time and voltage arrays for a signal with spikes.""" + voltage = np.full(duration, baseline) + # Add spikes to the voltage array + for i in range(0, duration, spike_interval): + voltage[i: i + spike_duration] = spike_value + time = np.arange(duration) + return time, voltage + + +def test_ISIs_single_spike(): + """Test the edge case where there are less than 2 spikes.""" + time, voltage = generate_spike_data(spike_interval=600) + trace = { + "T": time, + "V": voltage, + "stim_start": [0], + "stim_end": [1000], + } + + features = ["peak_time", "ISIs"] + feature_values = efel.getFeatureValues([trace], features)[0] + + assert len(feature_values["peak_time"]) == 1 + assert feature_values["ISIs"] is None + + +class TestRegularISI: + @pytest.fixture(autouse=True) + def setup_and_teardown(self): + # setup + efel.reset() + self.time, self.voltage = generate_spike_data() + self.trace = { + 'T': self.time, + 'V': self.voltage, + 'stim_start': [0], + 'stim_end': [1000], + } + self.features = ["single_burst_ratio", "ISIs", "irregularity_index", + "ISI_log_slope", "ISI_semilog_slope", "ISI_log_slope_skip", + "burst_ISI_indices" + ] + self.feature_values = get_feature_values( + [self.trace], + self.features, raise_warnings=False)[0] + # teardown + yield + efel.reset() + + def test_ISIs(self): + assert np.allclose(self.feature_values["ISIs"], 100.0) + + def test_single_burst_ratio(self): + assert "single_burst_ratio" in get_feature_names() + assert self.feature_values["single_burst_ratio"] == pytest.approx(1.0) + + def test_irregularity_index(self): + assert ( + self.feature_values["irregularity_index"] == pytest.approx(0.0, abs=1e-9) + ) + + def test_ISI_log_slope(self): + assert self.feature_values["ISI_log_slope"] == pytest.approx(0.0) + + def test_ISI_semilog_slope(self): + assert self.feature_values["ISI_semilog_slope"] == pytest.approx(0.0) + + def test_ISI_log_slope_skip(self): + assert self.feature_values["ISI_log_slope_skip"] == pytest.approx(0.0) + + def test_ISI_log_slope_skip_ValueError(self): + with pytest.raises(ValueError): + efel.set_double_setting("spike_skipf", 1.0) + get_feature_values( + [self.trace], + ["ISI_log_slope_skip"], raise_warnings=False)[0] + + def test_burst_ISI_indices(self): + assert self.feature_values["burst_ISI_indices"] is None + + +class TestThreeSpikes: + @pytest.fixture(autouse=True) + def setup_and_teardown(self): + # setup + efel.reset() + duration = 100 + spike_interval = 25 + self.time, self.voltage = generate_spike_data(duration, spike_interval) + self.trace = { + 'T': self.time, + 'V': self.voltage, + 'stim_start': [0], + 'stim_end': [duration], + } + self.features = ["single_burst_ratio", "ISIs", "peak_time"] + self.feature_values = get_feature_values( + [self.trace], + self.features, raise_warnings=False)[0] + # teardown + yield + efel.reset() + + def test_ISIs(self): + # duration/spike_interval is 4 but the beginning of the voltage + # trace is not a spike + assert len(self.feature_values["peak_time"]) == 3 + assert len(self.feature_values["ISIs"]) == 2 + assert np.allclose(self.feature_values["ISIs"], 25.0) + + def test_single_burst_ratio(self): + # none due to default ignore_first_ISI=True + assert self.feature_values["single_burst_ratio"] is None + + # set ignore_first_ISI=False + efel.set_int_setting("ignore_first_ISI", 0) + self.feature_values = get_feature_values( + [self.trace], + self.features, raise_warnings=False)[0] + assert self.feature_values["single_burst_ratio"] == pytest.approx(1.0) + + +class TestBurstISIIndices: + @pytest.fixture(autouse=True) + def setup_and_teardown(self): + efel.reset() + self.trace = { # bypassed values since ISI_values are defined in the mock + "T": [1, 2, 3], + "V": [4, 5, 6], + "stim_start": [0], + "stim_end": [1], + } + self.features = ["burst_ISI_indices"] + yield + efel.reset() + + @patch('efel.pyfeatures.isi.ISI_values') + def test_burst_ISI_indices(self, mock_ISI_values): + mock_ISI_values.return_value = np.array([50, 100, 50, 200, 50]) + + self.feature_values = efel.getFeatureValues( + [self.trace], self.features, raise_warnings=False + )[0] + + assert self.feature_values["burst_ISI_indices"].tolist() == [2, 4] diff --git a/tests/testdata/allfeatures/expectedresults.json b/tests/testdata/allfeatures/expectedresults.json index d7d7575b..06db917d 100644 --- a/tests/testdata/allfeatures/expectedresults.json +++ b/tests/testdata/allfeatures/expectedresults.json @@ -382,6 +382,13 @@ "interburst_voltage": [ -38.820613862535545 ], + "inv_ISI_values": [ + 4.918839153958547, + 2.021427127553367, + 3.2679738562121226, + 1.4803849000753657, + 3.995205753099918 + ], "inv_fifth_ISI": [ 3.995205753099918 ], @@ -404,7 +411,7 @@ 124.99999999855582 ], "irregularity_index": [ - 245.8499999998084 + 327.79999999974456 ], "is_not_stuck": [ 1 @@ -486,7 +493,7 @@ ], "sag_time_constant": null, "single_burst_ratio": [ - 1.204822211398266 + 1.1461337966987346 ], "spike_half_width": [ 1.563889523951499, diff --git a/tox.ini b/tox.ini index 110169ae..f458d4a0 100644 --- a/tox.ini +++ b/tox.ini @@ -3,7 +3,6 @@ envlist = docs,lint,py3-{test} minversion = 4 [gh-actions] python = - 3.8: py3 3.9: py3 3.10: py3 3.11: py3 @@ -14,7 +13,6 @@ python = envdir = {toxworkdir}/py3-test deps = pytest>=7.3.1 - scipy>=1.10.1 pytest-xdist>=3.3.1 extras = neo @@ -30,7 +28,6 @@ setenv = deps = setuptools>=69.0.3 pytest>=7.3.1 - scipy>=1.10.1 pytest-cov>=4.1.0 coverage>=7.3.0 gcovr @@ -52,7 +49,6 @@ commands = envdir={toxworkdir}/{envname} deps = pytest>=7.3.1 - scipy>=1.10.1 nbmake deap matplotlib