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

Fixed a bug in connected model blocks, assigning wrong column indices.

parent a38284dc
......@@ -927,9 +927,20 @@ struct qtl_pop_type {
}
ret->indices = filter_vector(indices, keep_indices);
/*MSG_DEBUG("Filtered indices " << ret->indices << " (" << ret->indices.size() << ')');*/
for (const auto& kv: LV) {
/*MSG_DEBUG("Filtering LV on " << kv.first << ", keeping " << keep_indices.size() << " elements");*/
ret->LV[kv.first] = filter_vector(kv.second, keep_indices);
}
/*MSG_DEBUG("Resulting LV (" << ret->LV.size() << ')');*/
/*for (const auto& lv: ret->LV) {*/
/*MSG_DEBUG("** " << lv.first << " (" << lv.second.size() << ')');*/
/*auto ri = ret->indices.begin();*/
/*for (const auto& v: lv.second) {*/
/*MSG_DEBUG("* " << (*ri++));*/
/*MSG_DEBUG("" << v);*/
/*}*/
/*}*/
return ret;
}
......
......@@ -1313,13 +1313,13 @@ void analysis_report::report_final_model(model_manager& mm)
/*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
<< "RSS " << mm.Mcurrent.rss() << std::endl
<< std::endl;
*/
//*/
if (mm.Mcurrent.n_obs() <= mm.Mcurrent.dof() - 1) {
report_file << "* There are not enough observations to compute a variance/covariance matrix." << std::endl;
......@@ -1874,7 +1874,7 @@ QTL::confidence_interval(const std::string& trait, const std::vector<QTL>& selec
auto lod = [&] (double l) { locus_set point = {{chr, l}}; return tm.results[point]; };
double threshold = lod(locus) - 1.5;
double threshold = lod(locus) - 1;
{
LOD.clear();
......
......@@ -167,15 +167,19 @@ init_connected_block_builder(
std::function<std::vector<L> (const qtl_pop_type&)> get_labels,
const collection<population_value>& all_pops)
{
std::set<L> uniq;
std::map<L, int> index;
for (auto& pop: all_pops) {
for (const L& l: get_labels(**pop)) {
int sz = index.size();
index.insert({l, sz});
}
uniq.insert(l);
}
}
for (const L& l: uniq) {
int sz = index.size();
index.insert({l, sz});
}
#if 0
#if 1
MSG_DEBUG("index size " << index.size());
std::stringstream s;
s << '{';
......@@ -188,9 +192,10 @@ init_connected_block_builder(
ret.n_rows = 0;
ret.n_columns = index.size();
ret.labels.reserve(ret.n_columns);
ret.elements.reserve(all_pops.size());
for (auto& kv: index) { ret.labels.push_back(kv.first); }
/*ret.labels.reserve(ret.n_columns);*/
/*for (auto& kv: index) { ret.labels.push_back(kv.first); }*/
ret.labels.assign(uniq.begin(), uniq.end());
for (auto& kv: all_pops) {
ret.elements.emplace_back();
......
......@@ -92,7 +92,7 @@ inline bool
arg_match(std::istream& is, const context_key& ck)
{
static std::string null("null");
return ((ck->pop ? arg_match(is, ck->pop->name) && arg_match(is, ck->pop->qtl_generation_name) && arg_match(is, ck->pop->size())
return ((ck->pop ? arg_match(is, ck->pop->name) && arg_match(is, ck->pop->qtl_generation_name) && arg_match(is, ck->pop->indices)
: arg_match(is, null))
&& (ck->gen ? arg_match(is, ck->gen) : arg_match(is, null))
&& (ck->chr ? arg_match(is, ck->chr->name) : arg_match(is, null))
......@@ -194,7 +194,7 @@ inline void
if (ck->pop) {
arg_write(os, ck->pop->name);
arg_write(os, ck->pop->qtl_generation_name);
arg_write(os, ck->pop->size());
arg_write(os, ck->pop->indices);
} else {
arg_write(os, null);
}
......
......@@ -146,7 +146,7 @@ l ocus_probabilities(qtl_chromosome_value qtl_chr, population_value pop,
optim_segment_computer_t
make_segment_computer(const chromosome* chr, population_value pop, const locus_key& lk, int ind)
make_segment_computer(const chromosome* chr, population_value pop, const locus_key& lk, int ind, size_t)
{
/*qtl_chromosome qc(ck->chr, lk);*/
const geno_matrix* gen = pop->gen.get();
......@@ -163,7 +163,7 @@ locus_probabilities(const context_key& ck, const locus_key& lk,
/*MSG_DEBUG("locus_probabilities(" << ck << ", " << lk << ", " << ind << ", " << loci << ")");*/
/*MSG_QUEUE_FLUSH();*/
/*const MatrixXd& LV = ck->pop->get_LV(ck->chr->name, ind);*/
auto ret = make_value<Mem|Sync>(make_segment_computer, as_value(ck->chr), as_value(ck->pop), as_value(lk), as_value(ind))->compute(loci);
auto ret = make_value<Mem|Sync>(make_segment_computer, as_value(ck->chr), as_value(ck->pop), as_value(lk), as_value(ind), as_value(ck->pop->indices[ind]))->compute(loci);
/*MSG_DEBUG("locus_probabilities(" << ck << ", " << lk << ", " << ind << ", " << loci << ")" << std::endl << ret);*/
return ret;
}
......
Supports Markdown
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