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

fixes [#3434].

parent 658049f5
......@@ -9,9 +9,15 @@
<file path="$PROJECT_DIR$/tests" />
</sourceRoots>
<excludeRoots>
<file path="$PROJECT_DIR$/PAG_XXV" />
<file path="$PROJECT_DIR$/attic" />
<file path="$PROJECT_DIR$/cmake-build-debug/CMakeFiles" />
<file path="$PROJECT_DIR$/cmake-build-release/CMakeFiles" />
<file path="$PROJECT_DIR$/data" />
<file path="$PROJECT_DIR$/include/attic" />
<file path="$PROJECT_DIR$/sample-data" />
<file path="$PROJECT_DIR$/sandbox" />
<file path="$PROJECT_DIR$/simulation" />
<file path="$PROJECT_DIR$/simulation/gen_BC2" />
<file path="$PROJECT_DIR$/simulation/gen_F9" />
<file path="$PROJECT_DIR$/simulation/gen_MAGIC4" />
......@@ -20,10 +26,12 @@
<file path="$PROJECT_DIR$/simulation/gen_RIL.old" />
<file path="$PROJECT_DIR$/simulation/gen_SIB" />
<file path="$PROJECT_DIR$/simulation/test/gen_RIL" />
<file path="$PROJECT_DIR$/simulator" />
<file path="$PROJECT_DIR$/src/bayes/gros_tests" />
<file path="$PROJECT_DIR$/src/bayes/test-data" />
<file path="$PROJECT_DIR$/src/bayes/tmp" />
<file path="$PROJECT_DIR$/src/experiments" />
<file path="$PROJECT_DIR$/test_ril_ril2" />
<file path="$PROJECT_DIR$/tests/TestCodeCoverage" />
</excludeRoots>
</component>
......
This diff is collapsed.
......@@ -18,8 +18,8 @@
#ifndef _SPELL_BETA_GAMMA_H_
#define _SPELL_BETA_GAMMA_H_
double ibeta(double a, double b, double c);
double gamma_q(double a, double b);
double log_ibeta(double a, double b, double c);
double log_gamma_q(double a, double b);
#endif
......
This diff is collapsed.
......@@ -269,7 +269,7 @@ init_disconnected_block_builder(
collection<model_block_type>
compute_parental_origins_multipop(const collection<population_value>& all_pops,
const value<chromosome_value>& chr,
const value<chromosome_value>& chr,
const value<locus_key>& lk,
const value<std::vector<double>>& loci,
const std::vector<double>& test_loci);
......@@ -281,6 +281,20 @@ compute_dominance_multipop(const collection<population_value>& all_pops,
const value<std::vector<double>>& loci,
const std::vector<double>& test_loci);
collection<model_block_type>
compute_parental_origins_multipop(const collection<population_value>& all_pops,
const value<chromosome_value>& chr,
const std::vector<locus_key>& lk,
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 std::vector<locus_key>& lk,
const value<std::vector<double>>& loci,
const std::vector<double>& test_loci);
/*collection<std::vector<char>>*/
/*contrast_groups_connected();*/
......@@ -472,34 +486,34 @@ void
init_computation(computation_along_chromosome& ret, ComputationType ct, ComputationResults cr, size_t n_steps, size_t Y_rows, size_t Y_cols)
{
if (cr & RSS) {
ret.rss.resize(Y_cols, n_steps);
ret.rss = MatrixXd::Zero(Y_cols, n_steps);
}
if (cr & Residuals) {
ret.residuals.resize(Y_rows, Y_cols * n_steps);
ret.residuals = MatrixXd::Zero(Y_rows, Y_cols * n_steps);
}
if (cr & Coefficients) {
ret.residuals.resize(Y_rows, Y_cols * n_steps);
ret.residuals = MatrixXd::Zero(Y_rows, Y_cols * n_steps);
}
if (cr & Rank) {
ret.rank.resize(n_steps);
ret.rank = VectorXd::Zero(n_steps);
}
if (ct & FTest) {
ret.ftest_pvalue.resize(Y_cols, n_steps);
ret.ftest_pvalue = MatrixXd::Zero(Y_cols, n_steps);
}
if (ct & FTestLOD) {
ret.ftest_lod.resize(Y_cols, n_steps);
ret.ftest_lod = MatrixXd::Zero(Y_cols, n_steps);
}
if (ct & Chi2) {
ret.chi2.resize(1, n_steps);
ret.chi2 = MatrixXd::Zero(1, n_steps);
}
if (ct & Chi2LOD) {
ret.chi2_lod.resize(1, n_steps);
ret.chi2_lod = MatrixXd::Zero(1, n_steps);
}
if (ct & Mahalanobis) {
ret.mahalanobis.resize(n_steps, 1);
ret.mahalanobis = MatrixXd::Zero(n_steps, 1);
}
if (ct & R2) {
ret.r2.resize(Y_cols, n_steps);
ret.r2 = MatrixXd::Zero(Y_cols, n_steps);
}
}
......@@ -507,13 +521,13 @@ init_computation(computation_along_chromosome& ret, ComputationType ct, Computat
template <CachingPolicy _Policy = CachingPolicy::Oneshot>
void
compute_along_interval(int i0, computation_along_chromosome& ret,
value<ComputationType> vct, value<ComputationResults> vcr,
size_t y_block_cols,
const value<model>& Mcurrent, const value<model>& M0,
const locus_key& base_key,
chromosome_value chr,
const std::vector<double>& loci,
const collection<parental_origin_per_locus_type>& popl)
value<ComputationType> vct, value<ComputationResults> vcr,
size_t y_block_cols,
const value<model>& Mcurrent, const value<model>& M0,
const locus_key& base_key,
chromosome_value chr,
const std::vector<double>& loci,
const collection<parental_origin_per_locus_type>& popl)
{
value<computation_along_chromosome*> pcac = &ret;
......@@ -534,13 +548,58 @@ compute_along_interval(int i0, computation_along_chromosome& ret,
as_value(y_block_cols),
Mcurrent, M0,
vmbk,
/*as_value(base_key + std::make_pair(chr, loci[i])),*/
/*as_value(base_key + std::make_pair(chr, loci[i])),*/
popl[i]));
}
for (auto& k: joiner) {
DUMP_FILE_LINE();
(void) *k; /* join them all */
MSG_DEBUG("Block" << std::endl << (*k));
MSG_QUEUE_FLUSH();
}
}
template <CachingPolicy _Policy = CachingPolicy::Oneshot>
void
compute_along_interval(int i0, computation_along_chromosome& ret,
value<ComputationType> vct, value<ComputationResults> vcr,
size_t y_block_cols,
const value<model>& Mcurrent, const value<model>& M0,
const std::vector<locus_key>& base_key,
chromosome_value chr,
const std::vector<double>& loci,
const collection<parental_origin_per_locus_type>& popl)
{
value<computation_along_chromosome*> pcac = &ret;
/*DUMP_FILE_LINE();*/
collection<int> joiner;
joiner.reserve(popl.size());
/* FIXME: create correct keys? FIXED */
for (auto i: range<int>(i0, i0 + popl.size(), 1)) {
/*DUMP_FILE_LINE();*/
/*value<model_block_key> vmbk = model_block_key(base_key);*/
/**vmbk += std::make_pair(chr, loci[*i]);*/
value<model_block_key> vmbk = model_block_key_struc::pop(chr, base_key[*i] + loci[*i]);
/*MSG_DEBUG("new block key " << vmbk);*/
joiner.push_back(make_value<_Policy>(compute_one,
i, pcac, /* flower */
vct, vcr,
as_value(y_block_cols),
Mcurrent, M0,
vmbk,
/*as_value(base_key + std::make_pair(chr, loci[i])),*/
popl[i]));
}
for (auto& k: joiner) {
// DUMP_FILE_LINE();
(void) *k; /* join them all */
// MSG_DEBUG("Block" << std::endl << (*k));
// MSG_QUEUE_FLUSH();
}
}
......
......@@ -105,7 +105,7 @@ struct chromosome {
void compute_haplo_sizes()
{
MSG_DEBUG("Computing haplotype sizes for chromosome " << name);
// MSG_DEBUG("Computing haplotype sizes for chromosome " << name);
haplo_sizes.clear();
condensed.marker_locus.clear();
condensed.marker_name.clear();
......
......@@ -1704,6 +1704,11 @@ struct model {
void
add_block_with_interactions(model_block_key mbk, const value<model_block_type>& block, const model* refptr=NULL)
{
if (!block) {
MSG_ERROR("Trying to add NULL block", "");
MSG_QUEUE_FLUSH();
abort();
}
if (mbk->order() > m_max_order) {
return;
}
......
......@@ -118,7 +118,8 @@ void f_test(const model& model_current, const model& model_new, int col_num, Mat
if (F <= 0) {
(*pvalue)(i, col_num) = 0;
} else {
(*pvalue)(i, col_num) = ibeta(dof_denom * .5, dof_num * .5, dof_denom / (dof_denom + dof_num * F));
(*pvalue)(i, col_num) = log_ibeta(dof_denom * .5, dof_num * .5,
dof_denom / (dof_denom + dof_num * F));
}
}
if (lod) {
......@@ -159,8 +160,8 @@ void f_test(const model& model_current, const model& model_new, int col_num, Mat
// (*ret)(col_num, 0) = 0;
// /*return 0;*/
// } else {
// (*ret)(col_num, 0) = gamma_q(df * .5, Chi2 * .5);
// /*return -log10(boost::math::gamma_q(df * .5, Chi2 * .5));*/
// (*ret)(col_num, 0) = log_gamma_q(df * .5, Chi2 * .5);
// /*return -log10(boost::math::log_gamma_q(df * .5, Chi2 * .5));*/
// }
//}
......@@ -190,12 +191,14 @@ VectorXd chi2(const model& model_current, const model& model_new, size_t block_s
Chi2 += log(rcur / rnew);
}
int df = block_size * (dof_new - dof_cur);
// MSG_DEBUG("chi2=" << Chi2);
if (Chi2 <= 0 || df <= 0) {
ret(c) = 0;
} else if (LOD) {
ret(c) = .5 / log(10) * Chi2;
ret(c) = .5 / log(10) * model_new.n_obs() * Chi2;
} else {
ret(c) = gamma_q(df * .5, Chi2 * .5);
ret(c) = log_gamma_q(df * .5, model_new.n_obs() * Chi2 * .5);
// ret(c) = model_new.n_obs() * Chi2;
}
}
return ret;
......
......@@ -19,12 +19,12 @@
#include <boost/math/special_functions/gamma.hpp>
#include <cmath>
double ibeta(double a, double b, double c)
double log_ibeta(double a, double b, double c)
{
return -log10(boost::math::ibeta(a, b, c));
}
double gamma_q(double a, double b)
double log_gamma_q(double a, double b)
{
return -log10(boost::math::gamma_q(a, b));
}
......
......@@ -419,28 +419,30 @@ int compute_one(int i, computation_along_chromosome* ret,
const parental_origin_per_locus_type& popl)
{
/*DUMP_FILE_LINE();*/
/*MSG_DEBUG("one popl" << std::endl << popl.transpose());*/
MSG_DEBUG("one popl" << std::endl << popl.transpose());
model Mext = M0.extend(k, popl);
Mext.compute();
/*MSG_DEBUG("Rank(Mext " << Mext.keys() << ") = " << Mext.rank());*/
/*MSG_DEBUG('[' << Mext.keys() << " rank=" << Mext.rank() << " rss=" << Mext.rss().transpose() << ']');*/
// MSG_DEBUG("Rank(Mext " << Mext.keys() << ") = " << Mext.rank());
// MSG_DEBUG('[' << Mext.keys() << " rank=" << Mext.rank() << " rss=" << Mext.rss().transpose() << ']');
/*DUMP_FILE_LINE();*/
// DUMP_FILE_LINE();
if (cr & RSS) {
one_computation_result(i, ret->rss, Mext.rss());
}
/*DUMP_FILE_LINE();*/
// DUMP_FILE_LINE();
if (cr & Residuals) {
one_computation_result(i, ret->residuals, Mext.residuals());
}
/*DUMP_FILE_LINE();*/
// DUMP_FILE_LINE();
if (cr & Coefficients) {
one_computation_result(i, ret->coefficients, Mext.coefficients());
}
// DUMP_FILE_LINE();
if (cr & Rank) {
one_computation_result(i, ret->rank, Mext.rank());
}
// DUMP_FILE_LINE();
if (ct & (FTest|FTestLOD)) {
f_test(Mcurrent, Mext, i,
ct & FTest ? &ret->ftest_pvalue : NULL,
......@@ -453,12 +455,15 @@ int compute_one(int i, computation_along_chromosome* ret,
}
#endif
}
// DUMP_FILE_LINE();
if (ct & Mahalanobis) {
mahalanobis(Mcurrent, Mext, i, &ret->mahalanobis);
}
// DUMP_FILE_LINE();
if (ct & R2) {
r2(Mcurrent, Mext, i, &ret->r2);
}
// DUMP_FILE_LINE();
if (ct & Chi2) {
// auto tmp = chi2(Mcurrent, Mext, y_block_size);
// MSG_DEBUG("HAVE CHI2 COLUMN " << tmp.rows() << 'x' << tmp.cols());
......@@ -466,6 +471,15 @@ int compute_one(int i, computation_along_chromosome* ret,
ret->chi2.col(i) = chi2(Mcurrent, Mext, y_block_size);
// ret->chi2.col(i) = chi2(Mcurrent, Mext, y_block_size).transpose();
}
if (ct & Chi2LOD) {
// auto tmp = chi2(Mcurrent, Mext, y_block_size);
// MSG_DEBUG("HAVE CHI2 COLUMN " << tmp.rows() << 'x' << tmp.cols());
// MSG_DEBUG("HAVE recipient " << ret->chi2.rows() << 'x' << ret->chi2.cols());
// MSG_DEBUG("Chi2LOD");
ret->chi2_lod.col(i) = chi2(Mcurrent, Mext, y_block_size, true);
// ret->chi2.col(i) = chi2(Mcurrent, Mext, y_block_size).transpose();
}
// DUMP_FILE_LINE();
return 0;
}
......@@ -479,7 +493,7 @@ init_model(const value<MatrixXd>& Y, const value<model_block_type>& indicator, c
// std::string s(decomp_method ? decomp_method : "");
M.m_max_order = 1;
M.add_block(model_block_key_struc::cross_indicator(), indicator);
MSG_DEBUG("covariables " << active_settings->with_covariables << " covars.cols=" << covars->data.cols());
// MSG_DEBUG("covariables " << active_settings->with_covariables << " covars.cols=" << covars->data.cols());
if (covars->data.cols() > 0) {
M.add_block(model_block_key_struc::covariables(), covars);
}
......@@ -548,76 +562,160 @@ assemble_parental_origins_multipop(
/* FIXME take into account the ACTUAL list of populations used! */
/*MSG_DEBUG("ALL_POPL " << (*chr)->name << " " << lk);*/
/*for (const auto& popl: all_popl) {*/
/*MSG_DEBUG("======================================");*/
/*MSG_DEBUG("" << popl);*/
/*MSG_DEBUG("======================================");*/
/*MSG_DEBUG("" << popl);*/
/*}*/
/*MSG_QUEUE_FLUSH();*/
static
std::function<int(const qtl_pop_type&)>
get_pop_size = [](const qtl_pop_type& pop)
{
return pop.size();
};
static
std::function<int(const qtl_pop_type&)>
get_pop_size = [](const qtl_pop_type& pop)
{
return pop.size();
};
std::function<std::vector<std::vector<char>>(const qtl_pop_type&)>
get_labels;
get_labels;
if (block_is_POP) {
get_labels
= [=](const qtl_pop_type& pop)
= [=](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;
if (active_settings->connected) {
/*MSG_DEBUG("CONNECTED!");*/
init_connected_block_builder(bb, get_pop_size, get_labels, all_pop);
} else {
/*MSG_DEBUG("DISCONNECTED!");*/
init_disconnected_block_builder(bb, get_pop_size, get_labels, all_pop);
}
/*DUMP_FILE_LINE();*/
collection<model_block_type> all_blocks;
all_blocks.reserve(all_popl.front().size());
size_t n = all_popl.front().size();
for (size_t i = 0; i < n; ++i) {
/*DUMP_FILE_LINE();*/
all_blocks.emplace_back(new immediate_value<model_block_type>(model_block_type()));
/*DUMP_FILE_LINE();*/
all_blocks.back()->column_labels = bb.labels;
auto builder = bb.begin(all_blocks.back()->data);
/*DUMP_FILE_LINE();*/
/*MSG_DEBUG("building at #" << i);*/
for (auto& popl: all_popl) {
/*MSG_DEBUG("popl" << std::endl << popl);*/
/*DUMP_FILE_LINE();*/
/*MSG_DEBUG(" block " << std::endl << popl[i]->data);*/
/*MSG_QUEUE_FLUSH();*/
builder += popl[i]->data;
}
/*MSG_DEBUG('#' << i << std::endl << all_blocks.back());*/
}
/*MSG_DEBUG("ALL BLOCKS" << std::endl << all_blocks);*/
return all_blocks;
}
collection<model_block_type>
assemble_parental_origins_multipop(
const value<chromosome_value>& chr,
const std::vector<locus_key>& lk,
const std::vector<collection<parental_origin_per_locus_type>>& all_popl,
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);*/
/*for (const auto& popl: all_popl) {*/
/*MSG_DEBUG("======================================");*/
/*MSG_DEBUG("" << popl);*/
/*}*/
/*MSG_QUEUE_FLUSH();*/
static
std::function<int(const qtl_pop_type&)>
get_pop_size = [](const qtl_pop_type& pop)
{
return pop.size();
};
std::function<std::vector<std::vector<char>>(const qtl_pop_type&)>
get_labels;
/*DUMP_FILE_LINE();*/
collection<model_block_type> all_blocks;
all_blocks.reserve(all_popl.front().size());
size_t n = all_popl.front().size();
for (size_t i = 0; i < n; ++i) {
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;
as_value(lk[i]))->row_labels;
};
} else {
get_labels
= [=](const qtl_pop_type& pop)
} 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;
as_value(lk[i]))->row_labels;
};
}
}
/*DUMP_FILE_LINE();*/
block_builder<std::vector<char>> bb;
if (active_settings->connected) {
/*MSG_DEBUG("CONNECTED!");*/
init_connected_block_builder(bb, get_pop_size, get_labels, all_pop);
} else {
/*MSG_DEBUG("DISCONNECTED!");*/
init_disconnected_block_builder(bb, get_pop_size, get_labels, all_pop);
}
/*DUMP_FILE_LINE();*/
collection<model_block_type> all_blocks;
all_blocks.reserve(all_popl.front().size());
size_t n = all_popl.front().size();
for (size_t i = 0; i < n; ++i) {
/*DUMP_FILE_LINE();*/
all_blocks.emplace_back(new immediate_value<model_block_type>(model_block_type()));
block_builder<std::vector<char>> bb;
if (active_settings->connected) {
/*MSG_DEBUG("CONNECTED!");*/
init_connected_block_builder(bb, get_pop_size, get_labels, all_pop);
} else {
/*MSG_DEBUG("DISCONNECTED!");*/
init_disconnected_block_builder(bb, get_pop_size, get_labels, all_pop);
}
/*DUMP_FILE_LINE();*/
all_blocks.emplace_back(new immediate_value<model_block_type>(model_block_type()));
/*DUMP_FILE_LINE();*/
all_blocks.back()->column_labels = bb.labels;
auto builder = bb.begin(all_blocks.back()->data);
auto builder = bb.begin(all_blocks.back()->data);
/*DUMP_FILE_LINE();*/
/*MSG_DEBUG("building at #" << i);*/
for (auto& popl: all_popl) {
/*MSG_DEBUG("building at #" << i);*/
for (auto& popl: all_popl) {
/*MSG_DEBUG("popl" << std::endl << popl);*/
/*DUMP_FILE_LINE();*/
/*MSG_DEBUG(" block " << std::endl << popl[i]->data);*/
/*MSG_DEBUG(" block " << std::endl << popl[i]->data);*/
/*MSG_QUEUE_FLUSH();*/
builder += popl[i]->data;
}
/*MSG_DEBUG('#' << i << std::endl << all_blocks.back());*/
}
builder += popl[i]->data;
}
/*MSG_DEBUG('#' << i << std::endl << all_blocks.back());*/
}
/*MSG_DEBUG("ALL BLOCKS" << std::endl << all_blocks);*/
return all_blocks;
return all_blocks;
}
......@@ -789,16 +887,16 @@ _compute_parental_origins_multipop(Ret (&Joint_POP_at_locus)(Args...),
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)
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;
std::vector<collection<parental_origin_per_locus_type>> all_popl;
/*DUMP_FILE_LINE();*/
all_popl.reserve(all_pop.size());
all_popl.reserve(all_pop.size());
/*DUMP_FILE_LINE();*/
/*MSG_DEBUG("loci are " << loci);*/
......@@ -810,15 +908,15 @@ compute_dominance_multipop(const collection<population_value>& all_pop,
lkeys.emplace_back(*lk + l);
}
for (auto& p: all_pop) {
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));*/
/*all_popl.emplace_back(compute_parental_origins(p, (*qtl_chr)->chr, loci));*/
}
}
for (auto& pl: all_popl) {
for (auto& l: pl) {
(void) *l;
......@@ -829,11 +927,53 @@ compute_dominance_multipop(const collection<population_value>& all_pop,
}
collection<model_block_type>
compute_dominance_multipop(const collection<population_value>& all_pop,
const value<chromosome_value>& chr,
const std::vector<locus_key>& lk,
/*const value<qtl_chromosome_value>& qtl_chr,*/
const value<std::vector<double>>& loci,