1 //
2 // Copyright (c) 2016, Guillaume GODIN.
3 // All rights reserved.
4 //
5 // Redistribution and use in source and binary forms, with or without
6 // modification, are permitted provided that the following conditions are
7 // met:
8 //
9 // * Redistributions of source code must retain the above copyright
10 // notice, this list of conditions and the following disclaimer.
11 // * Redistributions in binary form must reproduce the above
12 // copyright notice, this list of conditions and the following
13 // disclaimer in the documentation and/or other materials provided
14 // with the distribution.
15 // * Neither the name of Institue of Cancer Research.
16 // nor the names of its contributors may be used to endorse or promote
17 // products derived from this software without specific prior written
18 // permission.
19 //
20 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 //
32 // Guillaume GODIN access the AutoCorrelation 3D descriptors names in Dragon TDB
33
34 #include <GraphMol/RDKitBase.h>
35
36 #include "AUTOCORR3D.h"
37 #include "MolData3Ddescriptors.h"
38
39 #include <cmath>
40 #include <Eigen/Dense>
41 #include <Eigen/SVD>
42 #include <iostream>
43 #include <Eigen/Core>
44 #include <Eigen/QR>
45
46 using namespace Eigen;
47 namespace RDKit {
48 namespace Descriptors {
49
50 namespace {
51
52 MolData3Ddescriptors moldata3D;
53
getEigenVect(std::vector<double> v)54 VectorXd getEigenVect(std::vector<double> v) {
55 double* varray_ptr = &v[0];
56 Map<VectorXd> V(varray_ptr, v.size());
57 return V;
58 }
59
GetGeodesicMatrix(double * dist,int lag,int numAtoms)60 double* GetGeodesicMatrix(double* dist, int lag, int numAtoms) {
61 int sizeArray = numAtoms * numAtoms;
62 auto* Geodesic = new double[sizeArray];
63 for (int i = 0; i < sizeArray; i++) {
64 if (dist[i] == lag) {
65 Geodesic[i] = 1.0;
66 } else {
67 Geodesic[i] = 0.0;
68 }
69 }
70 return Geodesic;
71 }
72
73 // matrix prod that mimic the v2 version
74 // the code is in respect to Dragon 6 descriptors.
75 // replace the number of "Bicount" vertex per lag by a (numAtoms * (numAtoms -
76 // 1))!
77 // provided by Kobe team!
get3DautocorrelationDesc(double * dist3D,double * topologicaldistance,int numAtoms,const ROMol & mol,std::vector<double> & res)78 void get3DautocorrelationDesc(double* dist3D, double* topologicaldistance,
79 int numAtoms, const ROMol& mol,
80 std::vector<double>& res) {
81 Map<MatrixXd> dm(dist3D, numAtoms, numAtoms);
82 Map<MatrixXd> di(topologicaldistance, numAtoms, numAtoms);
83
84 std::vector<double> wp = moldata3D.GetRelativePol(mol);
85 VectorXd Wp = getEigenVect(wp);
86
87 std::vector<double> wm = moldata3D.GetRelativeMW(mol);
88 VectorXd Wm = getEigenVect(wm);
89
90 std::vector<double> wi = moldata3D.GetRelativeIonPol(mol);
91 VectorXd Wi = getEigenVect(wi);
92
93 std::vector<double> wv = moldata3D.GetRelativeVdW(mol);
94 VectorXd Wv = getEigenVect(wv);
95
96 std::vector<double> we = moldata3D.GetRelativeENeg(mol);
97 VectorXd We = getEigenVect(we);
98
99 std::vector<double> wu = moldata3D.GetUn(numAtoms);
100 VectorXd Wu = getEigenVect(wu);
101
102 std::vector<double> ws = moldata3D.GetIState(mol);
103 VectorXd Ws = getEigenVect(ws);
104
105 std::vector<double> wr = moldata3D.GetRelativeRcov(mol);
106 VectorXd Wr = getEigenVect(wr);
107
108 MatrixXd Bi;
109 MatrixXd tmp;
110 double TDBmat[8][10];
111 double dtmp;
112
113 for (int i = 0; i < 10; i++) {
114 double* Bimat = GetGeodesicMatrix(topologicaldistance, i + 1, numAtoms);
115 Map<MatrixXd> Bi(Bimat, numAtoms, numAtoms);
116 MatrixXd RBi = Bi.cwiseProduct(dm);
117 // double Bicount = (double)Bi.sum();
118
119 tmp = Wu.transpose() * RBi * Wu;
120 dtmp = (double)tmp(0);
121 if (std::isnan(dtmp)) {
122 dtmp = 0.0;
123 }
124 TDBmat[0][i] = dtmp;
125
126 tmp = Wm.transpose() * RBi * Wm;
127 dtmp = (double)tmp(0);
128 if (std::isnan(dtmp)) {
129 dtmp = 0.0;
130 }
131 TDBmat[1][i] = dtmp;
132
133 tmp = Wv.transpose() * RBi * Wv;
134 dtmp = (double)tmp(0);
135 if (std::isnan(dtmp)) {
136 dtmp = 0.0;
137 }
138 TDBmat[2][i] = dtmp;
139
140 tmp = We.transpose() * RBi * We;
141 dtmp = (double)tmp(0);
142 if (std::isnan(dtmp)) {
143 dtmp = 0.0;
144 }
145 TDBmat[3][i] = dtmp;
146
147 tmp = Wp.transpose() * RBi * Wp;
148 dtmp = (double)tmp(0);
149 if (std::isnan(dtmp)) {
150 dtmp = 0.0;
151 }
152 TDBmat[4][i] = dtmp;
153
154 tmp = Wi.transpose() * RBi * Wi;
155 dtmp = (double)tmp(0);
156 if (std::isnan(dtmp)) {
157 dtmp = 0.0;
158 }
159 TDBmat[5][i] = dtmp;
160
161 tmp = Ws.transpose() * RBi * Ws;
162 dtmp = (double)tmp(0);
163 if (std::isnan(dtmp)) {
164 dtmp = 0.0;
165 }
166 TDBmat[6][i] = dtmp;
167
168 tmp = Wr.transpose() * RBi * Wr;
169 dtmp = (double)tmp(0);
170 if (std::isnan(dtmp)) {
171 dtmp = 0.0;
172 }
173 TDBmat[7][i] = dtmp;
174 delete[] Bimat;
175 }
176
177 // update the Output vector!
178 for (unsigned int j = 0; j < 8; ++j) {
179 for (unsigned int i = 0; i < 10; ++i) {
180 res[j * 10 + i] =
181 std::round(1000 * TDBmat[j][i] / (numAtoms * (numAtoms - 1))) / 1000;
182 }
183 }
184 }
185
get3DautocorrelationDescCustom(double * dist3D,double * topologicaldistance,int numAtoms,const ROMol & mol,std::vector<double> & res,const std::string & customAtomPropName)186 void get3DautocorrelationDescCustom(double* dist3D, double* topologicaldistance,
187 int numAtoms, const ROMol& mol,
188 std::vector<double>& res,
189 const std::string& customAtomPropName) {
190 Map<MatrixXd> dm(dist3D, numAtoms, numAtoms);
191 Map<MatrixXd> di(topologicaldistance, numAtoms, numAtoms);
192
193 std::vector<double> customAtomArray =
194 moldata3D.GetCustomAtomProp(mol, customAtomPropName);
195 VectorXd Wc = getEigenVect(customAtomArray);
196
197 MatrixXd Bi;
198 MatrixXd tmp;
199 double TDBmat[10];
200 double dtmp;
201
202 for (int i = 0; i < 10; i++) {
203 double* Bimat = GetGeodesicMatrix(topologicaldistance, i + 1, numAtoms);
204 Map<MatrixXd> Bi(Bimat, numAtoms, numAtoms);
205 MatrixXd RBi = Bi.cwiseProduct(dm);
206
207 tmp = Wc.transpose() * RBi * Wc;
208 dtmp = (double)tmp(0);
209 if (std::isnan(dtmp)) {
210 dtmp = 0.0;
211 }
212 TDBmat[i] = dtmp;
213 delete[] Bimat;
214 }
215 // update the Output vector!
216 for (unsigned int i = 0; i < 10; ++i) {
217 res[i] = std::round(1000 * TDBmat[i] / (numAtoms * (numAtoms - 1))) / 1000;
218 }
219 }
220
Get3Dauto(double * dist3D,double * topologicaldistance,int numAtoms,const ROMol & mol,std::vector<double> & res)221 void Get3Dauto(double* dist3D, double* topologicaldistance, int numAtoms,
222 const ROMol& mol, std::vector<double>& res) {
223 // 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"};
224 get3DautocorrelationDesc(dist3D, topologicaldistance, numAtoms, mol, res);
225 }
226
Get3Dautoone(double * dist3D,double * topologicaldistance,int numAtoms,const ROMol & mol,std::vector<double> & res,const std::string & customAtomPropName)227 void Get3Dautoone(double* dist3D, double* topologicaldistance, int numAtoms,
228 const ROMol& mol, std::vector<double>& res,
229 const std::string& customAtomPropName) {
230 get3DautocorrelationDescCustom(dist3D, topologicaldistance, numAtoms, mol,
231 res, customAtomPropName);
232 }
233
234 } // end of anonymous namespace
235
AUTOCORR3D(const ROMol & mol,std::vector<double> & res,int confId,const std::string & customAtomPropName)236 void AUTOCORR3D(const ROMol& mol, std::vector<double>& res, int confId,
237 const std::string& customAtomPropName) {
238 PRECONDITION(mol.getNumConformers() >= 1, "molecule has no conformers")
239 int numAtoms = mol.getNumAtoms();
240
241 double* topologicaldistance =
242 MolOps::getDistanceMat(mol, false); // topological matrix
243 double* dist3D =
244 MolOps::get3DDistanceMat(mol, confId, false, true); // 3D distance matrix
245 if (customAtomPropName != "") {
246 res.clear();
247 res.resize(10);
248 Get3Dautoone(dist3D, topologicaldistance, numAtoms, mol, res,
249 customAtomPropName);
250 } else {
251 res.clear();
252 res.resize(80);
253 Get3Dauto(dist3D, topologicaldistance, numAtoms, mol, res);
254 }
255 }
256 } // namespace Descriptors
257 } // namespace RDKit
258