Commit 19e8ed44 authored by Damien Leroux's avatar Damien Leroux
Browse files

Promising developments regarding the overlumping.

parent ff92a1bd
......@@ -248,6 +248,14 @@ struct braille_grid {
line(int x0, int y0, int x1, int y1, int dash_length, int space_length, int r=255, int g=255, int b=255)
{
int counter = 0;
return line(x0, y0, x1, y1, dash_length, space_length, r, g, b, counter);
}
braille_grid&
line(int x0, int y0, int x1, int y1, int dash_length, int space_length, int r, int g, int b, int& counter)
{
/*int counter = 0;*/
space_length += dash_length;
if (abs(x1 - x0) >= abs(y1 - y0)) {
if (x1 < x0) {
......@@ -308,6 +316,8 @@ struct braille_grid {
}
}
}
counter += space_length - 1;
counter %= space_length;
return *this;
}
......@@ -545,8 +555,9 @@ struct braille_plot : public braille_grid {
}
braille_plot&
plot(const std::vector<double>& X, const std::vector<double>& Y, int R=255, int G=255, int B=255)
plot(const std::vector<double>& X, const std::vector<double>& Y, int R=255, int G=255, int B=255, int dash_length=1, int space_length=0)
{
int counter = 0;
if (X.size() == Y.size()) {
if (X.size() == 1) {
put_pixel(X[0], Y[0], R, G, B);
......@@ -554,7 +565,8 @@ struct braille_plot : public braille_grid {
double x0 = X[0];
double y0 = Y[0];
for (size_t i = 1; i < X.size(); ++i) {
line(x0, y0, X[i], Y[i], 1, 0, R, G, B);
/*MSG_DEBUG("plot line " << x0 << ',' << y0 << '-' << X[i] << ',' << Y[i]);*/
line(x0, y0, X[i], Y[i], dash_length, space_length, R, G, B, counter);
x0 = X[i];
y0 = Y[i];
}
......@@ -570,6 +582,13 @@ struct braille_plot : public braille_grid {
return *this;
}
braille_plot&
line(double x0, double y0, double x1, double y1, int dash_length, int space_length, int R, int G, int B, int& counter)
{
braille_grid::line(x2grid(x0), y2grid(y0), x2grid(x1), y2grid(y1), dash_length, space_length, R, G, B, counter);
return *this;
}
braille_plot&
hline(double y, int dash_length, int space_length, int R=255, int G=255, int B=255)
{
......
#ifndef _SPEL_GENOPROB_COMPUTER_H_
#define _SPEL_GENOPROB_COMPUTER_H_
#include "geno_matrix.h"
#include "labelled_matrix.h"
struct geno_prob_computer {
const std::vector<double>* locus;
const MatrixXd* LV;
const geno_matrix* this_gen;
std::unordered_map<double, MatrixXd> tr_cache;
std::vector<MatrixXd> TR;
MatrixXd init_kernel;
/*MatrixXd redux;*/
std::map<int, MatrixXd> left, right;
MatrixXd redux;
VectorXd inv_stat_dist;
std::vector<label_type> unique_labels;
geno_prob_computer()
: locus(0), LV(0), this_gen(0), TR(), init_kernel(), left(), right(), redux(), inv_stat_dist()
{}
geno_prob_computer(geno_prob_computer&& tmp)
: locus(tmp.locus)
, LV(tmp.LV)
, this_gen(tmp.this_gen)
, left(std::forward<std::map<int, MatrixXd>>(tmp.left))
, right(std::forward<std::map<int, MatrixXd>>(tmp.right))
, redux(tmp.redux)
, inv_stat_dist(tmp.inv_stat_dist)
{}
geno_prob_computer&
operator = (geno_prob_computer&& tmp)
{
locus = tmp.locus;
LV = tmp.LV;
this_gen = tmp.this_gen;
TR.swap(tmp.TR);
left.swap(tmp.left);
right.swap(tmp.right);
redux = tmp.redux;
inv_stat_dist = tmp.inv_stat_dist;
return *this;
}
MatrixXd& compute_tr(double dist)
{
auto it = tr_cache.find(dist);
if (it == tr_cache.end()) {
return tr_cache[dist] = this_gen->exp(dist);
}
return it->second;
}
MatrixXd normalize_by_col(const MatrixXd& m)
{
MatrixXd ret = MatrixXd::Zero(m.innerSize(), m.outerSize());
VectorXd sums = m.colwise().sum();
for (int i = 0; i < m.outerSize(); ++i) {
if (sums(i) != 0) {
ret.col(i) = m.col(i) / sums(i);
}
}
return ret;
}
MatrixXd lv(int i)
{
return LV->col(i).asDiagonal();
}
void init(const geno_matrix* gen)
{
/*tr_cache.clear();*/
this_gen = gen;
/*redux = make_redux(matrix_labels_to_cliques(break_matrix.column_labels));*/
experimental_lumper el(*this_gen);
redux = el.make_collect(el.partition_on_labels(), this_gen->inf_mat.cols());
inv_stat_dist = this_gen->stat_dist;
for (int i = 0; i < inv_stat_dist.size(); ++i) {
inv_stat_dist(i) = 1. / inv_stat_dist(i);
}
std::set<label_type> latemp(this_gen->labels.begin(), this_gen->labels.end());
unique_labels.assign(latemp.begin(), latemp.end());
TR.clear();
for (size_t i = 1; i < locus->size(); ++i) {
TR.emplace_back(compute_tr((*locus)[i] - (*locus)[i - 1]));
}
}
bool compute_partial_kernels()
{
std::map<int, std::string> lstring, rstring;
left.clear();
right.clear();
int l = 0, r = LV->outerSize() - 1;
MatrixXd kernel;
left[0] = kernel = (this_gen->stat_dist.array() * LV->col(0).array()).matrix();
lstring[0] = "LV0.";
/*MSG_DEBUG("0 " << lstring[0]);*/
/*MSG_DEBUG(left[0]);*/
for(++l; l <= r; ++l) {
kernel = LV->col(l).asDiagonal() * TR[l - 1] * kernel;
double s = kernel.sum();
if (s == 0) {
return false;
}
kernel /= s;
left[l] = kernel;
std::stringstream ss;
ss << "LV" << l << ".TR" << (l-1) << '.' << lstring[l-1];
lstring[l] = ss.str();
/*MSG_DEBUG(l << " " << lstring[l]);*/
/*MSG_DEBUG(left[l]);*/
}
right[r+1] = VectorXd::Ones(this_gen->rows()) / this_gen->rows();
/*right[r+1] = this_gen->stat_dist;*/
rstring[r+1] = "";
right[r-1] = right[r] = kernel = (this_gen->stat_dist.array() * LV->col(r).array()).matrix(); //LV->col(r) / LV->col(r).sum();
{ std::stringstream s; s << ".LV" << r; rstring[r-1] = rstring[r] = s.str(); }
for (--r; r > 0; --r) {
kernel = LV->col(r).asDiagonal() * TR[r] * kernel;
double s = kernel.sum();
if (s == 0) {
return false;
}
kernel /= s;
right[r-1] = kernel;
std::stringstream ss;
ss << rstring[r] << ".TR" << r << ".LV" << r;
rstring[r-1] = ss.str();
}
r = LV->outerSize() - 1;
kernel = init_kernel;
return true;
}
std::pair<size_t, size_t> find_lr(double loc)
{
size_t l = 0, r = LV->outerSize() - 1;
if (loc < locus->front() || loc > locus->back()) {
throw out_of_segment_exception();
}
for (; (l + 1) < r && (*locus)[l + 1] <= loc; ++l);
for (; r > 0 && (*locus)[r - 1] > loc; --r);
return {l, r};
}
VectorXd fast_compute_at(double loc, size_t l, size_t r)
{
#if 0
MSG_DEBUG("l=" << l);
MSG_DEBUG("r=" << r);
MSG_DEBUG("locus_l=" << (*locus)[l]);
MSG_DEBUG("locus_r=" << (*locus)[r]);
MSG_DEBUG("locus=" << loc);
#endif
MatrixXd& TRL = compute_tr(loc - (*locus)[l]);
MatrixXd& TRR = compute_tr((*locus)[r] - loc);
#if 0
MSG_DEBUG(MATRIX_SIZE(TRL));
MSG_DEBUG(TRL);
MSG_DEBUG(MATRIX_SIZE(TRR));
MSG_DEBUG(TRR);
MSG_DEBUG(MATRIX_SIZE(left[l]));
MSG_DEBUG(left[l]);
MSG_DEBUG(MATRIX_SIZE(right[l]));
MSG_DEBUG(right[l]);
#endif
auto L = (TRL * left[l]).array();
auto R = (TRR * right[l]).array();
#if 0
/*std::cout << loc << std::endl << (L * R).matrix().transpose() << std::endl;*/
MSG_DEBUG(MATRIX_SIZE(L));
MSG_DEBUG(L);
MSG_DEBUG(MATRIX_SIZE(R));
MSG_DEBUG(R);
MSG_DEBUG(MATRIX_SIZE(this_gen->get_redux()));
MSG_DEBUG(this_gen->get_redux());
#endif
VectorXd vv = redux * (L * R * inv_stat_dist.array()).matrix();
double s = vv.sum();
if (s == 0) {
return vv;
}
return vv / s;
}
labelled_matrix<MatrixXd, label_type, double>
fast_compute_over_segment(const std::vector<double>& steps)
{
if (!compute_partial_kernels()) {
labelled_matrix<MatrixXd, label_type, double>
ret(unique_labels, steps);
ret.data = MatrixXd::Zero(ret.row_labels.size(), ret.column_labels.size());
return ret;
}
#if 0
std::vector<typename GENERATION::col_label_type> labels;
size_t i = 0;
size_t j = 0;
while (i < redux.innerSize()) {
while (j < breakmat.column_labels.size() && redux(i, j) != 1) ++j;
labels.push_back(breakmat.column_labels[j]);
while (j < breakmat.column_labels.size() && redux(i, j) == 1) ++j;
++i;
}
labelled_matrix<MatrixXd, typename GENERATION::col_label_type, double>
ret(labels, steps);
#endif
/*MSG_DEBUG("LOCUS VECTORS" << std::endl << *LV);*/
labelled_matrix<MatrixXd, label_type, double>
ret(unique_labels, steps);
/*MSG_DEBUG(MATRIX_SIZE(ret));*/
size_t l = 0, last = steps.size() - 1, r = !!(LV->outerSize() - 1);
for (size_t i = 0; i < steps.size(); ++i) {
/* FIXME: optimize the CORRECT initialization of l & r */
std::tie(l, r) = find_lr(steps[i]);
ret.data.col(i) = fast_compute_at(steps[i], l, r);
if ((ret.data.col(i).array() < 0).any()) {
MSG_DEBUG("Bad value at column #" << i << std::endl << ret.data.col(i).transpose());
MSG_DEBUG((*this_gen));
#if 0
MSG_DEBUG("steps: " << steps);
MSG_DEBUG("LV:" << std::endl << (*LV));
MSG_DEBUG(*this_gen);
for (const auto& m: TR) {
MSG_DEBUG(m);
}
for (const auto& kv: left) {
MSG_DEBUG("left " << kv.first);
MSG_DEBUG(kv.second);
}
for (const auto& kv: right) {
MSG_DEBUG("right " << kv.first);
MSG_DEBUG(kv.second);
}
#endif
/*throw 0;*/
}
while (steps[i] >= (*locus)[r]) {
++l;
if (i < last) {
++r;
} else {
break;
}
/*r += i < last;*/
}
}
/*MSG_DEBUG("STATE PROBABILITIES" << std::endl << ret);*/
return ret;
}
};
#endif
......@@ -15,6 +15,7 @@
extern "C" {
#include <unistd.h>
#include <sys/ioctl.h>
#include <signal.h>
}
#define _WHITE "\x1b[37;1m"
......@@ -164,6 +165,9 @@ struct message_queue : public ostream_manager {
std::condition_variable m_condition;
std::condition_variable m_flush_condition;
static void catch_SIGSEG(int signo);
static sighandler_t& old_sighandler() { static sighandler_t _; return _; }
bool m_stop;
std::thread m_thread;
......@@ -176,6 +180,7 @@ struct message_queue : public ostream_manager {
, m_stop(false)
, m_thread([this] () { run(); })
{
old_sighandler() = signal(SIGSEGV, catch_SIGSEG);
}
~message_queue()
......@@ -206,7 +211,6 @@ struct message_queue : public ostream_manager {
while (!m_queue.empty()) {
m_flush_condition.wait(lock);
}
}
};
......@@ -427,6 +431,15 @@ void message_queue::run()
#define MSG_DEBUG_DEDENT CREATE_MESSAGE(msg_channel::Log, _DBG_DEDENT_MARK_S)
void message_queue::catch_SIGSEG(int)
{
MSG_DEBUG("**SEGFAULT DETECTED**");
MSG_QUEUE_FLUSH();
signal(SIGSEGV, old_sighandler());
raise(SIGSEGV);
}
struct scoped_indent {
scoped_indent() { MSG_DEBUG_INDENT; }
scoped_indent(const std::string& str) { MSG_DEBUG_INDENT_EXPR(str); }
......
This diff is collapsed.
......@@ -18,7 +18,7 @@
/*typedef std::set<int> subset;*/
typedef std::vector<int> subset;
#if 1
#if 0
#define _EPSILON (1.e-10)
namespace std {
......
......@@ -502,9 +502,10 @@ struct pedigree_type : public archivable {
size_t last_node_index() const { return tree.size() - 1; }
#endif
individual_index_type spawn_gamete(const std::string& gamete_name, int parent_node)
individual_index_type spawn_gamete(const std::string&, int parent_node)
{
int n = tree.add_node(parent_node);
node_generations.emplace_back(node_generations[parent_node]);
/*MSG_DEBUG_INDENT_EXPR("[compute " << gamete_name << " gamete] ");*/
/*compute_generation(n);*/
/*compute_LC(n);*/
......@@ -585,6 +586,7 @@ struct pedigree_type : public archivable {
void propagate_symmetries(int n, geno_matrix& gen)
{
MSG_DEBUG_INDENT_EXPR("[propagate symmetries #" << n << "] ");
std::vector<int> in, out;
auto expr = tree.extract_expression(n, in, out);
std::vector<pedigree_tree_type> input_trees;
......@@ -592,26 +594,42 @@ struct pedigree_type : public archivable {
for (int t: in) {
input_trees.emplace_back(tree.extract_subtree(t));
}
/*auto recompute = tree.get_deep_recompute_vec(n);*/
/*MSG_DEBUG("RECOMPUTE: " << recompute);*/
auto get_lumper
= [&, this] (int node) -> MatrixXb
{
/*if (tree[node].is_gamete() || recompute[node]) {*/
/*MSG_DEBUG("NIL lumper for node #" << node << " because" << (recompute[node] && tree[node].is_gamete() ? " recompute flag is set and it is a gamete" : tree[node].is_gamete() ? " it is a gamete" : " recompute flag is set"));*/
/*return {};*/
/*}*/
return get_node_gen(node)->collect.cast<bool>();
};
symmetry_propagator sp(expr);
MSG_DEBUG_INDENT_EXPR("[SYMMETRIES] ");
gen.symmetries = sp.compute_propagated_symmetries(
[&] (int i) -> const pedigree_tree_type& { MSG_DEBUG("request subtree for input #" << i); return input_trees[i]; },
[&] (int node) -> const symmetry_group_type& { MSG_DEBUG("request symmetries for node #" << node); return get_node_gen(node)->symmetries; },
[&] (int node) -> const symmetry_group_type& { MSG_DEBUG("request latent symmetries for node #" << node); return get_node_gen(node)->latent_symmetries; },
gen.labels, gen.inf_mat,
[this] (int node) -> MatrixXb { return get_node_gen(node)->collect.cast<bool>(); }
get_lumper
);
MSG_DEBUG_DEDENT;
MSG_DEBUG_INDENT_EXPR("[LATENT SYMMETRIES] ");
auto temp = sp.compute_propagated_latent_symmetries(
[&] (int i) -> const pedigree_tree_type& { MSG_DEBUG("request subtree for input #" << i); return input_trees[i]; },
[&] (int node) -> const symmetry_group_type& { MSG_DEBUG("request symmetries for node #" << node); return get_node_gen(node)->symmetries; },
[&] (int node) -> const symmetry_group_type& { MSG_DEBUG("request latent symmetries for node #" << node); return get_node_gen(node)->latent_symmetries; },
gen.labels, gen.inf_mat,
[this] (int node) -> MatrixXb { return get_node_gen(node)->collect.cast<bool>(); }
get_lumper
);
MSG_DEBUG_INDENT_EXPR("[AFTER SYMMETRY PROPAGATION] ");
MSG_DEBUG_DEDENT;
MSG_DEBUG(temp);
gen.latent_symmetries = temp - gen.symmetries;
MSG_DEBUG_INDENT_EXPR("[AFTER SYMMETRY PROPAGATION] ");
MSG_DEBUG(gen);
MSG_DEBUG_DEDENT;
MSG_DEBUG_DEDENT;
}
void compute_generation(int n)
......@@ -643,13 +661,13 @@ struct pedigree_type : public archivable {
/*MSG_DEBUG("TMP GAMETE GEN");*/
/*MSG_DEBUG(new_gen);*/
} else {
auto ngp1 = node_generations[np1];
auto ngp2 = node_generations[np2];
/*auto ngp1 = node_generations[np1];*/
/*auto ngp2 = node_generations[np2];*/
/*auto gp1 = generations[ngp1];*/
/*auto gp2 = generations[ngp2];*/
MSG_DEBUG("Child of " << tree.node2ind(tree.get_p1(np1)) << " and " << tree.node2ind(tree.get_p1(np2)));
cached_gen = &cache_geno[{ngp1, ngp2}];
cached_gen = &cache_geno[{np1, np2}];
if (*cached_gen) {
MSG_DEBUG("GENERATION HAS ALREADY BEEN COMPUTED");
node_generations[n] = *cached_gen;
......@@ -691,29 +709,29 @@ struct pedigree_type : public archivable {
MSG_QUEUE_FLUSH();
/*if (!(ind_number_to_node_number.size() == 9 && tree.size() == 23)) {*/
MSG_DEBUG("PROPAGATING SYMMETRIES");
/*MSG_DEBUG("PROPAGATING SYMMETRIES");*/
/*propagate_symmetries(new_gen, recompute, n);*/
/*study_expression_symmetries(new_gen);*/
/*complete_symmetries(new_gen);*/
MSG_DEBUG("COMPUTING LATENT SYMMETRY");
/*MSG_DEBUG("COMPUTING LATENT SYMMETRY");*/
/*}*/
}
node_generations[n] = generations.size();
generations.emplace_back(new geno_matrix());
MSG_DEBUG("BEFORE LUMPING");
MSG_DEBUG(new_gen);
/*MSG_DEBUG("BEFORE LUMPING");*/
/*MSG_DEBUG(new_gen);*/
*generations.back() = lump(new_gen, max_states);
if (tree[n].is_crossing()) {
propagate_symmetries(n, *generations.back());
} else if (tree[n].is_ancestor()) {
/*if (tree[n].is_crossing()) {*/
/*propagate_symmetries(n, *generations.back());*/
/*} else if (tree[n].is_ancestor()) {*/
generations.back()->symmetries = symmetry_group_type(generations.back()->labels);
}
/*}*/
/**node_generations[n] = lump(new_gen);*/
if (cached_gen) {
*cached_gen = node_generations[n];
}
MSG_DEBUG("DONE COMPUTING GENERATION FOR NODE #" << n);
MSG_DEBUG_INDENT_EXPR("[RESULT " << tree.make_node_label(n) << "] ");
MSG_DEBUG_INDENT_EXPR("[RESULT " << tree.make_node_label(n) << " gen#" << node_generations[n] << "] ");
MSG_DEBUG((*generations.back()));
MSG_DEBUG_DEDENT;
/*MSG_DEBUG((*generations[node_generations[n]]));*/
......@@ -844,7 +862,7 @@ struct pedigree_type : public archivable {
template <typename FIELD_TYPE>
FIELD_TYPE eval(size_t node, FIELD_TYPE geno_matrix::* field, KronFunc func, const std::vector<bool>& recompute, std::vector<bool>& visited) const
{
/*scoped_indent _;*/
scoped_indent _;
/*MSG_DEBUG("eval node " << node);*/
if (visited[node]) {
/*MSG_DEBUG("already visited => 1");*/
......@@ -988,8 +1006,8 @@ struct pedigree_type : public archivable {
} else {
node_to_iterator[node] = iterators.size();
iterators.emplace_back((this->*accessor)(true, node, gamete));
MSG_DEBUG("GAMETE");
MSG_DEBUG("" << (this->*accessor)(true, node, gamete));
/*MSG_DEBUG("GAMETE");*/
/*MSG_DEBUG("" << (this->*accessor)(true, node, gamete));*/
}
#if 0
switch (nodes[node].type) {
......@@ -1007,8 +1025,8 @@ struct pedigree_type : public archivable {
} else {
node_to_iterator[node] = iterators.size();
iterators.emplace_back((this->*accessor)(false, node, *generations[node_generations[node]]));
MSG_DEBUG("GENERATION");
MSG_DEBUG("" << (this->*accessor)(false, node, *generations[node_generations[node]]));
/*MSG_DEBUG("GENERATION");*/
/*MSG_DEBUG("" << (this->*accessor)(false, node, *generations[node_generations[node]]));*/
}
}
......
This diff is collapsed.
......@@ -284,6 +284,12 @@ struct permutation_type : public permutation_base<permutation_type> {
for (++i; i != j; ++i) {
if (b != bins[m_map[*i]]) {
/* whine */
MSG_DEBUG("CAN'T LUMP BECAUSE OF CONFLICT");
MSG_DEBUG("PERMUTATION");
MSG_DEBUG("" << (*this));
MSG_DEBUG("PARTITION");
MSG_DEBUG("" << P0);
MSG_QUEUE_FLUSH();
return {};
}
}
......@@ -294,6 +300,10 @@ struct permutation_type : public permutation_base<permutation_type> {
ret.m_map.push_back(bins[m_map[*S.begin()]]);
}
ret.compute_transposed();
MSG_DEBUG("LUMP");
MSG_DEBUG("" << (*this));
MSG_DEBUG("INTO");
MSG_DEBUG("" << ret);
return ret;
}
......@@ -360,7 +370,15 @@ template <typename OTHER_DERIVED>
/*MSG_DEBUG(other);*/
/*MSG_DEBUG_DEDENT;*/
/*MSG_QUEUE_FLUSH();*/
assert(size() == other.size());
if (size() != other.size()) {
MSG_DEBUG("SIZE MISMATCH " << size() << " vs " << other.size());
MSG_DEBUG("this:");
MSG_DEBUG("" << (*this));
MSG_DEBUG("other:");
MSG_DEBUG("" << other);
MSG_QUEUE_FLUSH();
assert(size() == other.size());
}
ret.m_map.resize(size(), -1);
ret.m_transposed.resize(size(), -1);
const std::vector<int>& m_map = map();
......
......@@ -260,6 +260,11 @@ struct symmetry_table_type {
MSG_DEBUG("* permutation");
MSG_DEBUG(permut);
MSG_QUEUE_FLUSH();
if (permut.size() != labels.size()) {
return {false, {}};
}
bool ok = true;
auto new_labels = permut * labels;
......@@ -771,10 +776,14 @@ struct symmetry_group_type {
symmetry_group_type
permute(const permutation_base<DERIVED>& permut) const
{
MSG_DEBUG_INDENT_EXPR("[group.permute] ");
MSG_DEBUG((*this));
MSG_DEBUG("" << permut);
symmetry_group_type ret;
for (const auto& sym: m_symmetries) {
ret.insert({permut * sym.table, sym.letters});