Commit cdc52247 authored by Damien Leroux's avatar Damien Leroux
Browse files

Implementing symmetries... Handling reentrants.

parent f3d8169e
......@@ -153,6 +153,42 @@ std::ostream& operator << (std::ostream& os, const std::map<char, int>& L)
return ret;
}
inline
MatrixXb generate_lozenge(size_t n0, size_t n1, const MatrixXb& permut)
{
int N = n0 * n1;
/*if (N == 2) {*/
/*return (MatrixXb::Ones(2, 2) - MatrixXb::Identity(2, 2));*/
/*}*/
std::vector<int> permut_col(permut.cols());
std::vector<int> permut_row(permut.rows());
for (int c = 0; c < permut.cols(); ++c) {
int r;
permut.col(c).maxCoeff(&r);
permut_col[r] = c;
permut_row[c] = r;
}
MatrixXb ret = MatrixXb::Zero(N, N);
for (size_t i0 = 0; i0 < n0; ++i0) {
for (size_t i1 = 0; i1 < n1; ++i1) {
int i = i1 * n0 + i0;
int j = i0 * n1 + i1;
/*ret(permut_col[i], permut_row[j]) = 1;*/
ret(permut_row[i], permut_col[j]) = 1;
}
}
/*MSG_DEBUG("lozenge(" << n0 << ", " << n1 << ',' << (antidiag ? "anti" : "diag") << ')');*/
/*MSG_DEBUG(ret);*/
/*MSG_DEBUG("-- self-adjoint?");*/
/*MSG_DEBUG((ret - ret.transpose()));*/
/*MSG_DEBUG("-- self-inverse?");*/
/*MSG_DEBUG((ret * ret));*/
return ret;
}
struct letter_permutation_type {
std::map<char, char> table;
......@@ -1011,25 +1047,25 @@ std::ostream& operator << (std::ostream& os, const geno_matrix& gm)
std::set<symmetry_table_type>& tmp,
const symmetry_table_type& s1, const symmetry_table_type& s2)
{
MSG_DEBUG_INDENT_EXPR("[combine_sym] ");
/*MSG_DEBUG_INDENT_EXPR("[combine_sym] ");*/
/*MSG_DEBUG("NEW SYM FROM");*/
/*MSG_DEBUG(s1.dump(m1.labels));*/
/*MSG_DEBUG("AND");*/
/*MSG_DEBUG(s2.dump(m2.labels));*/
/*MSG_DEBUG("...");*/
symmetry_table_type stmp = s1 * s2;
MSG_DEBUG(stmp.dump(ret.labels));
/*MSG_DEBUG(stmp.dump(ret.labels));*/
MatrixXd sm = stmp.matrix().cast<double>();
/*MSG_DEBUG(sm);*/
bool result = false;
if (stmp.is_consistent(ret.labels) && (sm.transpose() * ret.inf_mat * sm - ret.inf_mat).isZero(FLOAT_TOL)) {
/*if ((sm.transpose() * ret.inf_mat * sm - ret.inf_mat).isZero(FLOAT_TOL)) {*/
/*MSG_DEBUG("OK!" << std::endl << stmp);*/
MSG_DEBUG("OK!");
/*MSG_DEBUG("OK!");*/
tmp.emplace(stmp);
result = true;
} else {
MSG_DEBUG("NOT GOOD.");
/*MSG_DEBUG("NOT GOOD.");*/
/*MSG_DEBUG(ret.inf_mat);*/
/*MSG_DEBUG("--");*/
/*MSG_DEBUG((sm.transpose() * ret.inf_mat * sm - ret.inf_mat));*/
......@@ -1041,7 +1077,7 @@ std::ostream& operator << (std::ostream& os, const geno_matrix& gm)
/*MSG_DEBUG("switch map " << kv.first << " -> " << kv.second);*/
/*}*/
}
MSG_DEBUG_DEDENT;
/*MSG_DEBUG_DEDENT;*/
return result;
};
......@@ -1177,21 +1213,21 @@ geno_matrix kronecker(const geno_matrix& m1, const geno_matrix& m2)
MatrixXb loz = generate_lozenge(m1.cols(), m2.cols(), r0, r1);
symmetry_table_type solange(loz);
if (solange.is_consistent(ret.labels)) {
MSG_DEBUG("HAVE AN ACTUAL SYMMETRY:");
MSG_DEBUG(solange.dump(ret.labels));
/*MSG_DEBUG("HAVE AN ACTUAL SYMMETRY:");*/
/*MSG_DEBUG(solange.dump(ret.labels));*/
tmp.emplace(solange);
} else if (solange.is_consistent(ret.labels, true)) {
MSG_DEBUG("HAVE A LATENT SYMMETRY:");
MSG_DEBUG(solange.dump(ret.labels, true));
/*MSG_DEBUG("HAVE A LATENT SYMMETRY:");*/
/*MSG_DEBUG(solange.dump(ret.labels, true));*/
latent_tmp.emplace(solange);
} else {
MSG_DEBUG("NOT GOOD:");
MSG_DEBUG_INDENT_EXPR("[actual] ");
MSG_DEBUG(solange.dump(ret.labels, true));
MSG_DEBUG_DEDENT;
MSG_DEBUG_INDENT_EXPR("[latent] ");
MSG_DEBUG(solange.dump(ret.labels, true));
MSG_DEBUG_DEDENT;
/*MSG_DEBUG("NOT GOOD:");*/
/*MSG_DEBUG_INDENT_EXPR("[actual] ");*/
/*MSG_DEBUG(solange.dump(ret.labels, true));*/
/*MSG_DEBUG_DEDENT;*/
/*MSG_DEBUG_INDENT_EXPR("[latent] ");*/
/*MSG_DEBUG(solange.dump(ret.labels, true));*/
/*MSG_DEBUG_DEDENT;*/
/*latent_tmp.emplace(MatrixXb::Identity(m1.cols() * m2.cols(), m1.cols() * m2.cols()));*/
}
}
......@@ -1260,6 +1296,7 @@ geno_matrix ancestor_matrix(char a)
(MatrixXd(1, 1) << 1).finished(),
#ifdef WITH_OVERLUMPING
{{std::map<int, int>{{0, 0}}, letter_permutation_type{{a, a}}}},
/*{{std::map<int, int>{{0, 0}}, letter_permutation_type{{a, a}}}},*/
/*{}*/
#endif
};
......@@ -1609,10 +1646,10 @@ std::ostream& operator << (std::ostream& os, const std::forward_list<T>& l)
std::forward_list<const subset*> stack;
for (const auto& x: P) { stack.push_front(&x); }
size_t count = 0;
MSG_DEBUG_INDENT_EXPR("[derisavi] ");
MSG_DEBUG("+ + + + + + + + + + + + + + + + + + + + + + + + +");
MSG_DEBUG(M);
MSG_DEBUG("+ + + + + + + + + + + + + + + + + + + + + + + + +");
/*MSG_DEBUG_INDENT_EXPR("[derisavi] ");*/
/*MSG_DEBUG("+ + + + + + + + + + + + + + + + + + + + + + + + +");*/
/*MSG_DEBUG(M);*/
/*MSG_DEBUG("+ + + + + + + + + + + + + + + + + + + + + + + + +");*/
while (!stack.empty()) {
++count;
const subset& C = *stack.front();
......@@ -1629,23 +1666,23 @@ std::ostream& operator << (std::ostream& os, const std::forward_list<T>& l)
break;
}
if (0 && subsets.size() > 1) {
MSG_DEBUG("CONFLICT! " << B << ' ' << subsets);
MSG_QUEUE_FLUSH();
/*MSG_DEBUG("CONFLICT! " << B << ' ' << subsets);*/
/*MSG_QUEUE_FLUSH();*/
std::set<subset> sym_C = apply_symmetries(C);
MSG_DEBUG("REMOVING");
MSG_DEBUG("" << sym_C);
MSG_QUEUE_FLUSH();
/*MSG_DEBUG("REMOVING");*/
/*MSG_DEBUG("" << sym_C);*/
/*MSG_QUEUE_FLUSH();*/
stack.remove_if([&](const subset*& s) { return sym_C.find(*s) != sym_C.end(); });
for (const auto& c: sym_C) {
P.erase(c);
}
for (auto& s: subsets) {
for (const auto& sym_s: apply_symmetries(s.second)) {
MSG_DEBUG("ADDING " << sym_s);
MSG_QUEUE_FLUSH();
/*MSG_DEBUG("ADDING " << sym_s);*/
/*MSG_QUEUE_FLUSH();*/
auto io = P.insert(sym_s);
if (!io.second) {
MSG_DEBUG("PROUT SUBSET ALREADY IN P " << (*io.first));
/*MSG_DEBUG("PROUT SUBSET ALREADY IN P " << (*io.first));*/
}
if (io.first->size() > 1) {
stack.push_front(&*io.first);
......@@ -1657,17 +1694,17 @@ std::ostream& operator << (std::ostream& os, const std::forward_list<T>& l)
}
}
MSG_DEBUG("FINAL PARTITION after " << count << " iterations");
MSG_DEBUG("final base cost = " << compute_cost(P));
MSG_DEBUG("" << P);
/*MSG_DEBUG("FINAL PARTITION after " << count << " iterations");*/
/*MSG_DEBUG("final base cost = " << compute_cost(P));*/
/*MSG_DEBUG("" << P);*/
check_partition_against_symmetries(P);
MSG_DEBUG_INDENT_EXPR("[result] ");
/*MSG_DEBUG_INDENT_EXPR("[result] ");*/
experimental_lumper ret(lump_using_partition_weighted(M, P));
MSG_DEBUG(ret.M);
MSG_DEBUG_DEDENT;
MSG_DEBUG_DEDENT;
/*MSG_DEBUG(ret.M);*/
/*MSG_DEBUG_DEDENT;*/
/*MSG_DEBUG_DEDENT;*/
return ret;
}
......@@ -2011,13 +2048,13 @@ std::ostream& operator << (std::ostream& os, const std::forward_list<T>& l)
void split_and_apply_symmetries(std::set<subset>& P, MatrixXb& adj, const subset& splitted, const std::map<double, subset>& splits, std::forward_list<const subset*>& stack)
{
/*MSG_DEBUG("NEW SPLIT IMPLEMENTATION");*/
MSG_DEBUG("stack before " << stack);
MSG_DEBUG("SPLITTING " << splits << " FROM PARTITION");
MSG_DEBUG("" << P);
/*MSG_DEBUG("stack before " << stack);*/
/*MSG_DEBUG("SPLITTING " << splits << " FROM PARTITION");*/
/*MSG_DEBUG("" << P);*/
MSG_DEBUG_INDENT_EXPR("[adj BEFORE] ");
MSG_DEBUG(adj);
MSG_DEBUG_DEDENT;
/*MSG_DEBUG_INDENT_EXPR("[adj BEFORE] ");*/
/*MSG_DEBUG(adj);*/
/*MSG_DEBUG_DEDENT;*/
std::vector<bool> modified(M.cols(), false);
std::vector<bool> stacked(M.cols(), false);
......@@ -2028,11 +2065,11 @@ std::ostream& operator << (std::ostream& os, const std::forward_list<T>& l)
auto do_split = [&] (const subset& S1, const subset& S2)
{
MSG_DEBUG("Splitting " << S1 << " and " << S2);
/*MSG_DEBUG("Splitting " << S1 << " and " << S2);*/
for (int s1: S1) {
for (int s2: S2) {
MSG_DEBUG(s1 << "." << s2);
MSG_QUEUE_FLUSH();
/*MSG_DEBUG(s1 << "." << s2);*/
/*MSG_QUEUE_FLUSH();*/
adj(s2, s1) = adj(s1, s2) = 0;
}
modified[s1] = true;
......@@ -2106,10 +2143,10 @@ std::ostream& operator << (std::ostream& os, const std::forward_list<T>& l)
}
}
#endif
MSG_DEBUG_INDENT_EXPR("[adj AFTER] ");
MSG_DEBUG(adj);
MSG_DEBUG_DEDENT;
MSG_QUEUE_FLUSH();
/*MSG_DEBUG_INDENT_EXPR("[adj AFTER] ");*/
/*MSG_DEBUG(adj);*/
/*MSG_DEBUG_DEDENT;*/
/*MSG_QUEUE_FLUSH();*/
std::vector<bool> visited(M.cols(), false);
......@@ -2129,15 +2166,15 @@ std::ostream& operator << (std::ostream& os, const std::forward_list<T>& l)
enstack |= modified[i];
}
}
MSG_DEBUG("BUILT PART " << s);
/*MSG_DEBUG("BUILT PART " << s);*/
auto it = P.insert(s).first;
if (enstack || stacked[it->front()]) {
stack.push_front(&*it);
MSG_DEBUG("RESTACK " << (*stack.front()));
/*MSG_DEBUG("RESTACK " << (*stack.front()));*/
}
}
MSG_DEBUG("stack after " << stack);
/*MSG_DEBUG("stack after " << stack);*/
check_partition_against_symmetries(P);
}
......
......@@ -8,7 +8,7 @@ template <typename PARENT_TYPE, typename STATE_TYPE=size_t>
struct combination_type {
static const size_t nopar = (size_t) -1;
typedef combination_type<PARENT_TYPE> this_type;
typedef combination_type<PARENT_TYPE, STATE_TYPE> this_type;
typedef Eigen::Matrix<this_type, Eigen::Dynamic, 1> Vector;
......
This diff is collapsed.
......@@ -30,8 +30,8 @@ COV_OBJ=$(subst .cc,.cov.o,$(SRC))
ALL_BUT_MAIN_OBJ=$(subst main.o ,,$(OBJ))
DEBUG_OPTS=-ggdb
#DEBUG_OPTS=-ggdb -DNDEBUG
#DEBUG_OPTS=-ggdb
DEBUG_OPTS=-ggdb -DNDEBUG
#OPT_OPTS=-O3 -DEIGEN_NO_DEBUG -DNDEBUG
#OPT_OPTS=-O3 -DNDEBUG
#PROF_OPTS=-O2 -DEIGEN_NO_DEBUG -DNDEBUG -g -pg
......
......@@ -180,6 +180,7 @@ int main(int argc, char** argv)
auto test_case = [&] () { bool ret = arg & 1; arg >>= 1; return ret; };
/* 1 */
if (test_case())
{
pedigree_type ped(10);
......@@ -199,9 +200,11 @@ int main(int argc, char** argv)
size_t BCL = ped.crossing(A, F1); (void) BCL;
}
/* 2 */
if (test_case())
{
pedigree_type ped(10);
ped.n_alleles = 2;
MSG_DEBUG("#####################################################################");
MSG_DEBUG("A");
size_t A = ped.ancestor();
......@@ -225,11 +228,17 @@ int main(int argc, char** argv)
size_t F3 = ped.selfing(F2);
MSG_DEBUG((*ped.get_gen(F3)));
MSG_DEBUG("#####################################################################");
MSG_DEBUG("F4");
size_t F4 = ped.selfing(F3);
MSG_DEBUG((*ped.get_gen(F4)));
/*MSG_DEBUG("test F3 vs old impl");*/
/*geno_matrix old_F3 = lump(~~(*ped.get_gen(A) | *ped.get_gen(B)));*/
/*MSG_DEBUG((old_F3.inf_mat - ped.get_gen(F3)->inf_mat));*/
}
/* 4 */
if (test_case())
{
MSG_DEBUG("#####################################################################");
......@@ -238,6 +247,7 @@ int main(int argc, char** argv)
MSG_DEBUG("#####################################################################");
MSG_DEBUG("#####################################################################");
pedigree_type ped;
ped.n_alleles = 2;
MSG_DEBUG("A");
size_t A = ped.ancestor();
MSG_DEBUG("#####################################################################");
......@@ -266,6 +276,7 @@ int main(int argc, char** argv)
(void) F2_CP;
}
/* 8 */
if (test_case())
{
MSG_DEBUG("");
......@@ -301,6 +312,7 @@ int main(int argc, char** argv)
MSG_DEBUG(ped.get_gen(G)->lim_inf());
}
/* 16 */
if (test_case())
{
pedigree_type ped;
......@@ -313,9 +325,12 @@ int main(int argc, char** argv)
size_t F2 = ped.selfing(F1); (void) F2;
}
/* 32 */
if (test_case())
{
pedigree_type ped(10);
/*ped.max_states = 4;*/
/*ped.n_alleles = 2;*/
size_t A = ped.ancestor();
size_t B = ped.ancestor();
size_t F1_1 = ped.crossing(A, B);
......@@ -326,13 +341,13 @@ int main(int argc, char** argv)
size_t F3_1 = ped.crossing(F2_1, F2_2);
size_t F3_2 = ped.crossing(F2_1, F2_2);
MSG_DEBUG("F4s");
size_t F4_1 = ped.crossing(F3_1, F3_2);
size_t F4_2 = ped.crossing(F3_1, F3_2);
MSG_DEBUG("F5s");
size_t F5_1 = ped.crossing(F4_1, F4_2);
(void) F4_1; (void) F4_2;
size_t F4_1 = ped.crossing(F3_1, F3_2); (void) F4_1;
/*size_t F4_2 = ped.crossing(F3_1, F3_2); (void) F4_2;*/
/*MSG_DEBUG("F5s");*/
/*size_t F5_1 = ped.crossing(F4_1, F4_2); (void) F5_1;*/
}
/* 64 */
if (test_case())
{
pedigree_type ped(20);
......@@ -351,6 +366,7 @@ int main(int argc, char** argv)
size_t S1 = ped.selfing(S); (void) S1;
}
/* 128 */
if (test_case())
{
pedigree_type ped(20);
......@@ -372,9 +388,11 @@ int main(int argc, char** argv)
size_t S = ped.selfing(X); (void) S;
}
/* 256 */
if (test_case())
{
pedigree_type ped(20);
ped.n_alleles = 2;
size_t A = ped.ancestor();
size_t B = ped.ancestor();
size_t C = ped.ancestor();
......@@ -392,6 +410,7 @@ int main(int argc, char** argv)
/*size_t S = ped.selfing(X); (void) S;*/
}
/* 512 */
if (test_case())
{
pedigree_type ped(20);
......@@ -412,6 +431,7 @@ int main(int argc, char** argv)
size_t S = ped.selfing(X); (void) S;
}
/* 1024 */
if (test_case())
{
pedigree_type ped(20);
......@@ -424,6 +444,7 @@ int main(int argc, char** argv)
size_t S = ped.selfing(X); (void) S;
}
/* 2048 */
if (test_case())
{
pedigree_type ped(20);
......@@ -439,6 +460,7 @@ int main(int argc, char** argv)
size_t S = ped.selfing(X); (void) S;
}
/* 4096 */
if (test_case())
{
pedigree_type ped(20);
......@@ -454,6 +476,7 @@ int main(int argc, char** argv)
size_t S = ped.selfing(X); (void) S;
}
/* 8192 */
if (test_case())
{
auto A = ancestor_matrix('a');
......@@ -472,6 +495,7 @@ int main(int argc, char** argv)
MSG_DEBUG(kronecker(F1, gamete));
}
/* 16384 */
if (test_case())
{
pedigree_type ped(30);
......@@ -494,6 +518,7 @@ int main(int argc, char** argv)
size_t Y = ped.selfing(X); (void) Y;
}
/* 32768 */
if (test_case())
{
pedigree_type ped(30);
......
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