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
SymmMatrix(unsigned int N)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
SymmMatrix(unsigned int N,TYPE val)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
SymmMatrix(unsigned int N,DATA_SPTR data)50 SymmMatrix(unsigned int N, DATA_SPTR data)
51 : d_size(N), d_dataSize(N * (N + 1) / 2) {
52 d_data = data;
53 }
54
SymmMatrix(const SymmMatrix<TYPE> & other)55 SymmMatrix(const SymmMatrix<TYPE> &other)
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
~SymmMatrix()65 ~SymmMatrix() {}
66
67 //! returns the number of rows
numRows()68 inline unsigned int numRows() const { return d_size; }
69
70 //! returns the number of columns
numCols()71 inline unsigned int numCols() const { return d_size; }
72
getDataSize()73 inline unsigned int getDataSize() const { return d_dataSize; }
74
setToIdentity()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
getVal(unsigned int i,unsigned int j)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
setVal(unsigned int i,unsigned int j,TYPE val)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
getRow(unsigned int i,Vector<TYPE> & row)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
getCol(unsigned int i,Vector<TYPE> & col)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
getData()138 inline TYPE *getData() { return d_data.get(); }
139
140 //! returns a const pointer to our data array
getData()141 inline const TYPE *getData() const { return d_data.get(); }
142
143 SymmMatrix<TYPE> &operator*=(TYPE scale) {
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
151 SymmMatrix<TYPE> &operator/=(TYPE scale) {
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
159 SymmMatrix<TYPE> &operator+=(const SymmMatrix<TYPE> &other) {
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
170 SymmMatrix<TYPE> &operator-=(const SymmMatrix<TYPE> &other) {
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
182 SymmMatrix<TYPE> &operator*=(const SymmMatrix<TYPE> &B) {
183 CHECK_INVARIANT(d_size == B.numRows(),
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 */
transpose(SymmMatrix<TYPE> & transpose)219 SymmMatrix<TYPE> &transpose(SymmMatrix<TYPE> &transpose) const {
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
transposeInplace()230 SymmMatrix<TYPE> &transposeInplace() {
231 // nothing to be done we are symmetric
232 return (*this);
233 }
234
235 protected:
SymmMatrix()236 SymmMatrix() : d_data(0){};
237 unsigned int d_size{0};
238 unsigned int d_dataSize{0};
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>
multiply(const SymmMatrix<TYPE> & A,const SymmMatrix<TYPE> & B,SymmMatrix<TYPE> & C)262 SymmMatrix<TYPE> &multiply(const SymmMatrix<TYPE> &A, const SymmMatrix<TYPE> &B,
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>
multiply(const SymmMatrix<TYPE> & A,const Vector<TYPE> & x,Vector<TYPE> & y)313 Vector<TYPE> &multiply(const SymmMatrix<TYPE> &A, const Vector<TYPE> &x,
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
339 typedef SymmMatrix<double> DoubleSymmMatrix;
340 typedef SymmMatrix<int> IntSymmMatrix;
341 typedef SymmMatrix<unsigned int> UintSymmMatrix;
342 } // namespace RDNumeric
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
362