pub_soft/rdkit
3
0
Fork 0
mirror of https://github.com/rdkit/rdkit.git synced 2026-06-03 21:44:30 +08:00
Code Issues Packages Projects Releases Wiki Activity

dos2unix and untabify... sigh

Browse Source
This commit is contained in:
Greg Landrum
2008-05-18 05:06:57 +00:00
parent de146f7e28
commit a643a9913a
12 changed files with 871 additions and 871 deletions
Download Patch File Download Diff File Expand all files Collapse all files

228
Code/DistGeom/BoundsMatrix.h
View File

@@ -1,114 +1,114 @@
//
// Copyright (C) 2004-2006 Rational Discovery LLC
// All Rights Reserved
//
// @@ All Rights Reserved @@
//
#ifndef __RD_BOUNDS_MATRIX_H__
#define __RD_BOUNDS_MATRIX_H__
#include <RDGeneral/Invariant.h>
#include <boost/smart_ptr.hpp>
#include <iostream>
#include <iomanip>
#include <Numerics/SquareMatrix.h>
namespace DistGeom {
//! Class to store the distance bound
/*!
Basically a N by N matrix
with lower distance bounds on the lower traingle and upper bounds in the upper
triangle
*/
class BoundsMatrix : public RDNumeric::SquareMatrix<double> {
public:
typedef boost::shared_array<double> DATA_SPTR;
explicit BoundsMatrix(unsigned int N) : RDNumeric::SquareMatrix<double>(N,0.0) {};
BoundsMatrix(unsigned int N, DATA_SPTR data) :
RDNumeric::SquareMatrix<double>(N,data) {};
//! Get the upper bound between points i and j
inline double getUpperBound(unsigned int i, unsigned int j) const {
RANGE_CHECK(0, i, d_nRows-1);
RANGE_CHECK(0, j, d_nCols-1);
if (i < j) {
return getVal(i,j);
} else {
return getVal(j,i);
}
}
//! Set the lower bound between points i and j
inline void setUpperBound(unsigned int i, unsigned int j, double val) {
RANGE_CHECK(0, i, d_nRows-1);
RANGE_CHECK(0, j, d_nCols-1);
CHECK_INVARIANT(val >= 0.0, "Negative upper bound");
if (i < j) {
setVal(i,j,val);
} else {
setVal(j,i,val);
}
}
//! Set the upper bound between points i and j only if it is better than
//! previously existing value (i.e. the new value is smaller)
inline void setUpperBoundIfBetter(unsigned int i, unsigned int j, double val) {
if ((val < getUpperBound(i, j)) && (val > getLowerBound(i, j)) ) {
setUpperBound(i, j, val);
}
}
//! Set the lower bound between points i and j
inline void setLowerBound(unsigned int i, unsigned int j, double val) {
RANGE_CHECK(0, i, d_nRows-1);
RANGE_CHECK(0, j, d_nCols-1);
CHECK_INVARIANT(val >= 0.0, "Negative lower bound");
if (i < j) {
setVal(j,i,val);
} else {
setVal(i,j,val);
}
}
//! Set the lower bound between points i and j only if it is better than
//! previously existing value (i.e. the new value is larger)
inline void setLowerBoundIfBetter(unsigned int i, unsigned int j, double val) {
if ((val > getLowerBound(i,j)) && (val < getUpperBound(i,j))) {
setLowerBound(i,j, val);
}
}
//! Get the lower bound between points i and j
inline double getLowerBound(unsigned int i, unsigned int j) const {
RANGE_CHECK(0, i, d_nRows-1);
RANGE_CHECK(0, j, d_nCols-1);
if (i < j) {
return getVal(j,i);
} else {
return getVal(i,j);
}
}
//! Do a simple check of the current bounds - i.e. all lower bounds are
//! smaller than the existing upper bounds
inline bool checkValid() const {
unsigned int i, j;
for (i = 1; i < d_nRows; i++) {
for (j = 0; j < i; j++) {
if (getUpperBound(i,j) < getLowerBound(i,j)) {
return false;
}
}
}
return true;
}
};
typedef boost::shared_ptr<BoundsMatrix> BoundsMatPtr;
}
#endif
//
// Copyright (C) 2004-2006 Rational Discovery LLC
// All Rights Reserved
//
// @@ All Rights Reserved @@
//
#ifndef __RD_BOUNDS_MATRIX_H__
#define __RD_BOUNDS_MATRIX_H__
#include <RDGeneral/Invariant.h>
#include <boost/smart_ptr.hpp>
#include <iostream>
#include <iomanip>
#include <Numerics/SquareMatrix.h>
namespace DistGeom {
//! Class to store the distance bound
/*!
Basically a N by N matrix
with lower distance bounds on the lower traingle and upper bounds in the upper
triangle
*/
class BoundsMatrix : public RDNumeric::SquareMatrix<double> {
public:
typedef boost::shared_array<double> DATA_SPTR;
explicit BoundsMatrix(unsigned int N) : RDNumeric::SquareMatrix<double>(N,0.0) {};
BoundsMatrix(unsigned int N, DATA_SPTR data) :
RDNumeric::SquareMatrix<double>(N,data) {};
//! Get the upper bound between points i and j
inline double getUpperBound(unsigned int i, unsigned int j) const {
RANGE_CHECK(0, i, d_nRows-1);
RANGE_CHECK(0, j, d_nCols-1);
if (i < j) {
return getVal(i,j);
} else {
return getVal(j,i);
}
}
//! Set the lower bound between points i and j
inline void setUpperBound(unsigned int i, unsigned int j, double val) {
RANGE_CHECK(0, i, d_nRows-1);
RANGE_CHECK(0, j, d_nCols-1);
CHECK_INVARIANT(val >= 0.0, "Negative upper bound");
if (i < j) {
setVal(i,j,val);
} else {
setVal(j,i,val);
}
}
//! Set the upper bound between points i and j only if it is better than
//! previously existing value (i.e. the new value is smaller)
inline void setUpperBoundIfBetter(unsigned int i, unsigned int j, double val) {
if ((val < getUpperBound(i, j)) && (val > getLowerBound(i, j)) ) {
setUpperBound(i, j, val);
}
}
//! Set the lower bound between points i and j
inline void setLowerBound(unsigned int i, unsigned int j, double val) {
RANGE_CHECK(0, i, d_nRows-1);
RANGE_CHECK(0, j, d_nCols-1);
CHECK_INVARIANT(val >= 0.0, "Negative lower bound");
if (i < j) {
setVal(j,i,val);
} else {
setVal(i,j,val);
}
}
//! Set the lower bound between points i and j only if it is better than
//! previously existing value (i.e. the new value is larger)
inline void setLowerBoundIfBetter(unsigned int i, unsigned int j, double val) {
if ((val > getLowerBound(i,j)) && (val < getUpperBound(i,j))) {
setLowerBound(i,j, val);
}
}
//! Get the lower bound between points i and j
inline double getLowerBound(unsigned int i, unsigned int j) const {
RANGE_CHECK(0, i, d_nRows-1);
RANGE_CHECK(0, j, d_nCols-1);
if (i < j) {
return getVal(j,i);
} else {
return getVal(i,j);
}
}
//! Do a simple check of the current bounds - i.e. all lower bounds are
//! smaller than the existing upper bounds
inline bool checkValid() const {
unsigned int i, j;
for (i = 1; i < d_nRows; i++) {
for (j = 0; j < i; j++) {
if (getUpperBound(i,j) < getLowerBound(i,j)) {
return false;
}
}
}
return true;
}
};
typedef boost::shared_ptr<BoundsMatrix> BoundsMatPtr;
}
#endif

38
Code/DistGeom/ChiralViolationContrib.cpp
View File

@@ -10,7 +10,7 @@
namespace DistGeom {
ChiralViolationContrib::ChiralViolationContrib(ForceFields::ForceField *owner, const ChiralSet* cset,
double weight) {
double weight) {
PRECONDITION(owner,"bad owner");
PRECONDITION(cset, "bad chiral set")
@@ -41,17 +41,17 @@ namespace DistGeom {
// even if we are minimizing in higher dimension the chiral volume is
// calculated using only the first 3 dimensions
RDGeom::Point3D v1(pos[d_idx1*dim] - pos[d_idx4*dim],
pos[d_idx1*dim+1] - pos[d_idx4*dim+1],
pos[d_idx1*dim+2] - pos[d_idx4*dim+2]);
pos[d_idx1*dim+1] - pos[d_idx4*dim+1],
pos[d_idx1*dim+2] - pos[d_idx4*dim+2]);
RDGeom::Point3D v2(pos[d_idx2*dim] - pos[d_idx4*dim],
pos[d_idx2*dim+1] - pos[d_idx4*dim+1],
pos[d_idx2*dim+2] - pos[d_idx4*dim+2]);
pos[d_idx2*dim+1] - pos[d_idx4*dim+1],
pos[d_idx2*dim+2] - pos[d_idx4*dim+2]);
RDGeom::Point3D v3(pos[d_idx3*dim] - pos[d_idx4*dim],
pos[d_idx3*dim+1] - pos[d_idx4*dim+1],
pos[d_idx3*dim+2] - pos[d_idx4*dim+2]);
pos[d_idx3*dim+1] - pos[d_idx4*dim+1],
pos[d_idx3*dim+2] - pos[d_idx4*dim+2]);
RDGeom::Point3D v2xv3 = v2.crossProduct(v3);
@@ -76,17 +76,17 @@ namespace DistGeom {
// even if we are minimizing in higher dimension the chiral volume is
// calculated using only the first 3 dimensions
RDGeom::Point3D v1(pos[d_idx1*dim] - pos[d_idx4*dim],
pos[d_idx1*dim+1] - pos[d_idx4*dim+1],
pos[d_idx1*dim+2] - pos[d_idx4*dim+2]);
pos[d_idx1*dim+1] - pos[d_idx4*dim+1],
pos[d_idx1*dim+2] - pos[d_idx4*dim+2]);
RDGeom::Point3D v2(pos[d_idx2*dim] - pos[d_idx4*dim],
pos[d_idx2*dim+1] - pos[d_idx4*dim+1],
pos[d_idx2*dim+2] - pos[d_idx4*dim+2]);
pos[d_idx2*dim+1] - pos[d_idx4*dim+1],
pos[d_idx2*dim+2] - pos[d_idx4*dim+2]);
RDGeom::Point3D v3(pos[d_idx3*dim] - pos[d_idx4*dim],
pos[d_idx3*dim+1] - pos[d_idx4*dim+1],
pos[d_idx3*dim+2] - pos[d_idx4*dim+2]);
pos[d_idx3*dim+1] - pos[d_idx4*dim+1],
pos[d_idx3*dim+2] - pos[d_idx4*dim+2]);
RDGeom::Point3D v2xv3 = v2.crossProduct(v3);
@@ -117,16 +117,16 @@ namespace DistGeom {
grad[dim*d_idx3 + 2] += preFactor*((v2.y)*(v1.x) - (v2.x)*(v1.y));
grad[dim*d_idx4] += preFactor*(pos[d_idx1*dim+2]*(pos[d_idx2*dim+1] - pos[d_idx3*dim+1])
+ pos[d_idx2*dim+2]*(pos[d_idx3*dim+1] - pos[d_idx1*dim+1])
+ pos[d_idx3*dim+2]*(pos[d_idx1*dim+1] - pos[d_idx2*dim+1]));
+ pos[d_idx2*dim+2]*(pos[d_idx3*dim+1] - pos[d_idx1*dim+1])
+ pos[d_idx3*dim+2]*(pos[d_idx1*dim+1] - pos[d_idx2*dim+1]));
grad[dim*d_idx4+1] += preFactor*(pos[d_idx1*dim]*(pos[d_idx2*dim+2] - pos[d_idx3*dim+2])
+ pos[d_idx2*dim]*(pos[d_idx3*dim+2] - pos[d_idx1*dim+2])
+ pos[d_idx3*dim]*(pos[d_idx1*dim+2] - pos[d_idx2*dim+2]));
+ pos[d_idx2*dim]*(pos[d_idx3*dim+2] - pos[d_idx1*dim+2])
+ pos[d_idx3*dim]*(pos[d_idx1*dim+2] - pos[d_idx2*dim+2]));
grad[dim*d_idx4+2] += preFactor*(pos[d_idx1*dim+1]*(pos[d_idx2*dim] - pos[d_idx3*dim])
+ pos[d_idx2*dim+1]*(pos[d_idx3*dim] - pos[d_idx1*dim])
+ pos[d_idx3*dim+1]*(pos[d_idx1*dim] - pos[d_idx2*dim]));
+ pos[d_idx2*dim+1]*(pos[d_idx3*dim] - pos[d_idx1*dim])
+ pos[d_idx3*dim+1]*(pos[d_idx1*dim] - pos[d_idx2*dim]));
}
}

180
Code/DistGeom/DistGeomUtils.h
View File

@@ -1,90 +1,90 @@
//
// Copyright (C) 2004-2007 Greg Landrum and Rational Discovery LLC
//
// @@ All Rights Reserved @@
//
#ifndef _RD_DISTGEOMUTILS_H_
#define _RD_DISTGEOMUTILS_H_
#include "BoundsMatrix.h"
#include <Numerics/SymmMatrix.h>
#include <map>
#include <Geometry/point.h>
#include "EmbedObject.h"
namespace ForceFields {
class ForceField;
}
namespace DistGeom {
//! Pick a distance matrix at random such that the
//! distance satisfy the bounds in the BoundsMatrix
/*!
\param mmat Bounds matrix
\param distmat Storage for randomly chosen distances
\param seed Optional value to seed the random number generator
\return the largest element of the distance matrix
*/
double pickRandomDistMat(const BoundsMatrix &mmat,
RDNumeric::SymmMatrix<double> &distmat, int seed=-1);
//! Compute an initial embedded in 3D based on a distance matrix
/*!
This function follows the embed algorithm mentioned in
"Distance Geometry and Molecular Conformation" by G.M.Crippen and T.F.Havel
(pages 312-313)
\param distmat Distance matrix
\param positions A vector of pointers to Points to write out the resulting coordinates
\param randNegEig If set to true and if any of the eigen values are negative, we will
pick the corresponding components of the coordinates at random
\param numZeroFail Fail embedding is more this many (or more) eigen values are zero
\return true if the embedding was successful
*/
bool computeInitialCoords(const RDNumeric::SymmMatrix<double> &distmat,
RDGeom::PointPtrVect &positions, bool randNegEig=false,
unsigned int numZeroFail=2);
//! places atoms randomly in a box
/*!
\param positions A vector of pointers to Points to write out the resulting coordinates
\param boxSize If set to true and if any of the eigen values are negative, we will
pick the corresponding components of the coordinates at random
\return true if the embedding was successful
*/
bool computeRandomCoords(RDGeom::PointPtrVect &positions, double boxSize);
//! Setup the error function for violation of distance bounds as a forcefield
/*!
This is based on function E3 on page 311 of "Distance Geometry in Molecular
Modeling" Jeffrey M.Blaney and J.Scott Dixon, Review in Computational Chemistry,
Volume V
\param mmat Distance bounds matrix
\param positions A vector of pointers to Points to write out the resulting coordinates
\param csets The vector of chiral points (type: ChiralSet)
\param weightChiral weight to be used to enforce chirality
\param weightFourthDim another chiral weight
\param extraWeights an optional set of weights for distance bounds violations
\param basinSizeTol Optional: any distance bound with a basin (distance between max and
min bounds) larger than this value will not be included in the force
field used to cleanup the structure.
\return a pointer to a ForceField suitable for cleaning up the violations.
<b>NOTE:</b> the caller is responsible for deleting this force field.
*/
ForceFields::ForceField *constructForceField(const BoundsMatrix &mmat,
RDGeom::PointPtrVect &positions, const VECT_CHIRALSET &csets,
double weightChiral=1.0,
double weightFourthDim=0.1,
std::map< std::pair<int,int>,double> *extraWeights=0,
double basinSizeTol=5.0);
}
#endif
//
// Copyright (C) 2004-2007 Greg Landrum and Rational Discovery LLC
//
// @@ All Rights Reserved @@
//
#ifndef _RD_DISTGEOMUTILS_H_
#define _RD_DISTGEOMUTILS_H_
#include "BoundsMatrix.h"
#include <Numerics/SymmMatrix.h>
#include <map>
#include <Geometry/point.h>
#include "EmbedObject.h"
namespace ForceFields {
class ForceField;
}
namespace DistGeom {
//! Pick a distance matrix at random such that the
//! distance satisfy the bounds in the BoundsMatrix
/*!
\param mmat Bounds matrix
\param distmat Storage for randomly chosen distances
\param seed Optional value to seed the random number generator
\return the largest element of the distance matrix
*/
double pickRandomDistMat(const BoundsMatrix &mmat,
RDNumeric::SymmMatrix<double> &distmat, int seed=-1);
//! Compute an initial embedded in 3D based on a distance matrix
/*!
This function follows the embed algorithm mentioned in
"Distance Geometry and Molecular Conformation" by G.M.Crippen and T.F.Havel
(pages 312-313)
\param distmat Distance matrix
\param positions A vector of pointers to Points to write out the resulting coordinates
\param randNegEig If set to true and if any of the eigen values are negative, we will
pick the corresponding components of the coordinates at random
\param numZeroFail Fail embedding is more this many (or more) eigen values are zero
\return true if the embedding was successful
*/
bool computeInitialCoords(const RDNumeric::SymmMatrix<double> &distmat,
RDGeom::PointPtrVect &positions, bool randNegEig=false,
unsigned int numZeroFail=2);
//! places atoms randomly in a box
/*!
\param positions A vector of pointers to Points to write out the resulting coordinates
\param boxSize If set to true and if any of the eigen values are negative, we will
pick the corresponding components of the coordinates at random
\return true if the embedding was successful
*/
bool computeRandomCoords(RDGeom::PointPtrVect &positions, double boxSize);
//! Setup the error function for violation of distance bounds as a forcefield
/*!
This is based on function E3 on page 311 of "Distance Geometry in Molecular
Modeling" Jeffrey M.Blaney and J.Scott Dixon, Review in Computational Chemistry,
Volume V
\param mmat Distance bounds matrix
\param positions A vector of pointers to Points to write out the resulting coordinates
\param csets The vector of chiral points (type: ChiralSet)
\param weightChiral weight to be used to enforce chirality
\param weightFourthDim another chiral weight
\param extraWeights an optional set of weights for distance bounds violations
\param basinSizeTol Optional: any distance bound with a basin (distance between max and
min bounds) larger than this value will not be included in the force
field used to cleanup the structure.
\return a pointer to a ForceField suitable for cleaning up the violations.
<b>NOTE:</b> the caller is responsible for deleting this force field.
*/
ForceFields::ForceField *constructForceField(const BoundsMatrix &mmat,
RDGeom::PointPtrVect &positions, const VECT_CHIRALSET &csets,
double weightChiral=1.0,
double weightFourthDim=0.1,
std::map< std::pair<int,int>,double> *extraWeights=0,
double basinSizeTol=5.0);
}
#endif

164
Code/DistGeom/DistViolationContrib.cpp
View File

@@ -1,82 +1,82 @@
// $Id$
//
// Copyright (C) 2004-2006 Rational Discovery LLC
//
// @@ All Rights Reserved @@
//
#include "DistViolationContrib.h"
#include <ForceField/ForceField.h>
#include <RDGeneral/Invariant.h>
namespace DistGeom {
DistViolationContrib::DistViolationContrib(ForceFields::ForceField *owner,unsigned int idx1,
unsigned int idx2,
double ub, double lb, double weight) {
PRECONDITION(owner,"bad owner");
RANGE_CHECK(0,idx1,owner->positions().size()-1);
RANGE_CHECK(0,idx2,owner->positions().size()-1);
dp_forceField = owner;
d_end1Idx = idx1;
d_end2Idx = idx2;
d_ub = ub;
d_lb = lb;
d_weight = weight;
}
double DistViolationContrib::getEnergy(double *pos) const {
PRECONDITION(dp_forceField,"no owner");
PRECONDITION(pos,"bad vector");
double d = this->dp_forceField->distance(this->d_end1Idx,this->d_end2Idx,pos);
double val=0.0;
if (d > d_ub) {
val = ((d*d)/(d_ub*d_ub)) - 1.0;
} else if (d < d_lb) {
val = ((2*d_lb*d_lb)/(d_lb*d_lb + d*d)) - 1.0;
}
double res;
if(val>0.0){
res=d_weight*val*val;
}else{
res=0;
}
return res;
}
void DistViolationContrib::getGrad(double *pos, double *grad) const {
PRECONDITION(dp_forceField,"no owner");
PRECONDITION(pos,"bad vector");
PRECONDITION(grad,"bad vector");
unsigned int dim=this->dp_forceField->dimension();
double d = this->dp_forceField->distance(this->d_end1Idx,this->d_end2Idx,pos);
double preFactor = 0.0;
if (d > d_ub) {
double u2 = d_ub*d_ub;
preFactor = 4.*(((d*d)/u2) - 1.0)*(d/u2);
} else if (d < d_lb) {
double d2 = d*d;
double l2 = d_lb*d_lb;
preFactor = -8.*((l2-d2)/pow(l2+d2,3))*d*l2;
} else {
return;
}
double *end1Coords = &(pos[dim*this->d_end1Idx]);
double *end2Coords = &(pos[dim*this->d_end2Idx]);
for(unsigned int i=0;i<dim;i++){
double dGrad;
if(d>0.0){
dGrad= d_weight*preFactor * (end1Coords[i]-end2Coords[i])/d;
} else {
// FIX: this likely isn't right
dGrad= d_weight*preFactor * (end1Coords[i]-end2Coords[i]);
}
grad[dim*this->d_end1Idx+i] += dGrad;
grad[dim*this->d_end2Idx+i] -= dGrad;
}
}
}
// $Id$
//
// Copyright (C) 2004-2006 Rational Discovery LLC
//
// @@ All Rights Reserved @@
//
#include "DistViolationContrib.h"
#include <ForceField/ForceField.h>
#include <RDGeneral/Invariant.h>
namespace DistGeom {
DistViolationContrib::DistViolationContrib(ForceFields::ForceField *owner,unsigned int idx1,
unsigned int idx2,
double ub, double lb, double weight) {
PRECONDITION(owner,"bad owner");
RANGE_CHECK(0,idx1,owner->positions().size()-1);
RANGE_CHECK(0,idx2,owner->positions().size()-1);
dp_forceField = owner;
d_end1Idx = idx1;
d_end2Idx = idx2;
d_ub = ub;
d_lb = lb;
d_weight = weight;
}
double DistViolationContrib::getEnergy(double *pos) const {
PRECONDITION(dp_forceField,"no owner");
PRECONDITION(pos,"bad vector");
double d = this->dp_forceField->distance(this->d_end1Idx,this->d_end2Idx,pos);
double val=0.0;
if (d > d_ub) {
val = ((d*d)/(d_ub*d_ub)) - 1.0;
} else if (d < d_lb) {
val = ((2*d_lb*d_lb)/(d_lb*d_lb + d*d)) - 1.0;
}
double res;
if(val>0.0){
res=d_weight*val*val;
}else{
res=0;
}
return res;
}
void DistViolationContrib::getGrad(double *pos, double *grad) const {
PRECONDITION(dp_forceField,"no owner");
PRECONDITION(pos,"bad vector");
PRECONDITION(grad,"bad vector");
unsigned int dim=this->dp_forceField->dimension();
double d = this->dp_forceField->distance(this->d_end1Idx,this->d_end2Idx,pos);
double preFactor = 0.0;
if (d > d_ub) {
double u2 = d_ub*d_ub;
preFactor = 4.*(((d*d)/u2) - 1.0)*(d/u2);
} else if (d < d_lb) {
double d2 = d*d;
double l2 = d_lb*d_lb;
preFactor = -8.*((l2-d2)/pow(l2+d2,3))*d*l2;
} else {
return;
}
double *end1Coords = &(pos[dim*this->d_end1Idx]);
double *end2Coords = &(pos[dim*this->d_end2Idx]);
for(unsigned int i=0;i<dim;i++){
double dGrad;
if(d>0.0){
dGrad= d_weight*preFactor * (end1Coords[i]-end2Coords[i])/d;
} else {
// FIX: this likely isn't right
dGrad= d_weight*preFactor * (end1Coords[i]-end2Coords[i]);
}
grad[dim*this->d_end1Idx+i] += dGrad;
grad[dim*this->d_end2Idx+i] -= dGrad;
}
}
}

84
Code/DistGeom/DistViolationContrib.h
View File

@@ -1,42 +1,42 @@
//
// Copyright (C) 2004-2006 Rational Discovery LLC
//
// @@ All Rights Reserved @@
//
#ifndef __RD_DISTVIOLATIONCONTRIB_H__
#define __RD_DISTVIOLATIONCONTRIB_H__
#include <ForceField/Contrib.h>
namespace DistGeom {
//! A term to capture the violation of the upper and lower bounds by
//! distance between two points
class DistViolationContrib : public ForceFields::ForceFieldContrib {
public:
DistViolationContrib() : d_end1Idx(0), d_end2Idx(0), d_ub(1000.0), d_lb(0.0), d_weight(1.0) {};
//! Constructor
/*!
\param owner pointer to the owning ForceField
\param idx1 index of end1 in the ForceField's positions
\param idx2 index of end2 in the ForceField's positions
\param ub Upper bound on the distance
\param lb Lower bound on the distance
\param weight optional weight for this contribution
*/
DistViolationContrib(ForceFields::ForceField *owner,unsigned int idx1,unsigned int idx2,
double ub, double lb, double weight=1.0);
double getEnergy(double *pos) const;
void getGrad(double *pos, double *grad) const;
private:
unsigned int d_end1Idx,d_end2Idx; //!< indices of end points
double d_ub; //!< upper bound on the distance between d_end1Idx,d_end2Idx
double d_lb; //!< lower bound on the distance between d_end1Idx,d_end2Idx
double d_weight; //!< used to adjust relative contribution weights
};
}
#endif
//
// Copyright (C) 2004-2006 Rational Discovery LLC
//
// @@ All Rights Reserved @@
//
#ifndef __RD_DISTVIOLATIONCONTRIB_H__
#define __RD_DISTVIOLATIONCONTRIB_H__
#include <ForceField/Contrib.h>
namespace DistGeom {
//! A term to capture the violation of the upper and lower bounds by
//! distance between two points
class DistViolationContrib : public ForceFields::ForceFieldContrib {
public:
DistViolationContrib() : d_end1Idx(0), d_end2Idx(0), d_ub(1000.0), d_lb(0.0), d_weight(1.0) {};
//! Constructor
/*!
\param owner pointer to the owning ForceField
\param idx1 index of end1 in the ForceField's positions
\param idx2 index of end2 in the ForceField's positions
\param ub Upper bound on the distance
\param lb Lower bound on the distance
\param weight optional weight for this contribution
*/
DistViolationContrib(ForceFields::ForceField *owner,unsigned int idx1,unsigned int idx2,
double ub, double lb, double weight=1.0);
double getEnergy(double *pos) const;
void getGrad(double *pos, double *grad) const;
private:
unsigned int d_end1Idx,d_end2Idx; //!< indices of end points
double d_ub; //!< upper bound on the distance between d_end1Idx,d_end2Idx
double d_lb; //!< lower bound on the distance between d_end1Idx,d_end2Idx
double d_weight; //!< used to adjust relative contribution weights
};
}
#endif

202
Code/DistGeom/EmbedObject.h
View File

@@ -1,101 +1,101 @@
//
// Copyright (C) 2004-2006 Rational Discovery LLC
//
// @@ All Rights Reserved @@
//
#ifndef __RD_EMBED_OBJECT_H__
#define __RD_EMBED_OBJECT_H__
#include <RDGeneral/Invariant.h>
#include <RDGeneral/types.h>
#include <boost/smart_ptr.hpp>
namespace DistGeom {
/*! \brief Class used to store a quartet of points and chiral volume bounds on them
*
*/
class ChiralSet {
public:
unsigned int d_idx1;
unsigned int d_idx2;
unsigned int d_idx3;
unsigned int d_idx4;
double d_volumeLowerBound;
double d_volumeUpperBound;
ChiralSet(unsigned int pid1, unsigned int pid2, unsigned int pid3, unsigned int pid4,
double lowerVolBound, double upperVolBound) {
CHECK_INVARIANT(lowerVolBound <= upperVolBound, "Inconsistent bounds\n");
d_idx1 = pid1;
d_idx2 = pid2;
d_idx3 = pid3;
d_idx4 = pid4;
d_volumeLowerBound = lowerVolBound;
d_volumeUpperBound = upperVolBound;
}
inline double getUpperVolumeBound() const {
return d_volumeUpperBound;
}
inline double getLowerVolumeBound() const {
return d_volumeLowerBound;
}
};
typedef boost::shared_ptr<ChiralSet> ChiralSetPtr;
typedef std::vector<ChiralSetPtr> VECT_CHIRALSET;
/*! \brief Class that contains a an object that needs to be embedded.
*
* Basially a pointer to the MetricMatrix matric consisting of the upper and
* lower bounds on the distances, and a list of quartets of points to which chirality
* constraints need to be applied.
*/
/*
class EmbedObject {
public:
EmbedObject() {};
inline void setMetricMatrix(MetricMatrix *metricMat) {
PRE_CONDITION(metricMat, "");
d_metricMat = MetricMatPtr(metricMat);
}
inline void setMetricMatrix(MetricMatPtr metricMatPtr) {
d_metricMat = metricMatPtr;
}
inline const MetricMatPtr getMetricMatrix() const {
return d_metricMat;
}
inline MetricMatPtr getMetricMatrix() {
return d_metricMat;
}
inline void addChiralSet(ChiralSet *cset) {
PRE_CONDITION(cset, "");
d_chiralSets.push_back(ChiralSetPtr(cset));
}
inline void addChiralSet(ChiralSetPtr cset) {
d_chiralSets.push_back(cset);
}
const VECT_CHIRALSET &getChiralSets() const {
return d_chiralSets;
}
private:
MetricMatPtr d_metricMat;
VECT_CHIRALSET d_chiralSets;
};*/
}
#endif
//
// Copyright (C) 2004-2006 Rational Discovery LLC
//
// @@ All Rights Reserved @@
//
#ifndef __RD_EMBED_OBJECT_H__
#define __RD_EMBED_OBJECT_H__
#include <RDGeneral/Invariant.h>
#include <RDGeneral/types.h>
#include <boost/smart_ptr.hpp>
namespace DistGeom {
/*! \brief Class used to store a quartet of points and chiral volume bounds on them
*
*/
class ChiralSet {
public:
unsigned int d_idx1;
unsigned int d_idx2;
unsigned int d_idx3;
unsigned int d_idx4;
double d_volumeLowerBound;
double d_volumeUpperBound;
ChiralSet(unsigned int pid1, unsigned int pid2, unsigned int pid3, unsigned int pid4,
double lowerVolBound, double upperVolBound) {
CHECK_INVARIANT(lowerVolBound <= upperVolBound, "Inconsistent bounds\n");
d_idx1 = pid1;
d_idx2 = pid2;
d_idx3 = pid3;
d_idx4 = pid4;
d_volumeLowerBound = lowerVolBound;
d_volumeUpperBound = upperVolBound;
}
inline double getUpperVolumeBound() const {
return d_volumeUpperBound;
}
inline double getLowerVolumeBound() const {
return d_volumeLowerBound;
}
};
typedef boost::shared_ptr<ChiralSet> ChiralSetPtr;
typedef std::vector<ChiralSetPtr> VECT_CHIRALSET;
/*! \brief Class that contains a an object that needs to be embedded.
*
* Basially a pointer to the MetricMatrix matric consisting of the upper and
* lower bounds on the distances, and a list of quartets of points to which chirality
* constraints need to be applied.
*/
/*
class EmbedObject {
public:
EmbedObject() {};
inline void setMetricMatrix(MetricMatrix *metricMat) {
PRE_CONDITION(metricMat, "");
d_metricMat = MetricMatPtr(metricMat);
}
inline void setMetricMatrix(MetricMatPtr metricMatPtr) {
d_metricMat = metricMatPtr;
}
inline const MetricMatPtr getMetricMatrix() const {
return d_metricMat;
}
inline MetricMatPtr getMetricMatrix() {
return d_metricMat;
}
inline void addChiralSet(ChiralSet *cset) {
PRE_CONDITION(cset, "");
d_chiralSets.push_back(ChiralSetPtr(cset));
}
inline void addChiralSet(ChiralSetPtr cset) {
d_chiralSets.push_back(cset);
}
const VECT_CHIRALSET &getChiralSets() const {
return d_chiralSets;
}
private:
MetricMatPtr d_metricMat;
VECT_CHIRALSET d_chiralSets;
};*/
}
#endif

2
Code/DistGeom/FourthDimContrib.h
View File

@@ -22,7 +22,7 @@ namespace DistGeom {
*/
FourthDimContrib(ForceFields::ForceField *owner,unsigned int idx,
double weight) : d_idx(idx), d_weight(weight) {
double weight) : d_idx(idx), d_weight(weight) {
PRECONDITION(owner,"bad force field");
PRECONDITION(owner->dimension()==4,"force field has wrong dimension");
dp_forceField=owner;

108
Code/DistGeom/TriangleSmooth.cpp
View File

@@ -1,54 +1,54 @@
// $Id$
//
// Copyright (C) 2004-2006 Rational Discovery LLC
//
// @@ All Rights Reserved @@
//
#include "BoundsMatrix.h"
#include "TriangleSmooth.h"
namespace DistGeom {
bool triangleSmoothBounds(BoundsMatPtr boundsMat) {
return triangleSmoothBounds(boundsMat.get());
}
bool triangleSmoothBounds(BoundsMatrix *boundsMat) {
int npt = boundsMat->numRows();
int i, j, k;
double Uik, Lik, Ukj, sumUikUkj, diffLikUjk, diffLjkUik;
for (k = 0; k < npt; k++) {
for (i = 0; i < npt-1; i++) {
if (i == k) {
continue;
}
Uik = boundsMat->getUpperBound(i,k);
Lik = boundsMat->getLowerBound(i,k);
for (j = i+1; j < npt; j++) {
if (j == k) {
continue;
}
Ukj = boundsMat->getUpperBound(k,j);
sumUikUkj = Uik + Ukj;
if (boundsMat->getUpperBound(i,j) > sumUikUkj) {
boundsMat->setUpperBound(i,j, sumUikUkj);
}
diffLikUjk = Lik - Ukj;
diffLjkUik = boundsMat->getLowerBound(j,k) - Uik;
if (boundsMat->getLowerBound(i,j) < diffLikUjk) {
boundsMat->setLowerBound(i,j, diffLikUjk);
} else if (boundsMat->getLowerBound(i,j) < diffLjkUik) {
boundsMat->setLowerBound(i,j, diffLjkUik);
}
if (boundsMat->getLowerBound(i,j) > boundsMat->getUpperBound(i,j)) {
//std::cout << boundsMat->getLowerBound(i,j) << " " << boundsMat->getUpperBound(i,j) << "\n";
return false;
}
}
}
}
return true;
}
}
// $Id$
//
// Copyright (C) 2004-2006 Rational Discovery LLC
//
// @@ All Rights Reserved @@
//
#include "BoundsMatrix.h"
#include "TriangleSmooth.h"
namespace DistGeom {
bool triangleSmoothBounds(BoundsMatPtr boundsMat) {
return triangleSmoothBounds(boundsMat.get());
}
bool triangleSmoothBounds(BoundsMatrix *boundsMat) {
int npt = boundsMat->numRows();
int i, j, k;
double Uik, Lik, Ukj, sumUikUkj, diffLikUjk, diffLjkUik;
for (k = 0; k < npt; k++) {
for (i = 0; i < npt-1; i++) {
if (i == k) {
continue;
}
Uik = boundsMat->getUpperBound(i,k);
Lik = boundsMat->getLowerBound(i,k);
for (j = i+1; j < npt; j++) {
if (j == k) {
continue;
}
Ukj = boundsMat->getUpperBound(k,j);
sumUikUkj = Uik + Ukj;
if (boundsMat->getUpperBound(i,j) > sumUikUkj) {
boundsMat->setUpperBound(i,j, sumUikUkj);
}
diffLikUjk = Lik - Ukj;
diffLjkUik = boundsMat->getLowerBound(j,k) - Uik;
if (boundsMat->getLowerBound(i,j) < diffLikUjk) {
boundsMat->setLowerBound(i,j, diffLikUjk);
} else if (boundsMat->getLowerBound(i,j) < diffLjkUik) {
boundsMat->setLowerBound(i,j, diffLjkUik);
}
if (boundsMat->getLowerBound(i,j) > boundsMat->getUpperBound(i,j)) {
//std::cout << boundsMat->getLowerBound(i,j) << " " << boundsMat->getUpperBound(i,j) << "\n";
return false;
}
}
}
}
return true;
}
}

60
Code/DistGeom/TriangleSmooth.h
View File

@@ -1,30 +1,30 @@
//
// Copyright (C) 2004-2006 Rational Discovery LLC
//
// @@ All Rights Reserved @@
//
#ifndef __RD_TRIANGLE_SMOOTH_H__
#define __RD_TRIANGLE_SMOOTH_H__
#include "BoundsMatrix.h"
namespace DistGeom {
//! Smooth the upper and lower bound in a metric matrix so that triangle
//! inequality is not violated
/*!
This an implementation of the O(N^3) algorithm given on pages 252-253 of
"Distance Geometry and Molecular Conformation" by G.M.Crippen and T.F.Havel
Research Studies Press, 1988. There are other (slightly) more implementations
(see pages 301-302 in the above book), but that is for later
\param boundsMat A pointer to the distance bounds matrix
*/
bool triangleSmoothBounds(BoundsMatrix *boundsMat);
//! \overload
bool triangleSmoothBounds(BoundsMatPtr boundsMat);
}
#endif
//
// Copyright (C) 2004-2006 Rational Discovery LLC
//
// @@ All Rights Reserved @@
//
#ifndef __RD_TRIANGLE_SMOOTH_H__
#define __RD_TRIANGLE_SMOOTH_H__
#include "BoundsMatrix.h"
namespace DistGeom {
//! Smooth the upper and lower bound in a metric matrix so that triangle
//! inequality is not violated
/*!
This an implementation of the O(N^3) algorithm given on pages 252-253 of
"Distance Geometry and Molecular Conformation" by G.M.Crippen and T.F.Havel
Research Studies Press, 1988. There are other (slightly) more implementations
(see pages 301-302 in the above book), but that is for later
\param boundsMat A pointer to the distance bounds matrix
*/
bool triangleSmoothBounds(BoundsMatrix *boundsMat);
//! \overload
bool triangleSmoothBounds(BoundsMatPtr boundsMat);
}
#endif

432
Code/DistGeom/Wrap/DistGeom.cpp
View File

@@ -1,216 +1,216 @@
// $Id$
//
// Copyright (C) 2004-2006 Rational Discovery LLC
//
// @@ All Rights Reserved @@
//
#include <boost/python.hpp>
#define PY_ARRAY_UNIQUE_SYMBOL DistGeom_array_API
#include "Numeric/arrayobject.h"
#include <RDBoost/Wrap.h>
#include <Geometry/point.h>
#include <Numerics/Matrix.h>
#include <Numerics/SymmMatrix.h>
#include <DistGeom/BoundsMatrix.h>
#include <DistGeom/TriangleSmooth.h>
#include <DistGeom/DistGeomUtils.h>
#include <DistGeom/EmbedObject.h>
#include <ForceField/ForceField.h>
#include <vector>
#include <map>
namespace python = boost::python;
namespace RDKit {
bool doTriangleSmoothing(python::object boundsMatArg){
PyObject *boundsMatObj = boundsMatArg.ptr();
if(!PyArray_Check(boundsMatObj))
throw_value_error("Argument isn't an array");
PyArrayObject *boundsMat=reinterpret_cast<PyArrayObject *>(boundsMatObj);
// get the dimensions of the array
int nrows = boundsMat->dimensions[0];
int ncols = boundsMat->dimensions[1];
if(nrows!=ncols)
throw_value_error("The array has to be square");
if(nrows<=0)
throw_value_error("The array has to have a nonzero size");
if (boundsMat->descr->type_num != PyArray_DOUBLE)
throw_value_error("Only double arrays are currently supported");
unsigned int dSize = nrows*nrows;
double *cData = new double[dSize];
double *inData = reinterpret_cast<double *>(boundsMat->data);
memcpy(static_cast<void *>(cData),
static_cast<const void *>(inData),
dSize*sizeof(double));
DistGeom::BoundsMatrix::DATA_SPTR sdata(cData);
DistGeom::BoundsMatrix bm(nrows,sdata);
bool res=DistGeom::triangleSmoothBounds(&bm);
memcpy(static_cast<void *>(inData),
static_cast<const void *>(cData),
dSize*sizeof(double));
return (PyObject *) res;
}
PyObject *embedBoundsMatrix(python::object boundsMatArg,int maxIters=10,
bool randomizeOnFailure=false,int numZeroFail=2,
python::list weights=python::list(),
int randomSeed=-1){
PyObject *boundsMatObj = boundsMatArg.ptr();
if(!PyArray_Check(boundsMatObj))
throw_value_error("Argument isn't an array");
PyArrayObject *boundsMat=reinterpret_cast<PyArrayObject *>(boundsMatObj);
// get the dimensions of the array
unsigned int nrows = boundsMat->dimensions[0];
unsigned int ncols = boundsMat->dimensions[1];
if(nrows!=ncols)
throw_value_error("The array has to be square");
if(nrows<=0)
throw_value_error("The array has to have a nonzero size");
if (boundsMat->descr->type_num != PyArray_DOUBLE)
throw_value_error("Only double arrays are currently supported");
unsigned int dSize = nrows*nrows;
double *cData = new double[dSize];
double *inData = reinterpret_cast<double *>(boundsMat->data);
memcpy(static_cast<void *>(cData),
static_cast<const void *>(inData),
dSize*sizeof(double));
DistGeom::BoundsMatrix::DATA_SPTR sdata(cData);
DistGeom::BoundsMatrix bm(nrows,sdata);
RDGeom::Point3D *positions=new RDGeom::Point3D[nrows];
std::vector<RDGeom::Point *> posPtrs;
for (unsigned int i = 0; i < nrows; i++) {
posPtrs.push_back(&positions[i]);
}
RDNumeric::DoubleSymmMatrix distMat(nrows, 0.0);
// ---- ---- ---- ---- ---- ---- ---- ---- ----
// start the embedding:
bool gotCoords=false;
for(int iter=0;iter<maxIters && !gotCoords;iter++){
// pick a random distance matrix
DistGeom::pickRandomDistMat(bm,distMat,randomSeed);
// update the seed:
if(randomSeed>=0) randomSeed+=iter*999;
//std::cerr << "Random matrix:" << std::endl;
//std::cerr << distMat;
// and embed it:
gotCoords=DistGeom::computeInitialCoords(distMat,posPtrs,randomizeOnFailure,
numZeroFail);
}
if(gotCoords){
std::map<std::pair<int,int>,double> weightMap;
unsigned int nElems=PySequence_Size(weights.ptr());
for(unsigned int entryIdx=0;entryIdx<nElems;entryIdx++){
PyObject *entry=PySequence_GetItem(weights.ptr(),entryIdx);
if(!PySequence_Check(entry) || PySequence_Size(entry)!=3){
throw_value_error("weights argument must be a sequence of 3-sequences");
}
int idx1=PyInt_AsLong(PySequence_GetItem(entry,0));
int idx2=PyInt_AsLong(PySequence_GetItem(entry,1));
double w=PyFloat_AsDouble(PySequence_GetItem(entry,2));
weightMap[std::make_pair<int,int>(idx1,idx2)]=w;
}
DistGeom::VECT_CHIRALSET csets;
ForceFields::ForceField *field = DistGeom::constructForceField(bm,posPtrs,csets,0.0, 0.0,
&weightMap);
CHECK_INVARIANT(field,"could not build dgeom force field");
field->initialize();
if(field->calcEnergy()>1e-5){
int needMore=1;
while(needMore){
needMore=field->minimize();
}
}
delete field;
} else {
throw_value_error("could not embed matrix");
}
// ---- ---- ---- ---- ---- ---- ---- ---- ----
// construct the results matrix:
int dims[2];
dims[0] = nrows;
dims[1] = 3;
PyArrayObject *res = (PyArrayObject *)PyArray_FromDims(2,dims,PyArray_DOUBLE);
double *resData=reinterpret_cast<double *>(res->data);
for(unsigned int i=0;i<nrows;i++){
unsigned int iTab=i*3;
for (unsigned int j = 0; j < 3; ++j) {
resData[iTab + j]=positions[i][j]; //.x;
//resData[iTab+1]=positions[i].y;
//resData[iTab+2]=positions[i].z;
}
}
delete [] positions;
return (PyObject *) res;
}
}
BOOST_PYTHON_MODULE(DistGeom) {
python::scope().attr("__doc__") =
"Module containing functions for basic distance geometry operations"
;
import_array();
python::register_exception_translator<ValueErrorException>(&translate_value_error);
std::string docString;
docString = "Do triangle smoothing on a bounds matrix\n\n\
\n\
ARGUMENTS:\n\n\
- mat: a square Numeric array of doubles containing the bounds matrix, this matrix\n\
*is* modified by the smoothing\n\
\n\
RETURNS:\n\n\
a boolean indicating whether or not the smoothing worked.\n\
\n";
python::def("DoTriangleSmoothing", RDKit::doTriangleSmoothing,docString.c_str());
docString = "Embed a bounds matrix and return the coordinates\n\n\
\n\
ARGUMENTS:\n\n\
- boundsMatrix: a square Numeric array of doubles containing the bounds matrix, this matrix\n\
should already be smoothed\n\
- maxIters: (optional) the maximum number of random distance matrices to try\n\
- randomizeOnFailure: (optional) toggles using random coords if a matrix fails to embed\n\
- numZeroFail: (optional) sets the number of zero eigenvalues to be considered a failure\n\
- weights: (optional) a sequence of 3 sequences (i,j,weight) indicating elements of \n\
the bounds matrix whose weights should be adjusted\n\
- randomSeed: (optional) sets the random number seed used for embedding\n\
\n\
RETURNS:\n\n\
a Numeric array of doubles with the coordinates\n\
\n";
python::def("EmbedBoundsMatrix", RDKit::embedBoundsMatrix,
(python::arg("boundsMatrix"),python::arg("maxIters")=10,
python::arg("randomizeOnFailure")=false,
python::arg("numZeroFail")=2,
python::arg("weights")=python::list(),
python::arg("randomSeed")=-1),
docString.c_str());
}
// $Id$
//
// Copyright (C) 2004-2006 Rational Discovery LLC
//
// @@ All Rights Reserved @@
//
#include <boost/python.hpp>
#define PY_ARRAY_UNIQUE_SYMBOL DistGeom_array_API
#include "Numeric/arrayobject.h"
#include <RDBoost/Wrap.h>
#include <Geometry/point.h>
#include <Numerics/Matrix.h>
#include <Numerics/SymmMatrix.h>
#include <DistGeom/BoundsMatrix.h>
#include <DistGeom/TriangleSmooth.h>
#include <DistGeom/DistGeomUtils.h>
#include <DistGeom/EmbedObject.h>
#include <ForceField/ForceField.h>
#include <vector>
#include <map>
namespace python = boost::python;
namespace RDKit {
bool doTriangleSmoothing(python::object boundsMatArg){
PyObject *boundsMatObj = boundsMatArg.ptr();
if(!PyArray_Check(boundsMatObj))
throw_value_error("Argument isn't an array");
PyArrayObject *boundsMat=reinterpret_cast<PyArrayObject *>(boundsMatObj);
// get the dimensions of the array
int nrows = boundsMat->dimensions[0];
int ncols = boundsMat->dimensions[1];
if(nrows!=ncols)
throw_value_error("The array has to be square");
if(nrows<=0)
throw_value_error("The array has to have a nonzero size");
if (boundsMat->descr->type_num != PyArray_DOUBLE)
throw_value_error("Only double arrays are currently supported");
unsigned int dSize = nrows*nrows;
double *cData = new double[dSize];
double *inData = reinterpret_cast<double *>(boundsMat->data);
memcpy(static_cast<void *>(cData),
static_cast<const void *>(inData),
dSize*sizeof(double));
DistGeom::BoundsMatrix::DATA_SPTR sdata(cData);
DistGeom::BoundsMatrix bm(nrows,sdata);
bool res=DistGeom::triangleSmoothBounds(&bm);
memcpy(static_cast<void *>(inData),
static_cast<const void *>(cData),
dSize*sizeof(double));
return (PyObject *) res;
}
PyObject *embedBoundsMatrix(python::object boundsMatArg,int maxIters=10,
bool randomizeOnFailure=false,int numZeroFail=2,
python::list weights=python::list(),
int randomSeed=-1){
PyObject *boundsMatObj = boundsMatArg.ptr();
if(!PyArray_Check(boundsMatObj))
throw_value_error("Argument isn't an array");
PyArrayObject *boundsMat=reinterpret_cast<PyArrayObject *>(boundsMatObj);
// get the dimensions of the array
unsigned int nrows = boundsMat->dimensions[0];
unsigned int ncols = boundsMat->dimensions[1];
if(nrows!=ncols)
throw_value_error("The array has to be square");
if(nrows<=0)
throw_value_error("The array has to have a nonzero size");
if (boundsMat->descr->type_num != PyArray_DOUBLE)
throw_value_error("Only double arrays are currently supported");
unsigned int dSize = nrows*nrows;
double *cData = new double[dSize];
double *inData = reinterpret_cast<double *>(boundsMat->data);
memcpy(static_cast<void *>(cData),
static_cast<const void *>(inData),
dSize*sizeof(double));
DistGeom::BoundsMatrix::DATA_SPTR sdata(cData);
DistGeom::BoundsMatrix bm(nrows,sdata);
RDGeom::Point3D *positions=new RDGeom::Point3D[nrows];
std::vector<RDGeom::Point *> posPtrs;
for (unsigned int i = 0; i < nrows; i++) {
posPtrs.push_back(&positions[i]);
}
RDNumeric::DoubleSymmMatrix distMat(nrows, 0.0);
// ---- ---- ---- ---- ---- ---- ---- ---- ----
// start the embedding:
bool gotCoords=false;
for(int iter=0;iter<maxIters && !gotCoords;iter++){
// pick a random distance matrix
DistGeom::pickRandomDistMat(bm,distMat,randomSeed);
// update the seed:
if(randomSeed>=0) randomSeed+=iter*999;
//std::cerr << "Random matrix:" << std::endl;
//std::cerr << distMat;
// and embed it:
gotCoords=DistGeom::computeInitialCoords(distMat,posPtrs,randomizeOnFailure,
numZeroFail);
}
if(gotCoords){
std::map<std::pair<int,int>,double> weightMap;
unsigned int nElems=PySequence_Size(weights.ptr());
for(unsigned int entryIdx=0;entryIdx<nElems;entryIdx++){
PyObject *entry=PySequence_GetItem(weights.ptr(),entryIdx);
if(!PySequence_Check(entry) || PySequence_Size(entry)!=3){
throw_value_error("weights argument must be a sequence of 3-sequences");
}
int idx1=PyInt_AsLong(PySequence_GetItem(entry,0));
int idx2=PyInt_AsLong(PySequence_GetItem(entry,1));
double w=PyFloat_AsDouble(PySequence_GetItem(entry,2));
weightMap[std::make_pair<int,int>(idx1,idx2)]=w;
}
DistGeom::VECT_CHIRALSET csets;
ForceFields::ForceField *field = DistGeom::constructForceField(bm,posPtrs,csets,0.0, 0.0,
&weightMap);
CHECK_INVARIANT(field,"could not build dgeom force field");
field->initialize();
if(field->calcEnergy()>1e-5){
int needMore=1;
while(needMore){
needMore=field->minimize();
}
}
delete field;
} else {
throw_value_error("could not embed matrix");
}
// ---- ---- ---- ---- ---- ---- ---- ---- ----
// construct the results matrix:
int dims[2];
dims[0] = nrows;
dims[1] = 3;
PyArrayObject *res = (PyArrayObject *)PyArray_FromDims(2,dims,PyArray_DOUBLE);
double *resData=reinterpret_cast<double *>(res->data);
for(unsigned int i=0;i<nrows;i++){
unsigned int iTab=i*3;
for (unsigned int j = 0; j < 3; ++j) {
resData[iTab + j]=positions[i][j]; //.x;
//resData[iTab+1]=positions[i].y;
//resData[iTab+2]=positions[i].z;
}
}
delete [] positions;
return (PyObject *) res;
}
}
BOOST_PYTHON_MODULE(DistGeom) {
python::scope().attr("__doc__") =
"Module containing functions for basic distance geometry operations"
;
import_array();
python::register_exception_translator<ValueErrorException>(&translate_value_error);
std::string docString;
docString = "Do triangle smoothing on a bounds matrix\n\n\
\n\
ARGUMENTS:\n\n\
- mat: a square Numeric array of doubles containing the bounds matrix, this matrix\n\
*is* modified by the smoothing\n\
\n\
RETURNS:\n\n\
a boolean indicating whether or not the smoothing worked.\n\
\n";
python::def("DoTriangleSmoothing", RDKit::doTriangleSmoothing,docString.c_str());
docString = "Embed a bounds matrix and return the coordinates\n\n\
\n\
ARGUMENTS:\n\n\
- boundsMatrix: a square Numeric array of doubles containing the bounds matrix, this matrix\n\
should already be smoothed\n\
- maxIters: (optional) the maximum number of random distance matrices to try\n\
- randomizeOnFailure: (optional) toggles using random coords if a matrix fails to embed\n\
- numZeroFail: (optional) sets the number of zero eigenvalues to be considered a failure\n\
- weights: (optional) a sequence of 3 sequences (i,j,weight) indicating elements of \n\
the bounds matrix whose weights should be adjusted\n\
- randomSeed: (optional) sets the random number seed used for embedding\n\
\n\
RETURNS:\n\n\
a Numeric array of doubles with the coordinates\n\
\n";
python::def("EmbedBoundsMatrix", RDKit::embedBoundsMatrix,
(python::arg("boundsMatrix"),python::arg("maxIters")=10,
python::arg("randomizeOnFailure")=false,
python::arg("numZeroFail")=2,
python::arg("weights")=python::list(),
python::arg("randomSeed")=-1),
docString.c_str());
}

28
Code/DistGeom/Wrap/test_list.py
View File

@@ -1,14 +1,14 @@
tests=[
("python","rough_test.py",{})
]
longTests=[]
if __name__=='__main__':
import sys
import TestRunner
failed,tests = TestRunner.RunScript('test_list.py',0,1)
sys.exit(len(failed))
tests=[
("python","rough_test.py",{})
]
longTests=[]
if __name__=='__main__':
import sys
import TestRunner
failed,tests = TestRunner.RunScript('test_list.py',0,1)
sys.exit(len(failed))

216
Code/DistGeom/testDistGeom.cpp
View File

@@ -1,108 +1,108 @@
// $Id$
//
// Copyright (C) 2004-2006 Rational Discovery LLC
//
// @@ All Rights Reserved @@
//
#include "BoundsMatrix.h"
#include "TriangleSmooth.h"
#include <iostream>
#include <boost/smart_ptr.hpp>
#include <math.h>
#include <Numerics/SymmMatrix.h>
#include "DistGeomUtils.h"
#include <RDGeneral/utils.h>
using namespace DistGeom;
using namespace RDNumeric;
void test1() {
// test triangle smoothing
unsigned int npt = 5;
double x= sqrt(3.0);
BoundsMatrix * mmat = new BoundsMatrix(npt);
mmat->setUpperBound(0, 1, 1.0); mmat->setLowerBound(0, 1, 1.0);
mmat->setUpperBound(0, 2, x); mmat->setLowerBound(0, 2, x);
mmat->setUpperBound(0, 3, 10.0); mmat->setLowerBound(0, 3, 0.0);
mmat->setUpperBound(0, 4, 10.0); mmat->setLowerBound(0, 4, 0.0);
mmat->setUpperBound(1, 2, 1.0); mmat->setLowerBound(1, 2, 1.0);
mmat->setUpperBound(1, 3, x); mmat->setLowerBound(1, 3, x);
mmat->setUpperBound(1, 4, 10.0); mmat->setLowerBound(1, 4, 0.0);
mmat->setUpperBound(2, 3, 1.0); mmat->setLowerBound(2, 3, 1.0);
mmat->setUpperBound(2, 4, x); mmat->setLowerBound(2, 4, x);
mmat->setUpperBound(3, 4, 1.0); mmat->setLowerBound(3, 4, 1.0);
BoundsMatPtr mptr(mmat);
triangleSmoothBounds(mptr);
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(0, 1), 1.0, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(0, 1), 1.0, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(0, 2), 1.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(0, 2), 1.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(0, 3), 2.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(0, 3), 0.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(0, 4), 3.464, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(0, 4), 0.0, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(1, 2), 1.0, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(1, 2), 1.0, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(1, 3), 1.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(1, 3), 1.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(1, 4), 2.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(1, 4), 0.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(2, 3), 1.0, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(2, 3), 1.0, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(2, 4), 1.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(2, 4), 1.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(3, 4), 1.0, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(3, 4), 1.0, 0.001), "");
DoubleSymmMatrix dmat(npt, 0.0);
pickRandomDistMat(*mmat, dmat, 100);
double sumElem = 0.0;
for (unsigned int i = 0; i < dmat.getDataSize(); i++) {
sumElem += dmat.getData()[i];
}
CHECK_INVARIANT(RDKit::feq(sumElem, 14.3079, 0.001), "");
}
void testIssue216() {
RDNumeric::DoubleSymmMatrix dmat(4);
dmat.setVal(0,0, 0.0); dmat.setVal(0,1, 1.0); dmat.setVal(0,2, 1.0); dmat.setVal(0,3, 1.0);
dmat.setVal(1,1, 0.0); dmat.setVal(1,2, 1.0); dmat.setVal(1,3, 1.0);
dmat.setVal(2,2, 0.0); dmat.setVal(2,3, 1.0);
dmat.setVal(3,3, 0.0);
std::cout << dmat;
RDGeom::PointPtrVect pos;
for (int i = 0; i < 4; i++) {
RDGeom::Point3D *pt = new RDGeom::Point3D();
pos.push_back(pt);
}
bool gotCoords = DistGeom::computeInitialCoords(dmat, pos);
CHECK_INVARIANT(gotCoords, "");
for (int i = 1; i < 4; i++) {
RDGeom::Point3D pti = *(RDGeom::Point3D *)pos[i];
for (int j = 0; j < i; j++) {
RDGeom::Point3D ptj = *(RDGeom::Point3D *)pos[j];
ptj -= pti;
CHECK_INVARIANT(RDKit::feq(ptj.length(), 1.0, 0.02), "");
}
}
}
int main() {
std::cout << "***********************************************************\n";
std::cout << " test1 \n";
test1();
std::cout << "***********************************************************\n";
std::cout << " testIssue216 \n";
testIssue216();
std::cout << "***********************************************************\n\n";
return 0;
}
// $Id$
//
// Copyright (C) 2004-2006 Rational Discovery LLC
//
// @@ All Rights Reserved @@
//
#include "BoundsMatrix.h"
#include "TriangleSmooth.h"
#include <iostream>
#include <boost/smart_ptr.hpp>
#include <math.h>
#include <Numerics/SymmMatrix.h>
#include "DistGeomUtils.h"
#include <RDGeneral/utils.h>
using namespace DistGeom;
using namespace RDNumeric;
void test1() {
// test triangle smoothing
unsigned int npt = 5;
double x= sqrt(3.0);
BoundsMatrix * mmat = new BoundsMatrix(npt);
mmat->setUpperBound(0, 1, 1.0); mmat->setLowerBound(0, 1, 1.0);
mmat->setUpperBound(0, 2, x); mmat->setLowerBound(0, 2, x);
mmat->setUpperBound(0, 3, 10.0); mmat->setLowerBound(0, 3, 0.0);
mmat->setUpperBound(0, 4, 10.0); mmat->setLowerBound(0, 4, 0.0);
mmat->setUpperBound(1, 2, 1.0); mmat->setLowerBound(1, 2, 1.0);
mmat->setUpperBound(1, 3, x); mmat->setLowerBound(1, 3, x);
mmat->setUpperBound(1, 4, 10.0); mmat->setLowerBound(1, 4, 0.0);
mmat->setUpperBound(2, 3, 1.0); mmat->setLowerBound(2, 3, 1.0);
mmat->setUpperBound(2, 4, x); mmat->setLowerBound(2, 4, x);
mmat->setUpperBound(3, 4, 1.0); mmat->setLowerBound(3, 4, 1.0);
BoundsMatPtr mptr(mmat);
triangleSmoothBounds(mptr);
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(0, 1), 1.0, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(0, 1), 1.0, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(0, 2), 1.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(0, 2), 1.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(0, 3), 2.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(0, 3), 0.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(0, 4), 3.464, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(0, 4), 0.0, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(1, 2), 1.0, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(1, 2), 1.0, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(1, 3), 1.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(1, 3), 1.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(1, 4), 2.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(1, 4), 0.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(2, 3), 1.0, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(2, 3), 1.0, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(2, 4), 1.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(2, 4), 1.732, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getUpperBound(3, 4), 1.0, 0.001), "");
CHECK_INVARIANT(RDKit::feq(mmat->getLowerBound(3, 4), 1.0, 0.001), "");
DoubleSymmMatrix dmat(npt, 0.0);
pickRandomDistMat(*mmat, dmat, 100);
double sumElem = 0.0;
for (unsigned int i = 0; i < dmat.getDataSize(); i++) {
sumElem += dmat.getData()[i];
}
CHECK_INVARIANT(RDKit::feq(sumElem, 14.3079, 0.001), "");
}
void testIssue216() {
RDNumeric::DoubleSymmMatrix dmat(4);
dmat.setVal(0,0, 0.0); dmat.setVal(0,1, 1.0); dmat.setVal(0,2, 1.0); dmat.setVal(0,3, 1.0);
dmat.setVal(1,1, 0.0); dmat.setVal(1,2, 1.0); dmat.setVal(1,3, 1.0);
dmat.setVal(2,2, 0.0); dmat.setVal(2,3, 1.0);
dmat.setVal(3,3, 0.0);
std::cout << dmat;
RDGeom::PointPtrVect pos;
for (int i = 0; i < 4; i++) {
RDGeom::Point3D *pt = new RDGeom::Point3D();
pos.push_back(pt);
}
bool gotCoords = DistGeom::computeInitialCoords(dmat, pos);
CHECK_INVARIANT(gotCoords, "");
for (int i = 1; i < 4; i++) {
RDGeom::Point3D pti = *(RDGeom::Point3D *)pos[i];
for (int j = 0; j < i; j++) {
RDGeom::Point3D ptj = *(RDGeom::Point3D *)pos[j];
ptj -= pti;
CHECK_INVARIANT(RDKit::feq(ptj.length(), 1.0, 0.02), "");
}
}
}
int main() {
std::cout << "***********************************************************\n";
std::cout << " test1 \n";
test1();
std::cout << "***********************************************************\n";
std::cout << " testIssue216 \n";
testIssue216();
std::cout << "***********************************************************\n\n";
return 0;
}