Functions
LATest.cc File Reference
#include <iostream>
#include <cassert>
#include <cmath>
#include "lardata/RecoObjects/KalmanLinearAlgebra.h"
#include "boost/numeric/ublas/io.hpp"

Go to the source code of this file.

Functions

int main ()
 

Function Documentation

int main ( void  )

Definition at line 13 of file LATest.cc.

14 {
15  // Make sure assert is enabled.
16 
17  bool assert_flag = false;
18  assert((assert_flag = true, assert_flag));
19  if ( ! assert_flag ) {
20  std::cerr << "Assert is disabled" << std::endl;
21  return 1;
22  }
23 
24  // Test matrix inversion.
25 
26  // 1x1 KSymMatrix.
27 
29  for(unsigned int i = 0; i < m1.size1(); ++i) {
30  for(unsigned int j = 0; j <= i; ++j)
31  m1(i,j) = i+j+1.;
32  }
33  trkf::KSymMatrix<1>::type minv1(m1);
34  bool ok = trkf::syminvert(minv1);
35  assert(ok);
36  trkf::ublas::matrix<double> unit1 = prod(m1, minv1);
37  std::cout << m1 << std::endl;
38  std::cout << minv1 << std::endl;
39  std::cout << unit1 << std::endl;
40  for(unsigned int i = 0; i < m1.size1(); ++i) {
41  for(unsigned int j = 0; j <= i; ++j) {
42  double val = (i==j ? 1. : 0.);
43  assert(std::abs(unit1(i,j) - val) < 1.e-10);
44  }
45  }
46 
47  // 2x2 KSymMatrix.
48 
50  for(unsigned int i = 0; i < m2.size1(); ++i) {
51  for(unsigned int j = 0; j <= i; ++j)
52  m2(i,j) = i+j+1.;
53  }
55  ok = trkf::syminvert(minv2);
56  assert(ok);
57  trkf::ublas::matrix<double> unit2 = prod(m2, minv2);
58  std::cout << m2 << std::endl;
59  std::cout << minv2 << std::endl;
60  std::cout << unit2 << std::endl;
61  for(unsigned int i = 0; i < m2.size1(); ++i) {
62  for(unsigned int j = 0; j <= i; ++j) {
63  double val = (i==j ? 1. : 0.);
64  assert(std::abs(unit2(i,j) - val) < 1.e-10);
65  }
66  }
67 
68  // 3x3 KSymMatrix.
69 
71  for(unsigned int i = 0; i < m3.size1(); ++i) {
72  for(unsigned int j = 0; j <= i; ++j) {
73  m3(i,j) = i+j+1.;
74  if(i==j)
75  m3(i,j) += 1.;
76  }
77  }
79  ok = trkf::syminvert(minv3);
80  assert(ok);
81  trkf::ublas::matrix<double> unit3 = prod(m3, minv3);
82  std::cout << m3 << std::endl;
83  std::cout << minv3 << std::endl;
84  std::cout << unit3 << std::endl;
85  for(unsigned int i = 0; i < m3.size1(); ++i) {
86  for(unsigned int j = 0; j <= i; ++j) {
87  double val = (i==j ? 1. : 0.);
88  assert(std::abs(unit3(i,j) - val) < 1.e-10);
89  }
90  }
91 
92  // 4x4 KSymMatrix.
93 
95  for(unsigned int i = 0; i < m4.size1(); ++i) {
96  for(unsigned int j = 0; j <= i; ++j) {
97  m4(i,j) = i+j+1.;
98  if(i==j)
99  m4(i,j) += 1.;
100  }
101  }
102  trkf::KSymMatrix<4>::type minv4(m4);
103  ok = trkf::syminvert(minv4);
104  assert(ok);
105  trkf::ublas::matrix<double> unit4 = prod(m4, minv4);
106  std::cout << m4 << std::endl;
107  std::cout << minv4 << std::endl;
108  std::cout << unit4 << std::endl;
109  for(unsigned int i = 0; i < m4.size1(); ++i) {
110  for(unsigned int j = 0; j <= i; ++j) {
111  double val = (i==j ? 1. : 0.);
112  assert(std::abs(unit4(i,j) - val) < 1.e-10);
113  }
114  }
115 
116  // 5x5 TrackError.
117 
118  trkf::TrackError m5;
119  for(unsigned int i = 0; i < m5.size1(); ++i) {
120  for(unsigned int j = 0; j <= i; ++j) {
121  m5(i,j) = i+j+1.;
122  if(i==j)
123  m5(i,j) += 1.;
124  }
125  }
126  trkf::TrackError minv5(m5);
127  ok = trkf::syminvert(minv5);
128  assert(ok);
129  trkf::ublas::matrix<double> unit5 = prod(m5, minv5);
130  std::cout << m5 << std::endl;
131  std::cout << minv5 << std::endl;
132  std::cout << unit5 << std::endl;
133  for(unsigned int i = 0; i < m5.size1(); ++i) {
134  for(unsigned int j = 0; j <= i; ++j) {
135  double val = (i==j ? 1. : 0.);
136  assert(std::abs(unit5(i,j) - val) < 1.e-10);
137  }
138  }
139 
140  // 1x1 KMatrix.
141 
142  trkf::KMatrix<1,1>::type m6(1,1);
143  for(unsigned int i = 0; i < m6.size1(); ++i) {
144  for(unsigned int j = 0; j < m6.size2(); ++j)
145  m6(i,j) = i + 2*j + 1.;
146  }
147  trkf::KMatrix<1,1>::type minv6(m6);
148  ok = trkf::invert(minv6);
149  assert(ok);
150  trkf::ublas::matrix<double> unit6 = prod(m6, minv6);
151  std::cout << m6 << std::endl;
152  std::cout << minv6 << std::endl;
153  std::cout << unit6 << std::endl;
154  for(unsigned int i = 0; i < m6.size1(); ++i) {
155  for(unsigned int j = 0; j < m6.size2(); ++j) {
156  double val = (i==j ? 1. : 0.);
157  assert(std::abs(unit6(i,j) - val) < 1.e-10);
158  }
159  }
160 
161  // 2x2 KMatrix.
162 
163  trkf::KMatrix<2,2>::type m7(2,2);
164  for(unsigned int i = 0; i < m7.size1(); ++i) {
165  for(unsigned int j = 0; j < m7.size2(); ++j)
166  m7(i,j) = i + 2*j + 1.;
167  }
168  trkf::KMatrix<2,2>::type minv7(m7);
169  ok = trkf::invert(minv7);
170  assert(ok);
171  trkf::ublas::matrix<double> unit7 = prod(m7, minv7);
172  std::cout << m7 << std::endl;
173  std::cout << minv7 << std::endl;
174  std::cout << unit7 << std::endl;
175  for(unsigned int i = 0; i < m7.size1(); ++i) {
176  for(unsigned int j = 0; j < m7.size2(); ++j) {
177  double val = (i==j ? 1. : 0.);
178  assert(std::abs(unit7(i,j) - val) < 1.e-10);
179  }
180  }
181 
182  // 3x3 KMatrix.
183 
184  trkf::KMatrix<3,3>::type m8(3,3);
185  for(unsigned int i = 0; i < m8.size1(); ++i) {
186  for(unsigned int j = 0; j < m8.size2(); ++j)
187  m8(i,j) = i + 2*j + (i==j ? 1. : 0.);
188 
189  }
190  trkf::KMatrix<3,3>::type minv8(m8);
191  ok = trkf::invert(minv8);
192  assert(ok);
193  trkf::ublas::matrix<double> unit8 = prod(m8, minv8);
194  std::cout << m8 << std::endl;
195  std::cout << minv8 << std::endl;
196  std::cout << unit8 << std::endl;
197  for(unsigned int i = 0; i < m8.size1(); ++i) {
198  for(unsigned int j = 0; j < m8.size2(); ++j) {
199  double val = (i==j ? 1. : 0.);
200  assert(std::abs(unit8(i,j) - val) < 1.e-10);
201  }
202  }
203 
204  // 4x4 KMatrix.
205 
206  trkf::KMatrix<4,4>::type m9(4,4);
207  for(unsigned int i = 0; i < m9.size1(); ++i) {
208  for(unsigned int j = 0; j < m9.size2(); ++j)
209  m9(i,j) = i + 2*j + (i==j ? 1. : 0.);
210 
211  }
212  trkf::KMatrix<4,4>::type minv9(m9);
213  ok = trkf::invert(minv9);
214  assert(ok);
215  trkf::ublas::matrix<double> unit9 = prod(m9, minv9);
216  std::cout << m9 << std::endl;
217  std::cout << minv9 << std::endl;
218  std::cout << unit9 << std::endl;
219  for(unsigned int i = 0; i < m9.size1(); ++i) {
220  for(unsigned int j = 0; j < m9.size2(); ++j) {
221  double val = (i==j ? 1. : 0.);
222  assert(std::abs(unit9(i,j) - val) < 1.e-10);
223  }
224  }
225 
226  // 5x5 TrackMatrix.
227 
228  trkf::TrackMatrix m10(5,5);
229  for(unsigned int i = 0; i < m10.size1(); ++i) {
230  for(unsigned int j = 0; j < m10.size2(); ++j)
231  m10(i,j) = i + 2*j + (i==j ? 1. : 0.);
232 
233  }
234  trkf::TrackMatrix minv10(m10);
235  ok = trkf::invert(minv10);
236  assert(ok);
237  trkf::ublas::matrix<double> unit10 = prod(m10, minv10);
238  std::cout << m10 << std::endl;
239  std::cout << minv10 << std::endl;
240  std::cout << unit10 << std::endl;
241  for(unsigned int i = 0; i < m10.size1(); ++i) {
242  for(unsigned int j = 0; j < m10.size2(); ++j) {
243  double val = (i==j ? 1. : 0.);
244  assert(std::abs(unit10(i,j) - val) < 1.e-10);
245  }
246  }
247 
248  // Done (success).
249 
250  std::cout << "LATest: All tests passed." << std::endl;
251 
252  return 0;
253 }
ublas::symmetric_matrix< double, ublas::lower, ublas::row_major, ublas::bounded_array< double, N *(N+1)/2 > > type
ublas::matrix< double, ublas::row_major, ublas::bounded_array< double, N *M > > type
KSymMatrix< 5 >::type TrackError
Track error matrix, dimension 5x5.
T abs(T value)
const double e
static constexpr double m2
Definition: Units.h:72
bool syminvert(ublas::symmetric_matrix< T, TRI, L, A > &m)
bool invert(ublas::matrix< T, L, A > &m)
static constexpr double m3
Definition: Units.h:73
KMatrix< 5, 5 >::type TrackMatrix
General 5x5 matrix.
QTextStream & endl(QTextStream &s)