Rivet  1.8.3
MatrixN.hh
1 #ifndef RIVET_MATH_MATRIXN
2 #define RIVET_MATH_MATRIXN
3 
4 #include "Rivet/Math/MathHeader.hh"
5 #include "Rivet/Math/MathUtils.hh"
6 #include "Rivet/Math/Vectors.hh"
7 
8 #include "Rivet/Math/eigen/matrix.h"
9 
10 namespace Rivet {
11 
12 
13  template <size_t N>
14  class Matrix;
15  typedef Matrix<4> Matrix4;
16 
17  template <size_t N>
18  Matrix<N> multiply(const Matrix<N>& a, const Matrix<N>& b);
19  template <size_t N>
20  Matrix<N> divide(const Matrix<N>&, const double);
21  template <size_t N>
22  Matrix<N> operator*(const Matrix<N>& a, const Matrix<N>& b);
23 
24 
26 
27 
29  template <size_t N>
30  class Matrix {
31  template <size_t M>
32  friend Matrix<M> add(const Matrix<M>&, const Matrix<M>&);
33  template <size_t M>
34  friend Matrix<M> multiply(const double, const Matrix<M>&);
35  template <size_t M>
36  friend Matrix<M> multiply(const Matrix<M>&, const Matrix<M>&);
37  template <size_t M>
38  friend Vector<M> multiply(const Matrix<M>&, const Vector<M>&);
39  template <size_t M>
40  friend Matrix<M> divide(const Matrix<M>&, const double);
41 
42  public:
43  static Matrix<N> mkZero() {
44  Matrix<N> rtn;
45  return rtn;
46  }
47 
48  static Matrix<N> mkDiag(Vector<N> diag) {
49  Matrix<N> rtn;
50  for (size_t i = 0; i < N; ++i) {
51  rtn.set(i, i, diag[i]);
52  }
53  return rtn;
54  }
55 
56  static Matrix<N> mkIdentity() {
57  Matrix<N> rtn;
58  for (size_t i = 0; i < N; ++i) {
59  rtn.set(i, i, 1);
60  }
61  return rtn;
62  }
63 
64 
65  public:
66 
67  Matrix() {
68  _matrix.loadZero();
69  }
70 
71  Matrix(const Matrix<N>& other) {
72  _matrix = other._matrix;
73  }
74 
75  Matrix& set(const size_t i, const size_t j, const double value) {
76  if (i < N && j < N) {
77  _matrix(i, j) = value;
78  } else {
79  throw std::runtime_error("Attempted set access outside matrix bounds.");
80  }
81  return *this;
82  }
83 
84  double get(const size_t i, const size_t j) const {
85  if (i < N && j < N) {
86  return _matrix(i, j);
87  } else {
88  throw std::runtime_error("Attempted get access outside matrix bounds.");
89  }
90  }
91 
92  Vector<N> getRow(const size_t row) const {
93  Vector<N> rtn;
94  for (size_t i = 0; i < N; ++i) {
95  rtn.set(i, _matrix(row,i));
96  }
97  return rtn;
98  }
99 
100  Matrix<N>& setRow(const size_t row, const Vector<N>& r) {
101  for (size_t i = 0; i < N; ++i) {
102  _matrix(row,i) = r.get(i);
103  }
104  return *this;
105  }
106 
107  Vector<N> getColumn(const size_t col) const {
108  Vector<N> rtn;
109  for (size_t i = 0; i < N; ++i) {
110  rtn.set(i, _matrix(i,col));
111  }
112  return rtn;
113  }
114 
115  Matrix<N>& setColumn(const size_t col, const Vector<N>& c) {
116  for (size_t i = 0; i < N; ++i) {
117  _matrix(i,col) = c.get(i);
118  }
119  return *this;
120  }
121 
122  Matrix<N> transpose() const {
123  Matrix<N> tmp = *this;
124  tmp._matrix.replaceWithAdjoint();
125  return tmp;
126  }
127 
128  // Matrix<N>& transposeInPlace() {
129  // _matrix.replaceWithAdjoint();
130  // return *this;
131  // }
132 
134  Matrix<N> inverse() const {
135  Matrix<N> tmp;
136  tmp._matrix = _matrix.inverse();
137  return tmp;
138  }
139 
141  double det() const {
142  return _matrix.determinant();
143  }
144 
146  double trace() const {
147  double tr = 0.0;
148  for (size_t i = 0; i < N; ++i) {
149  tr += _matrix(i,i);
150  }
151  return tr;
152  // return _matrix.trace();
153  }
154 
157  Matrix<N> rtn;
158  rtn._matrix = -_matrix;
159  return rtn;
160  }
161 
163  size_t size() const {
164  return N;
165  }
166 
168  bool isZero(double tolerance=1E-5) const {
169  for (size_t i=0; i < N; ++i) {
170  for (size_t j=0; j < N; ++j) {
171  if (! Rivet::isZero(_matrix(i,j), tolerance) ) return false;
172  }
173  }
174  return true;
175  }
176 
178  bool isEqual(Matrix<N> other) const {
179  for (size_t i=0; i < N; ++i) {
180  for (size_t j=i; j < N; ++j) {
181  if (! Rivet::isZero(_matrix(i,j) - other._matrix(i,j)) ) return false;
182  }
183  }
184  return true;
185  }
186 
188  bool isSymm() const {
189  return isEqual(this->transpose());
190  }
191 
193  bool isDiag() const {
194  for (size_t i=0; i < N; ++i) {
195  for (size_t j=0; j < N; ++j) {
196  if (i == j) continue;
197  if (! Rivet::isZero(_matrix(i,j)) ) return false;
198  }
199  }
200  return true;
201  }
202 
203  bool operator==(const Matrix<N>& a) const {
204  return _matrix == a._matrix;
205  }
206 
207  bool operator!=(const Matrix<N>& a) const {
208  return _matrix != a._matrix;
209  }
210 
211  bool operator<(const Matrix<N>& a) const {
212  return _matrix < a._matrix;
213  }
214 
215  bool operator<=(const Matrix<N>& a) const {
216  return _matrix <= a._matrix;
217  }
218 
219  bool operator>(const Matrix<N>& a) const {
220  return _matrix > a._matrix;
221  }
222 
223  bool operator>=(const Matrix<N>& a) const {
224  return _matrix >= a._matrix;
225  }
226 
227  Matrix<N>& operator*=(const Matrix<N>& m) {
228  _matrix = _matrix * m._matrix;
229  return *this;
230  }
231 
232  Matrix<N>& operator*=(const double a) {
233  _matrix *= a;
234  return *this;
235  }
236 
237  Matrix<N>& operator/=(const double a) {
238  _matrix /= a;
239  return *this;
240  }
241 
242  Matrix<N>& operator+=(const Matrix<N>& m) {
243  _matrix += m._matrix;
244  return *this;
245  }
246 
247  Matrix<N>& operator-=(const Matrix<N>& m) {
248  _matrix -= m._matrix;
249  return *this;
250  }
251 
252  protected:
253  typedef Eigen::Matrix<double,N> EMatrix;
254  EMatrix _matrix;
255  };
256 
257 
259 
260 
261  template <size_t N>
262  inline Matrix<N> add(const Matrix<N>& a, const Matrix<N>& b) {
263  Matrix<N> result;
264  result._matrix = a._matrix + b._matrix;
265  return result;
266  }
267 
268  template <size_t N>
269  inline Matrix<N> subtract(const Matrix<N>& a, const Matrix<N>& b) {
270  return add(a, -b);
271  }
272 
273  template <size_t N>
274  inline Matrix<N> operator+(const Matrix<N> a, const Matrix<N>& b) {
275  return add(a, b);
276  }
277 
278  template <size_t N>
279  inline Matrix<N> operator-(const Matrix<N> a, const Matrix<N>& b) {
280  return subtract(a, b);
281  }
282 
283  template <size_t N>
284  inline Matrix<N> multiply(const double a, const Matrix<N>& m) {
285  Matrix<N> rtn;
286  rtn._matrix = a * m._matrix;
287  return rtn;
288  }
289 
290  template <size_t N>
291  inline Matrix<N> multiply(const Matrix<N>& m, const double a) {
292  return multiply(a, m);
293  }
294 
295  template <size_t N>
296  inline Matrix<N> divide(const Matrix<N>& m, const double a) {
297  return multiply(1/a, m);
298  }
299 
300  template <size_t N>
301  inline Matrix<N> operator*(const double a, const Matrix<N>& m) {
302  return multiply(a, m);
303  }
304 
305  template <size_t N>
306  inline Matrix<N> operator*(const Matrix<N>& m, const double a) {
307  return multiply(a, m);
308  }
309 
310  template <size_t N>
311  inline Matrix<N> multiply(const Matrix<N>& a, const Matrix<N>& b) {
312  Matrix<N> tmp;
313  tmp._matrix = a._matrix * b._matrix;
314  return tmp;
315  }
316 
317  template <size_t N>
318  inline Matrix<N> operator*(const Matrix<N>& a, const Matrix<N>& b) {
319  return multiply(a, b);
320  }
321 
322 
323  template <size_t N>
324  inline Vector<N> multiply(const Matrix<N>& a, const Vector<N>& b) {
325  Vector<N> tmp;
326  tmp._vec = a._matrix * b._vec;
327  return tmp;
328  }
329 
330  template <size_t N>
331  inline Vector<N> operator*(const Matrix<N>& a, const Vector<N>& b) {
332  return multiply(a, b);
333  }
334 
335  template <size_t N>
336  inline Matrix<N> transpose(const Matrix<N>& m) {
337  // Matrix<N> tmp;
338  // for (size_t i = 0; i < N; ++i) {
339  // for (size_t j = 0; j < N; ++j) {
340  // tmp.set(i, j, m.get(j, i));
341  // }
342  // }
343  // return tmp;
344  return m.transpose();
345  }
346 
347  template <size_t N>
348  inline Matrix<N> inverse(const Matrix<N>& m) {
349  return m.inverse();
350  }
351 
352  template <size_t N>
353  inline double det(const Matrix<N>& m) {
354  return m.determinant();
355  }
356 
357  template <size_t N>
358  inline double trace(const Matrix<N>& m) {
359  return m.trace();
360  }
361 
362 
364 
365 
367  template <size_t N>
368  inline string toString(const Matrix<N>& m) {
369  ostringstream ss;
370  ss << "[ ";
371  for (size_t i = 0; i < m.size(); ++i) {
372  ss << "( ";
373  for (size_t j = 0; j < m.size(); ++j) {
374  const double e = m.get(i, j);
375  ss << (Rivet::isZero(e) ? 0.0 : e) << " ";
376  }
377  ss << ") ";
378  }
379  ss << "]";
380  return ss.str();
381  }
382 
383 
385  template <size_t N>
386  inline ostream& operator<<(std::ostream& out, const Matrix<N>& m) {
387  out << toString(m);
388  return out;
389  }
390 
391 
393 
394 
396  template <size_t N>
397  inline bool fuzzyEquals(const Matrix<N>& ma, const Matrix<N>& mb, double tolerance=1E-5) {
398  for (size_t i = 0; i < N; ++i) {
399  for (size_t j = 0; j < N; ++j) {
400  const double a = ma.get(i, j);
401  const double b = mb.get(i, j);
402  if (!Rivet::fuzzyEquals(a, b, tolerance)) return false;
403  }
404  }
405  return true;
406  }
407 
408 
410  template <size_t N>
411  inline bool isZero(const Matrix<N>& m, double tolerance=1E-5) {
412  return m.isZero(tolerance);
413  }
414 
415 
416 }
417 
418 #endif
Definition: MC_JetAnalysis.hh:9
bool isZero(double tolerance=1E-5) const
Index-wise check for nullness, allowing for numerical precision.
Definition: MatrixN.hh:168
Matrix< N > operator-() const
Negate.
Definition: MatrixN.hh:156
bool isSymm() const
Check for symmetry under transposition.
Definition: MatrixN.hh:188
bool isDiag() const
Check that all off-diagonal elements are zero, allowing for numerical precision.
Definition: MatrixN.hh:193
double det() const
Calculate determinant.
Definition: MatrixN.hh:141
bool isEqual(Matrix< N > other) const
Check for index-wise equality, allowing for numerical precision.
Definition: MatrixN.hh:178
Matrix< N > inverse() const
Calculate inverse.
Definition: MatrixN.hh:134
double trace() const
Calculate trace.
Definition: MatrixN.hh:146
Vector< N > & set(const size_t index, const double value)
Set indexed value.
Definition: VectorN.hh:52
General -dimensional mathematical matrix object.
Definition: MatrixN.hh:14
bool isZero(double val, double tolerance=1E-8)
Definition: MathUtils.hh:17
A minimal base class for -dimensional vectors.
Definition: VectorN.hh:13
size_t size() const
Get dimensionality.
Definition: MatrixN.hh:163
std::string toString(const AnalysisInfo &ai)
String representation.
Definition: AnalysisInfo.cc:234
bool fuzzyEquals(double a, double b, double tolerance=1E-5)
Compare two floating point numbers for equality with a degree of fuzziness.
Definition: MathUtils.hh:34