Loading...
Searching...
No Matches
SymmMatrix.h
Go to the documentation of this file.
//
// Copyright (C) 2004-2006 Rational Discovery LLC
//
// @@ All Rights Reserved @@
// This file is part of the RDKit.
// The contents are covered by the terms of the BSD license
// which is included in the file license.txt, found at the root
// of the RDKit source tree.
//
#include <RDGeneral/export.h>
#ifndef __RD_SYMM_MATRIX_H__
#define __RD_SYMM_MATRIX_H__
#include "Matrix.h"
#include "SquareMatrix.h"
#include <cstring>
#include <boost/smart_ptr.hpp>
// #ifndef INVARIANT_SILENT_METHOD
// #define INVARIANT_SILENT_METHOD
// #endif
namespace RDNumeric {
//! A symmetric matrix class
/*!
The data is stored as the lower triangle, so
A[i,j] = data[i*(i+1) + j] when i >= j and
A[i,j] = data[j*(j+1) + i] when i < j
*/
template <class TYPE>
class SymmMatrix {
public:
typedef boost::shared_array<TYPE> DATA_SPTR;
explicit SymmMatrix(unsigned int N) : d_size(N), d_dataSize(N * (N + 1) / 2) {
TYPE *data = new TYPE[d_dataSize];
memset(static_cast<void *>(data), 0, d_dataSize * sizeof(TYPE));
d_data.reset(data);
}
SymmMatrix(unsigned int N, TYPE val)
: d_size(N), d_dataSize(N * (N + 1) / 2) {
TYPE *data = new TYPE[d_dataSize];
unsigned int i;
for (i = 0; i < d_dataSize; i++) {
data[i] = val;
}
d_data.reset(data);
}
SymmMatrix(unsigned int N, DATA_SPTR data)
: d_size(N), d_dataSize(N * (N + 1) / 2) {
d_data = data;
}
SymmMatrix(const SymmMatrix<TYPE> &other)
: d_size(other.numRows()), d_dataSize(other.getDataSize()) {
TYPE *data = new TYPE[d_dataSize];
const TYPE *otherData = other.getData();
memcpy(static_cast<void *>(data), static_cast<const void *>(otherData),
d_dataSize * sizeof(TYPE));
d_data.reset(data);
}
~SymmMatrix() = default;
//! returns the number of rows
inline unsigned int numRows() const { return d_size; }
//! returns the number of columns
inline unsigned int numCols() const { return d_size; }
inline unsigned int getDataSize() const { return d_dataSize; }
void setToIdentity() {
TYPE *data = d_data.get();
memset(static_cast<void *>(data), 0, d_dataSize * sizeof(TYPE));
for (unsigned int i = 0; i < d_size; i++) {
data[i * (i + 3) / 2] = (TYPE)1.0;
}
}
TYPE getVal(unsigned int i, unsigned int j) const {
URANGE_CHECK(i, d_size);
URANGE_CHECK(j, d_size);
unsigned int id;
if (i >= j) {
id = i * (i + 1) / 2 + j;
} else {
id = j * (j + 1) / 2 + i;
}
return d_data[id];
}
void setVal(unsigned int i, unsigned int j, TYPE val) {
URANGE_CHECK(i, d_size);
URANGE_CHECK(j, d_size);
unsigned int id;
if (i >= j) {
id = i * (i + 1) / 2 + j;
} else {
id = j * (j + 1) / 2 + i;
}
d_data[id] = val;
}
void getRow(unsigned int i, Vector<TYPE> &row) {
CHECK_INVARIANT(d_size == row.size(), "");
TYPE *rData = row.getData();
TYPE *data = d_data.get();
for (unsigned int j = 0; j < d_size; j++) {
unsigned int id;
if (j <= i) {
id = i * (i + 1) / 2 + j;
} else {
id = j * (j + 1) / 2 + i;
}
rData[j] = data[id];
}
}
void getCol(unsigned int i, Vector<TYPE> &col) {
CHECK_INVARIANT(d_size == col.size(), "");
TYPE *rData = col.getData();
TYPE *data = d_data.get();
for (unsigned int j = 0; j < d_size; j++) {
unsigned int id;
if (i <= j) {
id = j * (j + 1) / 2 + i;
} else {
id = i * (i + 1) / 2 + j;
}
rData[j] = data[id];
}
}
//! returns a pointer to our data array
inline TYPE *getData() { return d_data.get(); }
//! returns a const pointer to our data array
inline const TYPE *getData() const { return d_data.get(); }
SymmMatrix<TYPE> &operator*=(TYPE scale) {
TYPE *data = d_data.get();
for (unsigned int i = 0; i < d_dataSize; i++) {
data[i] *= scale;
}
return *this;
}
SymmMatrix<TYPE> &operator/=(TYPE scale) {
TYPE *data = d_data.get();
for (unsigned int i = 0; i < d_dataSize; i++) {
data[i] /= scale;
}
return *this;
}
SymmMatrix<TYPE> &operator+=(const SymmMatrix<TYPE> &other) {
CHECK_INVARIANT(d_size == other.numRows(),
"Sizes don't match in the addition");
const TYPE *oData = other.getData();
TYPE *data = d_data.get();
for (unsigned int i = 0; i < d_dataSize; i++) {
data[i] += oData[i];
}
return *this;
}
SymmMatrix<TYPE> &operator-=(const SymmMatrix<TYPE> &other) {
CHECK_INVARIANT(d_size == other.numRows(),
"Sizes don't match in the addition");
const TYPE *oData = other.getData();
TYPE *data = d_data.get();
for (unsigned int i = 0; i < d_dataSize; i++) {
data[i] -= oData[i];
}
return *this;
}
//! in-place matrix multiplication
SymmMatrix<TYPE> &operator*=(const SymmMatrix<TYPE> &B) {
CHECK_INVARIANT(d_size == B.numRows(),
"Size mismatch during multiplication");
TYPE *cData = new TYPE[d_dataSize];
const TYPE *bData = B.getData();
TYPE *data = d_data.get();
for (unsigned int i = 0; i < d_size; i++) {
unsigned int idC = i * (i + 1) / 2;
for (unsigned int j = 0; j < i + 1; j++) {
unsigned int idCt = idC + j;
cData[idCt] = (TYPE)0.0;
for (unsigned int k = 0; k < d_size; k++) {
unsigned int idA, idB;
if (k <= i) {
idA = i * (i + 1) / 2 + k;
} else {
idA = k * (k + 1) / 2 + i;
}
if (k <= j) {
idB = j * (j + 1) / 2 + k;
} else {
idB = k * (k + 1) / 2 + j;
}
cData[idCt] += (data[idA] * bData[idB]);
}
}
}
for (unsigned int i = 0; i < d_dataSize; i++) {
data[i] = cData[i];
}
delete[] cData;
return (*this);
}
/* Transpose will basically return a copy of itself
*/
SymmMatrix<TYPE> &transpose(SymmMatrix<TYPE> &transpose) const {
CHECK_INVARIANT(d_size == transpose.numRows(),
"Size mismatch during transposing");
TYPE *tData = transpose.getData();
TYPE *data = d_data.get();
for (unsigned int i = 0; i < d_dataSize; i++) {
tData[i] = data[i];
}
return transpose;
}
SymmMatrix<TYPE> &transposeInplace() {
// nothing to be done we are symmetric
return (*this);
}
protected:
SymmMatrix() : d_data(0) {}
unsigned int d_size{0};
unsigned int d_dataSize{0};
DATA_SPTR d_data;
private:
SymmMatrix<TYPE> &operator=(const SymmMatrix<TYPE> &other);
};
//! SymmMatrix-SymmMatrix multiplication
/*!
Multiply SymmMatrix A with a second SymmMatrix B
and write the result to C = A*B
\param A the first SymmMatrix
\param B the second SymmMatrix to multiply
\param C SymmMatrix to use for the results
\return the results of multiplying A by B.
This is just a reference to C.
This method is reimplemented here for efficiency reasons
(we basically don't want to use getter and setter functions)
*/
template <class TYPE>
SymmMatrix<TYPE> &multiply(const SymmMatrix<TYPE> &A, const SymmMatrix<TYPE> &B,
SymmMatrix<TYPE> &C) {
unsigned int aSize = A.numRows();
CHECK_INVARIANT(B.numRows() == aSize,
"Size mismatch in matric multiplication");
CHECK_INVARIANT(C.numRows() == aSize,
"Size mismatch in matric multiplication");
TYPE *cData = C.getData();
const TYPE *aData = A.getData();
const TYPE *bData = B.getData();
for (unsigned int i = 0; i < aSize; i++) {
unsigned int idC = i * (i + 1) / 2;
for (unsigned int j = 0; j < i + 1; j++) {
unsigned int idCt = idC + j;
cData[idCt] = (TYPE)0.0;
for (unsigned int k = 0; k < aSize; k++) {
unsigned int idA, idB;
if (k <= i) {
idA = i * (i + 1) / 2 + k;
} else {
idA = k * (k + 1) / 2 + i;
}
if (k <= j) {
idB = j * (j + 1) / 2 + k;
} else {
idB = k * (k + 1) / 2 + j;
}
cData[idCt] += (aData[idA] * bData[idB]);
}
}
}
return C;
}
//! SymmMatrix-Vector multiplication
/*!
Multiply a SymmMatrix A with a Vector x
so the result is y = A*x
\param A the SymmMatrix for multiplication
\param x the Vector by which to multiply
\param y Vector to use for the results
\return the results of multiplying x by A
This is just a reference to y.
This method is reimplemented here for efficiency reasons
(we basically don't want to use getter and setter functions)
*/
template <class TYPE>
Vector<TYPE> &multiply(const SymmMatrix<TYPE> &A, const Vector<TYPE> &x,
Vector<TYPE> &y) {
unsigned int aSize = A.numRows();
CHECK_INVARIANT(aSize == x.size(), "Size mismatch during multiplication");
CHECK_INVARIANT(aSize == y.size(), "Size mismatch during multiplication");
const TYPE *xData = x.getData();
const TYPE *aData = A.getData();
TYPE *yData = y.getData();
for (unsigned int i = 0; i < aSize; i++) {
yData[i] = (TYPE)(0.0);
unsigned int idA = i * (i + 1) / 2;
for (unsigned int j = 0; j < i + 1; j++) {
// idA = i*(i+1)/2 + j;
yData[i] += (aData[idA] * xData[j]);
idA++;
}
idA--;
for (unsigned int j = i + 1; j < aSize; j++) {
// idA = j*(j+1)/2 + i;
idA += j;
yData[i] += (aData[idA] * xData[j]);
}
}
return y;
}
typedef SymmMatrix<double> DoubleSymmMatrix;
typedef SymmMatrix<int> IntSymmMatrix;
typedef SymmMatrix<unsigned int> UintSymmMatrix;
} // namespace RDNumeric
//! ostream operator for Matrix's
template <class TYPE>
std::ostream &operator<<(std::ostream &target,
const RDNumeric::SymmMatrix<TYPE> &mat) {
unsigned int nr = mat.numRows();
unsigned int nc = mat.numCols();
target << "Rows: " << mat.numRows() << " Columns: " << mat.numCols() << "\n";
for (unsigned int i = 0; i < nr; i++) {
for (unsigned int j = 0; j < nc; j++) {
target << std::setw(7) << std::setprecision(3) << mat.getVal(i, j);
}
target << "\n";
}
return target;
}
#endif
CHECK_INVARIANT
#define CHECK_INVARIANT(expr, mess)
Definition Invariant.h:100
URANGE_CHECK
#define URANGE_CHECK(x, hi)
Definition Invariant.h:141
Matrix.h
SquareMatrix.h
operator<<
std::ostream & operator<<(std::ostream &target, const RDNumeric::SymmMatrix< TYPE > &mat)
ostream operator for Matrix's
Definition SymmMatrix.h:346
RDNumeric::SymmMatrix
A symmetric matrix class.
Definition SymmMatrix.h:30
RDNumeric::SymmMatrix::SymmMatrix
SymmMatrix(unsigned int N)
Definition SymmMatrix.h:34
RDNumeric::SymmMatrix::transposeInplace
SymmMatrix< TYPE > & transposeInplace()
Definition SymmMatrix.h:230
RDNumeric::SymmMatrix::transpose
SymmMatrix< TYPE > & transpose(SymmMatrix< TYPE > &transpose) const
Definition SymmMatrix.h:219
RDNumeric::SymmMatrix::getData
const TYPE * getData() const
returns a const pointer to our data array
Definition SymmMatrix.h:141
RDNumeric::SymmMatrix< double >::d_data
DATA_SPTR d_data
Definition SymmMatrix.h:239
RDNumeric::SymmMatrix::SymmMatrix
SymmMatrix()
Definition SymmMatrix.h:236
RDNumeric::SymmMatrix< double >::d_size
unsigned int d_size
Definition SymmMatrix.h:237
RDNumeric::SymmMatrix::operator/=
SymmMatrix< TYPE > & operator/=(TYPE scale)
Definition SymmMatrix.h:151
RDNumeric::SymmMatrix::operator*=
SymmMatrix< TYPE > & operator*=(TYPE scale)
Definition SymmMatrix.h:143
RDNumeric::SymmMatrix::getCol
void getCol(unsigned int i, Vector< TYPE > &col)
Definition SymmMatrix.h:122
RDNumeric::SymmMatrix::setVal
void setVal(unsigned int i, unsigned int j, TYPE val)
Definition SymmMatrix.h:95
RDNumeric::SymmMatrix::getRow
void getRow(unsigned int i, Vector< TYPE > &row)
Definition SymmMatrix.h:107
RDNumeric::SymmMatrix::numCols
unsigned int numCols() const
returns the number of columns
Definition SymmMatrix.h:71
RDNumeric::SymmMatrix::SymmMatrix
SymmMatrix(unsigned int N, DATA_SPTR data)
Definition SymmMatrix.h:50
RDNumeric::SymmMatrix< double >::numRows
unsigned int numRows() const
Definition SymmMatrix.h:68
RDNumeric::SymmMatrix< double >::getDataSize
unsigned int getDataSize() const
Definition SymmMatrix.h:73
RDNumeric::SymmMatrix::getVal
TYPE getVal(unsigned int i, unsigned int j) const
Definition SymmMatrix.h:83
RDNumeric::SymmMatrix::operator+=
SymmMatrix< TYPE > & operator+=(const SymmMatrix< TYPE > &other)
Definition SymmMatrix.h:159
RDNumeric::SymmMatrix::setToIdentity
void setToIdentity()
Definition SymmMatrix.h:75
RDNumeric::SymmMatrix< double >::d_dataSize
unsigned int d_dataSize
Definition SymmMatrix.h:238
RDNumeric::SymmMatrix::operator*=
SymmMatrix< TYPE > & operator*=(const SymmMatrix< TYPE > &B)
in-place matrix multiplication
Definition SymmMatrix.h:182
RDNumeric::SymmMatrix::SymmMatrix
SymmMatrix(const SymmMatrix< TYPE > &other)
Definition SymmMatrix.h:55
RDNumeric::SymmMatrix::~SymmMatrix
~SymmMatrix()=default
RDNumeric::SymmMatrix::getData
TYPE * getData()
returns a pointer to our data array
Definition SymmMatrix.h:138
RDNumeric::SymmMatrix::operator-=
SymmMatrix< TYPE > & operator-=(const SymmMatrix< TYPE > &other)
Definition SymmMatrix.h:170
RDNumeric::SymmMatrix::DATA_SPTR
boost::shared_array< TYPE > DATA_SPTR
Definition SymmMatrix.h:32
RDNumeric::SymmMatrix::SymmMatrix
SymmMatrix(unsigned int N, TYPE val)
Definition SymmMatrix.h:40
RDNumeric::Vector
A class to represent vectors of numbers.
Definition Vector.h:30
RDNumeric::Vector::getData
constexpr TYPE * getData()
returns a pointer to our data array
Definition Vector.h:100
RDNumeric::Vector::size
constexpr unsigned int size() const
return the size (dimension) of the vector
Definition Vector.h:75
export.h
RDNumeric
Definition AlignPoints.h:18
RDNumeric::UintSymmMatrix
SymmMatrix< unsigned int > UintSymmMatrix
Definition SymmMatrix.h:341
RDNumeric::DoubleSymmMatrix
SymmMatrix< double > DoubleSymmMatrix
Definition SymmMatrix.h:339
RDNumeric::multiply
Matrix< TYPE > & multiply(const Matrix< TYPE > &A, const Matrix< TYPE > &B, Matrix< TYPE > &C)
Matrix multiplication.
Definition Matrix.h:272
RDNumeric::IntSymmMatrix
SymmMatrix< int > IntSymmMatrix
Definition SymmMatrix.h:340