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