Commit 033edefa authored by Damien Leroux's avatar Damien Leroux
Browse files

Spell-marker works, now needs GMP.

parent 4a6d29c5
This diff is collapsed.
......@@ -93,7 +93,7 @@ add_executable(spell-pedigree ${SPELL_PEDIGREE_SRC})
add_executable(spell-marker ${SPELL_MARKER_SRC})
add_executable(spell-qtl ${SPELL_QTL_SRC})
target_link_libraries(spell-marker dl)
target_link_libraries(spell-marker dl gmp)
target_link_libraries(spell-qtl expat dl)
set(CMAKE_EXE_LINKER_FLAGS "-rdynamic")
......
......@@ -24,7 +24,7 @@
#include <cstring>
#include <unordered_map>
#include <boost/dynamic_bitset.hpp>
#include <gmp.h>
/* La spec est dans le code. */
......@@ -369,6 +369,7 @@ struct table_descr {
std::vector<state_index_type> offset_to_cursor;
uint64_t offset;
state_index_type mask;
bool unif;
inline
const genotype_comb_type::element_type&
......@@ -639,11 +640,14 @@ struct joint_variable_product_type {
tab.offset = 0;
coordinates.resize(tab.variable_names.size());
tab.offset_to_cursor.reserve(tab.data->size());
tab.unif = true;
double coef = tab.data->begin()->coef;
for (uint64_t i = 0; i < tab.data->size(); ++i) {
get_coords_at(tab, i, coordinates);
state_index_type cursor = compute_local_index(coordinates, tab);
tab.offset_to_cursor.push_back(cursor);
tab.domain_to_offset[cursor] = i;
tab.unif = tab.unif && std::abs(coef - tab.data->m_combination[i].coef) > _EPSILON;
}
if (0) {
std::stringstream ss;
......@@ -881,33 +885,40 @@ struct joint_variable_product_type {
#endif
inline
double
compute_coef()
void
compute_coef(mpf_t& accum)
{
auto ti = tables.begin();
auto tj = tables.end();
double accum = ti->current_element().coef;
mpf_t coef;
mpf_init_set_d(accum, ti->current_element().coef);
mpf_init(coef);
for (++ti; ti != tj; ++ti) {
accum *= ti->current_element().coef;
if (!ti->unif) {
mpf_set_d(coef, ti->current_element().coef);
mpf_mul(accum, accum, coef);
}
}
return accum;
}
inline
state_index_type
next_state(state_index_type from)
{
uint64_t valid = 0, all = (1 << tables.size()) - 1;
// uint64_t valid = 0, all = (1 << tables.size()) - 1;
boost::dynamic_bitset<> valid;
valid.resize(tables.size(), 0);
size_t ti = 0;
while (valid != all && !from.is_bad()) {
/*MSG_DEBUG("valid " << valid << " all " << all << " table " << (*tables[ti].data));*/
while (!valid.all() && !from.is_bad()) {
// MSG_DEBUG("valid " << valid << " all " << all << " table " << (*tables[ti].data));
if (!move_cursor_if_state_is_valid(tables[ti], from)) {
auto range = find_offsets_around_cursor(tables[ti], tables[ti].local(from));
/*MSG_DEBUG("offsets are " << range.first << " and " << range.second);*/
// MSG_DEBUG("offsets are " << range.first << " and " << range.second);
state_index_type
c1 = merge_cursors(tables[ti], from, tables[ti].offset_to_cursor[range.first]),
c2 = merge_cursors(tables[ti], from, tables[ti].offset_to_cursor[range.second]);
/*MSG_DEBUG("from = " << dump(from) << " c1 = " << dump(c1, tables[ti]) << " c2 = " << dump(c2, tables[ti]));*/
// MSG_DEBUG("from = " << dump(from) << " c1 = " << dump(c1, tables[ti]) << " c2 = " << dump(c2, tables[ti]));
state_index_type tmp; // = std::min(c1, c2);
if (c1 > from) {
if (c2 > from) {
......@@ -920,15 +931,16 @@ struct joint_variable_product_type {
} else {
return state_index_type::bad();
}
/*MSG_DEBUG("new lower bound is now " << dump(tmp));*/
// MSG_DEBUG("new lower bound is now " << dump(tmp));
if (tmp < from) {
return state_index_type::bad();
}
from = tmp;
valid = 0;
valid.reset();
// valid = 0;
/*ti = !ti;*/
} else {
valid |= 1 << ti;
valid.set(ti);
}
++ti;
ti %= tables.size();
......@@ -964,6 +976,8 @@ struct joint_variable_product_type {
genotype_comb_type
compute_generic()
{
if (invalid_product) { return {}; }
if (total_bits >= 32) {
std::vector<size_t> output_variable_indices;
for (size_t i = 0; i < output_variables_in_use.size(); ++i) {
......@@ -1010,7 +1024,6 @@ struct joint_variable_product_type {
}
key_computing_variant key_update;
if (invalid_product) { return {}; }
genotype_comb_type ret;
size_t counter = 0;
......@@ -1019,27 +1032,39 @@ struct joint_variable_product_type {
/*MSG_DEBUG(key.size()); MSG_QUEUE_FLUSH();*/
state_index_type z;
z.reset();
state_index_type cursor = next_state(z), last = cursor;
std::unordered_map<genotype_comb_type::key_list, double> values;
double norm = 0;
MSG_DEBUG("Initializing cursor");
state_index_type
cursor = next_state(z),
last = cursor;
std::unordered_map<genotype_comb_type::key_list, mpf_t> values;
mpf_t norm;
mpf_init(norm);
if (cursor.is_bad()) {
return ret;
}
for (;;)
{
mpf_t accum;
mpf_init(accum);
for (;;) {
/*if (++counter > 10) {*/
/*break;*/
/*}*/
/*jvp.compute_coordinates(cursor, coords);*/
/*MSG_DEBUG(std::bitset<7>(cursor) << ' ' << coords[0] << ' ' << coords[1] << ' ' << coords[2]);*/
/*MSG_DEBUG("COMPUTE @" << dump(cursor));*/
// MSG_DEBUG("COMPUTE @" << dump(cursor));
key_update(*this, key);
double coef = compute_coef();
if (coef != 0) {
norm += coef;
values[key] += coef;
mpf_clear(accum);
compute_coef(accum);
if (mpf_cmp_d(accum, 0)) {
mpf_add(norm, norm, accum);
auto it = values.find(key);
mpf_t& value = values[key];
if (it == values.end()) {
mpf_init_set(value, accum);
} else {
mpf_add(value, value, accum);
}
}
++counter;
......@@ -1062,15 +1087,17 @@ struct joint_variable_product_type {
}
}
ret.m_combination.reserve(values.size());
if (norm > 0) {
norm = 1. / norm;
if (mpf_cmp_d(norm, 0) > 0) {
mpf_ui_div(norm, 1, norm);
// norm = 1. / norm;
}
for (const auto& kv: values) {
ret.m_combination.emplace_back(kv.first, kv.second * norm);
mpf_mul(accum, kv.second, norm);
ret.m_combination.emplace_back(kv.first, mpf_get_d(accum));
/*ret.m_combination.emplace_back(kv.first, kv.second);*/
}
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");
return ret;
}
......
......@@ -1494,18 +1494,20 @@ struct instance_type {
return sub_instances[op.emitter]->compute(op.receiver, domains);
}
multiple_product_type mp;
/*MSG_DEBUG("evidence: " << evidence[op.emitter]);*/
MSG_DEBUG("evidence: " << 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]);
/*MSG_DEBUG("table: " << tables[op.emitter]);*/
// MSG_DEBUG("table: " << tables[op.emitter]);
/*MSG_DEBUG("incoming domains: " << op.incoming_domains);*/
MSG_DEBUG("Add emitter");
mp.add(tables[op.emitter]);
auto inci = op.incoming_messages.begin(), incj = op.incoming_messages.end();
auto domi = op.incoming_domains.begin();
std::vector<message_type> msgs;
msgs.reserve(op.incoming_messages.size());
for (; inci != incj; ++inci, ++domi) {
MSG_DEBUG("Add incoming");
msgs.emplace_back(messages[*inci] % *domi);
/*MSG_DEBUG("input #" << (*inci) << ": raw " << messages[*inci]);*/
/*MSG_DEBUG("input #" << (*inci) << ": dom " << (*domi));*/
......@@ -1513,6 +1515,7 @@ struct instance_type {
/*mp.add(msgs.back());*/
mp.add(messages[*inci]);
}
MSG_DEBUG("Compute.");
auto ret = mp.compute(op.output, domains);
MSG_DEBUG("" << ret);
MSG_DEBUG("");
......
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