RDKit
Open-source cheminformatics and machine learning.
SymmMatrix.h
Go to the documentation of this file.
1 //
2 // Copyright (C) 2004-2006 Rational Discovery LLC
3 //
4 // @@ All Rights Reserved @@
5 // This file is part of the RDKit.
6 // The contents are covered by the terms of the BSD license
7 // which is included in the file license.txt, found at the root
8 // of the RDKit source tree.
9 //
10 #ifndef __RD_SYMM_MATRIX_H__
11 #define __RD_SYMM_MATRIX_H__
12 
13 #include "Matrix.h"
14 #include "SquareMatrix.h"
15 #include <cstring>
16 #include <boost/smart_ptr.hpp>
17 
18 //#ifndef INVARIANT_SILENT_METHOD
19 //#define INVARIANT_SILENT_METHOD
20 //#endif
21 namespace RDNumeric {
22 //! A symmetric matrix class
23 /*!
24  The data is stored as the lower triangle, so
25  A[i,j] = data[i*(i+1) + j] when i >= j and
26  A[i,j] = data[j*(j+1) + i] when i < j
27 */
28 template <class TYPE>
29 class SymmMatrix {
30  public:
31  typedef boost::shared_array<TYPE> DATA_SPTR;
32 
33  explicit SymmMatrix(unsigned int N) : d_size(N), d_dataSize(N * (N + 1) / 2) {
34  TYPE *data = new TYPE[d_dataSize];
35  memset(static_cast<void *>(data), 0, d_dataSize * sizeof(TYPE));
36  d_data.reset(data);
37  }
38 
39  SymmMatrix(unsigned int N, TYPE val)
40  : d_size(N), d_dataSize(N * (N + 1) / 2) {
41  TYPE *data = new TYPE[d_dataSize];
42  unsigned int i;
43  for (i = 0; i < d_dataSize; i++) {
44  data[i] = val;
45  }
46  d_data.reset(data);
47  }
48 
49  SymmMatrix(unsigned int N, DATA_SPTR data)
50  : d_size(N), d_dataSize(N * (N + 1) / 2) {
51  d_data = data;
52  }
53 
55  : d_size(other.numRows()), d_dataSize(other.getDataSize()) {
56  TYPE *data = new TYPE[d_dataSize];
57  const TYPE *otherData = other.getData();
58 
59  memcpy(static_cast<void *>(data), static_cast<const void *>(otherData),
60  d_dataSize * sizeof(TYPE));
61  d_data.reset(data);
62  }
63 
65 
66  //! returns the number of rows
67  inline unsigned int numRows() const { return d_size; }
68 
69  //! returns the number of columns
70  inline unsigned int numCols() const { return d_size; }
71 
72  inline unsigned int getDataSize() const { return d_dataSize; }
73 
74  void setToIdentity() {
75  TYPE *data = d_data.get();
76  memset(static_cast<void *>(data), 0, d_dataSize * sizeof(TYPE));
77  for (unsigned int i = 0; i < d_size; i++) {
78  data[i * (i + 3) / 2] = (TYPE)1.0;
79  }
80  }
81 
82  TYPE getVal(unsigned int i, unsigned int j) const {
83  URANGE_CHECK(i, d_size - 1);
84  URANGE_CHECK(j, d_size - 1);
85  unsigned int id;
86  if (i >= j) {
87  id = i * (i + 1) / 2 + j;
88  } else {
89  id = j * (j + 1) / 2 + i;
90  }
91  return d_data[id];
92  }
93 
94  void setVal(unsigned int i, unsigned int j, TYPE val) {
95  URANGE_CHECK(i, d_size - 1);
96  URANGE_CHECK(j, d_size - 1);
97  unsigned int id;
98  if (i >= j) {
99  id = i * (i + 1) / 2 + j;
100  } else {
101  id = j * (j + 1) / 2 + i;
102  }
103  d_data[id] = val;
104  }
105 
106  void getRow(unsigned int i, Vector<TYPE> &row) {
107  CHECK_INVARIANT(d_size == row.size(), "");
108  TYPE *rData = row.getData();
109  TYPE *data = d_data.get();
110  for (unsigned int j = 0; j < d_size; j++) {
111  unsigned int id;
112  if (j <= i) {
113  id = i * (i + 1) / 2 + j;
114  } else {
115  id = j * (j + 1) / 2 + i;
116  }
117  rData[j] = data[id];
118  }
119  }
120 
121  void getCol(unsigned int i, Vector<TYPE> &col) {
122  CHECK_INVARIANT(d_size == col.size(), "");
123  TYPE *rData = col.getData();
124  TYPE *data = d_data.get();
125  for (unsigned int j = 0; j < d_size; j++) {
126  unsigned int id;
127  if (i <= j) {
128  id = j * (j + 1) / 2 + i;
129  } else {
130  id = i * (i + 1) / 2 + j;
131  }
132  rData[j] = data[id];
133  }
134  }
135 
136  //! returns a pointer to our data array
137  inline TYPE *getData() { return d_data.get(); }
138 
139  //! returns a const pointer to our data array
140  inline const TYPE *getData() const { return d_data.get(); }
141 
143  TYPE *data = d_data.get();
144  for (unsigned int i = 0; i < d_dataSize; i++) {
145  data[i] *= scale;
146  }
147  return *this;
148  }
149 
151  TYPE *data = d_data.get();
152  for (unsigned int i = 0; i < d_dataSize; i++) {
153  data[i] /= scale;
154  }
155  return *this;
156  }
157 
159  CHECK_INVARIANT(d_size == other.numRows(),
160  "Sizes don't match in the addition");
161  const TYPE *oData = other.getData();
162  TYPE *data = d_data.get();
163  for (unsigned int i = 0; i < d_dataSize; i++) {
164  data[i] += oData[i];
165  }
166  return *this;
167  }
168 
170  CHECK_INVARIANT(d_size == other.numRows(),
171  "Sizes don't match in the addition");
172  const TYPE *oData = other.getData();
173  TYPE *data = d_data.get();
174  for (unsigned int i = 0; i < d_dataSize; i++) {
175  data[i] -= oData[i];
176  }
177  return *this;
178  }
179 
180  //! in-place matrix multiplication
183  "Size mismatch during multiplication");
184  TYPE *cData = new TYPE[d_dataSize];
185  const TYPE *bData = B.getData();
186  TYPE *data = d_data.get();
187  for (unsigned int i = 0; i < d_size; i++) {
188  unsigned int idC = i * (i + 1) / 2;
189  for (unsigned int j = 0; j < i + 1; j++) {
190  unsigned int idCt = idC + j;
191  cData[idCt] = (TYPE)0.0;
192  for (unsigned int k = 0; k < d_size; k++) {
193  unsigned int idA, idB;
194  if (k <= i) {
195  idA = i * (i + 1) / 2 + k;
196  } else {
197  idA = k * (k + 1) / 2 + i;
198  }
199  if (k <= j) {
200  idB = j * (j + 1) / 2 + k;
201  } else {
202  idB = k * (k + 1) / 2 + j;
203  }
204  cData[idCt] += (data[idA] * bData[idB]);
205  }
206  }
207  }
208 
209  for (unsigned int i = 0; i < d_dataSize; i++) {
210  data[i] = cData[i];
211  }
212  delete[] cData;
213  return (*this);
214  }
215 
216  /* Transpose will basically return a copy of itself
217  */
219  CHECK_INVARIANT(d_size == transpose.numRows(),
220  "Size mismatch during transposing");
221  TYPE *tData = transpose.getData();
222  TYPE *data = d_data.get();
223  for (unsigned int i = 0; i < d_dataSize; i++) {
224  tData[i] = data[i];
225  }
226  return transpose;
227  }
228 
230  // nothing to be done we are symmetric
231  return (*this);
232  }
233 
234  protected:
235  SymmMatrix() : d_size(0), d_dataSize(0), d_data(0){};
236  unsigned int d_size;
237  unsigned int d_dataSize;
238  DATA_SPTR d_data;
239 
240  private:
241  SymmMatrix<TYPE> &operator=(const SymmMatrix<TYPE> &other);
242 };
243 
244 //! SymmMatrix-SymmMatrix multiplication
245 /*!
246  Multiply SymmMatrix A with a second SymmMatrix B
247  and write the result to C = A*B
248 
249  \param A the first SymmMatrix
250  \param B the second SymmMatrix to multiply
251  \param C SymmMatrix to use for the results
252 
253  \return the results of multiplying A by B.
254  This is just a reference to C.
255 
256  This method is reimplemented here for efficiency reasons
257  (we basically don't want to use getter and setter functions)
258 
259 */
260 template <class TYPE>
262  SymmMatrix<TYPE> &C) {
263  unsigned int aSize = A.numRows();
264  CHECK_INVARIANT(B.numRows() == aSize,
265  "Size mismatch in matric multiplication");
266  CHECK_INVARIANT(C.numRows() == aSize,
267  "Size mismatch in matric multiplication");
268  TYPE *cData = C.getData();
269  const TYPE *aData = A.getData();
270  const TYPE *bData = B.getData();
271  for (unsigned int i = 0; i < aSize; i++) {
272  unsigned int idC = i * (i + 1) / 2;
273  for (unsigned int j = 0; j < i + 1; j++) {
274  unsigned int idCt = idC + j;
275  cData[idCt] = (TYPE)0.0;
276  for (unsigned int k = 0; k < aSize; k++) {
277  unsigned int idA, idB;
278  if (k <= i) {
279  idA = i * (i + 1) / 2 + k;
280  } else {
281  idA = k * (k + 1) / 2 + i;
282  }
283  if (k <= j) {
284  idB = j * (j + 1) / 2 + k;
285  } else {
286  idB = k * (k + 1) / 2 + j;
287  }
288  cData[idCt] += (aData[idA] * bData[idB]);
289  }
290  }
291  }
292  return C;
293 }
294 
295 //! SymmMatrix-Vector multiplication
296 /*!
297  Multiply a SymmMatrix A with a Vector x
298  so the result is y = A*x
299 
300  \param A the SymmMatrix for multiplication
301  \param x the Vector by which to multiply
302  \param y Vector to use for the results
303 
304  \return the results of multiplying x by A
305  This is just a reference to y.
306 
307  This method is reimplemented here for efficiency reasons
308  (we basically don't want to use getter and setter functions)
309 
310 */
311 template <class TYPE>
313  Vector<TYPE> &y) {
314  unsigned int aSize = A.numRows();
315  CHECK_INVARIANT(aSize == x.size(), "Size mismatch during multiplication");
316  CHECK_INVARIANT(aSize == y.size(), "Size mismatch during multiplication");
317  const TYPE *xData = x.getData();
318  const TYPE *aData = A.getData();
319  TYPE *yData = y.getData();
320  for (unsigned int i = 0; i < aSize; i++) {
321  yData[i] = (TYPE)(0.0);
322  unsigned int idA = i * (i + 1) / 2;
323  for (unsigned int j = 0; j < i + 1; j++) {
324  // idA = i*(i+1)/2 + j;
325  yData[i] += (aData[idA] * xData[j]);
326  idA++;
327  }
328  idA--;
329  for (unsigned int j = i + 1; j < aSize; j++) {
330  // idA = j*(j+1)/2 + i;
331  idA += j;
332  yData[i] += (aData[idA] * xData[j]);
333  }
334  }
335  return y;
336 }
337 
341 }
342 
343 //! ostream operator for Matrix's
344 template <class TYPE>
345 std::ostream &operator<<(std::ostream &target,
346  const RDNumeric::SymmMatrix<TYPE> &mat) {
347  unsigned int nr = mat.numRows();
348  unsigned int nc = mat.numCols();
349  target << "Rows: " << mat.numRows() << " Columns: " << mat.numCols() << "\n";
350 
351  for (unsigned int i = 0; i < nr; i++) {
352  for (unsigned int j = 0; j < nc; j++) {
353  target << std::setw(7) << std::setprecision(3) << mat.getVal(i, j);
354  }
355  target << "\n";
356  }
357  return target;
358 }
359 
360 #endif
SymmMatrix< int > IntSymmMatrix
Definition: SymmMatrix.h:339
void getRow(unsigned int i, Vector< TYPE > &row)
Definition: SymmMatrix.h:106
SymmMatrix< TYPE > & operator/=(TYPE scale)
Definition: SymmMatrix.h:150
SymmMatrix(const SymmMatrix< TYPE > &other)
Definition: SymmMatrix.h:54
#define CHECK_INVARIANT(expr, mess)
Definition: Invariant.h:99
SymmMatrix< TYPE > & operator*=(TYPE scale)
Definition: SymmMatrix.h:142
A symmetric matrix class.
Definition: SymmMatrix.h:29
SymmMatrix(unsigned int N, TYPE val)
Definition: SymmMatrix.h:39
SymmMatrix< double > DoubleSymmMatrix
Definition: SymmMatrix.h:338
unsigned int size() const
return the size (dimension) of the vector
Definition: Vector.h:77
SymmMatrix< TYPE > & transposeInplace()
Definition: SymmMatrix.h:229
SymmMatrix< TYPE > & operator*=(const SymmMatrix< TYPE > &B)
in-place matrix multiplication
Definition: SymmMatrix.h:181
unsigned int d_dataSize
Definition: SymmMatrix.h:237
unsigned int numRows() const
returns the number of rows
Definition: SymmMatrix.h:67
const TYPE * getData() const
returns a const pointer to our data array
Definition: SymmMatrix.h:140
TYPE * getData()
returns a pointer to our data array
Definition: SymmMatrix.h:137
void setVal(unsigned int i, unsigned int j, TYPE val)
Definition: SymmMatrix.h:94
SymmMatrix< TYPE > & transpose(SymmMatrix< TYPE > &transpose) const
Definition: SymmMatrix.h:218
void getCol(unsigned int i, Vector< TYPE > &col)
Definition: SymmMatrix.h:121
boost::shared_array< TYPE > DATA_SPTR
Definition: SymmMatrix.h:31
unsigned int getDataSize() const
Definition: SymmMatrix.h:72
SymmMatrix< unsigned int > UintSymmMatrix
Definition: SymmMatrix.h:340
Matrix< TYPE > & multiply(const Matrix< TYPE > &A, const Matrix< TYPE > &B, Matrix< TYPE > &C)
Matrix multiplication.
Definition: Matrix.h:254
#define URANGE_CHECK(x, hi)
Definition: Invariant.h:140
SymmMatrix< TYPE > & operator+=(const SymmMatrix< TYPE > &other)
Definition: SymmMatrix.h:158
SymmMatrix(unsigned int N)
Definition: SymmMatrix.h:33
TYPE getVal(unsigned int i, unsigned int j) const
Definition: SymmMatrix.h:82
TYPE * getData()
returns a pointer to our data array
Definition: Vector.h:102
A class to represent vectors of numbers.
Definition: Vector.h:28
unsigned int numCols() const
returns the number of columns
Definition: SymmMatrix.h:70
std::ostream & operator<<(std::ostream &target, const RDNumeric::SymmMatrix< TYPE > &mat)
ostream operator for Matrix&#39;s
Definition: SymmMatrix.h:345
unsigned int d_size
Definition: SymmMatrix.h:235
SymmMatrix(unsigned int N, DATA_SPTR data)
Definition: SymmMatrix.h:49
SymmMatrix< TYPE > & operator-=(const SymmMatrix< TYPE > &other)
Definition: SymmMatrix.h:169