Skip to content
GitLab
Commit 4c1779a4 authored by Damien Leroux's avatar Damien Leroux
Browse files

Made lazy evaluations not lazy anymore.

parent 94ec6fe1
Hide whitespace changes
Inline Side-by-side
include/model.h
View file @ 4c1779a4
......@@ -2,6 +2,7 @@
#define _SPEL_MODEL_H_
#include <iostream>
#include "cache2.h"
#include "model/fwd.h"
#include "model/model.h"
#include "model/model_element.h"
......
include/model/model.h
View file @ 4c1779a4
#ifndef _SPEL_MODEL_MODEL_H_
#define _SPEL_MODEL_MODEL_H_
#include <stdexcept>
#include <Eigen/SVD>
#include <Eigen/QR>
......@@ -11,6 +13,72 @@
#define COMPONENT_EPSILON (1.e-10)
#include "labelled_matrix.h"
typedef labelled_matrix<Eigen::Matrix<double, -1, -1>, int, std::vector<char>> model_block_type;
namespace std {
template <typename T>
struct hash<std::vector<T>> {
size_t operator () (const std::vector<T>& v)
{
hash<T> h;
size_t accum = 0;
for (auto& x: v) {
accum ^= h(x);
}
return accum;
}
};
template <typename _Scalar, int A, int B, int C, int D, int E>
struct hash<Eigen::Matrix<_Scalar, A, B, C, D, E>> {
struct red_mat {
size_t accum;
void init(_Scalar s, int i, int j)
{
accum = hash<_Scalar>()(s);
(void)i; (void)j;
}
void operator () (_Scalar s, int i, int j)
{
accum = impl::ROTATE<7>(accum) ^ hash<_Scalar>()(s);
(void)i; (void)j;
}
};
size_t operator () (const value<Eigen::Matrix<_Scalar, A, B, C, D, E>>& m) const
{
red_mat rm;
m->visit(rm);
return rm.accum;
}
size_t operator () (const Eigen::Matrix<_Scalar, A, B, C, D, E>& m) const
{
red_mat rm;
m.visit(rm);
return rm.accum;
}
};
template <typename M, typename R, typename C>
struct hash<labelled_matrix<M, R, C>> {
size_t operator () (const labelled_matrix<M, R, C>& m) const
{
hash<M> hm;
hash<R> hr;
hash<C> hc;
size_t accum = hm(m.data);
for (auto& r: m.row_labels) { accum ^= hr(r); }
for (auto& c: m.column_labels) { accum ^= hc(c); }
return accum;
}
size_t operator () (const value<labelled_matrix<M, R, C>>& m) const
{
return operator () (*m);
}
};
}
static inline
constexpr bool
around_zero(double o)
......@@ -59,7 +127,7 @@ MatrixXd concat_right(const std::vector<MatrixXd>& mat_vec)
}
static inline
MatrixXd concat_right(const std::vector<const MatrixXd*>& mat_vec)
MatrixXd concat_right(const collection<model_block_type>& mat_vec)
{
size_t full_size = 0;
MatrixXd ret;
......@@ -72,7 +140,7 @@ MatrixXd concat_right(const std::vector<const MatrixXd*>& mat_vec)
for (auto m: mat_vec) {
/*MSG_DEBUG("concat_right in M(" << ret.innerSize() << ',' << ret.outerSize() << ") at col " << full_size << "matrix(" << m->innerSize() << ',' << m->outerSize() << ')');*/
/*ret.block(0, full_size, ret.innerSize(), m->outerSize()) = *m;*/
ret.middleCols(full_size, m->outerSize()) = *m;
ret.middleCols(full_size, m->outerSize()) = m->data;
full_size += m->outerSize();
}
return ret;
......@@ -142,45 +210,211 @@ MatrixXd components(const MatrixXd& M, const MatrixXd& P)
enum class SolverType { QR, SVD };
namespace std {
template <typename _Scalar, int A, int B, int C, int D, int E>
struct hash<Eigen::Matrix<_Scalar, A, B, C, D, E>> {
struct red_mat {
size_t accum;
void init(_Scalar s, int i, int j)
{
accum = hash<_Scalar>()(s);
(void)i; (void)j;
struct model {
model()
: m_Y(), m_blocs(), m_X(), m_rank(), m_rss(), m_coefficients(), m_solver_type(), m_computed(false)
{}
model(const value<MatrixXd>& y, SolverType st = SolverType::QR)
: m_Y(y)
, m_blocs(), m_X()
, m_rank(), m_rss(), m_coefficients()
, m_solver_type(st)
, m_computed(false)
/*{ MSG_DEBUG("new model " << __LINE__ << " with Y(" << y.innerSize() << ',' << y.outerSize() << ')'); }*/
{}
model(const model& mo)
: m_Y(mo.m_Y)
, m_blocs(mo.m_blocs), m_X(mo.m_X)
, m_rank(mo.m_rank), m_rss(mo.m_rss), m_coefficients(mo.m_coefficients)
, m_solver_type(mo.m_solver_type)
, m_computed(mo.m_computed)
/*{ MSG_DEBUG("new model " << __LINE__ << " with Y(" << m_Y->innerSize() << ',' << m_Y->outerSize() << ')'); }*/
{}
model& operator = (const model& mo)
{
m_Y = mo.m_Y;
m_blocs = mo.m_blocs;
m_computed = mo.m_computed;
m_X = mo.m_X;
m_rss = mo.m_rss;
m_coefficients = mo.m_coefficients;
return *this;
}
bool operator == (const model& m) const
{
return m_Y == m.m_Y && m_blocs == m.m_blocs && m_rss == m.m_rss;
}
void add_bloc(const value<model_block_type>& x)
{
m_computed = false;
m_blocs.push_back(x);
}
void remove_bloc(const value<model_block_type>& x)
{
m_computed = false;
m_blocs.erase(std::find(m_blocs.begin(), m_blocs.end(), x));
}
void use_SVD()
{
m_computed = false;
m_solver_type = SolverType::SVD;
}
void use_QR()
{
m_computed = false;
m_solver_type = SolverType::QR;
}
SolverType solver_type() const
{
return m_solver_type;
}
void compute()
{
m_computed = true;
m_X = new immediate_value<MatrixXd>(MatrixXd());
m_coefficients = new immediate_value<MatrixXd>(MatrixXd());
m_rss = new immediate_value<VectorXd>(VectorXd());
*m_X = concat_right(m_blocs);
if (m_solver_type == SolverType::QR) {
int m_columns = m_X->outerSize();
/*MSG_DEBUG("X(" << X().innerSize() << ',' << X().outerSize() << ')');*/
/*MSG_DEBUG("Y(" << Y().innerSize() << ',' << Y().outerSize() << ')');*/
FullPivHouseholderQR<MatrixXd> solver(*m_X);
solver.setThreshold(COMPONENT_EPSILON);
m_rank = solver.rank();
m_coefficients->resize(m_columns, Y().outerSize());
for (int i = 0; i < Y().outerSize(); ++i) {
m_coefficients->col(i) = solver.solve(Y().col(i));
}
} else {
JacobiSVD<MatrixXd> solver(*m_X);
m_rank = 0;
int nzsv = solver.nonzeroSingularValues();
for (int i = 0; i < nzsv; ++i) {
m_rank += !around_zero(solver.singularValues()(i));
}
*m_coefficients = solver.solve(Y());
}
MatrixXd tmp = Y() - X() * (*m_coefficients);
*m_rss = tmp.array().square().colwise().sum();
}
const MatrixXd& X() const
{
if (!m_computed) { throw std::runtime_error("Model not computed"); }
return *m_X;
}
const VectorXd& rss()
{
if (!m_computed) { throw std::runtime_error("Model not computed"); }
return *m_rss;
}
const MatrixXd& coefficients()
{
if (!m_computed) { throw std::runtime_error("Model not computed"); }
return *m_coefficients;
}
int rank()
{
if (!m_computed) { throw std::runtime_error("Model not computed"); }
return m_rank;
}
const MatrixXd& Y() const
{
if (!m_Y) {
MSG_ERROR("NULL Y!", "Fix the code");
throw 0;
}
return *m_Y;
}
model extend(const value<model_block_type>& m)
{
model ret(*this);
/*MSG_DEBUG("extend model " << __LINE__ << " with Y(" << m_Y->innerSize() << ',' << m_Y->outerSize() << ')');*/
/*model ret(m_Y, m_solver_type);*/
/*ret.add_bloc(X());*/
ret.add_bloc(m);
return ret;
}
MatrixXd XtX_pseudo_inverse()
{
JacobiSVD<MatrixXd> inverter(X().transpose() * X(), ComputeFullV);
auto& V = inverter.matrixV();
VectorXd inv_sv(inverter.singularValues());
for (int i = 0; i < inv_sv.innerSize(); ++i) {
if (!around_zero(inv_sv(i))) {
inv_sv(i) = 1. / inv_sv(i);
} else {
inv_sv(i) = 0.;
}
}
void operator () (_Scalar s, int i, int j)
return V * inv_sv.asDiagonal() * V.transpose();
}
private:
value<MatrixXd> m_Y;
collection<model_block_type> m_blocs;
value<MatrixXd> m_X;
int m_rank;
value<VectorXd> m_rss;
value<MatrixXd> m_coefficients;
SolverType m_solver_type;
/*decomposition_base* m_solver;*/
bool m_computed;
public:
friend
inline
md5_digest& operator << (md5_digest& md5, const model& m)
{
md5 << m.Y() << m.m_blocs;
return md5;
}
friend
inline
std::ostream& operator << (std::ostream& os, const model& m)
{
accum ^= hash<_Scalar>()(s);
(void)i; (void)j;
os << "<model Y(" << m.m_Y->innerSize() << ',' << m.m_Y->outerSize() << "), " << m.m_blocs.size() << " blocs>";
return os;
}
};
size_t operator () (const Eigen::Matrix<_Scalar, A, B, C, D, E>& m)
size_t hash() const
{
red_mat rm;
m.visit(rm);
return rm.accum;
std::hash<model_block_type> hm;
size_t accum = 0;
for (const auto& b: m_blocs) {
accum ^= hm(b);
}
return accum;
}
};
}
#if 0
struct model {
model()
: m_Y(), m_blocs(), m_X(), m_stamps(), m_solver_type(), m_solver()
{}
model(const MatrixXd* y, SolverType st = SolverType::QR)
: m_Y(*y)
, m_blocs(), m_X()
, m_stamps()
, m_solver_type(st)
, m_solver(0)
/*{ MSG_DEBUG("new model " << __LINE__ << " with Y(" << y->innerSize() << ',' << y->outerSize() << ')'); }*/
{}
model(const MatrixXd& y, SolverType st = SolverType::QR)
model(const value<MatrixXd>& y, SolverType st = SolverType::QR)
: m_Y(y)
, m_blocs(), m_X()
, m_stamps()
......@@ -214,7 +448,7 @@ struct model {
return m_Y == m.m_Y && m_blocs == m.m_blocs && m_solver_type == m.m_solver_type;
}
void add_bloc(const MatrixXd& x)
void add_bloc(const value<model_block_type>& x)
{
m_stamps.update_blocs();
m_blocs.push_back(x);
......@@ -222,7 +456,7 @@ struct model {
/*MSG_DEBUG("stamp_blocs=" << m_stamps.m_blocs << " stamp_X=" << m_stamps.m_X);*/
}
void remove_bloc(const MatrixXd& x)
void remove_bloc(const value<model_block_type>& x)
{
m_stamps.update_blocs();
m_blocs.erase(std::find(m_blocs.begin(), m_blocs.end(), x));
......@@ -233,7 +467,7 @@ struct model {
m_solver_type = SolverType::SVD;
m_stamps.m_solver = -1;
m_stamps.m_coefficients = -1;
m_stamps.m_residuals = -1;
m_stamps.m_rss = -1;
m_stamps.m_rank = -1;
}
......@@ -242,7 +476,7 @@ struct model {
m_solver_type = SolverType::QR;
m_stamps.m_solver = -1;
m_stamps.m_coefficients = -1;
m_stamps.m_residuals = -1;
m_stamps.m_rss = -1;
m_stamps.m_rank = -1;
}
......@@ -260,21 +494,36 @@ struct model {
return m_X;
}
const MatrixXd& residuals()
const ArrayXd residuals()
{
return ArrayXd();
#if 0
if (!m_stamps.residuals_is_uptodate()) {
//*
m_residuals = Y() - X() * coefficients();
m_rss = Y() - X() * coefficients();
/*/
m_residuals = Y() - coefficients() * X();
m_rss = Y() - coefficients() * X();
//*/
m_stamps.update_residuals();
}
return m_residuals;
return m_rss;
#endif
}
const VectorXd rss()
{
scoped_lock lg(m_rss_lock);
if (!m_stamps.rss_is_uptodate()) {
MatrixXd tmp = Y() - X() * solver()->solve(Y());
m_rss = tmp.array().square().colwise().sum();
m_stamps.update_rss();
}
return m_rss;
}
const MatrixXd& coefficients()
{
scoped_lock lg(m_rank_lock);
if (!m_stamps.coefficients_is_uptodate()) {
m_coefficients = solver()->solve(Y());
/*m_coefficients = solver()->solve(Y()).transpose();*/
......@@ -285,6 +534,7 @@ struct model {
int rank()
{
scoped_lock lg(m_rank_lock);
if (!m_stamps.rank_is_uptodate()) {
m_rank = solver()->rank();
m_stamps.update_rank();
......@@ -294,19 +544,24 @@ struct model {
const MatrixXd& Y() const
{
return m_Y;
if (!m_Y) {
MSG_ERROR("NULL Y!", "Fix the code");
throw 0;
}
return *m_Y;
}
model extend(const MatrixXd& m)
model extend(const value<model_block_type>& m)
{
/*model ret(*this);*/
model ret(*this);
/*MSG_DEBUG("extend model " << __LINE__ << " with Y(" << m_Y->innerSize() << ',' << m_Y->outerSize() << ')');*/
model ret(m_Y, m_solver_type);
ret.add_bloc(X());
/*model ret(m_Y, m_solver_type);*/
/*ret.add_bloc(X());*/
ret.add_bloc(m);
return ret;
}
#if 0
model extend(const std::vector<labelled_matrix<MatrixXd, int, char>>& mv)
{
/*MSG_DEBUG("extend model " << __LINE__ << " with Y(" << m_Y->innerSize() << ',' << m_Y->outerSize() << ')');*/
......@@ -318,6 +573,7 @@ struct model {
}
return ret;
}
#endif
MatrixXd XtX_pseudo_inverse()
{
......@@ -335,91 +591,26 @@ struct model {
}
private:
#if 0
struct stamps_type {
int m_X;
int m_blocs;
int m_rank;
int m_residuals;
int m_coefficients;
int m_solver;
#if 1
struct stamp_constraint {
typedef int stamps_type::* stamp_field;
stamp_field stamp;
std::vector<stamp_field> dependencies;
stamp_constraint(stamp_field s, std::initializer_list<stamp_field> deps)
: stamp(s), dependencies(deps)
{}
bool is_uptodate(const stamps_type* stamp_this) const
{
bool ok = true;
for (auto d: dependencies) { ok &= stamp_this->*stamp >= stamp_this->*d; }
return ok;
}
void update(stamps_type* stamp_this)
{
for (auto d: dependencies) {
stamp_this->*stamp = stamp_this->*stamp >= stamp_this->*d
? stamp_this->*stamp
: stamp_this->*d;
}
}
};
stamp_constraint cons_X, cons_rank, cons_resid, cons_coef, cons_solv;
#endif
stamps_type()
: m_X(-1), m_blocs(0), m_rank(-2), m_residuals(-2), m_coefficients(-2), m_solver(-1)
, cons_X(&stamps_type::m_X, {&stamps_type::m_blocs})
, cons_rank(&stamps_type::m_rank, {&stamps_type::m_blocs, &stamps_type::m_X, &stamps_type::m_solver})
, cons_resid(&stamps_type::m_residuals, {&stamps_type::m_blocs, &stamps_type::m_X, &stamps_type::m_solver})
, cons_coef(&stamps_type::m_coefficients, {&stamps_type::m_blocs, &stamps_type::m_X, &stamps_type::m_solver})
, cons_solv(&stamps_type::m_solver, {&stamps_type::m_blocs, &stamps_type::m_X})
{}
bool X_is_uptodate() const { return cons_X.is_uptodate(this); }
bool rank_is_uptodate() const { return cons_rank.is_uptodate(this); }
bool residuals_is_uptodate() const { return cons_resid.is_uptodate(this); }
bool coefficients_is_uptodate() const { return cons_coef.is_uptodate(this); }
bool solver_is_uptodate() const { return cons_solv.is_uptodate(this); }
/*bool X_is_uptodate() const { return m_X == m_blocs; }*/
/*bool rank_is_uptodate() const { return m_rank == m_blocs; }*/
/*bool residuals_is_uptodate() const { return m_residuals == m_blocs; }*/
/*bool coefficients_is_uptodate() const { return m_coefficients == m_blocs; }*/
/*bool solver_is_uptodate() const { return m_solver == m_blocs; }*/
void update_blocs() { ++m_blocs; }
void update_X() { cons_X.update(this); }
void update_rank() { cons_rank.update(this); }
void update_solver() { cons_solv.update(this); }
void update_coefficients() { cons_coef.update(this); }
void update_residuals() { cons_resid.update(this); }
};
#endif
struct stamps_type {
bool m_X;
bool m_rank;
bool m_residuals;
bool m_rss;
bool m_coefficients;
bool m_solver;
stamps_type() : m_X(false), m_rank(false), m_residuals(false), m_coefficients(false), m_solver(false) {}
stamps_type() : m_X(false), m_rank(false), m_rss(false), m_coefficients(false), m_solver(false) {}
void update_blocs() { m_X = false; m_rank = false; m_residuals = false; m_coefficients = false; m_solver = false; }
void update_X() { m_X = true; m_rank = false; m_residuals = false; m_coefficients = false; m_solver = false; }
void update_blocs() { m_X = false; m_rank = false; m_rss = false; m_coefficients = false; m_solver = false; }
void update_X() { m_X = true; m_rank = false; m_rss = false; m_coefficients = false; m_solver = false; }
void update_rank() { m_rank = true; }
void update_solver() { m_solver = true; m_rank = false; m_residuals = false; m_coefficients = false; }
void update_solver() { m_solver = true; m_rank = false; m_rss = false; m_coefficients = false; }
void update_coefficients() { m_coefficients = true; }
void update_residuals() { m_residuals = true; }
void update_rss() { m_rss = true; }
bool X_is_uptodate() const { return m_X; }
bool rank_is_uptodate() const { return m_rank; }
bool coefficients_is_uptodate() const { return m_coefficients; }
bool residuals_is_uptodate() const { return m_residuals; }
bool rss_is_uptodate() const { return m_rss; }
bool solver_is_uptodate() const { return m_solver; }
};
......@@ -497,16 +688,23 @@ private:
SolverType type() const { return SolverType::SVD; }
};
MatrixXd m_Y;
std::vector<MatrixXd> m_blocs;
value<MatrixXd> m_Y;
collection<model_block_type> m_blocs;
MatrixXd m_X;
int m_rank;
MatrixXd m_residuals;
VectorXd m_rss;
MatrixXd m_coefficients;
stamps_type m_stamps;
SolverType m_solver_type;
decomposition_base* m_solver;
spin_lock m_blocs_lock;
spin_lock m_X_lock;
spin_lock m_rank_lock;
spin_lock m_coefficients_lock;
spin_lock m_rss_lock;
spin_lock m_solver_lock;
public:
friend
inline
......@@ -520,13 +718,13 @@ public:
inline
std::ostream& operator << (std::ostream& os, const model& m)
{
os << "<model Y(" << m.m_Y.innerSize() << ',' << m.m_Y.outerSize() << "), " << m.m_blocs.size() << " blocs>";
os << "<model Y(" << m.m_Y->innerSize() << ',' << m.m_Y->outerSize() << "), " << m.m_blocs.size() << " blocs>";
return os;
}
size_t hash() const
{
std::hash<MatrixXd> hm;
std::hash<model_block_type> hm;
size_t accum = 0;
for (const auto& b: m_blocs) {
accum ^= hm(b);
......@@ -534,6 +732,7 @@ public:
return accum;
}
};
#endif
namespace std {
template <>
......
include/model/tests.h
View file @ 4c1779a4
......@@ -4,10 +4,12 @@
static inline
VectorXd f_test(model& model_current, model& model_new)
{
auto res_cur = model_current.residuals().array();
auto res_new = model_new.residuals().array();
VectorXd rss_cur = (res_cur * res_cur).colwise().sum();
VectorXd rss_new = (res_new * res_new).colwise().sum();
/*auto res_cur = model_current.residuals().array();*/
/*auto res_new = model_new.residuals().array();*/
/*VectorXd rss_cur = (res_cur * res_cur).colwise().sum();*/
/*VectorXd rss_new = (res_new * res_new).colwise().sum();*/
VectorXd rss_cur = model_current.rss();
VectorXd rss_new = model_new.rss();
int dof_cur = model_current.rank();
int dof_new = model_new.rank();