2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2013,2014, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 * Tests selection neighborhood searching.
40 * Increase coverage of these tests for different corner cases: other PBC cases
41 * than full 3D, large cutoffs (larger than half the box size), etc.
42 * At least some of these probably don't work correctly.
44 * \author Teemu Murtola <teemu.murtola@gmail.com>
45 * \ingroup module_selection
47 #include <gtest/gtest.h>
55 #include "gromacs/legacyheaders/pbc.h"
56 #include "gromacs/legacyheaders/vec.h"
58 #include "gromacs/selection/nbsearch.h"
59 #include "gromacs/random/random.h"
60 #include "gromacs/utility/smalloc.h"
62 #include "testutils/testasserts.h"
67 /********************************************************************
68 * NeighborhoodSearchTestData
71 class NeighborhoodSearchTestData
76 explicit TestPosition(const rvec x)
77 : refMinDist(0.0), refNearestPoint(-1)
79 copy_rvec(x, this->x);
85 std::set<int> refPairs;
87 typedef std::vector<TestPosition> TestPositionList;
89 NeighborhoodSearchTestData(int seed, real cutoff);
90 ~NeighborhoodSearchTestData();
92 gmx::AnalysisNeighborhoodPositions refPositions() const
94 return gmx::AnalysisNeighborhoodPositions(refPos_, refPosCount_);
96 gmx::AnalysisNeighborhoodPositions testPositions() const
100 snew(testPos_, testPositions_.size());
101 for (size_t i = 0; i < testPositions_.size(); ++i)
103 copy_rvec(testPositions_[i].x, testPos_[i]);
106 return gmx::AnalysisNeighborhoodPositions(testPos_,
107 testPositions_.size());
109 gmx::AnalysisNeighborhoodPositions testPosition(int index) const
111 return testPositions().selectSingleFromArray(index);
114 void addTestPosition(const rvec x)
116 GMX_RELEASE_ASSERT(testPos_ == NULL,
117 "Cannot add positions after testPositions() call");
118 testPositions_.push_back(TestPosition(x));
120 void generateRandomPosition(rvec x);
121 void generateRandomRefPositions(int count);
122 void generateRandomTestPositions(int count);
123 void computeReferences(t_pbc *pbc);
131 TestPositionList testPositions_;
134 mutable rvec *testPos_;
137 NeighborhoodSearchTestData::NeighborhoodSearchTestData(int seed, real cutoff)
138 : rng_(NULL), cutoff_(cutoff), refPosCount_(0), refPos_(NULL), testPos_(NULL)
140 // TODO: Handle errors.
141 rng_ = gmx_rng_init(seed);
143 set_pbc(&pbc_, epbcNONE, box_);
146 NeighborhoodSearchTestData::~NeighborhoodSearchTestData()
150 gmx_rng_destroy(rng_);
156 void NeighborhoodSearchTestData::generateRandomPosition(rvec x)
159 fx[XX] = gmx_rng_uniform_real(rng_);
160 fx[YY] = gmx_rng_uniform_real(rng_);
161 fx[ZZ] = gmx_rng_uniform_real(rng_);
163 // Add a small displacement to allow positions outside the box
164 x[XX] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
165 x[YY] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
166 x[ZZ] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
169 void NeighborhoodSearchTestData::generateRandomRefPositions(int count)
171 refPosCount_ = count;
172 snew(refPos_, refPosCount_);
173 for (int i = 0; i < refPosCount_; ++i)
175 generateRandomPosition(refPos_[i]);
179 void NeighborhoodSearchTestData::generateRandomTestPositions(int count)
181 testPositions_.reserve(count);
182 for (int i = 0; i < count; ++i)
185 generateRandomPosition(x);
190 void NeighborhoodSearchTestData::computeReferences(t_pbc *pbc)
192 real cutoff = cutoff_;
195 cutoff = std::numeric_limits<real>::max();
197 TestPositionList::iterator i;
198 for (i = testPositions_.begin(); i != testPositions_.end(); ++i)
200 i->refMinDist = cutoff;
201 i->refNearestPoint = -1;
203 for (int j = 0; j < refPosCount_; ++j)
208 pbc_dx(pbc, i->x, refPos_[j], dx);
212 rvec_sub(i->x, refPos_[j], dx);
214 const real dist = norm(dx);
215 if (dist < i->refMinDist)
217 i->refMinDist = dist;
218 i->refNearestPoint = j;
222 i->refPairs.insert(j);
228 /********************************************************************
229 * NeighborhoodSearchTest
232 class NeighborhoodSearchTest : public ::testing::Test
235 void testIsWithin(gmx::AnalysisNeighborhoodSearch *search,
236 const NeighborhoodSearchTestData &data);
237 void testMinimumDistance(gmx::AnalysisNeighborhoodSearch *search,
238 const NeighborhoodSearchTestData &data);
239 void testNearestPoint(gmx::AnalysisNeighborhoodSearch *search,
240 const NeighborhoodSearchTestData &data);
241 void testPairSearch(gmx::AnalysisNeighborhoodSearch *search,
242 const NeighborhoodSearchTestData &data);
244 gmx::AnalysisNeighborhood nb_;
247 void NeighborhoodSearchTest::testIsWithin(
248 gmx::AnalysisNeighborhoodSearch *search,
249 const NeighborhoodSearchTestData &data)
251 NeighborhoodSearchTestData::TestPositionList::const_iterator i;
252 for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
254 const bool bWithin = (i->refMinDist <= data.cutoff_);
255 EXPECT_EQ(bWithin, search->isWithin(i->x))
256 << "Distance is " << i->refMinDist;
260 void NeighborhoodSearchTest::testMinimumDistance(
261 gmx::AnalysisNeighborhoodSearch *search,
262 const NeighborhoodSearchTestData &data)
264 NeighborhoodSearchTestData::TestPositionList::const_iterator i;
265 for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
267 const real refDist = i->refMinDist;
268 EXPECT_REAL_EQ_TOL(refDist, search->minimumDistance(i->x),
269 gmx::test::ulpTolerance(20));
273 void NeighborhoodSearchTest::testNearestPoint(
274 gmx::AnalysisNeighborhoodSearch *search,
275 const NeighborhoodSearchTestData &data)
277 NeighborhoodSearchTestData::TestPositionList::const_iterator i;
278 for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
280 const gmx::AnalysisNeighborhoodPair pair = search->nearestPoint(i->x);
283 EXPECT_EQ(i->refNearestPoint, pair.refIndex());
284 EXPECT_EQ(0, pair.testIndex());
288 EXPECT_EQ(i->refNearestPoint, -1);
293 void NeighborhoodSearchTest::testPairSearch(
294 gmx::AnalysisNeighborhoodSearch *search,
295 const NeighborhoodSearchTestData &data)
297 NeighborhoodSearchTestData::TestPositionList::const_iterator i;
298 for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
300 std::set<int> checkSet = i->refPairs;
301 gmx::AnalysisNeighborhoodPairSearch pairSearch =
302 search->startPairSearch(i->x);
303 gmx::AnalysisNeighborhoodPair pair;
304 while (pairSearch.findNextPair(&pair))
306 EXPECT_EQ(0, pair.testIndex());
307 if (checkSet.erase(pair.refIndex()) == 0)
309 // TODO: Check whether the same pair was returned more than
310 // once and give a better error message if so.
312 << "Expected: Position " << pair.refIndex()
313 << " is within cutoff.\n"
314 << " Actual: It is not.";
317 EXPECT_TRUE(checkSet.empty()) << "Some positions were not returned by the pair search.";
321 /********************************************************************
322 * Test data generation
325 class TrivialTestData
328 static const NeighborhoodSearchTestData &get()
330 static TrivialTestData singleton;
331 return singleton.data_;
334 TrivialTestData() : data_(12345, 1.0)
336 data_.box_[XX][XX] = 5.0;
337 data_.box_[YY][YY] = 5.0;
338 data_.box_[ZZ][ZZ] = 5.0;
339 data_.generateRandomRefPositions(10);
340 data_.generateRandomTestPositions(5);
341 set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
342 data_.computeReferences(&data_.pbc_);
346 NeighborhoodSearchTestData data_;
349 class RandomBoxFullPBCData
352 static const NeighborhoodSearchTestData &get()
354 static RandomBoxFullPBCData singleton;
355 return singleton.data_;
358 RandomBoxFullPBCData() : data_(12345, 1.0)
360 data_.box_[XX][XX] = 10.0;
361 data_.box_[YY][YY] = 5.0;
362 data_.box_[ZZ][ZZ] = 7.0;
363 // TODO: Consider whether manually picking some positions would give better
365 data_.generateRandomRefPositions(1000);
366 data_.generateRandomTestPositions(100);
367 set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
368 data_.computeReferences(&data_.pbc_);
372 NeighborhoodSearchTestData data_;
375 class RandomTriclinicFullPBCData
378 static const NeighborhoodSearchTestData &get()
380 static RandomTriclinicFullPBCData singleton;
381 return singleton.data_;
384 RandomTriclinicFullPBCData() : data_(12345, 1.0)
386 data_.box_[XX][XX] = 5.0;
387 data_.box_[YY][XX] = 2.5;
388 data_.box_[YY][YY] = 2.5*sqrt(3.0);
389 data_.box_[ZZ][XX] = 2.5;
390 data_.box_[ZZ][YY] = 2.5*sqrt(1.0/3.0);
391 data_.box_[ZZ][ZZ] = 5.0*sqrt(2.0/3.0);
392 // TODO: Consider whether manually picking some positions would give better
394 data_.generateRandomRefPositions(1000);
395 data_.generateRandomTestPositions(100);
396 set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
397 data_.computeReferences(&data_.pbc_);
401 NeighborhoodSearchTestData data_;
404 class RandomBox2DPBCData
407 static const NeighborhoodSearchTestData &get()
409 static RandomBox2DPBCData singleton;
410 return singleton.data_;
413 RandomBox2DPBCData() : data_(12345, 1.0)
415 data_.box_[XX][XX] = 10.0;
416 data_.box_[YY][YY] = 7.0;
417 data_.box_[ZZ][ZZ] = 5.0;
418 // TODO: Consider whether manually picking some positions would give better
420 data_.generateRandomRefPositions(1000);
421 data_.generateRandomTestPositions(100);
422 set_pbc(&data_.pbc_, epbcXY, data_.box_);
423 data_.computeReferences(&data_.pbc_);
427 NeighborhoodSearchTestData data_;
430 /********************************************************************
434 TEST_F(NeighborhoodSearchTest, SimpleSearch)
436 const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
438 nb_.setCutoff(data.cutoff_);
439 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Simple);
440 gmx::AnalysisNeighborhoodSearch search =
441 nb_.initSearch(&data.pbc_, data.refPositions());
442 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
444 testIsWithin(&search, data);
445 testMinimumDistance(&search, data);
446 testNearestPoint(&search, data);
447 testPairSearch(&search, data);
450 TEST_F(NeighborhoodSearchTest, GridSearchBox)
452 const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
454 nb_.setCutoff(data.cutoff_);
455 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
456 gmx::AnalysisNeighborhoodSearch search =
457 nb_.initSearch(&data.pbc_, data.refPositions());
458 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
460 testIsWithin(&search, data);
461 testMinimumDistance(&search, data);
462 testNearestPoint(&search, data);
463 testPairSearch(&search, data);
466 TEST_F(NeighborhoodSearchTest, GridSearchTriclinic)
468 const NeighborhoodSearchTestData &data = RandomTriclinicFullPBCData::get();
470 nb_.setCutoff(data.cutoff_);
471 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
472 gmx::AnalysisNeighborhoodSearch search =
473 nb_.initSearch(&data.pbc_, data.refPositions());
474 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
476 testPairSearch(&search, data);
479 TEST_F(NeighborhoodSearchTest, GridSearch2DPBC)
481 const NeighborhoodSearchTestData &data = RandomBox2DPBCData::get();
483 nb_.setCutoff(data.cutoff_);
484 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
485 gmx::AnalysisNeighborhoodSearch search =
486 nb_.initSearch(&data.pbc_, data.refPositions());
487 // Currently, grid searching not supported with 2D PBC.
488 //ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
490 testIsWithin(&search, data);
491 testMinimumDistance(&search, data);
492 testNearestPoint(&search, data);
493 testPairSearch(&search, data);
496 TEST_F(NeighborhoodSearchTest, HandlesConcurrentSearches)
498 const NeighborhoodSearchTestData &data = TrivialTestData::get();
500 nb_.setCutoff(data.cutoff_);
501 gmx::AnalysisNeighborhoodSearch search1 =
502 nb_.initSearch(&data.pbc_, data.refPositions());
503 gmx::AnalysisNeighborhoodSearch search2 =
504 nb_.initSearch(&data.pbc_, data.refPositions());
506 gmx::AnalysisNeighborhoodPairSearch pairSearch1 =
507 search1.startPairSearch(data.testPosition(0));
508 gmx::AnalysisNeighborhoodPairSearch pairSearch2 =
509 search1.startPairSearch(data.testPosition(1));
511 testPairSearch(&search2, data);
513 gmx::AnalysisNeighborhoodPair pair;
514 pairSearch1.findNextPair(&pair);
515 EXPECT_EQ(0, pair.testIndex());
516 EXPECT_TRUE(data.testPositions_[0].refPairs.count(pair.refIndex()) == 1);
518 pairSearch2.findNextPair(&pair);
519 EXPECT_EQ(1, pair.testIndex());
520 EXPECT_TRUE(data.testPositions_[1].refPairs.count(pair.refIndex()) == 1);
523 TEST_F(NeighborhoodSearchTest, HandlesSkippingPairs)
525 const NeighborhoodSearchTestData &data = TrivialTestData::get();
527 nb_.setCutoff(data.cutoff_);
528 gmx::AnalysisNeighborhoodSearch search =
529 nb_.initSearch(&data.pbc_, data.refPositions());
530 gmx::AnalysisNeighborhoodPairSearch pairSearch =
531 search.startPairSearch(data.testPositions());
532 gmx::AnalysisNeighborhoodPair pair;
533 // TODO: This test needs to be adjusted if the grid search gets optimized
534 // to loop over the test positions in cell order (first, the ordering
535 // assumption here breaks, and second, it then needs to be tested
536 // separately for simple and grid searches).
537 int currentIndex = 0;
538 while (pairSearch.findNextPair(&pair))
540 while (currentIndex < pair.testIndex())
544 EXPECT_EQ(currentIndex, pair.testIndex());
545 EXPECT_TRUE(data.testPositions_[currentIndex].refPairs.count(pair.refIndex()) == 1);
546 pairSearch.skipRemainingPairsForTestPosition();