Commit 91141d66 authored by Damien Leroux's avatar Damien Leroux
Browse files

Integration in progress. About to trim old code from graphnode.

parent 10e3be25
......@@ -226,6 +226,13 @@ namespace std {
/* Permet d'itérer joliment sur la séquence inverse (rbegin->rend) avec la syntaxe for (v: container) {} */
template <typename CONTAINER> struct const_reversed_impl {
const CONTAINER& ref;
typedef typename std::decay<CONTAINER>::type::const_reverse_iterator iterator;
iterator begin() { return ref.rbegin(); }
iterator end() { return ref.rend(); }
};
template <typename CONTAINER> struct reversed_impl {
const CONTAINER& ref;
typedef typename std::decay<CONTAINER>::type::reverse_iterator iterator;
......@@ -234,6 +241,7 @@ template <typename CONTAINER> struct reversed_impl {
};
template <typename CONTAINER> reversed_impl<CONTAINER> reversed(CONTAINER&& x) { return {x}; }
template <typename CONTAINER> const_reversed_impl<CONTAINER> reversed(const CONTAINER& x) { return {x}; }
/* Descripteur du domaine d'une variable participant au produit */
......@@ -493,7 +501,7 @@ struct joint_variable_product_type {
tables.emplace_back();
tables.back().data = &t;
tables.back().variable_names.swap(vars);
MSG_DEBUG("[joint_product] Added table " << (&t) << " " << t);
/*MSG_DEBUG("[joint_product] Added table " << (&t) << " " << t);*/
}
inline
......@@ -519,8 +527,8 @@ struct joint_variable_product_type {
d.bit_width = 1;
while ((1ULL << d.bit_width) < d.size) { ++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_QUEUE_FLUSH();
/*MSG_DEBUG("Have domain for variable #" << d.variable_name << " size " << d.size << " bit_width " << d.bit_width);*/
/*MSG_QUEUE_FLUSH();*/
d.index_to_label.reserve(domain.size());
for (size_t i = 0; i < domain.size(); ++i) {
const auto& label = domain.m_combination[i].keys.keys.front().state;
......@@ -537,8 +545,9 @@ struct joint_variable_product_type {
/*d.shift = shift;*/
shift += d.bit_width;
d.assign(max_cursor, 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();*/
/*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();*/
if (0) {
std::stringstream ss;
for (const auto& kv: d.label_to_index) {
......@@ -570,10 +579,11 @@ struct joint_variable_product_type {
for (size_t v = 0; v < tab.variable_indices[0]; ++v) {
tab.variables_left[v] = true;
}
MSG_DEBUG("Table " << tab.data);
MSG_QUEUE_FLUSH();
MSG_DEBUG("Table variables " << tab.variable_names << " mask " << tab.mask);
MSG_QUEUE_FLUSH();
/*MSG_QUEUE_FLUSH();*/
/*MSG_DEBUG("Table " << tab.data);*/
/*MSG_QUEUE_FLUSH();*/
/*MSG_DEBUG("Table variables " << tab.variable_names << " mask " << tab.mask);*/
/*MSG_QUEUE_FLUSH();*/
tab.offset = 0;
coordinates.resize(tab.variable_names.size());
tab.offset_to_cursor.reserve(tab.data->size());
......@@ -583,7 +593,7 @@ struct joint_variable_product_type {
tab.offset_to_cursor.push_back(cursor);
tab.domain_to_offset[cursor] = i;
}
if (1) {
if (0) {
std::stringstream ss;
ss << "Table";
for (const auto& x: tab.domain_to_offset) {
......@@ -1000,8 +1010,8 @@ compute_product(Iterator tables_begin, Iterator tables_end, const std::vector<in
for (; tables_begin != tables_end; ++tables_begin) {
if (getter(*tables_begin).size()) {
jvp.add_table(getter(*tables_begin));
} else {
MSG_DEBUG("[compute_product] Didn't add table " << (*tables_begin) << " {" << getter(*tables_begin) << "} because it seems to be empty.");
/*} else {*/
/*MSG_DEBUG("[compute_product] Didn't add table " << (*tables_begin) << " {" << getter(*tables_begin) << "} because it seems to be empty.");*/
}
}
jvp.set_output(output_variables);
......
This diff is collapsed.
......@@ -6,6 +6,7 @@
#include <map>
#include <set>
#include <unordered_set>
#include <unordered_map>
#include <string>
#include <iostream>
#include "file.h"
......@@ -427,8 +428,6 @@ geno_matrix* read_geno_matrix(ifile& ifs)
template <class DERIVED>
struct rw_any {
bool fourcc(ifile& ifs, const char* cc)
......@@ -472,7 +471,21 @@ struct rw_any {
void operator () (ofile& ofs, ptrdiff_t i) { write_ptrdiff(ofs, i); }
void operator () (ifile& ifs, label_type& l) { l = {read_char(ifs), read_char(ifs)}; }
void operator () (ofile& ofs, const label_type& l) { write_char(ofs, l.first()); write_char(ofs, l.second()); }
void operator () (ofile& ofs, label_type& l) { write_char(ofs, l.first()); write_char(ofs, l.second()); }
template <typename K>
auto operator () (ifile& fs, K& has_file_io)
-> decltype(has_file_io.file_io(fs), void())
{
has_file_io.file_io(fs);
}
template <typename K>
auto operator () (ofile& fs, K& has_file_io)
-> decltype(has_file_io.file_io(fs), void())
{
has_file_io.file_io(fs);
}
template <typename V>
void operator () (ifile& ifs, std::set<V>& vec)
......@@ -488,11 +501,11 @@ struct rw_any {
}
template <typename V>
void operator () (ofile& ofs, const std::set<V>& vec)
void operator () (ofile& ofs, std::set<V>& vec)
{
if (fourcc(ofs, "OSET")) { return; }
write_size(ofs, vec.size());
for (const auto& e: vec) {
for (auto& e: vec) {
ref() (ofs, e);
}
}
......@@ -512,11 +525,11 @@ struct rw_any {
}
template <typename V>
void operator () (ofile& ofs, const std::unordered_set<V>& vec)
void operator () (ofile& ofs, std::unordered_set<V>& vec)
{
if (fourcc(ofs, "USET")) { return; }
write_size(ofs, vec.size());
for (const auto& e: vec) {
for (auto& e: vec) {
ref() (ofs, e);
}
}
......@@ -537,7 +550,16 @@ struct rw_any {
}
template <typename V, typename A>
void operator () (ofile& ofs, const std::vector<V, A>& vec)
void operator () (ofile& ofs, std::vector<V, A>& vec)
{
if (fourcc(ofs, "VECT")) { return; }
write_size(ofs, vec.size());
for (auto& e: vec) {
ref() (ofs, e);
}
}
void operator () (ofile& ofs, std::vector<bool>& vec)
{
if (fourcc(ofs, "VECT")) { return; }
write_size(ofs, vec.size());
......@@ -585,21 +607,54 @@ struct rw_any {
}
template <typename K, typename V, typename A, typename C>
void operator () (ofile& ofs, const std::map<K, V, A, C>& map)
void operator () (ifile& ifs, std::unordered_map<K, V, A, C>& map)
{
if (fourcc(ifs, "MAP ")) { return; }
size_t count = read_size(ifs);
map.clear();
for (size_t i = 0; i < count; ++i) {
K key;
V value;
ref() (ifs, key);
ref() (ifs, value);
map.emplace(std::move(key), std::move(value));
}
}
template <typename K, typename V, typename A, typename C>
void operator () (ofile& ofs, std::map<K, V, A, C>& map)
{
if (fourcc(ofs, "MAP ")) { return; }
write_size(ofs, map.size());
for (const auto& kv: map) {
for (auto& kv: map) {
ref() (ofs, const_cast<K&>(kv.first));
ref() (ofs, kv.second);
}
}
template <typename K, typename V, typename A, typename C>
void operator () (ofile& ofs, std::unordered_map<K, V, A, C>& map)
{
if (fourcc(ofs, "MAP ")) { return; }
write_size(ofs, map.size());
for (auto& kv: map) {
ref() (ofs, kv.first);
ref() (ofs, kv.second);
}
}
template <typename IO, typename A, typename B>
void operator () (IO& fs, std::pair<A, B>& p)
{
ref() (fs, p.first);
ref() (fs, p.second);
}
template <typename SCALAR, int ROW, int COL, int C, int D, int E>
void operator () (ifile& ifs, Eigen::Matrix<SCALAR, ROW, COL, C, D, E>& mat) { read_matrix(ifs, mat); }
template <typename SCALAR, int ROW, int COL, int C, int D, int E>
void operator () (ofile& ofs, const Eigen::Matrix<SCALAR, ROW, COL, C, D, E>& mat) { write_matrix(ofs, mat); }
void operator () (ofile& ofs, Eigen::Matrix<SCALAR, ROW, COL, C, D, E>& mat) { write_matrix(ofs, mat); }
void operator () (ifile& ifs, geno_matrix& mat) { read_geno_matrix(ifs, mat); }
void operator () (ofile& ofs, geno_matrix& mat) { write_geno_matrix(ofs, mat); }
......@@ -614,7 +669,7 @@ struct rw_any {
}
}
void operator () (ofile& ofs, const std::shared_ptr<geno_matrix>& ptr)
void operator () (ofile& ofs, std::shared_ptr<geno_matrix>& ptr)
{
if (ptr) {
ref() (ofs, *ptr.get());
......
......@@ -144,6 +144,8 @@ struct ostream_manager {
std::clog.rdbuf(&forbid);
std::cout.rdbuf(&forbid);
std::cerr.rdbuf(&forbid);
#else
std::clog.rdbuf(&hi);
#endif
}
......@@ -393,6 +395,16 @@ void message_queue::run()
}
}
#ifdef SPELL_UNSAFE_OUTPUT
#define MSG_ERROR(_msg_expr_, _workaround_expr_) do { std::cerr << _msg_expr_ << std::endl; } while (0)
#define MSG_WARNING(_msg_expr_) do { std::cout << _msg_expr_ << std::endl; } while (0)
#define MSG_INFO(_msg_expr_) do { std::cout << _msg_expr_ << std::endl; } while (0)
#define MSG_QUEUE_FLUSH()
#define MSG_DEBUG(_msg_expr_) do { std::clog << _msg_expr_ << std::endl; } while (0)
#else
#define MSG_ERROR(_msg_expr_, _workaround_expr_) \
do {\
CREATE_MESSAGE(msg_channel::Err, MESSAGE(msg_handler_t::e() << "[ERR] " << _msg_expr_ << msg_handler_t::n() << std::endl));\
......@@ -410,6 +422,11 @@ void message_queue::run()
CREATE_MESSAGE(msg_channel::Out, MESSAGE(msg_handler_t::i() << "[MSG] " << _msg_expr_ << msg_handler_t::n() << std::endl));\
} while (0)
#define MSG_QUEUE_FLUSH() do { \
CREATE_MESSAGE(msg_channel::Log, _DBG_SYNC_MARK_S); \
msg_handler_t::wait_for_flush(); \
} while (0)
#define MSG_DEBUG(_msg_expr_) \
do {\
if (msg_handler_t::debug_enabled()) {\
......@@ -417,10 +434,7 @@ void message_queue::run()
}\
} while (0);
#define MSG_QUEUE_FLUSH() do { \
CREATE_MESSAGE(msg_channel::Log, _DBG_SYNC_MARK_S); \
msg_handler_t::wait_for_flush(); \
} while (0)
#endif
#if 0
#define MSG_ERROR(_msg_expr_, _workaround_expr_) \
......
......@@ -297,16 +297,16 @@ struct rw_comb : public rw_any<rw_comb<PARENT_TYPE, STATE_TYPE>> {
using rw_any<rw_comb<PARENT_TYPE, STATE_TYPE>>::operator ();
void operator () (ifile& ifs, bn_label_type& l) { l.compact() = read_int(ifs); }
void operator () (ofile& ofs, const bn_label_type& l) { write_int(ofs, l.compact()); }
void operator () (ofile& ofs, bn_label_type& l) { write_int(ofs, l.compact()); }
void operator () (ifile& fs, typename comb_type::key_type& key) { ref() (fs, key.parent); ref() (fs, key.state); }
void operator () (ofile& fs, const typename comb_type::key_type& key) { ref() (fs, key.parent); ref() (fs, key.state); }
void operator () (ofile& fs, typename comb_type::key_type& key) { ref() (fs, key.parent); ref() (fs, key.state); }
void operator () (ifile& fs, typename comb_type::key_list& keys) { ref() (fs, keys.keys); }
void operator () (ofile& fs, const typename comb_type::key_list& keys) { ref() (fs, keys.keys); }
void operator () (ofile& fs, typename comb_type::key_list& keys) { ref() (fs, keys.keys); }
void operator () (ifile& fs, typename comb_type::element_type& elem) { ref() (fs, elem.keys); ref() (fs, elem.coef); }
void operator () (ofile& fs, const typename comb_type::element_type& elem) { ref() (fs, elem.keys); ref() (fs, elem.coef); }
void operator () (ofile& fs, typename comb_type::element_type& elem) { ref() (fs, elem.keys); ref() (fs, elem.coef); }
void operator () (ifile& fs, comb_type& comb)
{
......@@ -314,7 +314,7 @@ struct rw_comb : public rw_any<rw_comb<PARENT_TYPE, STATE_TYPE>> {
ref() (fs, comb.m_combination);
}
void operator () (ofile& fs, const comb_type& comb)
void operator () (ofile& fs, comb_type& comb)
{
if (fourcc(fs, "COMB")) { return; }
ref() (fs, comb.m_combination);
......@@ -330,7 +330,7 @@ struct rw_tree : public rw_any<rw_tree> {
using rw_any<rw_tree>::operator ();
void operator () (ifile& fs, pedigree_tree_type::node_type& node) { ref() (fs, node.p1); ref() (fs, node.p2); }
void operator () (ofile& fs, const pedigree_tree_type::node_type& node) { ref() (fs, node.p1); ref() (fs, node.p2); }
void operator () (ofile& fs, pedigree_tree_type::node_type& node) { ref() (fs, node.p1); ref() (fs, node.p2); }
template <typename STREAM_TYPE, typename TREE_TYPE>
void tree_io_impl(STREAM_TYPE& fs, TREE_TYPE&& tree)
......@@ -349,7 +349,7 @@ struct rw_tree : public rw_any<rw_tree> {
tree_io_impl(fs, tree);
}
void operator () (ofile& fs, const pedigree_tree_type& tree)
void operator () (ofile& fs, pedigree_tree_type& tree)
{
tree_io_impl(fs, tree);
}
......@@ -362,7 +362,7 @@ struct rw_tree : public rw_any<rw_tree> {
(*this)(fs, pi.p2);
}
void operator () (ofile& fs, const pedigree_item& pi)
void operator () (ofile& fs, pedigree_item& pi)
{
(*this)(fs, pi.gen_name);
(*this)(fs, pi.id);
......
#define SPELL_UNSAFE_OUTPUT
#include "bayes/generalized_product.h"
#include "graphnode.h"
......@@ -438,6 +439,13 @@ filter_pedigree(const std::vector<pedigree_item>& full_pedigree, const std::vect
inline bool ends_with(std::string const & value, std::string const & ending)
{
if (ending.size() > value.size()) return false;
return std::equal(ending.rbegin(), ending.rend(), value.rbegin());
}
int
main(int argc, char** argv)
......@@ -542,11 +550,11 @@ main(int argc, char** argv)
MSG_DEBUG("F4 * F5 (squeeze F4) size " << factors[100].size());
MSG_DEBUG("F4 * F5 * F6 (squeeze F4&F5) size " << factors[101].size());
#endif
graph_type g;
std::unique_ptr<graph_type> g;
if (argc == 1) {
g = graph_type::from_pedigree(read_csv("random.ped", ';'), 2, {}, {}, "boiboite");
g.to_image("boiboite", "png");
g->to_image("boiboite", "png");
/*g.dump_active();*/
/*g = graph_type::from_pedigree(read_csv("random12.ped", ';'), 2, {}, {}, "boiboite12");*/
......@@ -554,41 +562,69 @@ main(int argc, char** argv)
/*g.dump_active();*/
g = graph_type::from_pedigree(read_csv("/home/daleroux/devel/spel/sample-data/data_magic/pedigree_clean_DH.csv", ';'), 1, {"Cervil", "Levovil", "Criollo", "Stupicke", "Plovdiv", "LA1420", "Ferum", "LA0147", "ICS3"}, {"ICS3"});
g.to_image("magic-dh", "png");
g->to_image("magic-dh", "png");
/*g = graph_type::from_pedigree(read_csv("/home/daleroux/devel/spel/sample-data/data_magic/pedigree_clean-norepeat.csv", ';'), 1, {}, {});*/
/*g.to_image("magic-dh", "png");*/
/*g.dump_active();*/
} else {
std::string pedfile = argv[1];
std::string prefix = argv[2];
std::vector<std::string> in, out;
size_t pos;
{
std::string s(argv[3]);
while ((pos = s.find(",")) != std::string::npos) {
in.push_back(s.substr(0, pos));
s.erase(0, pos + 1);
if (ends_with(pedfile, ".ped")) {
std::string prefix = argv[2];
std::vector<std::string> in, out;
size_t pos;
{
std::string s(argv[3]);
while ((pos = s.find(",")) != std::string::npos) {
in.push_back(s.substr(0, pos));
s.erase(0, pos + 1);
}
in.push_back(s);
}
in.push_back(s);
}
{
std::string s(argv[4]);
while ((pos = s.find(",")) != std::string::npos) {
out.push_back(s.substr(0, pos));
s.erase(0, pos + 1);
{
std::string s(argv[4]);
while ((pos = s.find(",")) != std::string::npos) {
out.push_back(s.substr(0, pos));
s.erase(0, pos + 1);
}
out.push_back(s);
}
MSG_DEBUG("[args] Pedigree: " << pedfile);
MSG_DEBUG("[args] Prefix: " << prefix);
MSG_DEBUG("[args] In: " << in);
MSG_DEBUG("[args] Out: " << out);
/*g = graph_type::from_pedigree(read_csv(pedfile, ';'), 2, in, out);*/
/*g.to_image(MESSAGE(prefix << "-pre-optim"), "png");*/
/*g.optimize();*/
/*g.to_image(prefix, "png");*/
g = graph_type::from_pedigree(filter_pedigree(read_csv(pedfile, ';'), in, out), 2, in, out, MESSAGE("debug-" << prefix));
g->to_image(MESSAGE(prefix << "-filtered-pre-optim"), "png");
/*g->dump_active();*/
g->optimize();
MSG_DEBUG("Creating png...");
g->to_image(MESSAGE(prefix << "-filtered"), "png");
MSG_DEBUG("Creating instance...");
auto I = g->instance();
ofile of(MESSAGE(prefix << "-instance.data"));
I->file_io(of);
rw_comb<int, bn_label_type>()(of, g->domains);
} else {
std::unique_ptr<graph_type::instance_type> I(new graph_type::instance_type());
std::map<var_vec, genotype_comb_type> domains;
ifile ifs(pedfile);
I->file_io(ifs);
rw_comb<int, bn_label_type>()(ifs, domains);
MSG_DEBUG("Computing messages....");
auto ret = I->compute(0, domains);
MSG_DEBUG("MESSAGES");
for (const auto& kv: I->message_index) {
MSG_DEBUG(std::setw(4) << kv.first.first << " -> " << std::setw(4) << kv.first.second << " == " << I->messages[kv.second]);
}
MSG_DEBUG("OUTPUT");
for (const auto& t: ret) {
MSG_DEBUG("" << t);
}
out.push_back(s);
}
MSG_DEBUG("[args] Pedigree: " << pedfile);
MSG_DEBUG("[args] Prefix: " << prefix);
MSG_DEBUG("[args] In: " << in);
MSG_DEBUG("[args] Out: " << out);
g = graph_type::from_pedigree(read_csv(pedfile, ';'), 2, in, out);
g.to_image(prefix, "png");
g = graph_type::from_pedigree(filter_pedigree(read_csv(pedfile, ';'), in, out), 2, in, out);
g.to_image(MESSAGE(prefix << "-filtered"), "png");
g.compute_domains_and_factors();
}
return 0;
......
......@@ -77,6 +77,21 @@ dispatch_geno_probs(
std::vector<bn_label_type>
get_domain(std::map<var_vec, genotype_comb_type>& domains, variable_index_type var)
{
std::vector<bn_label_type> ret;
const auto& dom = domains[{var}];
ret.reserve(dom.size());
for (const auto& e: dom) {
ret.push_back(e.keys.keys.front().state);
}
return ret;
}
std::map<std::string, std::pair<count_jobs_type, job_type>>
job_registry = {
{"dummy-one", {
......@@ -142,8 +157,13 @@ job_registry = {
/*factor_graph fg(settings->pedigree, unique_n_alleles[n], settings->noise);*/
auto fg = graph_type::from_pedigree(settings->pedigree.items, unique_n_alleles[n], settings->input_generations, settings->output_generations);
/*MSG_DEBUG("COMPUTED FACTOR GRAPH" << std::endl << fg);*/
fg.finalize();
fg.save(settings->job_filename("factor-graph", unique_n_alleles[n]));
/*fg.finalize();*/
/*fg.save(settings->job_filename("factor-graph", unique_n_alleles[n]));*/
fg->optimize();
auto I = fg->instance();
ofile of(settings->job_filename("factor-graph", unique_n_alleles[n]));
I->file_io(of);
rw_comb<int, bn_label_type>()(of, fg->domains);
return true;
}
}},
......@@ -163,8 +183,13 @@ job_registry = {
/*MSG_DEBUG("ALLELE CONVERT " << c << " => " << ((int) allele_obs_to_idx[c]));*/
}
}
graph_type fg;
fg.load(settings->job_filename("factor-graph", used_alleles.size()));
/*graph_type fg;*/
/*fg.load(settings->job_filename("factor-graph", used_alleles.size()));*/
ifile ifs(settings->job_filename("factor-graph", used_alleles.size()));
std::unique_ptr<graph_type::instance_type> I(new graph_type::instance_type());
std::map<var_vec, genotype_comb_type> domains;
I->file_io(ifs);
rw_comb<int, bn_label_type>()(ifs, domains);
std::map<char, int> allele_obs2bn;
for (char c: used_alleles) { allele_obs2bn[c] = allele_obs2bn.size(); }
......
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