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/math/vec.h"
56 #include "gromacs/pbcutil/pbc.h"
57 #include "gromacs/random/random.h"
58 #include "gromacs/selection/nbsearch.h"
59 #include "gromacs/utility/smalloc.h"
61 #include "testutils/testasserts.h"
66 /********************************************************************
67 * NeighborhoodSearchTestData
70 class NeighborhoodSearchTestData
75 TestPosition() : refMinDist(0.0), refNearestPoint(-1)
79 explicit TestPosition(const rvec x)
80 : refMinDist(0.0), refNearestPoint(-1)
82 copy_rvec(x, this->x);
88 std::set<int> refPairs;
90 typedef std::vector<TestPosition> TestPositionList;
92 NeighborhoodSearchTestData(int seed, real cutoff);
93 ~NeighborhoodSearchTestData();
95 gmx::AnalysisNeighborhoodPositions refPositions() const
97 return gmx::AnalysisNeighborhoodPositions(refPos_, refPosCount_);
99 gmx::AnalysisNeighborhoodPositions testPositions() const
101 if (testPos_ == NULL)
103 snew(testPos_, testPositions_.size());
104 for (size_t i = 0; i < testPositions_.size(); ++i)
106 copy_rvec(testPositions_[i].x, testPos_[i]);
109 return gmx::AnalysisNeighborhoodPositions(testPos_,
110 testPositions_.size());
112 gmx::AnalysisNeighborhoodPositions testPosition(int index) const
114 return testPositions().selectSingleFromArray(index);
117 void addTestPosition(const rvec x)
119 GMX_RELEASE_ASSERT(testPos_ == NULL,
120 "Cannot add positions after testPositions() call");
121 testPositions_.push_back(TestPosition(x));
123 void generateRandomPosition(rvec x);
124 void generateRandomRefPositions(int count);
125 void generateRandomTestPositions(int count);
126 void computeReferences(t_pbc *pbc);
134 TestPositionList testPositions_;
137 mutable rvec *testPos_;
140 NeighborhoodSearchTestData::NeighborhoodSearchTestData(int seed, real cutoff)
141 : rng_(NULL), cutoff_(cutoff), refPosCount_(0), refPos_(NULL), testPos_(NULL)
143 // TODO: Handle errors.
144 rng_ = gmx_rng_init(seed);
146 set_pbc(&pbc_, epbcNONE, box_);
149 NeighborhoodSearchTestData::~NeighborhoodSearchTestData()
153 gmx_rng_destroy(rng_);
159 void NeighborhoodSearchTestData::generateRandomPosition(rvec x)
162 fx[XX] = gmx_rng_uniform_real(rng_);
163 fx[YY] = gmx_rng_uniform_real(rng_);
164 fx[ZZ] = gmx_rng_uniform_real(rng_);
166 // Add a small displacement to allow positions outside the box
167 x[XX] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
168 x[YY] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
169 x[ZZ] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
172 void NeighborhoodSearchTestData::generateRandomRefPositions(int count)
174 refPosCount_ = count;
175 snew(refPos_, refPosCount_);
176 for (int i = 0; i < refPosCount_; ++i)
178 generateRandomPosition(refPos_[i]);
182 void NeighborhoodSearchTestData::generateRandomTestPositions(int count)
184 testPositions_.reserve(count);
185 for (int i = 0; i < count; ++i)
188 generateRandomPosition(x);
193 void NeighborhoodSearchTestData::computeReferences(t_pbc *pbc)
195 real cutoff = cutoff_;
198 cutoff = std::numeric_limits<real>::max();
200 TestPositionList::iterator i;
201 for (i = testPositions_.begin(); i != testPositions_.end(); ++i)
203 i->refMinDist = cutoff;
204 i->refNearestPoint = -1;
206 for (int j = 0; j < refPosCount_; ++j)
211 pbc_dx(pbc, i->x, refPos_[j], dx);
215 rvec_sub(i->x, refPos_[j], dx);
217 const real dist = norm(dx);
218 if (dist < i->refMinDist)
220 i->refMinDist = dist;
221 i->refNearestPoint = j;
225 i->refPairs.insert(j);
231 /********************************************************************
232 * NeighborhoodSearchTest
235 class NeighborhoodSearchTest : public ::testing::Test
238 void testIsWithin(gmx::AnalysisNeighborhoodSearch *search,
239 const NeighborhoodSearchTestData &data);
240 void testMinimumDistance(gmx::AnalysisNeighborhoodSearch *search,
241 const NeighborhoodSearchTestData &data);
242 void testNearestPoint(gmx::AnalysisNeighborhoodSearch *search,
243 const NeighborhoodSearchTestData &data);
244 void testPairSearch(gmx::AnalysisNeighborhoodSearch *search,
245 const NeighborhoodSearchTestData &data);
247 gmx::AnalysisNeighborhood nb_;
250 void NeighborhoodSearchTest::testIsWithin(
251 gmx::AnalysisNeighborhoodSearch *search,
252 const NeighborhoodSearchTestData &data)
254 NeighborhoodSearchTestData::TestPositionList::const_iterator i;
255 for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
257 const bool bWithin = (i->refMinDist <= data.cutoff_);
258 EXPECT_EQ(bWithin, search->isWithin(i->x))
259 << "Distance is " << i->refMinDist;
263 void NeighborhoodSearchTest::testMinimumDistance(
264 gmx::AnalysisNeighborhoodSearch *search,
265 const NeighborhoodSearchTestData &data)
267 NeighborhoodSearchTestData::TestPositionList::const_iterator i;
268 for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
270 const real refDist = i->refMinDist;
271 EXPECT_REAL_EQ_TOL(refDist, search->minimumDistance(i->x),
272 gmx::test::ulpTolerance(20));
276 void NeighborhoodSearchTest::testNearestPoint(
277 gmx::AnalysisNeighborhoodSearch *search,
278 const NeighborhoodSearchTestData &data)
280 NeighborhoodSearchTestData::TestPositionList::const_iterator i;
281 for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
283 const gmx::AnalysisNeighborhoodPair pair = search->nearestPoint(i->x);
286 EXPECT_EQ(i->refNearestPoint, pair.refIndex());
287 EXPECT_EQ(0, pair.testIndex());
291 EXPECT_EQ(i->refNearestPoint, -1);
296 void NeighborhoodSearchTest::testPairSearch(
297 gmx::AnalysisNeighborhoodSearch *search,
298 const NeighborhoodSearchTestData &data)
300 NeighborhoodSearchTestData::TestPositionList::const_iterator i;
301 for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
303 std::set<int> checkSet = i->refPairs;
304 gmx::AnalysisNeighborhoodPairSearch pairSearch =
305 search->startPairSearch(i->x);
306 gmx::AnalysisNeighborhoodPair pair;
307 while (pairSearch.findNextPair(&pair))
309 EXPECT_EQ(0, pair.testIndex());
310 if (checkSet.erase(pair.refIndex()) == 0)
312 // TODO: Check whether the same pair was returned more than
313 // once and give a better error message if so.
315 << "Expected: Position " << pair.refIndex()
316 << " is within cutoff.\n"
317 << " Actual: It is not.";
320 EXPECT_TRUE(checkSet.empty()) << "Some positions were not returned by the pair search.";
324 /********************************************************************
325 * Test data generation
328 class TrivialTestData
331 static const NeighborhoodSearchTestData &get()
333 static TrivialTestData singleton;
334 return singleton.data_;
337 TrivialTestData() : data_(12345, 1.0)
339 data_.box_[XX][XX] = 5.0;
340 data_.box_[YY][YY] = 5.0;
341 data_.box_[ZZ][ZZ] = 5.0;
342 data_.generateRandomRefPositions(10);
343 data_.generateRandomTestPositions(5);
344 set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
345 data_.computeReferences(&data_.pbc_);
349 NeighborhoodSearchTestData data_;
352 class RandomBoxFullPBCData
355 static const NeighborhoodSearchTestData &get()
357 static RandomBoxFullPBCData singleton;
358 return singleton.data_;
361 RandomBoxFullPBCData() : data_(12345, 1.0)
363 data_.box_[XX][XX] = 10.0;
364 data_.box_[YY][YY] = 5.0;
365 data_.box_[ZZ][ZZ] = 7.0;
366 // TODO: Consider whether manually picking some positions would give better
368 data_.generateRandomRefPositions(1000);
369 data_.generateRandomTestPositions(100);
370 set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
371 data_.computeReferences(&data_.pbc_);
375 NeighborhoodSearchTestData data_;
378 class RandomTriclinicFullPBCData
381 static const NeighborhoodSearchTestData &get()
383 static RandomTriclinicFullPBCData singleton;
384 return singleton.data_;
387 RandomTriclinicFullPBCData() : data_(12345, 1.0)
389 data_.box_[XX][XX] = 5.0;
390 data_.box_[YY][XX] = 2.5;
391 data_.box_[YY][YY] = 2.5*sqrt(3.0);
392 data_.box_[ZZ][XX] = 2.5;
393 data_.box_[ZZ][YY] = 2.5*sqrt(1.0/3.0);
394 data_.box_[ZZ][ZZ] = 5.0*sqrt(2.0/3.0);
395 // TODO: Consider whether manually picking some positions would give better
397 data_.generateRandomRefPositions(1000);
398 data_.generateRandomTestPositions(100);
399 set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
400 data_.computeReferences(&data_.pbc_);
404 NeighborhoodSearchTestData data_;
407 class RandomBox2DPBCData
410 static const NeighborhoodSearchTestData &get()
412 static RandomBox2DPBCData singleton;
413 return singleton.data_;
416 RandomBox2DPBCData() : data_(12345, 1.0)
418 data_.box_[XX][XX] = 10.0;
419 data_.box_[YY][YY] = 7.0;
420 data_.box_[ZZ][ZZ] = 5.0;
421 // TODO: Consider whether manually picking some positions would give better
423 data_.generateRandomRefPositions(1000);
424 data_.generateRandomTestPositions(100);
425 set_pbc(&data_.pbc_, epbcXY, data_.box_);
426 data_.computeReferences(&data_.pbc_);
430 NeighborhoodSearchTestData data_;
433 /********************************************************************
437 TEST_F(NeighborhoodSearchTest, SimpleSearch)
439 const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
441 nb_.setCutoff(data.cutoff_);
442 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Simple);
443 gmx::AnalysisNeighborhoodSearch search =
444 nb_.initSearch(&data.pbc_, data.refPositions());
445 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
447 testIsWithin(&search, data);
448 testMinimumDistance(&search, data);
449 testNearestPoint(&search, data);
450 testPairSearch(&search, data);
453 TEST_F(NeighborhoodSearchTest, GridSearchBox)
455 const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
457 nb_.setCutoff(data.cutoff_);
458 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
459 gmx::AnalysisNeighborhoodSearch search =
460 nb_.initSearch(&data.pbc_, data.refPositions());
461 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
463 testIsWithin(&search, data);
464 testMinimumDistance(&search, data);
465 testNearestPoint(&search, data);
466 testPairSearch(&search, data);
469 TEST_F(NeighborhoodSearchTest, GridSearchTriclinic)
471 const NeighborhoodSearchTestData &data = RandomTriclinicFullPBCData::get();
473 nb_.setCutoff(data.cutoff_);
474 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
475 gmx::AnalysisNeighborhoodSearch search =
476 nb_.initSearch(&data.pbc_, data.refPositions());
477 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
479 testPairSearch(&search, data);
482 TEST_F(NeighborhoodSearchTest, GridSearch2DPBC)
484 const NeighborhoodSearchTestData &data = RandomBox2DPBCData::get();
486 nb_.setCutoff(data.cutoff_);
487 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
488 gmx::AnalysisNeighborhoodSearch search =
489 nb_.initSearch(&data.pbc_, data.refPositions());
490 // Currently, grid searching not supported with 2D PBC.
491 //ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
493 testIsWithin(&search, data);
494 testMinimumDistance(&search, data);
495 testNearestPoint(&search, data);
496 testPairSearch(&search, data);
499 TEST_F(NeighborhoodSearchTest, HandlesConcurrentSearches)
501 const NeighborhoodSearchTestData &data = TrivialTestData::get();
503 nb_.setCutoff(data.cutoff_);
504 gmx::AnalysisNeighborhoodSearch search1 =
505 nb_.initSearch(&data.pbc_, data.refPositions());
506 gmx::AnalysisNeighborhoodSearch search2 =
507 nb_.initSearch(&data.pbc_, data.refPositions());
509 gmx::AnalysisNeighborhoodPairSearch pairSearch1 =
510 search1.startPairSearch(data.testPosition(0));
511 gmx::AnalysisNeighborhoodPairSearch pairSearch2 =
512 search1.startPairSearch(data.testPosition(1));
514 testPairSearch(&search2, data);
516 gmx::AnalysisNeighborhoodPair pair;
517 pairSearch1.findNextPair(&pair);
518 EXPECT_EQ(0, pair.testIndex());
519 EXPECT_TRUE(data.testPositions_[0].refPairs.count(pair.refIndex()) == 1);
521 pairSearch2.findNextPair(&pair);
522 EXPECT_EQ(1, pair.testIndex());
523 EXPECT_TRUE(data.testPositions_[1].refPairs.count(pair.refIndex()) == 1);
526 TEST_F(NeighborhoodSearchTest, HandlesSkippingPairs)
528 const NeighborhoodSearchTestData &data = TrivialTestData::get();
530 nb_.setCutoff(data.cutoff_);
531 gmx::AnalysisNeighborhoodSearch search =
532 nb_.initSearch(&data.pbc_, data.refPositions());
533 gmx::AnalysisNeighborhoodPairSearch pairSearch =
534 search.startPairSearch(data.testPositions());
535 gmx::AnalysisNeighborhoodPair pair;
536 // TODO: This test needs to be adjusted if the grid search gets optimized
537 // to loop over the test positions in cell order (first, the ordering
538 // assumption here breaks, and second, it then needs to be tested
539 // separately for simple and grid searches).
540 int currentIndex = 0;
541 while (pairSearch.findNextPair(&pair))
543 while (currentIndex < pair.testIndex())
547 EXPECT_EQ(currentIndex, pair.testIndex());
548 EXPECT_TRUE(data.testPositions_[currentIndex].refPairs.count(pair.refIndex()) == 1);
549 pairSearch.skipRemainingPairsForTestPosition();