1 //  $Id$
2 //
3 //   Copyright (C) 2005-2013 Greg Landrum and Rational Discovery LLC
4 //
5 //   @@ All Rights Reserved @@
6 //  This file is part of the RDKit.
7 //  The contents are covered by the terms of the BSD license
8 //  which is included in the file license.txt, found at the root
9 //  of the RDKit source tree.
10 //
11 #include <RDGeneral/test.h>
12 #include "UniformGrid3D.h"
13 #include <DataStructs/DiscreteValueVect.h>
14 #include "point.h"
15 #include <RDGeneral/Invariant.h>
16 #include <RDGeneral/utils.h>
17 #include "GridUtils.h"
18 #include <iostream>
19 #include <sstream>
20 #include <fstream>
21 #include <ios>
22 #include <cstdlib>
23 
24 using namespace RDGeom;
25 using namespace RDKit;
26 
testUniformGrid1()27 void testUniformGrid1() {
28   UniformGrid3D grd(6.0, 5.0, 4.0);
29   CHECK_INVARIANT(grd.getSize() == 960, "");
30   CHECK_INVARIANT(RDKit::feq(grd.getSpacing(), .5), "");
31   CHECK_INVARIANT(grd.getNumX() == 12, "");
32   CHECK_INVARIANT(grd.getNumY() == 10, "");
33   CHECK_INVARIANT(grd.getNumZ() == 8, "");
34 
35   grd.setSphereOccupancy(Point3D(0.0, 0.0, 0.0), 1.5, 0.25);
36   CHECK_INVARIANT(grd.getOccupancyVect()->getTotalVal() == 523, "");
37   // writeGridToFile(grd, "junk.grd");
38 
39   UniformGrid3D grd2(grd);
40   CHECK_INVARIANT(grd2.getSize() == 960, "");
41   CHECK_INVARIANT(RDKit::feq(grd2.getSpacing(), .5), "");
42   CHECK_INVARIANT(grd2.getNumX() == 12, "");
43   CHECK_INVARIANT(grd2.getNumY() == 10, "");
44   CHECK_INVARIANT(grd2.getNumZ() == 8, "");
45   CHECK_INVARIANT(grd2.getOccupancyVect()->getTotalVal() == 523, "");
46 
47   // make sure the data are actually decoupled:
48   grd.setSphereOccupancy(Point3D(1.0, 1.0, 0.0), 1.5, 0.25);
49   CHECK_INVARIANT(grd.getOccupancyVect()->getTotalVal() > 523, "");
50   CHECK_INVARIANT(grd2.getOccupancyVect()->getTotalVal() == 523, "");
51 }
52 
testUniformGrid2()53 void testUniformGrid2() {
54   // test distance metrics:
55   UniformGrid3D grd(10.0, 10.0, 10.0);
56   grd.setSphereOccupancy(Point3D(-2.0, -2.0, 0.0), 1.5, 0.25);
57   grd.setSphereOccupancy(Point3D(-2.0, 2.0, 0.0), 1.5, 0.25);
58   grd.setSphereOccupancy(Point3D(2.0, -2.0, 0.0), 1.5, 0.25);
59   grd.setSphereOccupancy(Point3D(2.0, 2.0, 0.0), 1.5, 0.25);
60   writeGridToFile(grd, "junk.grd");
61   double dist = tanimotoDistance(grd, grd);
62   CHECK_INVARIANT(RDKit::feq(dist, 0.0), "");
63   dist = tverskyIndex(grd, grd, 1.0, 1.0);
64   CHECK_INVARIANT(RDKit::feq(dist, 1.0), "");
65   dist = tverskyIndex(grd, grd, 1.0, 0.0);
66   CHECK_INVARIANT(RDKit::feq(dist, 1.0), "");
67   dist = tverskyIndex(grd, grd, 0.0, 1.0);
68   CHECK_INVARIANT(RDKit::feq(dist, 1.0), "");
69   dist = tverskyIndex(grd, grd, 0.25, 0.75);
70   CHECK_INVARIANT(RDKit::feq(dist, 1.0), "");
71   dist = tverskyIndex(grd, grd, 0.75, 0.25);
72   CHECK_INVARIANT(RDKit::feq(dist, 1.0), "");
73 
74   UniformGrid3D grd2(10.0, 10.0, 10.0);
75   grd2.setSphereOccupancy(Point3D(-2.0, -2.0, 0.0), 1.5, 0.25);
76   grd2.setSphereOccupancy(Point3D(-2.0, 2.0, 0.0), 1.5, 0.25);
77   grd2.setSphereOccupancy(Point3D(2.0, -2.0, 0.0), 1.5, 0.25);
78   dist = tanimotoDistance(grd, grd2);
79   CHECK_INVARIANT(RDKit::feq(dist, 0.25), "");
80   dist = tverskyIndex(grd, grd2, 1.0, 1.0);
81   CHECK_INVARIANT(RDKit::feq(dist, 0.75), "");
82   dist = tverskyIndex(grd, grd2, 1.0, 0.0);
83   CHECK_INVARIANT(RDKit::feq(dist, 0.75), "");
84   dist = tverskyIndex(grd, grd2, 0.0, 1.0);
85   CHECK_INVARIANT(RDKit::feq(dist, 1.0), "");
86   dist = tverskyIndex(grd, grd2, 0.25, 0.75);
87   CHECK_INVARIANT(RDKit::feq(dist, 0.923), "");
88   dist = tverskyIndex(grd, grd2, 0.75, 0.25);
89   CHECK_INVARIANT(RDKit::feq(dist, 0.8), "");
90   dist = protrudeDistance(grd, grd2);
91   CHECK_INVARIANT(RDKit::feq(dist, 0.25), "");
92   dist = protrudeDistance(grd2, grd);
93   CHECK_INVARIANT(RDKit::feq(dist, 0.0), "");
94 
95   UniformGrid3D grd3(10.0, 10.0, 10.0);
96   grd3.setSphereOccupancy(Point3D(-2.0, -2.0, 0.0), 1.5, 0.25);
97   grd3.setSphereOccupancy(Point3D(-2.0, 2.0, 0.0), 1.5, 0.25);
98   dist = tanimotoDistance(grd, grd3);
99   CHECK_INVARIANT(RDKit::feq(dist, 0.5), "");
100   dist = tverskyIndex(grd, grd3, 1.0, 1.0);
101   CHECK_INVARIANT(RDKit::feq(dist, 0.5), "");
102   dist = tverskyIndex(grd, grd3, 1.0, 0.0);
103   CHECK_INVARIANT(RDKit::feq(dist, 0.5), "");
104   dist = tverskyIndex(grd, grd3, 0.0, 1.0);
105   CHECK_INVARIANT(RDKit::feq(dist, 1.0), "");
106   dist = tverskyIndex(grd, grd3, 0.25, 0.75);
107   CHECK_INVARIANT(RDKit::feq(dist, 0.8), "");
108   dist = tverskyIndex(grd, grd3, 0.75, 0.25);
109   CHECK_INVARIANT(RDKit::feq(dist, 0.5714), "");
110   dist = protrudeDistance(grd, grd3);
111   CHECK_INVARIANT(RDKit::feq(dist, 0.5), "");
112   dist = protrudeDistance(grd3, grd);
113   CHECK_INVARIANT(RDKit::feq(dist, 0.0), "");
114 
115   UniformGrid3D grd4(10.0, 10.0, 10.0);
116   grd4.setSphereOccupancy(Point3D(2.0, 2.0, 0.0), 1.5, 0.25);
117   grd4.setSphereOccupancy(Point3D(2.0, -2.0, 0.0), 1.5, 0.25);
118   dist = tanimotoDistance(grd3, grd4);
119   CHECK_INVARIANT(RDKit::feq(dist, 1.0), "");
120   dist = tverskyIndex(grd3, grd4, 1.0, 1.0);
121   CHECK_INVARIANT(RDKit::feq(dist, 0.0), "");
122   dist = tverskyIndex(grd3, grd4, 1.0, 0.0);
123   CHECK_INVARIANT(RDKit::feq(dist, 0.0), "");
124   dist = tverskyIndex(grd3, grd4, 0.0, 1.0);
125   CHECK_INVARIANT(RDKit::feq(dist, 0.0), "");
126   dist = tverskyIndex(grd3, grd4, 0.25, 0.75);
127   CHECK_INVARIANT(RDKit::feq(dist, 0.0), "");
128   dist = tverskyIndex(grd3, grd4, 0.75, 0.25);
129   CHECK_INVARIANT(RDKit::feq(dist, 0.0), "");
130 
131   UniformGrid3D grd5(10.0, 10.0, 10.0);
132   grd5.setSphereOccupancy(Point3D(-2.0, -2.0, 0.0), 1.5, 0.25);
133   dist = tanimotoDistance(grd, grd5);
134   CHECK_INVARIANT(RDKit::feq(dist, 0.75), "");
135   dist = tverskyIndex(grd, grd5, 1.0, 1.0);
136   CHECK_INVARIANT(RDKit::feq(dist, 0.25), "");
137   dist = tverskyIndex(grd, grd5, 1.0, 0.0);
138   CHECK_INVARIANT(RDKit::feq(dist, 0.25), "");
139   dist = tverskyIndex(grd, grd5, 0.0, 1.0);
140   CHECK_INVARIANT(RDKit::feq(dist, 1.0), "");
141   dist = tverskyIndex(grd, grd5, 0.25, 0.75);
142   CHECK_INVARIANT(RDKit::feq(dist, 0.5714), "");
143   dist = tverskyIndex(grd, grd5, 0.75, 0.25);
144   CHECK_INVARIANT(RDKit::feq(dist, 0.3077), "");
145   dist = protrudeDistance(grd, grd5);
146   CHECK_INVARIANT(RDKit::feq(dist, 0.75), "");
147   dist = protrudeDistance(grd5, grd);
148   CHECK_INVARIANT(RDKit::feq(dist, 0.00), "");
149 }
150 
testUniformGridPickling()151 void testUniformGridPickling() {
152   {
153     // test tanimoto distance
154     UniformGrid3D grd(10.0, 10.0, 10.0);
155     grd.setSphereOccupancy(Point3D(-2.0, -2.0, 0.0), 1.5, 0.25);
156     grd.setSphereOccupancy(Point3D(-2.0, 2.0, 0.0), 1.5, 0.25);
157     grd.setSphereOccupancy(Point3D(2.0, -2.0, 0.0), 1.5, 0.25);
158     grd.setSphereOccupancy(Point3D(2.0, 2.0, 0.0), 1.5, 0.25);
159     UniformGrid3D grd2(grd.toString());
160     double dist = tanimotoDistance(grd, grd2);
161     CHECK_INVARIANT(RDKit::feq(dist, 0.0), "");
162   }
163 
164   {
165     std::string dirName = getenv("RDBASE");
166     dirName += "/Code/Geometry/testData/";
167     std::string pklName = dirName + "grid1.bin";
168     std::ifstream inS;
169     inS.open(pklName.c_str(), std::ios_base::binary);
170     unsigned int length;
171     inS >> length;
172     auto *buff = new char[length];
173     unsigned int nRead = 0;
174     while (nRead < length) {
175       nRead += inS.readsome(buff + nRead, length - nRead);
176     }
177     inS.close();
178     std::string pkl(buff, length);
179     delete[] buff;
180     UniformGrid3D grd(pkl);
181 
182     UniformGrid3D grd2(10.0, 10.0, 10.0);
183     grd2.setSphereOccupancy(Point3D(-2.0, -2.0, 0.0), 1.5, 0.25);
184     grd2.setSphereOccupancy(Point3D(-2.0, 2.0, 0.0), 1.5, 0.25);
185     grd2.setSphereOccupancy(Point3D(2.0, -2.0, 0.0), 1.5, 0.25);
186     grd2.setSphereOccupancy(Point3D(2.0, 2.0, 0.0), 1.5, 0.25);
187 
188     std::string pkl2 = grd2.toString();
189     TEST_ASSERT(pkl2.length() == pkl.length());
190     TEST_ASSERT(pkl2 == pkl);
191 
192     TEST_ASSERT(grd.getSize() == grd2.getSize());
193     TEST_ASSERT(grd.getNumX() == grd2.getNumX());
194     TEST_ASSERT(grd.getNumY() == grd2.getNumY());
195     TEST_ASSERT(grd.getNumZ() == grd2.getNumZ());
196     TEST_ASSERT(grd.compareParams(grd2));
197     double dist = tanimotoDistance(grd, grd2);
198     TEST_ASSERT(RDKit::feq(dist, 0.0));
199   }
200 }
201 
testUniformGridOps()202 void testUniformGridOps() {
203   UniformGrid3D grd(10.0, 10.0, 10.0);
204   grd.setSphereOccupancy(Point3D(-2.0, -2.0, 0.0), 1.0, 0.25);
205   grd.setSphereOccupancy(Point3D(-2.0, 2.0, 0.0), 1.0, 0.25);
206 
207   UniformGrid3D grd2(10.0, 10.0, 10.0);
208   grd2.setSphereOccupancy(Point3D(2.0, -2.0, 0.0), 1.0, 0.25);
209   grd2.setSphereOccupancy(Point3D(2.0, 2.0, 0.0), 1.0, 0.25);
210 
211   double dist = tanimotoDistance(grd, grd2);
212   CHECK_INVARIANT(RDKit::feq(dist, 1.0), "");
213 
214   UniformGrid3D grd3(grd);
215   grd3 |= grd2;
216 
217   dist = tanimotoDistance(grd3, grd);
218   CHECK_INVARIANT(RDKit::feq(dist, 0.5), "");
219   dist = tanimotoDistance(grd3, grd2);
220   CHECK_INVARIANT(RDKit::feq(dist, 0.5), "");
221 
222   UniformGrid3D grd4(10.0, 10.0, 10.0);
223   grd4.setSphereOccupancy(Point3D(-2.0, -2.0, 0.0), 1.0, 0.25);
224   grd4.setSphereOccupancy(Point3D(-2.0, 2.0, 0.0), 1.0, 0.25);
225   grd4.setSphereOccupancy(Point3D(2.0, -2.0, 0.0), 1.0, 0.25);
226 
227   UniformGrid3D grd5(grd4);
228   grd5 &= grd2;
229 
230   dist = tanimotoDistance(grd5, grd);
231   CHECK_INVARIANT(RDKit::feq(dist, 1.0), "");
232   dist = tanimotoDistance(grd5, grd2);
233   CHECK_INVARIANT(RDKit::feq(dist, 0.5), "");
234 }
235 
testUniformGridIndexing()236 void testUniformGridIndexing() {
237   // test distance metrics:
238   UniformGrid3D grd(5.0, 5.0, 5.0);
239 
240   {
241     unsigned int xi = 3, yi = 3, zi = 3;
242     unsigned int idx = grd.getGridIndex(xi, yi, zi);
243     unsigned int nxi, nyi, nzi;
244     grd.getGridIndices(idx, nxi, nyi, nzi);
245     TEST_ASSERT(nxi == xi);
246     TEST_ASSERT(nyi == yi);
247     TEST_ASSERT(nzi == zi);
248   }
249   {
250     unsigned int xi = 3, yi = 3, zi = 5;
251     unsigned int idx = grd.getGridIndex(xi, yi, zi);
252     unsigned int nxi, nyi, nzi;
253     grd.getGridIndices(idx, nxi, nyi, nzi);
254     TEST_ASSERT(nxi == xi);
255     TEST_ASSERT(nyi == yi);
256     TEST_ASSERT(nzi == zi);
257   }
258   {
259     unsigned int xi = 3, yi = 6, zi = 3;
260     unsigned int idx = grd.getGridIndex(xi, yi, zi);
261     unsigned int nxi, nyi, nzi;
262     grd.getGridIndices(idx, nxi, nyi, nzi);
263     TEST_ASSERT(nxi == xi);
264     TEST_ASSERT(nyi == yi);
265     TEST_ASSERT(nzi == zi);
266   }
267   {
268     unsigned int xi = 0, yi = 0, zi = 0;
269     unsigned int idx = grd.getGridIndex(xi, yi, zi);
270     unsigned int nxi, nyi, nzi;
271     grd.getGridIndices(idx, nxi, nyi, nzi);
272     TEST_ASSERT(nxi == xi);
273     TEST_ASSERT(nyi == yi);
274     TEST_ASSERT(nzi == zi);
275   }
276   {
277     unsigned int xi = 8, yi = 2, zi = 1;
278     unsigned int idx = grd.getGridIndex(xi, yi, zi);
279     unsigned int nxi, nyi, nzi;
280     grd.getGridIndices(idx, nxi, nyi, nzi);
281     TEST_ASSERT(nxi == xi);
282     TEST_ASSERT(nyi == yi);
283     TEST_ASSERT(nzi == zi);
284   }
285 }
286 
main()287 int main() {
288   std::cout << "***********************************************************\n";
289   std::cout << "Testing Grid\n";
290 
291   std::cout << "\t---------------------------------\n";
292   std::cout << "\t testUniformGrid1 \n\n";
293   testUniformGrid1();
294 
295   std::cout << "\t---------------------------------\n";
296   std::cout << "\t testUniformGrid2 \n\n";
297   testUniformGrid2();
298 
299   std::cout << "\t---------------------------------\n";
300   std::cout << "\t testUniformGridPickling \n\n";
301   testUniformGridPickling();
302 
303   std::cout << "\t---------------------------------\n";
304   std::cout << "\t testGridOps \n\n";
305   testUniformGridOps();
306 
307   std::cout << "\t---------------------------------\n";
308   std::cout << "\t testUniformGridIndexing \n\n";
309   testUniformGridIndexing();
310 
311   return 0;
312 }
313