mathkit  1.5
 All Classes Namespaces Functions Variables Typedefs
Matrix.cpp
1 #include <sstream>
2 #include <stdexcept>
3 #include "Matrix.hpp"
4 
5 using std::abs;
6 using std::ostringstream;
7 using std::invalid_argument;
8 
9 namespace mathkit {
10 
11  Matrix::Matrix() : _data(), _dimen(0, 0) {}
12 
13  Matrix::Matrix(const Table & data) {
14  _data = data;
15  _dimen = Dimen(data[0].size(), data.size());
16  }
17 
18  Matrix::Matrix(Dimen dimen, double val) {
19  _data = Table(dimen.second, Vector(dimen.first, val));
20  _dimen = dimen;
21  }
22 
23  Matrix::Matrix(const Vector & data, Dimen dimen, bool byrow) {
24  setData(data, dimen, byrow);
25  }
26 
27  void Matrix::setData(const Table & data) {
28  _data = data;
29  _dimen = Dimen(data[0].size(), data.size());
30  }
31 
32  void Matrix::setData(Dimen dimen, double val) {
33  _data = Table(dimen.second, Vector(dimen.first, val));
34  _dimen = dimen;
35  }
36 
37  void Matrix::setData(const Vector & data, Dimen dimen, bool byrow) {
38  _dimen = dimen;
39  _data = Table(dimen.second);
40  if (byrow) {
41  int k = 0;
42  for (int i = 0; i < dimen.first; ++i)
43  for (int j = 0; j < dimen.second; ++j)
44  _data[j].push_back(data[k++]);
45  }
46  else {
47  int k = 0;
48  for (int i = 0; i < dimen.second; ++i)
49  for (int j = 0; j < dimen.first; ++j)
50  _data[i].push_back(data[k++]);
51  }
52  }
53 
55  return _data;
56  }
57 
58  Dimen Matrix::dimen() const {
59  return _dimen;
60  }
61 
62  bool Matrix::dimen(Dimen dimen) {
63  if (!_data.empty()) return false;
64  _dimen = dimen;
65  return true;
66  }
67 
68  int Matrix::nrow() const {
69  return _dimen.first;
70  }
71 
72  int Matrix::ncol() const {
73  return _dimen.second;
74  }
75 
76  Matrix Matrix::round(int ndec) const {
77  Matrix mat = *this;
78  int m = mat._dimen.first, n = mat._dimen.second;
79  for (int i = 0; i < m; ++i)
80  for (int j = 0; j < n; ++j)
81  mat(i, j) = prec(mat(i, j), ndec);
82  return mat;
83  }
84 
85  bool Matrix::square() const {
86  return _dimen.first == _dimen.second;
87  }
88 
89  bool Matrix::symmetric() const {
90  if (!square()) return false;
91  return *this == transpose();
92  }
93 
94  Vector Matrix::getRow(int m) const {
95  Vector row(_dimen.second);
96  for (int i = 0; i < _dimen.second; ++i) row[i] = _data[i][m];
97  return row;
98  }
99 
100  Vector Matrix::getCol(int n) const {
101  return _data[n];
102  }
103 
104  void Matrix::setRow(int m, const Vector & row) {
105  for (int i = 0; i < _dimen.second; ++i) _data[i][m] = row[i];
106  }
107 
108  void Matrix::setCol(int n, const Vector & col) {
109  _data[n] = col;
110  }
111 
112  void Matrix::exchange(int i, int j, bool row) {
113  Vector temp;
114  if (row) {
115  temp = getRow(i);
116  setRow(i, getRow(j));
117  setRow(j, temp);
118  }
119  else {
120  temp = getCol(i);
121  setCol(i, getCol(j));
122  setCol(j, temp);
123  }
124  }
125 
127  Table tbl;
128  for (int i = 0; i < _dimen.first; ++i) tbl.push_back(getRow(i));
129  return Matrix(tbl);
130  }
131 
132  int Matrix::rank(const Epsilon & eps) const {
133  Matrix mat = *this;
134  int m = mat._dimen.first;
135  int n = mat._dimen.second;
136  int limit = m < n ? m : n;
137  int r = 0;
138  for (int i = 0; i < n; ++i) {
139  if (eps(mat(r, i))) {
140  bool exchanged = false;
141  for (int j = r + 1; j < m; ++j) {
142  if (!eps(mat(j, i))) {
143  mat.exchange(r, j);
144  exchanged = true;
145  break;
146  }
147  }
148  if (!exchanged) continue;
149  }
150  for (int k = r + 1; k < m; ++k) {
151  if (!eps(mat(k, i))) {
152  double s = mat(k, i) / mat(r, i);
153  mat.setRow(k, -s * mat.getRow(r) + mat.getRow(k));
154  }
155  }
156  r += 1;
157  if (r == limit) break;
158  }
159  return r;
160  }
161 
162  Matrix Matrix::submat(int frow, int lrow, int fcol, int lcol) {
163  int m = lrow - frow + 1, n = lcol - fcol + 1;
164  Table tbl(n, Vector(m));
165  for (int i = 0; i < m; ++i)
166  for (int j = 0; j < n; ++j)
167  tbl[j][i] = _data[fcol + j][frow + i];
168  return Matrix(tbl);
169  }
170 
171  Matrix Matrix::resid(int i, int j) const {
172  Table tbl;
173  for (int n = 0; n < _dimen.second; ++n) {
174  if (n != j) {
175  Vector col = _data[n];
176  Vector::iterator iter = col.begin();
177  for (int m = 0; m < i; ++m) ++iter;
178  col.erase(iter);
179  tbl.push_back(col);
180  }
181  }
182  return Matrix(tbl);
183  }
184 
185  double & Matrix::operator()(int m, int n) {
186  return _data[n][m];
187  }
188 
189  double Matrix::operator()(int m, int n) const {
190  return _data[n][m];
191  }
192 
193  Matrix Matrix::operator+(const Matrix & mat) const {
194  if (_dimen != mat._dimen) throw invalid_argument("Incompatible matrices.");
195  Table tbl;
196  for (int i = 0; i < _dimen.second; ++i) tbl.push_back(_data[i] + mat._data[i]);
197  return Matrix(tbl);
198  }
199 
200  Matrix Matrix::operator-(const Matrix & mat) const {
201  if (_dimen != mat._dimen) throw invalid_argument("Incompatible matrices.");
202  Table tbl;
203  for (int i = 0; i < _dimen.second; ++i) tbl.push_back(_data[i] - mat._data[i]);
204  return Matrix(tbl);
205  }
206 
207  Matrix Matrix::operator*(const Matrix & mat) const {
208  if (_dimen.second != mat._dimen.first) throw invalid_argument("Incompatible matrices.");
209  Table tbl(mat._dimen.second, Vector(_dimen.first, 0));
210  for (int i = 0; i < mat._dimen.second; ++i)
211  for (int j = 0; j < _dimen.first; ++j)
212  tbl[i] = tbl[i] + mat._data[i][j] * _data[j];
213  return Matrix(tbl);
214  }
215 
216  bool Matrix::operator==(const Matrix & mat) const {
217  return _dimen == mat._dimen && _data == mat._data;
218  }
219 
220  bool Matrix::operator!=(const Matrix & mat) const {
221  return !(_dimen == mat._dimen && _data == mat._data);
222  }
223 
224  Matrix operator*(const Matrix & mat, double scalar) {
225  Table tbl;
226  int n = mat.ncol();
227  for (int i = 0; i < n; ++i) tbl.push_back(scalar * mat.getCol(i));
228  return Matrix(tbl);
229  }
230 
231  Matrix operator*(double scalar, const Matrix & mat) {
232  return mat * scalar;
233  }
234 
235  ostream & operator<<(ostream & outs, const Matrix & mat) {
236  save(outs, mat.getData());
237  return outs;
238  }
239 
240  istream & operator>>(istream & ins, Matrix & mat) {
241  if (mat.dimen() != Dimen(0, 0))
242  mat.setData(load(ins, mat.nrow(), mat.ncol()));
243  return ins;
244  }
245 
246  Matrix t(const Matrix & mat) {
247  return mat.transpose();
248  }
249 
250  Matrix eye(int n) {
251  Matrix mat(Dimen(n, n));
252  for (int i = 0; i < n; ++i) mat(i, i) = 1;
253  return mat;
254  }
255 
256  Vector diag(const Matrix & mat) {
257  Vector vec;
258  int n = mat.dimen().second;
259  for (int i = 0; i < n; ++i) vec.push_back(mat(i, i));
260  return vec;
261  }
262 
263  Matrix diag(const Vector & vec) {
264  int n = vec.size();
265  Matrix mat(Dimen(n, n));
266  for (int i = 0; i < n; ++i) mat(i, i) = vec[i];
267  return mat;
268  }
269 
270  double det(const Matrix & mat, const Epsilon & eps) {
271  if (!mat.square()) throw invalid_argument("Square matrix needed.");
272  int n = mat.ncol();
273  if (n == 1) return eps(mat(0, 0)) ? 0 : mat(0, 0);
274  double result = 0;
275  for (int i = 0; i < n; ++i) {
276  if (!eps(mat(0, i))) {
277  double resid = mat(0, i) * det(mat.resid(0, i));
278  if (i % 2 != 0) resid = -resid;
279  result += resid;
280  }
281  }
282  return result;
283  }
284 
285  Matrix inv(const Matrix & mat, const Epsilon & eps) {
286  if (!mat.square()) throw invalid_argument("Square matrix needed.");
287  int n = mat.ncol();
288  Matrix I = eye(n);
289  Table tbl;
290  Matrices lup = lu(mat, eps);
291  for (int i = 0; i < n; ++i) {
292  Vector result = lusolve(lup, I.getCol(i), eps);
293  if (result.empty()) throw invalid_argument("Singular matrix.");
294  tbl.push_back(result);
295  }
296  return Matrix(tbl);
297  }
298 
299  double trace(const Matrix & mat, const Epsilon & eps) {
300  if (!mat.square()) throw invalid_argument("Square matrix needed.");
301  int n = mat.ncol();
302  double val = 0;
303  for (int i = 0; i < n; ++i)
304  if (!eps(mat(i, i))) val += mat(i, i);
305  return val;
306  }
307 
308  bool positive(const Matrix & mat, const Epsilon & eps) {
309  if (!mat.symmetric()) throw invalid_argument("Symmetric matrix needed.");
310  Matrices lup = lu(mat, eps);
311  if (lup.empty()) return false;
312  const Matrix &U = lup[1];
313  int n = U.ncol();
314  for (int i = 0; i < n; ++i)
315  if (U(i, i) < 0) return false;
316  return true;
317  }
318 
319  Table table(const Vector & data, Dimen dimen, bool byrow) {
320  Table tbl(dimen.second);
321  if (byrow) {
322  int k = 0;
323  for (int i = 0; i < dimen.first; ++i)
324  for (int j = 0; j < dimen.second; ++j)
325  tbl[j].push_back(data[k++]);
326  }
327  else {
328  int k = 0;
329  for (int i = 0; i < dimen.second; ++i)
330  for (int j = 0; j < dimen.first; ++j)
331  tbl[i].push_back(data[k++]);
332  }
333  return tbl;
334  }
335 
336  Matrix rowmat(const Vector & data) {
337  return Matrix(data, Dimen(1, data.size()));
338  }
339 
340  Matrix colmat(const Vector & data) {
341  return Matrix(data, Dimen(data.size(), 1));
342  }
343 
344  Matrix cbind(const Matrix & mat1, const Matrix & mat2) {
345  if (mat1.nrow() != mat2.nrow()) throw invalid_argument("Incompatible matrices.");
346  Table tbl;
347  for (int i = 0; i < mat1.ncol(); ++i) tbl.push_back(mat1.getCol(i));
348  for (int i = 0; i < mat2.ncol(); ++i) tbl.push_back(mat2.getCol(i));
349  return Matrix(tbl);
350  }
351 
352  Matrix rbind(const Matrix & mat1, const Matrix & mat2) {
353  int n = mat1.ncol();
354  if (n != mat2.ncol()) throw invalid_argument("Incompatible matrices.");
355  Table tbl(n);
356  for (int i = 0; i < mat1.nrow(); ++i) {
357  Vector row = mat1.getRow(i);
358  for (int j = 0; j < n; ++j) tbl[j].push_back(row[j]);
359  }
360  for (int i = 0; i < mat2.nrow(); ++i) {
361  Vector row = mat2.getRow(i);
362  for (int j = 0; j < n; ++j) tbl[j].push_back(row[j]);
363  }
364  return Matrix(tbl);
365  }
366 
367  Matrices lu(const Matrix & mat, const Epsilon & eps) {
368  if (!mat.square()) throw invalid_argument("Square matrix needed.");
369  int n = mat.ncol();
370  Matrix L = eye(n), P = eye(n);
371  Matrix U = mat;
372  for (int i = 0; i < n; ++i) {
373  if (eps(U(i, i))) {
374  bool exchanged = false;
375  for (int k = i + 1; k < n; ++k) {
376  if (!eps(U(k, i))) {
377  U.exchange(i, k);
378  P.exchange(i, k);
379  for (int j = 0; j < i; ++j) {
380  double temp = L(i, j);
381  L(i, j) = L(k, j);
382  L(k, j) = temp;
383  }
384  exchanged = true;
385  break;
386  }
387  }
388  if (!exchanged) return Matrices();
389  }
390  for (int j = i + 1; j < n; ++j) {
391  if (!eps(U(j, i))) {
392  double s = U(j, i) / U(i, i);
393  L(j, i) = s;
394  U.setRow(j, -s * U.getRow(i) + U.getRow(j));
395  }
396  }
397  }
398  Matrices result;
399  result.push_back(L);
400  result.push_back(U);
401  result.push_back(P);
402  return result;
403  }
404 
405  Matrices qr(const Matrix & mat, const Epsilon & eps) {
406  int n = mat.ncol();
407  if (mat.rank(eps) < n) return Matrices();
408  Table tbl;
409  for (int i = 0; i < n; ++i) {
410  Vector vec = mat.getCol(i);
411  for (int j = 0; j < i; ++j) vec = vec - tbl[j] * dot_prod(tbl[j], vec);
412  tbl.push_back(vec / norm(vec));
413  }
414  Matrix Q(tbl);
415  Matrix R(Dimen(n, n));
416  for (int i = 0; i < n; ++i)
417  for (int j = i; j < n; ++j)
418  R(i, j) = dot_prod(tbl[i], mat.getCol(j));
419  Matrices result;
420  result.push_back(Q);
421  result.push_back(R);
422  return result;
423  }
424 
425  Vector linsolve(const Matrix & mat, const Vector & vec, const Epsilon & eps) {
426  return lusolve(lu(mat, eps), vec, eps);
427  }
428 
429  Vector lusolve(const Matrices & lup, const Vector & vec, const Epsilon & eps) {
430  if (lup.empty()) return Vector();
431  const Matrix &L = lup[0], &U = lup[1], &P = lup[2];
432  unsigned int n = U.nrow();
433  if (n != vec.size()) throw invalid_argument("Incompatible arguments.");
434  Matrix b = colmat(vec);
435  if (P != eye(n)) b = P * b;
436  Vector result;
437  for (unsigned int i = 0; i < n; ++i) {
438  double val = 0;
439  for (unsigned int j = 0; j < i; ++j)
440  if (!eps(L(i, j)))val += L(i, j) * result[j];
441  result.push_back(b(i, 0) - val);
442  }
443  b = colmat(result);
444  for (int i = n - 1; i >= 0; --i) {
445  double val = 0;
446  for (int j = n - 1; j > i; --j)
447  if (!eps(U(i, j))) val += U(i, j) * result[j];
448  result[i] = (b(i, 0) - val) / U(i, i);
449  }
450  return result;
451  }
452 
453  Vector qrsolve(const Matrices & qrmat, const Vector & vec, const Epsilon & eps) {
454  if (qrmat.empty()) return Vector();
455  const Matrix &Q = qrmat[0], &R = qrmat[1];
456  unsigned int m = Q.nrow();
457  if (vec.size() != m) throw invalid_argument("Incompatible arguments.");
458  Matrix b = colmat(vec);
459  b = t(Q) * b;
460  int n = R.nrow();
461  Vector result(n);
462  for (int i = n - 1; i >= 0; --i) {
463  double val = 0;
464  for (int j = n - 1; j > i; --j)
465  if (!eps(R(i, j))) val += R(i, j) * result[j];
466  result[i] = (b(i, 0) - val) / R(i, i);
467  }
468  return result;
469  }
470 
471  string to_str(const Matrix & mat, bool byrow, string delim) {
472  ostringstream oss;
473  if (byrow) {
474  int m = mat.nrow();
475  for (int i = 0; i < m; ++i) oss << to_str(mat.getRow(i));
476  }
477  else {
478  int n = mat.ncol();
479  for (int i = 0; i < n; ++i) oss << to_str(mat.getCol(i));
480  }
481  string result = oss.str();
482  if (delim.size() == 2) result = delim[0] + result + delim[1];
483  return result;
484  }
485 }