Commit 43ca8be5 authored by Damien Leroux's avatar Damien Leroux
Browse files

On our way to implement efficient permutations.

parent cdc52247
......@@ -189,6 +189,64 @@ std::ostream& operator << (std::ostream& os, const std::map<char, int>& L)
return ret;
}
inline
MatrixXb generate_lozenge(const std::vector<size_t> sizes&, const std::vector<size_t>& idx_permut, const std::vector<bool>& reverse_gamete)
{
std::vector<size_t> cursors(sizes.size(), 0);
size_t sz = 1;
for (size_t s: sizes) { sz *= s; }
int start = cursors.size() - 1;
auto idx1 = [&] ()
{
size_t ret = 0;
for (int i = start; i >= 0; --i) {
ret *= sizes[i];
ret += cursors[i];
}
return ret;
};
auto idx2 = [&] () {
size_t ret = 0;
for (int i = start; i >= 0; --i) {
ret *= sizes[idx_permut[i]];
size_t curs = cursors[idx_permut[i]];
if (reverse_gamete[i]) {
curs = sizes[i] - 1 - curs;
}
ret += curs;
}
return ret;
};
auto next = [&] ()
{
auto i = cursors.rbegin();
auto j = cursors.rend();
auto si = sizes.rbegin();
for (; i != j; ++i) {
++*i;
if (*i == *si) {
*i = 0;
} else {
break;
}
}
return i != j;
};
MatrixXb ret = MatrixXb::Zero(sz, sz);
do {
ret(idx2(), idx1()) = 1;
} while (next());
return ret;
}
struct letter_permutation_type {
std::map<char, char> table;
......@@ -218,6 +276,19 @@ std::ostream& operator << (std::ostream& os, const std::map<char, int>& L)
}
}
static
letter_permutation_type identity(const std::vector<label_type>& labels)
{
letter_permutation_type ret;
for (const auto& l: labels) {
ret.table[l.first] = l.first;
ret.table[l.second] = l.second;
}
return ret;
}
bool empty() const { return table.size() == 0; }
MatrixXb matrix() const
{
MatrixXb ret = MatrixXb::Zero(table.size(), table.size());
......@@ -227,8 +298,8 @@ std::ostream& operator << (std::ostream& os, const std::map<char, int>& L)
indices[kv.first] = sz;
}
for (const auto& kv: table) {
MSG_DEBUG(kv.second << '[' << indices[kv.second] << "] " << kv.first << '[' << indices[kv.first] << ']');
MSG_QUEUE_FLUSH();
/*MSG_DEBUG(kv.second << '[' << indices[kv.second] << "] " << kv.first << '[' << indices[kv.first] << ']');*/
/*MSG_QUEUE_FLUSH();*/
ret(indices[kv.second], indices[kv.first]) = 1;
}
return ret;
......@@ -320,9 +391,14 @@ std::ostream& operator << (std::ostream& os, const std::map<char, int>& L)
MatrixXb m = matrix() * other.matrix();
if (!(m.transpose() * m - MatrixXb::Identity(m.rows(), m.cols())).isZero()) {
MSG_DEBUG("FOIRURE MATRICE DE PERMUTATION DE LETTRES COMBINEE");
MSG_DEBUG(matrix());
MSG_DEBUG("--");
MSG_DEBUG(other.matrix());
MSG_DEBUG("--");
MSG_DEBUG(m);
MSG_QUEUE_FLUSH();
abort();
/*abort();*/
return {};
}
return {*this, m};
}
......@@ -350,6 +426,25 @@ std::ostream& operator << (std::ostream& os, const std::map<char, int>& L)
/*}*/
}
symmetry_table_type(const MatrixXb& state_matrix, const std::vector<int>& permut_col, const std::vector<int>& permut_row)
: table(), letters()
{
/*std::set<std::pair<char, char>> tmp_switches;*/
for (int j = 0; j < state_matrix.cols(); ++j) {
for (int i = 0; i < state_matrix.rows(); ++i) {
if (state_matrix(i, j)) {
/*MSG_DEBUG("" << labels[j] << " -> " << labels[i]);*/
/*tmp_switches.emplace(labels[j].first, labels[i].first);*/
/*if (NOT_GAMETE(labels[i].second) && NOT_GAMETE(labels[j].second)) {*/
/*tmp_switches.emplace(labels[j].second, labels[i].second);*/
/*}*/
table[permut_col[j]] = permut_row[i];
/*MSG_DEBUG(i << ',' << j << " => " << permut_col[j] << ',' << permut_row[i]);*/
}
}
}
}
symmetry_table_type(/*const std::vector<label_type>& labels, */const MatrixXb& state_matrix, const letter_permutation_type& l)
: /*switches(),*/ table(), letters(l)
{
......@@ -441,6 +536,11 @@ std::ostream& operator << (std::ostream& os, const std::map<char, int>& L)
MSG_DEBUG("letters");
MSG_DEBUG(letters);
MSG_DEBUG_DEDENT;
for (const auto& s: smap) {
if (s.second == '!') {
return false;
}
}
return smap == letters.table;
/*MSG_DEBUG("switches.size() = " << switches.size());*/
/*MSG_DEBUG("switch_map().size() = " << smap.size());*/
......@@ -1042,6 +1142,7 @@ std::ostream& operator << (std::ostream& os, const geno_matrix& gm)
#ifdef WITH_OVERLUMPING
#if 0
inline
bool combine_sym(const geno_matrix& ret, const geno_matrix& m1, const geno_matrix& m2,
std::set<symmetry_table_type>& tmp,
......@@ -1080,8 +1181,9 @@ std::ostream& operator << (std::ostream& os, const geno_matrix& gm)
/*MSG_DEBUG_DEDENT;*/
return result;
};
#endif /* 0 */
#endif
#endif /* WITH_OVERLUMPING */
inline
......@@ -1297,7 +1399,7 @@ geno_matrix ancestor_matrix(char a)
#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
};
}
......@@ -1311,11 +1413,11 @@ inline
geno_matrix lump_using_partition_weighted(const geno_matrix& m, const std::set<subset>& lumping_partition)
{
m.check_not_nan();
MSG_DEBUG("|| LUMPING");
MSG_DEBUG("" << m);
MSG_DEBUG("|| WITH PARTITION");
MSG_DEBUG("" << lumping_partition);
MSG_QUEUE_FLUSH();
/*MSG_DEBUG("|| LUMPING");*/
/*MSG_DEBUG("" << m);*/
/*MSG_DEBUG("|| WITH PARTITION");*/
/*MSG_DEBUG("" << lumping_partition);*/
/*MSG_QUEUE_FLUSH();*/
geno_matrix ret_lump;
......@@ -1396,17 +1498,17 @@ geno_matrix lump_using_partition_weighted(const geno_matrix& m, const std::set<s
{
std::set<symmetry_table_type> tmp_sym;
MatrixXb shrink = L1.cast<bool>();
MSG_DEBUG("SHRINK . SHRINK:transpose");
MSG_DEBUG((shrink * shrink.transpose()));
/*MSG_DEBUG("SHRINK . SHRINK:transpose");*/
/*MSG_DEBUG((shrink * shrink.transpose()));*/
for (const auto& S: m.symmetries) {
/*tmp_sym.emplace(ret_lump.labels, shrink * S.matrix() * shrink.transpose());*/
MSG_DEBUG("SHRINKING SYMMETRY");
MSG_DEBUG(S.matrix());
MSG_DEBUG("WITH");
MSG_DEBUG(shrink);
/*MSG_DEBUG("SHRINKING SYMMETRY");*/
/*MSG_DEBUG(S.matrix());*/
/*MSG_DEBUG("WITH");*/
/*MSG_DEBUG(shrink);*/
auto it_ok = tmp_sym.emplace(shrink * S.matrix() * shrink.transpose(), S.letters);
MSG_DEBUG("RESULT");
MSG_DEBUG(it_ok.first->matrix());
/*MSG_DEBUG("RESULT");*/
/*MSG_DEBUG(it_ok.first->matrix());*/
}
ret_lump.symmetries.assign(tmp_sym.begin(), tmp_sym.end());
}
......@@ -1421,10 +1523,10 @@ geno_matrix lump_using_partition_weighted(const geno_matrix& m, const std::set<s
MatrixXb ret = latent_lump(shrink, S.matrix());
auto it_ok = tmp_sym.emplace(ret, S.letters);
MSG_DEBUG("LUMPING LATENT SYMMETRY");
MSG_DEBUG(S.dump(m.labels, true));
MSG_DEBUG("INTO");
MSG_DEBUG(it_ok.first->dump(ret_lump.labels, true));
/*MSG_DEBUG("LUMPING LATENT SYMMETRY");*/
/*MSG_DEBUG(S.dump(m.labels, true));*/
/*MSG_DEBUG("INTO");*/
/*MSG_DEBUG(it_ok.first->dump(ret_lump.labels, true));*/
/*tmp_sym.emplace(shrink * S.matrix() * shrink.transpose(), S.letters);*/
/*tmp_sym.emplace(S);*/
......@@ -1567,7 +1669,7 @@ std::ostream& operator << (std::ostream& os, const std::forward_list<T>& l)
#endif
}
experimental_lumper(const geno_matrix& m, const symmetry_list_type& sym_table)
experimental_lumper(const geno_matrix& m, const symmetry_list_type& /*sym_table*/)
: M(m)
{
/*complete_symmetries();*/
......@@ -1715,7 +1817,7 @@ std::ostream& operator << (std::ostream& os, const std::forward_list<T>& l)
std::map<double, subset> keys;
double cost_;
conflict_type(const experimental_lumper& el, const std::set<subset>& P0,
conflict_type(const experimental_lumper& /*el*/, const std::set<subset>& /*P0*/,
const subset& c, const subset& b, const std::map<double, subset>& k)
: C(c), B(b), keys(k)
/*, cost_(opportunity_cost(el, P0))*/
......@@ -1836,6 +1938,7 @@ std::ostream& operator << (std::ostream& os, const std::forward_list<T>& l)
}
#endif
return ret;
(void) el; (void) P0;
}
};
......@@ -2106,6 +2209,7 @@ std::ostream& operator << (std::ostream& os, const std::forward_list<T>& l)
}
}
(void) splitted;
#if 0
/*for (int s1: splitted) {*/
/*for (const auto& split: splits) {*/
......
This diff is collapsed.
#ifndef _SPEL_PERMUT_H
#define _SPEL_PERMUT_H
#include <map>
#include <iostream>
#include "eigen.h"
struct permutation_type;
struct transposed_permutation_type;
template <typename TYPE> struct transposition_result;
template <> struct transposition_result<permutation_type> {
typedef transposed_permutation_type type;
typedef transposed_permutation_type return_type;
static return_type result(const permutation_type* p);
};
template <> struct transposition_result<transposed_permutation_type> {
typedef permutation_type type;
typedef const permutation_type& return_type;
static return_type result(const transposed_permutation_type* p);
};
template <typename DERIVED>
struct permutation_base {
virtual const std::vector<int>& map() const = 0;
virtual const std::vector<int>& transposed_map() const = 0;
typename transposition_result<DERIVED>::return_type transpose() const;
size_t size() const { return map().size(); }
template <typename OTHER_DERIVED>
permutation_type operator * (const permutation_base<OTHER_DERIVED>& other);
const DERIVED* derived() const { return dynamic_cast<const DERIVED*>(this); }
template <typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>
operator % (const Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& other)
{
Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> ret(other.rows(), other.cols());
const std::vector<int>& m_map = map();
const std::vector<int>& m_transposed = transposed_map();
for (int i = 0; i < other.rows(); ++i) {
for (int j = 0; j < other.cols(); ++j) {
ret(m_map[i], m_transposed[j]) = other(i, j);
}
}
return ret;
}
template <typename VALUE_TYPE>
std::vector<VALUE_TYPE>
operator * (const std::vector<VALUE_TYPE>& vec) const
{
std::vector<VALUE_TYPE> ret;
ret.reserve(vec.size());
const std::vector<int>& m_transposed = transposed_map();
for (size_t i = 0; i < vec.size(); ++i) {
ret.push_back(vec[m_transposed[i]]);
}
return ret;
}
template <typename Scalar>
Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> matrix() const
{
const std::vector<int>& m_map = map();
Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> ret = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>::Zero(m_map.size(), m_map.size());
for (size_t i = 0; i < m_map.size(); ++i) {
ret(m_map[i], i) = 1;
}
return ret;
}
template <typename OTHER_DERIVED>
bool operator == (const permutation_base<OTHER_DERIVED>& other) const
{
return map() == other.map();
}
template <typename OTHER_DERIVED>
bool operator < (const permutation_base<OTHER_DERIVED>& other) const
{
return map() < other.map();
}
friend
std::ostream& operator << (std::ostream& os, const permutation_base<DERIVED>& perm)
{
return os << ('.' + ('@' - '.') * perm.matrix<char>().array());
}
};
struct permutation_type : public permutation_base<permutation_type> {
std::vector<int> m_map;
std::vector<int> m_transposed;
static
permutation_type identity(int n)
{
permutation_type ret;
ret.m_map.reserve(n);
ret.m_transposed.reserve(n);
for (int i = 0; i < n; ++i) {
ret.m_map.push_back(i);
ret.m_transposed.push_back(i);
}
return ret;
}
permutation_type() : m_map(), m_transposed() {}
permutation_type(const permutation_type& p) : m_map(p.m_map), m_transposed(p.m_transposed) {}
permutation_type(permutation_type&& p) : m_map(std::move(p.m_map)), m_transposed(std::move(p.m_transposed)) {}
permutation_type(const std::vector<int>& m) : m_map(m), m_transposed() { compute_transposed(); }
permutation_type(std::vector<int>&& m) : m_map(std::move(m)), m_transposed() { compute_transposed(); }
permutation_type(const std::vector<int>& m, const std::vector<int>& t) : m_map(m), m_transposed(t) {}
template <typename Scalar>
permutation_type(const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>& mat)
: m_map(mat.rows(), -1), m_transposed(mat.rows(), -1)
{
assert(mat.rows() == mat.cols());
for (int col = 0; col < mat.cols(); ++col) {
int row;
mat.col(col).maxCoeff(&row);
m_map[col] = row;
m_transposed[row] = col;
}
}
const std::vector<int>& map() const { return m_map; }
const std::vector<int>& transposed_map() const { return m_transposed; }
template <typename OTHER_DERIVED>
permutation_type& operator = (const permutation_base<OTHER_DERIVED>& p) { m_map = p.map(); m_transposed = p.transposed_map(); return *this; }
permutation_type& operator = (permutation_type&& p) { m_map.swap(p.m_map); m_transposed.swap(p.m_transposed); return *this; }
permutation_type& compute_transposed()
{
m_transposed.clear();
m_transposed.resize(m_map.size());
for (size_t i = 0; i < m_map.size(); ++i) {
m_transposed[m_map[i]] = i;
}
return *this;
}
};
struct transposed_permutation_type : public permutation_base<transposed_permutation_type> {
const permutation_type* ptr;
transposed_permutation_type(const permutation_type* viewed) : ptr(viewed) {}
const std::vector<int>& map() const { return ptr->m_transposed; }
const std::vector<int>& transposed_map() const { return ptr->m_map; }
};
/*static*/
typename transposition_result<permutation_type>::return_type
transposition_result<permutation_type>::
result(const permutation_type* p) { return {p}; }
/*static*/
typename transposition_result<transposed_permutation_type>::return_type
transposition_result<transposed_permutation_type>::
result(const transposed_permutation_type* p) { return *p->ptr; }
template <typename DERIVED>
typename transposition_result<DERIVED>::return_type
permutation_base<DERIVED>::
transpose() const
{
return transposition_result<DERIVED>::result(derived());
}
template <typename DERIVED>
template <typename OTHER_DERIVED>
permutation_type
permutation_base<DERIVED>::
operator * (const permutation_base<OTHER_DERIVED>& other)
{
permutation_type ret;
assert(size() == other.size());
ret.m_map.resize(size(), -1);
ret.m_transposed.resize(size(), -1);
const std::vector<int>& m_map = map();
const std::vector<int>& m_transposed = transposed_map();
const std::vector<int>& other_map = other.map();
const std::vector<int>& other_transposed = other.transposed_map();
for (size_t i = 0; i < other.size(); ++i) {
ret.m_map[i] = m_map[other_map[i]];
ret.m_transposed[i] = other_transposed[m_transposed[i]];
}
#ifndef NDEBUG
for (size_t i = 0; i < other.size(); ++i) {
assert(ret.m_transposed[m_map[other_map[i]]] == i);
}
#endif
return ret;
/*return ret.compute_transposed();*/
}
#endif
......@@ -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
......@@ -164,5 +164,8 @@ test_bayes_dai: ../sandbox/bayes.cc bayes.d $(ALL_BUT_MAIN_OBJ)
test_sib: pedigree_analysis.cc static_data.o
$C $< static_data.o -o $@
matrixexp: geno_matrix.cc ../include/lumping2.h ../include/geno_matrix.h Makefile
matrixexp: geno_matrix.cc ../include/lumping2.h ../include/geno_matrix.h ../include/pedigree.h Makefile
$C -DORDER=$(ORDER) $< -o $@ -DTEST_GENO_MATRIX
test_permutation: test_permutation.cc ../include/permutation.h Makefile
$C -DORDER=$(ORDER) $< -o $@
......@@ -545,6 +545,88 @@ int main(int argc, char** argv)
size_t Y = ped.selfing(X); (void) Y;
}
/* 65536 */
if (test_case())
{
geno_matrix A = ancestor_matrix('a');
geno_matrix B = ancestor_matrix('b');
geno_matrix F1 = lump(A|B);
geno_matrix F2_1 = lump(F1|F1);
geno_matrix F2_2 = F2_1;
geno_matrix ones_4;
ones_4.variant = Geno;
ones_4.stat_dist = (VectorXd(4) << .25, .25, .25, .25).finished();
ones_4.inf_mat = MatrixXd::Ones(4, 4);
ones_4.labels = F2_1.labels;
ones_4.collect = ones_4.dispatch = MatrixXd::Identity(4, 4);
ones_4 = lump_using_partition_weighted(ones_4, {{0}, {1}, {2}, {3}});
MSG_DEBUG("" << ones_4);
geno_matrix F3_1 = lump(F2_1 | ones_4);
geno_matrix F3_2 = lump(F2_2 | ones_4);
}
/* 131072 */
if (test_case())
{
pedigree_type ped(30);
size_t A = ped.ancestor();
size_t B = ped.ancestor();
size_t C = ped.ancestor();
size_t D = ped.ancestor();
size_t AB = ped.crossing(A, B);
size_t CD = ped.crossing(C, D);
size_t F2_1 = ped.crossing(AB, CD);
size_t F2_2 = ped.crossing(AB, CD);
size_t F2_rev = ped.crossing(CD, AB);
size_t F3_1 = ped.crossing(F2_1, F2_2);
size_t F3_rev = ped.crossing(F2_1, F2_rev);
auto g1 = ped.get_gen(F3_1);
auto g2 = ped.get_gen(F3_rev);
/*size_t sz = ped.get_gen(F2_1)->cols() * 2;*/
MatrixXb prout_1 = MatrixXb::Identity(8, 8);
MatrixXb prout_2 = generate_lozenge(2, 2);
MatrixXb prout_3 = MatrixXb::Ones(2, 2) - MatrixXb::Identity(2, 2);
MatrixXb sub = kroneckerProduct(prout_1, prout_2);
MatrixXb lz = kroneckerProduct(sub, prout_3);
MSG_DEBUG("--");
MSG_DEBUG("Compare g1->inf_mat and g2->inf_mat: " << std::boolalpha << (g1->inf_mat - g2->inf_mat).isZero());
/*MSG_DEBUG((g1->inf_mat - lz * g2->inf_mat * lz.transpose()));*/
MSG_DEBUG("Compare g1->collect and g2->collect: " << std::boolalpha << (g1->collect - g2->collect).isZero());
lz = g1->collect.cast<bool>() * lz * g1->collect.cast<bool>().transpose();
MatrixXd im = lz.cast<double>().transpose() * g2->inf_mat * lz.cast<double>();
MSG_DEBUG("Compare g1->inf_mat and lz o g2->inf_mat: " << std::boolalpha << (g1->inf_mat - im).isZero());
MSG_DEBUG(('.' + (('@' - '.') * lz.array().cast<char>())));
std::vector<label_type> g2l(g2->labels.size());
for (int i = 0; i < lz.cols(); ++i) {
int r;
lz.col(i).maxCoeff(&r);
g2l[r] = g2->labels[i];
}
MSG_DEBUG("Compare labels: " << std::boolalpha << (g2l == g1->labels));
MSG_DEBUG("Compare (wrong) labels: " << std::boolalpha << (g2->labels == g1->labels));
}
/* 262144 */
if (test_case())
{
pedigree_type ped;
size_t A = ped.ancestor();
size_t B = ped.ancestor();
size_t C = ped.ancestor();
size_t D = ped.ancestor();
size_t E = ped.ancestor();
size_t F = ped.ancestor();
size_t AB = ped.crossing(A, B);
size_t CD = ped.crossing(C, D);
size_t EF = ped.crossing(E, F);
size_t BCa = ped.crossing(A, AB);
size_t BCc = ped.crossing(C, CD);
size_t BCf = ped.crossing(EF, F);
size_t F1_1 = ped.crossing(BCa, BCc);
size_t F1_2 = ped.crossing(BCc, BCf);
size_t F2 = ped.crossing(F1_1, F1_2); (void) F2;
}
return 0;
(void) argc; (void) argv;
}
......
#include "permutation.h"
#include <iomanip>
int main(int argc, char** argv)
{
permutation_type p0 = permutation_type::identity(10);
std::cout << "p0" << std::endl << p0 << std::endl << std::endl;
std::vector<int> perm = {0, 2, 4, 1, 3, 5, 6, 8, 7, 9};
permutation_type p1(perm);
std::cout << "p1" << std::endl << p1 << std::endl << std::endl;
std::cout << "p1.transpose" << std::endl << (p1.transpose()) << std::endl << std::endl;
std::cout << "p1 * p1" << std::endl << (p1 * p1) << std::endl << std::endl;
std::cout << "p1 * p1.transpose" << std::endl << (p1 * p1.transpose()) << std::endl << std::endl;
std::vector<std::string> oups = {"pouet", "plop", "coin", "prout", "plof", "hop", "bouh", "chou", "genou", "caillou"};