Belle II Software  release-08-01-10
matrix.test.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 /*
10 This file contains some basic investigation on numerically matrix algorithms.
11 */
12 
13 #include <gtest/gtest.h>
14 
15 #include <vector>
16 #include <Eigen/Dense>
17 
18 
19 // 1. Study of matrix inversions
20 
21 // TL;DR use a rank revealing matrix inversion (usually involving pivoting)
22 // * pivoted householder QR decomposition
23 // * full pivoted LU
24 
25 // no good for inversion if rank of the matrix is not full.
26 // TEST(TrackFindingCDCTest, TMatrixDSym_invert)
27 // {
28 // TMatrixDSym weight(2);
29 // weight(0,0) = 2;
30 // weight(0,1) = 0;
31 // weight(1,0) = 0;
32 // weight(1,1) = 0;
33 
34 // EXPECT_EQ(2, weight.GetNrows());
35 // EXPECT_EQ(2, weight.GetNcols());
36 
37 // TMatrixDSym cov = weight;
38 // cov.Invert();
39 
40 // EXPECT_EQ(2, weight.GetNrows());
41 // EXPECT_EQ(2, weight.GetNcols());
42 
43 // // Already fails here
44 // EXPECT_EQ(2, cov.GetNrows());
45 // EXPECT_EQ(2, cov.GetNcols());
46 
47 // EXPECT_NEAR(1.0 / 2.0, cov(0,0), 10e-7);
48 // EXPECT_NEAR(0.0, cov(0,1), 10e-7);
49 // EXPECT_NEAR(0.0, cov(1,0), 10e-7);
50 // EXPECT_NEAR(0.0, cov(1,1), 10e-7);
51 // }
52 
53 // no good for inversion if rank of the matrix is not full.
54 // TEST(TrackFindingCDCTest, SMatrixDSym_inverse)
55 // {
56 // using CovarianceMatrix = ROOT::Math::SMatrix< double, 2, 2, ROOT::Math::MatRepSym< double, 2>
57 // >;
58 // using WeightMatrix = ROOT::Math::SMatrix< double, 2, 2, ROOT::Math::MatRepSym< double, 2> >;
59 
60 // WeightMatrix weight;
61 // weight(0,0) = 2;
62 // weight(0,1) = 0;
63 // weight(1,0) = 0;
64 // weight(1,1) = 0;
65 
66 // int state = 0;
67 // CovarianceMatrix cov = weight.Inverse(state);
68 
69 // EXPECT_NEAR(1.0 / 2.0, cov(0,0), 10e-7);
70 // EXPECT_NEAR(0.0, cov(0,1), 10e-7);
71 // EXPECT_NEAR(0.0, cov(1,0), 10e-7);
72 // EXPECT_NEAR(0.0, cov(1,1), 10e-7);
73 // }
74 
75 // no good for inversion if rank of the matrix is not full.
76 // TEST(TrackFindingCDCTest, SMatrixDSym_inverseChol)
77 // {
78 // using CovarianceMatrix = ROOT::Math::SMatrix< double, 2, 2, ROOT::Math::MatRepSym< double, 2>
79 // >;
80 // using WeightMatrix = ROOT::Math::SMatrix< double, 2, 2, ROOT::Math::MatRepSym< double, 2> >;
81 
82 // WeightMatrix weight;
83 // weight(0,0) = 2;
84 // weight(0,1) = 0;
85 // weight(1,0) = 0;
86 // weight(1,1) = 0;
87 
88 // int state = 0;
89 // CovarianceMatrix cov = weight.InverseChol(state);
90 
91 // EXPECT_NEAR(1.0 / 2.0, cov(0,0), 10e-7);
92 // EXPECT_NEAR(0.0, cov(0,1), 10e-7);
93 // EXPECT_NEAR(0.0, cov(1,0), 10e-7);
94 // EXPECT_NEAR(0.0, cov(1,1), 10e-7);
95 // }
96 
97 // no good for inversion if rank of the matrix is not full.
98 // TEST(TrackFindingCDCTest, Eigen_inverse)
99 // {
100 // Eigen::Matrix<double, 2, 2> weight;
101 // weight(0,0) = 2;
102 // weight(0,1) = 0;
103 // weight(1,0) = 0;
104 // weight(1,1) = 0;
105 
106 // Eigen::Matrix<double, 2, 2> cov = weight.inverse();
107 
108 // EXPECT_NEAR(1.0 / 2.0, cov(0,0), 10e-7);
109 // EXPECT_NEAR(0.0, cov(0,1), 10e-7);
110 // EXPECT_NEAR(0.0, cov(1,0), 10e-7);
111 // EXPECT_NEAR(0.0, cov(1,1), 10e-7);
112 // }
113 
114 // no good for inversion if rank of the matrix is not full.
115 // TEST(TrackFindingCDCTest, Eigen_llt_inverse)
116 // {
117 // Eigen::Matrix<double, 2, 2> weight;
118 // weight(0,0) = 2;
119 // weight(0,1) = 0;
120 // weight(1,0) = 0;
121 // weight(1,1) = 0;
122 
123 // Eigen::Matrix<double, 2, 2> cov = weight.llt().inverse();
124 
125 // EXPECT_NEAR(1.0 / 2.0, cov(0,0), 10e-7);
126 // EXPECT_NEAR(0.0, cov(0,1), 10e-7);
127 // EXPECT_NEAR(0.0, cov(1,0), 10e-7);
128 // EXPECT_NEAR(0.0, cov(1,1), 10e-7);
129 // }
130 
131 // No good for inversion if rank of the matrix is not full.
132 // TEST(TrackFindingCDCTest, Eigen_lldt_inverse)
133 // {
134 // Eigen::Matrix<double, 2, 2> weight;
135 // weight(0,0) = 2;
136 // weight(0,1) = 0;
137 // weight(1,0) = 0;
138 // weight(1,1) = 0;
139 
140 // Eigen::Matrix<double, 2, 2> cov = weight.lldt().inverse();
141 
142 // EXPECT_NEAR(1.0 / 2.0, cov(0,0), 10e-7);
143 // EXPECT_NEAR(0.0, cov(0,1), 10e-7);
144 // EXPECT_NEAR(0.0, cov(1,0), 10e-7);
145 // EXPECT_NEAR(0.0, cov(1,1), 10e-7);
146 // }
147 
148 TEST(TrackFindingCDCTest, Eigen_fullPivLu_inverse)
149 {
150  Eigen::Matrix<double, 2, 2> weight;
151  weight(0, 0) = 2;
152  weight(0, 1) = 0;
153  weight(1, 0) = 0;
154  weight(1, 1) = 0;
155 
156  Eigen::Matrix<double, 2, 2> cov = weight.fullPivLu().inverse();
157 
158  EXPECT_NEAR(1.0 / 2.0, cov(0, 0), 10e-7);
159  EXPECT_NEAR(0.0, cov(0, 1), 10e-7);
160  EXPECT_NEAR(0.0, cov(1, 0), 10e-7);
161  EXPECT_NEAR(0.0, cov(1, 1), 10e-7);
162 }
163 
164 // No good for inversion if rank of the matrix is not full.
165 // TEST(TrackFindingCDCTest, Eigen_partialPivLu_inverse)
166 // {
167 // Eigen::Matrix<double, 2, 2> weight;
168 // weight(0,0) = 2;
169 // weight(0,1) = 0;
170 // weight(1,0) = 0;
171 // weight(1,1) = 0;
172 
173 // Eigen::Matrix<double, 2, 2> cov = weight.partialPivLu().inverse();
174 
175 // EXPECT_NEAR(1.0 / 2.0, cov(0,0), 10e-7);
176 // EXPECT_NEAR(0.0, cov(0,1), 10e-7);
177 // EXPECT_NEAR(0.0, cov(1,0), 10e-7);
178 // EXPECT_NEAR(0.0, cov(1,1), 10e-7);
179 // }
180 
181 // Can not do inversion
182 // TEST(TrackFindingCDCTest, Eigen_householderQr_inverse)
183 // {
184 // Eigen::Matrix<double, 2, 2> weight;
185 // weight(0,0) = 2;
186 // weight(0,1) = 0;
187 // weight(1,0) = 0;
188 // weight(1,1) = 0;
189 
190 // Eigen::Matrix<double, 2, 2> cov = weight.householderQr().inverse();
191 
192 // EXPECT_NEAR(1.0 / 2.0, cov(0,0), 10e-7);
193 // EXPECT_NEAR(0.0, cov(0,1), 10e-7);
194 // EXPECT_NEAR(0.0, cov(1,0), 10e-7);
195 // EXPECT_NEAR(0.0, cov(1,1), 10e-7);
196 // }
197 
198 /* Method of choice for covariance <-> weight matrices */
199 TEST(TrackFindingCDCTest, Eigen_colPivHouseholderQr_inverse)
200 {
201  Eigen::Matrix<double, 2, 2> weight;
202  weight(0, 0) = 2;
203  weight(0, 1) = 0;
204  weight(1, 0) = 0;
205  weight(1, 1) = 0;
206 
207  Eigen::Matrix<double, 2, 2> cov = weight.colPivHouseholderQr().inverse();
208 
209  EXPECT_NEAR(1.0 / 2.0, cov(0, 0), 10e-7);
210  EXPECT_NEAR(0.0, cov(0, 1), 10e-7);
211  EXPECT_NEAR(0.0, cov(1, 0), 10e-7);
212  EXPECT_NEAR(0.0, cov(1, 1), 10e-7);
213 }
214 
215 TEST(TrackFindingCDCTest, Eigen_fullPivHouseholderQr_inverse)
216 {
217  Eigen::Matrix<double, 2, 2> weight;
218  weight(0, 0) = 2;
219  weight(0, 1) = 0;
220  weight(1, 0) = 0;
221  weight(1, 1) = 0;
222 
223  Eigen::Matrix<double, 2, 2> cov = weight.fullPivHouseholderQr().inverse();
224 
225  EXPECT_NEAR(1.0 / 2.0, cov(0, 0), 10e-7);
226  EXPECT_NEAR(0.0, cov(0, 1), 10e-7);
227  EXPECT_NEAR(0.0, cov(1, 0), 10e-7);
228  EXPECT_NEAR(0.0, cov(1, 1), 10e-7);
229 }
Belle2::TEST
TEST(TestgetDetectorRegion, TestgetDetectorRegion)
Test Constructors.
Definition: utilityFunctions.cc:25