Tulip 5.7.1
Large graphs analysis and drawing
Loading...
Searching...
No Matches
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
23namespace tlp {
24template <typename Obj>
25class Matrix<Obj, 1> : public Array<Vector<Obj, 1>, 1> {
26public:
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//===================================================================
50template <typename Obj, size_t SIZE>
51MATRIX::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//===================================================================
59template <typename Obj, size_t SIZE>
60MATRIX &MATRIX::fill(Obj obj) {
61 for (size_t i = 0; i < SIZE; ++i)
62 (*this)[i].fill(obj);
63
64 return (*this);
65}
66//======================================================
67template <typename Obj, size_t SIZE>
68MATRIX &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//======================================================
75template <typename Obj, size_t SIZE>
76MATRIX &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//======================================================
83template <typename Obj, size_t SIZE>
84bool 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//======================================================
93template <typename Obj, size_t SIZE>
94bool 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//===================================================================
103template <typename Obj, size_t SIZE>
104MATRIX &MATRIX::operator*=(const MATRIX &mat) {
105 (*this) = (*this) * mat;
106 return (*this);
107}
108//=====================================================================================
109template <typename Obj, size_t SIZE>
110MATRIX &MATRIX::operator*=(const Obj &obj) {
111 for (size_t i = 0; i < SIZE; ++i)
112 (*this)[i] *= obj;
113
114 return (*this);
115}
116//=====================================================================================
117template <typename Obj, size_t SIZE>
118MATRIX &MATRIX::operator/=(const MATRIX &mat) {
119 MATRIX tmpMat(mat);
120 tmpMat.inverse();
121 (*this) *= tmpMat;
122 return (*this);
123}
124//=====================================================================================
125template <typename Obj, size_t SIZE>
126MATRIX &MATRIX::operator/=(const Obj &obj) {
127 for (size_t i = 0; i < SIZE; ++i)
128 (*this)[i] /= obj;
129
130 return (*this);
131}
132//=====================================================================================
133template <typename Obj, size_t SIZE>
134Obj 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//=====================================================================================
175template <typename Obj, size_t SIZE>
176MATRIX 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//=====================================================================================
237template <typename Obj, size_t SIZE>
238MATRIX &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//=====================================================================================
252template <typename Obj, size_t SIZE>
253MATRIX &MATRIX::inverse() {
254 (*this) = (*this).cofactor().transpose() /= (*this).determinant();
255 return (*this);
256}
257//=====================================================================================
258template <typename Obj, size_t SIZE>
259MATRIX tlp::operator+(const MATRIX &mat1, const MATRIX &mat2) {
260 return MATRIX(mat1) += mat2;
261}
262//=====================================================================================
263template <typename Obj, size_t SIZE>
264MATRIX tlp::operator-(const MATRIX &mat1, const MATRIX &mat2) {
265 return MATRIX(mat1) -= mat2;
266}
267//=====================================================================================
268template <typename Obj, size_t SIZE>
269MATRIX 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//=====================================================================================
283template <typename Obj, size_t SIZE>
284MATRIX MATRIX::operator/(const MATRIX &mat2) const {
285 return MATRIX(*this) /= mat2;
286}
287//=====================================================================================
288template <typename Obj, size_t SIZE>
289MATRIX MATRIX::operator/(const Obj &obj) const {
290 return MATRIX(*this) /= obj;
291}
292//=====================================================================================
293template <typename Obj, size_t SIZE>
294MATRIX tlp::operator*(const MATRIX &mat, const Obj &obj) {
295 return MATRIX(mat) *= obj;
296}
297//=====================================================================================
298template <typename Obj, size_t SIZE>
299tlp::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//=====================================================================================
315template <typename Obj, size_t SIZE>
316tlp::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//=====================================================================================
333template <typename Obj, size_t SIZE>
334tlp::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
351template <typename Obj, size_t SIZE>
352bool 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
410template <typename Obj, size_t SIZE>
411bool 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}