model.h 5.83 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
#ifndef _SPEL_MODEL_MODEL_H_
#define _SPEL_MODEL_MODEL_H_

#include <Eigen/Core>
#include <Eigen/SVD>
#include <Eigen/QR>


#define COMPONENT_EPSILON (1.e-10)


using namespace Eigen;

static inline
MatrixXd concat_right(const std::vector<MatrixXd>& mat_vec)
{
    size_t full_size = 0;
    MatrixXd ret;
    for (const auto& m: mat_vec) {
        full_size += m.outerSize();
    }
    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();
    }
    return ret;
}

static inline
MatrixXd concat_down(const std::vector<MatrixXd>& mat_vec)
{
    size_t full_size = 0;
    MatrixXd ret;
    for (const auto& m: mat_vec) {
        full_size += m.innerSize();
    }
    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();
    }
    return ret;
}


static inline
std::pair<int, MatrixXd>
rank_and_components(const MatrixXd& M)
{
    JacobiSVD<MatrixXd> svd(M, ComputeThinU, COMPONENT_EPSILON);

    std::cout << "Singular values " << svd.singularValues().transpose() << std::endl;
    int nzsv = svd.nonzeroSingularValues();

    return {nzsv, svd.matrixU().leftCols(nzsv)};
}


static inline
MatrixXd components(const MatrixXd& M, const MatrixXd& P)
{
    MatrixXd pnorm(P.innerSize(), P.outerSize());
    for (int i = 0; i < P.outerSize(); ++i) {
        pnorm.col(i) = P.col(i).normalized();
    }
    MatrixXd orth = M - pnorm * pnorm.transpose() * M; /* feu ! */
    return rank_and_components(orth).second;
}


enum class SolverType { Null, QR, SVD };

struct model {
    model(const MatrixXd& y)
        : m_Y(y)
        , m_blocs(), m_X()
        , m_stamps()
        , m_solver(0)
    {}

    void add_bloc(const MatrixXd& x)
    {
        m_stamps.update_blocs();
        m_blocs.push_back(x);
    }

    void use_SVD()
    {
        if (m_solver) { delete m_solver; }
        m_solver = new decomposition_SVD(X());
        m_stamps.update_solver();
    }

    void use_QR()
    {
        if (m_solver) { delete m_solver; }
        m_solver = new decomposition_QR(X());
        m_stamps.update_solver();
    }

    SolverType solver_type() const
    {
        if (!m_solver) {
            return SolverType::Null;
        }
        return m_solver->type();
    }

    const MatrixXd& X()
    {
        if (!m_stamps.X_is_uptodate()) {
            m_X = concat_right(m_blocs);
            m_stamps.update_X();
        }
        return m_X;
    }

    const MatrixXd& residuals()
    {
        if (!m_stamps.residuals_is_uptodate()) {
            //*
            m_residuals = m_Y - X() * coefficients();
            /*/
            m_residuals = m_Y - coefficients() * X();
            //*/
            m_stamps.update_residuals();
        }
        return m_residuals;
    }

    const MatrixXd& coefficients()
    {
        if (!m_stamps.coefficients_is_uptodate()) {
            m_coefficients = solver()->solve(m_Y);
            /*m_coefficients = solver()->solve(m_Y).transpose();*/
            m_stamps.update_coefficients();
        }
        return m_coefficients;
    }

    int rank()
    {
        if (!m_stamps.rank_is_uptodate()) {
            m_rank = solver()->rank();
            m_stamps.update_rank();
        }
        return m_rank;
    }

    const MatrixXd& Y() const
    {
        return m_Y;
    }

private:
    struct stamps_type {
        int m_X;
        int m_blocs;
        int m_rank;
        int m_residuals;
        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; }
        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; }
    };

    struct decomposition_base {
        virtual MatrixXd solve(const MatrixXd& column) = 0;
        virtual int rank() const = 0;
        virtual SolverType type() const = 0;
    };

    decomposition_base* solver()
    {
        if (!(m_solver != NULL && m_stamps.solver_is_uptodate())) {
            /* whine */
            throw 0;
        }
        return m_solver;
    }

    struct decomposition_QR : decomposition_base {
        ColPivHouseholderQR<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 ret(rhs.innerSize(), m_columns);
            for (int i = 0; i < m_columns; ++i) {
                ret.col(i) = m_solver.solve(rhs.col(i));
            }
            return ret;
        }
        SolverType type() const { return SolverType::QR; }
    };

    struct decomposition_SVD : decomposition_base {
        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); }
        SolverType type() const { return SolverType::SVD; }
    };

    MatrixXd m_Y;
    std::vector<MatrixXd> m_blocs;
    MatrixXd m_X;
    int m_rank;
    MatrixXd m_residuals;
    MatrixXd m_coefficients;
    stamps_type m_stamps;
    decomposition_base* m_solver;
};


#endif