Engauge Digitizer  2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Matrix.cpp
Go to the documentation of this file.
1 /******************************************************************************************************
2  * (C) 2016 markummitchell@github.com. This file is part of Engauge Digitizer, which is released *
3  * under GNU General Public License version 2 (GPLv2) or (at your option) any later version. See file *
4  * LICENSE or go to gnu.org/licenses for details. Distribution requires prior written permission. *
5  ******************************************************************************************************/
6 
7 #include "EngaugeAssert.h"
8 #include "Matrix.h"
9 #include <qmath.h>
10 #include <QTextStream>
11 
13 {
14  initialize (N, N);
15 }
16 
17 Matrix::Matrix (int rows,
18  int cols)
19 {
20  initialize (rows, cols);
21 }
22 
23 Matrix::Matrix (const Matrix &other)
24 {
25  m_rows = other.rows();
26  m_cols = other.cols();
27  m_vector.resize (m_rows * m_cols);
28  for (int row = 0; row < m_rows; row++) {
29  for (int col = 0; col < m_cols; col++) {
30  set (row, col, other.get (row, col));
31  }
32  }
33 }
34 
36 {
37  m_rows = other.rows();
38  m_cols = other.cols();
39  m_vector.resize (m_rows * m_cols);
40  for (int row = 0; row < m_rows; row++) {
41  for (int col = 0; col < m_cols; col++) {
42  set (row, col, other.get (row, col));
43  }
44  }
45 
46  return *this;
47 }
48 
49 void Matrix::addRowToAnotherWithScaling (int rowFrom,
50  int rowTo,
51  double factor)
52 {
53  for (int col = 0; col < cols (); col++) {
54  double oldValueFrom = get (rowFrom, col);
55  double oldValueTo = get (rowTo, col);
56  double newValueTo = oldValueFrom * factor + oldValueTo;
57  set (rowTo, col, newValueTo);
58  }
59 }
60 
61 int Matrix::cols () const
62 {
63  return m_cols;
64 }
65 
66 double Matrix::determinant () const
67 {
68  ENGAUGE_ASSERT (m_rows == m_cols);
69 
70  double rtn;
71 
72  if (m_rows == 1) {
73 
74  rtn = m_vector [0];
75 
76  } else {
77 
78  const int COL = 0; // We arbitrarily iterate through the first column
79 
80  // This is a recursive algorithm
81  rtn = 0.0;
82  double multiplier = +1;
83  for (int row = 0; row < m_rows; row++) {
84  Matrix min = minorReduced (row, COL);
85  rtn += multiplier * get (row, COL) * min.determinant ();
86  multiplier *= -1.0;
87  }
88  }
89 
90  return rtn;
91 }
92 
93 int Matrix::fold2dIndexes (int row, int col) const
94 {
95  return row * m_cols + col;
96 }
97 
98 double Matrix::get (int row, int col) const
99 {
100  int foldedIndex = fold2dIndexes (row, col);
101  return m_vector [foldedIndex];
102 }
103 
104 void Matrix::initialize (int rows,
105  int cols)
106 {
107  m_rows = rows;
108  m_cols = cols;
109  m_vector.resize (rows * cols);
110  for (int row = 0; row < m_rows; row++) {
111  for (int col = 0; col < m_cols; col++) {
112 
113  // Identity matrix for square matrix, otherwise zero matrix
114  if (row == col && m_rows == m_cols) {
115  set (row, col, 1.0);
116  } else {
117  set (row, col, 0.0);
118  }
119  }
120  }
121 }
122 
123 Matrix Matrix::inverse (int significantDigits,
124  MatrixConsistent &matrixConsistent) const
125 {
126  // Set epsilon threshold for valueFailsEpsilonTest
127  double maxValue = 0;
128  for (int row = 0; row < m_rows; row++) {
129  for (int col = 0; col < m_cols; col++) {
130  double value = qAbs (get (row, col));
131  if (value > maxValue) {
132  maxValue = value;
133  }
134  }
135  }
136 
137  // Available algorithms are inverseCramersRule and inverseGaussianElimination
138  matrixConsistent = MATRIX_CONSISTENT;
139  return inverseGaussianElimination (significantDigits,
140  matrixConsistent);
141 }
142 
143 Matrix Matrix::inverseCramersRule (MatrixConsistent & matrixConsistent,
144  double epsilonThreshold) const
145 {
146  ENGAUGE_ASSERT (m_rows == m_cols);
147 
148  Matrix inv (m_rows);
149  int row, col;
150 
151  if (m_rows > 1) {
152 
153  // Compute cofactor matrix
154  double multiplierStartForRow = 1.0;
155  Matrix cofactor (m_rows);
156  for (row = 0; row < m_rows; row++) {
157  double multiplier = multiplierStartForRow;
158  for (col = 0; col < m_cols; col++) {
159  Matrix min = minorReduced (row, col);
160  double element = multiplier * min.determinant ();
161  multiplier *= -1.0;
162  cofactor.set (row, col, element);
163  }
164  multiplierStartForRow *= -1.0;
165  }
166 
167  // Compute adjoint
168  Matrix adjoint = cofactor.transpose ();
169 
170  double determ = determinant ();
171  if (valueFailsEpsilonTest (determ,
172  epsilonThreshold)) {
173  matrixConsistent = MATRIX_INCONSISTENT;
174  return inv;
175  }
176 
177  // Inverse is the adjoint divided by the determinant
178  for (row = 0; row < m_rows; row++) {
179  for (col = 0; col < m_cols; col++) {
180  inv.set (row, col, adjoint.get (row, col) / determ);
181  }
182  }
183  } else {
184  double denominator = get (0, 0);
185  if (valueFailsEpsilonTest (denominator,
186  epsilonThreshold)) {
187  matrixConsistent = MATRIX_INCONSISTENT;
188  return inv;
189  }
190  inv.set (0, 0, 1.0 / denominator);
191  }
192 
193  return inv;
194 }
195 
196 Matrix Matrix::inverseGaussianElimination (int significantDigits,
197  MatrixConsistent &matrixConsistent) const
198 {
199  // From https://en.wikipedia.org/wiki/Gaussian_elimination
200 
201  int row, col, rowFrom, rowTo;
202  Matrix inv (rows ());
203 
204  // Track the row switches that may or may not be performed below
205  QVector<int> rowIndexes (rows ());
206  for (row = 0; row < rows (); row++) {
207  rowIndexes [row] = row;
208  }
209 
210  // Create the working matrix and populate the left half. Number of columns in the working matrix is twice the number
211  // of cols in this matrix, but we will not populate the right half until after the bubble sort below
212  Matrix working (rows (), 2 * cols ());
213  for (row = 0; row < rows (); row++) {
214  for (col = 0; col < cols (); col++) {
215  double value = get (row, col);
216  working.set (row, col, value);
217  }
218  }
219 
220  // Simple bubble sort to rearrange the rows with 0 leading zeros at start, followed by rows with 1 leading zero at start, ...
221  for (int row1 = 0; row1 < rows (); row1++) {
222  for (int row2 = row1 + 1; row2 < rows (); row2++) {
223  if (working.leadingZeros (row1) > working.leadingZeros (row2)) {
224  working.switchRows (row1, row2);
225 
226  // Remember this switch
227  int row1Index = rowIndexes [row1];
228  int row2Index = rowIndexes [row2];
229  rowIndexes [row1] = row2Index;
230  rowIndexes [row2] = row1Index;
231  }
232  }
233  }
234 
235  // Populate the right half now
236  for (row = 0; row < rows (); row++) {
237  for (col = cols (); col < 2 * cols (); col++) {
238  int colIdentitySubmatrix = col - cols ();
239  double value = (row == colIdentitySubmatrix) ? 1 : 0;
240  working.set (row, col, value);
241  }
242  }
243 
244  // Loop through the "from" row going down. This results in the lower off-diagonal terms becoming zero, in the left half
245  for (rowFrom = 0; rowFrom < rows (); rowFrom++) {
246  int colFirstWithNonZero = rowFrom;
247 
248  // In pathological situations we have (rowFrom, colFirstWithNonzero) = 0 in which case the solution cannot be obtained
249  // so we exit
250  if (qAbs (working.get (rowFrom, colFirstWithNonZero)) <= 0) {
251  matrixConsistent = MATRIX_INCONSISTENT;
252  return inv;
253  }
254 
255  // Normalize the 'from' row with first nonzero term set to 1
256  working.normalizeRow (rowFrom, colFirstWithNonZero, significantDigits, matrixConsistent);
257  if (matrixConsistent == MATRIX_INCONSISTENT) {
258  return inv;
259  }
260 
261  // Apply the 'from' row to all the 'to' rows
262  for (rowTo = rowFrom + 1; rowTo < rows (); rowTo++) {
263 
264  if (qAbs (working.get (rowTo, colFirstWithNonZero)) > 0) {
265 
266  // We need to merge rowFrom and rowTo into rowTo
267  double denominator = working.get (rowFrom, colFirstWithNonZero);
268  double factor = -1.0 * working.get (rowTo, colFirstWithNonZero) / denominator;
269  working.addRowToAnotherWithScaling (rowFrom, rowTo, factor);
270  }
271  }
272  }
273 
274  // Loop through the "from" row doing up. This results in the upper off-diagonal terms becoming zero, in the left half
275  for (rowFrom = rows () - 1; rowFrom >= 0; rowFrom--) {
276  int colFirstWithNonZero = rowFrom; // This is true since we should have 1s all down the diagonal at this point
277 
278  // Normalize the 'from' row with diagonal term set to 1. The first term should be like 0.9999 or 1.0001 but we want exactly one
279  MatrixConsistent matrixConsistentIter;
280  working.normalizeRow (rowFrom, colFirstWithNonZero, significantDigits, matrixConsistentIter);
281  if (matrixConsistentIter == MATRIX_INCONSISTENT) {
282  return inv;
283  }
284 
285  // Apply the 'from' row to all the 'to' rows
286  for (rowTo = rowFrom - 1; rowTo >= 0; rowTo--) {
287 
288  if (qAbs (working.get (rowTo, colFirstWithNonZero)) > 0) {
289 
290  // We need to merge rowFro and rowTo into rowTo
291  double denominator = working.get (rowFrom, colFirstWithNonZero);
292  double factor = -1.0 * working.get (rowTo, colFirstWithNonZero) / denominator;
293  working.addRowToAnotherWithScaling (rowFrom, rowTo, factor);
294  }
295  }
296  }
297 
298  // Extract right half of rectangular matrix which is the inverse
299 
300  for (row = 0; row < working.rows (); row++) {
301 
302  int rowBeforeSort = rowIndexes [row];
303 
304  for (col = 0; col < working.rows (); col++) {
305  int colRightHalf = col + inv.cols ();
306  double value = working.get (rowBeforeSort, colRightHalf);
307  inv.set (row, col, value);
308  }
309  }
310 
311  return inv;
312 }
313 
314 unsigned int Matrix::leadingZeros (int row) const
315 {
316  unsigned int sum = 0;
317  for (int col = 0; col < cols (); col++) {
318  if (qAbs (get (row, col)) > 0) {
319  ++sum;
320  }
321  }
322 
323  return sum;
324 }
325 
326 Matrix Matrix::minorReduced (int rowOmit, int colOmit) const
327 {
328  ENGAUGE_ASSERT (m_rows == m_cols);
329 
330  Matrix outMinor (m_rows - 1);
331  int rowMinor = 0;
332  for (int row = 0; row < m_rows; row++) {
333 
334  if (row != rowOmit) {
335 
336  int colMinor = 0;
337  for (int col = 0; col < m_cols; col++) {
338 
339  if (col != colOmit) {
340 
341  outMinor.set (rowMinor, colMinor, get (row, col));
342  ++colMinor;
343  }
344  }
345  ++rowMinor;
346  }
347  }
348 
349  return outMinor;
350 }
351 
352 void Matrix::normalizeRow (int rowToNormalize,
353  int colToNormalize,
354  int significantDigits,
355  MatrixConsistent &matrixConsistent)
356 {
357  double denominator = get (rowToNormalize, colToNormalize);
358 
359  // Epsilon is computed from smallest value in row
360  double smallestAbsValue = 0;
361  for (int col = 0; col < cols (); col++) {
362  double absValue = qAbs (get (rowToNormalize, 0));
363  if (col == 0 || absValue < smallestAbsValue) {
364  smallestAbsValue = absValue;
365  }
366  }
367  double epsilonThreshold = smallestAbsValue / qPow (10.0, significantDigits);
368 
369  if (valueFailsEpsilonTest (denominator,
370  epsilonThreshold)) {
371 
372  matrixConsistent = MATRIX_INCONSISTENT;
373 
374  } else {
375 
376  matrixConsistent = MATRIX_CONSISTENT;
377 
378  double factor = 1.0 / denominator;
379  for (int col = 0; col < cols (); col++) {
380  double value = get (rowToNormalize, col);
381  set (rowToNormalize, col, factor * value);
382  }
383  }
384 }
385 
386 Matrix Matrix::operator* (const Matrix &other) const
387 {
388  ENGAUGE_ASSERT (m_cols == other.rows ());
389 
390  Matrix out (m_rows, other.cols ());
391 
392  for (int row = 0; row < m_rows; row++) {
393  for (int col = 0; col < other.cols (); col++) {
394  double sum = 0;
395  for (int index = 0; index < m_cols; index++) {
396  sum += get (row, index) * other.get (index, col);
397  }
398  out.set (row, col, sum);
399  }
400  }
401 
402  return out;
403 }
404 
405 QVector<double> Matrix::operator* (const QVector<double> other) const
406 {
407  ENGAUGE_ASSERT (m_cols == other.size ());
408 
409  QVector<double> out;
410  out.resize (m_rows);
411  for (int row = 0; row < m_rows; row++) {
412  double sum = 0;
413  for (int col = 0; col < m_cols; col++) {
414  sum += get (row, col) * other [col];
415  }
416 
417  out [row] = sum;
418  }
419 
420  return out;
421 }
422 
423 int Matrix::rows () const
424 {
425  return m_rows;
426 }
427 
428 void Matrix::set (int row, int col, double value)
429 {
430  m_vector [fold2dIndexes (row, col)] = value;
431 }
432 
433 void Matrix::switchRows (int row1,
434  int row2)
435 {
436  // Switch rows
437  for (int col = 0; col < cols (); col++) {
438  double temp1 = get (row1, col);
439  double temp2 = get (row2, col);
440  set (row2, col, temp2);
441  set (row1, col, temp1);
442  }
443 }
444 
445 QString Matrix::toString () const
446 {
447  QString out;
448  QTextStream str (&out);
449 
450  str << "(";
451  for (int row = 0; row < rows (); row++) {
452  if (row > 0) {
453  str << ", ";
454  }
455  str << "(";
456  for (int col = 0; col < cols (); col++) {
457  if (col > 0) {
458  str << ", ";
459  }
460  str << get (row, col);
461  }
462  str << ")";
463  }
464  str << ")";
465 
466  return out;
467 }
468 
470 {
471  Matrix out (m_cols, m_rows);
472 
473  for (int row = 0; row < m_rows; row++) {
474  for (int col = 0; col < m_cols; col++) {
475  out.set (col, row, get (row, col));
476  }
477  }
478 
479  return out;
480 }
481 
482 bool Matrix::valueFailsEpsilonTest (double value,
483  double epsilonThreshold) const
484 {
485  return qAbs (value) < qAbs (epsilonThreshold);
486 }
int rows() const
Height of matrix.
Definition: Matrix.cpp:423
Matrix minorReduced(int rowOmit, int colOmit) const
Return minor matrix which is the original with the specified row and column omitted. The name &#39;minor&#39; is a reserved word.
Definition: Matrix.cpp:326
Matrix inverse(int significantDigits, MatrixConsistent &matrixConsistent) const
Return the inverse of this matrix.
Definition: Matrix.cpp:123
Matrix operator*(const Matrix &other) const
Multiplication operator with a matrix.
Definition: Matrix.cpp:386
double determinant() const
Return the determinant of this matrix.
Definition: Matrix.cpp:66
MatrixConsistent
Indicates if matrix is consistent (i.e. has at least one solution)
Definition: Matrix.h:14
Matrix & operator=(const Matrix &matrix)
Assignment operator.
Definition: Matrix.cpp:35
Matrix(int N)
Simple constructor of square matrix with initialization to identity matrix.
Definition: Matrix.cpp:12
Matrix transpose() const
Return the transpose of the current matrix.
Definition: Matrix.cpp:469
Matrix class that supports arbitrary NxN size.
Definition: Matrix.h:20
QString toString() const
Dump matrix to a string.
Definition: Matrix.cpp:445
int cols() const
Width of matrix.
Definition: Matrix.cpp:61
void set(int row, int col, double value)
Set (row, col) element.
Definition: Matrix.cpp:428
#define ENGAUGE_ASSERT(cond)
Drop in replacement for Q_ASSERT if defined(QT_NO_DEBUG) &amp;&amp; !defined(QT_FORCE_ASSERTS) define ENGAUGE...
Definition: EngaugeAssert.h:20
double get(int row, int col) const
Return (row, col) element.
Definition: Matrix.cpp:98