Commit 4211ab81 authored by Damien Leroux's avatar Damien Leroux
Browse files

Now able to compute factors once again.

parent 22f21721
......@@ -26,6 +26,100 @@ operator - (const std::vector<V>& u, const std::vector<V>& v)
struct compute_labels {
bn_label_type
find_label(size_t n, const genotype_comb_type::element_type& labels)
{
auto it
= std::find_if(labels.keys.begin(), labels.keys.end(),
[=] (const genotype_comb_type::key_type& k) { return k.parent == n; }
);
if (it == labels.keys.end()) {
MSG_ERROR("COULDN'T FIND LABEL FOR " << n << " IN " << labels, "");
MSG_QUEUE_FLUSH();
return {};
}
return it->state;
}
bn_label_type
operator () (const pedigree_tree_type& tree, size_t n, const genotype_comb_type::element_type& labels, const std::vector<bool>& recompute)
{
if (tree.get_p2(n) == NONE) {
/* gamete or ancestor */
if (tree.get_p1(n) == NONE) {
/* ancestor */
return find_label(n, labels);
} else {
auto gl = find_label(n, labels);
auto sub = operator () (tree, tree.get_p1(n), labels, recompute);
if (gl.first == GAMETE_L) {
return {sub.first, 0, sub.first_allele, 0};
} else {
return {sub.second, 0, sub.second_allele, 0};
}
}
} else if (recompute[n]) {
auto subl = operator () (tree, tree.get_p1(n), labels, recompute);
auto subr = operator () (tree, tree.get_p2(n), labels, recompute);
return {subl.first, subr.first, subl.first_allele, subr.first_allele};
} else {
return find_label(n, labels);
}
}
std::vector<bn_label_type>
operator () (const pedigree_tree_type& tree, size_t n, const genotype_comb_type& comb, const std::vector<bool>& recompute)
{
std::vector<bn_label_type> ret;
ret.reserve(comb.m_combination.size());
for (const auto& e: comb) {
ret.emplace_back(operator () (tree, n, e, recompute));
}
return ret;
}
};
#if 0
struct factor_set_type {
size_t n_alleles;
std::map<size_t, genotype_comb_type>& joint_prob_tables;
std::map<size_t, std::vector<bn_label_type>>& variable_domains;
void
compute_factor(const graph_type& G, size_t id, size_t p1, size_t p2)
{
if (p1 == 0 && p2 == 0) {
/* ancestor */
std::vector<bn_label_type> labels;
char letter = ped.ancestor_letters.find(ind_node)->second;
for (size_t i = 0; i < n_alleles; ++i) {
labels.emplace_back(letter, letter, i, i);
}
variable_domains[ind_node] = labels;
} else if (p1 == p2) {
/* selfing */
} else if (p2 == 0) {
/* doubled haploid */
} else {
/* standard cross */
}
}
};
#endif
struct bayesian_node_type;
typedef std::shared_ptr<bayesian_node_type> bayesian_node;
struct bn_interface_type;
......
This diff is collapsed.
......@@ -62,7 +62,7 @@ struct bn_label_type {
}
};
typedef combination_type<size_t, bn_label_type> genotype_comb_type;
typedef combination_type<int, bn_label_type> genotype_comb_type;
template <typename Arg>
......@@ -1133,7 +1133,7 @@ struct pedigree_type {
for (size_t i = 0; i < m.labels.size(); ++i) {
for (char a1 = 0; a1 < (char) n_alleles; ++a1) {
for (char a2 = 0; a2 < (char) n_alleles; ++a2) {
ret[n++].m_combination.emplace_back(genotype_comb_type::key_type{node, {m.labels[i].first(), m.labels[i].second(), a1, a2}}, norm);
ret[n++].m_combination.emplace_back(genotype_comb_type::key_type{(int) node, {m.labels[i].first(), m.labels[i].second(), a1, a2}}, norm);
}
}
}
......@@ -1144,7 +1144,7 @@ struct pedigree_type {
double norm = 1.; // / n_alleles;
for (size_t i = 0; i < m.labels.size(); ++i) {
for (char a1 = 0; a1 < (char) n_alleles; ++a1) {
ret[n++].m_combination.emplace_back(genotype_comb_type::key_type{node, {m.labels[i].first(), m.labels[i].second(), a1, 0}}, norm);
ret[n++].m_combination.emplace_back(genotype_comb_type::key_type{(int) node, {m.labels[i].first(), m.labels[i].second(), a1, 0}}, norm);
}
}
return ret;
......
This diff is collapsed.
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