model.h 19.7 KB
Newer Older
1
2
3
#ifndef _SPEL_MODEL_MODEL_H_
#define _SPEL_MODEL_MODEL_H_

4
5
#include <stdexcept>

6
7
8
#include <Eigen/SVD>
#include <Eigen/QR>

9
10
#include <cmath>
#include <boost/math/special_functions/beta.hpp>
11
#include <boost/math/special_functions/gamma.hpp>
12

13
14
15

#define COMPONENT_EPSILON (1.e-10)

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

82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
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;
    }
}

106
107
108

using namespace Eigen;

109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
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;
}

129
static inline
130
MatrixXd concat_right(const collection<model_block_type>& mat_vec)
131
132
133
{
    size_t full_size = 0;
    MatrixXd ret;
134
135
    for (auto m: mat_vec) {
        full_size += m->outerSize();
136
        /*MSG_DEBUG("preparing concat_right with matrix(" << m->innerSize() << ',' << m->outerSize() << ')');*/
137
    }
138
    ret.resize(mat_vec.front()->innerSize(), full_size);
139
    full_size = 0;
140
    for (auto m: mat_vec) {
141
142
        /*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;*/
143
        ret.middleCols(full_size, m->outerSize()) = m->data;
144
        full_size += m->outerSize();
145
146
147
148
    }
    return ret;
}

Damien Leroux's avatar
Damien Leroux committed
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
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;
}

166
static inline
167
MatrixXd concat_down(const std::vector<const MatrixXd*>& mat_vec)
168
169
170
{
    size_t full_size = 0;
    MatrixXd ret;
171
172
    for (auto m: mat_vec) {
        full_size += m->innerSize();
173
    }
174
    ret.resize(full_size, mat_vec.front()->outerSize());
175
    full_size = 0;
176
    for (auto m: mat_vec) {
177
178
        /*ret.block(full_size, 0, m->innerSize(), ret.outerSize()) = *m;*/
        ret.middleRows(full_size, m->innerSize()) = *m;
179
        full_size += m->innerSize();
180
181
182
183
184
185
186
187
188
    }
    return ret;
}


static inline
std::pair<int, MatrixXd>
rank_and_components(const MatrixXd& M)
{
189
    JacobiSVD<MatrixXd> svd(M, ComputeThinU);
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209

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


210
enum class SolverType { QR, SVD };
211

Damien Leroux's avatar
Damien Leroux committed
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
248
249
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
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
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.;
            }
Damien Leroux's avatar
Damien Leroux committed
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
394
        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)
Damien Leroux's avatar
Damien Leroux committed
395
        {
396
397
            os << "<model Y(" << m.m_Y->innerSize() << ',' << m.m_Y->outerSize() << "), " << m.m_blocs.size() << " blocs>";
            return os;
Damien Leroux's avatar
Damien Leroux committed
398
        }
399
400

    size_t hash() const
Damien Leroux's avatar
Damien Leroux committed
401
    {
402
403
404
405
406
407
        std::hash<model_block_type> hm;
        size_t accum = 0;
        for (const auto& b: m_blocs) {
            accum ^= hm(b);
        }
        return accum;
Damien Leroux's avatar
Damien Leroux committed
408
409
410
    }
};

411
#if 0
412
struct model {
413
414
415
    model()
        : m_Y(), m_blocs(), m_X(), m_stamps(), m_solver_type(), m_solver()
    {}
416

417
    model(const value<MatrixXd>& y, SolverType st = SolverType::QR)
418
        : m_Y(y)
419
420
421
        , m_blocs(), m_X()
        , m_stamps()
        , m_solver_type(st)
422
        , m_solver(0)
423
    /*{ MSG_DEBUG("new model " << __LINE__ << " with Y(" << y.innerSize() << ',' << y.outerSize() << ')'); }*/
424
425
    {}

426
427
428
429
430
431
    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)
432
    /*{ MSG_DEBUG("new model " << __LINE__ << " with Y(" << m_Y->innerSize() << ',' << m_Y->outerSize() << ')'); }*/
433
434
    {}

435
436
437
438
439
440
441
442
443
444
    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;
    }
445

Damien Leroux's avatar
Damien Leroux committed
446
447
448
449
450
    bool operator == (const model& m) const
    {
        return m_Y == m.m_Y && m_blocs == m.m_blocs && m_solver_type == m.m_solver_type;
    }

451
    void add_bloc(const value<model_block_type>& x)
452
453
454
    {
        m_stamps.update_blocs();
        m_blocs.push_back(x);
455
456
        /*MSG_DEBUG("added bloc to model (" << x->innerSize() << ',' << x->outerSize() << ')');*/
        /*MSG_DEBUG("stamp_blocs=" << m_stamps.m_blocs << " stamp_X=" << m_stamps.m_X);*/
457
458
    }

459
    void remove_bloc(const value<model_block_type>& x)
460
461
462
463
464
    {
        m_stamps.update_blocs();
        m_blocs.erase(std::find(m_blocs.begin(), m_blocs.end(), x));
    }

465
466
    void use_SVD()
    {
467
468
469
        m_solver_type = SolverType::SVD;
        m_stamps.m_solver = -1;
        m_stamps.m_coefficients = -1;
470
        m_stamps.m_rss = -1;
471
        m_stamps.m_rank = -1;
472
473
474
475
    }

    void use_QR()
    {
476
477
478
        m_solver_type = SolverType::QR;
        m_stamps.m_solver = -1;
        m_stamps.m_coefficients = -1;
479
        m_stamps.m_rss = -1;
480
        m_stamps.m_rank = -1;
481
482
483
484
    }

    SolverType solver_type() const
    {
485
        return m_solver_type;
486
487
488
489
490
491
492
493
494
495
496
    }

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

497
    const ArrayXd residuals()
498
    {
499
500
		return ArrayXd();
#if 0
501
502
        if (!m_stamps.residuals_is_uptodate()) {
            //*
503
            m_rss = Y() - X() * coefficients();
504
            /*/
505
            m_rss = Y() - coefficients() * X();
506
507
508
            //*/
            m_stamps.update_residuals();
        }
509
510
        return m_rss;
#endif
511
512
    }

513
514
515
516
517
518
519
520
521
522
523
	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;
	}

524
525
    const MatrixXd& coefficients()
    {
526
		scoped_lock lg(m_rank_lock);
527
        if (!m_stamps.coefficients_is_uptodate()) {
528
529
            m_coefficients = solver()->solve(Y());
            /*m_coefficients = solver()->solve(Y()).transpose();*/
530
531
532
533
534
535
536
            m_stamps.update_coefficients();
        }
        return m_coefficients;
    }

    int rank()
    {
537
		scoped_lock lg(m_rank_lock);
538
539
540
541
542
543
544
545
546
        if (!m_stamps.rank_is_uptodate()) {
            m_rank = solver()->rank();
            m_stamps.update_rank();
        }
        return m_rank;
    }

    const MatrixXd& Y() const
    {
547
548
549
550
551
		if (!m_Y) {
			MSG_ERROR("NULL Y!", "Fix the code");
			throw 0;
		}
        return *m_Y;
552
553
    }

554
    model extend(const value<model_block_type>& m)
555
    {
556
        model ret(*this);
557
        /*MSG_DEBUG("extend model " << __LINE__ << " with Y(" << m_Y->innerSize() << ',' << m_Y->outerSize() << ')');*/
558
559
        /*model ret(m_Y, m_solver_type);*/
        /*ret.add_bloc(X());*/
560
561
562
563
        ret.add_bloc(m);
        return ret;
    }

564
#if 0
565
    model extend(const std::vector<labelled_matrix<MatrixXd, int, char>>& mv)
566
    {
567
568
569
570
571
572
        /*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);
573
        }
574
575
        return ret;
    }
576
#endif
577
578
579
580
581
582
583
584
585

    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);
586
            } else {
587
                inv_sv(i) = 0.;
588
589
            }
        }
590
        return V * inv_sv.asDiagonal() * V.transpose();
591
592
    }

593
private:
594
595
596
    struct stamps_type {
        bool m_X;
        bool m_rank;
597
        bool m_rss;
598
599
600
        bool m_coefficients;
        bool m_solver;

601
        stamps_type() : m_X(false), m_rank(false), m_rss(false), m_coefficients(false), m_solver(false) {}
602

603
604
        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; }
605
        void update_rank() { m_rank = true; }
606
        void update_solver() { m_solver = true; m_rank = false; m_rss = false; m_coefficients = false; }
607
        void update_coefficients() { m_coefficients = true; }
608
        void update_rss() { m_rss = true; }
609
610
611
612

        bool X_is_uptodate() const { return m_X; }
        bool rank_is_uptodate() const { return m_rank; }
        bool coefficients_is_uptodate() const { return m_coefficients; }
613
        bool rss_is_uptodate() const { return m_rss; }
614
615
        bool solver_is_uptodate() const { return m_solver; }
    };
616
617
618
619
620

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

624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
    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;
    }

641
642
643
    decomposition_base* solver()
    {
        if (!(m_solver != NULL && m_stamps.solver_is_uptodate())) {
644
            return _new_solver();
645
646
647
648
649
        }
        return m_solver;
    }

    struct decomposition_QR : decomposition_base {
650
        /*
651
        ColPivHouseholderQR<MatrixXd> m_solver;
652
653
654
655
        /*/
        // Doesn't seem to perform well when there are many more rows than columns
        FullPivHouseholderQR<MatrixXd> m_solver;
        //*/
656
        int m_columns;
657
        decomposition_QR(const MatrixXd& m) : m_solver(m), m_columns(m.outerSize()) { m_solver.setThreshold(COMPONENT_EPSILON); }
658
        int rank() const { return m_solver.rank(); }
659
        MatrixXd solve(const MatrixXd& lhs)
660
        {
661
662
663
664
            /*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;*/
665
                /*std::cout << "lhs = " << lhs.col(i).transpose() << std::endl;*/
666
                ret.col(i) = m_solver.solve(lhs.col(i));
667
            }
668
            /*std::cout << "QR done" << std::endl;*/
669
670
671
672
673
674
675
676
            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) {}
677
678
679
680
681
682
683
684
685
686
        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;
        }
687
        MatrixXd solve(const MatrixXd& lhs) { return m_solver.solve(lhs); }
688
689
690
        SolverType type() const { return SolverType::SVD; }
    };

691
692
    value<MatrixXd> m_Y;
    collection<model_block_type> m_blocs;
693
694
    MatrixXd m_X;
    int m_rank;
695
    VectorXd m_rss;
696
697
    MatrixXd m_coefficients;
    stamps_type m_stamps;
698
    SolverType m_solver_type;
699
700
    decomposition_base* m_solver;

701
702
703
704
705
706
707
	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;

Damien Leroux's avatar
Damien Leroux committed
708
709
710
711
712
713
714
public:
    friend
    inline
    md5_digest& operator << (md5_digest& md5, const model& m)
    {
        md5 << m.Y() << m.m_blocs;
        return md5;
715
716
    }

Damien Leroux's avatar
Damien Leroux committed
717
718
719
720
    friend
        inline
        std::ostream& operator << (std::ostream& os, const model& m)
        {
721
            os << "<model Y(" << m.m_Y->innerSize() << ',' << m.m_Y->outerSize() << "), " << m.m_blocs.size() << " blocs>";
Damien Leroux's avatar
Damien Leroux committed
722
            return os;
723
724
        }

Damien Leroux's avatar
Damien Leroux committed
725
726
    size_t hash() const
    {
727
        std::hash<model_block_type> hm;
Damien Leroux's avatar
Damien Leroux committed
728
729
730
        size_t accum = 0;
        for (const auto& b: m_blocs) {
            accum ^= hm(b);
731
        }
Damien Leroux's avatar
Damien Leroux committed
732
        return accum;
733
    }
Damien Leroux's avatar
Damien Leroux committed
734
};
735
#endif
Damien Leroux's avatar
Damien Leroux committed
736
737
738
739
740
741

namespace std {
template <>
struct hash<model> {
    size_t operator () (const model& m) { return m.hash(); }
};
742
743
744

}

745
746
#endif