Engauge Digitizer  2
Matrix.cpp
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  return m_vector [fold2dIndexes (row, col)];
101 }
102 
103 void Matrix::initialize (int rows,
104  int cols)
105 {
106  m_rows = rows;
107  m_cols = cols;
108  m_vector.resize (rows * cols);
109  for (int row = 0; row < m_rows; row++) {
110  for (int col = 0; col < m_cols; col++) {
111 
112  // Identity matrix for square matrix, otherwise zero matrix
113  if (row == col && m_rows == m_cols) {
114  set (row, col, 1.0);
115  } else {
116  set (row, col, 0.0);
117  }
118  }
119  }
120 }
121 
122 Matrix Matrix::inverse (int significantDigits,
123  MatrixConsistent &matrixConsistent) const
124 {
125  // Set epsilon threshold for valueFailsEpsilonTest
126  double maxValue = 0;
127  for (int row = 0; row < m_rows; row++) {
128  for (int col = 0; col < m_cols; col++) {
129  double value = qAbs (get (row, col));
130  if (value > maxValue) {
131  maxValue = value;
132  }
133  }
134  }
135 
136  double epsilonThreshold = maxValue / qPow (10.0, significantDigits);
137 
138  // Available algorithms are inverseCramersRule and inverseGaussianElimination
139  matrixConsistent = MATRIX_CONSISTENT;
140  return inverseGaussianElimination (matrixConsistent,
141  epsilonThreshold);
142 }
143 
144 Matrix Matrix::inverseCramersRule (MatrixConsistent & matrixConsistent,
145  double epsilonThreshold) const
146 {
147  ENGAUGE_ASSERT (m_rows == m_cols);
148 
149  Matrix inv (m_rows);
150  int row, col;
151 
152  if (m_rows > 1) {
153 
154  // Compute cofactor matrix
155  double multiplierStartForRow = 1.0;
156  Matrix cofactor (m_rows);
157  for (row = 0; row < m_rows; row++) {
158  double multiplier = multiplierStartForRow;
159  for (col = 0; col < m_cols; col++) {
160  Matrix min = minorReduced (row, col);
161  double element = multiplier * min.determinant ();
162  multiplier *= -1.0;
163  cofactor.set (row, col, element);
164  }
165  multiplierStartForRow *= -1.0;
166  }
167 
168  // Compute adjoint
169  Matrix adjoint = cofactor.transpose ();
170 
171  double determ = determinant ();
172  if (valueFailsEpsilonTest (determ,
173  epsilonThreshold)) {
174  matrixConsistent = MATRIX_INCONSISTENT;
175  return inv;
176  }
177 
178  // Inverse is the adjoint divided by the determinant
179  for (row = 0; row < m_rows; row++) {
180  for (col = 0; col < m_cols; col++) {
181  inv.set (row, col, adjoint.get (row, col) / determ);
182  }
183  }
184  } else {
185  double denominator = get (0, 0);
186  if (valueFailsEpsilonTest (denominator,
187  epsilonThreshold)) {
188  matrixConsistent = MATRIX_INCONSISTENT;
189  return inv;
190  }
191  inv.set (0, 0, 1.0 / denominator);
192  }
193 
194  return inv;
195 }
196 
197 Matrix Matrix::inverseGaussianElimination (MatrixConsistent &matrixConsistent,
198  double epsilonThreshold) const
199 {
200  // From https://en.wikipedia.org/wiki/Gaussian_elimination
201 
202  int row, col, rowFrom, rowTo;
203  Matrix inv (rows ());
204 
205  // Track the row switches that may or may not be performed below
206  QVector<int> rowIndexes (rows ());
207  for (row = 0; row < rows (); row++) {
208  rowIndexes [row] = row;
209  }
210 
211  // Create the working matrix and populate the left half. Number of columns in the working matrix is twice the number
212  // of cols in this matrix, but we will not populate the right half until after the bubble sort below
213  Matrix working (rows (), 2 * cols ());
214  for (row = 0; row < rows (); row++) {
215  for (col = 0; col < cols (); col++) {
216  double value = get (row, col);
217  working.set (row, col, value);
218  }
219  }
220 
221  // Simple bubble sort to rearrange the rows with 0 leading zeros at start, followed by rows with 1 leading zero at start, ...
222  for (int row1 = 0; row1 < rows (); row1++) {
223  for (int row2 = row1 + 1; row2 < rows (); row2++) {
224  if (working.leadingZeros (row1) > working.leadingZeros (row2)) {
225  working.switchRows (row1, row2);
226 
227  // Remember this switch
228  int row1Index = rowIndexes [row1];
229  int row2Index = rowIndexes [row2];
230  rowIndexes [row1] = row2Index;
231  rowIndexes [row2] = row1Index;
232  }
233  }
234  }
235 
236  // Populate the right half now
237  for (row = 0; row < rows (); row++) {
238  for (col = cols (); col < 2 * cols (); col++) {
239  int colIdentitySubmatrix = col - cols ();
240  double value = (row == colIdentitySubmatrix) ? 1 : 0;
241  working.set (row, col, value);
242  }
243  }
244 
245  // Loop through the "from" row going down. This results in the lower off-diagonal terms becoming zero, in the left half
246  for (rowFrom = 0; rowFrom < rows (); rowFrom++) {
247  int colFirstWithNonZero = rowFrom;
248 
249  // In pathological situations we have (rowFrom, colFirstWithNonzero) = 0 in which case the solution cannot be obtained
250  // so we exit
251  if (working.get (rowFrom, colFirstWithNonZero) == 0) {
252  matrixConsistent = MATRIX_INCONSISTENT;
253  return inv;
254  }
255 
256  // Normalize the 'from' row with first nonzero term set to 1
257  working.normalizeRow (rowFrom, colFirstWithNonZero, matrixConsistent, epsilonThreshold);
258  if (matrixConsistent == MATRIX_INCONSISTENT) {
259  return inv;
260  }
261 
262  // Apply the 'from' row to all the 'to' rows
263  for (rowTo = rowFrom + 1; rowTo < rows (); rowTo++) {
264 
265  if (working.get (rowTo, colFirstWithNonZero) != 0) {
266 
267  // We need to merge rowFrom and rowTo into rowTo
268  double denominator = working.get (rowFrom, colFirstWithNonZero);
269  if (valueFailsEpsilonTest (denominator, epsilonThreshold)) {
270  matrixConsistent = MATRIX_INCONSISTENT;
271  return inv;
272  }
273  double factor = -1.0 * working.get (rowTo, colFirstWithNonZero) / denominator;
274  working.addRowToAnotherWithScaling (rowFrom, rowTo, factor);
275  }
276  }
277  }
278 
279  // Loop through the "from" row doing up. This results in the upper off-diagonal terms becoming zero, in the left half
280  for (rowFrom = rows () - 1; rowFrom >= 0; rowFrom--) {
281  int colFirstWithNonZero = rowFrom; // This is true since we should have 1s all down the diagonal at this point
282 
283  // 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
284  MatrixConsistent matrixConsistent;
285  working.normalizeRow (rowFrom, colFirstWithNonZero, matrixConsistent, epsilonThreshold);
286  if (matrixConsistent == MATRIX_INCONSISTENT) {
287  return inv;
288  }
289 
290  // Apply the 'from' row to all the 'to' rows
291  for (rowTo = rowFrom - 1; rowTo >= 0; rowTo--) {
292 
293  if (working.get (rowTo, colFirstWithNonZero) != 0) {
294 
295  // We need to merge rowFro and rowTo into rowTo
296  double denominator = working.get (rowFrom, colFirstWithNonZero);
297  if (valueFailsEpsilonTest (denominator, epsilonThreshold)) {
298  matrixConsistent = MATRIX_INCONSISTENT;
299  return inv;
300  }
301  double factor = -1.0 * working.get (rowTo, colFirstWithNonZero) / denominator;
302  working.addRowToAnotherWithScaling (rowFrom, rowTo, factor);
303  }
304  }
305  }
306 
307  // Extract right half of rectangular matrix which is the inverse
308 
309  for (row = 0; row < working.rows (); row++) {
310 
311  int rowBeforeSort = rowIndexes [row];
312 
313  for (col = 0; col < working.rows (); col++) {
314  int colRightHalf = col + inv.cols ();
315  double value = working.get (rowBeforeSort, colRightHalf);
316  inv.set (row, col, value);
317  }
318  }
319 
320  return inv;
321 }
322 
323 unsigned int Matrix::leadingZeros (int row) const
324 {
325  unsigned int sum = 0;
326  for (int col = 0; col < cols (); col++) {
327  if (get (row, col) != 0) {
328  ++sum;
329  }
330  }
331 
332  return sum;
333 }
334 
335 Matrix Matrix::minorReduced (int rowOmit, int colOmit) const
336 {
337  ENGAUGE_ASSERT (m_rows == m_cols);
338 
339  Matrix outMinor (m_rows - 1);
340  int rowMinor = 0;
341  for (int row = 0; row < m_rows; row++) {
342 
343  if (row != rowOmit) {
344 
345  int colMinor = 0;
346  for (int col = 0; col < m_cols; col++) {
347 
348  if (col != colOmit) {
349 
350  outMinor.set (rowMinor, colMinor, get (row, col));
351  ++colMinor;
352  }
353  }
354  ++rowMinor;
355  }
356  }
357 
358  return outMinor;
359 }
360 
361 void Matrix::normalizeRow (int rowToNormalize,
362  int colToNormalize,
363  MatrixConsistent &matrixConsistent,
364  double epsilonThreshold)
365 {
366  double denominator = get (rowToNormalize, colToNormalize);
367 
368  if (valueFailsEpsilonTest (denominator,
369  epsilonThreshold)) {
370 
371  matrixConsistent = MATRIX_INCONSISTENT;
372 
373  } else {
374 
375  matrixConsistent = MATRIX_CONSISTENT;
376 
377  double factor = 1.0 / denominator;
378  for (int col = 0; col < cols (); col++) {
379  double value = get (rowToNormalize, col);
380  set (rowToNormalize, col, factor * value);
381  }
382  }
383 }
384 
385 Matrix Matrix::operator* (const Matrix &other) const
386 {
387  ENGAUGE_ASSERT (m_cols == other.rows ());
388 
389  Matrix out (m_rows, other.cols ());
390 
391  for (int row = 0; row < m_rows; row++) {
392  for (int col = 0; col < other.cols (); col++) {
393  double sum = 0;
394  for (int index = 0; index < m_cols; index++) {
395  sum += get (row, index) * other.get (index, col);
396  }
397  out.set (row, col, sum);
398  }
399  }
400 
401  return out;
402 }
403 
404 QVector<double> Matrix::operator* (const QVector<double> other) const
405 {
406  ENGAUGE_ASSERT (m_cols == other.size ());
407 
408  QVector<double> out;
409  out.resize (m_rows);
410  for (int row = 0; row < m_rows; row++) {
411  double sum = 0;
412  for (int col = 0; col < m_cols; col++) {
413  sum += get (row, col) * other [col];
414  }
415 
416  out [row] = sum;
417  }
418 
419  return out;
420 }
421 
422 int Matrix::rows () const
423 {
424  return m_rows;
425 }
426 
427 void Matrix::set (int row, int col, double value)
428 {
429  m_vector [fold2dIndexes (row, col)] = value;
430 }
431 
432 void Matrix::switchRows (int row1,
433  int row2)
434 {
435  // Switch rows
436  for (int col = 0; col < cols (); col++) {
437  double temp1 = get (row1, col);
438  double temp2 = get (row2, col);
439  set (row2, col, temp2);
440  set (row1, col, temp1);
441  }
442 }
443 
444 QString Matrix::toString () const
445 {
446  QString out;
447  QTextStream str (&out);
448 
449  str << "(";
450  for (int row = 0; row < rows (); row++) {
451  if (row > 0) {
452  str << ", ";
453  }
454  str << "(";
455  for (int col = 0; col < cols (); col++) {
456  if (col > 0) {
457  str << ", ";
458  }
459  str << get (row, col);
460  }
461  str << ")";
462  }
463  str << ")";
464 
465  return out;
466 }
467 
469 {
470  Matrix out (m_cols, m_rows);
471 
472  for (int row = 0; row < m_rows; row++) {
473  for (int col = 0; col < m_cols; col++) {
474  out.set (col, row, get (row, col));
475  }
476  }
477 
478  return out;
479 }
480 
481 bool Matrix::valueFailsEpsilonTest (double value,
482  double epsilonThreshold) const
483 {
484  return (qAbs (value) < qAbs (epsilonThreshold));
485 }
Matrix operator*(const Matrix &other) const
Multiplication operator with a matrix.
Definition: Matrix.cpp:385
Matrix transpose() const
Return the transpose of the current matrix.
Definition: Matrix.cpp:468
Matrix & operator=(const Matrix &matrix)
Assignment operator.
Definition: Matrix.cpp:35
Matrix inverse(int significantDigits, MatrixConsistent &matrixConsistent) const
Return the inverse of this matrix.
Definition: Matrix.cpp:122
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:335
double determinant() const
Return the determinant of this matrix.
Definition: Matrix.cpp:66
Matrix(int N)
Simple constructor of square matrix with initialization to identity matrix.
Definition: Matrix.cpp:12
int cols() const
Width of matrix.
Definition: Matrix.cpp:61
Matrix class that supports arbitrary NxN size.
Definition: Matrix.h:20
int rows() const
Height of matrix.
Definition: Matrix.cpp:422
double get(int row, int col) const
Return (row, col) element.
Definition: Matrix.cpp:98
void set(int row, int col, double value)
Set (row, col) element.
Definition: Matrix.cpp:427
QString toString() const
Dump matrix to a string.
Definition: Matrix.cpp:444