Commit 4a6d29c5 authored by Damien Leroux's avatar Damien Leroux
Browse files

Fixed issue with locus_keys getting empty when challenging a QTL.

parent 4b160a3f
This diff is collapsed.
......@@ -1276,20 +1276,21 @@ struct search_lg_type {
deselect(double position, const collection<population_value>& all_pops, value<model> M)
{
auto i = M->m_blocks.begin(), j = M->m_blocks.end();
for (; i != j;) {
std::vector<decltype(M->m_blocks.begin())> to_remove;
for (; i != j; ++i) {
auto& key = i->first;
auto& value = i->second;
if (key->has_locus(chrom, position)) {
if (!key->remove_position(chrom, position)) {
auto tmp = i;
++i;
M->m_blocks.erase(tmp);
} else {
reduce(all_pops, position, value);
++i;
}
to_remove.emplace_back(i);
}
}
for (auto it: to_remove) {
auto& key = it->first;
if (!key->remove_position(chrom, position)) {
M->m_blocks.erase(it);
} else {
++i;
reduce(all_pops, position, it->second);
}
}
//
......@@ -1995,10 +1996,12 @@ struct model_manager {
void
restore_chromosome_selection(chromosome_value chrom, const std::pair<model_block_collection, std::vector<locus_key>>& chrom_sel)
{
MSG_DEBUG("before restore_chromosome_selection " << vMcurrent->keys());
auto& si = search_intervals.find(chrom)->second;
si.selection = chrom_sel.second.front();
si.recompute_modes();
vMcurrent->add_blocks(chrom_sel.first);
MSG_DEBUG("after restore_chromosome_selection " << vMcurrent->keys());
// auto sii = search_intervals.find(chrom)->second.begin();
// for (const auto& lk: chrom_sel.second) {
// sii->selection = lk;
......@@ -2184,6 +2187,7 @@ struct model_manager {
challenge_qtl(chromosome_value chrom, double position)
{
std::pair<bool, test_result> ret;
MSG_DEBUG("before challenge_qtl " << vMcurrent->keys());
#if 0
model_manager mm(*this);
mm.deselect(chrom, position);
......@@ -2193,12 +2197,13 @@ struct model_manager {
ret.first = (ret.second.locus == position && ret.second.over_threshold);
#else
deselect(chrom, position);
MSG_DEBUG("Using model " << vMcurrent->keys() << " to challenge " << chrom << ':' << position);
MSG_DEBUG("Using model " << vMcurrent->keys() << " to challenge " << chrom->name << ':' << position);
MSG_QUEUE_FLUSH();
ret.second = search_new_best(false, chrom);
ret.first = (ret.second.locus == position && ret.second.over_threshold);
select(chrom, position);
#endif
MSG_DEBUG("after challenge_qtl " << vMcurrent->keys());
return ret;
}
};
......
......@@ -208,6 +208,19 @@ struct model_block_key_struc {
}
}
bool
is_valid()
{
switch (type) {
case mbk_Covar:
case mbk_CI: return true;
case mbk_POP: return !loci->is_empty();
case mbk_Dominance: return left->is_valid();
case mbk_Interaction: return left->is_valid() && right->is_valid();
};
return false; /* make compiler happy. */
}
friend
std::ostream&
operator << (std::ostream& os, const model_block_key_struc& mbk)
......@@ -366,6 +379,10 @@ struct model_block_key_struc {
bool b = right->remove_position(chr, locus);
return a && b;
};
if (!is_valid()) {
MSG_ERROR("Created an invalid key!", "");
throw 0;
}
return false;
}
......@@ -975,40 +992,40 @@ struct model {
return ret;
/*} else if (mbk.selection.size() == 1) {*/
} else if (mbk->type == mbk_POP) {
size_t n_qtl = mb.column_labels.front().size();
// size_t n_qtl = mb.column_labels.front().size();
// if (n_qtl == 1) {
if (1) {
/* only one chromosome: do the contrast groups trick */
ret = {{{mbk, contrast_groups(m_all_pops, mbk->loci)}}};
} else {
MSG_DEBUG("Creating constraints for " << mbk);
std::vector<std::set<char>> uniq_letters_per_qtl(mb.column_labels.front().size());
/* need the epistasis magic here */
for (const auto& vec: mb.column_labels) {
auto i = uniq_letters_per_qtl.begin();
for (char c: vec) { i->insert(c); ++i; }
}
std::vector<size_t> letter_counts(uniq_letters_per_qtl.size());
for (size_t i = 0; i < letter_counts.size(); ++i) { letter_counts[i] = uniq_letters_per_qtl[i].size(); }
size_t start = 1;
size_t finish = (1 << letter_counts.size()) - 1;
for (size_t variant = start; variant < finish; ++variant) {
MatrixXd constraint = MatrixXd::Identity(1, 1);
for (size_t i = 0; i < letter_counts.size(); ++i) {
bool flat = (variant >> i) & 1;
MatrixXd tmp;
if (flat) {
tmp = kroneckerProduct(constraint, MatrixXd::Ones(1, letter_counts[i]));
} else {
tmp = kroneckerProduct(constraint, MatrixXd::Identity(letter_counts[i], letter_counts[i]));
}
constraint = tmp;
}
ret.emplace_back();
MSG_DEBUG("Created constraint for " << mbk);
MSG_DEBUG(constraint);
ret.back().insert({{mbk, constraint}});
}
// } else {
// MSG_DEBUG("Creating constraints for " << mbk);
// std::vector<std::set<char>> uniq_letters_per_qtl(mb.column_labels.front().size());
// /* need the epistasis magic here */
// for (const auto& vec: mb.column_labels) {
// auto i = uniq_letters_per_qtl.begin();
// for (char c: vec) { i->insert(c); ++i; }
// }
// std::vector<size_t> letter_counts(uniq_letters_per_qtl.size());
// for (size_t i = 0; i < letter_counts.size(); ++i) { letter_counts[i] = uniq_letters_per_qtl[i].size(); }
// size_t start = 1;
// size_t finish = (1 << letter_counts.size()) - 1;
// for (size_t variant = start; variant < finish; ++variant) {
// MatrixXd constraint = MatrixXd::Identity(1, 1);
// for (size_t i = 0; i < letter_counts.size(); ++i) {
// bool flat = (variant >> i) & 1;
// MatrixXd tmp;
// if (flat) {
// tmp = kroneckerProduct(constraint, MatrixXd::Ones(1, letter_counts[i]));
// } else {
// tmp = kroneckerProduct(constraint, MatrixXd::Identity(letter_counts[i], letter_counts[i]));
// }
// constraint = tmp;
// }
// ret.emplace_back();
// MSG_DEBUG("Created constraint for " << mbk);
// MSG_DEBUG(constraint);
// ret.back().insert({{mbk, constraint}});
// }
}
} else if (mbk->type == mbk_Interaction) {
/*MSG_DEBUG("Computing constraint for interaction " << mbk);*/
......@@ -1093,6 +1110,14 @@ struct model {
add_block(const model_block_key& key, const value<model_block_type>& x)
{
/*MSG_DEBUG("add_block(" << key << ", <some block>)");*/
if (!x) {
MSG_ERROR("Attempting to insert a null block in the model.", "");
abort();
}
if (!key->is_valid()) {
MSG_ERROR("Attempting to insert an invalid block key in the model.", "");
abort();
}
if (m_blocks.find(key) != m_blocks.end()) {
/*MSG_DEBUG("Key already exists");*/
return;
......@@ -1134,6 +1159,10 @@ struct model {
MSG_ERROR("Attempting to insert a null block in the model.", "");
abort();
}
if (!kv.first->is_valid()) {
MSG_ERROR("Attempting to insert an invalid block key in the model.", "");
abort();
}
}
m_blocks.insert(mbc.begin(), mbc.end());
}
......@@ -1255,7 +1284,7 @@ struct model {
void compute()
{
MSG_DEBUG('[' << std::this_thread::get_id() << "] Recomputing model with selection " << keys());
// MSG_DEBUG('[' << std::this_thread::get_id() << "] Recomputing model with selection " << keys());
if (m_computed) {
return;
}
......@@ -1272,7 +1301,7 @@ struct model {
}
}
}
MSG_DEBUG('[' << std::this_thread::get_id() << "] Filtered selection " << keys());
// MSG_DEBUG('[' << std::this_thread::get_id() << "] Filtered selection " << keys());
m_computed = true;
m_X = new immediate_value<MatrixXd>(MatrixXd());
......@@ -1316,9 +1345,9 @@ struct model {
n_cols += mb.outerSize();
}
n_rows = n_pop_rows;
MSG_DEBUG("n_rows=" << n_rows << " n_cols=" << n_cols);
// MSG_DEBUG("n_rows=" << n_rows << " n_cols=" << n_cols);
// MSG_DEBUG("DEBUG X" << std::endl << (*m_X));
{
if (0) {
/*os << "<model Y(" << m.m_Y->innerSize() << ',' << m.m_Y->outerSize() << "), " << m.m_blocks.size() << " blocks>";*/
MatrixXd big = MatrixXd::Zero(1, X().cols());
model_print::matrix_with_sections<std::string, void, std::string, std::vector<char>> mws(big);
......@@ -1335,8 +1364,8 @@ struct model {
auto it = cmap.find(b.first);
if (it != cmap.end()) {
const MatrixXd& constraint = it->second;
MSG_DEBUG('[' << std::this_thread::get_id() << "] Adding (" << constraint.innerSize() << ',' << constraint.outerSize() << ") constraint at " << n_rows << "," << n_cols << " within (" << m_X->rows() << ',' << m_X->cols() << ')');
MSG_DEBUG(constraint);
// MSG_DEBUG('[' << std::this_thread::get_id() << "] Adding (" << constraint.innerSize() << ',' << constraint.outerSize() << ") constraint at " << n_rows << "," << n_cols << " within (" << m_X->rows() << ',' << m_X->cols() << ')');
// MSG_DEBUG(constraint);
if (constraint.size()) {
m_X->block(n_rows, n_cols, constraint.innerSize(), constraint.outerSize()) = constraint;
row_incr = constraint.innerSize();
......
......@@ -177,9 +177,10 @@ VectorXd chi2(const model& model_current, const model& model_new, size_t block_s
int dof_cur = model_current.rank();
int dof_new = model_new.rank();
double Chi2;
assert(!(model_current.Y().outerSize() % block_size) && "STOOPID. BLOCK SIZE MUST BE A DIVISOR OF THE NUMBER OF COLUMNS IN Y.");
VectorXd ret(model_current.Y().outerSize() / block_size);
for (int c = 0, ofs = 0; ofs < model_current.Y().outerSize(); ofs += block_size, ++c) {
MatrixXd Y = model_current.Y();
assert(!(Y.outerSize() % block_size) && "STOOPID. BLOCK SIZE MUST BE A DIVISOR OF THE NUMBER OF COLUMNS IN Y.");
VectorXd ret(Y.outerSize() / block_size);
for (int c = 0, ofs = 0; ofs < Y.outerSize(); ofs += block_size, ++c) {
Chi2 = 0;
for (size_t i = ofs; i < (ofs + block_size); ++i) {
double rcur = rss_cur(i);
......
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