Commit 554cc32e authored by Damien Leroux's avatar Damien Leroux
Browse files

New implementation for symmetries. Promising. Almost complete.

parent 26a28c6f
This diff is collapsed.
......@@ -24,6 +24,13 @@ extern "C" {
#define _CYAN "\x1b[36;1m"
#define _NORMAL "\x1b[0;m"
#define _DBG_INDENT_MARK '\x11'
#define _DBG_DEDENT_MARK '\x12'
#define _DBG_INDENT_MARK_S "\x11"
#define _DBG_DEDENT_MARK_S "\x12"
#define _DBG_SYNC_MARK '\x13'
#define _DBG_SYNC_MARK_S "\x13"
#define MSG_HANDLER_IS_SYNCED
enum msg_channel { Out, Err, Log };
......@@ -281,7 +288,7 @@ inline
int ostream_manager::HeaderInserter::overflow(int ch)
{
int retval = 0;
if (ch == 1) {
if (ch == _DBG_INDENT_MARK) {
/*indent += 3;*/
prefix_insert_mode = !prefix_insert_mode;
if (prefix_insert_mode) {
......@@ -290,9 +297,9 @@ int ostream_manager::HeaderInserter::overflow(int ch)
} else {
prefix.push_back(prefix_ss.str());
}
} else if (ch == 2) {
} else if (ch == _DBG_DEDENT_MARK) {
if (prefix_insert_mode) {
overflow('\x1');
overflow(_DBG_INDENT_MARK);
}
if (prefix.size()) {
prefix.pop_back();
......@@ -340,7 +347,7 @@ void message_queue::run()
m_queue.pop_front();
}
if (next->message == "\x3") {
if (next->message == _DBG_SYNC_MARK_S) {
m_flush_condition.notify_one();
continue;
} else if (next->message.size() == 0) {
......@@ -379,7 +386,7 @@ void message_queue::run()
} while (0)
#define MSG_QUEUE_FLUSH() do { \
CREATE_MESSAGE(msg_channel::Log, "\x3"); \
CREATE_MESSAGE(msg_channel::Log, _DBG_SYNC_MARK_S); \
msg_handler_t::wait_for_flush(); \
} while (0)
......@@ -415,9 +422,9 @@ void message_queue::run()
} while(0)
#endif
#define MSG_DEBUG_INDENT CREATE_MESSAGE(msg_channel::Log, "\x1 \x1")
#define MSG_DEBUG_INDENT_EXPR(_str_) CREATE_MESSAGE(msg_channel::Log, MESSAGE('\x1' << _str_ << '\x1'))
#define MSG_DEBUG_DEDENT CREATE_MESSAGE(msg_channel::Log, "\x2")
#define MSG_DEBUG_INDENT CREATE_MESSAGE(msg_channel::Log, _DBG_INDENT_MARK_S " " _DBG_INDENT_MARK_S)
#define MSG_DEBUG_INDENT_EXPR(_str_) CREATE_MESSAGE(msg_channel::Log, MESSAGE(_DBG_INDENT_MARK << _str_ << _DBG_INDENT_MARK))
#define MSG_DEBUG_DEDENT CREATE_MESSAGE(msg_channel::Log, _DBG_DEDENT_MARK_S)
struct scoped_indent {
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
#ifndef _SPEL_PERMUT_H
#define _SPEL_PERMUT_H
#ifndef _SPEL_PERMUT_H_
#define _SPEL_PERMUT_H_
#include <map>
#include <iostream>
#include "eigen.h"
#include "braille_plot.h"
#include "lumping2.h"
struct permutation_type;
struct transposed_permutation_type;
......@@ -47,22 +48,27 @@ struct permutation_base {
typename transposition_result<DERIVED>::return_type transpose() const;
size_t size() const { return map().size(); }
size_t cols() const { return size(); }
size_t rows() const { return size(); }
size_t innerSize() const { return size(); }
size_t outerSize() const { return size(); }
template <typename OTHER_DERIVED>
permutation_type operator * (const permutation_base<OTHER_DERIVED>& other);
permutation_type operator * (const permutation_base<OTHER_DERIVED>& other) const;
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)
operator % (const Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& other) const
{
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);
/*ret(m_transposed[i], m_map[j]) = other(i, j);*/
ret(i, m_map[j]) = other(m_transposed[i], j);
}
}
return ret;
......@@ -98,6 +104,12 @@ struct permutation_base {
return map() == other.map();
}
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
{
......@@ -109,6 +121,7 @@ struct permutation_base {
{
/*return os << ('.' + ('@' - '.') * perm.matrix<char>().array());*/
braille_grid grid(perm.size(), perm.size());
grid.set_background(32, 20, 32);
const auto& map = perm.map();
for (size_t i = 0; i < map.size(); ++i) {
grid.put_pixel(i, map[i]);
......@@ -208,11 +221,14 @@ struct permutation_type : public permutation_base<permutation_type> {
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(); }
explicit permutation_type(const std::vector<int>& m) : m_map(m), m_transposed() { compute_transposed(); }
explicit permutation_type(std::vector<int>&& m) : m_map(std::move(m)), m_transposed() { compute_transposed(); }
explicit permutation_type(std::initializer_list<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 DERIVED>
permutation_type(const permutation_base<DERIVED>& p) : m_map(p.map()), m_transposed(p.transposed_map()) {}
template <typename Scalar>
permutation_type(const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>& mat)
explicit 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());
......@@ -250,6 +266,51 @@ struct permutation_type : public permutation_base<permutation_type> {
}
return *this;
}
permutation_type lump(const std::set<subset>& P0) const
{
std::map<int, int> bins;
int b = 0;
for (const auto& S: P0) {
for (int s: S) {
bins[s] = b;
}
++b;
}
/* check consistency */
for (const auto& S: P0) {
auto i = S.begin(), j = S.end();
b = bins[m_map[*i]];
for (++i; i != j; ++i) {
if (b != bins[m_map[*i]]) {
/* whine */
return {};
}
}
}
/* perform lump */
permutation_type ret;
for (const auto& S: P0) {
ret.m_map.push_back(bins[m_map[*S.begin()]]);
}
ret.compute_transposed();
return ret;
}
permutation_type lump(const Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic>& mat) const
{
std::set<subset> P;
for (int r = 0; r < mat.rows(); ++r) {
subset S;
for (int c = 0; c < mat.cols(); ++c) {
if (mat(r, c)) {
S.push_back(c);
}
}
P.insert(S);
}
return lump(P);
}
};
......@@ -289,9 +350,16 @@ template <typename DERIVED>
template <typename OTHER_DERIVED>
permutation_type
permutation_base<DERIVED>::
operator * (const permutation_base<OTHER_DERIVED>& other)
operator * (const permutation_base<OTHER_DERIVED>& other) const
{
permutation_type ret;
/*MSG_DEBUG_INDENT_EXPR("[product A] ");*/
/*MSG_DEBUG((*this));*/
/*MSG_DEBUG_DEDENT;*/
/*MSG_DEBUG_INDENT_EXPR("[product B] ");*/
/*MSG_DEBUG(other);*/
/*MSG_DEBUG_DEDENT;*/
/*MSG_QUEUE_FLUSH();*/
assert(size() == other.size());
ret.m_map.resize(size(), -1);
ret.m_transposed.resize(size(), -1);
......@@ -306,16 +374,18 @@ template <typename OTHER_DERIVED>
#ifndef NDEBUG
for (size_t i = 0; i < other.size(); ++i) {
assert(ret.m_transposed[m_map[other_map[i]]] == i);
assert(ret.m_transposed[m_map[other_map[i]]] == (int) i);
}
#endif
/*MSG_DEBUG_INDENT_EXPR("[product R] ");*/
/*MSG_DEBUG(ret);*/
/*MSG_DEBUG_DEDENT;*/
return ret;
/*return ret.compute_transposed();*/
}
#endif
template <typename DERIVED1, typename DERIVED2>
permutation_type
kronecker(const permutation_base<DERIVED1>& p1, const permutation_base<DERIVED2>& p2)
......@@ -336,3 +406,67 @@ template <typename DERIVED1, typename DERIVED2>
}
return ret;
}
#if 0
inline
permutation_type
kronecker(const std::vector<permutation_type>& pvec, const std::vector<int>& refs)
{
std::vector<size_t> iterators(pvec.size(), 0);
std::vector<size_t> sizes(pvec.size(), 1);
for (size_t i = 0; i < pvec.size(); ++i) {
if (refs[i] != -1) {
sizes[i] = pvec[i].size();
}
}
auto get_i =
[&] (size_t i)
{
if (refs[i] == -1) {
return iterators[i];
} else {
return 0;
}
};
auto get_pi =
[&] (size_t i)
{
if (refs[i] == -1) {
return pvec[i].m_map[iterators[i]];
} else {
return pvec[i].m_map[iterators[refs[i]]];
}
};
auto next =
[&] ()
{
size_t i = 0;
while (i < iterators.size()) {
if (refs[i] == -1) {
++iterators[i];
if (iterators[i] == sizes[i]) {
iterators[i] = 0;
++i;
} else {
break;
}
} else {
++i;
}
}
return i != iterators.size();
};
auto get_col =
[&] ()
{
};
}
#endif
#endif
......@@ -164,8 +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 ../include/pedigree.h Makefile
matrixexp: geno_matrix.cc ../include/lumping2.h ../include/geno_matrix.h ../include/pedigree.h ../include/permutation.h ../include/braille_plot.h Makefile ../include/pedigree_tree.h
$C -DORDER=$(ORDER) $< -o $@ -DTEST_GENO_MATRIX
test_permutation: test_permutation.cc ../include/permutation.h Makefile
test_permutation: test_permutation.cc ../include/permutation.h ../include/braille_plot.h ../include/pedigree_tree.h Makefile
$C -DORDER=$(ORDER) $< -o $@
This diff is collapsed.
This diff is collapsed.
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