From a643a9913aab14f6a70da5addf1c8db626175f0f Mon Sep 17 00:00:00 2001 From: Greg Landrum Date: Sun, 18 May 2008 05:06:57 +0000 Subject: [PATCH] dos2unix and untabify... sigh --- Code/DistGeom/BoundsMatrix.h | 228 ++++++------ Code/DistGeom/ChiralViolationContrib.cpp | 38 +- Code/DistGeom/DistGeomUtils.h | 180 +++++----- Code/DistGeom/DistViolationContrib.cpp | 164 ++++----- Code/DistGeom/DistViolationContrib.h | 84 ++--- Code/DistGeom/EmbedObject.h | 202 +++++------ Code/DistGeom/FourthDimContrib.h | 2 +- Code/DistGeom/TriangleSmooth.cpp | 108 +++--- Code/DistGeom/TriangleSmooth.h | 60 ++-- Code/DistGeom/Wrap/DistGeom.cpp | 432 +++++++++++------------ Code/DistGeom/Wrap/test_list.py | 28 +- Code/DistGeom/testDistGeom.cpp | 216 ++++++------ 12 files changed, 871 insertions(+), 871 deletions(-) diff --git a/Code/DistGeom/BoundsMatrix.h b/Code/DistGeom/BoundsMatrix.h index e8310f57d..77c216c43 100644 --- a/Code/DistGeom/BoundsMatrix.h +++ b/Code/DistGeom/BoundsMatrix.h @@ -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 -#include -#include -#include -#include - -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 { - public: - typedef boost::shared_array DATA_SPTR; - - explicit BoundsMatrix(unsigned int N) : RDNumeric::SquareMatrix(N,0.0) {}; - BoundsMatrix(unsigned int N, DATA_SPTR data) : - RDNumeric::SquareMatrix(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 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 +#include +#include +#include +#include + +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 { + public: + typedef boost::shared_array DATA_SPTR; + + explicit BoundsMatrix(unsigned int N) : RDNumeric::SquareMatrix(N,0.0) {}; + BoundsMatrix(unsigned int N, DATA_SPTR data) : + RDNumeric::SquareMatrix(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 BoundsMatPtr; +} + +#endif + diff --git a/Code/DistGeom/ChiralViolationContrib.cpp b/Code/DistGeom/ChiralViolationContrib.cpp index e7a87b263..3bcb3d291 100644 --- a/Code/DistGeom/ChiralViolationContrib.cpp +++ b/Code/DistGeom/ChiralViolationContrib.cpp @@ -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])); } } diff --git a/Code/DistGeom/DistGeomUtils.h b/Code/DistGeom/DistGeomUtils.h index 65663f698..edf74665b 100644 --- a/Code/DistGeom/DistGeomUtils.h +++ b/Code/DistGeom/DistGeomUtils.h @@ -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 -#include -#include -#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 &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 &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. - NOTE: 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,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 +#include +#include +#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 &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 &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. + NOTE: 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,double> *extraWeights=0, + double basinSizeTol=5.0); + +} + +#endif diff --git a/Code/DistGeom/DistViolationContrib.cpp b/Code/DistGeom/DistViolationContrib.cpp index 6e7b7beac..7ce1e90ab 100644 --- a/Code/DistGeom/DistViolationContrib.cpp +++ b/Code/DistGeom/DistViolationContrib.cpp @@ -1,82 +1,82 @@ -// $Id$ -// -// Copyright (C) 2004-2006 Rational Discovery LLC -// -// @@ All Rights Reserved @@ -// - -#include "DistViolationContrib.h" -#include -#include - -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;i0.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 +#include + +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;i0.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; + } + } +} diff --git a/Code/DistGeom/DistViolationContrib.h b/Code/DistGeom/DistViolationContrib.h index 899b8e199..00f91aace 100644 --- a/Code/DistGeom/DistViolationContrib.h +++ b/Code/DistGeom/DistViolationContrib.h @@ -1,42 +1,42 @@ -// -// Copyright (C) 2004-2006 Rational Discovery LLC -// -// @@ All Rights Reserved @@ -// -#ifndef __RD_DISTVIOLATIONCONTRIB_H__ -#define __RD_DISTVIOLATIONCONTRIB_H__ - -#include - -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 + +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 diff --git a/Code/DistGeom/EmbedObject.h b/Code/DistGeom/EmbedObject.h index a746f61fc..0bfe36959 100644 --- a/Code/DistGeom/EmbedObject.h +++ b/Code/DistGeom/EmbedObject.h @@ -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 -#include -#include - -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 ChiralSetPtr; - typedef std::vector 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 +#include +#include + +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 ChiralSetPtr; + typedef std::vector 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 diff --git a/Code/DistGeom/FourthDimContrib.h b/Code/DistGeom/FourthDimContrib.h index 373e07fc8..264dfa02d 100644 --- a/Code/DistGeom/FourthDimContrib.h +++ b/Code/DistGeom/FourthDimContrib.h @@ -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; diff --git a/Code/DistGeom/TriangleSmooth.cpp b/Code/DistGeom/TriangleSmooth.cpp index a14fcaddf..c2d3d2842 100644 --- a/Code/DistGeom/TriangleSmooth.cpp +++ b/Code/DistGeom/TriangleSmooth.cpp @@ -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; + } +} + diff --git a/Code/DistGeom/TriangleSmooth.h b/Code/DistGeom/TriangleSmooth.h index 75979230b..55548ece9 100644 --- a/Code/DistGeom/TriangleSmooth.h +++ b/Code/DistGeom/TriangleSmooth.h @@ -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 + diff --git a/Code/DistGeom/Wrap/DistGeom.cpp b/Code/DistGeom/Wrap/DistGeom.cpp index 72e8060f7..0c1564c1e 100644 --- a/Code/DistGeom/Wrap/DistGeom.cpp +++ b/Code/DistGeom/Wrap/DistGeom.cpp @@ -1,216 +1,216 @@ -// $Id$ -// -// Copyright (C) 2004-2006 Rational Discovery LLC -// -// @@ All Rights Reserved @@ -// - -#include - -#define PY_ARRAY_UNIQUE_SYMBOL DistGeom_array_API -#include "Numeric/arrayobject.h" -#include - -#include -#include -#include - -#include -#include -#include -#include - -#include - -#include -#include - -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(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(boundsMat->data); - memcpy(static_cast(cData), - static_cast(inData), - dSize*sizeof(double)); - DistGeom::BoundsMatrix::DATA_SPTR sdata(cData); - DistGeom::BoundsMatrix bm(nrows,sdata); - - bool res=DistGeom::triangleSmoothBounds(&bm); - memcpy(static_cast(inData), - static_cast(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(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(boundsMat->data); - memcpy(static_cast(cData), - static_cast(inData), - dSize*sizeof(double)); - - DistGeom::BoundsMatrix::DATA_SPTR sdata(cData); - DistGeom::BoundsMatrix bm(nrows,sdata); - - RDGeom::Point3D *positions=new RDGeom::Point3D[nrows]; - std::vector 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=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,double> weightMap; - unsigned int nElems=PySequence_Size(weights.ptr()); - for(unsigned int entryIdx=0;entryIdx(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(res->data); - for(unsigned int i=0;i(&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 + +#define PY_ARRAY_UNIQUE_SYMBOL DistGeom_array_API +#include "Numeric/arrayobject.h" +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include + +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(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(boundsMat->data); + memcpy(static_cast(cData), + static_cast(inData), + dSize*sizeof(double)); + DistGeom::BoundsMatrix::DATA_SPTR sdata(cData); + DistGeom::BoundsMatrix bm(nrows,sdata); + + bool res=DistGeom::triangleSmoothBounds(&bm); + memcpy(static_cast(inData), + static_cast(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(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(boundsMat->data); + memcpy(static_cast(cData), + static_cast(inData), + dSize*sizeof(double)); + + DistGeom::BoundsMatrix::DATA_SPTR sdata(cData); + DistGeom::BoundsMatrix bm(nrows,sdata); + + RDGeom::Point3D *positions=new RDGeom::Point3D[nrows]; + std::vector 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=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,double> weightMap; + unsigned int nElems=PySequence_Size(weights.ptr()); + for(unsigned int entryIdx=0;entryIdx(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(res->data); + for(unsigned int i=0;i(&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()); + + +} diff --git a/Code/DistGeom/Wrap/test_list.py b/Code/DistGeom/Wrap/test_list.py index 319e754e5..77c6dab1e 100644 --- a/Code/DistGeom/Wrap/test_list.py +++ b/Code/DistGeom/Wrap/test_list.py @@ -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)) diff --git a/Code/DistGeom/testDistGeom.cpp b/Code/DistGeom/testDistGeom.cpp index 28f840d42..95ec4dfec 100644 --- a/Code/DistGeom/testDistGeom.cpp +++ b/Code/DistGeom/testDistGeom.cpp @@ -1,108 +1,108 @@ -// $Id$ -// -// Copyright (C) 2004-2006 Rational Discovery LLC -// -// @@ All Rights Reserved @@ -// - - -#include "BoundsMatrix.h" -#include "TriangleSmooth.h" -#include -#include -#include -#include -#include "DistGeomUtils.h" -#include - -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 +#include +#include +#include +#include "DistGeomUtils.h" +#include + +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; +}