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