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

Now featuring dominance blocks in model. Not optional yet.

parent cfdec1f8
......@@ -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
......
......@@ -5,7 +5,7 @@
#include "computations/basic_data.h"
#include "computations/probabilities.h"
#include "computations/model.h"
#include "computations/frontends2.h"
#include "computations/frontends4.h"
#endif
......@@ -84,6 +84,9 @@ test_loci(chromosome_value);
labelled_matrix<MatrixXd, std::vector<char>, label_type>
compute_state_to_parental_origin(const context_key& ck, const locus_key& lk);
labelled_matrix<MatrixXd, std::vector<char>, label_type>
compute_state_to_dominance(const context_key& ck, const locus_key& lk);
MatrixXb
get_contrast_groups(generation_value gen, const locus_key& lk);
......
This diff is collapsed.
......@@ -252,6 +252,13 @@ compute_parental_origins_multipop(const collection<population_value>& all_pops,
const value<std::vector<double>>& loci,
const std::vector<double>& test_loci);
collection<model_block_type>
compute_dominance_multipop(const collection<population_value>& all_pops,
const value<chromosome_value>& chr,
const value<locus_key>& lk,
const value<std::vector<double>>& loci,
const std::vector<double>& test_loci);
/*collection<std::vector<char>>*/
/*contrast_groups_connected();*/
......@@ -280,7 +287,8 @@ assemble_parental_origins_multipop(
const value<chromosome_value>& chr,
const value<locus_key>& lk,
const std::vector<collection<parental_origin_per_locus_type>>& all_popl,
const collection<population_value>& all_pops);
const collection<population_value>& all_pops,
bool block_is_POP);
collection<parental_origin_per_locus_type>
disassemble_parental_origins_multipop(
......
......@@ -30,6 +30,9 @@ parental_origin(const locus_probabilities_type&, generation_value, const selecte
parental_origin_per_locus_type
joint_parental_origin_at_locus(const context_key& ck, const locus_key& lk);
parental_origin_per_locus_type
joint_dominance_at_locus(const context_key& ck, const locus_key& lk, bool);
collection<parental_origin_per_locus_type>
parental_origin_per_locus(const collection<parental_origin_type>&);
......
......@@ -166,6 +166,8 @@ struct labelled_matrix {
size_t innerSize() const { return data.innerSize(); }
size_t outerSize() const { return data.outerSize(); }
int cols() const { return data.cols(); }
int rows() const { return data.rows(); }
struct formatted {
const labelled_matrix<MATRIX, row_label_type, col_label_type>* m;
......
This diff is collapsed.
......@@ -166,6 +166,9 @@ struct settings_t {
std::unordered_map<std::thread::id, std::vector<std::string>> thread_stacks;
size_t max_order;
bool cross_indicator_can_interact;
settings_t()
: notes()
, map_filename(), map()
......@@ -204,6 +207,8 @@ struct settings_t {
, LV()
, linked_pops()
, thread_stacks()
, max_order(1)
, cross_indicator_can_interact(false)
{}
std::vector<std::pair<const chromosome*, double>>
......
......@@ -77,7 +77,7 @@ qtl_generation(const population_value& pop)
#else
generation_value
generation_by_name(const population_value& pop, const std::string& name)
generation_by_name(const population_value& pop, const std::string& /*name*/)
{
return pop->gen.get();
}
......@@ -115,7 +115,7 @@ qtl_chromosome_by_name(const std::string& name)
collection<std::string>
all_observed_generation_names(population_value pop)
all_observed_generation_names(population_value /*pop*/)
{
collection<std::string> ret;
/* FIXME */
......@@ -291,6 +291,41 @@ compute_state_to_parental_origin(const context_key& ck, const locus_key& lk)
}
labelled_matrix<MatrixXd, std::vector<char>, label_type>
compute_state_to_dominance(const context_key& ck, const locus_key& lk)
{
auto stdom = compute_state_to_parental_origin(ck, lk);
std::map<std::vector<char>, VectorXd> dominance_vectors;
int nrows = stdom.rows();
for (int r1 = 0; r1 < nrows; ++r1) {
for (int r2 = r1 + 1; r2 < nrows; ++r2) {
VectorXd prod = (stdom.data.row(r1).array() * stdom.data.row(r2).array()).matrix();
if (prod.sum() > 0) {
std::vector<char> label;
label.push_back('(');
label.insert(label.end(), stdom.row_labels[r1].begin(), stdom.row_labels[r1].end());
label.push_back('/');
label.insert(label.end(), stdom.row_labels[r2].begin(), stdom.row_labels[r2].end());
label.push_back(')');
dominance_vectors[label] = prod;
}
}
}
stdom.row_labels.clear();
stdom.row_labels.reserve(dominance_vectors.size());
stdom.data.resize(dominance_vectors.size(), stdom.data.cols());
for (const auto& lv: dominance_vectors) {
stdom.data.row(stdom.row_labels.size()) = lv.second;
stdom.row_labels.emplace_back(lv.first);
}
MSG_DEBUG("Computed stdom matrix" << std::endl << stdom);
return stdom;
}
stpom_data
compute_state_to_parental_origin_haplo(const context_key& ck, const locus_key& lk)
{
......
......@@ -179,11 +179,22 @@ iqtlm_backward_one(model_manager& mm)
for (const auto& chr_sel: mm.get_selection()) {
for (double d: chr_sel.second) {
results.emplace_back(mm.challenge_qtl(chr_sel.first, d));
std::pair<bool, test_result> result = mm.challenge_qtl(chr_sel.first, d);
MSG_DEBUG("Challenging QTL at " << chr_sel.first << ':' << d << " gives " << result.first << ' ' << result.second);
results.emplace_back(result);
not_modified &= results.back().first;
}
}
{
std::stringstream dbg;
dbg << "backward";
for (const auto& bt: results) {
dbg << " {" << (bt.first ? "true" : "false") << ' ' << bt.second << '}';
}
MSG_DEBUG(dbg.str());
}
if (not_modified) {
return false;
}
......@@ -210,11 +221,15 @@ in_history(std::set<std::map<chromosome_value, locus_key>>& history, const std::
bool
iqtlm_backward(model_manager& mm, std::set<std::map<chromosome_value, locus_key>>& history)
{
DUMP_FILE_LINE();
while (iqtlm_backward_one(mm)) {
DUMP_FILE_LINE();
if (in_history(history, mm.get_selection())) {
DUMP_FILE_LINE();
return false;
}
}
DUMP_FILE_LINE();
return true;
}
......@@ -222,19 +237,27 @@ iqtlm_backward(model_manager& mm, std::set<std::map<chromosome_value, locus_key>
bool
iqtlm_forward(model_manager& mm, std::set<std::map<chromosome_value, locus_key>>& history)
{
DUMP_FILE_LINE();
while (true) {
DUMP_FILE_LINE();
auto results = mm.search_new_best_per_chromosome();
bool modified = false;
for (const auto& result: results) {
DUMP_FILE_LINE();
modified |= result.second.over_threshold;
if (result.second.over_threshold) {
DUMP_FILE_LINE();
result.second.select(mm);
MSG_DEBUG("added " << result);
}
}
DUMP_FILE_LINE();
if (!modified) {
DUMP_FILE_LINE();
return true;
}
if (in_history(history, mm.get_selection())) {
DUMP_FILE_LINE();
return false;
}
}
......@@ -266,14 +289,14 @@ qtl_detect_iqtlm(model_manager& mm)
std::set<std::map<chromosome_value, locus_key>> history;
std::map<chromosome_value, locus_key> last_sel, new_sel;
new_sel = mm.get_selection();
while (new_sel != last_sel) {
do {
if (iqtlm_backward(mm, history) && iqtlm_forward(mm, history)) {
last_sel = new_sel;
new_sel = mm.get_selection();
} else {
break;
}
}
} while (new_sel != last_sel);
return mm;
}
......
......@@ -167,7 +167,7 @@ cross_indicator(const std::string& trait)
}
}
model_block_type val;
int p = 0;
/*int p = 0;*/
std::vector<int> columns;
const auto& lp = active_settings->linked_pops;
std::vector<int> pop_col(lp.size(), -1);
......@@ -184,7 +184,13 @@ cross_indicator(const std::string& trait)
}
val.column_labels.clear();
val.column_labels.resize(val.data.outerSize(), {' '});
val.column_labels.reserve(val.data.outerSize());
for (const auto& pop_vec: pops) {
for (const auto& pop: pop_vec) {
val.column_labels.emplace_back(pop->qtl_generation_name.begin(), pop->qtl_generation_name.end());
}
}
/*for (const auto& pop: pops) {*/
/*int sz = (*pop)->observed_traits.front().values.size();*/
/*val.data.block(n_rows, p, sz, 1) = VectorXd::Ones(sz);*/
......@@ -410,7 +416,7 @@ init_model(const value<MatrixXd>& Y, const value<model_block_type>& indicator)
model M(Y, all_populations());
const char* decomp_method = getenv("DECOMP_METHOD");
std::string s(decomp_method ? decomp_method : "");
M.add_block({}, indicator);
M.add_block(model_block_key_struc::cross_indicator(), indicator);
/*std::cout << "indicator" << std::endl << indicator->transpose() << std::endl;*/
/*std::cout << "Xt" << std::endl << M.X().transpose() << std::endl;*/
/*std::cout << "Yt" << std::endl << M.Y().transpose() << std::endl;*/
......@@ -470,7 +476,8 @@ assemble_parental_origins_multipop(
const value<chromosome_value>& chr,
const value<locus_key>& lk,
const std::vector<collection<parental_origin_per_locus_type>>& all_popl,
const collection<population_value>& all_pop)
const collection<population_value>& all_pop,
bool block_is_POP)
{
/* FIXME take into account the ACTUAL list of populations used! */
/*MSG_DEBUG("ALL_POPL " << (*chr)->name << " " << lk);*/
......@@ -487,14 +494,29 @@ assemble_parental_origins_multipop(
};
std::function<std::vector<std::vector<char>>(const qtl_pop_type&)>
get_labels = [=](const qtl_pop_type& pop)
{
context_key ck(new context_key_struc(&pop, *chr, std::vector<double>()));
/*MSG_DEBUG(ck);*/
return make_value<Sync|Mem>(compute_state_to_parental_origin,
value<context_key>{ck},
lk)->row_labels;
};
get_labels;
if (block_is_POP) {
get_labels
= [=](const qtl_pop_type& pop)
{
context_key ck(new context_key_struc(&pop, *chr, std::vector<double>()));
/*MSG_DEBUG(ck);*/
return make_value<Sync|Mem>(compute_state_to_parental_origin,
value<context_key>{ck},
lk)->row_labels;
};
} else {
get_labels
= [=](const qtl_pop_type& pop)
{
context_key ck(new context_key_struc(&pop, *chr, std::vector<double>()));
/*MSG_DEBUG(ck);*/
return make_value<Sync|Mem>(compute_state_to_dominance,
value<context_key>{ck},
lk)->row_labels;
};
}
/*DUMP_FILE_LINE();*/
block_builder<std::vector<char>> bb;
......@@ -656,6 +678,91 @@ compute_effective_loci(const locus_key& lk, const std::vector<double>& loci)
}
template <typename Ret, typename... Args>
collection<model_block_type>
_compute_parental_origins_multipop(Ret (&Joint_POP_at_locus)(Args...),
const collection<population_value>& all_pop,
const value<chromosome_value>& chr,
const value<locus_key>& lk,
const value<std::vector<double>>& loci,
const std::vector<double>& test_loci)
{
/*DUMP_FILE_LINE();*/
std::vector<collection<parental_origin_per_locus_type>> all_popl;
/*DUMP_FILE_LINE();*/
all_popl.reserve(all_pop.size());
/*DUMP_FILE_LINE();*/
/*MSG_DEBUG("loci are " << loci);*/
/*MSG_DEBUG("given " << lk << ", test loci are " << test_loci);*/
collection<locus_key> lkeys;
lkeys.reserve(test_loci.size());
for (double l: test_loci) {
lkeys.emplace_back(*lk + l);
}
for (auto& p: all_pop) {
/*DUMP_FILE_LINE();*/
context_key ck(new context_key_struc(*p, *chr, *loci));
all_popl.emplace_back(make_collection<Disk>(Joint_POP_at_locus,
value<context_key>{ck}, lkeys));
/*std::cout << all_popl << std::endl;*/
/*all_popl.emplace_back(compute_parental_origins(p, (*qtl_chr)->chr, loci));*/
}
for (auto& pl: all_popl) {
for (auto& l: pl) {
(void) *l;
}
}
return assemble_parental_origins_multipop(chr, lk, all_popl, all_pop, &Joint_POP_at_locus == &joint_parental_origin_at_locus);
}
collection<model_block_type>
compute_dominance_multipop(const collection<population_value>& all_pop,
const value<chromosome_value>& chr,
const value<locus_key>& lk,
/*const value<qtl_chromosome_value>& qtl_chr,*/
const value<std::vector<double>>& loci,
const std::vector<double>& test_loci)
{
/*DUMP_FILE_LINE();*/
std::vector<collection<parental_origin_per_locus_type>> all_popl;
/*DUMP_FILE_LINE();*/
all_popl.reserve(all_pop.size());
/*DUMP_FILE_LINE();*/
/*MSG_DEBUG("loci are " << loci);*/
/*MSG_DEBUG("given " << lk << ", test loci are " << test_loci);*/
collection<locus_key> lkeys;
lkeys.reserve(test_loci.size());
for (double l: test_loci) {
lkeys.emplace_back(*lk + l);
}
for (auto& p: all_pop) {
/*DUMP_FILE_LINE();*/
context_key ck(new context_key_struc(*p, *chr, *loci));
all_popl.emplace_back(make_collection<Disk>(joint_dominance_at_locus,
value<context_key>{ck}, lkeys, as_value(false)));
/*std::cout << all_popl << std::endl;*/
/*all_popl.emplace_back(compute_parental_origins(p, (*qtl_chr)->chr, loci));*/
}
for (auto& pl: all_popl) {
for (auto& l: pl) {
(void) *l;
}
}
return assemble_parental_origins_multipop(chr, lk, all_popl, all_pop, false);
}
collection<model_block_type>
compute_parental_origins_multipop(const collection<population_value>& all_pop,
const value<chromosome_value>& chr,
......@@ -664,6 +771,7 @@ compute_parental_origins_multipop(const collection<population_value>& all_pop,
const value<std::vector<double>>& loci,
const std::vector<double>& test_loci)
{
#if 1
/*locus_key lk(new locus_key_struc((*qtl_chr)->qtl));*/
/*int n_qtl = (*qtl_chr)->qtl.size();*/
......@@ -720,7 +828,10 @@ compute_parental_origins_multipop(const collection<population_value>& all_pop,
(void) *l;
}
}
return assemble_parental_origins_multipop(chr, lk, all_popl, all_pop);
return assemble_parental_origins_multipop(chr, lk, all_popl, all_pop, true);
#else
return _compute_parental_origins_multipop(joint_parental_origin_at_locus, all_pop, chr, lk, loci, test_loci);
#endif
}
......
......@@ -315,6 +315,29 @@ joint_geno_prob_at_locus(const context_key& ck, const locus_key& lk)
}
/* Compute the joint probability of selected loci UNION new locus */
parental_origin_per_locus_type
joint_dominance_at_locus(const context_key& ck, const locus_key& lk, bool)
{
value<context_key> vck{ck};
value<locus_key> vlk{lk};
value<locus_key> vlkp{lk->parent};
auto pop = make_value<Disk>(joint_geno_prob_at_locus, vck, vlk);
auto stdom = make_value<Mem>(compute_state_to_dominance, vck, vlkp);
/*(void) *pop;*/
/*(void) *stfopom;*/
/*MSG_DEBUG("joint_dominance_at_locus");*/
/*MSG_DEBUG("******* " << MATRIX_SIZE(*pop));*/
/*MSG_DEBUG(pop);*/
/*MSG_DEBUG("******* " << MATRIX_SIZE(*stdom));*/
/*MSG_DEBUG(stdom);*/
/*MSG_QUEUE_FLUSH();*/
return ((*stdom) * (*pop)).transpose();
}
/* Compute the joint probability of selected loci UNION new locus */
parental_origin_per_locus_type
joint_parental_origin_at_locus(const context_key& ck, const locus_key& lk)
......
......@@ -927,6 +927,17 @@ arguments = {
SAFE_IGNORE_CALLBACK_ARGS;
}},
{{"max-order"},
{"order"},
"Maximum block order allowed in the model",
false,
{1},
[](CALLBACK_ARGS)
{
ensure(target)->max_order = to<size_t>(*++ai);
SAFE_IGNORE_CALLBACK_ARGS;
}},
{{"epistasis"},
{},
"Detect epistasis",
......
......@@ -147,10 +147,11 @@ int main(int argc, const char** argv)
model_manager mm(trait_pops.first, pops,
trait_matrix(trait_pops.first, pops),
active_settings->get_threshold(trait_pops.first) * .9, /* FIXME, get_threshold needs pops_by_trait, so maybe move pops_by_trait to settings */
ComputationType::FTest,
//active_settings->get_threshold(trait_pops.first) * .9, /* FIXME, get_threshold needs pops_by_trait, so maybe move pops_by_trait to settings */
/*ComputationType::FTest,*/
active_settings->max_order, /* FIXME make max_order configurable in settings_t */
SolverType::QR);
mm.threshold = active_settings->get_threshold(trait_pops.first) * .9;
ar.attach_model_manager(mm);
#if 1
......@@ -174,7 +175,7 @@ int main(int argc, const char** argv)
} else if (active_settings->cofactor_algorithm == "backward") {
backward(mm);
} else if (active_settings->cofactor_algorithm == "all") {
mm.add_all_loci();
mm.select_all_loci();
}
MSG_QUEUE_FLUSH();
//*/
......@@ -211,8 +212,8 @@ int main(int argc, const char** argv)
mm.add(active_settings->find_chromosome("chr7"), 104.87);
mm = mm.model_for_estimation();
mm.Mcurrent.enable_constraints();
mm.Mcurrent.compute();
mm.vMcurrent->enable_constraints();
mm.vMcurrent->compute();
ar.report_final_model(mm);
......@@ -221,9 +222,14 @@ int main(int argc, const char** argv)
mm.add(active_settings->find_chromosome("chr7"), 96.89);
//*/
mm = mm.model_for_estimation();
mm.Mcurrent.enable_constraints();
mm.Mcurrent.compute();
/*mm = mm.model_for_estimation();*/
mm.vMcurrent->compute();
MSG_DEBUG("MODEL BEFORE CONSTRAINTS" << std::endl << (*mm.vMcurrent));
mm.vMcurrent->enable_constraints();
mm.vMcurrent->compute();
MSG_DEBUG("MODEL AFTER CONSTRAINTS" << std::endl << (*mm.vMcurrent));
ar.report_final_model(mm);
......@@ -243,7 +249,7 @@ int main(int argc, const char** argv)
MSG_DEBUG("TEST MANAGER SAYS");
MSG_DEBUG("" << tm.results);
auto ri = tm.results.begin(), rj = tm.results.end();
auto ri = tm.results.begin(); //, rj = tm.results.end();
double v0, v1, l0, l1;
for (const chromosome& c: active_settings->map) {
v0 = v1 = 0;
......@@ -267,12 +273,12 @@ int main(int argc, const char** argv)
}
#if 0
mm.Mcurrent.output_X_to_file();
mm.Mcurrent.output_XtX_inv_to_file();
mm.vMcurrent->output_X_to_file();
mm.vMcurrent->output_XtX_inv_to_file();
/*mm.QTLs_to_cofactors();*/
/*mm.Mcurrent.compute();*/
mm.Mcurrent.output_X_to_file();
mm.Mcurrent.output_XtX_inv_to_file();
/*mm.vMcurrent->compute();*/
mm.vMcurrent->output_X_to_file();
mm.vMcurrent->output_XtX_inv_to_file();
/*MSG_DEBUG("Trait " << trait_pops.first << std::endl << mm.Mcurrent);*/
MSG_DEBUG("Final model");
MSG_DEBUG("" << mm.Mcurrent);
......@@ -280,9 +286,9 @@ int main(int argc, const char** argv)
MSG_DEBUG("pmode = " << mm.pmode);
model::xtx_printable xtx(mm.Mcurrent, true);
MSG_DEBUG("XtX^-1" << std::endl << xtx);
MSG_DEBUG("Coefficients " << mm.Mcurrent.coefficients().transpose());
/*MSG_DEBUG("Residuals " << mm.Mcurrent.residuals().transpose());*/
MSG_DEBUG("RSS " << mm.Mcurrent.rss());
MSG_DEBUG("Coefficients " << mm.vMcurrent->coefficients().transpose());
/*MSG_DEBUG("Residuals " << mm.vMcurrent->residuals().transpose());*/
MSG_DEBUG("RSS " << mm.vMcurrent->rss());
/*MSG_DEBUG("Model " << std::endl << mm.Mcurrent);*/
MSG_QUEUE_FLUSH();
......@@ -298,7 +304,7 @@ int main(int argc, const char** argv)
}
#if 0
for (const auto& mbk: mm.Mcurrent.m_blocks) {
for (const auto& mbk: mm.vMcurrent->m_blocks) {
for (const auto& chr_lk: mbk.first.selection) {
for (double l: chr_lk.second) {
if (poi[chr_lk.first->name][l].size()) {
......
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