Skip to content

GitLab

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

closes [#3437]

parent 9f008bbc
Hide whitespace changes
Inline Side-by-side
include/bayes/output.h
View file @ cbe99d99
......@@ -802,7 +802,13 @@ struct LV_database {
std::map<std::string, std::vector<MatrixXd>> ret;
for (const auto& chr_gen_lv_vec: data) {
const std::string& chr = chr_gen_lv_vec.first;
const auto& lv_vec = chr_gen_lv_vec.second.find(gen)->second;
auto it = chr_gen_lv_vec.second.find(gen);
if (it == chr_gen_lv_vec.second.end()) {
/* whine */
MSG_ERROR("Generation " << gen << " doesn't exist!", "Double-check the generation names given in commandline.");
continue;
}
const auto& lv_vec = it->second;
auto& ret_lv_vec = ret[chr];
ret_lv_vec.reserve(ind_vec.size());
for (size_t i: ind_vec) {
......
include/settings.h
View file @ cbe99d99
......@@ -218,6 +218,9 @@ struct settings_t {
std::map<std::string, trait_metadata_type> trait_metadata;
double Rjoin;
double Rskip;
settings_t()
: notes()
, map_filename(), map()
......@@ -269,6 +272,8 @@ struct settings_t {
, with_lg()
, covar_per_pop()
, lod_support(1.)
, Rjoin(10.)
, Rskip(1.)
{}
std::vector<std::pair<const chromosome*, double>>
......@@ -596,7 +601,11 @@ settings_t::finalize()
MatrixXd xtx = MatrixXd::Zero(pleio.traits.size(), pleio.traits.size());
for (const auto& vp: pops) {
for (const auto& p: vp) {
xtx += p->observed_traits.front().XtX();
if (!p->observed_traits.front().raw.size()) {
MSG_ERROR("Invalid dataset", "");
} else {
xtx += p->observed_traits.front().XtX();
}
}
}
auto P = trait::do_PCA(xtx, pleio.tolerance);
......@@ -633,6 +642,10 @@ settings_t::finalize()
}
}
} else {
if (!pptr->gen) {
MSG_ERROR("PROUT", "");
continue;
}
context_key ck(new context_key_struc(pptr->gen.get()));
locus_key lk;
size_t n_par = get_n_parents(ck, lk);
......
src/bayes/cli.cc
View file @ cbe99d99
......@@ -560,6 +560,14 @@ bn_settings_t* bn_settings_t::from_args(int argc, const char** argv)
if (ret->output_mode == 0) {
ret->output_mode = bn_settings_t::OutputPopData;
}
if (ret->output_generations.size() == 0) {
for (const auto& kv: ret->pedigree.generation_names) {
ret->output_generations.emplace_back(kv.second);
}
}
std::sort(ret->output_generations.begin(), ret->output_generations.end());
auto it = std::unique(ret->output_generations.begin(), ret->output_generations.end());
ret->output_generations.resize(it - ret->output_generations.begin());
}
return ret;
......
src/main.cc
View file @ cbe99d99
......@@ -146,6 +146,7 @@ int main(int argc, const char** argv)
active_settings->set_title("Checking the validity of the configuration");
if (!active_settings->sanity_check()) {
MSG_QUEUE_FLUSH();
exit(-1);
}
......@@ -165,7 +166,7 @@ int main(int argc, const char** argv)
collection<int> all_lp;
for (const auto& chr_interval: all_loci) {
chromosome_value chr = chr_interval.first;
const auto& loci = chr_interval.second.front().all_positions;
const auto& loci = chr_interval.second./*front().*/all_positions;
locus_key empty_lk;
ensure_directories_exist(MESSAGE(active_settings->work_directory << '/' << active_settings->name << ".n-point"));
// for (const std::string& gen: active_settings->npoint_gen) {
......@@ -297,6 +298,11 @@ int main(int argc, const char** argv)
MSG_DEBUG("initial selection " << active_settings->initial_selection << " qtl detection algorithm " << active_settings->qtl_algorithm);
MSG_INFO(__FILE__ << ':' << __LINE__ << " mm.selection = " << mm.get_selection());
auto sel = mm.get_selection();
MSG_INFO(__FILE__ << ':' << __LINE__ << " mm.selection = " << mm.get_selection());
mm.search_intervals = full_search_intervals(sel);
MSG_INFO(__FILE__ << ':' << __LINE__ << " mm.selection = " << mm.get_selection());
if (active_settings->qtl_algorithm != "none") {
active_settings->set_title(
MESSAGE("QTL detection (" << active_settings->qtl_algorithm << ')'));
......@@ -306,8 +312,6 @@ int main(int argc, const char** argv)
MSG_INFO(__FILE__ << ':' << __LINE__ << " mm.selection = " << mm.get_selection());
MSG_INFO("QTL acceptance threshold set to " << mm.threshold);
MSG_INFO(__FILE__ << ':' << __LINE__ << " mm.selection = " << mm.get_selection());
auto sel = mm.get_selection();
mm.search_intervals = full_search_intervals(sel);
MSG_INFO(__FILE__ << ':' << __LINE__ << " mm.selection = " << mm.get_selection());
//*
......@@ -339,16 +343,42 @@ int main(int argc, const char** argv)
mm.add(active_settings->find_chromosome("chr7"), 96.89);
//*/
/*mm = mm.model_for_estimation();*/
MSG_INFO(__FILE__ << ':' << __LINE__ << " mm.selection = " << mm.get_selection());
mm.vMcurrent->compute();
/*MSG_DEBUG("MODEL BEFORE CONSTRAINTS" << std::endl << (*mm.vMcurrent));*/
mm.vMcurrent->enable_constraints();
mm.vMcurrent->compute();
/*mm = mm.model_for_estimation();*/
if (active_settings->initial_selection == "") {
MSG_INFO(__FILE__ << ':' << __LINE__ << " mm.selection = " << mm.get_selection());
// MSG_DEBUG("MODEL BEFORE CONSTRAINTS" << std::endl << (*mm.vMcurrent));
/*MSG_DEBUG("MODEL AFTER CONSTRAINTS" << std::endl << (*mm.vMcurrent));*/
MSG_INFO(__FILE__ << ':' << __LINE__ << " mm.selection = " << mm.get_selection());
mm.vMcurrent->enable_constraints();
MSG_INFO(__FILE__ << ':' << __LINE__ << " mm.selection = " << mm.get_selection());
mm.vMcurrent->compute();
MSG_INFO(__FILE__ << ':' << __LINE__ << " mm.selection = " << mm.get_selection());
ar.report_final_model(mm);
// MSG_DEBUG("MODEL AFTER CONSTRAINTS" << std::endl << (*mm.vMcurrent));
MSG_INFO(__FILE__ << ':' << __LINE__ << " mm.selection = " << mm.get_selection());
ar.report_final_model(mm);
MSG_INFO(__FILE__ << ':' << __LINE__ << " mm.selection = " << mm.get_selection());
} else {
model raw(*mm.vMcurrent);
for (const auto& chr: active_settings->map) {
raw.remove_blocks_with_chromosome(&chr);
}
raw.compute();
double test = 0;
if (mm.vMcurrent->Y().cols() > 1) {
VectorXd tmp = chi2(raw, *mm.vMcurrent, mm.vMcurrent->Y().cols());
test = tmp(0);
} else {
MatrixXd tmp(1, 1);
f_test(raw, *mm.vMcurrent, 0, &tmp, NULL);
test = tmp(0, 0);
}
MSG_DEBUG("TEST: " << test);
}
ar.detach_model_manager(mm);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment