diff --git a/efel/cppcore/CMakeLists.txt b/efel/cppcore/CMakeLists.txt index 8743f97f..3803c306 100644 --- a/efel/cppcore/CMakeLists.txt +++ b/efel/cppcore/CMakeLists.txt @@ -19,7 +19,7 @@ cmake_minimum_required(VERSION 2.6) set(CMAKE_BUILD_TYPE Debug) -set(FEATURESRCS Utils.cpp LibV1.cpp LibV2.cpp LibV3.cpp LibV4.cpp LibV5.cpp +set(FEATURESRCS Utils.cpp LibV1.cpp LibV2.cpp LibV3.cpp LibV5.cpp FillFptrTable.cpp DependencyTree.cpp efel.cpp cfeature.cpp mapoperations.cpp) @@ -33,6 +33,6 @@ add_library(efel SHARED ${FEATURESRCS}) install(TARGETS efel LIBRARY DESTINATION lib) install(FILES efel.h cfeature.h FillFptrTable.h LibV1.h LibV2.h LibV3.h - LibV4.h LibV5.h mapoperations.h Utils.h DependencyTree.h eFELLogger.h + LibV5.h mapoperations.h Utils.h DependencyTree.h eFELLogger.h types.h DESTINATION include) diff --git a/efel/cppcore/FillFptrTable.cpp b/efel/cppcore/FillFptrTable.cpp index a31acd50..284d16ab 100644 --- a/efel/cppcore/FillFptrTable.cpp +++ b/efel/cppcore/FillFptrTable.cpp @@ -20,23 +20,17 @@ int FillFptrTable() { //****************** for FptrTableV1 ***************************** FptrTableV1["interpolate"] = &LibV1::interpolate; - FptrTableV1["peak_indices"] = &LibV1::peak_indices; 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["min_AHP_indices"] = &LibV1::min_AHP_indices; - FptrTableV1["min_AHP_values"] = &LibV1::min_AHP_values; - FptrTableV1["voltage_base"] = &LibV1::rest_voltage_value; FptrTableV1["burst_ISI_indices"] = &LibV1::burst_ISI_indices; FptrTableV1["adaptation_index"] = &LibV1::adaptation_index; FptrTableV1["trace_check"] = &LibV1::trace_check; - FptrTableV1["spike_half_width"] = &LibV1::spike_width1; FptrTableV1["spike_width2"] = &LibV1::spike_width2; FptrTableV1["burst_mean_freq"] = &LibV1::burst_mean_freq; FptrTableV1["interburst_voltage"] = &LibV1::interburst_voltage; - FptrTableV1["AHP_depth_abs"] = &LibV1::AHP_depth_abs; // passive properties FptrTableV1["time_constant"] = &LibV1::time_constant; FptrTableV1["voltage_deflection"] = &LibV1::voltage_deflection; @@ -65,17 +59,9 @@ int FillFptrTable() { //****************** for FptrTableV2 ***************************** /* - FptrTableV2["peak_indices"] = &LibV2::peak_indices; - FptrTableV2["ISI_values"] = &LibV2::ISI_values; - FptrTableV2["peak_voltage"] = &LibV2::peak_voltage; - FptrTableV2["firing_rate"] = &LibV2::firing_rate; - FptrTableV2["peak_time"] = &LibV2::peak_time; - FptrTableV2["first_spike_time"] = &LibV2::first_spike_time; FptrTableV2["AHP_min_indices"] = &LibV2::AHP_min_indices; FptrTableV2["AHP_values"] = &LibV2::AHP_values; - FptrTableV2["rest_voltage_value"] = &LibV2::rest_voltage_value; FptrTableV2["burst_vector"] = &LibV2::burst_vector; - FptrTableV2["adaptation_index"] = &LibV2::adaptation_index; */ // AP parameter @@ -137,18 +123,11 @@ int FillFptrTable() { //****************** end of FptrTableV2 ***************************** //****************** FptrTableV3 ***************************** - FptrTableV3["AP_amplitude"] = &LibV3::AP_amplitude; - FptrTableV3["AP_width"] = &LibV3::AP_width; - - FptrTableV3["AP_begin_indices"] = &LibV3::AP_begin_indices; - FptrTableV3["AP_end_indices"] = &LibV3::AP_end_indices; // eFeatures FptrTableV3["depolarized_base"] = &LibV3::depolarized_base; //****************** end of FptrTableV3 ***************************** - FptrTableV4["peak_indices"] = &LibV4::peak_indices; - //****************** FptrTableV5 ***************************** FptrTableV5["ISI_log_slope"] = &LibV5::ISI_log_slope; diff --git a/efel/cppcore/FillFptrTable.h b/efel/cppcore/FillFptrTable.h index 1dac33b6..7817c07a 100644 --- a/efel/cppcore/FillFptrTable.h +++ b/efel/cppcore/FillFptrTable.h @@ -24,13 +24,11 @@ #include "LibV1.h" #include "LibV2.h" #include "LibV3.h" -#include "LibV4.h" #include "LibV5.h" extern feature2function FptrTableV1; extern feature2function FptrTableV2; extern feature2function FptrTableV3; -extern feature2function FptrTableV4; extern feature2function FptrTableV5; int FillFptrTable(); diff --git a/efel/cppcore/LibV1.cpp b/efel/cppcore/LibV1.cpp index a3952534..7da591a6 100644 --- a/efel/cppcore/LibV1.cpp +++ b/efel/cppcore/LibV1.cpp @@ -85,66 +85,6 @@ int LibV1::interpolate(mapStr2intVec& IntFeatureData, return retVal; } -static int __peak_indices(double dThreshold, vector& V, - vector& PeakIndex) { - vector upVec, dnVec; - double dtmp; - int itmp; - - for (size_t i = 1; i < V.size(); i++) { - if (V[i] > dThreshold && V[i - 1] < dThreshold) { - upVec.push_back(i); - } else if (V[i] < dThreshold && V[i - 1] > dThreshold) { - dnVec.push_back(i); - } - } - if (dnVec.size() == 0) { - GErrorStr += - "\nVoltage never goes below or above threshold in spike detection.\n"; - return 0; - } - - if (dnVec.size() != upVec.size()) { - GErrorStr += "\nVoltage never goes below threshold after last spike.\n"; - return 0; - } - - PeakIndex.clear(); - int j = 0; - for (size_t i = 0; i < upVec.size(); i++) { - dtmp = -1e9; - itmp = -1; - for (j = upVec[i]; j <= dnVec[i]; j++) { - if (dtmp < V[j]) { - dtmp = V[j]; - itmp = j; - } - } - if (itmp != -1) PeakIndex.push_back(itmp); - } - return PeakIndex.size(); -} -int LibV1::peak_indices(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - // printf("\n LibV1 This is inside peak_indices()"); - int retVal, nSize; - retVal = - CheckInMap(IntFeatureData, StringData, "peak_indices", nSize); - if (retVal) - return nSize; - vector PeakIndex; - vector v, Th; - retVal = getVec(DoubleFeatureData, StringData, "V", v); - if (retVal <= 0) return -1; - retVal = getDoubleParam(DoubleFeatureData, "Threshold", Th); - if (retVal <= 0) return -1; - int retval = __peak_indices(Th[0], v, PeakIndex); - if (retval >= 0) - setVec(IntFeatureData, StringData, "peak_indices", PeakIndex); - return retval; -} - // *** Spikecount *** int LibV1::Spikecount(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, @@ -358,67 +298,6 @@ int LibV1::first_spike_time(mapStr2intVec& IntFeatureData, return first_spike.size(); } -// min_AHP_indices -// find the minimum between two spikes, -// and the minimum between the last spike and the time the stimulus ends -int LibV1::min_AHP_indices(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - int retVal, nSize; - retVal = CheckInMap(IntFeatureData, StringData, "min_AHP_indices", nSize); - if (retVal) - return nSize; - - vector peak_indices_plus; - vector min_ahp_indices; - vector v; - vector min_ahp_values; - vector stim_end; - vector t; - retVal = getVec(DoubleFeatureData, StringData, "V", v); - if (retVal <= 0) return -1; - retVal = getVec(IntFeatureData, StringData, "peak_indices", - peak_indices_plus); - if (retVal < 1) { - GErrorStr += "\n At least one spike required for calculation of " - "min_AHP_indices.\n"; - return -1; - } - retVal = getVec(DoubleFeatureData, StringData, "stim_end", stim_end); - if (retVal <= 0) return -1; - retVal = getVec(DoubleFeatureData, StringData, "T", t); - if (retVal <= 0) return -1; - - int end_index = distance( - t.begin(), find_if(t.begin(), t.end(), - bind2nd(greater_equal(), stim_end[0]))); - // if the last spike happens to be close to the end of the stimulus - // there will not be a proper AHP, this case is not properly dealt with here - if (end_index > peak_indices_plus.back() + 5) { - peak_indices_plus.push_back(end_index); - } - for (size_t i = 0; i < peak_indices_plus.size() - 1; i++) { - int ahpindex = distance( - v.begin(), min_element(v.begin() + peak_indices_plus[i], - v.begin() + peak_indices_plus[i + 1])); - min_ahp_indices.push_back(ahpindex); - min_ahp_values.push_back(v[ahpindex]); - } - setVec(IntFeatureData, StringData, "min_AHP_indices", min_ahp_indices); - setVec(DoubleFeatureData, StringData, "min_AHP_values", min_ahp_values); - return min_ahp_indices.size(); -} - -int LibV1::min_AHP_values(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - int retVal, nSize; - retVal = CheckInMap(DoubleFeatureData, StringData, - "min_AHP_values", nSize); - if (retVal) return nSize; - return -1; -} - int LibV1::AP_height(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData) { @@ -517,29 +396,6 @@ int LibV1::AP_amplitude(mapStr2intVec& IntFeatureData, return apamplitude.size(); } -// AHP_depth_abs -// naming conflict here AHP_depth_abs does the same as min_AHP_values. -// In my opinion the should not be a feature called 'AHP_depth_abs', -// use min_AHP_values instead. -// A more interesting feature would be 'AHP_depth' anyways, which calculates the -// depth of the AHP relative to the voltage base -int LibV1::AHP_depth_abs(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - int retVal, nSize; - retVal = CheckInMap(DoubleFeatureData, StringData, - "AHP_depth_abs", nSize); - if (retVal) - return nSize; - - vector vAHP; - retVal = getVec(DoubleFeatureData, StringData, - "min_AHP_values", vAHP); - if (retVal <= 0) return -1; - setVec(DoubleFeatureData, StringData, "AHP_depth_abs", vAHP); - return vAHP.size(); -} - // *** AHP_depth_abs_slow *** // same as AHP_depth_abs but the minimum search starts // 5 ms (or custom duration) after the spike, @@ -625,43 +481,6 @@ int LibV1::AHP_slow_time(mapStr2intVec& IntFeatureData, return -1; } -int LibV1::rest_voltage_value(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - int retVal, nSize; - retVal = CheckInMap(DoubleFeatureData, StringData, - "voltage_base", nSize); - if (retVal) - return nSize; - - vector v, t, stimStart, vRest; - double startTime, endTime; - retVal = getVec(DoubleFeatureData, StringData, "V", v); - if (retVal < 0) return -1; - retVal = getVec(DoubleFeatureData, StringData, "T", t); - if (retVal < 0) return -1; - retVal = getVec(DoubleFeatureData, StringData, "stim_start", stimStart); - if (retVal < 0) return -1; - startTime = stimStart[0] * .25; // It is 25% from start (0), so if stimulus - // starts at 100ms then StartTime will be - // 25mS - // as in above case end time will be 25% less than startTime - endTime = stimStart[0] * .75; - int nCount = 0; - double vSum = 0; - // calculte the mean of voltage between startTime and endTime - for (size_t i = 0; i < t.size(); i++) { - if (t[i] >= startTime) { - vSum = vSum + v[i]; - nCount++; - } - if (t[i] > endTime) break; - } - vRest.push_back(vSum / nCount); - setVec(DoubleFeatureData, StringData, "voltage_base", vRest); - return 1; -} - static int __burst_ISI_indices(double BurstFactor, vector& PeakIndex, vector& ISIValues, vector& BurstIndex) { @@ -1220,88 +1039,6 @@ int LibV1::spike_width2(mapStr2intVec& IntFeatureData, return retVal; } -// spike half width -// for spike amplitude = v_peak - v_AHP -static int __spike_width1(const vector& t, const vector& v, - const vector& peak_indices, - const vector& min_ahp_indices, double stim_start, - vector& spike_width1) { - int start_index = distance( - t.begin(), find_if(t.begin(), t.end(), - bind2nd(greater_equal(), stim_start))); - vector min_ahp_indices_plus(min_ahp_indices.size() + 1, start_index); - copy(min_ahp_indices.begin(), min_ahp_indices.end(), - min_ahp_indices_plus.begin() + 1); - for (size_t i = 1; i < min_ahp_indices_plus.size(); i++) { - double v_half = (v[peak_indices[i - 1]] + v[min_ahp_indices_plus[i]]) / 2.; - // interpolate this one time step where the voltage is close to v_half in - // the rising and in the falling edge - double v_dev; - double delta_v; - double t_dev_rise; - double t_dev_fall; - double delta_t; - int rise_index = - distance(v.begin(), find_if(v.begin() + min_ahp_indices_plus[i - 1], - v.begin() + peak_indices[i - 1], - bind2nd(greater_equal(), v_half))); - v_dev = v_half - v[rise_index]; - delta_v = v[rise_index] - v[rise_index - 1]; - delta_t = t[rise_index] - t[rise_index - 1]; - t_dev_rise = delta_t * v_dev / delta_v; - int fall_index = - distance(v.begin(), find_if(v.begin() + peak_indices[i - 1], - v.begin() + min_ahp_indices_plus[i], - bind2nd(less_equal(), v_half))); - v_dev = v_half - v[fall_index]; - delta_v = v[fall_index] - v[fall_index - 1]; - delta_t = t[fall_index] - t[fall_index - 1]; - t_dev_fall = delta_t * v_dev / delta_v; - spike_width1.push_back(t[fall_index] + t_dev_rise - t[rise_index] + - t_dev_fall); - } - return spike_width1.size(); -} - -int LibV1::spike_width1(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - int retVal, nSize; - retVal = CheckInMap(DoubleFeatureData, StringData, - "spike_half_width", nSize); - if (retVal) - return nSize; - - vector PeakIndex, minAHPIndex; - vector V, t, dv1, dv2, spike_width1; - vector stim_start; - retVal = getVec(DoubleFeatureData, StringData, "V", V); - if (retVal < 0) return -1; - retVal = getVec(DoubleFeatureData, StringData, "T", t); - if (retVal < 0) return -1; - retVal = getVec(DoubleFeatureData, StringData, "stim_start", stim_start); - if (retVal < 0) return -1; - retVal = getVec(IntFeatureData, StringData, "min_AHP_indices", minAHPIndex); - if (retVal < 0) return -1; - retVal = getVec(IntFeatureData, StringData, "peak_indices", PeakIndex); - if (retVal < 0) return -1; - if (PeakIndex.size() <= 1) { - GErrorStr += "\nError: More than one spike is needed for spikewidth " - "calculation.\n"; - return -1; - } - // - // Take derivative of voltage from 1st AHPmin to the peak of the spike - // Using Central difference derivative vec1[i] = ((vec[i+1]+vec[i-1])/2)/dx - retVal = __spike_width1(t, V, PeakIndex, minAHPIndex, stim_start[0], - spike_width1); - if (retVal >= 0) { - setVec(DoubleFeatureData, StringData, "spike_half_width", - spike_width1); - } - return retVal; -} - // passive properties implementation // // *** timeconstant *** diff --git a/efel/cppcore/LibV1.h b/efel/cppcore/LibV1.h index a13ce9bf..4443dcf7 100644 --- a/efel/cppcore/LibV1.h +++ b/efel/cppcore/LibV1.h @@ -33,9 +33,6 @@ int printVectorD(char* strName, vector vec); int interpolate(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); -int peak_indices(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); - int ISI_values(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); @@ -52,26 +49,9 @@ int first_spike_time(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); -int spike_width1(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); - int spike_width2(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); -int min_AHP_indices(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData); - -int AHP_depth_abs(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); - -int min_AHP_values(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); - -int rest_voltage_value(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData); - int burst_ISI_indices(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); diff --git a/efel/cppcore/LibV3.cpp b/efel/cppcore/LibV3.cpp index c088c0c1..06fe5b9a 100644 --- a/efel/cppcore/LibV3.cpp +++ b/efel/cppcore/LibV3.cpp @@ -31,276 +31,6 @@ using std::list; using std::min_element; -// spike amplitude: peak_voltage - v[AP_begin_indices] -int LibV3::AP_amplitude(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - int retVal, nSize; - retVal = - CheckInMap(DoubleFeatureData, StringData, "AP_amplitude", nSize); - if (retVal > 0) - return nSize; - - vector peakvoltage; - vector apbeginindices; - vector v; - retVal = getVec(DoubleFeatureData, StringData, "V", v); - if (retVal <= 0) { return -1; } - retVal = getVec(DoubleFeatureData, StringData, "peak_voltage", - peakvoltage); - if (retVal <= 0) { return -1; } - retVal = getVec(IntFeatureData, StringData, "AP_begin_indices", - apbeginindices); - if (retVal <= 0) { return -1; } - - vector apamplitude; - apamplitude.resize(peakvoltage.size()); - for (size_t i = 0; i < apamplitude.size(); i++) { - apamplitude[i] = peakvoltage[i] - v[apbeginindices[i]]; - } - setVec(DoubleFeatureData, StringData, "AP_amplitude", apamplitude); - return apamplitude.size(); -} - -// *** AP_width *** -// -// spike width calculation according to threshold value for the first spike -// unfortunately spike width means the width of the spike on onset not at half -// maximum -static int __AP_width(const vector& t, const vector& v, - double stimstart, double threshold, - const vector& peakindices, - const vector& minahpindices, - vector& apwidth) { - // printf("\n Inside AP_width...\n"); - // printVectorD("t", t); - // printVectorD("v", v); - // printVectorI("peakindices", peakindices); - // printVectorI("minahpindices", minahpindices); - // printf("\nStimStart = %f , thereshold = %f ", stimstart, threshold); - vector indices(minahpindices.size() + 1); - int start_index = distance( - t.begin(), - find_if(t.begin(), t.end(), bind2nd(greater_equal(), stimstart))); - indices[0] = start_index; - copy(minahpindices.begin(), minahpindices.end(), indices.begin() + 1); - for (size_t i = 0; i < indices.size() - 1; i++) { - /* - // FWHM (not used): - // half maximum - double v_hm = (v[peakindices[i]] + threshold) / 2.; - // half maximum indices - int hm_index1 = distance(v.begin(), find_if(v.begin() + indices[i], - v.begin() + indices[i + 1], bind2nd(greater_equal(), v_hm))); - int hm_index2 = distance(v.begin(), find_if(v.begin() + peakindices[i], - v.begin() + indices[i + 1], bind2nd(less_equal(), v_hm))); - apwidth.push_back(t[hm_index2] - t[hm_index1]); - */ - int onset_index = distance( - v.begin(), find_if(v.begin() + indices[i], v.begin() + indices[i + 1], - bind2nd(greater_equal(), threshold))); - // int end_index = distance(v.begin(), find_if(v.begin() + peakindices[i], - // v.begin() + indices[i + 1], bind2nd(less_equal(), threshold))); - int end_index = distance( - v.begin(), find_if(v.begin() + onset_index, v.begin() + indices[i + 1], - bind2nd(less_equal(), threshold))); - apwidth.push_back(t[end_index] - t[onset_index]); - } - return apwidth.size(); -} - -int LibV3::AP_width(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - int retval; - int nsize; - retval = CheckInMap(DoubleFeatureData, StringData, "AP_width", - nsize); - if (retval) { - return nsize; - } - vector t; - retval = getVec(DoubleFeatureData, StringData, "T", t); - if (retval < 0) return -1; - vector v; - retval = getVec(DoubleFeatureData, StringData, "V", v); - if (retval < 0) return -1; - vector threshold; - retval = getDoubleParam(DoubleFeatureData, "Threshold", threshold); - if (retval < 0) return -1; - vector stimstart; - retval = getVec(DoubleFeatureData, StringData, "stim_start", stimstart); - if (retval < 0) return -1; - vector peakindices; - retval = getVec(IntFeatureData, StringData, "peak_indices", peakindices); - if (retval <= 0) { - GErrorStr += "\nNo spike in trace.\n"; - return -1; - } - vector minahpindices; - retval = getVec(IntFeatureData, StringData, "min_AHP_indices", minahpindices); - if (retval < 0) return -1; - vector apwidth; - retval = __AP_width(t, v, stimstart[0], threshold[0], peakindices, - minahpindices, apwidth); - if (retval >= 0) { - setVec(DoubleFeatureData, StringData, "AP_width", apwidth); - } - return retval; -} -// end of AP_width - -// AP parameters -// -// *** AP begin indices *** -// -static int __AP_begin_indices(const vector& t, const vector& v, - double stimstart, double stimend, - const vector& ahpi, vector& apbi) { - // derivative at peak start according to eCode specification 10mV/ms - // according to Shaul 12mV/ms - const double derivativethreshold = 12.; - vector dvdt(v.size()); - vector dv; - vector dt; - getCentralDifferenceDerivative(1., v, dv); - getCentralDifferenceDerivative(1., t, dt); - transform(dv.begin(), dv.end(), dt.begin(), dvdt.begin(), std::divides()); - - // restrict to time interval where stimulus is applied - vector minima; - int stimbeginindex = distance( - t.begin(), - find_if(t.begin(), t.end(), bind2nd(greater_equal(), stimstart))); - minima.push_back(stimbeginindex); - for (size_t i = 0; i < ahpi.size(); i++) { - if (ahpi[i] > stimbeginindex) { - minima.push_back(ahpi[i]); - } - if (t[ahpi[i]] > stimend) { - break; - } - } - // if the AHP_indices are already restricted make sure that we do not miss - // the last spike - if (t[minima.back()] < stimend) { - int stimendindex = - distance(t.begin(), find_if(t.begin() + minima.back(), t.end(), - bind2nd(greater_equal(), stimend))); - minima.push_back(stimendindex); - } - for (size_t i = 0; i < minima.size() - 1; i++) { - // assure that the width of the slope is bigger than 4 - int newbegin = minima[i]; - int begin = minima[i]; - int width = 5; - bool skip = false; - do { - begin = distance( - dvdt.begin(), - find_if(dvdt.begin() + newbegin, dvdt.begin() + minima[i + 1], - bind2nd(greater_equal(), derivativethreshold))); - if (begin == minima[i + 1]) { - // could not find a spike in between these minima - skip = true; - break; - } - newbegin = begin + 1; - } while (find_if(dvdt.begin() + begin, dvdt.begin() + begin + width, - bind2nd(std::less(), derivativethreshold)) != - dvdt.begin() + begin + width); - if (skip) { - continue; - } - apbi.push_back(begin); - } - return apbi.size(); -} -int LibV3::AP_begin_indices(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - int retVal; - int nSize; - retVal = CheckInMap(IntFeatureData, StringData, "AP_begin_indices", - nSize); - if (retVal) { - return nSize; - } - vector t; - retVal = getVec(DoubleFeatureData, StringData, "T", t); - if (retVal < 0) return -1; - vector v; - retVal = getVec(DoubleFeatureData, StringData, "V", v); - if (retVal < 0) return -1; - vector stimstart; - retVal = getVec(DoubleFeatureData, StringData, "stim_start", stimstart); - if (retVal < 0) return -1; - vector stimend; - retVal = getVec(DoubleFeatureData, StringData, "stim_end", stimend); - if (retVal < 0) return -1; - vector ahpi; - retVal = getVec(IntFeatureData, StringData, "min_AHP_indices", ahpi); - if (retVal < 0) return -1; - vector apbi; - retVal = __AP_begin_indices(t, v, stimstart[0], stimend[0], ahpi, apbi); - if (retVal >= 0) { - setVec(IntFeatureData, StringData, "AP_begin_indices", apbi); - } - return retVal; -} - -// *** AP end indices *** -// -static int __AP_end_indices(const vector& t, const vector& v, - const vector& pi, vector& apei) { - // derivative at peak end according to eCode specification -10mV/ms - // according to Shaul -12mV/ms - const double derivativethreshold = -12.; - // assume constant time steps - double timestep = t[1] - t[0]; - vector dvdt; - getCentralDifferenceDerivative(timestep, v, dvdt); - - apei.resize(pi.size()); - vector picopy(pi.begin(), pi.end()); - picopy.push_back(v.size() - 1); - for (size_t i = 0; i < apei.size(); i++) { - // assure that the width of the slope is bigger than 4 - apei[i] = distance( - dvdt.begin(), - find_if(dvdt.begin() + picopy[i] + 1, dvdt.begin() + picopy[i + 1], - bind2nd(greater_equal(), derivativethreshold))); - } - return apei.size(); -} -int LibV3::AP_end_indices(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - int retVal; - int nSize; - retVal = CheckInMap(IntFeatureData, StringData, "AP_end_indices", - nSize); - if (retVal) { - return nSize; - } - - vector t; - retVal = getVec(DoubleFeatureData, StringData, "T", t); - if (retVal < 0) return -1; - vector v; - retVal = getVec(DoubleFeatureData, StringData, "V", v); - if (retVal < 0) return -1; - vector pi; - retVal = getVec(IntFeatureData, StringData, "peak_indices", pi); - if (retVal < 0) return -1; - vector apei; - retVal = __AP_end_indices(t, v, pi, apei); - if (retVal >= 0) { - setVec(IntFeatureData, StringData, "AP_end_indices", apei); - } - return retVal; -} - static int __depolarized_base(const vector& t, const vector& v, double stimstart, double stimend, const vector& apbi, diff --git a/efel/cppcore/LibV3.h b/efel/cppcore/LibV3.h index d14113c8..16e4daf1 100644 --- a/efel/cppcore/LibV3.h +++ b/efel/cppcore/LibV3.h @@ -27,20 +27,6 @@ using std::vector; namespace LibV3 { -int AP_amplitude(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); - -int AP_width(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData); - -int AP_begin_indices(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData); - -int AP_end_indices(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); - -// eFeatures int depolarized_base(mapStr2intVec& IntFeatureData, mapStr2doubleVec& DoubleFeatureData, diff --git a/efel/cppcore/LibV4.cpp b/efel/cppcore/LibV4.cpp deleted file mode 100644 index 3e857a44..00000000 --- a/efel/cppcore/LibV4.cpp +++ /dev/null @@ -1,97 +0,0 @@ -/* Copyright (c) 2015, EPFL/Blue Brain Project - * - * This file is part of eFEL - * - * This library is free software; you can redistribute it and/or modify it under - * the terms of the GNU Lesser General Public License version 3.0 as published - * by the Free Software Foundation. - * - * This library is distributed in the hope that it will be useful, but WITHOUT - * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS - * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more - * details. - * - * You should have received a copy of the GNU Lesser General Public License - * along with this library; if not, write to the Free Software Foundation, Inc., - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - */ - -#include "LibV4.h" -#include -#include -#include -#include - -// new algorithm to find the spikes in the voltage trace -// based on minimum spike height -// -// get all local minima -// assume spikes between local minima -// check if the assumed spike is bigger than min_spike_height -static int __peak_indices(const vector& v, double min_spike_height, - double threshold, vector& peakindices) { - vector dv; - vector minimum_indices; - getCentralDifferenceDerivative(1., v, dv); - int falling_index = 0; - int minimum_index = 0; - do { - minimum_indices.push_back(minimum_index); - falling_index = - distance(dv.begin(), find_if(dv.begin() + minimum_index, dv.end(), - std::bind2nd(std::less(), 0.))); - minimum_index = - distance(dv.begin(), find_if(dv.begin() + falling_index, dv.end(), - std::bind2nd(std::greater_equal(), 0.))); - } while (dv.begin() + minimum_index != dv.end()); - minimum_indices.push_back(dv.size() - 1); - - for (unsigned i = 0; i < minimum_indices.size() - 1; i++) { - int maximum_index = - distance(v.begin(), max_element(v.begin() + minimum_indices[i], - v.begin() + minimum_indices[i + 1])); - double spike_height1 = v[maximum_index] - v[minimum_indices[i + 1]]; - double spike_height2 = v[maximum_index] - v[minimum_indices[i]]; - if ((v[maximum_index] > threshold && (spike_height1 > min_spike_height || - spike_height2 > min_spike_height)) || - (spike_height1 > min_spike_height && - spike_height2 > min_spike_height)) { - // if(spike_height1 > min_spike_height && spike_height2 > - // min_spike_height) { - peakindices.push_back(maximum_index); - } - } - return peakindices.size(); -} - -int LibV4::peak_indices(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, - mapStr2Str& StringData) { - int size; - if (CheckInMap(IntFeatureData, StringData, "peak_indices", size)) { - return size; - } - - vector peakindices; - vector v; - vector min_spike_height; - vector threshold; - if (getVec(DoubleFeatureData, StringData, "V", v) <= 0) { - return -1; - } - if (getDoubleParam(DoubleFeatureData, "min_spike_height", min_spike_height) <= - 0) { - return -1; - } - if (getDoubleParam(DoubleFeatureData, "Threshold", threshold) <= 0) { - return -1; - } - - int retval = __peak_indices(v, min_spike_height[0], threshold[0], peakindices); - if (retval >= 0) { - setVec(IntFeatureData, StringData, "peak_indices", peakindices); - return peakindices.size(); - } - - return retval; -} diff --git a/efel/cppcore/LibV4.h b/efel/cppcore/LibV4.h deleted file mode 100644 index 11ed36a8..00000000 --- a/efel/cppcore/LibV4.h +++ /dev/null @@ -1,33 +0,0 @@ -/* Copyright (c) 2015, EPFL/Blue Brain Project - * - * This file is part of eFEL - * - * This library is free software; you can redistribute it and/or modify it under - * the terms of the GNU Lesser General Public License version 3.0 as published - * by the Free Software Foundation. - * - * This library is distributed in the hope that it will be useful, but WITHOUT - * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS - * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more - * details. - * - * You should have received a copy of the GNU Lesser General Public License - * along with this library; if not, write to the Free Software Foundation, Inc., - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - */ - -#ifndef __LIBV4 -#define __LIBV4 - -#include "mapoperations.h" -#include "Utils.h" - -#include - -using std::vector; - -namespace LibV4 { -int peak_indices(mapStr2intVec& IntFeatureData, - mapStr2doubleVec& DoubleFeatureData, mapStr2Str& StringData); -} -#endif diff --git a/efel/cppcore/cfeature.cpp b/efel/cppcore/cfeature.cpp index 5a2d9e69..985fe580 100644 --- a/efel/cppcore/cfeature.cpp +++ b/efel/cppcore/cfeature.cpp @@ -36,7 +36,6 @@ cFeature::cFeature(const string& strDepFile, const string& outdir) mapFptrLib["LibV1"] = &FptrTableV1; mapFptrLib["LibV2"] = &FptrTableV2; mapFptrLib["LibV3"] = &FptrTableV3; - mapFptrLib["LibV4"] = &FptrTableV4; mapFptrLib["LibV5"] = &FptrTableV5; fillfeaturetypes(); diff --git a/setup.py b/setup.py index a66cfae7..95fcf704 100644 --- a/setup.py +++ b/setup.py @@ -33,7 +33,6 @@ 'LibV1.cpp', 'LibV2.cpp', 'LibV3.cpp', - 'LibV4.cpp', 'LibV5.cpp', 'FillFptrTable.cpp', 'DependencyTree.cpp', @@ -44,7 +43,6 @@ 'LibV1.h', 'LibV2.h', 'LibV3.h', - 'LibV4.h', 'LibV5.h', 'FillFptrTable.h', 'DependencyTree.h', diff --git a/tests/DependencyV5_LibV4peakindices.txt b/tests/DependencyV5_LibV5peakindices.txt similarity index 85% rename from tests/DependencyV5_LibV4peakindices.txt rename to tests/DependencyV5_LibV5peakindices.txt index 9309e755..7933ea70 100644 --- a/tests/DependencyV5_LibV4peakindices.txt +++ b/tests/DependencyV5_LibV5peakindices.txt @@ -1,22 +1,22 @@ -LibV1:peak_indices #LibV1:interpolate +LibV5:peak_indices #LibV1:interpolate LibV1:trace_check #LibV1:peak_time LibV1:ISI_values #LibV1:peak_time LibV1:doublet_ISI #LibV1:peak_time -LibV1:peak_voltage #LibV1:peak_indices -LibV1:burst_ISI_indices #LibV1:peak_indices #LibV1:ISI_values +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 #LibV1:peak_indices +LibV1:peak_time #LibV5:peak_indices LibV1:time_to_first_spike #LibV1:peak_time LibV1:adaptation_index #LibV1:peak_time LibV1:adaptation_index2 #LibV1:peak_time LibV1:spike_width2 #LibV5:min_AHP_indices -LibV1:AP_width #LibV1:peak_indices #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:burst_number #LibV1:burst_mean_freq 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 #LibV1:peak_indices +LibV1:AHP_depth_abs_slow #LibV5:peak_indices LibV1:AHP_slow_time #LibV1:AHP_depth_abs_slow LibV1:time_constant LibV1:voltage_deflection @@ -27,17 +27,17 @@ LibV1:interpolate LibV1:steady_state_voltage LibV3:depolarized_base LibV1:ISI_CV #LibV1:ISI_values -LibV1:Spikecount #LibV4:peak_indices +LibV1:Spikecount #LibV5:peak_indices LibV1:AHP_depth #LibV5:voltage_base #LibV5:min_AHP_values -LibV2:AP_rise_indices #LibV1:peak_indices #LibV5:AP_begin_indices -LibV2:AP_end_indices #LibV1:peak_indices -LibV2:AP_fall_indices #LibV1:peak_indices #LibV5:AP_begin_indices #LibV2:AP_end_indices +LibV2:AP_rise_indices #LibV5:peak_indices #LibV5:AP_begin_indices +LibV2:AP_end_indices #LibV5:peak_indices +LibV2:AP_fall_indices #LibV5:peak_indices #LibV5:AP_begin_indices #LibV2:AP_end_indices LibV2:AP_duration #LibV5:AP_begin_indices #LibV2:AP_end_indices LibV2:AP_duration_half_width #LibV2:AP_rise_indices #LibV2:AP_fall_indices -LibV2:AP_rise_time #LibV5:AP_begin_indices #LibV1:peak_indices -LibV2:AP_fall_time #LibV1:peak_indices #LibV2:AP_end_indices -LibV2:AP_rise_rate #LibV5:AP_begin_indices #LibV1:peak_indices -LibV2:AP_fall_rate #LibV1:peak_indices #LibV2:AP_end_indices +LibV2:AP_rise_time #LibV5:AP_begin_indices #LibV5:peak_indices +LibV2:AP_fall_time #LibV5:peak_indices #LibV2:AP_end_indices +LibV2:AP_rise_rate #LibV5:AP_begin_indices #LibV5:peak_indices +LibV2:AP_fall_rate #LibV5:peak_indices #LibV2:AP_end_indices LibV2:fast_AHP #LibV5:AP_begin_indices #LibV5:min_AHP_indices LibV2:AP_amplitude_change #LibV1:AP_amplitude LibV2:AP_duration_change #LibV2:AP_duration @@ -86,7 +86,7 @@ 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 #LibV1:peak_indices +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 @@ -99,8 +99,8 @@ LibV5:AP2_AP1_diff #LibV1:AP_amplitude LibV5:AP2_AP1_peak_diff #LibV1:peak_voltage LibV5:AP1_width #LibV5:spike_half_width LibV5:AP2_width #LibV5:spike_half_width -LibV5:AHP_depth_from_peak #LibV1:peak_indices #LibV5:min_AHP_indices -LibV5:AHP_time_from_peak #LibV1:peak_indices #LibV5:min_AHP_indices +LibV5:AHP_depth_from_peak #LibV5:peak_indices #LibV5:min_AHP_indices +LibV5:AHP_time_from_peak #LibV5:peak_indices #LibV5:min_AHP_indices 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 @@ -112,7 +112,7 @@ 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 #LibV1:peak_indices +LibV5:spike_half_width #LibV5:min_AHP_indices #LibV5:peak_indices LibV5:AP_begin_indices #LibV5:min_AHP_indices LibV5:AHP_depth_abs #LibV5:min_AHP_values LibV5:AP_begin_width #LibV5:min_AHP_indices #LibV5:AP_begin_indices @@ -138,7 +138,7 @@ LibV5:BAC_width #LibV1:AP_width;location_epsp LibV5:BAC_maximum_voltage #LibV1:maximum_voltage;location_epsp LibV5:all_ISI_values #LibV1:peak_time LibV5:AP_amplitude_from_voltagebase #LibV5:voltage_base #LibV1:peak_voltage -LibV5:min_voltage_between_spikes #LibV1:peak_indices +LibV5:min_voltage_between_spikes #LibV5:peak_indices LibV5:voltage LibV5:steady_state_voltage_stimend LibV5:voltage_base diff --git a/tests/test_basic.py b/tests/test_basic.py index bcad538d..9751ec68 100644 --- a/tests/test_basic.py +++ b/tests/test_basic.py @@ -159,6 +159,17 @@ def test_setDependencyFileLocation_wrongpath(): efel.setDependencyFileLocation, "thisfiledoesntexist") +def test_setDependencyFileLocation(): + """basic: Test if setDependencyFileLocation works""" + import efel + efel.reset() + dep_file = os.path.join(os.path.dirname(os.path.abspath(__file__)), + 'DependencyV5_LibV5peakindices.txt') + efel.setDependencyFileLocation(dep_file) + result = efel.getDependencyFileLocation() + assert result == dep_file + + def test_nonexisting_feature(): """basic: Test nonexisting feature""" @@ -1435,42 +1446,6 @@ def test_spikecount_stimint1(): spikecount_stimint + 2) -def test_spikecount_libv4peakindices(): - """basic: Test Spikecount in combination with LibV4 peak_indices""" - - import efel - efel.reset() - - stim_start = 500.0 - stim_end = 900.0 - - time = efel.io.load_fragment('%s#col=1' % meanfrequency1_url) - voltage = efel.io.load_fragment('%s#col=2' % meanfrequency1_url) - - trace = {} - - trace['T'] = time - trace['V'] = voltage - trace['stim_start'] = [stim_start] - trace['stim_end'] = [stim_end] - - features = ['peak_indices', 'Spikecount'] - - test_peak = os.path.join(os.path.dirname(os.path.abspath(__file__)), - 'DependencyV5_LibV4peakindices.txt') - efel.setDependencyFileLocation(test_peak) - - feature_values = \ - efel.getFeatureValues( - [trace], - features) - - peak_indices = feature_values[0]['peak_indices'] - spikecount = feature_values[0]['Spikecount'][0] - assert len(peak_indices) == 5 - assert len(peak_indices) == spikecount - - def test_ohmic_inputresistance(): """basic: Test ohmic_input_resistance"""