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

Implemented -log10(pvalue(F-test)).

parent 6a2dabd3
Hide whitespace changes
Inline Side-by-side
include/model/model.h
View file @ dd689f2f
......@@ -5,6 +5,9 @@
#include <Eigen/SVD>
#include <Eigen/QR>
#include <cmath>
#include <boost/math/special_functions/beta.hpp>
#define COMPONENT_EPSILON (1.e-10)
......@@ -12,35 +15,35 @@
using namespace Eigen;
static inline
MatrixXd concat_right(const std::vector<MatrixXd>& mat_vec)
MatrixXd concat_right(const std::vector<const MatrixXd*>& mat_vec)
{
size_t full_size = 0;
MatrixXd ret;
for (const auto& m: mat_vec) {
full_size += m.outerSize();
for (auto m: mat_vec) {
full_size += m->outerSize();
}
ret.resize(mat_vec.front().innerSize(), full_size);
ret.resize(mat_vec.front()->innerSize(), full_size);
full_size = 0;
for (const auto& m: mat_vec) {
ret.block(0, full_size, ret.innerSize(), m.outerSize()) = m;
full_size += m.outerSize();
for (auto m: mat_vec) {
ret.block(0, full_size, ret.innerSize(), m->outerSize()) = *m;
full_size += m->outerSize();
}
return ret;
}
static inline
MatrixXd concat_down(const std::vector<MatrixXd>& mat_vec)
MatrixXd concat_down(const std::vector<const MatrixXd*>& mat_vec)
{
size_t full_size = 0;
MatrixXd ret;
for (const auto& m: mat_vec) {
full_size += m.innerSize();
for (auto m: mat_vec) {
full_size += m->innerSize();
}
ret.resize(full_size, mat_vec.front().outerSize());
ret.resize(full_size, mat_vec.front()->outerSize());
full_size = 0;
for (const auto& m: mat_vec) {
ret.block(full_size, 0, m.innerSize(), ret.outerSize()) = m;
full_size += m.innerSize();
for (auto m: mat_vec) {
ret.block(full_size, 0, m->innerSize(), ret.outerSize()) = *m;
full_size += m->innerSize();
}
return ret;
}
......@@ -71,42 +74,69 @@ MatrixXd components(const MatrixXd& M, const MatrixXd& P)
}
enum class SolverType { Null, QR, SVD };
enum class SolverType { QR, SVD };
struct model {
model(const MatrixXd& y)
model(const MatrixXd* y, SolverType st = SolverType::QR)
: m_Y(y)
, m_blocs(), m_X()
, m_stamps()
, m_solver_type(st)
, m_solver(0)
{}
model(const MatrixXd& y, SolverType st = SolverType::QR)
: m_Y(&y)
, m_blocs(), m_X()
, m_stamps()
, m_solver_type(st)
, m_solver(0)
{}
void add_bloc(const MatrixXd& x)
model(const model& mo)
: m_Y(mo.m_Y)
, m_blocs(mo.m_blocs), m_X()
, m_stamps()
, m_solver_type(mo.m_solver_type)
, m_solver(0)
{}
void add_bloc(const MatrixXd& x) { add_bloc(&x); }
void remove_bloc(const MatrixXd& x) { remove_bloc(&x); }
void add_bloc(const MatrixXd* x)
{
m_stamps.update_blocs();
m_blocs.push_back(x);
}
void remove_bloc(const MatrixXd* x)
{
m_stamps.update_blocs();
m_blocs.erase(std::find(m_blocs.begin(), m_blocs.end(), x));
}
void use_SVD()
{
if (m_solver) { delete m_solver; }
m_solver = new decomposition_SVD(X());
m_stamps.update_solver();
m_solver_type = SolverType::SVD;
m_stamps.m_solver = -1;
m_stamps.m_coefficients = -1;
m_stamps.m_residuals = -1;
m_stamps.m_rank = -1;
}
void use_QR()
{
if (m_solver) { delete m_solver; }
m_solver = new decomposition_QR(X());
m_stamps.update_solver();
m_solver_type = SolverType::QR;
m_stamps.m_solver = -1;
m_stamps.m_coefficients = -1;
m_stamps.m_residuals = -1;
m_stamps.m_rank = -1;
}
SolverType solver_type() const
{
if (!m_solver) {
return SolverType::Null;
}
return m_solver->type();
return m_solver_type;
}
const MatrixXd& X()
......@@ -122,9 +152,9 @@ struct model {
{
if (!m_stamps.residuals_is_uptodate()) {
//*
m_residuals = m_Y - X() * coefficients();
m_residuals = Y() - X() * coefficients();
/*/
m_residuals = m_Y - coefficients() * X();
m_residuals = Y() - coefficients() * X();
//*/
m_stamps.update_residuals();
}
......@@ -134,8 +164,8 @@ struct model {
const MatrixXd& coefficients()
{
if (!m_stamps.coefficients_is_uptodate()) {
m_coefficients = solver()->solve(m_Y);
/*m_coefficients = solver()->solve(m_Y).transpose();*/
m_coefficients = solver()->solve(Y());
/*m_coefficients = solver()->solve(Y()).transpose();*/
m_stamps.update_coefficients();
}
return m_coefficients;
......@@ -152,11 +182,98 @@ struct model {
const MatrixXd& Y() const
{
return m_Y;
return *m_Y;
}
model extend(const MatrixXd& m) { return extend(&m); }
model extend(const MatrixXd* m)
{
/*model ret(*this);*/
model ret(m_Y, m_solver_type);
ret.add_bloc(&X());
ret.add_bloc(m);
return ret;
}
MatrixXd f_test(const MatrixXd& m) { return f_test(&m); }
#if 0
F=((RSS_current - RSS_new)/(dof_new-dof_current))
F/= rss_new /( nind - dof_new)
#endif
VectorXd f_test(const MatrixXd* m)
{
model ext = extend(m);
auto res_cur = residuals().array();
auto res_new = ext.residuals().array();
VectorXd rss_cur = (res_cur * res_cur).colwise().sum();
VectorXd rss_new = (res_new * res_new).colwise().sum();
int dof_cur = rank();
int dof_new = ext.rank();
int nind = m_Y->innerSize();
if (dof_new <= dof_cur) {
std::cout << "dof_new[" << dof_new << "] <= dof_cur[" << dof_cur << ']' << std::endl;
return VectorXd::Zero(m_Y->innerSize());
}
if (nind <= dof_new) {
std::cout << "nind[" << nind << "] <= dof_new[" << dof_new << ']' << std::endl;
/* FIXME: HALT! */
return VectorXd::Zero(m_Y->innerSize());
}
VectorXd F(rss_new.innerSize());
double dof_num = dof_new - dof_cur;
double dof_denom = nind - dof_new;
double dof_num_inv = 1. / dof_num;
double dof_denom_inv = 1. / dof_denom;
for (int i = 0; i < rss_new.innerSize(); ++i) {
if (rss_new(i) == 0) {
F(i) = 100.; /* FIXME: magic number */
} else if (dof_num > 0) {
F(i) = ((rss_cur(i) - rss_new(i)) * dof_num_inv) / (rss_new(i) * dof_denom_inv);
if (F(i) < 0) {
F(i) = 0;
} else {
/* NdDL : c'est immonde. */
F(i) = -log10(boost::math::ibeta(dof_denom * .5, dof_num * .5, dof_denom / (dof_denom + dof_num * F(i))));
}
} else {
F(i) = 0;
}
}
return F;
}
private:
struct stamps_type {
#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
int m_X;
int m_blocs;
int m_rank;
......@@ -164,17 +281,32 @@ private:
int m_coefficients;
int m_solver;
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; }
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() { m_X = m_blocs; }
void update_rank() { m_rank = m_blocs; }
void update_solver() { m_solver = m_blocs; }
void update_coefficients() { m_coefficients = m_blocs; }
void update_residuals() { m_residuals = 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); }
};
struct decomposition_base {
......@@ -183,26 +315,50 @@ private:
virtual SolverType type() const = 0;
};
decomposition_base* _new_solver()
{
if (m_solver) {
delete m_solver;
}
switch (m_solver_type) {
case SolverType::QR:
m_solver = new decomposition_QR(X());
break;
case SolverType::SVD:
m_solver = new decomposition_SVD(X());
break;
};
m_stamps.update_solver();
return m_solver;
}
decomposition_base* solver()
{
if (!(m_solver != NULL && m_stamps.solver_is_uptodate())) {
/* whine */
throw 0;
return _new_solver();
}
return m_solver;
}
struct decomposition_QR : decomposition_base {
//*
ColPivHouseholderQR<MatrixXd> m_solver;
/*/
// Doesn't seem to perform well when there are many more rows than columns
FullPivHouseholderQR<MatrixXd> m_solver;
//*/
int m_columns;
decomposition_QR(const MatrixXd& m) : m_solver(m), m_columns(m.outerSize()) {}
int rank() const { return m_solver.rank(); }
MatrixXd solve(const MatrixXd& rhs)
MatrixXd solve(const MatrixXd& lhs)
{
MatrixXd ret(rhs.innerSize(), m_columns);
for (int i = 0; i < m_columns; ++i) {
ret.col(i) = m_solver.solve(rhs.col(i));
/*std::cout << "QR " << lhs.outerSize() << " columns" << std::endl;*/
MatrixXd ret(m_columns, lhs.outerSize());
for (int i = 0; i < lhs.outerSize(); ++i) {
/*std::cout << "QR processing column #" << i << std::endl;*/
ret.col(i) = m_solver.solve(lhs.col(i));
}
/*std::cout << "QR done" << std::endl;*/
return ret;
}
SolverType type() const { return SolverType::QR; }
......@@ -212,17 +368,18 @@ private:
JacobiSVD<MatrixXd> m_solver;
decomposition_SVD(const MatrixXd& m) : m_solver(m, ComputeThinU|ComputeThinV) {}
int rank() const { return m_solver.nonzeroSingularValues(); }
MatrixXd solve(const MatrixXd& rhs) { return m_solver.solve(rhs); }
MatrixXd solve(const MatrixXd& lhs) { return m_solver.solve(lhs); }
SolverType type() const { return SolverType::SVD; }
};
MatrixXd m_Y;
std::vector<MatrixXd> m_blocs;
const MatrixXd* m_Y;
std::vector<const MatrixXd*> m_blocs;
MatrixXd m_X;
int m_rank;
MatrixXd m_residuals;
MatrixXd m_coefficients;
stamps_type m_stamps;
SolverType m_solver_type;
decomposition_base* m_solver;
};
......
src/test_model.cc
View file @ dd689f2f
#include "model.h"
#include "chrono.h"
MatrixXd test_comp(const MatrixXd& M, const MatrixXd& P)
{
......@@ -10,10 +11,15 @@ MatrixXd test_comp(const MatrixXd& M, const MatrixXd& P)
}
bool closeTo(const MatrixXd& a, const MatrixXd& b, double prec = 1.e-10)
{
return (a - b).isMuchSmallerThan(1., prec) || a.isApprox(b, prec);
}
int main(int argc, char** argv)
{
#if 0
/*MatrixXd P = MatrixXd::Ones(5, 1);*/
/*MatrixXd M = MatrixXd::Random(5, 2);*/
MatrixXd tmp(5, 2);
......@@ -46,8 +52,8 @@ int main(int argc, char** argv)
MatrixXd gros = concat_right({bouh, P});
test_comp(M, gros);
/*model m(tmp);*/
model m(MatrixXd::Ones(5, 1));
model m(tmp);
/*model m(MatrixXd::Ones(5, 1));*/
m.add_bloc(zob);
std::cout << "MODEL Y" << std::endl << m.Y() << std::endl;
std::cout << "MODEL X" << std::endl << m.X() << std::endl;
......@@ -59,7 +65,84 @@ int main(int argc, char** argv)
std::cout << "MODEL X RANK " << m.rank() << std::endl;
std::cout << "MODEL COEFFICIENTS" << std::endl << m.coefficients() << std::endl;
std::cout << "MODEL RESIDUALS" << std::endl << m.residuals() << std::endl;
#endif
#define NIND 1000
#define NTRUC 210
#define XRANK 200
#define Hij(_x_) (1./(_x_))
MatrixXd Y = MatrixXd::Random(NIND, 20);
model bigmodel(Y);
#if 0
MatrixXd H(NIND, NIND);
for (int i = 0; i < H.innerSize(); ++i) {
for (int j = 0; j < H.outerSize(); ++j) {
H(i, j) = 1. / (1 + i + j);
}
}
bigmodel.add_bloc(H);
#else
MatrixXd X(NIND, NTRUC);
X.leftCols(XRANK) = MatrixXd::Random(NIND, XRANK);
for (int i = XRANK; i < NTRUC; i += XRANK) {
X.block(0, i, NIND, (i + XRANK) <= NTRUC ? XRANK : (NTRUC - i - XRANK)) = X.leftCols(XRANK);
}
/*MatrixXd X = MatrixXd::Random(NIND, NTRUC);*/
bigmodel.add_bloc(X);
#endif
bigmodel.X(); /* pour initialiser X avant le chrono */
std::cout << "Model rank = " << bigmodel.rank() << std::endl;
MatrixXd coefQR, coefSVD, residQR, residSVD;
chrono::start("SVD gros modèle");
bigmodel.use_SVD();
coefSVD = bigmodel.coefficients();
std::cout << "SVD coefficients (" << coefSVD.innerSize() << ',' << coefSVD.outerSize() << ')' << std::endl;
residSVD = bigmodel.residuals();
std::cout << "SVD residuals (" << residSVD.innerSize() << ',' << residSVD.outerSize() << ')' << std::endl;
chrono::stop("SVD gros modèle");
chrono::start("QR gros modèle");
bigmodel.use_QR();
coefQR = bigmodel.coefficients();
std::cout << "QR coefficients (" << coefQR.innerSize() << ',' << coefQR.outerSize() << ')' << std::endl;
residQR = bigmodel.residuals();
std::cout << "QR residuals (" << residQR.innerSize() << ',' << residQR.outerSize() << ')' << std::endl;
chrono::stop("QR gros modèle");
std::cout << "coef QR == coef SVD ? " << closeTo(coefQR, coefSVD, 1.e-15) << closeTo(coefQR, coefSVD, 1.e-10) << closeTo(coefQR, coefSVD, 1.e-5) << closeTo(coefQR, coefSVD, 1.e-1) << std::endl;
std::cout << "resid QR == resid SVD ? " << closeTo(residQR, residSVD, 1.e-15) << closeTo(residQR, residSVD, 1.e-10) << closeTo(residQR, residSVD, 1.e-5) << closeTo(residQR, residSVD, 1.e-1) << std::endl;
/*std::cout << "resid QR" << std::endl << residQR << std::endl << std::endl;*/
/*std::cout << "resid SVD" << std::endl << residSVD << std::endl << std::endl;*/
/*std::cout << "coef QR" << std::endl << coefQR << std::endl << std::endl;*/
/*std::cout << "coef SVD" << std::endl << coefSVD << std::endl << std::endl;*/
{
MatrixXd X1(6, 3);
MatrixXd X2(6, 2);
MatrixXd Y(6, 1);
X1 << 1, 2, 3,
2, 2, 2,
1, 3, 4,
2, 3, 5,
1, 2, 1,
0, 1, 2;
X2 << 5, 1,
3, 1,
2, 2,
-1, 1,
1, 0,
1, 1;
Y << 12, 10, 12, 10, 5, 5;
model M(Y);
M.add_bloc(X1);
M.use_QR();
std::cout << "M1 residuals " << M.residuals().transpose() << std::endl;
auto F = M.f_test(X2);
auto M2 = M.extend(X2);
std::cout << "F(M1, M2) " << F.transpose() << std::endl;
std::cout << "M2.residuals " << M2.residuals().transpose() << std::endl;
}
return 0;
}
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