model.h 11.5 KB
Newer Older
1
2
3
4
5
6
7
#ifndef _SPEL_MODEL_MODEL_H_
#define _SPEL_MODEL_MODEL_H_

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

8
9
10
#include <cmath>
#include <boost/math/special_functions/beta.hpp>

11
12
13
14
15
16
17

#define COMPONENT_EPSILON (1.e-10)


using namespace Eigen;

static inline
18
MatrixXd concat_right(const std::vector<const MatrixXd*>& mat_vec)
19
20
21
{
    size_t full_size = 0;
    MatrixXd ret;
22
23
    for (auto m: mat_vec) {
        full_size += m->outerSize();
24
    }
25
    ret.resize(mat_vec.front()->innerSize(), full_size);
26
    full_size = 0;
27
28
29
    for (auto m: mat_vec) {
        ret.block(0, full_size, ret.innerSize(), m->outerSize()) = *m;
        full_size += m->outerSize();
30
31
32
33
34
    }
    return ret;
}

static inline
35
MatrixXd concat_down(const std::vector<const MatrixXd*>& mat_vec)
36
37
38
{
    size_t full_size = 0;
    MatrixXd ret;
39
40
    for (auto m: mat_vec) {
        full_size += m->innerSize();
41
    }
42
    ret.resize(full_size, mat_vec.front()->outerSize());
43
    full_size = 0;
44
45
46
    for (auto m: mat_vec) {
        ret.block(full_size, 0, m->innerSize(), ret.outerSize()) = *m;
        full_size += m->innerSize();
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
    }
    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;
}


77
enum class SolverType { QR, SVD };
78
79

struct model {
80
    model(const MatrixXd* y, SolverType st = SolverType::QR)
81
82
83
        : m_Y(y)
        , m_blocs(), m_X()
        , m_stamps()
84
85
86
87
88
89
90
91
92
        , 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)
93
94
95
        , m_solver(0)
    {}

96
97
98
99
100
101
102
103
104
105
106
107
    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)
108
109
110
111
112
    {
        m_stamps.update_blocs();
        m_blocs.push_back(x);
    }

113
114
115
116
117
118
    void remove_bloc(const MatrixXd* x)
    {
        m_stamps.update_blocs();
        m_blocs.erase(std::find(m_blocs.begin(), m_blocs.end(), x));
    }

119
120
    void use_SVD()
    {
121
122
123
124
125
        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;
126
127
128
129
    }

    void use_QR()
    {
130
131
132
133
134
        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;
135
136
137
138
    }

    SolverType solver_type() const
    {
139
        return m_solver_type;
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
    }

    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()) {
            //*
155
            m_residuals = Y() - X() * coefficients();
156
            /*/
157
            m_residuals = Y() - coefficients() * X();
158
159
160
161
162
163
164
165
166
            //*/
            m_stamps.update_residuals();
        }
        return m_residuals;
    }

    const MatrixXd& coefficients()
    {
        if (!m_stamps.coefficients_is_uptodate()) {
167
168
            m_coefficients = solver()->solve(Y());
            /*m_coefficients = solver()->solve(Y()).transpose();*/
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
            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
    {
185
        return *m_Y;
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
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
    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;
    }


248
249
private:
    struct stamps_type {
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
#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

277
278
279
280
281
282
283
        int m_X;
        int m_blocs;
        int m_rank;
        int m_residuals;
        int m_coefficients;
        int m_solver;

284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
        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; }*/
304
        void update_blocs() { ++m_blocs; }
305
306
307
308
309
        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); }
310
311
312
313
314
315
316
317
    };

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

318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
    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;
    }

335
336
337
    decomposition_base* solver()
    {
        if (!(m_solver != NULL && m_stamps.solver_is_uptodate())) {
338
            return _new_solver();
339
340
341
342
343
        }
        return m_solver;
    }

    struct decomposition_QR : decomposition_base {
344
        //*
345
        ColPivHouseholderQR<MatrixXd> m_solver;
346
347
348
349
        /*/
        // Doesn't seem to perform well when there are many more rows than columns
        FullPivHouseholderQR<MatrixXd> m_solver;
        //*/
350
351
352
        int m_columns;
        decomposition_QR(const MatrixXd& m) : m_solver(m), m_columns(m.outerSize()) {}
        int rank() const { return m_solver.rank(); }
353
        MatrixXd solve(const MatrixXd& lhs)
354
        {
355
356
357
358
359
            /*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));
360
            }
361
            /*std::cout << "QR done" << std::endl;*/
362
363
364
365
366
367
368
369
370
            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(); }
371
        MatrixXd solve(const MatrixXd& lhs) { return m_solver.solve(lhs); }
372
373
374
        SolverType type() const { return SolverType::SVD; }
    };

375
376
    const MatrixXd* m_Y;
    std::vector<const MatrixXd*> m_blocs;
377
378
379
380
381
    MatrixXd m_X;
    int m_rank;
    MatrixXd m_residuals;
    MatrixXd m_coefficients;
    stamps_type m_stamps;
382
    SolverType m_solver_type;
383
384
385
386
387
388
    decomposition_base* m_solver;
};


#endif