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

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

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

11
12
13

#define COMPONENT_EPSILON (1.e-10)

14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
static inline
constexpr bool
around_zero(double o)
{
    return o < COMPONENT_EPSILON && o > -COMPONENT_EPSILON;
}

static inline
constexpr bool
much_smaller_than(double a, double b)
{
    return a < (COMPONENT_EPSILON * b);
}

static inline
void
set_if_much_smaller_than(double& a, double b)
{
    double tmp = COMPONENT_EPSILON * b;
    if (a < tmp) {
        a = tmp;
    }
}

38
39
40

using namespace Eigen;

41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
static inline
MatrixXd concat_right(const std::vector<MatrixXd>& mat_vec)
{
    size_t full_size = 0;
    MatrixXd ret;
    for (auto& m: mat_vec) {
        full_size += m.outerSize();
        /*MSG_DEBUG("preparing concat_right with matrix(" << m->innerSize() << ',' << m->outerSize() << ')');*/
    }
    ret.resize(mat_vec.front().innerSize(), full_size);
    full_size = 0;
    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;
        full_size += m.outerSize();
    }
    return ret;
}

61
static inline
62
MatrixXd concat_right(const std::vector<const MatrixXd*>& mat_vec)
63
64
65
{
    size_t full_size = 0;
    MatrixXd ret;
66
67
    for (auto m: mat_vec) {
        full_size += m->outerSize();
68
        /*MSG_DEBUG("preparing concat_right with matrix(" << m->innerSize() << ',' << m->outerSize() << ')');*/
69
    }
70
    ret.resize(mat_vec.front()->innerSize(), full_size);
71
    full_size = 0;
72
    for (auto m: mat_vec) {
73
74
75
        /*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;
76
        full_size += m->outerSize();
77
78
79
80
    }
    return ret;
}

Damien Leroux's avatar
Damien Leroux committed
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
static inline
MatrixXd concat_down(const std::vector<MatrixXd>& mat_vec)
{
    size_t full_size = 0;
    MatrixXd ret;
    for (auto& m: mat_vec) {
        full_size += m.innerSize();
    }
    ret.resize(full_size, mat_vec.front().outerSize());
    full_size = 0;
    for (auto& m: mat_vec) {
        ret.middleRows(full_size, m.innerSize()) = m;
        full_size += m.innerSize();
    }
    return ret;
}

98
static inline
99
MatrixXd concat_down(const std::vector<const MatrixXd*>& mat_vec)
100
101
102
{
    size_t full_size = 0;
    MatrixXd ret;
103
104
    for (auto m: mat_vec) {
        full_size += m->innerSize();
105
    }
106
    ret.resize(full_size, mat_vec.front()->outerSize());
107
    full_size = 0;
108
    for (auto m: mat_vec) {
109
110
        /*ret.block(full_size, 0, m->innerSize(), ret.outerSize()) = *m;*/
        ret.middleRows(full_size, m->innerSize()) = *m;
111
        full_size += m->innerSize();
112
113
114
115
116
117
118
119
120
    }
    return ret;
}


static inline
std::pair<int, MatrixXd>
rank_and_components(const MatrixXd& M)
{
121
    JacobiSVD<MatrixXd> svd(M, ComputeThinU);
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141

    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;
}


142
enum class SolverType { QR, SVD };
143

Damien Leroux's avatar
Damien Leroux committed
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

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;
        }
        void operator () (_Scalar s, int i, int j)
        {
            accum ^= hash<_Scalar>()(s);
            (void)i; (void)j;
        }
    };
    size_t operator () (const Eigen::Matrix<_Scalar, A, B, C, D, E>& m)
    {
        red_mat rm;
        m.visit(rm);
        return rm.accum;
    }
};
}

170
struct model {
171
172
173
    model()
        : m_Y(), m_blocs(), m_X(), m_stamps(), m_solver_type(), m_solver()
    {}
174
    model(const MatrixXd* y, SolverType st = SolverType::QR)
175
        : m_Y(*y)
176
177
        , m_blocs(), m_X()
        , m_stamps()
178
179
        , m_solver_type(st)
        , m_solver(0)
180
    /*{ MSG_DEBUG("new model " << __LINE__ << " with Y(" << y->innerSize() << ',' << y->outerSize() << ')'); }*/
181
182
183
    {}

    model(const MatrixXd& y, SolverType st = SolverType::QR)
184
        : m_Y(y)
185
186
187
        , m_blocs(), m_X()
        , m_stamps()
        , m_solver_type(st)
188
        , m_solver(0)
189
    /*{ MSG_DEBUG("new model " << __LINE__ << " with Y(" << y.innerSize() << ',' << y.outerSize() << ')'); }*/
190
191
    {}

192
193
194
195
196
197
    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)
198
    /*{ MSG_DEBUG("new model " << __LINE__ << " with Y(" << m_Y->innerSize() << ',' << m_Y->outerSize() << ')'); }*/
199
200
    {}

201
202
203
204
205
206
207
208
209
210
    model& operator = (const model& mo)
    {
        m_Y = mo.m_Y;
        m_blocs = mo.m_blocs;
        m_stamps = mo.m_stamps;
        m_solver_type = mo.m_solver_type;
        m_solver = NULL;
        m_stamps.m_solver = -1;
        return *this;
    }
211

Damien Leroux's avatar
Damien Leroux committed
212
213
214
215
216
    bool operator == (const model& m) const
    {
        return m_Y == m.m_Y && m_blocs == m.m_blocs && m_solver_type == m.m_solver_type;
    }

217
    void add_bloc(const MatrixXd& x)
218
219
220
    {
        m_stamps.update_blocs();
        m_blocs.push_back(x);
221
222
        /*MSG_DEBUG("added bloc to model (" << x->innerSize() << ',' << x->outerSize() << ')');*/
        /*MSG_DEBUG("stamp_blocs=" << m_stamps.m_blocs << " stamp_X=" << m_stamps.m_X);*/
223
224
    }

225
    void remove_bloc(const MatrixXd& x)
226
227
228
229
230
    {
        m_stamps.update_blocs();
        m_blocs.erase(std::find(m_blocs.begin(), m_blocs.end(), x));
    }

231
232
    void use_SVD()
    {
233
234
235
236
237
        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;
238
239
240
241
    }

    void use_QR()
    {
242
243
244
245
246
        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;
247
248
249
250
    }

    SolverType solver_type() const
    {
251
        return m_solver_type;
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
    }

    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()) {
            //*
267
            m_residuals = Y() - X() * coefficients();
268
            /*/
269
            m_residuals = Y() - coefficients() * X();
270
271
272
273
274
275
276
277
278
            //*/
            m_stamps.update_residuals();
        }
        return m_residuals;
    }

    const MatrixXd& coefficients()
    {
        if (!m_stamps.coefficients_is_uptodate()) {
279
280
            m_coefficients = solver()->solve(Y());
            /*m_coefficients = solver()->solve(Y()).transpose();*/
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
            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
    {
297
        return m_Y;
298
299
    }

300
    model extend(const MatrixXd& m)
301
302
    {
        /*model ret(*this);*/
303
        /*MSG_DEBUG("extend model " << __LINE__ << " with Y(" << m_Y->innerSize() << ',' << m_Y->outerSize() << ')');*/
304
        model ret(m_Y, m_solver_type);
305
        ret.add_bloc(X());
306
307
308
309
        ret.add_bloc(m);
        return ret;
    }

310
    model extend(const std::vector<labelled_matrix<MatrixXd, int, char>>& mv)
311
    {
312
313
314
315
316
317
        /*MSG_DEBUG("extend model " << __LINE__ << " with Y(" << m_Y->innerSize() << ',' << m_Y->outerSize() << ')');*/
        /*model ret(*this);*/
        model ret(m_Y, m_solver_type);
        ret.add_bloc(X());
        for (auto& m: mv) {
            ret.add_bloc(m.data);
318
        }
319
320
321
322
323
324
325
326
327
328
329
        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);
330
            } else {
331
                inv_sv(i) = 0.;
332
333
            }
        }
334
        return V * inv_sv.asDiagonal() * V.transpose();
335
336
    }

337
private:
338
#if 0
339
    struct stamps_type {
340
341
342
343
344
345
346
        int m_X;
        int m_blocs;
        int m_rank;
        int m_residuals;
        int m_coefficients;
        int m_solver;

347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
#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; }*/
394
        void update_blocs() { ++m_blocs; }
395
396
397
398
399
        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); }
400
    };
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
#endif

    struct stamps_type {
        bool m_X;
        bool m_rank;
        bool m_residuals;
        bool m_coefficients;
        bool m_solver;

        stamps_type() : m_X(false), m_rank(false), m_residuals(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_rank() { m_rank = true; }
        void update_solver() { m_solver = true; m_rank = false; m_residuals = false; m_coefficients = false; }
        void update_coefficients() { m_coefficients = true; }
        void update_residuals() { m_residuals = 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 solver_is_uptodate() const { return m_solver; }
    };
425
426
427
428
429

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

433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
    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;
    }

450
451
452
    decomposition_base* solver()
    {
        if (!(m_solver != NULL && m_stamps.solver_is_uptodate())) {
453
            return _new_solver();
454
455
456
457
458
        }
        return m_solver;
    }

    struct decomposition_QR : decomposition_base {
459
        /*
460
        ColPivHouseholderQR<MatrixXd> m_solver;
461
462
463
464
        /*/
        // Doesn't seem to perform well when there are many more rows than columns
        FullPivHouseholderQR<MatrixXd> m_solver;
        //*/
465
        int m_columns;
466
        decomposition_QR(const MatrixXd& m) : m_solver(m), m_columns(m.outerSize()) { m_solver.setThreshold(COMPONENT_EPSILON); }
467
        int rank() const { return m_solver.rank(); }
468
        MatrixXd solve(const MatrixXd& lhs)
469
        {
470
471
472
473
            /*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;*/
474
                /*std::cout << "lhs = " << lhs.col(i).transpose() << std::endl;*/
475
                ret.col(i) = m_solver.solve(lhs.col(i));
476
            }
477
            /*std::cout << "QR done" << std::endl;*/
478
479
480
481
482
483
484
485
            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) {}
486
487
488
489
490
491
492
493
494
495
        int rank() const
        {
            int nz_accum = 0;
            int nzsv = m_solver.nonzeroSingularValues();

            for (int i = 0; i < nzsv; ++i) {
                nz_accum += !around_zero(m_solver.singularValues()(i));
            }
            return nz_accum;
        }
496
        MatrixXd solve(const MatrixXd& lhs) { return m_solver.solve(lhs); }
497
498
499
        SolverType type() const { return SolverType::SVD; }
    };

500
501
    MatrixXd m_Y;
    std::vector<MatrixXd> m_blocs;
502
503
504
505
506
    MatrixXd m_X;
    int m_rank;
    MatrixXd m_residuals;
    MatrixXd m_coefficients;
    stamps_type m_stamps;
507
    SolverType m_solver_type;
508
509
    decomposition_base* m_solver;

Damien Leroux's avatar
Damien Leroux committed
510
511
512
513
514
515
516
public:
    friend
    inline
    md5_digest& operator << (md5_digest& md5, const model& m)
    {
        md5 << m.Y() << m.m_blocs;
        return md5;
517
518
    }

Damien Leroux's avatar
Damien Leroux committed
519
520
521
522
523
524
    friend
        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>";
            return os;
525
526
        }

Damien Leroux's avatar
Damien Leroux committed
527
528
529
530
531
532
    size_t hash() const
    {
        std::hash<MatrixXd> hm;
        size_t accum = 0;
        for (const auto& b: m_blocs) {
            accum ^= hm(b);
533
        }
Damien Leroux's avatar
Damien Leroux committed
534
        return accum;
535
    }
Damien Leroux's avatar
Damien Leroux committed
536
537
538
539
540
541
542
};

namespace std {
template <>
struct hash<model> {
    size_t operator () (const model& m) { return m.hash(); }
};
543
544
545

}

546
547
#endif