2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
5 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
38 * Tests selection neighborhood searching.
41 * Increase coverage of these tests for different corner cases: other PBC cases
42 * than full 3D, large cutoffs (larger than half the box size), etc.
43 * At least some of these probably don't work correctly.
45 * \author Teemu Murtola <teemu.murtola@gmail.com>
46 * \ingroup module_selection
50 #include "gromacs/selection/nbsearch.h"
60 #include <gtest/gtest.h>
62 #include "gromacs/math/functions.h"
63 #include "gromacs/math/vec.h"
64 #include "gromacs/pbcutil/pbc.h"
65 #include "gromacs/random/threefry.h"
66 #include "gromacs/random/uniformrealdistribution.h"
67 #include "gromacs/topology/block.h"
68 #include "gromacs/utility/listoflists.h"
69 #include "gromacs/utility/smalloc.h"
70 #include "gromacs/utility/stringutil.h"
72 #include "testutils/testasserts.h"
78 /********************************************************************
79 * NeighborhoodSearchTestData
82 class NeighborhoodSearchTestData
87 RefPair(int refIndex, real distance) :
96 bool operator<(const RefPair& other) const { return refIndex < other.refIndex; }
100 // The variables below are state variables that are only used
101 // during the actual testing after creating a copy of the reference
102 // pair list, not as part of the reference data.
103 // Simpler to have just a single structure for both purposes.
111 explicit TestPosition(const rvec x) : refMinDist(0.0), refNearestPoint(-1)
113 copy_rvec(x, this->x);
119 std::vector<RefPair> refPairs;
122 typedef std::vector<TestPosition> TestPositionList;
124 NeighborhoodSearchTestData(uint64_t seed, real cutoff);
126 gmx::AnalysisNeighborhoodPositions refPositions() const
128 return gmx::AnalysisNeighborhoodPositions(refPos_);
130 gmx::AnalysisNeighborhoodPositions testPositions() const
132 if (testPos_.empty())
134 testPos_.reserve(testPositions_.size());
135 for (size_t i = 0; i < testPositions_.size(); ++i)
137 testPos_.emplace_back(testPositions_[i].x);
140 return gmx::AnalysisNeighborhoodPositions(testPos_);
142 gmx::AnalysisNeighborhoodPositions testPosition(int index) const
144 return testPositions().selectSingleFromArray(index);
147 void addTestPosition(const rvec x)
149 GMX_RELEASE_ASSERT(testPos_.empty(), "Cannot add positions after testPositions() call");
150 testPositions_.emplace_back(x);
152 gmx::RVec generateRandomPosition();
153 static std::vector<int> generateIndex(int count, uint64_t seed);
154 void generateRandomRefPositions(int count);
155 void generateRandomTestPositions(int count);
156 void useRefPositionsAsTestPositions();
157 void computeReferences(t_pbc* pbc) { computeReferencesInternal(pbc, false); }
158 void computeReferencesXY(t_pbc* pbc) { computeReferencesInternal(pbc, true); }
160 bool containsPair(int testIndex, const RefPair& pair) const
162 const std::vector<RefPair>& refPairs = testPositions_[testIndex].refPairs;
163 std::vector<RefPair>::const_iterator foundRefPair =
164 std::lower_bound(refPairs.begin(), refPairs.end(), pair);
165 return !(foundRefPair == refPairs.end() || foundRefPair->refIndex != pair.refIndex);
168 // Return a tolerance that accounts for the magnitudes of the coordinates
169 // when doing subtractions, so that we set the ULP tolerance relative to the
170 // coordinate values rather than their difference.
171 // i.e., 10.0-9.9999999 will achieve a few ULP accuracy relative
172 // to 10.0, but a much larger error relative to the difference.
173 gmx::test::FloatingPointTolerance relativeTolerance() const
175 real magnitude = std::max(box_[XX][XX], std::max(box_[YY][YY], box_[ZZ][ZZ]));
176 return gmx::test::relativeToleranceAsUlp(magnitude, 4);
179 gmx::DefaultRandomEngine rng_;
184 std::vector<gmx::RVec> refPos_;
185 TestPositionList testPositions_;
188 void computeReferencesInternal(t_pbc* pbc, bool bXY);
190 mutable std::vector<gmx::RVec> testPos_;
193 //! Shorthand for a collection of reference pairs.
194 typedef std::vector<NeighborhoodSearchTestData::RefPair> RefPairList;
196 NeighborhoodSearchTestData::NeighborhoodSearchTestData(uint64_t seed, real cutoff) :
202 set_pbc(&pbc_, PbcType::No, box_);
205 gmx::RVec NeighborhoodSearchTestData::generateRandomPosition()
207 gmx::UniformRealDistribution<real> dist;
213 // Add a small displacement to allow positions outside the box
214 x[XX] += 0.2 * dist(rng_) - 0.1;
215 x[YY] += 0.2 * dist(rng_) - 0.1;
216 x[ZZ] += 0.2 * dist(rng_) - 0.1;
220 std::vector<int> NeighborhoodSearchTestData::generateIndex(int count, uint64_t seed)
222 gmx::DefaultRandomEngine rngIndex(seed);
223 gmx::UniformRealDistribution<real> dist;
224 std::vector<int> result;
226 for (int i = 0; i < count; ++i)
228 if (dist(rngIndex) > 0.5)
236 void NeighborhoodSearchTestData::generateRandomRefPositions(int count)
238 refPosCount_ = count;
239 refPos_.reserve(count);
240 for (int i = 0; i < count; ++i)
242 refPos_.push_back(generateRandomPosition());
246 void NeighborhoodSearchTestData::generateRandomTestPositions(int count)
248 testPositions_.reserve(count);
249 for (int i = 0; i < count; ++i)
251 addTestPosition(generateRandomPosition());
255 void NeighborhoodSearchTestData::useRefPositionsAsTestPositions()
257 testPositions_.reserve(refPosCount_);
258 for (const auto& refPos : refPos_)
260 addTestPosition(refPos);
264 void NeighborhoodSearchTestData::computeReferencesInternal(t_pbc* pbc, bool bXY)
266 real cutoff = cutoff_;
269 cutoff = std::numeric_limits<real>::max();
271 for (TestPosition& testPos : testPositions_)
273 testPos.refMinDist = cutoff;
274 testPos.refNearestPoint = -1;
275 testPos.refPairs.clear();
276 for (int j = 0; j < refPosCount_; ++j)
281 pbc_dx(pbc, testPos.x, refPos_[j], dx);
285 rvec_sub(testPos.x, refPos_[j], dx);
287 // TODO: This may not work intuitively for 2D with the third box
288 // vector not parallel to the Z axis, but neither does the actual
289 // neighborhood search.
290 const real dist = !bXY ? norm(dx) : std::hypot(dx[XX], dx[YY]);
291 if (dist < testPos.refMinDist)
293 testPos.refMinDist = dist;
294 testPos.refNearestPoint = j;
296 if (dist > 0 && dist <= cutoff)
298 RefPair pair(j, dist);
299 GMX_RELEASE_ASSERT(testPos.refPairs.empty() || testPos.refPairs.back() < pair,
300 "Reference pairs should be generated in sorted order");
301 testPos.refPairs.push_back(pair);
307 /********************************************************************
311 class ExclusionsHelper
314 static void markExcludedPairs(RefPairList* refPairs, int testIndex, const gmx::ListOfLists<int>* excls);
316 ExclusionsHelper(int refPosCount, int testPosCount);
318 void generateExclusions();
320 const gmx::ListOfLists<int>* exclusions() const { return &excls_; }
322 gmx::ArrayRef<const int> refPosIds() const
324 return gmx::makeArrayRef(exclusionIds_).subArray(0, refPosCount_);
326 gmx::ArrayRef<const int> testPosIds() const
328 return gmx::makeArrayRef(exclusionIds_).subArray(0, testPosCount_);
334 std::vector<int> exclusionIds_;
335 gmx::ListOfLists<int> excls_;
339 void ExclusionsHelper::markExcludedPairs(RefPairList* refPairs, int testIndex, const gmx::ListOfLists<int>* excls)
342 for (const int excludedIndex : (*excls)[testIndex])
344 NeighborhoodSearchTestData::RefPair searchPair(excludedIndex, 0.0);
345 RefPairList::iterator excludedRefPair =
346 std::lower_bound(refPairs->begin(), refPairs->end(), searchPair);
347 if (excludedRefPair != refPairs->end() && excludedRefPair->refIndex == excludedIndex)
349 excludedRefPair->bFound = true;
350 excludedRefPair->bExcluded = true;
356 ExclusionsHelper::ExclusionsHelper(int refPosCount, int testPosCount) :
357 refPosCount_(refPosCount),
358 testPosCount_(testPosCount)
360 // Generate an array of 0, 1, 2, ...
361 // TODO: Make the tests work also with non-trivial exclusion IDs,
363 exclusionIds_.resize(std::max(refPosCount, testPosCount), 1);
364 exclusionIds_[0] = 0;
365 std::partial_sum(exclusionIds_.begin(), exclusionIds_.end(), exclusionIds_.begin());
368 void ExclusionsHelper::generateExclusions()
370 // TODO: Consider a better set of test data, where the density of the
371 // particles would be higher, or where the exclusions would not be random,
372 // to make a higher percentage of the exclusions to actually be within the
374 for (int i = 0; i < testPosCount_; ++i)
376 excls_.pushBackListOfSize(20);
377 gmx::ArrayRef<int> exclusionsForAtom = excls_.back();
378 for (int j = 0; j < 20; ++j)
380 exclusionsForAtom[j] = i + j * 3;
385 /********************************************************************
386 * NeighborhoodSearchTest
389 class NeighborhoodSearchTest : public ::testing::Test
392 static void testIsWithin(gmx::AnalysisNeighborhoodSearch* search,
393 const NeighborhoodSearchTestData& data);
394 static void testMinimumDistance(gmx::AnalysisNeighborhoodSearch* search,
395 const NeighborhoodSearchTestData& data);
396 static void testNearestPoint(gmx::AnalysisNeighborhoodSearch* search,
397 const NeighborhoodSearchTestData& data);
398 static void testPairSearch(gmx::AnalysisNeighborhoodSearch* search,
399 const NeighborhoodSearchTestData& data);
400 static void testPairSearchIndexed(gmx::AnalysisNeighborhood* nb,
401 const NeighborhoodSearchTestData& data,
403 static void testPairSearchFull(gmx::AnalysisNeighborhoodSearch* search,
404 const NeighborhoodSearchTestData& data,
405 const gmx::AnalysisNeighborhoodPositions& pos,
406 const gmx::ListOfLists<int>* excls,
407 const gmx::ArrayRef<const int>& refIndices,
408 const gmx::ArrayRef<const int>& testIndices,
411 gmx::AnalysisNeighborhood nb_;
414 void NeighborhoodSearchTest::testIsWithin(gmx::AnalysisNeighborhoodSearch* search,
415 const NeighborhoodSearchTestData& data)
417 NeighborhoodSearchTestData::TestPositionList::const_iterator i;
418 for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
420 const bool bWithin = (i->refMinDist <= data.cutoff_);
421 EXPECT_EQ(bWithin, search->isWithin(i->x)) << "Distance is " << i->refMinDist;
425 void NeighborhoodSearchTest::testMinimumDistance(gmx::AnalysisNeighborhoodSearch* search,
426 const NeighborhoodSearchTestData& data)
428 NeighborhoodSearchTestData::TestPositionList::const_iterator i;
430 for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
432 const real refDist = i->refMinDist;
433 EXPECT_REAL_EQ_TOL(refDist, search->minimumDistance(i->x), data.relativeTolerance());
437 void NeighborhoodSearchTest::testNearestPoint(gmx::AnalysisNeighborhoodSearch* search,
438 const NeighborhoodSearchTestData& data)
440 NeighborhoodSearchTestData::TestPositionList::const_iterator i;
441 for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
443 const gmx::AnalysisNeighborhoodPair pair = search->nearestPoint(i->x);
446 EXPECT_EQ(i->refNearestPoint, pair.refIndex());
447 EXPECT_EQ(0, pair.testIndex());
448 EXPECT_REAL_EQ_TOL(i->refMinDist, std::sqrt(pair.distance2()), data.relativeTolerance());
452 EXPECT_EQ(i->refNearestPoint, -1);
457 //! Helper function for formatting test failure messages.
458 std::string formatVector(const rvec x)
460 return gmx::formatString("[%.3f, %.3f, %.3f]", x[XX], x[YY], x[ZZ]);
464 * Helper function to check that all expected pairs were found.
466 void checkAllPairsFound(const RefPairList& refPairs,
467 const std::vector<gmx::RVec>& refPos,
471 // This could be elegantly expressed with Google Mock matchers, but that
472 // has a significant effect on the runtime of the tests...
474 RefPairList::const_iterator first;
475 for (RefPairList::const_iterator i = refPairs.begin(); i != refPairs.end(); ++i)
485 ADD_FAILURE() << "Some pairs (" << count << "/" << refPairs.size() << ") "
486 << "within the cutoff were not found. First pair:\n"
487 << " Ref: " << first->refIndex << " at "
488 << formatVector(refPos[first->refIndex]) << "\n"
489 << "Test: " << testPosIndex << " at " << formatVector(testPos) << "\n"
490 << "Dist: " << first->distance;
494 void NeighborhoodSearchTest::testPairSearch(gmx::AnalysisNeighborhoodSearch* search,
495 const NeighborhoodSearchTestData& data)
497 testPairSearchFull(search, data, data.testPositions(), nullptr, {}, {}, false);
500 void NeighborhoodSearchTest::testPairSearchIndexed(gmx::AnalysisNeighborhood* nb,
501 const NeighborhoodSearchTestData& data,
504 std::vector<int> refIndices(data.generateIndex(data.refPos_.size(), seed++));
505 std::vector<int> testIndices(data.generateIndex(data.testPositions_.size(), seed++));
506 gmx::AnalysisNeighborhoodSearch search =
507 nb->initSearch(&data.pbc_, data.refPositions().indexed(refIndices));
508 testPairSearchFull(&search, data, data.testPositions(), nullptr, refIndices, testIndices, false);
511 void NeighborhoodSearchTest::testPairSearchFull(gmx::AnalysisNeighborhoodSearch* search,
512 const NeighborhoodSearchTestData& data,
513 const gmx::AnalysisNeighborhoodPositions& pos,
514 const gmx::ListOfLists<int>* excls,
515 const gmx::ArrayRef<const int>& refIndices,
516 const gmx::ArrayRef<const int>& testIndices,
519 std::map<int, RefPairList> refPairs;
520 // TODO: Some parts of this code do not work properly if pos does not
521 // initially contain all the test positions.
522 if (testIndices.empty())
524 for (size_t i = 0; i < data.testPositions_.size(); ++i)
526 refPairs[i] = data.testPositions_[i].refPairs;
531 for (int index : testIndices)
533 refPairs[index] = data.testPositions_[index].refPairs;
536 if (excls != nullptr)
538 GMX_RELEASE_ASSERT(!selfPairs, "Self-pairs testing not implemented with exclusions");
539 for (auto& entry : refPairs)
541 const int testIndex = entry.first;
542 ExclusionsHelper::markExcludedPairs(&entry.second, testIndex, excls);
545 if (!refIndices.empty())
547 GMX_RELEASE_ASSERT(!selfPairs, "Self-pairs testing not implemented with indexing");
548 for (auto& entry : refPairs)
550 for (auto& refPair : entry.second)
552 refPair.bIndexed = false;
554 for (int index : refIndices)
556 NeighborhoodSearchTestData::RefPair searchPair(index, 0.0);
557 auto refPair = std::lower_bound(entry.second.begin(), entry.second.end(), searchPair);
558 if (refPair != entry.second.end() && refPair->refIndex == index)
560 refPair->bIndexed = true;
563 for (auto& refPair : entry.second)
565 if (!refPair.bIndexed)
567 refPair.bFound = true;
573 gmx::AnalysisNeighborhoodPositions posCopy(pos);
574 if (!testIndices.empty())
576 posCopy.indexed(testIndices);
578 gmx::AnalysisNeighborhoodPairSearch pairSearch =
579 selfPairs ? search->startSelfPairSearch() : search->startPairSearch(posCopy);
580 gmx::AnalysisNeighborhoodPair pair;
581 while (pairSearch.findNextPair(&pair))
583 const int testIndex = (testIndices.empty() ? pair.testIndex() : testIndices[pair.testIndex()]);
584 const int refIndex = (refIndices.empty() ? pair.refIndex() : refIndices[pair.refIndex()]);
586 if (refPairs.count(testIndex) == 0)
588 ADD_FAILURE() << "Expected: No pairs are returned for test position " << testIndex << ".\n"
589 << " Actual: Pair with ref " << refIndex << " is returned.";
592 NeighborhoodSearchTestData::RefPair searchPair(refIndex, std::sqrt(pair.distance2()));
593 const auto foundRefPair =
594 std::lower_bound(refPairs[testIndex].begin(), refPairs[testIndex].end(), searchPair);
595 if (foundRefPair == refPairs[testIndex].end() || foundRefPair->refIndex != refIndex)
597 ADD_FAILURE() << "Expected: Pair (ref: " << refIndex << ", test: " << testIndex
598 << ") is not within the cutoff.\n"
599 << " Actual: It is returned.";
601 else if (foundRefPair->bExcluded)
603 ADD_FAILURE() << "Expected: Pair (ref: " << refIndex << ", test: " << testIndex
604 << ") is excluded from the search.\n"
605 << " Actual: It is returned.";
607 else if (!foundRefPair->bIndexed)
609 ADD_FAILURE() << "Expected: Pair (ref: " << refIndex << ", test: " << testIndex
610 << ") is not part of the indexed set.\n"
611 << " Actual: It is returned.";
613 else if (foundRefPair->bFound)
615 ADD_FAILURE() << "Expected: Pair (ref: " << refIndex << ", test: " << testIndex
616 << ") is returned only once.\n"
617 << " Actual: It is returned multiple times.";
622 foundRefPair->bFound = true;
624 EXPECT_REAL_EQ_TOL(foundRefPair->distance, searchPair.distance, data.relativeTolerance())
625 << "Distance computed by the neighborhood search does not match.";
628 searchPair = NeighborhoodSearchTestData::RefPair(testIndex, 0.0);
629 const auto otherRefPair = std::lower_bound(
630 refPairs[refIndex].begin(), refPairs[refIndex].end(), searchPair);
631 GMX_RELEASE_ASSERT(otherRefPair != refPairs[refIndex].end(),
632 "Precomputed reference data is not symmetric");
633 otherRefPair->bFound = true;
638 for (auto& entry : refPairs)
640 const int testIndex = entry.first;
641 checkAllPairsFound(entry.second, data.refPos_, testIndex, data.testPositions_[testIndex].x);
645 /********************************************************************
646 * Test data generation
649 class TrivialTestData
652 static const NeighborhoodSearchTestData& get()
654 static TrivialTestData singleton;
655 return singleton.data_;
658 TrivialTestData() : data_(12345, 1.0)
660 // Make the box so small we are virtually guaranteed to have
661 // several neighbors for the five test positions
662 data_.box_[XX][XX] = 3.0;
663 data_.box_[YY][YY] = 3.0;
664 data_.box_[ZZ][ZZ] = 3.0;
665 data_.generateRandomRefPositions(10);
666 data_.generateRandomTestPositions(5);
667 set_pbc(&data_.pbc_, PbcType::Xyz, data_.box_);
668 data_.computeReferences(&data_.pbc_);
672 NeighborhoodSearchTestData data_;
675 class TrivialSelfPairsTestData
678 static const NeighborhoodSearchTestData& get()
680 static TrivialSelfPairsTestData singleton;
681 return singleton.data_;
684 TrivialSelfPairsTestData() : data_(12345, 1.0)
686 data_.box_[XX][XX] = 3.0;
687 data_.box_[YY][YY] = 3.0;
688 data_.box_[ZZ][ZZ] = 3.0;
689 data_.generateRandomRefPositions(20);
690 data_.useRefPositionsAsTestPositions();
691 set_pbc(&data_.pbc_, PbcType::Xyz, data_.box_);
692 data_.computeReferences(&data_.pbc_);
696 NeighborhoodSearchTestData data_;
699 class TrivialNoPBCTestData
702 static const NeighborhoodSearchTestData& get()
704 static TrivialNoPBCTestData singleton;
705 return singleton.data_;
708 TrivialNoPBCTestData() : data_(12345, 1.0)
710 data_.generateRandomRefPositions(10);
711 data_.generateRandomTestPositions(5);
712 data_.computeReferences(nullptr);
716 NeighborhoodSearchTestData data_;
719 class RandomBoxFullPBCData
722 static const NeighborhoodSearchTestData& get()
724 static RandomBoxFullPBCData singleton;
725 return singleton.data_;
728 RandomBoxFullPBCData() : data_(12345, 1.0)
730 data_.box_[XX][XX] = 10.0;
731 data_.box_[YY][YY] = 5.0;
732 data_.box_[ZZ][ZZ] = 7.0;
733 // TODO: Consider whether manually picking some positions would give better
735 data_.generateRandomRefPositions(1000);
736 data_.generateRandomTestPositions(100);
737 set_pbc(&data_.pbc_, PbcType::Xyz, data_.box_);
738 data_.computeReferences(&data_.pbc_);
742 NeighborhoodSearchTestData data_;
745 class RandomBoxSelfPairsData
748 static const NeighborhoodSearchTestData& get()
750 static RandomBoxSelfPairsData singleton;
751 return singleton.data_;
754 RandomBoxSelfPairsData() : data_(12345, 1.0)
756 data_.box_[XX][XX] = 10.0;
757 data_.box_[YY][YY] = 5.0;
758 data_.box_[ZZ][ZZ] = 7.0;
759 data_.generateRandomRefPositions(1000);
760 data_.useRefPositionsAsTestPositions();
761 set_pbc(&data_.pbc_, PbcType::Xyz, data_.box_);
762 data_.computeReferences(&data_.pbc_);
766 NeighborhoodSearchTestData data_;
769 class RandomBoxXYFullPBCData
772 static const NeighborhoodSearchTestData& get()
774 static RandomBoxXYFullPBCData singleton;
775 return singleton.data_;
778 RandomBoxXYFullPBCData() : data_(54321, 1.0)
780 data_.box_[XX][XX] = 10.0;
781 data_.box_[YY][YY] = 5.0;
782 data_.box_[ZZ][ZZ] = 7.0;
783 // TODO: Consider whether manually picking some positions would give better
785 data_.generateRandomRefPositions(1000);
786 data_.generateRandomTestPositions(100);
787 set_pbc(&data_.pbc_, PbcType::Xyz, data_.box_);
788 data_.computeReferencesXY(&data_.pbc_);
792 NeighborhoodSearchTestData data_;
795 class RandomTriclinicFullPBCData
798 static const NeighborhoodSearchTestData& get()
800 static RandomTriclinicFullPBCData singleton;
801 return singleton.data_;
804 RandomTriclinicFullPBCData() : data_(12345, 1.0)
806 data_.box_[XX][XX] = 5.0;
807 data_.box_[YY][XX] = 2.5;
808 data_.box_[YY][YY] = 2.5 * std::sqrt(3.0);
809 data_.box_[ZZ][XX] = 2.5;
810 data_.box_[ZZ][YY] = 2.5 * std::sqrt(1.0 / 3.0);
811 data_.box_[ZZ][ZZ] = 5.0 * std::sqrt(2.0 / 3.0);
812 // TODO: Consider whether manually picking some positions would give better
814 data_.generateRandomRefPositions(1000);
815 data_.generateRandomTestPositions(100);
816 set_pbc(&data_.pbc_, PbcType::Xyz, data_.box_);
817 data_.computeReferences(&data_.pbc_);
821 NeighborhoodSearchTestData data_;
824 class RandomBox2DPBCData
827 static const NeighborhoodSearchTestData& get()
829 static RandomBox2DPBCData singleton;
830 return singleton.data_;
833 RandomBox2DPBCData() : data_(12345, 1.0)
835 data_.box_[XX][XX] = 10.0;
836 data_.box_[YY][YY] = 7.0;
837 data_.box_[ZZ][ZZ] = 5.0;
838 // TODO: Consider whether manually picking some positions would give better
840 data_.generateRandomRefPositions(1000);
841 data_.generateRandomTestPositions(100);
842 set_pbc(&data_.pbc_, PbcType::XY, data_.box_);
843 data_.computeReferences(&data_.pbc_);
847 NeighborhoodSearchTestData data_;
850 class RandomBoxNoPBCData
853 static const NeighborhoodSearchTestData& get()
855 static RandomBoxNoPBCData singleton;
856 return singleton.data_;
859 RandomBoxNoPBCData() : data_(12345, 1.0)
861 data_.box_[XX][XX] = 10.0;
862 data_.box_[YY][YY] = 5.0;
863 data_.box_[ZZ][ZZ] = 7.0;
864 // TODO: Consider whether manually picking some positions would give better
866 data_.generateRandomRefPositions(1000);
867 data_.generateRandomTestPositions(100);
868 set_pbc(&data_.pbc_, PbcType::No, data_.box_);
869 data_.computeReferences(nullptr);
873 NeighborhoodSearchTestData data_;
876 /********************************************************************
880 TEST_F(NeighborhoodSearchTest, SimpleSearch)
882 const NeighborhoodSearchTestData& data = RandomBoxFullPBCData::get();
884 nb_.setCutoff(data.cutoff_);
885 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Simple);
886 gmx::AnalysisNeighborhoodSearch search = nb_.initSearch(&data.pbc_, data.refPositions());
887 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
889 testIsWithin(&search, data);
890 testMinimumDistance(&search, data);
891 testNearestPoint(&search, data);
892 testPairSearch(&search, data);
895 testPairSearchIndexed(&nb_, data, 123);
898 TEST_F(NeighborhoodSearchTest, SimpleSearchXY)
900 const NeighborhoodSearchTestData& data = RandomBoxXYFullPBCData::get();
902 nb_.setCutoff(data.cutoff_);
903 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Simple);
905 gmx::AnalysisNeighborhoodSearch search = nb_.initSearch(&data.pbc_, data.refPositions());
906 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
908 testIsWithin(&search, data);
909 testMinimumDistance(&search, data);
910 testNearestPoint(&search, data);
911 testPairSearch(&search, data);
914 TEST_F(NeighborhoodSearchTest, GridSearchBox)
916 const NeighborhoodSearchTestData& data = RandomBoxFullPBCData::get();
918 nb_.setCutoff(data.cutoff_);
919 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
920 gmx::AnalysisNeighborhoodSearch search = nb_.initSearch(&data.pbc_, data.refPositions());
921 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
923 testIsWithin(&search, data);
924 testMinimumDistance(&search, data);
925 testNearestPoint(&search, data);
926 testPairSearch(&search, data);
929 testPairSearchIndexed(&nb_, data, 456);
932 TEST_F(NeighborhoodSearchTest, GridSearchTriclinic)
934 const NeighborhoodSearchTestData& data = RandomTriclinicFullPBCData::get();
936 nb_.setCutoff(data.cutoff_);
937 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
938 gmx::AnalysisNeighborhoodSearch search = nb_.initSearch(&data.pbc_, data.refPositions());
939 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
941 testPairSearch(&search, data);
944 TEST_F(NeighborhoodSearchTest, GridSearch2DPBC)
946 const NeighborhoodSearchTestData& data = RandomBox2DPBCData::get();
948 nb_.setCutoff(data.cutoff_);
949 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
950 gmx::AnalysisNeighborhoodSearch search = nb_.initSearch(&data.pbc_, data.refPositions());
951 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
953 testIsWithin(&search, data);
954 testMinimumDistance(&search, data);
955 testNearestPoint(&search, data);
956 testPairSearch(&search, data);
959 TEST_F(NeighborhoodSearchTest, GridSearchNoPBC)
961 const NeighborhoodSearchTestData& data = RandomBoxNoPBCData::get();
963 nb_.setCutoff(data.cutoff_);
964 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
965 gmx::AnalysisNeighborhoodSearch search = nb_.initSearch(&data.pbc_, data.refPositions());
966 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
968 testPairSearch(&search, data);
971 TEST_F(NeighborhoodSearchTest, GridSearchXYBox)
973 const NeighborhoodSearchTestData& data = RandomBoxXYFullPBCData::get();
975 nb_.setCutoff(data.cutoff_);
976 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
978 gmx::AnalysisNeighborhoodSearch search = nb_.initSearch(&data.pbc_, data.refPositions());
979 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
981 testIsWithin(&search, data);
982 testMinimumDistance(&search, data);
983 testNearestPoint(&search, data);
984 testPairSearch(&search, data);
987 TEST_F(NeighborhoodSearchTest, SimpleSelfPairsSearch)
989 const NeighborhoodSearchTestData& data = TrivialSelfPairsTestData::get();
991 nb_.setCutoff(data.cutoff_);
992 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Simple);
993 gmx::AnalysisNeighborhoodSearch search = nb_.initSearch(&data.pbc_, data.refPositions());
994 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
996 testPairSearchFull(&search, data, data.testPositions(), nullptr, {}, {}, true);
999 TEST_F(NeighborhoodSearchTest, GridSelfPairsSearch)
1001 const NeighborhoodSearchTestData& data = RandomBoxSelfPairsData::get();
1003 nb_.setCutoff(data.cutoff_);
1004 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
1005 gmx::AnalysisNeighborhoodSearch search = nb_.initSearch(&data.pbc_, data.refPositions());
1006 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
1008 testPairSearchFull(&search, data, data.testPositions(), nullptr, {}, {}, true);
1011 TEST_F(NeighborhoodSearchTest, HandlesConcurrentSearches)
1013 const NeighborhoodSearchTestData& data = TrivialTestData::get();
1015 nb_.setCutoff(data.cutoff_);
1016 gmx::AnalysisNeighborhoodSearch search1 = nb_.initSearch(&data.pbc_, data.refPositions());
1017 gmx::AnalysisNeighborhoodSearch search2 = nb_.initSearch(&data.pbc_, data.refPositions());
1019 // These checks are fragile, and unfortunately depend on the random
1020 // engine used to create the test positions. There is no particular reason
1021 // why exactly particles 0 & 2 should have neighbors, but in this case they do.
1022 gmx::AnalysisNeighborhoodPairSearch pairSearch1 = search1.startPairSearch(data.testPosition(0));
1023 gmx::AnalysisNeighborhoodPairSearch pairSearch2 = search1.startPairSearch(data.testPosition(2));
1025 testPairSearch(&search2, data);
1027 gmx::AnalysisNeighborhoodPair pair;
1028 ASSERT_TRUE(pairSearch1.findNextPair(&pair))
1029 << "Test data did not contain any pairs for position 0 (problem in the test).";
1030 EXPECT_EQ(0, pair.testIndex());
1032 NeighborhoodSearchTestData::RefPair searchPair(pair.refIndex(), std::sqrt(pair.distance2()));
1033 EXPECT_TRUE(data.containsPair(0, searchPair));
1036 ASSERT_TRUE(pairSearch2.findNextPair(&pair))
1037 << "Test data did not contain any pairs for position 2 (problem in the test).";
1038 EXPECT_EQ(2, pair.testIndex());
1040 NeighborhoodSearchTestData::RefPair searchPair(pair.refIndex(), std::sqrt(pair.distance2()));
1041 EXPECT_TRUE(data.containsPair(2, searchPair));
1045 TEST_F(NeighborhoodSearchTest, HandlesNoPBC)
1047 const NeighborhoodSearchTestData& data = TrivialNoPBCTestData::get();
1049 nb_.setCutoff(data.cutoff_);
1050 gmx::AnalysisNeighborhoodSearch search = nb_.initSearch(&data.pbc_, data.refPositions());
1051 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
1053 testIsWithin(&search, data);
1054 testMinimumDistance(&search, data);
1055 testNearestPoint(&search, data);
1056 testPairSearch(&search, data);
1059 TEST_F(NeighborhoodSearchTest, HandlesNullPBC)
1061 const NeighborhoodSearchTestData& data = TrivialNoPBCTestData::get();
1063 nb_.setCutoff(data.cutoff_);
1064 gmx::AnalysisNeighborhoodSearch search = nb_.initSearch(nullptr, data.refPositions());
1065 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
1067 testIsWithin(&search, data);
1068 testMinimumDistance(&search, data);
1069 testNearestPoint(&search, data);
1070 testPairSearch(&search, data);
1073 TEST_F(NeighborhoodSearchTest, HandlesSkippingPairs)
1075 const NeighborhoodSearchTestData& data = TrivialTestData::get();
1077 nb_.setCutoff(data.cutoff_);
1078 gmx::AnalysisNeighborhoodSearch search = nb_.initSearch(&data.pbc_, data.refPositions());
1079 gmx::AnalysisNeighborhoodPairSearch pairSearch = search.startPairSearch(data.testPositions());
1080 gmx::AnalysisNeighborhoodPair pair;
1081 // TODO: This test needs to be adjusted if the grid search gets optimized
1082 // to loop over the test positions in cell order (first, the ordering
1083 // assumption here breaks, and second, it then needs to be tested
1084 // separately for simple and grid searches).
1085 int currentIndex = 0;
1086 while (pairSearch.findNextPair(&pair))
1088 while (currentIndex < pair.testIndex())
1092 EXPECT_EQ(currentIndex, pair.testIndex());
1093 NeighborhoodSearchTestData::RefPair searchPair(pair.refIndex(), std::sqrt(pair.distance2()));
1094 EXPECT_TRUE(data.containsPair(currentIndex, searchPair));
1095 pairSearch.skipRemainingPairsForTestPosition();
1100 TEST_F(NeighborhoodSearchTest, SimpleSearchExclusions)
1102 const NeighborhoodSearchTestData& data = RandomBoxFullPBCData::get();
1104 ExclusionsHelper helper(data.refPosCount_, data.testPositions_.size());
1105 helper.generateExclusions();
1107 nb_.setCutoff(data.cutoff_);
1108 nb_.setTopologyExclusions(helper.exclusions());
1109 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Simple);
1110 gmx::AnalysisNeighborhoodSearch search =
1111 nb_.initSearch(&data.pbc_, data.refPositions().exclusionIds(helper.refPosIds()));
1112 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
1114 testPairSearchFull(&search,
1116 data.testPositions().exclusionIds(helper.testPosIds()),
1117 helper.exclusions(),
1123 TEST_F(NeighborhoodSearchTest, GridSearchExclusions)
1125 const NeighborhoodSearchTestData& data = RandomBoxFullPBCData::get();
1127 ExclusionsHelper helper(data.refPosCount_, data.testPositions_.size());
1128 helper.generateExclusions();
1130 nb_.setCutoff(data.cutoff_);
1131 nb_.setTopologyExclusions(helper.exclusions());
1132 nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
1133 gmx::AnalysisNeighborhoodSearch search =
1134 nb_.initSearch(&data.pbc_, data.refPositions().exclusionIds(helper.refPosIds()));
1135 ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
1137 testPairSearchFull(&search,
1139 data.testPositions().exclusionIds(helper.testPosIds()),
1140 helper.exclusions(),