Commit 15f72f0b authored by Damien Leroux's avatar Damien Leroux
Browse files

Integration in progress. Now computing LVs.

parent 714f8fc8
......@@ -2,7 +2,7 @@
#define _SPEL_MARKER_OBS_BAYES_H_
#include "error.h"
#include "generation_rs_fwd.h"
/*#include "generation_rs_fwd.h"*/
/*#include "bayes/graph.h"*/
#if 0
......@@ -40,7 +40,8 @@ size_t size(const std::vector<size_t>& dim) {
#include "pedigree.h"
#include "bayes/bn.h"
/*#include "bayes/bn.h"*/
#include "bayes/factor_var3.h"
#include "bayes/cli.h"
#include "bayes/dispatch.h"
......
......@@ -4,9 +4,10 @@
#include "error.h"
#include <map>
#include <vector>
#include "generation_rs.h"
/*#include "generation_rs.h"*/
#include "commandline.h"
#include "bn.h"
/*#include "bn.h"*/
#include "factor_var3.h"
enum JobDispatchScheme { JDS_None, JDS_MT, JDS_SGE, JDS_SSH };
......@@ -31,20 +32,21 @@ struct bn_settings_t {
std::string pop_name;
std::string qtl_generation_name;
std::string work_directory;
std::vector<chromosome> map;
std::string map_filename;
std::map<std::string, marker_observation_spec> marker_observation_specs;
std::vector<std::string> marker_observation_specs_filenames;
std::map<std::string, population_marker_observation> observed_mark;
double noise;
double tolerance;
/* Inputs */
pedigree_type pedigree;
std::string pedigree_filename;
std::string design_filename;
/* BN */
pedigree_bayesian_network* bn;
/*pedigree_bayesian_network* bn;*/
/* Management */
std::map<std::string, std::map<char, VectorXd>> compiled_obs_specs;
std::vector<std::map<char, int>> alleles_per_marker;
std::vector<size_t> unique_n_alleles;
/*std::map<std::string, std::map<char, VectorXd>> compiled_obs_specs;*/
std::vector<std::string> marker_names;
std::map<std::string, size_t> marker_index;
......@@ -61,17 +63,16 @@ struct bn_settings_t {
, pop_name()
, qtl_generation_name()
, work_directory(".")
, map()
, map_filename()
, marker_observation_specs()
, marker_observation_specs_filenames()
, observed_mark()
, noise(0.)
, tolerance(1.e-10)
, pedigree()
, pedigree_filename()
, design_filename()
, bn(NULL)
, compiled_obs_specs()
/*, bn(NULL)*/
, alleles_per_marker()
, marker_names()
, marker_index()
......
This diff is collapsed.
......@@ -45,6 +45,22 @@ size_t read_size(std::ifstream& ifs)
return ret;
}
/** PTRDIFF_T **/
inline
void write_ptrdiff(std::ofstream& ofs, ptrdiff_t sz)
{
ofs.write((const char*) &sz, sizeof sz);
}
inline
ptrdiff_t read_ptrdiff(std::ifstream& ifs)
{
ptrdiff_t ret;
ifs.read((char*) &ret, sizeof ret);
return ret;
}
/** STRING **/
inline
......@@ -405,6 +421,215 @@ geno_matrix* read_geno_matrix(std::ifstream& ifs)
}
template <class DERIVED>
struct rw_any {
bool fourcc(std::ifstream& ifs, const char* cc)
{
if (check_fourcc(ifs, cc)) {
MSG_ERROR("File is not valid or has been corrupted (expected " << cc << ')', "");
return true;
}
return false;
}
bool fourcc(std::ofstream& ofs, const char* cc)
{
write_fourcc(ofs, cc);
return false;
}
virtual ~rw_any() {}
DERIVED& ref() { return *dynamic_cast<DERIVED*>(this); }
void operator () (std::ifstream& ifs, std::string& s) { s = read_str(ifs); }
void operator () (std::ofstream& ofs, const std::string& s) { write_str(ofs, s); }
void operator () (std::ifstream& ifs, char& i) { i = read_char(ifs); }
void operator () (std::ofstream& ofs, char i) { write_char(ofs, i); }
void operator () (std::ifstream& ifs, int& i) { i = read_int(ifs); }
void operator () (std::ofstream& ofs, int i) { write_int(ofs, i); }
void operator () (std::ifstream& ifs, bool& i) { i = !!read_char(ifs); }
void operator () (std::ofstream& ofs, bool i) { write_char(ofs, i); }
void operator () (std::ifstream& ifs, double& i) { i = read_double(ifs); }
void operator () (std::ofstream& ofs, double i) { write_double(ofs, i); }
void operator () (std::ifstream& ifs, size_t& i) { i = read_size(ifs); }
void operator () (std::ofstream& ofs, size_t i) { write_size(ofs, i); }
void operator () (std::ifstream& ifs, ptrdiff_t& i) { i = read_ptrdiff(ifs); }
void operator () (std::ofstream& ofs, ptrdiff_t i) { write_ptrdiff(ofs, i); }
void operator () (std::ifstream& ifs, label_type& l) { l = {read_char(ifs), read_char(ifs)}; }
void operator () (std::ofstream& ofs, const label_type& l) { write_char(ofs, l.first()); write_char(ofs, l.second()); }
template <typename V>
void operator () (std::ifstream& ifs, std::set<V>& vec)
{
if (fourcc(ifs, "OSET")) { return; }
size_t sz = read_size(ifs);
vec.clear();
for (size_t i = 0; i < sz; ++i) {
V tmp;
ref() (ifs, tmp);
vec.emplace(tmp);
}
}
template <typename V>
void operator () (std::ofstream& ofs, const std::set<V>& vec)
{
if (fourcc(ofs, "OSET")) { return; }
write_size(ofs, vec.size());
for (const auto& e: vec) {
ref() (ofs, e);
}
}
template <typename V>
void operator () (std::ifstream& ifs, std::unordered_set<V>& vec)
{
if (fourcc(ifs, "USET")) { return; }
size_t sz = read_size(ifs);
vec.clear();
for (size_t i = 0; i < sz; ++i) {
V tmp;
ref() (ifs, tmp);
vec.emplace(tmp);
}
}
template <typename V>
void operator () (std::ofstream& ofs, const std::unordered_set<V>& vec)
{
if (fourcc(ofs, "USET")) { return; }
write_size(ofs, vec.size());
for (const auto& e: vec) {
ref() (ofs, e);
}
}
void operator () (std::ifstream& ifs, std::vector<bool>::reference i) { i = !!read_char(ifs); }
template <typename V, typename A>
void operator () (std::ifstream& ifs, std::vector<V, A>& vec)
{
if (fourcc(ifs, "VECT")) { return; }
size_t sz = read_size(ifs);
vec.clear();
vec.reserve(sz);
for (size_t i = 0; i < sz; ++i) {
vec.emplace_back();
ref() (ifs, vec.back());
}
}
template <typename V, typename A>
void operator () (std::ofstream& ofs, const std::vector<V, A>& vec)
{
if (fourcc(ofs, "VECT")) { return; }
write_size(ofs, vec.size());
for (const auto& e: vec) {
ref() (ofs, e);
}
}
template <typename V, typename A, typename RWElem>
void operator () (std::ifstream& ifs, std::vector<V, A>& vec, RWElem rw_elem)
{
if (fourcc(ifs, "VECT")) { return; }
size_t sz = read_size(ifs);
vec.clear();
vec.reserve(sz);
for (size_t i = 0; i < sz; ++i) {
vec.emplace_back();
rw_elem(ifs, vec.back());
}
}
template <typename V, typename A, typename RWElem>
void operator () (std::ofstream& ofs, std::vector<V, A>& vec, RWElem rw_elem)
{
if (fourcc(ofs, "VECT")) { return; }
write_size(ofs, vec.size());
for (auto& e: vec) {
rw_elem(ofs, e);
}
}
template <typename K, typename V, typename A, typename C>
void operator () (std::ifstream& ifs, std::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 () (std::ofstream& ofs, const std::map<K, V, A, C>& map)
{
if (fourcc(ofs, "MAP ")) { return; }
write_size(ofs, map.size());
for (const auto& kv: map) {
ref() (ofs, kv.first);
ref() (ofs, kv.second);
}
}
template <typename SCALAR, int ROW, int COL, int C, int D, int E>
void operator () (std::ifstream& 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 () (std::ofstream& ofs, const Eigen::Matrix<SCALAR, ROW, COL, C, D, E>& mat) { write_matrix(ofs, mat); }
void operator () (std::ifstream& ifs, geno_matrix& mat) { read_geno_matrix(ifs, mat); }
void operator () (std::ofstream& ofs, geno_matrix& mat) { write_geno_matrix(ofs, mat); }
void operator () (std::ifstream& ifs, std::shared_ptr<geno_matrix>& ptr)
{
ptr.reset();
ptr = std::make_shared<geno_matrix>();
ref() (ifs, *ptr);
if (!ptr->size()) {
ptr.reset();
}
}
void operator () (std::ofstream& ofs, const std::shared_ptr<geno_matrix>& ptr)
{
if (ptr) {
ref() (ofs, *ptr.get());
} else {
geno_matrix _;
ref() (ofs, _);
}
}
};
struct rw_base : public rw_any<rw_base> {
virtual ~rw_base() {}
using rw_any<rw_base>::fourcc;
using rw_any<rw_base>::ref;
using rw_any<rw_base>::operator ();
};
/** **/
......
......@@ -143,6 +143,29 @@ struct combination_type {
return ret;
}
friend
key_list
operator / (const key_list& k1, const std::vector<PARENT_TYPE>& discard)
{
key_list ret;
ret.keys.reserve(k1.keys.size() - discard.size());
auto ki = k1.begin();
auto kj = k1.end();
auto di = discard.begin();
auto dj = discard.end();
for (; ki != kj && di != dj; ++ki) {
if (ki->parent == *di) {
++di;
} else if (ki->parent > *di) {
++di;
} else {
ret.keys.push_back(*ki);
}
}
ret.keys.insert(ret.keys.end(), ki, kj);
return ret;
}
std::pair<key_list, key_list>
extract(const std::vector<PARENT_TYPE>& parents) const
{
......@@ -640,7 +663,7 @@ normalize_over(const combination_type<PARENT_TYPE, STATE_TYPE>& comb, const std:
auto k = e.keys % variables;
norm_coef[k] += e.coef;
}
MSG_DEBUG("NORM FACTORS: " << norm_coef);
/*MSG_DEBUG("NORM FACTORS: " << norm_coef);*/
for (auto& kv: norm_coef) { kv.second = 1. / kv.second; }
for (const auto& e: comb) {
auto k = e.keys % variables;
......@@ -663,7 +686,7 @@ sum_over_dual(const combination_type<PARENT_TYPE, STATE_TYPE>& comb, const std::
ret[ke.second.keys] += {ke.first, ke.second.coef};
}
MSG_DEBUG("sum_over_dual(" << comb << ", " << parents << ") => " << ret);
/*MSG_DEBUG("sum_over_dual(" << comb << ", " << parents << ") => " << ret);*/
return ret;
}
......
......@@ -284,187 +284,6 @@ template <> struct make_one<VectorXd> {
};
template <class DERIVED>
struct rw_any {
bool fourcc(std::ifstream& ifs, const char* cc)
{
if (check_fourcc(ifs, cc)) {
MSG_ERROR("File is not valid or has been corrupted", "");
return true;
}
return false;
}
bool fourcc(std::ofstream& ofs, const char* cc)
{
write_fourcc(ofs, cc);
return false;
}
virtual ~rw_any() {}
DERIVED& ref() { return *dynamic_cast<DERIVED*>(this); }
void operator () (std::ifstream& ifs, std::string& s) { s = read_str(ifs); }
void operator () (std::ofstream& ofs, const std::string& s) { write_str(ofs, s); }
void operator () (std::ifstream& ifs, char& i) { i = read_char(ifs); }
void operator () (std::ofstream& ofs, char i) { write_char(ofs, i); }
void operator () (std::ifstream& ifs, int& i) { i = read_int(ifs); }
void operator () (std::ofstream& ofs, int i) { write_int(ofs, i); }
void operator () (std::ifstream& ifs, bool& i) { i = !!read_char(ifs); }
void operator () (std::ofstream& ofs, bool i) { write_char(ofs, i); }
void operator () (std::ifstream& ifs, double& i) { i = read_double(ifs); }
void operator () (std::ofstream& ofs, double i) { write_double(ofs, i); }
void operator () (std::ifstream& ifs, size_t& i) { i = read_size(ifs); }
void operator () (std::ofstream& ofs, size_t i) { write_size(ofs, i); }
void operator () (std::ifstream& ifs, label_type& l) { l = {read_char(ifs), read_char(ifs)}; }
void operator () (std::ofstream& ofs, const label_type& l) { write_char(ofs, l.first()); write_char(ofs, l.second()); }
void operator () (std::ifstream& ifs, bn_label_type& l) { l.compact() = read_int(ifs); }
void operator () (std::ofstream& ofs, const bn_label_type& l) { write_int(ofs, l.compact()); }
template <typename V>
void operator () (std::ifstream& ifs, std::set<V>& vec)
{
if (fourcc(ifs, "OSET")) { return; }
size_t sz = read_size(ifs);
vec.clear();
for (size_t i = 0; i < sz; ++i) {
V tmp;
ref() (ifs, tmp);
vec.emplace(tmp);
}
}
template <typename V>
void operator () (std::ofstream& ofs, const std::set<V>& vec)
{
if (fourcc(ofs, "OSET")) { return; }
write_size(ofs, vec.size());
for (const auto& e: vec) {
ref() (ofs, e);
}
}
template <typename V>
void operator () (std::ifstream& ifs, std::unordered_set<V>& vec)
{
if (fourcc(ifs, "USET")) { return; }
size_t sz = read_size(ifs);
vec.clear();
for (size_t i = 0; i < sz; ++i) {
V tmp;
ref() (ifs, tmp);
vec.emplace(tmp);
}
}
template <typename V>
void operator () (std::ofstream& ofs, const std::unordered_set<V>& vec)
{
if (fourcc(ofs, "USET")) { return; }
write_size(ofs, vec.size());
for (const auto& e: vec) {
ref() (ofs, e);
}
}
void operator () (std::ifstream& ifs, std::vector<bool>::reference i) { i = !!read_char(ifs); }
template <typename V, typename A>
void operator () (std::ifstream& ifs, std::vector<V, A>& vec)
{
if (fourcc(ifs, "VECT")) { return; }
size_t sz = read_size(ifs);
vec.clear();
vec.reserve(sz);
for (size_t i = 0; i < sz; ++i) {
vec.emplace_back();
ref() (ifs, vec.back());
}
}
template <typename V, typename A>
void operator () (std::ofstream& ofs, const std::vector<V, A>& vec)
{
if (fourcc(ofs, "VECT")) { return; }
write_size(ofs, vec.size());
for (const auto& e: vec) {
ref() (ofs, e);
}
}
template <typename K, typename V, typename A, typename C>
void operator () (std::ifstream& ifs, std::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 () (std::ofstream& ofs, const std::map<K, V, A, C>& map)
{
if (fourcc(ofs, "MAP ")) { return; }
write_size(ofs, map.size());
for (const auto& kv: map) {
ref() (ofs, kv.first);
ref() (ofs, kv.second);
}
}
template <typename SCALAR, int ROW, int COL, int C, int D, int E>
void operator () (std::ifstream& 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 () (std::ofstream& ofs, const Eigen::Matrix<SCALAR, ROW, COL, C, D, E>& mat) { write_matrix(ofs, mat); }
void operator () (std::ifstream& ifs, geno_matrix& mat) { read_geno_matrix(ifs, mat); }
void operator () (std::ofstream& ofs, geno_matrix& mat) { write_geno_matrix(ofs, mat); }
void operator () (std::ifstream& ifs, std::shared_ptr<geno_matrix>& ptr)
{
ptr.reset();
ptr = std::make_shared<geno_matrix>();
ref() (ifs, *ptr);
if (!ptr->size()) {
ptr.reset();
}
}
void operator () (std::ofstream& ofs, const std::shared_ptr<geno_matrix>& ptr)
{
if (ptr) {
ref() (ofs, *ptr.get());
} else {
geno_matrix _;
ref() (ofs, _);
}
}
};
struct rw_base : public rw_any<rw_base> {
virtual ~rw_base() {}
using rw_any<rw_base>::fourcc;
using rw_any<rw_base>::ref;
using rw_any<rw_base>::operator ();
};
template <typename PARENT_TYPE, typename STATE_TYPE>
struct rw_comb : public rw_any<rw_comb<PARENT_TYPE, STATE_TYPE>> {
typedef combination_type<PARENT_TYPE, STATE_TYPE> comb_type;
......@@ -475,16 +294,16 @@ struct rw_comb : public rw_any<rw_comb<PARENT_TYPE, STATE_TYPE>> {
using rw_any<rw_comb<PARENT_TYPE, STATE_TYPE>>::ref;
using rw_any<rw_comb<PARENT_TYPE, STATE_TYPE>>::operator ();
void operator () (std::ifstream& fs, typename comb_type::key_type& key) { ref() (fs, key.parent); ref() (fs, key.state); }
void operator () (std::ifstream& ifs, bn_label_type& l) { l.compact() = read_int(ifs); }
void operator () (std::ofstream& ofs, const bn_label_type& l) { write_int(ofs, l.compact()); }
void operator () (std::ifstream& fs, typename comb_type::key_type& key) { ref() (fs, key.parent); ref() (fs, key.state); }
void operator () (std::ofstream& fs, const typename comb_type::key_type& key) { ref() (fs, key.parent); ref() (fs, key.state); }
void operator () (std::ifstream& fs, typename comb_type::key_list& keys) { ref() (fs, keys.keys); }
void operator () (std::ofstream& fs, const typename comb_type::key_list& keys) { ref() (fs, keys.keys); }
void operator () (std::ifstream& fs, typename comb_type::element_type& elem) { ref() (fs, elem.keys); ref() (fs, elem.coef); }
void operator () (std::ofstream& fs, const typename comb_type::element_type& elem) { ref() (fs, elem.keys); ref() (fs, elem.coef); }
void operator () (std::ifstream& fs, comb_type& comb)
......@@ -511,24 +330,26 @@ struct rw_tree : public rw_any<rw_tree> {
void operator () (std::ifstream& fs, pedigree_tree_type::node_type& node) { ref() (fs, node.p1); ref() (fs, node.p2); }
void operator () (std::ofstream& fs, const 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)
{
ref() (fs, tree.m_leaves);
ref() (fs, tree.m_roots);
ref() (fs, tree.m_nodes);
ref() (fs, tree.m_must_recompute);
ref() (fs, tree.m_node_number_to_ind_number);
ref() (fs, tree.m_ind_number_to_node_number);
ref() (fs, tree.m_original_ordering);
}
void operator () (std::ifstream& fs, pedigree_tree_type& tree)
{