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