test_model.cc 5.08 KB
Newer Older
1
#include "model.h"
2
#include "chrono.h"
3
4
5
6
7
8
9
10
11
12
13

MatrixXd test_comp(const MatrixXd& M, const MatrixXd& P)
{
    std::cout << "M" << std::endl << M << std::endl;
    std::cout << "P" << std::endl << P << std::endl;
    MatrixXd ret = components(M, P);
    std::cout << "components(M, P)" << std::endl << ret << std::endl;
    return ret;
}


14
15
16
17
bool closeTo(const MatrixXd& a, const MatrixXd& b, double prec = 1.e-10)
{
    return (a - b).isMuchSmallerThan(1., prec) || a.isApprox(b, prec);
}
18
19
20
21


int main(int argc, char** argv)
{
22
#if 0
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
    /*MatrixXd P = MatrixXd::Ones(5, 1);*/
    /*MatrixXd M = MatrixXd::Random(5, 2);*/
    MatrixXd tmp(5, 2);
    tmp << 1, 0,
           1, 3,
           2, 4,
           1, 1,
           0, 1;
    test_comp(tmp, MatrixXd::Ones(5, 1));
    MatrixXd zob(5, 3);
    zob << 0, 1, 2,
           0, 1, 2,
           3, 0, 0,
           0, 1, 3,
           0, 2, 3;
    test_comp(zob, MatrixXd::Ones(5, 1));
    MatrixXd test(5, 2);
    test.col(0) = MatrixXd::Random(5, 1);
    test.col(1) = MatrixXd::Ones(5, 1) - test.col(0);
    test_comp(test, MatrixXd::Ones(5, 1));
    MatrixXd P(5, 1);
    P << 1, 0, 0, 0, 0;
    MatrixXd M(5, 3);
    M << 1, 0, 0,
         0, 1, 0,
         0, 0, 1,
         1, 0, 0,
         0, 0, 1;
    MatrixXd bouh = test_comp(M, P);
    MatrixXd gros = concat_right({bouh, P});
    test_comp(M, gros);

55
56
    model m(tmp);
    /*model m(MatrixXd::Ones(5, 1));*/
57
58
59
60
61
62
63
64
65
66
67
    m.add_bloc(zob);
    std::cout << "MODEL Y" << std::endl << m.Y() << std::endl;
    std::cout << "MODEL X" << std::endl << m.X() << std::endl;
    m.use_SVD();
    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;
    m.use_QR();
    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;
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
#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;*/
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
    {
        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;
    }
146
147
148
    return 0;
}