Tulip  5.7.4
Large graphs analysis and drawing
Matrix.cxx
1 /*
2  *
3  * This file is part of Tulip (https://tulip.labri.fr)
4  *
5  * Authors: David Auber and the Tulip development Team
6  * from LaBRI, University of Bordeaux
7  *
8  * Tulip is free software; you can redistribute it and/or modify
9  * it under the terms of the GNU Lesser General Public License
10  * as published by the Free Software Foundation, either version 3
11  * of the License, or (at your option) any later version.
12  *
13  * Tulip is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
16  * See the GNU General Public License for more details.
17  *
18  */
19 #include <cmath>
20 
21 //===================================================================
22 // Specialisation
23 namespace tlp {
24 template <typename Obj>
25 class Matrix<Obj, 1> : public Array<Vector<Obj, 1>, 1> {
26 public:
27  Obj determinant() {
28  return (*this)[0][0];
29  }
30  // Matrix<Obj,1>& fill(Obj obj) {return *this;}
31  Matrix<Obj, 1> &inverse() {
32  (*this)[0][0] = 1.0 / (*this)[0][0];
33  return *this;
34  }
35  Matrix<Obj, 1> &transpose() {
36  return *this;
37  }
38  Matrix<Obj, 1> &operator*=(const Matrix<Obj, 1> &mat) {
39  (*this)[0][0] *= mat[0][0];
40  return *this;
41  }
42  // Matrix<Obj,1>& operator/=(const Obj &obj){return *this;}
43  Matrix<Obj, 1> cofactor() {
44  return *this;
45  }
46 };
47 } // namespace tlp
48 
49 //===================================================================
50 template <typename Obj, size_t SIZE>
51 MATRIX::Matrix(const std::vector<std::vector<Obj>> &covarianceMatrix) {
52  for (size_t i = 0; i < SIZE; i++)
53  for (size_t j = 0; j < SIZE; j++)
54  (*this)[i][j] =
55  covarianceMatrix[i][j] / (sqrt(covarianceMatrix[i][i] * covarianceMatrix[j][j]));
56 }
57 
58 //===================================================================
59 template <typename Obj, size_t SIZE>
60 MATRIX &MATRIX::fill(Obj obj) {
61  for (size_t i = 0; i < SIZE; ++i)
62  (*this)[i].fill(obj);
63 
64  return (*this);
65 }
66 //======================================================
67 template <typename Obj, size_t SIZE>
68 MATRIX &MATRIX::operator+=(const MATRIX &mat) {
69  for (size_t i = 0; i < SIZE; ++i)
70  (*this)[i] += mat[i];
71 
72  return (*this);
73 }
74 //======================================================
75 template <typename Obj, size_t SIZE>
76 MATRIX &MATRIX::operator-=(const MATRIX &mat) {
77  for (size_t i = 0; i < SIZE; ++i)
78  (*this)[i] -= mat[i];
79 
80  return (*this);
81 }
82 //======================================================
83 template <typename Obj, size_t SIZE>
84 bool MATRIX::operator==(const MATRIX &mat) const {
85  for (size_t i = 0; i < SIZE; ++i) {
86  if (((*this)[i] != mat[i]))
87  return false;
88  }
89 
90  return true;
91 }
92 //======================================================
93 template <typename Obj, size_t SIZE>
94 bool MATRIX::operator!=(const MATRIX &mat) const {
95  for (size_t i = 0; i < SIZE; ++i) {
96  if (((*this)[i] != mat[i]))
97  return true;
98  }
99 
100  return false;
101 }
102 //===================================================================
103 template <typename Obj, size_t SIZE>
104 MATRIX &MATRIX::operator*=(const MATRIX &mat) {
105  (*this) = (*this) * mat;
106  return (*this);
107 }
108 //=====================================================================================
109 template <typename Obj, size_t SIZE>
110 MATRIX &MATRIX::operator*=(const Obj &obj) {
111  for (size_t i = 0; i < SIZE; ++i)
112  (*this)[i] *= obj;
113 
114  return (*this);
115 }
116 //=====================================================================================
117 template <typename Obj, size_t SIZE>
118 MATRIX &MATRIX::operator/=(const MATRIX &mat) {
119  MATRIX tmpMat(mat);
120  tmpMat.inverse();
121  (*this) *= tmpMat;
122  return (*this);
123 }
124 //=====================================================================================
125 template <typename Obj, size_t SIZE>
126 MATRIX &MATRIX::operator/=(const Obj &obj) {
127  for (size_t i = 0; i < SIZE; ++i)
128  (*this)[i] /= obj;
129 
130  return (*this);
131 }
132 //=====================================================================================
133 template <typename Obj, size_t SIZE>
134 Obj MATRIX::determinant() const {
135  switch (SIZE) {
136  case 2:
137  return (*this)[0][0] * (*this)[1][1] - (*this)[1][0] * (*this)[0][1];
138  break;
139 
140  case 3:
141  return (*this)[0][0] * ((*this)[1][1] * (*this)[2][2] - (*this)[1][2] * (*this)[2][1]) -
142  (*this)[0][1] * ((*this)[1][0] * (*this)[2][2] - (*this)[1][2] * (*this)[2][0]) +
143  (*this)[0][2] * ((*this)[1][0] * (*this)[2][1] - (*this)[1][1] * (*this)[2][0]);
144  break;
145 
146  default:
147  int j2;
148  Obj det = 0;
149 
150  for (size_t j1 = 0; j1 < SIZE; ++j1) {
151  tlp::Matrix<Obj, SIZE - 1> m;
152 
153  for (size_t i = 1; i < SIZE; i++) {
154  j2 = 0;
155 
156  for (size_t j = 0; j < SIZE; ++j) {
157  if (j == j1)
158  continue;
159 
160  m[i - 1][j2] = (*this)[i][j];
161  ++j2;
162  }
163  }
164 
165  if (j1 & 1)
166  det += (*this)[0][j1] * m.determinant();
167  else
168  det -= (*this)[0][j1] * m.determinant();
169  }
170 
171  return (det);
172  }
173 }
174 //=====================================================================================
175 template <typename Obj, size_t SIZE>
176 MATRIX MATRIX::cofactor() const {
177  MATRIX result;
178 
179  switch (SIZE) {
180  case 2:
181  (result)[0][0] = (*this)[1][1];
182  (result)[0][1] = -(*this)[1][0];
183  (result)[1][0] = -(*this)[0][1];
184  (result)[1][1] = (*this)[0][0];
185  break;
186 
187  case 3:
188  (result)[0][0] = (*this)[1][1] * (*this)[2][2] - (*this)[1][2] * (*this)[2][1];
189  (result)[0][1] = -((*this)[1][0] * (*this)[2][2] - (*this)[2][0] * (*this)[1][2]);
190  (result)[0][2] = (*this)[1][0] * (*this)[2][1] - (*this)[1][1] * (*this)[2][0];
191  (result)[1][0] = -((*this)[0][1] * (*this)[2][2] - (*this)[0][2] * (*this)[2][1]);
192  (result)[1][1] = (*this)[0][0] * (*this)[2][2] - (*this)[0][2] * (*this)[2][0];
193  (result)[1][2] = -((*this)[0][0] * (*this)[2][1] - (*this)[0][1] * (*this)[2][0]);
194  (result)[2][0] = (*this)[0][1] * (*this)[1][2] - (*this)[0][2] * (*this)[1][1];
195  (result)[2][1] = -((*this)[0][0] * (*this)[1][2] - (*this)[0][2] * (*this)[1][0]);
196  (result)[2][2] = (*this)[0][0] * (*this)[1][1] - (*this)[0][1] * (*this)[1][0];
197  break;
198 
199  default:
200  int i1, j1;
201  tlp::Matrix<Obj, SIZE - 1> c;
202 
203  for (size_t j = 0; j < SIZE; ++j) {
204  for (size_t i = 0; i < SIZE; ++i) {
205  i1 = 0;
206 
207  for (size_t ii = 0; ii < SIZE; ++ii) {
208  if (ii == i)
209  continue;
210 
211  j1 = 0;
212 
213  for (size_t jj = 0; jj < SIZE; jj++) {
214  if (jj == j)
215  continue;
216 
217  c[i1][j1] = (*this)[ii][jj];
218  ++j1;
219  }
220 
221  ++i1;
222  }
223 
224  if ((i + j) & 1)
225  result[i][j] = c.determinant();
226  else
227  result[i][j] = -c.determinant();
228  }
229  }
230 
231  break;
232  }
233 
234  return result;
235 }
236 //=====================================================================================
237 template <typename Obj, size_t SIZE>
238 MATRIX &MATRIX::transpose() {
239  Obj tmp;
240 
241  for (size_t i = 1; i < SIZE; ++i) {
242  for (size_t j = 0; j < i; ++j) {
243  tmp = (*this)[i][j];
244  (*this)[i][j] = (*this)[j][i];
245  (*this)[j][i] = tmp;
246  }
247  }
248 
249  return (*this);
250 }
251 //=====================================================================================
252 template <typename Obj, size_t SIZE>
253 MATRIX &MATRIX::inverse() {
254  (*this) = (*this).cofactor().transpose() /= (*this).determinant();
255  return (*this);
256 }
257 //=====================================================================================
258 template <typename Obj, size_t SIZE>
259 MATRIX tlp::operator+(const MATRIX &mat1, const MATRIX &mat2) {
260  return MATRIX(mat1) += mat2;
261 }
262 //=====================================================================================
263 template <typename Obj, size_t SIZE>
264 MATRIX tlp::operator-(const MATRIX &mat1, const MATRIX &mat2) {
265  return MATRIX(mat1) -= mat2;
266 }
267 //=====================================================================================
268 template <typename Obj, size_t SIZE>
269 MATRIX tlp::operator*(const MATRIX &mat1, const MATRIX &mat2) {
270  MATRIX result;
271 
272  for (size_t i = 0; i < SIZE; ++i)
273  for (size_t j = 0; j < SIZE; ++j) {
274  result[i][j] = mat1[i][0] * mat2[0][j];
275 
276  for (size_t k = 1; k < SIZE; ++k)
277  result[i][j] += mat1[i][k] * mat2[k][j];
278  }
279 
280  return result;
281 }
282 //=====================================================================================
283 template <typename Obj, size_t SIZE>
284 MATRIX MATRIX::operator/(const MATRIX &mat2) const {
285  return MATRIX(*this) /= mat2;
286 }
287 //=====================================================================================
288 template <typename Obj, size_t SIZE>
289 MATRIX MATRIX::operator/(const Obj &obj) const {
290  return MATRIX(*this) /= obj;
291 }
292 //=====================================================================================
293 template <typename Obj, size_t SIZE>
294 MATRIX tlp::operator*(const MATRIX &mat, const Obj &obj) {
295  return MATRIX(mat) *= obj;
296 }
297 //=====================================================================================
298 template <typename Obj, size_t SIZE>
299 tlp::Vector<Obj, SIZE> tlp::operator*(const MATRIX &mat, const tlp::Vector<Obj, SIZE> &vec) {
300  tlp::Vector<Obj, SIZE> result;
301 
302  for (size_t row = 0; row < SIZE; ++row) {
303  result[row] = mat[row][0] * vec[0];
304  }
305 
306  for (size_t col = 1; col < SIZE; ++col) {
307  for (size_t row = 0; row < SIZE; ++row) {
308  result[row] += mat[row][col] * vec[col];
309  }
310  }
311 
312  return result;
313 }
314 //=====================================================================================
315 template <typename Obj, size_t SIZE>
316 tlp::Vector<Obj, SIZE> tlp::operator*(const tlp::Vector<Obj, SIZE> &vec, const MATRIX &mat) {
317  tlp::Vector<Obj, SIZE> result;
318 
319  for (size_t row = 0; row < SIZE; ++row) {
320  result[row] = mat[0][row] * vec[0];
321  }
322 
323  for (size_t col = 1; col < SIZE; ++col) {
324  for (size_t row = 0; row < SIZE; ++row) {
325  result[row] += mat[col][row] * vec[col];
326  }
327  }
328 
329  return result;
330 }
331 
332 //=====================================================================================
333 template <typename Obj, size_t SIZE>
334 tlp::Vector<Obj, SIZE> MATRIX::powerIteration(const unsigned int nIterations) const {
335  tlp::Vector<Obj, SIZE> iteration;
336 
337  for (size_t i = 0; i < SIZE; i++)
338  iteration[i] = 1;
339 
340  for (unsigned int i = 0; i < nIterations; i++) {
341  iteration = (*this) * iteration;
342 
343  iteration /= iteration.norm();
344  }
345 
346  return iteration;
347 }
348 
349 //=====================================================================================
350 
351 template <typename Obj, size_t SIZE>
352 bool MATRIX::simplify(tlp::Matrix<Obj, 2> &simplifiedMatrix) const {
353  if (SIZE != 3) {
354  tlp::warning() << "Computation allowed only for 3x3 Matrices. Yours sizes : " << SIZE << "x"
355  << SIZE << std::endl;
356 
357  return false;
358  }
359 
360  // We start with a matrix representing an equation system under the following form :
361  //
362  // ax + by + cz = 0
363  // dx + ey + fz = 0
364  // gx + hy + iz = 0
365  //
366  // So M looks like :
367  //
368  // ( ax by cz ) *(e1)*
369  // M = ( dx ey fz ) *(e2)*
370  // ( gx hy iz ) *(e3)*
371  //
372  // What we want is something like that :
373  //
374  // jx + ky = 0
375  // lx + mz = 0
376  //
377  // So to reduce the matrix, we will use the Gaussian Elimination.
378  // For the first line we will apply a Gaussian Elimination between (e1) and (e2)
379  // For the second line we will apply a Gaussian Elimination between (e1) and (e3)
380 
381  float coeff;
382 
383  // First Gaussian Elimination :
384  // The pivot is z
385 
386  coeff = (*this)[1][2] / (*this)[0][2]; // fz / cz
387 
388  // After that:
389  // jx = dx - (coeff * ax)
390  // ky = ey - (coeff * by)
391  simplifiedMatrix[0][0] = (*this)[1][0] - (coeff * (*this)[0][0]);
392  simplifiedMatrix[0][1] = (*this)[1][1] - (coeff * (*this)[0][1]);
393 
394  // Second Gaussian Elimination :
395  // The pivot is y
396 
397  coeff = (*this)[2][1] / (*this)[0][1]; // hy / by
398 
399  // Idem :
400  // lx = gx - (coeff * ax)
401  // mz = iz - (coeff * cz)
402  simplifiedMatrix[1][0] = (*this)[2][0] - (coeff * (*this)[0][0]);
403  simplifiedMatrix[1][1] = (*this)[2][2] - (coeff * (*this)[0][2]);
404 
405  return true;
406 }
407 
408 //=====================================================================================
409 
410 template <typename Obj, size_t SIZE>
411 bool MATRIX::computeEigenVector(const float x, tlp::Vector<Obj, 3> &eigenVector) const {
412  if (SIZE != 2) {
413  tlp::warning() << "Computation allowed only for 2x2 Matrices. Yours sizes : " << SIZE << "x"
414  << SIZE << std::endl;
415 
416  return false;
417  }
418 
419  eigenVector[0] = x; // Fixed by user
420 
421  // We know that the matrix we are using is under that form :
422  //
423  // ( ax by )
424  // M = ( )
425  // ( cx dz )
426  //
427  // Since we have a fixed x, we can compute y and z :
428  //
429  // y = -a / b
430  // z = -c / d
431 
432  float a, b, c, d;
433 
434  a = (*this)[0][0];
435  b = (*this)[0][1];
436  c = (*this)[1][0];
437  d = (*this)[1][1];
438 
439  eigenVector[1] = (-a * x) / b;
440  eigenVector[2] = (-c * x) / d;
441 
442  return true;
443 }