Tulip  4.2.0
Better Visualization Through Research
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Properties Groups Pages
Matrix.cxx
1 /*
2  *
3  * This file is part of Tulip (www.tulip-software.org)
4  *
5  * Authors: David Auber and the Tulip development Team
6  * from LaBRI, University of Bordeaux 1 and Inria Bordeaux - Sud Ouest
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 }
48 
49 
50 //===================================================================
51 template<typename Obj,unsigned int SIZE>
52 MATRIX::Matrix(const std::vector<std::vector<Obj> > &covarianceMatrix) {
53  for(unsigned int i=0; i < SIZE; i++)
54  for(unsigned int j=0; j < SIZE; j++)
55  (*this)[i][j] = covarianceMatrix[i][j] / (sqrt(covarianceMatrix[i][i] * covarianceMatrix[j][j]));
56 }
57 
58 //===================================================================
59 template<typename Obj,unsigned int SIZE>
60 MATRIX & MATRIX::fill(Obj obj) {
61  for (unsigned int i=0; i<SIZE; ++i)
62  (*this)[i].fill(obj);
63 
64  return (*this);
65 }
66 //======================================================
67 template<typename Obj,unsigned int SIZE>
68 MATRIX& MATRIX::operator+=(const MATRIX &mat) {
69  for (unsigned int i=0; i<SIZE; ++i)
70  (*this)[i] += mat[i];
71 
72  return (*this);
73 }
74 //======================================================
75 template<typename Obj,unsigned int SIZE>
76 MATRIX& MATRIX::operator-=(const MATRIX &mat) {
77  for (unsigned int i=0; i<SIZE; ++i)
78  (*this)[i] -= mat[i];
79 
80  return (*this);
81 }
82 //======================================================
83 template<typename Obj,unsigned int SIZE>
84 bool MATRIX::operator==(const MATRIX &mat) const {
85  for (unsigned int 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,unsigned int SIZE>
94 bool MATRIX::operator!=(const MATRIX &mat) const {
95  for (unsigned int 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,unsigned int SIZE>
104 MATRIX & MATRIX::operator*=(const MATRIX &mat) {
105  (*this) = (*this) * mat;
106  return (*this);
107 }
108 //=====================================================================================
109 template<typename Obj,unsigned int SIZE>
110 MATRIX & MATRIX::operator*=(const Obj &obj) {
111  for (unsigned int i=0; i<SIZE; ++i)
112  (*this)[i] *= obj;
113 
114  return (*this);
115 }
116 //=====================================================================================
117 template<typename Obj,unsigned int 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,unsigned int SIZE>
126 MATRIX & MATRIX::operator/=(const Obj &obj) {
127  for (unsigned int i=0; i<SIZE; ++i)
128  (*this)[i] /= obj;
129 
130  return (*this);
131 }
132 //=====================================================================================
133 template<typename Obj,unsigned int 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 (unsigned int j1=0; j1<SIZE; ++j1) {
151  tlp::Matrix<Obj, SIZE - 1> m;
152 
153  for (unsigned int i=1; i<SIZE; i++) {
154  j2 = 0;
155 
156  for (unsigned int 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,unsigned int 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 (unsigned int j=0; j<SIZE; ++j) {
204  for (unsigned int i=0; i<SIZE; ++i) {
205  i1 = 0;
206 
207  for (unsigned int ii=0; ii<SIZE; ++ii) {
208  if (ii == i)
209  continue;
210 
211  j1 = 0;
212 
213  for (unsigned int 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) result[i][j]=c.determinant();
225  else result[i][j]=-c.determinant();
226  }
227  }
228 
229  break;
230  }
231 
232  return result;
233 }
234 //=====================================================================================
235 template<typename Obj,unsigned int SIZE>
236 MATRIX & MATRIX::transpose() {
237  register Obj tmp;
238 
239  for (unsigned int i=1; i<SIZE; ++i) {
240  for (unsigned int j=0; j<i; ++j) {
241  tmp = (*this)[i][j];
242  (*this)[i][j] = (*this)[j][i];
243  (*this)[j][i] = tmp;
244  }
245  }
246 
247  return (*this);
248 }
249 //=====================================================================================
250 template<typename Obj,unsigned int SIZE>
251 MATRIX & MATRIX::inverse() {
252  (*this) = (*this).cofactor().transpose() /= (*this).determinant();
253  return (*this);
254 }
255 //=====================================================================================
256 template<typename Obj,unsigned int SIZE>
257 MATRIX tlp::operator+(const MATRIX &mat1 ,const MATRIX &mat2) {
258  return MATRIX(mat1)+=mat2;
259 }
260 //=====================================================================================
261 template<typename Obj,unsigned int SIZE>
262 MATRIX tlp::operator-(const MATRIX &mat1 ,const MATRIX &mat2) {
263  return MATRIX(mat1)-=mat2;
264 }
265 //=====================================================================================
266 template<typename Obj,unsigned int SIZE>
267 MATRIX tlp::operator*(const MATRIX &mat1 ,const MATRIX &mat2) {
268  MATRIX result;
269 
270  for (unsigned int i=0; i<SIZE; ++i)
271  for (unsigned int j=0; j<SIZE; ++j) {
272  result[i][j] = mat1[i][0] * mat2[0][j];
273 
274  for (unsigned int k=1; k<SIZE; ++k)
275  result[i][j] += mat1[i][k] * mat2[k][j];
276  }
277 
278  return result;
279 }
280 //=====================================================================================
281 template<typename Obj,unsigned int SIZE>
282 MATRIX MATRIX::operator/(const MATRIX &mat2) const {
283  return MATRIX(*this)/=mat2;
284 }
285 //=====================================================================================
286 template<typename Obj,unsigned int SIZE>
287 MATRIX MATRIX::operator/(const Obj &obj) const {
288  return MATRIX(*this) /= obj;
289 }
290 //=====================================================================================
291 template<typename Obj,unsigned int SIZE>
292 MATRIX tlp::operator*(const MATRIX &mat ,const Obj &obj) {
293  return MATRIX(mat) *= obj;
294 }
295 //=====================================================================================
296 template<typename Obj, unsigned int SIZE>
297 tlp::Vector<Obj, SIZE> tlp::operator*( const MATRIX &mat , const tlp::Vector<Obj, SIZE> &vec) {
298  tlp::Vector<Obj,SIZE> result;
299 
300  for (unsigned int row=0; row<SIZE; ++row) {
301  result[row] = mat[row][0] * vec[0];
302  }
303 
304  for (unsigned int col=1; col<SIZE; ++col) {
305  for (unsigned int row=0; row<SIZE; ++row) {
306  result[row] += mat[row][col] * vec[col];
307  }
308  }
309 
310  return result;
311 }
312 //=====================================================================================
313 template<typename Obj,unsigned int SIZE>
314 tlp::Vector<Obj,SIZE> tlp::operator*( const tlp::Vector<Obj,SIZE> &vec, const MATRIX &mat) {
315  tlp::Vector<Obj,SIZE> result;
316 
317  for (unsigned int row=0; row<SIZE; ++row) {
318  result[row] = mat[0][row] * vec[0];
319  }
320 
321  for (unsigned int col=1; col<SIZE; ++col) {
322  for (unsigned int row=0; row<SIZE; ++row) {
323  result[row] += mat[col][row] * vec[col];
324  }
325  }
326 
327  return result;
328 }
329 
330 //=====================================================================================
331 template<typename Obj, unsigned int SIZE>
332 tlp::Vector<Obj, SIZE> MATRIX::powerIteration(const unsigned int nIterations) const {
333  tlp::Vector<Obj, SIZE> iteration;
334 
335  for(unsigned int i=0; i < SIZE; i++)
336  iteration[i] = 1;
337 
338  for(unsigned int i=0; i < nIterations; i++) {
339  iteration = (*this) * iteration;
340 
341  iteration /= iteration.norm();
342  }
343 
344  return iteration;
345 }
346 
347 //=====================================================================================
348 
349 template<typename Obj, unsigned int SIZE>
350 bool MATRIX::simplify(tlp::Matrix<Obj, 2> &simplifiedMatrix) const {
351  if (SIZE != 3) {
352  tlp::warning() << "Computation allowed only for 3x3 Matrices. Yours sizes : " << SIZE << "x" << SIZE << std::endl;
353 
354  return false;
355  }
356 
357  // We start with a matrix representing an equation system under the following form :
358  //
359  // ax + by + cz = 0
360  // dx + ey + fz = 0
361  // gx + hy + iz = 0
362  //
363  // So M looks like :
364  //
365  // ( ax by cz ) *(e1)*
366  // M = ( dx ey fz ) *(e2)*
367  // ( gx hy iz ) *(e3)*
368  //
369  // What we want is something like that :
370  //
371  // jx + ky = 0
372  // lx + mz = 0
373  //
374  // So to reduce the matrix, we will use the Gaussian Elimination.
375  // For the first line we will apply a Gaussian Elimination between (e1) and (e2)
376  // For the second line we will apply a Gaussian Elimination between (e1) and (e3)
377 
378  float coeff;
379 
380  // First Gaussian Elimination :
381  // The pivot is z
382 
383  coeff = (*this)[1][2] / (*this)[0][2]; // fz / cz
384 
385  // After that:
386  // jx = dx - (coeff * ax)
387  // ky = ey - (coeff * by)
388  simplifiedMatrix[0][0] = (*this)[1][0] - (coeff * (*this)[0][0]);
389  simplifiedMatrix[0][1] = (*this)[1][1] - (coeff * (*this)[0][1]);
390 
391  // Second Gaussian Elimination :
392  // The pivot is y
393 
394  coeff = (*this)[2][1] / (*this)[0][1]; // hy / by
395 
396  // Idem :
397  // lx = gx - (coeff * ax)
398  // mz = iz - (coeff * cz)
399  simplifiedMatrix[1][0] = (*this)[2][0] - (coeff * (*this)[0][0]);
400  simplifiedMatrix[1][1] = (*this)[2][2] - (coeff * (*this)[0][2]);
401 
402  return true;
403 }
404 
405 //=====================================================================================
406 
407 template<typename Obj, unsigned int SIZE>
408 bool MATRIX::computeEigenVector(const float x, tlp::Vector<Obj, 3> &eigenVector) const {
409  if (SIZE != 2) {
410  tlp::warning() << "Computation allowed only for 2x2 Matrices. Yours sizes : " << SIZE << "x" << SIZE << std::endl;
411 
412  return false;
413  }
414 
415  eigenVector[0] = x; // Fixed by user
416 
417  // We know that the matrix we are using is under that form :
418  //
419  // ( ax by )
420  // M = ( )
421  // ( cx dz )
422  //
423  // Since we have a fixed x, we can compute y and z :
424  //
425  // y = -a / b
426  // z = -c / d
427 
428  float a, b, c, d;
429 
430  a = (*this)[0][0];
431  b = (*this)[0][1];
432  c = (*this)[1][0];
433  d = (*this)[1][1];
434 
435  eigenVector[1] = (-a * x) / b;
436  eigenVector[2] = (-c * x) / d;
437 
438  return true;
439 }