Commit 26a28c6f authored by Damien Leroux's avatar Damien Leroux
Browse files

Integration in progress.

parent a8208909
......@@ -5,45 +5,33 @@
#include "eigen.h"
struct braille_plot {
struct braille_grid {
int m_width;
int m_height;
typedef Eigen::Matrix<unsigned long, Eigen::Dynamic, Eigen::Dynamic> plot_data_type;
plot_data_type m_data;
double m_xmin;
double m_xmax;
double m_ymin;
double m_ymax;
double m_xrange_inv;
double m_yrange_inv;
braille_plot(int w, int h, double min_x, double max_x, double min_y, double max_y)
: m_width((w << 1)), m_height((h << 2)), m_data(plot_data_type::Constant(h + 1, w + 1, 0x0080A0E2)),
m_xmin(min_x), m_xmax(max_x), m_ymin(min_y), m_ymax(max_y),
m_xrange_inv(1. / (m_xmax - m_xmin)),
m_yrange_inv(1. / (m_ymax - m_ymin))
braille_grid(int dot_width, int dot_height)
: m_width(dot_width), m_height(dot_height)
, m_data(plot_data_type::Constant((dot_height + 3) >> 2, (dot_width + 1) >> 1, 0x0080A0E2))
{}
braille_plot&
put_pixel(double x, double y)
braille_grid&
put_pixel(int ix, int iy)
{
int ix = x2plot(x);
int iy = y2plot(y);
unsigned long bit = xy2bit(ix, iy);
int r = y2data(iy);
int c = x2data(ix);
if (r >= 0 && r < m_data.rows() && c >= 0 && c < m_data.cols()) {
m_data(r, c) |= bit;
/*std::cout << "put_pixel(" << std::dec << x << ", " << y << ") ix=" << ix << " iy=" << iy << " bit=" << std::hex << bit << std::dec << " r=" << r << " c=" << c << " data=" << std::hex << m_data(r, c) << std::endl;*/
/*std::cout << "put_pixel(" << std::dec << ix << ", " << iy << ") bit=" << std::hex << bit << std::dec << " r=" << r << " c=" << c << " data=" << std::hex << m_data(r, c) << std::endl;*/
}
return *this;
}
braille_plot&
clear_pixel(double x, double y)
braille_grid&
clear_pixel(int ix, int iy)
{
int ix = x2plot(x);
int iy = y2plot(y);
unsigned char bit = xy2bit(ix, iy);
int r = y2data(iy);
int c = x2data(ix);
......@@ -54,10 +42,8 @@ struct braille_plot {
}
bool
get_pixel(double x, double y) const
get_pixel(int ix, int iy) const
{
int ix = x2plot(x);
int iy = y2plot(y);
unsigned char bit = xy2bit(ix, iy);
int r = y2data(iy);
int c = x2data(ix);
......@@ -67,30 +53,9 @@ struct braille_plot {
return false;
}
template <typename FUN>
braille_plot&
plot(FUN fun, double step)
{
for (double x = m_xmin; x < m_xmax; x += step) {
put_pixel(x, fun(x));
}
return *this;
}
friend
std::ostream& operator << (std::ostream& os, const braille_plot& plot)
std::ostream& operator << (std::ostream& os, const braille_grid& plot)
{
/*static Eigen::IOFormat fmt(Eigen::StreamPrecision, Eigen::DontAlignCols,*/
/*"", "\n",*/
/*"", "",*/
/*"", "");*/
/*return os << plot.m_data.format(fmt);*/
/*static Eigen::IOFormat fmt(Eigen::StreamPrecision, Eigen::DontAlignCols,*/
/*"\x28", "\n",*/
/*"\x28", "",*/
/*"", "");*/
/*return os << plot.m_data.format(fmt);*/
for (int r = 0; r < plot.m_data.rows(); ++r) {
for (int c = 0; c < plot.m_data.cols(); ++c) {
os.write((char*) &plot.m_data(r, c), 3);
......@@ -115,17 +80,181 @@ private:
return bitval[ ((y & 0x3) << 1) | (x & 1) ];
}
};
struct braille_plot : public braille_grid {
double m_xmin;
double m_xmax;
double m_ymin;
double m_ymax;
double m_xrange_inv;
double m_yrange_inv;
double m_xres;
double m_yres;
braille_plot(int w, int h, double min_x, double max_x, double min_y, double max_y)
: braille_grid(w << 1, h << 2),
/*m_width((w << 1)), m_height((h << 2)), m_data(plot_data_type::Constant(h + 1, w + 1, 0x0080A0E2)),*/
m_xmin(min_x), m_xmax(max_x), m_ymin(min_y), m_ymax(max_y),
m_xrange_inv(1. / (m_xmax - m_xmin)),
m_yrange_inv(1. / (m_ymax - m_ymin)),
m_xres((m_xmax - m_xmin) / (m_width - 1)),
m_yres((m_ymax - m_ymin) / (m_height - 1))
{
/*std::cout << "m_xres=" << m_xres << " m_yres=" << m_yres << std::endl;*/
}
template <typename FUN>
braille_plot&
plot(FUN fun, double step)
{
for (double x = m_xmin; x < m_xmax; x += step) {
put_pixel(x, fun(x));
}
return *this;
}
braille_plot&
line(double x0, double y0, double x1, double y1, int dash_length, int space_length)
{
int gx0 = x2grid(x0);
int gy0 = y2grid(y0);
int gx1 = x2grid(x1);
int gy1 = y2grid(y1);
int gdx = abs(gx0 - gx1);
int gdy = abs(gy0 - gy1);
double dx, dy;
int counter = 0;
space_length += dash_length;
double x = x0, y = y0;
if (gdx < gdy) {
dy = y0 > y1 ? -m_yres : m_yres;
dx = ((x0 > x1 ? -m_xres : m_xres) * gdx) / gdy;
/*std::cout << "dx=" << dx << " dy=" << dy << std::endl;*/
if (dy > 0) {
for (; y <= y1; y += dy, x += dx) {
if (counter++ < dash_length) {
put_pixel(x, y);
}
counter %= space_length;
}
} else {
for (; y >= y1; y += dy, x += dx) {
if (counter++ < dash_length) {
put_pixel(x, y);
}
counter %= space_length;
}
}
} else {
dx = x0 > x1 ? -m_xres : m_xres;
dy = ((y0 > y1 ? -m_yres : m_yres) * gdy) / gdx;
/*std::cout << "dx=" << dx << " dy=" << dy << std::endl;*/
if (dx > 0) {
for (; x <= x1; y += dy, x += dx) {
if (counter++ < dash_length) {
put_pixel(x, y);
}
counter %= space_length;
}
} else {
for (; x >= x1; y += dy, x += dx) {
if (counter++ < dash_length) {
put_pixel(x, y);
}
counter %= space_length;
}
}
}
return *this;
}
braille_plot&
hline(double y, int dash_length, int space_length)
{
return line(m_xmin, y, m_xmax, y, dash_length, space_length);
}
braille_plot&
vline(double x, int dash_length, int space_length)
{
return line(x, m_ymin, x, m_ymax, dash_length, space_length);
}
braille_plot&
haxis(double y, double tick_origin, double tick_major, double tick_minor)
{
std::vector<double> ticks_maj, ticks_min;
for (double x = tick_origin; x > m_xmin; x -= tick_major) { ticks_maj.push_back(x); }
for (double x = tick_origin; x < m_xmax; x += tick_major) { ticks_maj.push_back(x); }
for (double x = tick_origin; x > m_xmin; x -= tick_minor) { ticks_min.push_back(x); }
for (double x = tick_origin; x < m_xmax; x += tick_minor) { ticks_min.push_back(x); }
return haxis(y, ticks_maj, ticks_min);
}
braille_plot&
vaxis(double x, double tick_origin, double tick_major, double tick_minor)
{
std::vector<double> ticks_maj, ticks_min;
for (double y = tick_origin; y > m_ymin; y -= tick_major) { ticks_maj.push_back(y); }
for (double y = tick_origin; y < m_ymax; y += tick_major) { ticks_maj.push_back(y); }
for (double y = tick_origin; y > m_ymin; y -= tick_minor) { ticks_min.push_back(y); }
for (double y = tick_origin; y < m_ymax; y += tick_minor) { ticks_min.push_back(y); }
return vaxis(x, ticks_maj, ticks_min);
}
braille_plot&
vaxis(double x, const std::vector<double>& ticks_major, const std::vector<double>& ticks_minor)
{
vline(x, 1, 0);
for (double y: ticks_minor) {
put_pixel(x - m_xres, y);
put_pixel(x + m_xres, y);
}
for (double y: ticks_major) {
put_pixel(x - m_xres - m_xres, y);
put_pixel(x + m_xres, y);
put_pixel(x - m_xres, y);
put_pixel(x + m_xres + m_xres, y);
}
return *this;
}
braille_plot&
haxis(double y, const std::vector<double>& ticks_major, const std::vector<double>& ticks_minor)
{
hline(y, 1, 0);
for (double x: ticks_minor) {
put_pixel(x, y - m_yres);
put_pixel(x, y + m_yres);
}
for (double x: ticks_major) {
put_pixel(x, y - m_yres - m_yres);
put_pixel(x, y + m_yres);
put_pixel(x, y - m_yres);
put_pixel(x, y + m_yres + m_yres);
}
return *this;
}
braille_plot& put_pixel(double x, double y)
{
braille_grid::put_pixel(x2grid(x), y2grid(y));
return *this;
}
int x2plot(double x) const
private:
int x2grid(double x) const
{
double rel_x = (x - m_xmin) * m_xrange_inv;
return (int) round(m_width * rel_x);
return (int) round((m_width - 1) * rel_x);
}
int y2plot(double y) const
int y2grid(double y) const
{
double rel_y = 1 - (y - m_ymin) * m_yrange_inv;
return (int) round(m_height * rel_y);
return (int) round((m_height - 1) * rel_y);
}
};
......
......@@ -5,6 +5,7 @@
#include "error.h"
#include "labelled_matrix.h"
#include "lumping2.h"
#include "permutation.h"
#include <limits>
#include <cmath>
......@@ -190,7 +191,7 @@ std::ostream& operator << (std::ostream& os, const std::map<char, int>& L)
}
inline
MatrixXb generate_lozenge(const std::vector<size_t> sizes&, const std::vector<size_t>& idx_permut, const std::vector<bool>& reverse_gamete)
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);
......@@ -409,7 +410,8 @@ std::ostream& operator << (std::ostream& os, const std::map<char, int>& L)
struct symmetry_table_type {
//std::vector<std::pair<char, char>> switches; /* over letters */
std::map<int, int> table; /* over states */
//std::map<int, int> table; /* over states */
permutation_type table;
letter_permutation_type letters;
symmetry_table_type(/*const std::vector<std::pair<char, char>>& _s, */const std::map<int, int>& _t, const letter_permutation_type& _l)
......@@ -426,7 +428,8 @@ 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)
#if 0
symmetry_table_type(const permutation_type& state_matrix, const std::vector<int>& permut_col, const std::vector<int>& permut_row)
: table(), letters()
{
/*std::set<std::pair<char, char>> tmp_switches;*/
......@@ -444,44 +447,11 @@ std::ostream& operator << (std::ostream& os, const std::map<char, int>& L)
}
}
}
symmetry_table_type(/*const std::vector<label_type>& labels, */const MatrixXb& state_matrix, const letter_permutation_type& l)
: /*switches(),*/ table(), letters(l)
{
#if 0
MSG_DEBUG_INDENT_EXPR("[symmetry_table_type ctor] ");
MSG_DEBUG("from matrix");
MSG_DEBUG(state_matrix);
#endif
/*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[j] = i;
}
}
}
/*switches.assign(tmp_switches.begin(), tmp_switches.end());*/
/*MSG_DEBUG(switch_matrix());*/
#if 0
for (const auto& s: switches) {
MSG_DEBUG("switch " << s);
}
MSG_DEBUG("--");
for (const auto& t: table) {
MSG_DEBUG(t.first << " -> " << t.second);
}
MSG_DEBUG("---------------");
MSG_DEBUG_DEDENT;
#endif
}
symmetry_table_type(/*const std::vector<label_type>& labels, */const permutation_type& state_matrix, const letter_permutation_type& l)
: /*switches(),*/ table(state_matrix), letters(l)
{}
symmetry_table_type() : table() {}
symmetry_table_type& operator = (const symmetry_table_type& other) { table = other.table; letters = other.letters; return *this; }
......@@ -1309,6 +1279,28 @@ geno_matrix kronecker(const geno_matrix& m1, const geno_matrix& m2)
}
if (ret.variant == Geno) {
permutation_type loz = generate_lozenge(m1.cols(), m2.cols());
symmetry_table_type solange(loz);
if (solange.is_consistent(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));*/
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;*/
/*latent_tmp.emplace(MatrixXb::Identity(m1.cols() * m2.cols(), m1.cols() * m2.cols()));*/
}
#if 1
#else
std::vector<bool> iter = {false, true};
for (auto r0: iter) {
for (auto r1: iter) {
......@@ -1333,6 +1325,7 @@ geno_matrix kronecker(const geno_matrix& m1, const geno_matrix& m2)
/*latent_tmp.emplace(MatrixXb::Identity(m1.cols() * m2.cols(), m1.cols() * m2.cols()));*/
}
}
#endif
}
#if 0
MatrixXb loz = generate_lozenge(m1.cols(), m2.cols());
......
......@@ -4,6 +4,7 @@
#include <map>
#include <iostream>
#include "eigen.h"
#include "braille_plot.h"
struct permutation_type;
struct transposed_permutation_type;
......@@ -22,9 +23,27 @@ template <> struct transposition_result<transposed_permutation_type> {
template <typename DERIVED>
struct permutation_base {
typedef DERIVED derived_type;
virtual const std::vector<int>& map() const = 0;
virtual const std::vector<int>& transposed_map() const = 0;
struct const_iterator {
std::vector<int>::const_iterator m_i, m_beg;
typedef std::pair<int, int> value_type;
value_type tmp;
const_iterator(std::vector<int>::const_iterator m, std::vector<int>::const_iterator b) : m_i(m), m_beg(b), tmp() {}
bool operator == (const_iterator c) { return m_i == c.m_i; }
bool operator != (const_iterator c) { return m_i != c.m_i; }
bool operator < (const_iterator c) { return m_i < c.m_i; }
value_type operator * () const { return {m_i - m_beg, *m_i}; }
const_iterator& operator ++ () { ++m_i; return *this; }
const_iterator& operator ++ (int) { ++m_i; return *this; }
const value_type* operator -> () { tmp = **this; return &tmp; }
};
const_iterator begin() const { return {map().begin(), map().begin()}; }
const_iterator end() const { return {map().end(), map().begin()}; }
typename transposition_result<DERIVED>::return_type transpose() const;
size_t size() const { return map().size(); }
......@@ -88,8 +107,16 @@ struct permutation_base {
friend
std::ostream& operator << (std::ostream& os, const permutation_base<DERIVED>& perm)
{
return os << ('.' + ('@' - '.') * perm.matrix<char>().array());
/*return os << ('.' + ('@' - '.') * perm.matrix<char>().array());*/
braille_grid grid(perm.size(), perm.size());
const auto& map = perm.map();
for (size_t i = 0; i < map.size(); ++i) {
grid.put_pixel(i, map[i]);
}
return os << grid;
}
int operator [] (int state) const { return map()[state]; }
};
......@@ -110,6 +137,74 @@ struct permutation_type : public permutation_base<permutation_type> {
return ret;
}
static
permutation_type anti_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(n - 1 - i);
ret.m_transposed.push_back(n - 1 - i);
}
return ret;
}
static
permutation_type lozenge(size_t n0, size_t n1)
{
int N = n0 * n1;
permutation_type ret;
ret.m_map.resize(N);
ret.m_transposed.resize(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.m_map[j] = i;
ret.m_transposed[i] = j;
}
}
return ret;
}
template <typename DERIVED>
static
permutation_type lozenge(size_t n0, size_t n1, const permutation_base<DERIVED>& 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());*/
const std::vector<int>& permut_col = permut.transposed_map();
const std::vector<int>& permut_row = permut.map();
permutation_type ret;
ret.m_map.resize(N);
ret.m_transposed.resize(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;*/
ret.m_map[permut_col[j]] = permut_row[i];
ret.m_transposed[permut_row[i]] = permut_col[j];
}
}
/*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;
}
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)) {}
......@@ -129,12 +224,21 @@ struct permutation_type : public permutation_base<permutation_type> {
m_transposed[row] = col;
}
}
template <typename Alloc>
permutation_type(const std::map<int, int, Alloc>& map)
: m_map(map.size(), -1), m_transposed(map.size(), -1)
{
for (const auto& kv: map) {
m_map[kv.first] = kv.second;
m_transposed[kv.second] = kv.first;
}
}
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 = (const permutation_type& p) { m_map = p.map(); m_transposed = p.transposed_map(); return *this; }
permutation_type& operator = (const transposed_permutation_type& p);
permutation_type& operator = (permutation_type&& p) { m_map.swap(p.m_map); m_transposed.swap(p.m_transposed); return *this; }
permutation_type& compute_transposed()
......@@ -158,6 +262,7 @@ struct transposed_permutation_type : public permutation_base<transposed_permutat
const std::vector<int>& transposed_map() const { return ptr->m_map; }
};
permutation_type& permutation_type::operator = (const transposed_permutation_type& p) { m_map = p.map(); m_transposed = p.transposed_map(); return *this; }
......@@ -211,3 +316,23 @@ template <typename OTHER_DERIVED>
#endif
template <typename DERIVED1, typename DERIVED2>
permutation_type
kronecker(const permutation_base<DERIVED1>& p1, const permutation_base<DERIVED2>& p2)
{
permutation_type ret;
size_t sz1 = p1.size();
size_t sz2 = p2.size();
auto compute_state = [&] (int st1, int st2) { return st1 * sz2 + st2; };
ret.m_map.resize(p1.size() * sz2);
ret.m_transposed.resize(sz1 * sz2);
size_t s = 0;
for (size_t s1 = 0; s1 < sz1; ++s1) {
for (size_t s2 = 0; s2 < sz2; ++s2, ++s) {
size_t t = compute_state(p1[s1], p2[s2]);
ret.m_map[s] = t;
ret.m_transposed[t] = s;
}
}
return ret;
}
#include "permutation.h"
#include "geno_matrix.h"
#include <iomanip>
......@@ -36,6 +36,129 @@ int main(int argc, char** argv)
std::cout << (p2 % mat) << std::endl;
permutation_type p3({2, 1, 0});
std::cout << "kronecker of" << std::endl << p2 << std::endl << "and" << std::endl << p3 << std::endl << "=" << std::endl;
std::cout << kronecker(p2, p3) << std::endl;
std::cout << "kronecker of" << std::endl << p3 << std::endl << "and" << std::endl << p2 << std::endl << "=" << std::endl;
std::cout << kronecker(p3, p2) << std::endl;
{