CluE  1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
matrix.cpp
Go to the documentation of this file.
1 #include "../point/matrix.h"
2 
3 #include "../exception/invalidargumentexception.h"
4 #include "../exception/invalidruntimeconfigurationexception.h"
5 
6 #include <math.h>
7 #include <iostream>
8 #include <numeric>
9 #include <vector>
10 #include <algorithm>
11 
12 using namespace CluE;
13 
14 Matrix::Matrix(size_t r, size_t c, std::vector<double> const & e):rows(r),cols(c)
15 {
16  if (e.size() != r*c)
17  throw InvalidArgumentException(0, "Incompatible dimensions!", "e");
18 
19  this->entries = e;
20 }
21 
22 Matrix::Matrix(Point const & p, bool isRow /*= false*/)
23 {
24  size_t d = p.dimension();
25 
26  if (isRow)
27  {
28  this->rows = 1;
29  this->cols = d;
30  } else {
31  this->rows = d;
32  this->cols = 1;
33  }
34 
35  this->entries = std::vector<double>(d);
36  for (size_t i=0; i<d; ++i)
37  this->entries[i] = p[i];
38 }
39 
41 {
42  Matrix m(this->cols, this->rows);
43  for (size_t i=0; i<this->rows; ++i)
44  for (size_t j=0; j<this->cols; ++j)
45  m(j,i) = (*this)(i,j);
46  return m;
47 }
48 
50 {
51  if (this->rows!=this->cols)
52  throw InvalidRuntimeConfigurationException(0, "Matrix is not quadratic.");
53 
54  size_t n = this->rows;
55  Matrix c(n,n);
56 
57  for (size_t i=0; i<n; ++i)
58  {
59  for (size_t j=0; j<i; ++j)
60  {
61  double sum = (*this)(i,j);
62  for (size_t k=0; k<j; ++k)
63  sum -= c(i,k)*c(j,k);
64  c(i,j) = sum/c(j,j);
65  }
66 
67  double sum = (*this)(i,i);
68  for (size_t k=0; k<i; ++k)
69  sum -= c(i,k)*c(i,k);
70  if (sum > 0)
71  c(i,i) = sqrt(sum);
72  else
73  {
74  std::cout << "cholesky(): i = " << i << " sum = " << sum << std::endl;
75  throw InvalidRuntimeConfigurationException(1, "Matrix is not positive definite.");
76  }
77  }
78 
79  return c;
80 }
81 
82 bool Matrix::isSPD() const
83 {
84  if (this->rows!=this->cols)
85  return false;
86 
87  size_t n = this->rows;
88  Matrix c(n,n);
89 
90  for (size_t i=0; i<n; ++i)
91  {
92  for (size_t j=0; j<i; ++j)
93  {
94  double sum = (*this)(i,j);
95  for (size_t k=0; k<j; ++k)
96  sum -= c(i,k)*c(j,k);
97  c(i,j) = sum/c(j,j);
98  }
99 
100  double sum = (*this)(i,i);
101  for (size_t k=0; k<i; ++k)
102  sum -= c(i,k)*c(i,k);
103  if (sum > 0)
104  c(i,i) = sqrt(sum);
105  else
106  return false;
107  }
108 
109  return true;
110 }
111 
113 {
114  Matrix m = this->cholesky();
115 
116  size_t n = m.numRows();
117  double det = 1;
118  for (size_t i=0; i<n; ++i)
119  det *= m(i,i)*m(i,i);
120 
121  return det;
122 }
123 
125 {
126  Matrix ic = ltInverse(this->cholesky());
127 
128  // multiplication of inverse of cholesky decomposition with its transpose
129  size_t n = this->rows;
130  Matrix inv(n,n);
131  for (size_t i=0; i<n; ++i)
132  for (size_t j=0; j<n; ++j)
133  {
134  double sum = 0;
135  for (size_t k=0; k<=i && k<=j; ++k)
136  {
137  sum += ic(i,k)*ic(j,k);
138  }
139  inv(i,j) = sum;
140  }
141 
142  return inv;
143 }
144 
145 
147 {
148  size_t r = this->numRows();
149  size_t c = this->numColumns();
150 
151  Matrix gs(*this);
152  for (size_t i=0; i<c; ++i)
153  {
154  for (size_t j=0; j<i; ++j)
155  {
156  double s = 0;
157  for (size_t k=0; k<r; ++k)
158  s += gs(k,j)*(*this)(k,i);
159  for (size_t k=0; k<r; ++k)
160  gs(k,i) -= s*gs(k,j);
161  }
162  double sum = 0;
163  for (size_t k=0; k<r; ++k)
164  sum += gs(k,i)*gs(k,i);
165  for (size_t k=0; k<r; ++k)
166  gs(k,i) /= sqrt(sum);
167  }
168 
169  return gs;
170 }
172 {
173  size_t r = this->numRows();
174  size_t c = this->numColumns();
175 
176  if(m.numRows() != r || m.numColumns() != c)
177  throw InvalidArgumentException(0, "Incompatible dimensions!", "m");
178 
179  for(size_t i = 0; i < r; ++i)
180  for(size_t j = 0; j < c; ++j)
181  (*this)(i,j) += m(i,j);
182 
183  return *this;
184 }
185 
187 {
188  size_t r = this->numRows();
189  size_t c = this->numColumns();
190 
191  if(m.numRows() != r || m.numColumns() != c)
192  throw InvalidArgumentException(0, "Incompatible dimensions!", "m");
193 
194  for(size_t i = 0; i < r; ++i)
195  for(size_t j = 0; j < c; ++j)
196  (*this)(i,j) -= m(i,j);
197 
198  return *this;
199 }
200 
202 {
203  return Matrix(*this) += m;
204 }
205 
207 {
208  return Matrix(*this) -= m;
209 }
210 
212 {
213  (*this)=(*this)*m;
214  return *this;
215 }
216 
218 {
219  size_t r = this->numRows();
220  size_t n = this->numColumns();
221  size_t c = m.numColumns();
222 
223  if(m.numRows() != n)
224  throw InvalidArgumentException(0, "Incompatible dimensions!", "m");
225 
226  Matrix ret(r,c);
227  for(size_t i = 0; i < r; ++i)
228  for(size_t j = 0; j < c; ++j)
229  for(size_t k = 0; k < n; ++k)
230  ret(i,j) += (*this)(i,k)*m(k,j);
231 
232  return ret;
233 }
234 
235 std::ostream& CluE::operator<<(std::ostream& os, Matrix const& m)
236 {
237  size_t r = m.numRows();
238  size_t c = m.numColumns();
239  os << "(";
240  for(size_t i=0; i<r; ++i)
241  {
242  os << "(";
243  for(size_t j=0; j<c; ++j)
244  {
245  os << m(i,j);
246  if(j<c-1)
247  os << ", ";
248  }
249  os << ")";
250  if(i<r-1)
251  os << ", ";
252  }
253  os << ")";
254  return os;
255 }
256 
257 Matrix CluE::operator*(double scalar, Matrix const & m)
258 {
259  Matrix ret(m);
260 
261  size_t r = ret.numRows();
262  size_t c = ret.numColumns();
263 
264  for(size_t i = 0; i < r; ++i)
265  for(size_t j = 0; j < c; ++j)
266  ret(i,j) = scalar * m(i,j);
267 
268  return ret;
269 }
270 
271 Point CluE::operator*(Point const& row, Matrix const& m)
272 {
273  size_t r = m.numRows();
274  size_t c = m.numColumns();
275 
276  if(row.dimension() != r)
277  throw InvalidArgumentException(0, "Incompatible dimensions!", "row, m");
278 
279  Point ret(c);
280  for(size_t j = 0; j < c; ++j)
281  for(size_t k = 0; k < r; ++k)
282  ret[j] += row[k]*m(k,j);
283 
284  return ret;
285 }
286 
287 Point CluE::operator*(Matrix const& m, Point const& column)
288 {
289  size_t r = m.numRows();
290  size_t c = m.numColumns();
291 
292  if(column.dimension() != c)
293  throw InvalidArgumentException(0, "Incompatible dimensions!", "column, m");
294 
295  Point ret(r);
296  for(size_t i = 0; i < r; ++i)
297  for(size_t k = 0; k < c; ++k)
298  ret[i] += m(i,k)*column[k];
299 
300  return ret;
301 }
302 
304 {
305  Matrix m(d, d);
306  for (size_t i=0; i<d; ++i)
307  m(i,i)=1;
308 
309  return m;
310 }
311 
313 {
314  size_t n = lt.rows;
315  Matrix inv(n,n);
316 
317  // compute inverse of lt (backwards substitution)
318  for (size_t i=n; i-- > 0; )
319  {
320  inv(i,i) = 1/lt(i,i);
321  for (size_t j=i+1; j<n; ++j)
322  {
323  double sum = 0;
324  for (size_t k=i; k<j; ++k)
325  {
326  sum += lt(j,k)*inv(k,i);
327  }
328  inv(j,i) = -sum*inv(j,j);
329  }
330  }
331 
332  return inv;
333 }
Matrix operator-(Matrix const &m) const
Definition: matrix.cpp:206
Matrix spdInverse() const
Definition: matrix.cpp:124
size_t numColumns() const
Definition: matrix.h:110
size_t dimension() const
Definition: point.h:87
static Matrix ltInverse(Matrix const &lt)
Definition: matrix.cpp:312
size_t rows
Definition: matrix.h:120
Matrix & operator-=(Matrix const &m)
Definition: matrix.cpp:186
Matrix operator*(Matrix const &m) const
Definition: matrix.cpp:217
Weighted matrix of arbitrary dimension.
Definition: matrix.h:17
std::vector< double > entries
Definition: matrix.h:121
Matrix & operator*=(Matrix const &m)
Definition: matrix.cpp:211
Matrix(size_t r=0, size_t c=0)
Constructs a matrix.
Definition: matrix.h:23
Matrix cholesky() const
Definition: matrix.cpp:49
Matrix & operator+=(Matrix const &m)
Definition: matrix.cpp:171
size_t numRows() const
Definition: matrix.h:102
double spdDeterminant() const
Definition: matrix.cpp:112
bool isSPD() const
Definition: matrix.cpp:82
Matrix operator*(double scalar, Matrix const &m)
Definition: matrix.cpp:257
static Matrix identity(double dimension)
Definition: matrix.cpp:303
Weighted point of arbitrary dimension.
Definition: point.h:17
Indicates invalid values of arguments.
Matrix transpose() const
Definition: matrix.cpp:40
Matrix gramSchmidt() const
Definition: matrix.cpp:146
Matrix operator+(Matrix const &m) const
Definition: matrix.cpp:201
Indicates that a computation entered an invalid configuration state.
std::ostream & operator<<(std::ostream &, FrequencyDistribution &)
size_t cols
Definition: matrix.h:120