Skip to content

Commit

Permalink
Merge pull request #353 from AurelienJaquier/fix-AP-end-indices
Browse files Browse the repository at this point in the history 
Fix AP_begin_indices-dependent features #minor
  • Loading branch information
@AurelienJaquier
AurelienJaquier authored Jan 11, 2024
2 parents 2bd34e6 + cb35767 commit 3ffb816
Show file tree
Hide file tree
Showing 8 changed files with 24,295 additions and 31 deletions.
7 changes: 7 additions & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,13 @@ All notable changes to this project will be documented in this file.
The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).

## [5.5.0] - 2024-01

### C++ changes
- AP_end_indices, AP_rise_time, AP_fall_time, AP_rise_rate, AP_fall_rate do not take into account peaks before stim_start anymore
- New test and test data for spontaneous firing case.
The data is provided by github user SzaBoglarka using cell https://modeldb.science/114047


## [5.4.0] - 2024-01

Expand Down
2 changes: 1 addition & 1 deletion efel/DependencyV5.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ LibV5:AP_end_indices #LibV5:peak_indices #LibV1:interpolate
LibV2:AP_fall_indices #LibV5:peak_indices #LibV5:AP_begin_indices #LibV5:AP_end_indices #LibV1:interpolate
LibV2:AP_duration #LibV5:AP_begin_indices #LibV5:AP_end_indices #LibV1:interpolate
LibV2:AP_duration_half_width #LibV2:AP_rise_indices #LibV2:AP_fall_indices #LibV1:interpolate
LibV2:AP_rise_time #LibV5:AP_begin_indices #LibV5:peak_indices #LibV1:AP_amplitude #LibV1:interpolate
LibV2:AP_rise_time #LibV5:AP_begin_indices #LibV5:peak_indices #LibV1:interpolate
LibV2:AP_fall_time #LibV5:peak_indices #LibV5:AP_end_indices #LibV1:interpolate
LibV2:AP_rise_rate #LibV5:AP_begin_indices #LibV5:peak_indices #LibV1:interpolate
LibV2:AP_fall_rate #LibV5:peak_indices #LibV5:AP_end_indices #LibV1:interpolate
Expand Down
57 changes: 40 additions & 17 deletions efel/cppcore/LibV2.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -168,28 +168,38 @@ int LibV2::AP_duration_half_width(mapStr2intVec& IntFeatureData,
static int __AP_rise_time(const vector<double>& t, const vector<double>& v,
const vector<int>& apbeginindices,
const vector<int>& peakindices,
const vector<double>& apamplitude, double beginperc,
double endperc, vector<double>& aprisetime) {
double beginperc, double endperc,
vector<double>& aprisetime) {
aprisetime.resize(std::min(apbeginindices.size(), peakindices.size()));
// Make sure that we do not use peaks starting before the 1st AP_begin_index
// Because AP_begin_indices only takes into account peaks after stimstart
vector<int> newpeakindices;
if (apbeginindices.size() > 0) {
newpeakindices = peaks_after_stim_start(apbeginindices[0], peakindices);
}
double begin_v;
double end_v;
double begin_indice;
double end_indice;
double apamplitude;
for (size_t i = 0; i < aprisetime.size(); i++) {
begin_v = v[apbeginindices[i]] + beginperc * apamplitude[i];
end_v = v[apbeginindices[i]] + endperc * apamplitude[i];
// do not use AP_amplitude feature because it does not take into account
// peaks after stim_end
apamplitude = v[newpeakindices[i]] - v[apbeginindices[i]];
begin_v = v[apbeginindices[i]] + beginperc * apamplitude;
end_v = v[apbeginindices[i]] + endperc * apamplitude;

// Get begin indice
size_t j = apbeginindices[i];
// change slightly begin_v for almost equal case
// truncature error can change begin_v even when beginperc == 0.0
while (j < peakindices[i] && v[j] < begin_v - 0.0000000000001) {
while (j < newpeakindices[i] && v[j] < begin_v - 0.0000000000001){
j++;
}
begin_indice = j;

// Get end indice
j = peakindices[i];
j = newpeakindices[i];
// change slightly end_v for almost equal case
// truncature error can change end_v even when beginperc == 0.0
while (j > apbeginindices[i] && v[j] > end_v + 0.0000000000001) {
Expand All @@ -207,14 +217,13 @@ int LibV2::AP_rise_time(mapStr2intVec& IntFeatureData,
// Fetching all required features in one go.
const auto& doubleFeatures = getFeatures(
DoubleFeatureData,
{"T", "V", "AP_amplitude", "rise_start_perc", "rise_end_perc"});
{"T", "V", "rise_start_perc", "rise_end_perc"});
const auto& intFeatures =
getFeatures(IntFeatureData, {"AP_begin_indices", "peak_indices"});
vector<double> aprisetime;
int retval = __AP_rise_time(
doubleFeatures.at("T"), doubleFeatures.at("V"),
intFeatures.at("AP_begin_indices"), intFeatures.at("peak_indices"),
doubleFeatures.at("AP_amplitude"),
doubleFeatures.at("rise_start_perc").empty()
? 0.0
: doubleFeatures.at("rise_start_perc").front(),
Expand All @@ -230,28 +239,34 @@ int LibV2::AP_rise_time(mapStr2intVec& IntFeatureData,

// *** AP_fall_time according to E10 and E18 ***
static int __AP_fall_time(const vector<double>& t,
const double stimstart,
const vector<int>& peakindices,
const vector<int>& apendindices,
vector<double>& apfalltime) {
apfalltime.resize(std::min(peakindices.size(), apendindices.size()));
// Make sure that we do not use peaks starting before stim start
// Because AP_end_indices only takes into account peaks after stim start
vector<int> newpeakindices = peaks_after_stim_start(stimstart, peakindices, t);

for (size_t i = 0; i < apfalltime.size(); i++) {
apfalltime[i] = t[apendindices[i]] - t[peakindices[i]];
apfalltime[i] = t[apendindices[i]] - t[newpeakindices[i]];
}
return apfalltime.size();
}
int LibV2::AP_fall_time(mapStr2intVec& IntFeatureData,
mapStr2doubleVec& DoubleFeatureData,
mapStr2Str& StringData) {
const auto& doubleFeatures = getFeatures(DoubleFeatureData, {"T"});
const auto& doubleFeatures = getFeatures(DoubleFeatureData, {"T", "stim_start"});
const auto& intFeatures =
getFeatures(IntFeatureData, {"peak_indices", "AP_end_indices"});

const vector<double>& t = doubleFeatures.at("T");
const double stim_start = doubleFeatures.at("stim_start")[0];
const vector<int>& peakindices = intFeatures.at("peak_indices");
const vector<int>& apendindices = intFeatures.at("AP_end_indices");

vector<double> apfalltime;
int retval = __AP_fall_time(t, peakindices, apendindices, apfalltime);
int retval = __AP_fall_time(t, stim_start, peakindices, apendindices, apfalltime);

if (retval > 0) {
setVec(DoubleFeatureData, StringData, "AP_fall_time", apfalltime);
Expand All @@ -265,9 +280,13 @@ static int __AP_rise_rate(const vector<double>& t, const vector<double>& v,
const vector<int>& peakindices,
vector<double>& apriserate) {
apriserate.resize(std::min(peakindices.size(), apbeginindices.size()));
vector<int> newpeakindices;
if (apbeginindices.size() > 0) {
newpeakindices = peaks_after_stim_start(apbeginindices[0], peakindices);
}
for (size_t i = 0; i < apriserate.size(); i++) {
apriserate[i] = (v[peakindices[i]] - v[apbeginindices[i]]) /
(t[peakindices[i]] - t[apbeginindices[i]]);
apriserate[i] = (v[newpeakindices[i]] - v[apbeginindices[i]]) /
(t[newpeakindices[i]] - t[apbeginindices[i]]);
}
return apriserate.size();
}
Expand All @@ -294,30 +313,34 @@ int LibV2::AP_rise_rate(mapStr2intVec& IntFeatureData,

// *** AP_fall_rate according to E12 and E20 ***
static int __AP_fall_rate(const vector<double>& t, const vector<double>& v,
const double stimstart,
const vector<int>& peakindices,
const vector<int>& apendindices,
vector<double>& apfallrate) {
apfallrate.resize(std::min(apendindices.size(), peakindices.size()));
vector<int> newpeakindices = peaks_after_stim_start(stimstart, peakindices, t);

for (size_t i = 0; i < apfallrate.size(); i++) {
apfallrate[i] = (v[apendindices[i]] - v[peakindices[i]]) /
(t[apendindices[i]] - t[peakindices[i]]);
apfallrate[i] = (v[apendindices[i]] - v[newpeakindices[i]]) /
(t[apendindices[i]] - t[newpeakindices[i]]);
}
return apfallrate.size();
}
int LibV2::AP_fall_rate(mapStr2intVec& IntFeatureData,
mapStr2doubleVec& DoubleFeatureData,
mapStr2Str& StringData) {
const auto& doubleFeatures = getFeatures(DoubleFeatureData, {"T", "V"});
const auto& doubleFeatures = getFeatures(DoubleFeatureData, {"T", "V", "stim_start"});
const auto& intFeatures =
getFeatures(IntFeatureData, {"peak_indices", "AP_end_indices"});

const vector<double>& t = doubleFeatures.at("T");
const vector<double>& v = doubleFeatures.at("V");
const double stim_start = doubleFeatures.at("stim_start")[0];
const vector<int>& peakindices = intFeatures.at("peak_indices");
const vector<int>& apendindices = intFeatures.at("AP_end_indices");

vector<double> apfallrate;
int retval = __AP_fall_rate(t, v, peakindices, apendindices, apfallrate);
int retval = __AP_fall_rate(t, v, stim_start, peakindices, apendindices, apfallrate);

if (retval > 0) {
setVec(DoubleFeatureData, StringData, "AP_fall_rate", apfallrate);
Expand Down
34 changes: 22 additions & 12 deletions efel/cppcore/LibV5.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -536,31 +536,40 @@ int LibV5::AP_begin_indices(mapStr2intVec& IntFeatureData,
}

static int __AP_end_indices(const vector<double>& t, const vector<double>& v,
const vector<int>& pi, vector<int>& apei,
double derivativethreshold) {
const double stimstart, const vector<int>& pi,
vector<int>& apei, double derivativethreshold) {

vector<double> dvdt(v.size());
vector<double> dv;
vector<double> dt;
vector<int> peak_indices;
int max_slope;
getCentralDifferenceDerivative(1., v, dv);
getCentralDifferenceDerivative(1., t, dt);
transform(dv.begin(), dv.end(), dt.begin(), dvdt.begin(),
std::divides<double>());

apei.resize(pi.size());
vector<int> picopy(pi.begin(), pi.end());
picopy.push_back(v.size() - 1);
auto stimbeginindex = distance(t.begin(),
find_if(t.begin(), t.end(),
[stimstart](double time){ return time >= stimstart; }));

for (size_t i = 0; i < pi.size(); i++) {
if (pi[i] > stimbeginindex) {
peak_indices.push_back(pi[i]);
}
}
peak_indices.push_back(v.size() - 1);

for (size_t i = 0; i < apei.size(); i++) {
for (size_t i = 0; i < peak_indices.size() - 1; i++) {
max_slope =
distance(dvdt.begin(), std::min_element(dvdt.begin() + picopy[i] + 1,
dvdt.begin() + picopy[i + 1]));
distance(dvdt.begin(), std::min_element(dvdt.begin() + peak_indices[i] + 1,
dvdt.begin() + peak_indices[i + 1]));
// assure that the width of the slope is bigger than 4
apei[i] = distance(dvdt.begin(), find_if(dvdt.begin() + max_slope,
dvdt.begin() + picopy[i + 1],
apei.push_back(distance(dvdt.begin(), find_if(dvdt.begin() + max_slope,
dvdt.begin() + peak_indices[i + 1],
[derivativethreshold](double x) {
return x >= derivativethreshold;
}));
})));
}
return apei.size();
}
Expand All @@ -570,14 +579,15 @@ int LibV5::AP_end_indices(mapStr2intVec& IntFeatureData,
mapStr2Str& StringData) {
const auto& T = getFeature(DoubleFeatureData, "T");
const auto& V = getFeature(DoubleFeatureData, "V");
const auto& stim_start = getFeature(DoubleFeatureData, "stim_start");
const auto& peak_indices = getFeature(IntFeatureData, "peak_indices");

vector<double> dTh;
int retVal = getParam(DoubleFeatureData, "DownDerivativeThreshold", dTh);
double downDerivativeThreshold = (retVal <= 0) ? -12.0 : dTh[0];

vector<int> AP_end_indices;
retVal = __AP_end_indices(T, V, peak_indices, AP_end_indices,
retVal = __AP_end_indices(T, V, stim_start[0], peak_indices, AP_end_indices,
downDerivativeThreshold);
if (retVal >= 0) {
setVec(IntFeatureData, StringData, "AP_end_indices", AP_end_indices);
Expand Down
25 changes: 25 additions & 0 deletions efel/cppcore/Utils.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -200,6 +200,31 @@ std::pair<size_t, size_t> get_time_index(const std::vector<double>& t,
return std::pair<size_t, size_t>(startIndex, endIndex);
}

// expects an empty newpeakindices
vector<int> peaks_after_stim_start(const double stimstart,
const vector<int>& peakindices,
const vector<double>& t) {
vector<int> newpeakindices;
for (size_t i = 0; i < peakindices.size(); i++) {
if (t[peakindices[i]] > stimstart) {
newpeakindices.push_back(peakindices[i]);
}
}
return newpeakindices;
}

// expects an empty newpeakindices
vector<int> peaks_after_stim_start(const double stimstartindex,
const vector<int>& peakindices){
vector<int> newpeakindices;
for (size_t i = 0; i < peakindices.size(); i++) {
if (peakindices[i] > stimstartindex) {
newpeakindices.push_back(peakindices[i]);
}
}
return newpeakindices;
}

template <class T>
double vec_mean(const vector<T>& v) {
/*
Expand Down
9 changes: 8 additions & 1 deletion efel/cppcore/Utils.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,14 @@ double vec_mean(const vector<T> &v);

std::pair<size_t, size_t> get_time_index(const std::vector<double> &t, double startTime,
double endTime, double precisionThreshold);


vector<int> peaks_after_stim_start(const double stimstart,
const vector<int>& peakindices,
const vector<double>& t);

vector<int> peaks_after_stim_start(const double stimstartindex,
const vector<int>& peakindices);

template <class ForwardIterator>
ForwardIterator first_min_element(ForwardIterator first, ForwardIterator last) {
ForwardIterator lowest = first;
Expand Down
Loading

0 comments on commit 3ffb816

Please sign in to comment.