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 <QTextStream>
10 
12 {
13  initialize (N, N);
14 }
15 
17  int cols)
18 {
19  initialize (rows, cols);
20 }
21 
22 Matrix::Matrix (const Matrix &other)
23 {
24  m_rows = other.rows();
25  m_cols = other.cols();
26  m_vector.resize (m_rows * m_cols);
27  for (int row = 0; row < m_rows; row++) {
28  for (int col = 0; col < m_cols; col++) {
29  set (row, col, other.get (row, col));
30  }
31  }
32 }
33 
35 {
36  m_rows = other.rows();
37  m_cols = other.cols();
38  m_vector.resize (m_rows * m_cols);
39  for (int row = 0; row < m_rows; row++) {
40  for (int col = 0; col < m_cols; col++) {
41  set (row, col, other.get (row, col));
42  }
43  }
44 
45  return *this;
46 }
47 
48 void Matrix::addRowToAnotherWithScaling (int rowFrom,
49  int rowTo,
50  double factor)
51 {
52  for (int col = 0; col < cols (); col++) {
53  double oldValueFrom = get (rowFrom, col);
54  double oldValueTo = get (rowTo, col);
55  double newValueTo = oldValueFrom * factor + oldValueTo;
56  set (rowTo, col, newValueTo);
57  }
58 }
59 
60 int Matrix::cols () const
61 {
62  return m_cols;
63 }
64 
65 double Matrix::determinant () const
66 {
67  ENGAUGE_ASSERT (m_rows == m_cols);
68 
69  double rtn;
70 
71  if (m_rows == 1) {
72 
73  rtn = m_vector [0];
74 
75  } else {
76 
77  const int COL = 0; // We arbitrarily iterate through the first column
78 
79  // This is a recursive algorithm
80  rtn = 0.0;
81  double multiplier = +1;
82  for (int row = 0; row < m_rows; row++) {
83  Matrix min = minorReduced (row, COL);
84  rtn += multiplier * get (row, COL) * min.determinant ();
85  multiplier *= -1.0;
86  }
87  }
88 
89  return rtn;
90 }
91 
92 int Matrix::fold2dIndexes (int row, int col) const
93 {
94  return row * m_cols + col;
95 }
96 
97 double Matrix::get (int row, int col) const
98 {
99  return m_vector [fold2dIndexes (row, col)];
100 }
101 
102 void Matrix::initialize (int rows,
103  int cols)
104 {
105  m_rows = rows;
106  m_cols = cols;
107  m_vector.resize (rows * cols);
108  for (int row = 0; row < m_rows; row++) {
109  for (int col = 0; col < m_cols; col++) {
110 
111  // Identity matrix for square matrix, otherwise zero matrix
112  if (row == col && m_rows == m_cols) {
113  set (row, col, 1.0);
114  } else {
115  set (row, col, 0.0);
116  }
117  }
118  }
119 }
120 
122 {
123  // Available algorithms are inverseCramersRule and inverseGaussianElimination
124  return inverseGaussianElimination ();
125 }
126 
127 Matrix Matrix::inverseCramersRule () const
128 {
129  ENGAUGE_ASSERT (m_rows == m_cols);
130 
131  Matrix inv (m_rows);
132  int row, col;
133 
134  if (m_rows > 1) {
135 
136  // Compute cofactor matrix
137  double multiplierStartForRow = 1.0;
138  Matrix cofactor (m_rows);
139  for (row = 0; row < m_rows; row++) {
140  double multiplier = multiplierStartForRow;
141  for (col = 0; col < m_cols; col++) {
142  Matrix min = minorReduced (row, col);
143  double element = multiplier * min.determinant ();
144  multiplier *= -1.0;
145  cofactor.set (row, col, element);
146  }
147  multiplierStartForRow *= -1.0;
148  }
149 
150  // Compute adjoint
151  Matrix adjoint = cofactor.transpose ();
152 
153  double determ = determinant ();
154 
155  // Inverse is the adjoint divided by the determinant
156  for (row = 0; row < m_rows; row++) {
157  for (col = 0; col < m_cols; col++) {
158  inv.set (row, col, adjoint.get (row, col) / determ);
159  }
160  }
161  } else {
162  inv.set (0, 0, 1.0 / get (0, 0));
163  }
164 
165  return inv;
166 }
167 
168 Matrix Matrix::inverseGaussianElimination () const
169 {
170  // From https://en.wikipedia.org/wiki/Gaussian_elimination
171 
172  int row, col, rowFrom, rowTo;
173 
174  // Track the row switches that may or may not be performed below
175  QVector<int> rowIndexes (rows ());
176  for (row = 0; row < rows (); row++) {
177  rowIndexes [row] = row;
178  }
179 
180  // Creat the working matrix and populate the left half. Number of columns in the working matrix is twice the number
181  // of cols in this matrix, but we will not populate the right half until after the bubble sort below
182  Matrix working (rows (), 2 * cols ());
183  for (row = 0; row < rows (); row++) {
184  for (col = 0; col < cols (); col++) {
185  double value = get (row, col);
186  working.set (row, col, value);
187  }
188  }
189 
190  // Simple bubble sort to rearrange the rows with 0 leading zeros at start, followed by rows with 1 leading zero at start, ...
191  for (int row1 = 0; row1 < rows (); row1++) {
192  for (int row2 = row1 + 1; row2 < rows (); row2++) {
193  if (working.leadingZeros (row1) > working.leadingZeros (row2)) {
194  working.switchRows (row1, row2);
195 
196  // Remember this switch
197  int row1Index = rowIndexes [row1];
198  int row2Index = rowIndexes [row2];
199  rowIndexes [row1] = row2Index;
200  rowIndexes [row2] = row1Index;
201  }
202  }
203  }
204 
205  // Populate the right half now
206  for (row = 0; row < rows (); row++) {
207  for (col = cols (); col < 2 * cols (); col++) {
208  int colIdentitySubmatrix = col - cols ();
209  double value = (row == colIdentitySubmatrix) ? 1 : 0;
210  working.set (row, col, value);
211  }
212  }
213 \
214  // Loop through the "from" row going down. This results in the lower off-diagonal terms becoming zero, in the left half
215  for (rowFrom = 0; rowFrom < rows (); rowFrom++) {
216  int colFirstWithNonZero = rowFrom;
217 
218  // In pathological situations we have (rowFrom, colFirstWithNonzero) = 0 in which case the solution cannot be obtained
219  // so we exit
220  if (working.get (rowFrom, colFirstWithNonZero) == 0) {
221 
222  // Fail
223  break;
224  }
225 
226  // Normalize the 'from' row with first nonzero term set to 1
227  working.normalizeRow (rowFrom, colFirstWithNonZero);
228 
229  // Apply the 'from' row to all the 'to' rows
230  for (rowTo = rowFrom + 1; rowTo < rows (); rowTo++) {
231 
232  if (working.get (rowTo, colFirstWithNonZero) != 0) {
233 
234  // We need to merge rowFrom and rowTo into rowTo
235  double factor = -1.0 * working.get (rowTo, colFirstWithNonZero) / working.get (rowFrom, colFirstWithNonZero);
236  working.addRowToAnotherWithScaling (rowFrom, rowTo, factor);
237  }
238  }
239  }
240 
241  // Loop through the "from" row doing up. This results in the upper off-diagonal terms becoming zero, in the left half
242  for (rowFrom = rows () - 1; rowFrom >= 0; rowFrom--) {
243  int colFirstWithNonZero = rowFrom; // This is true since we should have 1s all down the diagonal at this point
244 
245  // 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
246  working.normalizeRow (rowFrom, colFirstWithNonZero);
247 
248  // Apply the 'from' row to all the 'to' rows
249  for (rowTo = rowFrom - 1; rowTo >= 0; rowTo--) {
250 
251  if (working.get (rowTo, colFirstWithNonZero) != 0) {
252 
253  // We need to merge rowFro and rowTo into rowTo
254  double factor = -1.0 * working.get (rowTo, colFirstWithNonZero) / working.get (rowFrom, colFirstWithNonZero);
255  working.addRowToAnotherWithScaling (rowFrom, rowTo, factor);
256  }
257  }
258  }
259 
260  // Extract right half of rectangular matrix which is the inverse
261  Matrix inv (working.rows ());
262 
263  for (row = 0; row < working.rows (); row++) {
264 
265  int rowBeforeSort = rowIndexes [row];
266 
267  for (col = 0; col < working.rows (); col++) {
268  int colRightHalf = col + inv.cols ();
269  double value = working.get (rowBeforeSort, colRightHalf);
270  inv.set (row, col, value);
271  }
272  }
273 
274  return inv;
275 }
276 
277 unsigned int Matrix::leadingZeros (int row) const
278 {
279  unsigned int sum = 0;
280  for (int col = 0; col < cols (); col++) {
281  if (get (row, col) != 0) {
282  ++sum;
283  }
284  }
285 
286  return sum;
287 }
288 
289 Matrix Matrix::minorReduced (int rowOmit, int colOmit) const
290 {
291  ENGAUGE_ASSERT (m_rows == m_cols);
292 
293  Matrix outMinor (m_rows - 1);
294  int rowMinor = 0;
295  for (int row = 0; row < m_rows; row++) {
296 
297  if (row != rowOmit) {
298 
299  int colMinor = 0;
300  for (int col = 0; col < m_cols; col++) {
301 
302  if (col != colOmit) {
303 
304  outMinor.set (rowMinor, colMinor, get (row, col));
305  ++colMinor;
306  }
307  }
308  ++rowMinor;
309  }
310  }
311 
312  return outMinor;
313 }
314 
315 void Matrix::normalizeRow (int rowToNormalize,
316  int colToNormalize)
317 {
318  ENGAUGE_ASSERT (get (rowToNormalize, colToNormalize) != 0);
319 
320  double factor = 1.0 / get (rowToNormalize, colToNormalize);
321  for (int col = 0; col < cols (); col++) {
322  double value = get (rowToNormalize, col);
323  set (rowToNormalize, col, factor * value);
324  }
325 }
326 
327 Matrix Matrix::operator* (const Matrix &other) const
328 {
329  ENGAUGE_ASSERT (m_cols == other.rows ());
330 
331  Matrix out (m_rows, other.cols ());
332 
333  for (int row = 0; row < m_rows; row++) {
334  for (int col = 0; col < other.cols (); col++) {
335  double sum = 0;
336  for (int index = 0; index < m_cols; index++) {
337  sum += get (row, index) * other.get (index, col);
338  }
339  out.set (row, col, sum);
340  }
341  }
342 
343  return out;
344 }
345 
346 QVector<double> Matrix::operator* (const QVector<double> other) const
347 {
348  ENGAUGE_ASSERT (m_cols == other.size ());
349 
350  QVector<double> out;
351  out.resize (m_rows);
352  for (int row = 0; row < m_rows; row++) {
353  double sum = 0;
354  for (int col = 0; col < m_cols; col++) {
355  sum += get (row, col) * other [col];
356  }
357 
358  out [row] = sum;
359  }
360 
361  return out;
362 }
363 
364 int Matrix::rows () const
365 {
366  return m_rows;
367 }
368 
369 void Matrix::set (int row, int col, double value)
370 {
371  m_vector [fold2dIndexes (row, col)] = value;
372 }
373 
374 void Matrix::switchRows (int row1,
375  int row2)
376 {
377  // Switch rows
378  for (int col = 0; col < cols (); col++) {
379  double temp1 = get (row1, col);
380  double temp2 = get (row2, col);
381  set (row2, col, temp2);
382  set (row1, col, temp1);
383  }
384 }
385 
386 QString Matrix::toString () const
387 {
388  QString out;
389  QTextStream str (&out);
390 
391  str << "(";
392  for (int row = 0; row < rows (); row++) {
393  if (row > 0) {
394  str << ", ";
395  }
396  str << "(";
397  for (int col = 0; col < cols (); col++) {
398  if (col > 0) {
399  str << ", ";
400  }
401  str << get (row, col);
402  }
403  str << ")";
404  }
405  str << ")";
406 
407  return out;
408 }
409 
411 {
412  Matrix out (m_cols, m_rows);
413 
414  for (int row = 0; row < m_rows; row++) {
415  for (int col = 0; col < m_cols; col++) {
416  out.set (col, row, get (row, col));
417  }
418  }
419 
420  return out;
421 }
Matrix inverse() const
Return the inverse of this matrix.
Definition: Matrix.cpp:121
Matrix operator*(const Matrix &other) const
Multiplication operator with a matrix.
Definition: Matrix.cpp:327
Matrix transpose() const
Return the transpose of the current matrix.
Definition: Matrix.cpp:410
Matrix & operator=(const Matrix &matrix)
Assignment operator.
Definition: Matrix.cpp:34
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:289
double determinant() const
Return the determinant of this matrix.
Definition: Matrix.cpp:65
Matrix(int N)
Simple constructor of square matrix with initialization to identity matrix.
Definition: Matrix.cpp:11
int cols() const
Width of matrix.
Definition: Matrix.cpp:60
Matrix class that supports arbitrary NxN size.
Definition: Matrix.h:14
int rows() const
Height of matrix.
Definition: Matrix.cpp:364
double get(int row, int col) const
Return (row, col) element.
Definition: Matrix.cpp:97
void set(int row, int col, double value)
Set (row, col) element.
Definition: Matrix.cpp:369
QString toString() const
Dump matrix to a string.
Definition: Matrix.cpp:386