Commit 27ec4c15 authored by damien's avatar damien
Browse files

Generalized product seems to be working now to compute factors with and without gametes.

parent f3153233
......@@ -134,7 +134,7 @@ struct state_index_type {
std::ostream&
operator << (std::ostream& os, const state_index_type& s)
{
MSG_DEBUG("printing state_index " << s.data);
// MSG_DEBUG("printing state_index " << s.data);
for (int i = n_blocks - 1; i >= 0; --i) {
uint64_t mask = 1ULL << 63;
while (mask && !(s.data[i] & mask)) { mask >>= 1; }
......@@ -694,7 +694,7 @@ struct joint_variable_product_type {
}
}
}
if (0) {
if (1) {
std::stringstream ss;
ss << "key extractors";
for (const auto& t_ki: key_extractors) {
......@@ -861,6 +861,7 @@ struct joint_variable_product_type {
auto ti = jvp.tables.begin();
bn_label_type G = ti->current_element().keys.keys.front().state;
MSG_DEBUG("Have gamete label " << G); MSG_QUEUE_FLUSH();
MSG_DEBUG("Have output variable names " << jvp.output_variable_names);
for (const auto& t_ki: jvp.key_extractors) {
*ki++ = jvp.tables[t_ki.first].current_element().keys.keys[t_ki.second];
/**ki++ = ti->current_element().keys.keys.front();*/
......@@ -871,6 +872,7 @@ struct joint_variable_product_type {
bool f1 = G.first == GAMETE_L;
bool f2 = G.second == GAMETE_L;
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,
f2 ? p2.state.first : p2.state.second,
......@@ -881,6 +883,37 @@ struct joint_variable_product_type {
}
};
struct update_key_build_factor_keep_gametes {
size_t
key_size(const joint_variable_product_type& jvp) const { return jvp.output_variable_names.size() + 1; } /* parents + gameteS, where gameteS will be replaced by the output variable */
void
operator () (joint_variable_product_type& jvp, genotype_comb_type::key_list& key) const
{
auto ki = key.begin();
auto ti = jvp.tables.begin();
bn_label_type G = ti->current_element().keys.keys.front().state;
MSG_DEBUG("Have gamete label " << G); MSG_QUEUE_FLUSH();
MSG_DEBUG("Have output variable names " << jvp.output_variable_names);
for (const auto& t_ki: jvp.key_extractors) {
*ki++ = jvp.tables[t_ki.first].current_element().keys.keys[t_ki.second];
/**ki++ = ti->current_element().keys.keys.front();*/
}
MSG_DEBUG("Have temp key " << key); MSG_QUEUE_FLUSH();
const auto& p1 = key.keys[1 + G.first_allele];
const auto& p2 = key.keys[1 + G.second_allele];
bool f1 = G.first == GAMETE_L;
bool f2 = G.second == GAMETE_L;
MSG_DEBUG("Have p1 " << p1 << " p2 " << p2 << " f1 " << f1 << " f2 " << f2); MSG_QUEUE_FLUSH();
*ki = {
-jvp.output_variable_names.front(), {f1 ? p1.state.first : p1.state.second,
f2 ? p2.state.first : p2.state.second,
f1 ? p1.state.first_allele : p1.state.second_allele,
f2 ? p2.state.first_allele : p2.state.second_allele}
};
MSG_DEBUG("Have complete key " << key); MSG_QUEUE_FLUSH();
}
};
#if 0
inline
void
......@@ -1115,6 +1148,7 @@ struct joint_variable_product_type {
inline genotype_comb_type compute() { return compute_generic<update_key>(); }
inline genotype_comb_type build_factor() { return compute_generic<update_key_build_factor>(); }
inline genotype_comb_type build_factor_with_gametes() { return compute_generic<update_key_build_factor_keep_gametes>(); }
/* Algorithm adapted from wikipedia page https://en.wikipedia.org/wiki/Stoer%E2%80%93Wagner_algorithm */
std::pair<int, std::vector<int>>
......@@ -1233,15 +1267,17 @@ generate_factor(const std::vector<int>& parents, bool dh, int spawnling, const g
}
int gamete = -spawnling;
genotype_comb_type G, S;
joint_variable_product_type jvp;
genotype_comb_type
if (keep_gametes) {
G = state_to_combination(gamete, label_g) * (1. / label_g.size());
MSG_DEBUG("before adding gamete domain " << domains);
domains[{-spawnling}] = G;
MSG_DEBUG("after adding gamete domain " << domains);
joint_variable_product_type jvp;
domains[{gamete}] = G;
} else {
G = state_to_combination(spawnling, label_g) * (1. / label_g.size());
domains[{spawnling}] = G;
}
MSG_DEBUG("joint parent domain " << parent_domain);
MSG_DEBUG("all domains " << domains);
jvp.add_table(G);
......@@ -1252,7 +1288,7 @@ generate_factor(const std::vector<int>& parents, bool dh, int spawnling, const g
jvp.set_output(sorted_parents + std::vector<int>{spawnling});
}
jvp.compile(domains);
return jvp.build_factor();
return keep_gametes ? jvp.build_factor_with_gametes() : jvp.build_factor();
}
inline
......
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