mirror of
https://github.com/rdkit/rdkit.git
synced 2026-06-03 21:44:30 +08:00
91008ff11d
* Address compile warnings & trivial improvements * revert unwanted initializers; use RDUNUSED_PARAM for unused params * revert fix in testRDFcustom; marked with 'TO DO' comment
236 lines
8.5 KiB
C++
236 lines
8.5 KiB
C++
//
|
|
// Copyright (c) 2016, Guillaume GODIN.
|
|
// All rights reserved.
|
|
//
|
|
// Redistribution and use in source and binary forms, with or without
|
|
// modification, are permitted provided that the following conditions are
|
|
// met:
|
|
//
|
|
// * Redistributions of source code must retain the above copyright
|
|
// notice, this list of conditions and the following disclaimer.
|
|
// * Redistributions in binary form must reproduce the above
|
|
// copyright notice, this list of conditions and the following
|
|
// disclaimer in the documentation and/or other materials provided
|
|
// with the distribution.
|
|
// * Neither the name of Institue of Cancer Research.
|
|
// nor the names of its contributors may be used to endorse or promote
|
|
// products derived from this software without specific prior written
|
|
// permission.
|
|
//
|
|
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
|
|
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
|
|
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
|
|
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
|
|
// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
|
|
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
|
|
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
|
|
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
|
|
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
|
|
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
|
|
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
|
//
|
|
// Guillaume GODIN access the AutoCorrelation 3D descriptors names in Dragon TDB
|
|
|
|
#include <GraphMol/RDKitBase.h>
|
|
|
|
#include "AUTOCORR3D.h"
|
|
#include "MolData3Ddescriptors.h"
|
|
|
|
#include <math.h>
|
|
#include <Eigen/Dense>
|
|
#include <Eigen/SVD>
|
|
#include <iostream>
|
|
#include <Eigen/Core>
|
|
#include <Eigen/QR>
|
|
#include <boost/math/special_functions/fpclassify.hpp>
|
|
|
|
using namespace Eigen;
|
|
namespace RDKit {
|
|
namespace Descriptors {
|
|
|
|
namespace {
|
|
|
|
MolData3Ddescriptors moldata3D;
|
|
|
|
VectorXd getEigenVect(std::vector<double> v) {
|
|
double* varray_ptr = &v[0];
|
|
Map<VectorXd> V(varray_ptr, v.size());
|
|
return V;
|
|
}
|
|
|
|
double* GetGeodesicMatrix(double* dist, int lag, int numAtoms) {
|
|
int sizeArray = numAtoms * numAtoms;
|
|
double* Geodesic = new double[sizeArray];
|
|
for (int i = 0; i < sizeArray; i++) {
|
|
if (dist[i] == lag)
|
|
Geodesic[i] = 1.0;
|
|
else
|
|
Geodesic[i] = 0.0;
|
|
}
|
|
return Geodesic;
|
|
}
|
|
|
|
// matrix prod that mimic the v2 version
|
|
// the code is in respect to Dragon 6 descriptors.
|
|
// replace the number of "Bicount" vertex per lag by a (numAtoms * (numAtoms -
|
|
// 1))!
|
|
// provided by Kobe team!
|
|
void get3DautocorrelationDesc(double* dist3D, double* topologicaldistance,
|
|
int numAtoms, const ROMol& mol,
|
|
std::vector<double>& res) {
|
|
Map<MatrixXd> dm(dist3D, numAtoms, numAtoms);
|
|
Map<MatrixXd> di(topologicaldistance, numAtoms, numAtoms);
|
|
|
|
std::vector<double> wp = moldata3D.GetRelativePol(mol);
|
|
VectorXd Wp = getEigenVect(wp);
|
|
|
|
std::vector<double> wm = moldata3D.GetRelativeMW(mol);
|
|
VectorXd Wm = getEigenVect(wm);
|
|
|
|
std::vector<double> wi = moldata3D.GetRelativeIonPol(mol);
|
|
VectorXd Wi = getEigenVect(wi);
|
|
|
|
std::vector<double> wv = moldata3D.GetRelativeVdW(mol);
|
|
VectorXd Wv = getEigenVect(wv);
|
|
|
|
std::vector<double> we = moldata3D.GetRelativeENeg(mol);
|
|
VectorXd We = getEigenVect(we);
|
|
|
|
std::vector<double> wu = moldata3D.GetUn(numAtoms);
|
|
VectorXd Wu = getEigenVect(wu);
|
|
|
|
std::vector<double> ws = moldata3D.GetIState(mol);
|
|
VectorXd Ws = getEigenVect(ws);
|
|
|
|
std::vector<double> wr = moldata3D.GetRelativeRcov(mol);
|
|
VectorXd Wr = getEigenVect(wr);
|
|
|
|
MatrixXd Bi;
|
|
MatrixXd tmp;
|
|
double TDBmat[8][10];
|
|
double dtmp;
|
|
|
|
for (int i = 0; i < 10; i++) {
|
|
double* Bimat = GetGeodesicMatrix(topologicaldistance, i + 1, numAtoms);
|
|
Map<MatrixXd> Bi(Bimat, numAtoms, numAtoms);
|
|
MatrixXd RBi = Bi.cwiseProduct(dm);
|
|
// double Bicount = (double)Bi.sum();
|
|
|
|
tmp = Wu.transpose() * RBi * Wu;
|
|
dtmp = (double)tmp(0);
|
|
if (boost::math::isnan(dtmp)) dtmp = 0.0;
|
|
TDBmat[0][i] = dtmp;
|
|
|
|
tmp = Wm.transpose() * RBi * Wm;
|
|
dtmp = (double)tmp(0);
|
|
if (boost::math::isnan(dtmp)) dtmp = 0.0;
|
|
TDBmat[1][i] = dtmp;
|
|
|
|
tmp = Wv.transpose() * RBi * Wv;
|
|
dtmp = (double)tmp(0);
|
|
if (boost::math::isnan(dtmp)) dtmp = 0.0;
|
|
TDBmat[2][i] = dtmp;
|
|
|
|
tmp = We.transpose() * RBi * We;
|
|
dtmp = (double)tmp(0);
|
|
if (boost::math::isnan(dtmp)) dtmp = 0.0;
|
|
TDBmat[3][i] = dtmp;
|
|
|
|
tmp = Wp.transpose() * RBi * Wp;
|
|
dtmp = (double)tmp(0);
|
|
if (boost::math::isnan(dtmp)) dtmp = 0.0;
|
|
TDBmat[4][i] = dtmp;
|
|
|
|
tmp = Wi.transpose() * RBi * Wi;
|
|
dtmp = (double)tmp(0);
|
|
if (boost::math::isnan(dtmp)) dtmp = 0.0;
|
|
TDBmat[5][i] = dtmp;
|
|
|
|
tmp = Ws.transpose() * RBi * Ws;
|
|
dtmp = (double)tmp(0);
|
|
if (boost::math::isnan(dtmp)) dtmp = 0.0;
|
|
TDBmat[6][i] = dtmp;
|
|
|
|
tmp = Wr.transpose() * RBi * Wr;
|
|
dtmp = (double)tmp(0);
|
|
if (boost::math::isnan(dtmp)) dtmp = 0.0;
|
|
TDBmat[7][i] = dtmp;
|
|
delete[] Bimat;
|
|
}
|
|
|
|
// update the Output vector!
|
|
for (unsigned int j = 0; j < 8; ++j) {
|
|
for (unsigned int i = 0; i < 10; ++i) {
|
|
res[j * 10 + i] =
|
|
round(1000 * TDBmat[j][i] / (numAtoms * (numAtoms - 1))) / 1000;
|
|
}
|
|
}
|
|
}
|
|
|
|
void get3DautocorrelationDescCustom(double* dist3D, double* topologicaldistance,
|
|
int numAtoms, const ROMol& mol,
|
|
std::vector<double>& res,
|
|
const std::string &customAtomPropName) {
|
|
Map<MatrixXd> dm(dist3D, numAtoms, numAtoms);
|
|
Map<MatrixXd> di(topologicaldistance, numAtoms, numAtoms);
|
|
|
|
std::vector<double> customAtomArray = moldata3D.GetCustomAtomProp(mol,customAtomPropName);
|
|
VectorXd Wc = getEigenVect(customAtomArray);
|
|
|
|
MatrixXd Bi;
|
|
MatrixXd tmp;
|
|
double TDBmat[10];
|
|
double dtmp;
|
|
|
|
for (int i = 0; i < 10; i++) {
|
|
double* Bimat = GetGeodesicMatrix(topologicaldistance, i + 1, numAtoms);
|
|
Map<MatrixXd> Bi(Bimat, numAtoms, numAtoms);
|
|
MatrixXd RBi = Bi.cwiseProduct(dm);
|
|
|
|
tmp = Wc.transpose() * RBi * Wc;
|
|
dtmp = (double)tmp(0);
|
|
if (std::isnan(dtmp)) dtmp = 0.0;
|
|
TDBmat[i] = dtmp;
|
|
delete[] Bimat;
|
|
}
|
|
// update the Output vector!
|
|
for (unsigned int i = 0; i < 10; ++i) {
|
|
res[i] = round(1000 * TDBmat[i] / (numAtoms * (numAtoms - 1))) / 1000;
|
|
}
|
|
}
|
|
|
|
void Get3Dauto(double* dist3D, double* topologicaldistance, int numAtoms,
|
|
const ROMol& mol, std::vector<double>& res) {
|
|
// AUTOCORRNAMES={"TDB01u","TDB02u","TDB03u","TDB04u","TDB05u","TDB06u","TDB07u","TDB08u","TDB09u","TDB10u","TDB01m","TDB02m","TDB03m","TDB04m","TDB05m","TDB06m","TDB07m","TDB08m","TDB09m","TDB10m","TDB01v","TDB02v","TDB03v","TDB04v","TDB05v","TDB06v","TDB07v","TDB08v","TDB09v","TDB10v","TDB01e","TDB02e","TDB03e","TDB04e","TDB05e","TDB06e","TDB07e","TDB08e","TDB09e","TDB10e","TDB01p","TDB02p","TDB03p","TDB04p","TDB05p","TDB06p","TDB07p","TDB08p","TDB09p","TDB10p","TDB01i","TDB02i","TDB03i","TDB04i","TDB05i","TDB06i","TDB07i","TDB08i","TDB09i","TDB10i","TDB01s","TDB02s","TDB03s","TDB04s","TDB05s","TDB06s","TDB07s","TDB08s","TDB09s","TDB10s","TDB01r","TDB02r","TDB03r","TDB04r","TDB05r","TDB06r","TDB07r","TDB08r","TDB09r","TDB10r"};
|
|
get3DautocorrelationDesc(dist3D, topologicaldistance, numAtoms, mol, res);
|
|
}
|
|
|
|
void Get3Dautoone(double* dist3D, double* topologicaldistance, int numAtoms,
|
|
const ROMol& mol, std::vector<double>& res, const std::string &customAtomPropName) {
|
|
get3DautocorrelationDescCustom(dist3D, topologicaldistance, numAtoms, mol, res, customAtomPropName);
|
|
}
|
|
|
|
} // end of anonymous namespace
|
|
|
|
void AUTOCORR3D(const ROMol& mol, std::vector<double>& res, int confId,
|
|
const std::string &customAtomPropName) {
|
|
PRECONDITION(mol.getNumConformers() >= 1, "molecule has no conformers")
|
|
int numAtoms = mol.getNumAtoms();
|
|
|
|
double* topologicaldistance =
|
|
MolOps::getDistanceMat(mol, false); // topological matrix
|
|
double* dist3D =
|
|
MolOps::get3DDistanceMat(mol, confId, false, true); // 3D distance matrix
|
|
if (customAtomPropName != "") {
|
|
res.clear();
|
|
res.resize(10);
|
|
Get3Dautoone(dist3D, topologicaldistance, numAtoms, mol, res, customAtomPropName);
|
|
} else {
|
|
res.clear();
|
|
res.resize(80);
|
|
Get3Dauto(dist3D, topologicaldistance, numAtoms, mol, res);
|
|
}
|
|
}
|
|
} // end of Descriptors namespace
|
|
} // end of RDKit namespace
|