1 /**********************************************************************
2
3 Audacity: A Digital Audio Editor
4
5 Matrix.cpp
6
7 Dominic Mazzoni
8
9 **********************************************************************/
10
11 #include "Matrix.h"
12
13 #include <stdlib.h>
14 #include <math.h>
15
16 #include <wx/defs.h>
17
Vector()18 Vector::Vector()
19 {
20 }
21
Vector(unsigned len,double * data)22 Vector::Vector(unsigned len, double *data)
23 : mN{ len }
24 , mData(len)
25 {
26 if (data)
27 std::copy(data, data + len, mData.get());
28 else
29 std::fill(mData.get(), mData.get() + len, 0.0);
30 }
31
Vector(unsigned len,float * data)32 Vector::Vector(unsigned len, float *data)
33 : mN{ len }
34 , mData{ len }
35 {
36 if (data)
37 std::copy(data, data + len, mData.get());
38 else
39 std::fill(mData.get(), mData.get() + len, 0.0);
40 }
41
operator =(const Vector & other)42 Vector& Vector::operator=(const Vector &other)
43 {
44 wxASSERT(Len() == other.Len());
45 std::copy(other.mData.get(), other.mData.get() + mN, mData.get());
46 return *this;
47 }
48
Vector(const Vector & other)49 Vector::Vector(const Vector &other)
50 : mN{ other.Len() }
51 , mData{ mN }
52 {
53 std::copy(other.mData.get(), other.mData.get() + mN, mData.get());
54 }
55
~Vector()56 Vector::~Vector()
57 {
58 }
59
Reinit(unsigned len)60 void Vector::Reinit(unsigned len)
61 {
62 Vector temp(len);
63 Swap(temp);
64 }
65
Swap(Vector & that)66 void Vector::Swap(Vector &that)
67 {
68 std::swap(mN, that.mN);
69 mData.swap(that.mData);
70 }
71
Sum() const72 double Vector::Sum() const
73 {
74 double sum = 0.0;
75 for(unsigned i = 0; i < Len(); i++)
76 sum += mData[i];
77 return sum;
78 }
79
Matrix(unsigned rows,unsigned cols,double ** data)80 Matrix::Matrix(unsigned rows, unsigned cols, double **data)
81 : mRows{ rows }
82 , mCols{ cols }
83 , mRowVec{ mRows }
84 {
85 for(unsigned i = 0; i < mRows; i++) {
86 mRowVec[i].Reinit( mCols );
87 for(unsigned j = 0; j < mCols; j++) {
88 if (data)
89 (*this)[i][j] = data[i][j];
90 else
91 (*this)[i][j] = 0.0;
92 }
93 }
94 }
95
operator =(const Matrix & other)96 Matrix& Matrix::operator=(const Matrix &other)
97 {
98 CopyFrom(other);
99 return *this;
100 }
101
Matrix(const Matrix & other)102 Matrix::Matrix(const Matrix &other)
103 {
104 CopyFrom(other);
105 }
106
CopyFrom(const Matrix & other)107 void Matrix::CopyFrom(const Matrix &other)
108 {
109 mRows = other.mRows;
110 mCols = other.mCols;
111 mRowVec.reinit(mRows);
112 for (unsigned i = 0; i < mRows; i++) {
113 mRowVec[i].Reinit( mCols );
114 mRowVec[i] = other.mRowVec[i];
115 }
116 }
117
~Matrix()118 Matrix::~Matrix()
119 {
120 }
121
SwapRows(unsigned i,unsigned j)122 void Matrix::SwapRows(unsigned i, unsigned j)
123 {
124 mRowVec[i].Swap(mRowVec[j]);
125 }
126
IdentityMatrix(unsigned N)127 Matrix IdentityMatrix(unsigned N)
128 {
129 Matrix M(N, N);
130 for(unsigned i = 0; i < N; i++)
131 M[i][i] = 1.0;
132 return M;
133 }
134
operator +(const Vector & left,const Vector & right)135 Vector operator+(const Vector &left, const Vector &right)
136 {
137 wxASSERT(left.Len() == right.Len());
138 Vector v(left.Len());
139 for(unsigned i = 0; i < left.Len(); i++)
140 v[i] = left[i] + right[i];
141 return v;
142 }
143
operator -(const Vector & left,const Vector & right)144 Vector operator-(const Vector &left, const Vector &right)
145 {
146 wxASSERT(left.Len() == right.Len());
147 Vector v(left.Len());
148 for(unsigned i = 0; i < left.Len(); i++)
149 v[i] = left[i] - right[i];
150 return v;
151 }
152
operator *(const Vector & left,const Vector & right)153 Vector operator*(const Vector &left, const Vector &right)
154 {
155 wxASSERT(left.Len() == right.Len());
156 Vector v(left.Len());
157 for(unsigned i = 0; i < left.Len(); i++)
158 v[i] = left[i] * right[i];
159 return v;
160 }
161
operator *(const Vector & left,double right)162 Vector operator*(const Vector &left, double right)
163 {
164 Vector v(left.Len());
165 for(unsigned i = 0; i < left.Len(); i++)
166 v[i] = left[i] * right;
167 return v;
168 }
169
VectorSubset(const Vector & other,unsigned start,unsigned len)170 Vector VectorSubset(const Vector &other, unsigned start, unsigned len)
171 {
172 Vector v(len);
173 for(unsigned i = 0; i < len; i++)
174 v[i] = other[start+i];
175 return v;
176 }
177
VectorConcatenate(const Vector & left,const Vector & right)178 Vector VectorConcatenate(const Vector& left, const Vector& right)
179 {
180 Vector v(left.Len() + right.Len());
181 for(unsigned i = 0; i < left.Len(); i++)
182 v[i] = left[i];
183 for(unsigned i = 0; i < right.Len(); i++)
184 v[i + left.Len()] = right[i];
185 return v;
186 }
187
operator *(const Vector & left,const Matrix & right)188 Vector operator*(const Vector &left, const Matrix &right)
189 {
190 wxASSERT(left.Len() == right.Rows());
191 Vector v(right.Cols());
192 for(unsigned i = 0; i < right.Cols(); i++) {
193 v[i] = 0.0;
194 for(unsigned j = 0; j < right.Rows(); j++)
195 v[i] += left[j] * right[j][i];
196 }
197 return v;
198 }
199
operator *(const Matrix & left,const Vector & right)200 Vector operator*(const Matrix &left, const Vector &right)
201 {
202 wxASSERT(left.Cols() == right.Len());
203 Vector v(left.Rows());
204 for(unsigned i = 0; i < left.Rows(); i++) {
205 v[i] = 0.0;
206 for(unsigned j = 0; j < left.Cols(); j++)
207 v[i] += left[i][j] * right[j];
208 }
209 return v;
210 }
211
operator +(const Matrix & left,const Matrix & right)212 Matrix operator+(const Matrix &left, const Matrix &right)
213 {
214 wxASSERT(left.Rows() == right.Rows());
215 wxASSERT(left.Cols() == right.Cols());
216 Matrix M(left.Rows(), left.Cols());
217 for(unsigned i = 0; i < left.Rows(); i++)
218 for(unsigned j = 0; j < left.Cols(); j++)
219 M[i][j] = left[i][j] + right[i][j];
220 return M;
221 }
222
operator *(const Matrix & left,const double right)223 Matrix operator*(const Matrix &left, const double right)
224 {
225 Matrix M(left.Rows(), left.Cols());
226 for(unsigned i = 0; i < left.Rows(); i++)
227 for(unsigned j = 0; j < left.Cols(); j++)
228 M[i][j] = left[i][j] * right;
229 return M;
230 }
231
ScalarMultiply(const Matrix & left,const Matrix & right)232 Matrix ScalarMultiply(const Matrix &left, const Matrix &right)
233 {
234 wxASSERT(left.Rows() == right.Rows());
235 wxASSERT(left.Cols() == right.Cols());
236 Matrix M(left.Rows(), left.Cols());
237 for(unsigned i = 0; i < left.Rows(); i++)
238 for(unsigned j = 0; j < left.Cols(); j++)
239 M[i][j] = left[i][j] * right[i][j];
240 return M;
241 }
242
MatrixMultiply(const Matrix & left,const Matrix & right)243 Matrix MatrixMultiply(const Matrix &left, const Matrix &right)
244 {
245 wxASSERT(left.Cols() == right.Rows());
246 Matrix M(left.Rows(), right.Cols());
247 for(unsigned i = 0; i < left.Rows(); i++)
248 for(unsigned j = 0; j < right.Cols(); j++) {
249 M[i][j] = 0.0;
250 for(unsigned k = 0; k < left.Cols(); k++)
251 M[i][j] += left[i][k] * right[k][j];
252 }
253 return M;
254 }
255
MatrixSubset(const Matrix & input,unsigned startRow,unsigned numRows,unsigned startCol,unsigned numCols)256 Matrix MatrixSubset(const Matrix &input,
257 unsigned startRow, unsigned numRows,
258 unsigned startCol, unsigned numCols)
259 {
260 Matrix M(numRows, numCols);
261 for(unsigned i = 0; i < numRows; i++)
262 for(unsigned j = 0; j < numCols; j++)
263 M[i][j] = input[startRow+i][startCol+j];
264 return M;
265 }
266
MatrixConcatenateCols(const Matrix & left,const Matrix & right)267 Matrix MatrixConcatenateCols(const Matrix& left, const Matrix& right)
268 {
269 wxASSERT(left.Rows() == right.Rows());
270 Matrix M(left.Rows(), left.Cols() + right.Cols());
271 for(unsigned i = 0; i < left.Rows(); i++) {
272 for(unsigned j = 0; j < left.Cols(); j++)
273 M[i][j] = left[i][j];
274 for(unsigned j = 0; j < right.Cols(); j++)
275 M[i][j+left.Cols()] = right[i][j];
276 }
277 return M;
278 }
279
TransposeMatrix(const Matrix & other)280 Matrix TransposeMatrix(const Matrix& other)
281 {
282 Matrix M(other.Cols(), other.Rows());
283 for(unsigned i = 0; i < other.Rows(); i++)
284 for(unsigned j = 0; j < other.Cols(); j++)
285 M[j][i] = other[i][j];
286 return M;
287 }
288
InvertMatrix(const Matrix & input,Matrix & Minv)289 bool InvertMatrix(const Matrix& input, Matrix& Minv)
290 {
291 // Very straightforward implementation of
292 // Gauss-Jordan elimination to invert a matrix.
293 // Returns true if successful
294
295 wxASSERT(input.Rows() == input.Cols());
296 auto N = input.Rows();
297
298 Matrix M = input;
299 Minv = IdentityMatrix(N);
300
301 // Do the elimination one column at a time
302 for(unsigned i = 0; i < N; i++) {
303 // Pivot the row with the largest absolute value in
304 // column i, into row i
305 double absmax = 0.0;
306 unsigned int argmax = 0;
307
308 for(unsigned j = i; j < N; j++)
309 if (fabs(M[j][i]) > absmax) {
310 absmax = fabs(M[j][i]);
311 argmax = j;
312 }
313
314 // If no row has a nonzero value in that column,
315 // the matrix is singular and we have to give up.
316 if (absmax == 0)
317 return false;
318
319 if (i != argmax) {
320 M.SwapRows(i, argmax);
321 Minv.SwapRows(i, argmax);
322 }
323
324 // Divide this row by the value of M[i][i]
325 double factor = 1.0 / M[i][i];
326 M[i] = M[i] * factor;
327 Minv[i] = Minv[i] * factor;
328
329 // Eliminate the rest of the column
330 for(unsigned j = 0; j < N; j++) {
331 if (j == i)
332 continue;
333 if (fabs(M[j][i]) > 0) {
334 // Subtract a multiple of row i from row j
335 factor = M[j][i];
336 for(unsigned k = 0; k < N; k++) {
337 M[j][k] -= (M[i][k] * factor);
338 Minv[j][k] -= (Minv[i][k] * factor);
339 }
340 }
341 }
342 }
343
344 return true;
345 }
346