Commit a38284dc authored by Damien Leroux's avatar Damien Leroux
Browse files

Multiple bugfixes. Now handling Euralis data properly.

parent 1507a980
......@@ -9,13 +9,13 @@ depend:
cd src/pedigree/ && make .depend
spell-qtl: depend
cd src && $(MAKE)
+cd src && $(MAKE)
spell-pedigree: depend
cd src/pedigree && $(MAKE)
+cd src/pedigree && $(MAKE)
spell-marker: depend
cd src/bayes && $(MAKE)
+cd src/bayes && $(MAKE)
clean:
cd src && make clean
......
......@@ -21,11 +21,11 @@ VERSION_DEFINES=$(shell git tag | tail -1 | sed 's/\([0-9]*\)[.]\([0-9]*\)\([.][
CXXARGS=-Wextra -Wall -pthread -fPIC -DEIGEN_NO_DEPRECATED_WARNING $(VERSION_DEFINES)
#CXXARGS=-Wextra -Wall -pthread -fPIC -DEIGEN_NO_DEPRECATED_WARNING $(VERSION_DEFINES)
DEBUG_OPTS=-ggdb -O0
#DEBUG_OPTS=-ggdb -O0
#DEBUG_OPTS=-ggdb -O0 -DNDEBUG
#DEBUG_OPTS=-ggdb -O2 -DNDEBUG
#OPT_OPTS=-O2 -DEIGEN_NO_DEBUG -DNDEBUG
#OPT_OPTS=-O3 -DEIGEN_NO_DEBUG -DNDEBUG
OPT_OPTS=-O3 -DEIGEN_NO_DEBUG -DNDEBUG
#OPT_OPTS=-ggdb -O2 -DEIGEN_NO_DEBUG -DNDEBUG
#OPT_OPTS=-ggdb -O2 -DNDEBUG
#DEBUG_OPTS=-ggdb -DNDEBUG
......
......@@ -655,17 +655,44 @@ struct LV_database {
const auto& mark_ind_lv = gen_dat.second;
size_t n_states = mark_ind_lv.begin()->second.front().size();
size_t n_mark = chr.raw.marker_name.size();
size_t n_condensed_mark = chr.condensed.marker_name.size();
size_t n_ind = mark_ind_lv.begin()->second.size();
auto& ind_LVmat = data[chr.name][gen];
ind_LVmat.resize(n_ind, {n_states, n_mark});
/*MSG_DEBUG("[assemble.using.map] #haplo_size=" << chr.haplo_sizes.size());*/
/*MSG_QUEUE_FLUSH();*/
for (size_t ind = 0; ind < n_ind; ++ind) {
auto& mat = ind_LVmat[ind];
mat.resize(n_states, n_mark);
mat = MatrixXd::Ones(n_states, n_condensed_mark);
size_t m = 0;
int haplo_counter = 0;
auto haplo_size_i = chr.haplo_sizes.begin();
for (const auto& mark: chr.raw.marker_name) {
const auto& svec = mark_ind_lv.find(mark)->second;
mat.col(m) = svec[ind];
++m;
/*MSG_DEBUG("[assemble.using.map] haplo_size=" << (*haplo_size_i) << " column=" << m);*/
/*MSG_QUEUE_FLUSH();*/
mat.col(m).array() *= svec[ind].array();
++haplo_counter;
if (haplo_counter == *haplo_size_i) {
++m;
++haplo_size_i;
haplo_counter = 0;
}
}
auto hi = chr.begin();
for (int c = 0; c < mat.cols(); ++c, ++hi) {
double sum = mat.col(c).sum();
if (sum == 0) {
std::stringstream ss;
ss << "Detected inconsistent observations! Generation is " << gen << ", chromosome " << chr.name << ", individual #" << (ind + 1) << ", involving markers";
auto haplo = *hi;
for (size_t i = haplo.first; i < haplo.second; ++i) {
ss << ' ' << chr.raw.marker_name[i];
}
MSG_ERROR(ss.str(), "Fix the genotype data");
} else {
mat.col(c) /= sum;
}
}
}
}
......@@ -853,14 +880,20 @@ struct qtl_pop_type {
std::map<std::string, std::vector<MatrixXd>> LV;
std::string observed_traits_filename;
std::vector<trait> observed_traits;
std::map<char, std::string> ancestor_names;
qtl_pop_type()
: name(), qtl_generation_name(), indices(), gen(), LV(), observed_traits_filename(), observed_traits()
{}
qtl_pop_type(const std::string& n, const std::string& qgn, const std::vector<size_t>& i, std::shared_ptr<geno_matrix> g,
const std::map<std::string, std::vector<MatrixXd>>& lv, const std::string& otf, const std::vector<trait>& ot)
: name(n), qtl_generation_name(qgn), indices(i), gen(g), LV(lv), observed_traits_filename(otf), observed_traits(ot)
{}
const std::map<std::string, std::vector<MatrixXd>>& lv, const std::string& otf, const std::vector<trait>& ot, const std::map<char, std::string>& anam)
: name(n), qtl_generation_name(qgn), indices(i), gen(g), LV(lv), observed_traits_filename(otf), observed_traits(ot), ancestor_names(anam)
{
/*MSG_DEBUG("new qtl_pop with " << ancestor_names.size() << " ancestor names");*/
/*for (const auto& kv: ancestor_names) {*/
/*MSG_DEBUG(" * " << kv.first << ": " << kv.second);*/
/*}*/
}
std::shared_ptr<qtl_pop_type>
filter_clone(const trait& this_trait) const
......@@ -872,6 +905,7 @@ struct qtl_pop_type {
ret->gen = gen;
ret->observed_traits_filename = observed_traits_filename;
ret->observed_traits.emplace_back(this_trait);
ret->ancestor_names = ancestor_names;
auto i = indices.begin();
auto j = indices.end();
auto filt = keep.begin();
......@@ -923,6 +957,7 @@ struct pop_data_type {
LV_database LV;
std::map<std::string, std::vector<size_t>> families;
std::map<size_t, std::shared_ptr<geno_matrix>> gen_by_id;
std::map<char, std::string> ancestor_names;
void
assemble_chromosomes(const std::vector<chromosome>& map)
......@@ -991,6 +1026,12 @@ struct pop_data_type {
write_size(ofs, kv.first);
write_size(ofs, std::find(generations.begin(), generations.end(), kv.second) - generations.begin());
}
write_fourcc(ofs, "ANAM");
write_size(ofs, ancestor_names.size());
for (const auto& kv: ancestor_names) {
write_char(ofs, kv.first);
write_str(ofs, kv.second);
}
}
static
......@@ -1044,6 +1085,12 @@ struct pop_data_type {
size_t k = read_size(ifs);
ret.gen_by_id[k] = ret.generations[read_size(ifs)];
}
CHECK_4CC("ANAM");
size_t n_an = read_size(ifs);
for (size_t a = 0; a < n_an; ++a) {
char c = read_char(ifs);
ret.ancestor_names[c] = read_str(ifs);
}
return ret;
}
......@@ -1110,7 +1157,7 @@ struct pop_data_type {
}
auto aqg = all_qtl_geno_matrix();
for (const auto& kv: aqg) {
std::shared_ptr<qtl_pop_type> this_pop = std::make_shared<qtl_pop_type>(name, qtl_generation_name, kv.second, kv.first, LV.extract(qtl_generation_name, kv.second), traits_filename, extract_traits(kv.second));
std::shared_ptr<qtl_pop_type> this_pop = std::make_shared<qtl_pop_type>(name, qtl_generation_name, kv.second, kv.first, LV.extract(qtl_generation_name, kv.second), traits_filename, extract_traits(kv.second), ancestor_names);
for (const auto& this_trait: traits) {
auto& lp = linked_pops[this_trait.name];
auto pop_ptr = this_pop->filter_clone(this_trait);
......
......@@ -6,6 +6,15 @@
#include "model.h"
/*#include "model/tests.h"*/
#include <boost/math/distributions/normal.hpp> // for normal_distribution
using boost::math::normal; // typedef provides default type is double.
using boost::math::cdf;
using boost::math::mean;
using boost::math::variance;
using boost::math::quantile;
using boost::math::complement;
std::pair<bool, double>
detect_strongest_qtl(chromosome_value chr, const locus_key& lk,
const model& M0, const std::vector<double> pos);
......@@ -172,7 +181,8 @@ struct QTL {
}
std::pair<double, double>
confidence_interval() const
confidence_interval(const std::string& trait, const std::vector<QTL>& selection);
#if 0
{
/*MSG_DEBUG_INDENT_EXPR("[Confidence interval] ");*/
/*MSG_DEBUG("LOD: " << LOD);*/
......@@ -200,6 +210,7 @@ struct QTL {
/*MSG_DEBUG_DEDENT;*/
return {left, right};
}
#endif
};
......@@ -247,18 +258,30 @@ struct analysis_report {
report_file.close();
report_file.open(MESSAGE(report_path << "/full_map.txt"), std::fstream::out);
for (const auto& chr: active_settings->map) {
report_file << chr.pretty_print(80, poi[chr.name], roi[chr.name]) << std::endl;
report_file << chr.pretty_print(200, poi[chr.name], roi[chr.name]) << std::endl;
}
}
void attach_model_manager(model_manager& mm);
void detach_model_manager(model_manager& mm);
void report_trait(const std::string& name, const MatrixXd& values)
void report_trait(const std::string& /*name*/, const MatrixXd& values)
{
static Eigen::IOFormat trait_format(Eigen::FullPrecision, Eigen::DontAlignCols, "\t", "\n", "", "", "", "");
std::string filename = MESSAGE(full_path << '/' << "trait_values.txt");
ofile(filename) << values.format(trait_format);
ofile of(filename);
of << values.format(trait_format);
of.close();
}
void report_lod(const QTL& qtl)
{
std::string filename = MESSAGE(full_path << '/' << qtl.chromosome << ':' << qtl.locus << "_LOD.txt");
ofile of(filename);
for (size_t i = 0; i < qtl.LOD.size(); ++i) {
of << qtl.LOD_loci[i] << '\t' << qtl.LOD[i] << std::endl;
}
of.close();
}
void report_model(const model& Mcurrent)
......@@ -282,7 +305,7 @@ struct analysis_report {
return ret.str();
}
void report_computation(const model& Mcurrent, const chromosome* chrom_under_study, const computation_along_chromosome& cac, ComputationType ct, ComputationResults cr, const std::vector<double>& testpos, probability_mode pmode)
void report_computation(const model& Mcurrent, const chromosome* chrom_under_study, const computation_along_chromosome& cac, ComputationType ct, ComputationResults /*cr*/, const std::vector<double>& testpos, probability_mode pmode)
{
if (output_test | output_rss | output_rank) {
/*MSG_DEBUG(MATRIX_SIZE(cac.ftest_pvalue));*/
......@@ -330,12 +353,13 @@ struct analysis_report {
}
f << std::endl;
}
f.close();
}
}
void report_final_model(model_manager& mm);
void report_qtls(const std::vector<QTL>& qtls);
void report_qtls(std::vector<QTL>& qtls);
};
......@@ -376,7 +400,7 @@ struct model_manager {
/* slight hack. */
std::vector<double> testpos;
/* additional data for output and postprocessing */
/* additiona, colpops.front()->ancestor_namesl data for output and postprocessing */
std::map<std::string, std::map<double, std::pair<double, double>>> qtl_confidence_interval;
analysis_report* reporter;
......@@ -391,7 +415,7 @@ struct model_manager {
, all_pops(colpops)
, pop_blocks()
, n_par()
, Mcurrent(y, colpops, st)
, Mcurrent(y, colpops, (*colpops.front())->ancestor_names, st)
, threshold(th)
, test_type(ct)
, all_loci()
......@@ -1182,32 +1206,246 @@ void analysis_report::detach_model_manager(model_manager& mm)
mm.reporter = NULL;
}
inline
const std::string&
significance_string(double pvalue)
{
static std::string
threestar = "***",
twostar = "**",
onestar = "*",
zerostar = ".",
bad = " ";
if (pvalue <= .001) {
return threestar;
} else if (pvalue <= .01) {
return twostar;
} else if (pvalue <= .05) {
return onestar;
} else if (pvalue <= .1) {
return zerostar;
} else {
return bad;
}
}
/* FIXME: define some kind of global array so significance strings and legend are always consistent. */
inline
const std::string&
significance_legend() { static std::string _ = "Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1"; return _; }
inline
std::string
format_pvalue(double pvalue)
{
std::stringstream pv;
double d = 100 * pvalue;
pv << std::setw(7);
if (d < .1) {
pv.setf(std::ios::scientific, std::ios::floatfield);
pv.precision(1);
} else {
pv.setf(std::ios::fixed, std::ios::floatfield);
pv.precision(2);
}
pv << d << "% " << std::setw(3) << std::left << significance_string(pvalue);
return pv.str();
}
inline std::string empty_cell() { return "____________"; }
inline
std::pair<std::vector<size_t>, std::vector<std::pair<model_block_key, std::vector<std::vector<char>>>>>
filter_columns(const MatrixXd& conrow, const model_block_collection& blocks)
{
std::pair<std::vector<size_t>, std::vector<std::pair<model_block_key, std::vector<std::vector<char>>>>> ret;
auto& columns = ret.first;
auto& block_parts = ret.second;
size_t j = 0;
for (const auto& kv: blocks) {
const auto& labels = kv.second->column_labels;
for (size_t i = 0; i < labels.size(); ++i, ++j) {
if (conrow(j)) {
if (!block_parts.size() || block_parts.back().first != kv.first) {
MSG_DEBUG("first used column in block " << kv.first);
block_parts.emplace_back();
block_parts.back().first = kv.first;
block_parts.back().second.push_back(labels[i]);
} else {
block_parts.back().second.push_back(labels[i]);
}
MSG_DEBUG("added column #" << j);
MSG_QUEUE_FLUSH();
columns.push_back(j);
}
}
}
return ret;
}
inline
file& header(file& f, const std::string& str)
{
size_t w = str.size() + 1;
auto fill = [&] (char c) -> file& { char prev = f.m_impl.fill(c); f << std::setw(w) << ""; f.m_impl.fill(prev); return f; };
f << std::endl;
fill(' ') << '|' << std::endl
<< str << " |" << std::endl;
fill('_') << '|' << std::endl << std::endl;
return f;
}
inline
void analysis_report::report_final_model(model_manager& mm)
{
mm.Mcurrent.output_X_to_file(full_path);
mm.Mcurrent.output_XtX_inv_to_file(full_path);
model::xtx_printable xtx(mm.Mcurrent, true);
const auto& rss = mm.Mcurrent.rss();
/*const auto& XtX_inv = mm.Mcurrent.XtX_pseudo_inverse();*/
const auto& X = mm.Mcurrent.X();
/*const auto& XtX = mm.Mcurrent.XtX();*/
/*size_t dof = mm.Mcurrent.rank();*/
/*size_t n = mm.Mcurrent.Y().rows();*/
auto qtls = mm.QTLs();
report_qtls(qtls);
/*
report_file
<< "Final model" << std::endl << mm.Mcurrent << std::endl
<< "XtX^-1" << std::endl << xtx << std::endl
<< "Coefficients " << mm.Mcurrent.coefficients().transpose() << std::endl
<< "RSS " << mm.Mcurrent.rss() << std::endl
<< std::endl;
report_qtls(qtls);
*/
if (mm.Mcurrent.n_obs() <= mm.Mcurrent.dof() - 1) {
report_file << "* There are not enough observations to compute a variance/covariance matrix." << std::endl;
} else {
static std::string coefstr;
for (int i = 0; i < rss.size(); ++i) {
report_file
<< "##########################" << std::endl
<< "## Variable #" << (i + 1) << std::endl;
MatrixXd coef = mm.Mcurrent.coefficients().col(i).transpose();
MSG_DEBUG(MATRIX_SIZE(coef));
{
model_print::matrix_with_sections<std::string, void, model_block_key, std::vector<char>> mws(coef);
for (const auto& kv: mm.Mcurrent.m_blocks) {
mws.add_column_section(kv.first, mm.Mcurrent.labels_to_ancestor_names(*kv.second));
}
mws.add_row_section(coefstr, 1);
header(report_file, "Coefficients") << mws << std::endl;
}
/*report_file*/
/*<< "Coefficients " << std::endl << mws << std::endl;*/
/*double var_resid_est = rss(i) / (n - dof - 1);*/
MatrixXd vcov = mm.Mcurrent.vcov(i);
model::xtx_printable vc(mm.Mcurrent, vcov);
header(report_file, "Covariance matrix") << vc << std::endl;
auto conmat = X.bottomRows(X.rows() - mm.Mcurrent.n_obs());
for (int r = 0; r < conmat.rows(); ++r) {
auto subset = filter_columns(conmat.row(r), mm.Mcurrent.m_blocks);
const auto& columns = subset.first;
const auto& blocks = subset.second;
/*Eigen::Matrix<std::string, Eigen::Dynamic, Eigen::Dynamic> pvmat = Eigen::Matrix<std::string, Eigen::Dynamic, Eigen::Dynamic>::Constant(columns.size(), columns.size(), empty_cell());*/
Eigen::Matrix<std::string, Eigen::Dynamic, Eigen::Dynamic> pvmat(columns.size(), columns.size());
Eigen::Matrix<std::string, Eigen::Dynamic, Eigen::Dynamic> ctrmat(columns.size(), columns.size());
/*MatrixXd ctrmat = MatrixXd::Constant(columns.size(), columns.size(), std::numeric_limits<double>::quiet_NaN());*/
for (int i1 = 0; i1 < pvmat.cols(); ++i1) {
for (int i2 = i1 + 1; i2 < pvmat.cols(); ++i2) {
int a1 = columns[i1];
int a2 = columns[i2];
double sigma = vcov(a1, a1) + vcov(a2, a2) - 2 * vcov(a1, a2);
double delta = coef(a1) - coef(a2);
normal s(0, sqrt(sigma));
double pvalue = 2 * cdf(complement(s, fabs(delta)));
ctrmat(i1, i2) = MESSAGE(std::setprecision(3) << -delta << ' ' << std::setw(3) << std::left << significance_string(pvalue));
ctrmat(i2, i1) = MESSAGE(std::setprecision(3) << delta << ' ' << std::setw(3) << std::left << significance_string(pvalue));
pvmat(i1, i2) = pvmat(i2, i1) = format_pvalue(pvalue);
}
}
model_print::matrix_with_sections<void, std::vector<char>, model_block_key, std::vector<char>, Eigen::Matrix<std::string, Eigen::Dynamic, Eigen::Dynamic>> ctrsig(pvmat);
model_print::matrix_with_sections<void, std::vector<char>, model_block_key, std::vector<char>, Eigen::Matrix<std::string, Eigen::Dynamic, Eigen::Dynamic>> ctr(ctrmat);
for (const auto& kl: blocks) {
auto names = mm.Mcurrent.labels_to_ancestor_names(kl.second);
for (auto& n: names) { n.push_back(' '); n.push_back(' '); n.push_back(' '); n.push_back(' '); }
ctr.add_row_section(names);
ctr.add_column_section(kl.first, names);
ctrsig.add_row_section(names);
ctrsig.add_column_section(kl.first, names);
}
header(report_file, "Contrasts") << ctr << std::endl << significance_legend() << std::endl;
header(report_file, "Contrast significance") << ctrsig << std::endl << significance_legend() << std::endl;
}
#if 0
std::vector<std::vector<int>> groups(conmat.rows());
for (int i = 0; i < conmat.rows(); ++i) {
for (int j = 0; j < conmat.cols(); ++j) {
if (conmat(i, j)) {
groups[i].push_back(j);
}
}
}
std::vector<std::string> colnames;
colnames.reserve(X.cols());
for (const auto& k_blk: mm.Mcurrent.m_blocks) {
for (const auto& v: mm.Mcurrent.labels_to_ancestor_names(*k_blk.second)) {
colnames.emplace_back(v.begin(), v.end());
}
}
for (const auto& grp: groups) {
report_file << "Contrast group";
for (int j: grp) {
report_file << ' ' << mm.Mcurrent.m_ancestor_names['a' + j];
}
report_file << std::endl;
for (size_t i1 = 0; i1 < grp.size(); ++i1) {
for (size_t i2 = i1 + 1; i2 < grp.size(); ++i2) {
int a1 = grp[i1];
int a2 = grp[i2];
std::string name1 = colnames[a1];
std::string name2 = colnames[a2];
double sigma = vcov(a1, a1) + vcov(a2, a2) - 2 * vcov(a1, a2);
double delta = coef(a1) - coef(a2);
normal s(0, sqrt(sigma));
/*MSG_DEBUG("Normal distribution, mean = " << s.mean() << ", variance = " << variance(s));*/
double q = quantile(s, .975);
bool out = fabs(delta) > q;
double pvalue = 2 * cdf(complement(s, fabs(delta)));
/*MSG_DEBUG("quantile(cdf(.975)) == x ? " << (fabs(quantile(s, cdf(s, .975)) - .975) < 1.e-5));*/
/*report_file << "** " << name1 << " -- " << name2 << " µ=" << mean(s) << " sigma=" << variance(s) << std::endl;*/
/*report_file << name1 << " -- " << name2 << ": " << (out ? "true" : "false") << " quantile(.975)=" << q << " delta=" << delta << " cdf(" << fabs(delta) << ")=" << cdf(s, fabs(delta)) << std::endl;*/
/*report_file << name1 << " -- " << name2 << ": " << pvalue << " (-log10 = " << (-log10(pvalue)) << ") " << (out ? "true" : "false") << std::endl;*/
std::string pvstr = format_pvalue(pvalue);
report_file << name1 << " -- " << name2 << ": " << pvstr << std::endl;
}
}
}
#endif
}
}
}
inline
void analysis_report::report_qtls(const std::vector<QTL>& qtls)
void analysis_report::report_qtls(std::vector<QTL>& qtls)
{
for (const auto& qtl: qtls) {
for (auto& qtl: qtls) {
const auto& ci = roi[qtl.chromosome][qtl.locus] = qtl.confidence_interval(trait_name, qtls);
report_lod(qtl);
std::string name = MESSAGE(trait_name << " @ " << qtl.locus << " [" << ci.first << ':' << ci.second << ']');
if (poi[qtl.chromosome][qtl.locus].size()) {
poi[qtl.chromosome][qtl.locus] = MESSAGE(poi[qtl.chromosome][qtl.locus] << ',' << trait_name);
poi[qtl.chromosome][qtl.locus] = MESSAGE(poi[qtl.chromosome][qtl.locus] << ',' << name);
} else {
poi[qtl.chromosome][qtl.locus] = trait_name;
poi[qtl.chromosome][qtl.locus] = name;
}
roi[qtl.chromosome][qtl.locus] = qtl.confidence_interval();
report_file << "QTL detected on chromosome " << qtl.chromosome << " at " << qtl.locus << "cM with confidence interval {" << roi[qtl.chromosome][qtl.locus].first << ':' << roi[qtl.chromosome][qtl.locus].second << '}' << std::endl;
}
}
......@@ -1569,7 +1807,7 @@ struct test_all_domain : public test_domain_search_method {
test_all_domain() : sd_begin(), sd_end(), sd_cur() {}
void
init(const test_manager* tm, const std::vector<genome_search_domain>& dom, locus_set& point)
init(const test_manager*, const std::vector<genome_search_domain>& dom, locus_set& point)
{
sd_begin.clear();
sd_cur.clear();
......@@ -1598,7 +1836,7 @@ struct test_all_domain : public test_domain_search_method {
point.reserve(sd_cur.size());
bool carry;
do {
MSG_DEBUG('#' << counter);
/*MSG_DEBUG('#' << counter);*/
carry = false;
++sd_cur[i];
if (sd_cur[i] == sd_end[i]) {
......@@ -1618,5 +1856,59 @@ struct test_all_domain : public test_domain_search_method {
}
};
inline
std::pair<double, double>
QTL::confidence_interval(const std::string& trait, const std::vector<QTL>& selection)
{
test_all_domain tad;
test_manager tm(trait, 1, 0, true);
tm.test_type = ComputationType::FTestLOD;
std::vector<genome_search_domain> vgsd(1);
const ::chromosome* chr = active_settings->find_chromosome(chromosome);
auto steps = compute_steps(chr->condensed.marker_locus, active_settings->step);
vgsd.back().emplace_back(chr, steps);
tm.compute(vgsd, tad);
size_t first = 0, last = steps.size() - 1;
double step_min, step_max;
auto lod = [&] (double l) { locus_set point = {{chr, l}}; return tm.results[point]; };
double threshold = lod(locus) - 1.5;
{
LOD.clear();
LOD.reserve(tm.results.size());
LOD_loci = steps;
for (const auto& r: tm.results) {
LOD.push_back(r.second);
}
}
{
auto i = tm.results.begin();
while (i->second < threshold) { ++i; ++first; }
if (first == 0) {
step_min = 0;
} else {
step_min = interpolate(steps[first - 1], lod(steps[first - 1]), steps[first], lod(steps[first]), threshold);
}
}
{
/* TODO add selection \ {qtl} */
(void) selection;
auto i = tm.results.rbegin();
while (i->second < threshold) { ++i; --last; }
if (last == steps.size() - 1) {
step_max = steps.back();
} else {
step_max = interpolate(steps[last], lod(steps[last]), steps[last + 1], lod(steps[last + 1]), threshold);
}
}