Commit f3153233 authored by damien's avatar damien
Browse files
parents dce8adf2 b211f878
#!/bin/bash
BUILD_ROOT=$1
cd ../doc/user_manual && make BINDIR=$BUILD_ROOT/bin SRCDIR=$2 && cp user_manual.pdf $BUILD_ROOT/doc
...@@ -77,10 +77,11 @@ struct bn_settings_t { ...@@ -77,10 +77,11 @@ struct bn_settings_t {
static const int OutputPopData = 1; static const int OutputPopData = 1;
static const int OutputOnePoint = 2; static const int OutputOnePoint = 2;
static const int OutputGametes = 4;
int output_mode; int output_mode;
std::vector<std::string> closing_messages; std::vector<std::string> closing_messages;
bn_settings_t() bn_settings_t()
: scheme(JDS_None) : scheme(JDS_None)
, command_line() , command_line()
......
...@@ -134,7 +134,8 @@ struct state_index_type { ...@@ -134,7 +134,8 @@ struct state_index_type {
std::ostream& std::ostream&
operator << (std::ostream& os, const state_index_type& s) operator << (std::ostream& os, const state_index_type& s)
{ {
for (int i = n_blocks - 1; i >= 0; ++i) { MSG_DEBUG("printing state_index " << s.data);
for (int i = n_blocks - 1; i >= 0; --i) {
uint64_t mask = 1ULL << 63; uint64_t mask = 1ULL << 63;
while (mask && !(s.data[i] & mask)) { mask >>= 1; } while (mask && !(s.data[i] & mask)) { mask >>= 1; }
while (mask) { os << ((s.data[i] & mask) == 0); mask >>= 1; } while (mask) { os << ((s.data[i] & mask) == 0); mask >>= 1; }
...@@ -525,13 +526,13 @@ struct joint_variable_product_type { ...@@ -525,13 +526,13 @@ struct joint_variable_product_type {
bool bool
move_cursor_if_state_is_valid(table_descr& t, state_index_type cursor) move_cursor_if_state_is_valid(table_descr& t, state_index_type cursor)
{ {
/*MSG_DEBUG("table @" << t.data << " cursor=" << dump(cursor));*/ MSG_DEBUG("table @" << t.data << " cursor=" << dump(cursor));
auto it = t.domain_to_offset.find(t.local(cursor)); auto it = t.domain_to_offset.find(t.local(cursor));
if (it == t.domain_to_offset.end()) { if (it == t.domain_to_offset.end()) {
/*MSG_DEBUG("not found");*/ MSG_DEBUG("not found");
return false; return false;
} }
/*MSG_DEBUG("found");*/ MSG_DEBUG("found");
t.offset = it->second; t.offset = it->second;
return true; return true;
} }
...@@ -560,6 +561,7 @@ struct joint_variable_product_type { ...@@ -560,6 +561,7 @@ struct joint_variable_product_type {
void void
compile(const std::map<std::vector<int>, genotype_comb_type>& all_variable_domains) compile(const std::map<std::vector<int>, genotype_comb_type>& all_variable_domains)
{ {
scoped_indent _("[jvp compile] ");
compile_domains(all_variable_domains); compile_domains(all_variable_domains);
compile_tables(); compile_tables();
compile_output(); compile_output();
...@@ -569,6 +571,8 @@ struct joint_variable_product_type { ...@@ -569,6 +571,8 @@ struct joint_variable_product_type {
void void
compile_domains(const std::map<std::vector<int>, genotype_comb_type>& all_variable_domains) compile_domains(const std::map<std::vector<int>, genotype_comb_type>& all_variable_domains)
{ {
scoped_indent _("[domain] ");
MSG_DEBUG(all_variable_domains);
all_domains = &all_variable_domains; all_domains = &all_variable_domains;
total_bits = 0; total_bits = 0;
for (int v: all_variable_names) { for (int v: all_variable_names) {
...@@ -580,7 +584,7 @@ struct joint_variable_product_type { ...@@ -580,7 +584,7 @@ struct joint_variable_product_type {
d.bit_width = 1; d.bit_width = 1;
while ((1ULL << d.bit_width) < d.size) { ++d.bit_width; } while ((1ULL << d.bit_width) < d.size) { ++d.bit_width; }
total_bits += d.bit_width; total_bits += d.bit_width;
/*MSG_DEBUG("Have domain for variable #" << d.variable_name << " size " << d.size << " bit_width " << d.bit_width);*/ MSG_DEBUG("Have domain for variable #" << d.variable_name << " size " << d.size << " bit_width " << d.bit_width);
/*MSG_QUEUE_FLUSH();*/ /*MSG_QUEUE_FLUSH();*/
d.index_to_label.reserve(domain.size()); d.index_to_label.reserve(domain.size());
for (size_t i = 0; i < domain.size(); ++i) { for (size_t i = 0; i < domain.size(); ++i) {
...@@ -599,9 +603,9 @@ struct joint_variable_product_type { ...@@ -599,9 +603,9 @@ struct joint_variable_product_type {
shift += d.bit_width; shift += d.bit_width;
d.assign(max_cursor, d.max); d.assign(max_cursor, d.max);
/*MSG_QUEUE_FLUSH();*/ /*MSG_QUEUE_FLUSH();*/
/*MSG_DEBUG("Have domain for variable #" << d.variable_name << " size " << d.size << " bit_width " << d.bit_width << " mask " << boost::dynamic_bitset<uint64_t>(total_bits, d.mask) << " max " << boost::dynamic_bitset<uint64_t>(total_bits, d.max));*/ MSG_DEBUG("Have domain for variable #" << d.variable_name << " size " << d.size << " bit_width " << d.bit_width << " mask " << boost::dynamic_bitset<uint64_t>(total_bits, d.mask) << " max " << boost::dynamic_bitset<uint64_t>(total_bits, d.max));
/*MSG_QUEUE_FLUSH();*/ /*MSG_QUEUE_FLUSH();*/
if (0) { if (1) {
std::stringstream ss; std::stringstream ss;
for (const auto& kv: d.label_to_index) { for (const auto& kv: d.label_to_index) {
ss << ' ' << kv.first << ':' << kv.second; ss << ' ' << kv.first << ':' << kv.second;
...@@ -616,6 +620,7 @@ struct joint_variable_product_type { ...@@ -616,6 +620,7 @@ struct joint_variable_product_type {
void void
compile_tables() compile_tables()
{ {
scoped_indent _("[tables] ");
std::vector<size_t> coordinates; std::vector<size_t> coordinates;
coordinates.reserve(all_variable_names.size()); coordinates.reserve(all_variable_names.size());
for (auto& tab: tables) { for (auto& tab: tables) {
...@@ -633,9 +638,9 @@ struct joint_variable_product_type { ...@@ -633,9 +638,9 @@ struct joint_variable_product_type {
tab.variables_left[v] = true; tab.variables_left[v] = true;
} }
/*MSG_QUEUE_FLUSH();*/ /*MSG_QUEUE_FLUSH();*/
/*MSG_DEBUG("Table " << tab.data);*/ MSG_DEBUG("Table " << *tab.data);
/*MSG_QUEUE_FLUSH();*/ /*MSG_QUEUE_FLUSH();*/
/*MSG_DEBUG("Table variables " << tab.variable_names << " mask " << tab.mask);*/ MSG_DEBUG("Table variables " << tab.variable_names << " mask " << tab.mask);
/*MSG_QUEUE_FLUSH();*/ /*MSG_QUEUE_FLUSH();*/
tab.offset = 0; tab.offset = 0;
coordinates.resize(tab.variable_names.size()); coordinates.resize(tab.variable_names.size());
...@@ -649,7 +654,7 @@ struct joint_variable_product_type { ...@@ -649,7 +654,7 @@ struct joint_variable_product_type {
tab.domain_to_offset[cursor] = i; tab.domain_to_offset[cursor] = i;
tab.unif = tab.unif && std::abs(coef - tab.data->m_combination[i].coef) > _EPSILON; tab.unif = tab.unif && std::abs(coef - tab.data->m_combination[i].coef) > _EPSILON;
} }
if (0) { if (1) {
std::stringstream ss; std::stringstream ss;
ss << "Table"; ss << "Table";
for (const auto& x: tab.domain_to_offset) { for (const auto& x: tab.domain_to_offset) {
...@@ -669,21 +674,23 @@ struct joint_variable_product_type { ...@@ -669,21 +674,23 @@ struct joint_variable_product_type {
void void
compile_output() compile_output()
{ {
/*MSG_DEBUG("all output variables " << output_variable_names);*/ scoped_indent _("[output] ");
MSG_DEBUG("all output variables " << output_variable_names);
output_variable_names = output_variable_names % all_variable_names; output_variable_names = output_variable_names % all_variable_names;
/*MSG_DEBUG("all_variable_names " << all_variable_names << " output_variable_names " << output_variable_names);*/ MSG_DEBUG("all_variable_names " << all_variable_names << " output_variable_names " << output_variable_names);
output_variables_in_use.resize(all_variable_names.size(), false); output_variables_in_use.resize(all_variable_names.size(), false);
for (int v: output_variable_names) { for (int v: output_variable_names) {
size_t index = std::find(all_variable_names.begin(), all_variable_names.end(), v) - all_variable_names.begin(); size_t index = std::find(all_variable_names.begin(), all_variable_names.end(), v) - all_variable_names.begin();
output_variables_in_use[index] = true; output_variables_in_use[index] = true;
MSG_DEBUG("variable " << v << " has index " << index);
for (size_t t = 0; t < tables.size(); ++t) { for (size_t t = 0; t < tables.size(); ++t) {
auto iter = std::find(tables[t].variable_names.begin(), tables[t].variable_names.end(), v); auto iter = std::find(tables[t].variable_names.begin(), tables[t].variable_names.end(), v);
if (iter != tables[t].variable_names.end()) { if (iter != tables[t].variable_names.end()) {
/*MSG_DEBUG("variable " << v << " found in table " << t << " at index " << (iter - tables[t].variable_names.begin()));*/ MSG_DEBUG("variable " << v << " found in table " << t << " at index " << (iter - tables[t].variable_names.begin()));
key_extractors.emplace_back(t, iter - tables[t].variable_names.begin()); key_extractors.emplace_back(t, iter - tables[t].variable_names.begin());
break; break;
/*} else {*/ } else {
/*MSG_DEBUG("variable " << v << " not found in table " << t);*/ MSG_DEBUG("variable " << v << " not found in table " << t);
} }
} }
} }
...@@ -836,9 +843,11 @@ struct joint_variable_product_type { ...@@ -836,9 +843,11 @@ struct joint_variable_product_type {
operator () (const joint_variable_product_type& jvp, genotype_comb_type::key_list& key) const operator () (const joint_variable_product_type& jvp, genotype_comb_type::key_list& key) const
{ {
auto ki = key.begin(); auto ki = key.begin();
MSG_DEBUG("extractors " << jvp.key_extractors);
for (const auto& t_ki: jvp.key_extractors) { for (const auto& t_ki: jvp.key_extractors) {
*ki++ = jvp.tables[t_ki.first].current_element().keys.keys[t_ki.second]; *ki++ = jvp.tables[t_ki.first].current_element().keys.keys[t_ki.second];
} }
MSG_DEBUG("key: " << key);
} }
}; };
...@@ -851,24 +860,24 @@ struct joint_variable_product_type { ...@@ -851,24 +860,24 @@ struct joint_variable_product_type {
auto ki = key.begin(); auto ki = key.begin();
auto ti = jvp.tables.begin(); auto ti = jvp.tables.begin();
bn_label_type G = ti->current_element().keys.keys.front().state; bn_label_type G = ti->current_element().keys.keys.front().state;
/*MSG_DEBUG("Have gamete label " << G); MSG_QUEUE_FLUSH();*/ MSG_DEBUG("Have gamete label " << G); MSG_QUEUE_FLUSH();
for (const auto& t_ki: jvp.key_extractors) { for (const auto& t_ki: jvp.key_extractors) {
*ki++ = jvp.tables[t_ki.first].current_element().keys.keys[t_ki.second]; *ki++ = jvp.tables[t_ki.first].current_element().keys.keys[t_ki.second];
/**ki++ = ti->current_element().keys.keys.front();*/ /**ki++ = ti->current_element().keys.keys.front();*/
} }
/*MSG_DEBUG("Have temp key " << key); MSG_QUEUE_FLUSH();*/ MSG_DEBUG("Have temp key " << key); MSG_QUEUE_FLUSH();
const auto& p1 = key.keys[G.first_allele]; const auto& p1 = key.keys[G.first_allele];
const auto& p2 = key.keys[G.second_allele]; const auto& p2 = key.keys[G.second_allele];
bool f1 = G.first == GAMETE_L; bool f1 = G.first == GAMETE_L;
bool f2 = G.second == GAMETE_L; bool f2 = G.second == GAMETE_L;
/*MSG_DEBUG("Have p1 " << p1 << " p2 " << p2 << " f1 " << f1 << " f2 " << f2); MSG_QUEUE_FLUSH();*/ MSG_DEBUG("Have p1 " << p1 << " p2 " << p2 << " f1 " << f1 << " f2 " << f2); MSG_QUEUE_FLUSH();
*ki = { *ki = {
jvp.output_variable_names.back(), {f1 ? p1.state.first : p1.state.second, jvp.output_variable_names.back(), {f1 ? p1.state.first : p1.state.second,
f2 ? p2.state.first : p2.state.second, f2 ? p2.state.first : p2.state.second,
f1 ? p1.state.first_allele : p1.state.second_allele, f1 ? p1.state.first_allele : p1.state.second_allele,
f2 ? p2.state.first_allele : p2.state.second_allele} f2 ? p2.state.first_allele : p2.state.second_allele}
}; };
/*MSG_DEBUG("Have complete key " << key); MSG_QUEUE_FLUSH();*/ MSG_DEBUG("Have complete key " << key); MSG_QUEUE_FLUSH();
} }
}; };
...@@ -976,6 +985,8 @@ struct joint_variable_product_type { ...@@ -976,6 +985,8 @@ struct joint_variable_product_type {
genotype_comb_type genotype_comb_type
compute_generic() compute_generic()
{ {
scoped_indent _("[compute jvp] ");
MSG_DEBUG("output variable names " << output_variable_names);
if (invalid_product) { return {}; } if (invalid_product) { return {}; }
if (total_bits >= 32) { if (total_bits >= 32) {
...@@ -1098,6 +1109,7 @@ struct joint_variable_product_type { ...@@ -1098,6 +1109,7 @@ struct joint_variable_product_type {
} }
std::sort(ret.m_combination.begin(), ret.m_combination.end()); std::sort(ret.m_combination.begin(), ret.m_combination.end());
MSG_DEBUG("Computed " << counter << " coefficients, projected onto " << ret.size() << " elements"); MSG_DEBUG("Computed " << counter << " coefficients, projected onto " << ret.size() << " elements");
MSG_DEBUG(ret);
return ret; return ret;
} }
...@@ -1184,9 +1196,9 @@ compute_product(Iterator tables_begin, Iterator tables_end, const std::vector<in ...@@ -1184,9 +1196,9 @@ compute_product(Iterator tables_begin, Iterator tables_end, const std::vector<in
inline inline
genotype_comb_type genotype_comb_type
generate_factor(const std::vector<int>& parents, bool dh, int spawnling, const genotype_comb_type& parent_domain, std::map<std::vector<int>, genotype_comb_type>& domains) generate_factor(const std::vector<int>& parents, bool dh, int spawnling, const genotype_comb_type& parent_domain, std::map<std::vector<int>, genotype_comb_type>& domains, bool keep_gametes)
{ {
/*scoped_indent _(MESSAGE("[generate_factor " << (dh ? "dh " : "") << spawnling << " from parents " << parents << " domain " << parent_domain << "] "));*/ scoped_indent _(MESSAGE("[generate_factor " << (dh ? "dh " : "") << spawnling << " from parents " << parents << " domain " << parent_domain << "] "));
std::vector<bn_label_type> label_g; std::vector<bn_label_type> label_g;
std::vector<int> sorted_parents(parents); std::vector<int> sorted_parents(parents);
if (parents.size() == 2) { if (parents.size() == 2) {
...@@ -1219,19 +1231,27 @@ generate_factor(const std::vector<int>& parents, bool dh, int spawnling, const g ...@@ -1219,19 +1231,27 @@ generate_factor(const std::vector<int>& parents, bool dh, int spawnling, const g
{GAMETE_R, GAMETE_R, 0, 0} {GAMETE_R, GAMETE_R, 0, 0}
}; };
} }
int gamete = -spawnling;
genotype_comb_type genotype_comb_type
G = state_to_combination(-1, label_g) * (1. / label_g.size()); G = state_to_combination(gamete, label_g) * (1. / label_g.size());
domains[{-1}] = G; MSG_DEBUG("before adding gamete domain " << domains);
domains[{-spawnling}] = G;
MSG_DEBUG("after adding gamete domain " << domains);
joint_variable_product_type jvp; joint_variable_product_type jvp;
/*MSG_DEBUG("joint parent domain " << parent_domain);*/ MSG_DEBUG("joint parent domain " << parent_domain);
MSG_DEBUG("all domains " << domains);
jvp.add_table(G); jvp.add_table(G);
jvp.add_table(parent_domain); jvp.add_table(parent_domain);
jvp.set_output(sorted_parents); if (keep_gametes) {
jvp.set_output(sorted_parents + std::vector<int>{gamete, spawnling});
} else {
jvp.set_output(sorted_parents + std::vector<int>{spawnling});
}
jvp.compile(domains); jvp.compile(domains);
jvp.output_variable_names.push_back(spawnling);
return jvp.build_factor(); return jvp.build_factor();
} }
......
...@@ -430,12 +430,14 @@ struct factor_graph_type : public recursive_graph_type<factor_graph_type> { ...@@ -430,12 +430,14 @@ struct factor_graph_type : public recursive_graph_type<factor_graph_type> {
std::vector<message_type> node_domains; std::vector<message_type> node_domains;
std::map<node_index_type, node_index_type> anchor_points; std::map<node_index_type, node_index_type> anchor_points;
std::vector<bool> node_domain_computed; std::vector<bool> node_domain_computed;
bool keep_gametes;
std::vector<int> gametes;
graph_type& operator = (graph_type&& other) = delete; graph_type& operator = (graph_type&& other) = delete;
factor_graph_type(const graph_type& other) = delete; factor_graph_type(const graph_type& other) = delete;
factor_graph_type(size_t n_al=1) factor_graph_type(bool keepgam, size_t n_al=1)
: io(), domains(), /*parent(nullptr), index_in_parent(0),*/ n_alleles(n_al), is_dh(), node_domains(), anchor_points(), node_domain_computed() : io(), domains(), /*parent(nullptr), index_in_parent(0),*/ n_alleles(n_al), is_dh(), node_domains(), anchor_points(), node_domain_computed(), keep_gametes(keepgam)
{} {}
std::shared_ptr<factor_graph_type> std::shared_ptr<factor_graph_type>
...@@ -472,9 +474,9 @@ struct factor_graph_type : public recursive_graph_type<factor_graph_type> { ...@@ -472,9 +474,9 @@ struct factor_graph_type : public recursive_graph_type<factor_graph_type> {
static static
std::unique_ptr<factor_graph_type> std::unique_ptr<factor_graph_type>
from_pedigree(const std::vector<pedigree_item>& ped, size_t n_al, const std::vector<std::string>& input_gen, const std::vector<std::string>& output_gen/*, const std::string& debug_prefix = ""*/) from_pedigree(const std::vector<pedigree_item>& ped, size_t n_al, bool keep_gametes, const std::vector<std::string>& input_gen, const std::vector<std::string>& output_gen/*, const std::string& debug_prefix = ""*/)
{ {
std::unique_ptr<factor_graph_type> g(new factor_graph_type(n_al)); std::unique_ptr<factor_graph_type> g(new factor_graph_type(keep_gametes, n_al));
/*system("rm -vf test_*.png");*/ /*system("rm -vf test_*.png");*/
for (const auto& p: filter_pedigree(ped, input_gen, output_gen)) { for (const auto& p: filter_pedigree(ped, input_gen, output_gen)) {
/*MSG_DEBUG("-------------------------------------------------------------------------------");*/ /*MSG_DEBUG("-------------------------------------------------------------------------------");*/
...@@ -1023,22 +1025,24 @@ struct factor_graph_type : public recursive_graph_type<factor_graph_type> { ...@@ -1023,22 +1025,24 @@ struct factor_graph_type : public recursive_graph_type<factor_graph_type> {
} else if (node_is_subgraph(n)) { } else if (node_is_subgraph(n)) {
subgraph(n)->compute_domains_and_factors(); subgraph(n)->compute_domains_and_factors();
} else /* factor */ { } else /* factor */ {
/*scoped_indent _(MESSAGE("[compute_factor {" << rule_of(n) << "} => " << own_variable_of(n) << "] "));*/ scoped_indent _(MESSAGE("[compute_factor {" << rule_of(n) << "} => " << own_variable_of(n) << "] "));
/*MSG_DEBUG(" ### ### COMPUTING FACTOR {" << rule_of(n) << "} => " << own_variable_of(n) << " ### ###");*/ MSG_DEBUG(" ### ### COMPUTING FACTOR {" << rule_of(n) << "} => " << own_variable_of(n) << " ### ###");
/*MSG_QUEUE_FLUSH();*/ MSG_QUEUE_FLUSH();
joint_variable_product_type jvp; joint_variable_product_type jvp;
variable_index_type spawnling = own_variable_of(n);
var_vec varset = variables_of(n); var_vec varset = variables_of(n);
std::sort(varset.begin(), varset.end()); std::sort(varset.begin(), varset.end());
varset = varset - var_vec{spawnling};
/*for (node_index_type o: nei_out(n)) {*/ /*for (node_index_type o: nei_out(n)) {*/
/*varset = varset + variables_of(o);*/ /*varset = varset + variables_of(o);*/
/*}*/ /*}*/
message_type stack; message_type stack;
stack.reserve(64); /* FIXME arbitrary value */ stack.reserve(64); /* FIXME arbitrary value */
for (node_index_type i: inputs) { for (node_index_type i: inputs) {
/*MSG_DEBUG("[jpar_dom] varset " << varset);*/ MSG_DEBUG("[jpar_dom] varset " << varset);
auto idom = get_node_domain(i, varset); auto idom = get_node_domain(i, varset);
/*MSG_DEBUG("[jpar_dom] from domain " << idom);*/ MSG_DEBUG("[jpar_dom] from domain " << idom);
for (const auto& t: idom) { for (const auto& t: idom) {
stack.emplace_back(t); stack.emplace_back(t);
} }
...@@ -1055,16 +1059,21 @@ struct factor_graph_type : public recursive_graph_type<factor_graph_type> { ...@@ -1055,16 +1059,21 @@ struct factor_graph_type : public recursive_graph_type<factor_graph_type> {
genotype_comb_type jpar_dom; genotype_comb_type jpar_dom;
/*MSG_QUEUE_FLUSH();*/ /*MSG_QUEUE_FLUSH();*/
if (stack.size() > 1) { if (stack.size() > 1) {
/*MSG_DEBUG("Domains " << (parent() ? MESSAGE(parent() << "->" << index_in_parent()) : std::string("top-level")) << ' ' << domains);*/ MSG_DEBUG("Domains " << (parent() ? MESSAGE(parent() << "->" << index_in_parent()) : std::string("top-level")) << ' ' << domains);
jpar_dom = compute_product(stack.begin(), stack.end(), varset, domains); jpar_dom = compute_product(stack.begin(), stack.end(), varset, domains);
} else { } else {
jpar_dom = stack.front(); jpar_dom = stack.front();
} }
/*MSG_DEBUG("[jpar_dom] stack " << stack);*/ MSG_DEBUG("[jpar_dom] stack " << stack);
/*MSG_DEBUG("[jpar_dom] result " << jpar_dom);*/ MSG_DEBUG("[jpar_dom] result " << jpar_dom);
variable_index_type spawnling = own_variable_of(n); node_domains[n].emplace_back(generate_factor(rule_of(n), is_dh[n], spawnling, jpar_dom, domains, keep_gametes));
node_domains[n].emplace_back(generate_factor(rule_of(n), is_dh[n], spawnling, jpar_dom, domains)); MSG_DEBUG("keep_gametes " << std::boolalpha << keep_gametes);
if (keep_gametes) {
gametes.emplace_back(-spawnling);
}
/* Add domain for spawnling */ /* Add domain for spawnling */
MSG_DEBUG("node_domains[n] " << node_domains[n]);
MSG_DEBUG("spawnling vec " << var_vec{spawnling});
auto msg = (node_domains[n] % var_vec{spawnling}); auto msg = (node_domains[n] % var_vec{spawnling});
auto dom = msg.front(); auto dom = msg.front();
for (auto& e: dom) { for (auto& e: dom) {
...@@ -1498,8 +1507,8 @@ struct instance_type { ...@@ -1498,8 +1507,8 @@ struct instance_type {
mp.add(evidence[op.emitter]); mp.add(evidence[op.emitter]);
MSG_DEBUG("internal evidence: " << internal_evidence[op.emitter]); MSG_DEBUG("internal evidence: " << internal_evidence[op.emitter]);
mp.add(internal_evidence[op.emitter]); mp.add(internal_evidence[op.emitter]);
// MSG_DEBUG("table: " << tables[op.emitter]); MSG_DEBUG("table: " << tables[op.emitter]);
/*MSG_DEBUG("incoming domains: " << op.incoming_domains);*/ MSG_DEBUG("incoming domains: " << op.incoming_domains);
MSG_DEBUG("Add emitter"); MSG_DEBUG("Add emitter");
mp.add(tables[op.emitter]); mp.add(tables[op.emitter]);
auto inci = op.incoming_messages.begin(), incj = op.incoming_messages.end(); auto inci = op.incoming_messages.begin(), incj = op.incoming_messages.end();
...@@ -1535,12 +1544,12 @@ struct instance_type { ...@@ -1535,12 +1544,12 @@ struct instance_type {
} }
for (size_t o: variant.operations) { for (size_t o: variant.operations) {
const auto& op = operations[o]; const auto& op = operations[o];
/*scoped_indent _(MESSAGE("[" << get_message_index(op.emitter, op.receiver) << "] "));*/ scoped_indent _(MESSAGE("[" << get_message_index(op.emitter, op.receiver) << "] "));
messages[get_message_index(op.emitter, op.receiver)] = compute_message(op, domains); messages[get_message_index(op.emitter, op.receiver)] = compute_message(op, domains);
/*MSG_DEBUG("MESSAGES");*/ MSG_DEBUG("MESSAGES");
/*for (const auto& kv: message_index) {*/ for (const auto& kv: message_index) {
/*MSG_DEBUG(std::setw(4) << kv.first.first << " -> " << std::setw(4) << kv.first.second << " ==" << kv.second << "== " << messages[kv.second]);*/ MSG_DEBUG(std::setw(4) << kv.first.first << " -> " << std::setw(4) << kv.first.second << " ==" << kv.second << "== " << messages[kv.second]);
/*}*/ }
} }
message_type ret; message_type ret;
for (const auto& op: variant.output_operations) { for (const auto& op: variant.output_operations) {
......
...@@ -712,6 +712,114 @@ struct rw_base : public rw_any<rw_base> { ...@@ -712,6 +712,114 @@ struct rw_base : public rw_any<rw_base> {
/** **/ /** **/
struct gamete_LV_database {
/* IND.ped / MARK => 1 or 2 gametes (.second will be empty when cross type is DH) */
std::map<int, std::map<std::string, std::pair<Eigen::Vector2d, Eigen::Vector2d>>> data;
std::map<double, Eigen::Matrix2d> cache;
gamete_LV_database() : data() {}
Eigen::Matrix2d
get_TR(double dist)
{
auto it = cache.find(dist);
if (it != cache.end()) {
double s = .5 + exp(-dist * .02) * .5;
double r = 1. - s;
Eigen::Matrix2d ret;
ret << s, r, r, s;
cache[dist] = ret;
return ret;
}
return it->second;
}
void
add_gametes(const std::string& mark, int ind, const VectorXd& lv, bool dh)
{
std::pair<Eigen::Vector2d, Eigen::Vector2d> g;
if (!dh) {
g.first(0) = lv(0) + lv(1);
g.first(1) = lv(2) + lv(3);
g.second(0) = lv(0) + lv(2);
g.second(1) = lv(1) + lv(3);
} else {
g.first = lv;
g.second(0) = g.second(1) = 0;
}
data[ind][mark] = g;
}
double
map_likelihood(const std::vector<std::string>& marker_names, const std::vector<double>& distances)
{
if (marker_names.size() != distances.size() + 1) {
MSG_ERROR("Need to have exactly one more marker name than the number of inter-marker distances", "");
return NAN;
}
// implicitly, we just garanteed there is at least one marker name in the list
Eigen::Vector2d all_accum;
all_accum(0) = all_accum(1) = 0;
std::vector<Eigen::Matrix2d> TR;
TR.reserve(distances.size());
for (double d: distances) {
TR.push_back(get_TR(d));
}
double ret = 0;
for (auto& kv: data) {
int ind = kv.first;
auto& LV = kv.second;
auto name_i = marker_names.begin();
auto& gam0 = LV[*name_i++];
Eigen::Vector2d accum;
accum << 1, 1;
std::vector<double> norms;
double norm_accum = 0;
norms.reserve(distances.size());