src/gromacs/selection/tests/nbsearch.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, 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.
   8  *
   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.
  13  *
  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.
  18  *
  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.
  23  *
  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.
  31  *
  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.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Tests selection neighborhood searching.
  38  *
  39  * \todo
  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.
  43  *
  44  * \author Teemu Murtola <teemu.murtola@gmail.com>
  45  * \ingroup module_selection
  46  */
  47 #include "gmxpre.h"
  48
  49 #include "gromacs/selection/nbsearch.h"
  50
  51 #include <cmath>
  52
  53 #include <algorithm>
  54 #include <limits>
  55 #include <map>
  56 #include <numeric>
  57 #include <vector>
  58
  59 #include <gtest/gtest.h>
  60
  61 #include "gromacs/math/functions.h"
  62 #include "gromacs/math/vec.h"
  63 #include "gromacs/pbcutil/pbc.h"
  64 #include "gromacs/random/threefry.h"
  65 #include "gromacs/random/uniformrealdistribution.h"
  66 #include "gromacs/topology/block.h"
  67 #include "gromacs/utility/smalloc.h"
  68 #include "gromacs/utility/stringutil.h"
  69
  70 #include "testutils/testasserts.h"
  71
  72
  73 namespace
  74 {
  75
  76 /********************************************************************
  77  * NeighborhoodSearchTestData
  78  */
  79
  80 class NeighborhoodSearchTestData
  81 {
  82     public:
  83         struct RefPair
  84         {
  85             RefPair(int refIndex, real distance)
  86                 : refIndex(refIndex), distance(distance), bFound(false),
  87                   bExcluded(false), bIndexed(true)
  88             {
  89             }
  90
  91             bool operator<(const RefPair &other) const
  92             {
  93                 return refIndex < other.refIndex;
  94             }
  95
  96             int                 refIndex;
  97             real                distance;
  98             // The variables below are state variables that are only used
  99             // during the actual testing after creating a copy of the reference
 100             // pair list, not as part of the reference data.
 101             // Simpler to have just a single structure for both purposes.
 102             bool                bFound;
 103             bool                bExcluded;
 104             bool                bIndexed;
 105         };
 106
 107         struct TestPosition
 108         {
 109             explicit TestPosition(const rvec x)
 110                 : refMinDist(0.0), refNearestPoint(-1)
 111             {
 112                 copy_rvec(x, this->x);
 113             }
 114
 115             rvec                 x;
 116             real                 refMinDist;
 117             int                  refNearestPoint;
 118             std::vector<RefPair> refPairs;
 119         };
 120
 121         typedef std::vector<TestPosition> TestPositionList;
 122
 123         NeighborhoodSearchTestData(uint64_t seed, real cutoff);
 124
 125         gmx::AnalysisNeighborhoodPositions refPositions() const
 126         {
 127             return gmx::AnalysisNeighborhoodPositions(refPos_);
 128         }
 129         gmx::AnalysisNeighborhoodPositions testPositions() const
 130         {
 131             if (testPos_.empty())
 132             {
 133                 testPos_.reserve(testPositions_.size());
 134                 for (size_t i = 0; i < testPositions_.size(); ++i)
 135                 {
 136                     testPos_.emplace_back(testPositions_[i].x);
 137                 }
 138             }
 139             return gmx::AnalysisNeighborhoodPositions(testPos_);
 140         }
 141         gmx::AnalysisNeighborhoodPositions testPosition(int index) const
 142         {
 143             return testPositions().selectSingleFromArray(index);
 144         }
 145
 146         void addTestPosition(const rvec x)
 147         {
 148             GMX_RELEASE_ASSERT(testPos_.empty(),
 149                                "Cannot add positions after testPositions() call");
 150             testPositions_.emplace_back(x);
 151         }
 152         gmx::RVec generateRandomPosition();
 153         std::vector<int> generateIndex(int count, uint64_t seed) const;
 154         void generateRandomRefPositions(int count);
 155         void generateRandomTestPositions(int count);
 156         void useRefPositionsAsTestPositions();
 157         void computeReferences(t_pbc *pbc)
 158         {
 159             computeReferencesInternal(pbc, false);
 160         }
 161         void computeReferencesXY(t_pbc *pbc)
 162         {
 163             computeReferencesInternal(pbc, true);
 164         }
 165
 166         bool containsPair(int testIndex, const RefPair &pair) const
 167         {
 168             const std::vector<RefPair>          &refPairs = testPositions_[testIndex].refPairs;
 169             std::vector<RefPair>::const_iterator foundRefPair
 170                 = std::lower_bound(refPairs.begin(), refPairs.end(), pair);
 171             return !(foundRefPair == refPairs.end() || foundRefPair->refIndex != pair.refIndex);
 172         }
 173
 174         // Return a tolerance that accounts for the magnitudes of the coordinates
 175         // when doing subtractions, so that we set the ULP tolerance relative to the
 176         // coordinate values rather than their difference.
 177         // i.e., 10.0-9.9999999 will achieve a few ULP accuracy relative
 178         // to 10.0, but a much larger error relative to the difference.
 179         gmx::test::FloatingPointTolerance relativeTolerance() const
 180         {
 181             real magnitude = std::max(box_[XX][XX], std::max(box_[YY][YY], box_[ZZ][ZZ]));
 182             return gmx::test::relativeToleranceAsUlp(magnitude, 4);
 183         }
 184
 185         gmx::DefaultRandomEngine         rng_;
 186         real                             cutoff_;
 187         matrix                           box_;
 188         t_pbc                            pbc_;
 189         int                              refPosCount_;
 190         std::vector<gmx::RVec>           refPos_;
 191         TestPositionList                 testPositions_;
 192
 193     private:
 194         void computeReferencesInternal(t_pbc *pbc, bool bXY);
 195
 196         mutable std::vector<gmx::RVec>   testPos_;
 197 };
 198
 199 //! Shorthand for a collection of reference pairs.
 200 typedef std::vector<NeighborhoodSearchTestData::RefPair> RefPairList;
 201
 202 NeighborhoodSearchTestData::NeighborhoodSearchTestData(uint64_t seed, real cutoff)
 203     : rng_(seed), cutoff_(cutoff), refPosCount_(0)
 204 {
 205     clear_mat(box_);
 206     set_pbc(&pbc_, epbcNONE, box_);
 207 }
 208
 209 gmx::RVec NeighborhoodSearchTestData::generateRandomPosition()
 210 {
 211     gmx::UniformRealDistribution<real>  dist;
 212     rvec fx, x;
 213     fx[XX] = dist(rng_);
 214     fx[YY] = dist(rng_);
 215     fx[ZZ] = dist(rng_);
 216     mvmul(box_, fx, x);
 217     // Add a small displacement to allow positions outside the box
 218     x[XX] += 0.2 * dist(rng_) - 0.1;
 219     x[YY] += 0.2 * dist(rng_) - 0.1;
 220     x[ZZ] += 0.2 * dist(rng_) - 0.1;
 221     return x;
 222 }
 223
 224 std::vector<int> NeighborhoodSearchTestData::generateIndex(int count, uint64_t seed) const
 225 {
 226     gmx::DefaultRandomEngine             rngIndex(seed);
 227     gmx::UniformRealDistribution<real>   dist;
 228     std::vector<int>                     result;
 229
 230     for (int i = 0; i < count; ++i)
 231     {
 232         if (dist(rngIndex) > 0.5)
 233         {
 234             result.push_back(i);
 235         }
 236     }
 237     return result;
 238 }
 239
 240 void NeighborhoodSearchTestData::generateRandomRefPositions(int count)
 241 {
 242     refPosCount_ = count;
 243     refPos_.reserve(count);
 244     for (int i = 0; i < count; ++i)
 245     {
 246         refPos_.push_back(generateRandomPosition());
 247     }
 248 }
 249
 250 void NeighborhoodSearchTestData::generateRandomTestPositions(int count)
 251 {
 252     testPositions_.reserve(count);
 253     for (int i = 0; i < count; ++i)
 254     {
 255         addTestPosition(generateRandomPosition());
 256     }
 257 }
 258
 259 void NeighborhoodSearchTestData::useRefPositionsAsTestPositions()
 260 {
 261     testPositions_.reserve(refPosCount_);
 262     for (const auto &refPos : refPos_)
 263     {
 264         addTestPosition(refPos);
 265     }
 266 }
 267
 268 void NeighborhoodSearchTestData::computeReferencesInternal(t_pbc *pbc, bool bXY)
 269 {
 270     real cutoff = cutoff_;
 271     if (cutoff <= 0)
 272     {
 273         cutoff = std::numeric_limits<real>::max();
 274     }
 275     for (TestPosition &testPos : testPositions_)
 276     {
 277         testPos.refMinDist      = cutoff;
 278         testPos.refNearestPoint = -1;
 279         testPos.refPairs.clear();
 280         for (int j = 0; j < refPosCount_; ++j)
 281         {
 282             rvec dx;
 283             if (pbc != nullptr)
 284             {
 285                 pbc_dx(pbc, testPos.x, refPos_[j], dx);
 286             }
 287             else
 288             {
 289                 rvec_sub(testPos.x, refPos_[j], dx);
 290             }
 291             // TODO: This may not work intuitively for 2D with the third box
 292             // vector not parallel to the Z axis, but neither does the actual
 293             // neighborhood search.
 294             const real dist =
 295                 !bXY ? norm(dx) : std::hypot(dx[XX], dx[YY]);
 296             if (dist < testPos.refMinDist)
 297             {
 298                 testPos.refMinDist      = dist;
 299                 testPos.refNearestPoint = j;
 300             }
 301             if (dist > 0 && dist <= cutoff)
 302             {
 303                 RefPair pair(j, dist);
 304                 GMX_RELEASE_ASSERT(testPos.refPairs.empty() || testPos.refPairs.back() < pair,
 305                                    "Reference pairs should be generated in sorted order");
 306                 testPos.refPairs.push_back(pair);
 307             }
 308         }
 309     }
 310 }
 311
 312 /********************************************************************
 313  * ExclusionsHelper
 314  */
 315
 316 class ExclusionsHelper
 317 {
 318     public:
 319         static void markExcludedPairs(RefPairList *refPairs, int testIndex,
 320                                       const t_blocka *excls);
 321
 322         ExclusionsHelper(int refPosCount, int testPosCount);
 323
 324         void generateExclusions();
 325
 326         const t_blocka *exclusions() const { return &excls_; }
 327
 328         gmx::ArrayRef<const int> refPosIds() const
 329         {
 330             return gmx::makeArrayRef(exclusionIds_).subArray(0, refPosCount_);
 331         }
 332         gmx::ArrayRef<const int> testPosIds() const
 333         {
 334             return gmx::makeArrayRef(exclusionIds_).subArray(0, testPosCount_);
 335         }
 336
 337     private:
 338         int              refPosCount_;
 339         int              testPosCount_;
 340         std::vector<int> exclusionIds_;
 341         std::vector<int> exclsIndex_;
 342         std::vector<int> exclsAtoms_;
 343         t_blocka         excls_;
 344 };
 345
 346 // static
 347 void ExclusionsHelper::markExcludedPairs(RefPairList *refPairs, int testIndex,
 348                                          const t_blocka *excls)
 349 {
 350     int count = 0;
 351     for (int i = excls->index[testIndex]; i < excls->index[testIndex + 1]; ++i)
 352     {
 353         const int                           excludedIndex = excls->a[i];
 354         NeighborhoodSearchTestData::RefPair searchPair(excludedIndex, 0.0);
 355         RefPairList::iterator               excludedRefPair
 356             = std::lower_bound(refPairs->begin(), refPairs->end(), searchPair);
 357         if (excludedRefPair != refPairs->end()
 358             && excludedRefPair->refIndex == excludedIndex)
 359         {
 360             excludedRefPair->bFound    = true;
 361             excludedRefPair->bExcluded = true;
 362             ++count;
 363         }
 364     }
 365 }
 366
 367 ExclusionsHelper::ExclusionsHelper(int refPosCount, int testPosCount)
 368     : refPosCount_(refPosCount), testPosCount_(testPosCount)
 369 {
 370     // Generate an array of 0, 1, 2, ...
 371     // TODO: Make the tests work also with non-trivial exclusion IDs,
 372     // and test that.
 373     exclusionIds_.resize(std::max(refPosCount, testPosCount), 1);
 374     exclusionIds_[0] = 0;
 375     std::partial_sum(exclusionIds_.begin(), exclusionIds_.end(),
 376                      exclusionIds_.begin());
 377
 378     excls_.nr           = 0;
 379     excls_.index        = nullptr;
 380     excls_.nra          = 0;
 381     excls_.a            = nullptr;
 382     excls_.nalloc_index = 0;
 383     excls_.nalloc_a     = 0;
 384 }
 385
 386 void ExclusionsHelper::generateExclusions()
 387 {
 388     // TODO: Consider a better set of test data, where the density of the
 389     // particles would be higher, or where the exclusions would not be random,
 390     // to make a higher percentage of the exclusions to actually be within the
 391     // cutoff.
 392     exclsIndex_.reserve(testPosCount_ + 1);
 393     exclsAtoms_.reserve(testPosCount_ * 20);
 394     exclsIndex_.push_back(0);
 395     for (int i = 0; i < testPosCount_; ++i)
 396     {
 397         for (int j = 0; j < 20; ++j)
 398         {
 399             exclsAtoms_.push_back(i + j*3);
 400         }
 401         exclsIndex_.push_back(exclsAtoms_.size());
 402     }
 403     excls_.nr    = exclsIndex_.size();
 404     excls_.index = exclsIndex_.data();
 405     excls_.nra   = exclsAtoms_.size();
 406     excls_.a     = exclsAtoms_.data();
 407 }
 408
 409 /********************************************************************
 410  * NeighborhoodSearchTest
 411  */
 412
 413 class NeighborhoodSearchTest : public ::testing::Test
 414 {
 415     public:
 416         void testIsWithin(gmx::AnalysisNeighborhoodSearch  *search,
 417                           const NeighborhoodSearchTestData &data);
 418         void testMinimumDistance(gmx::AnalysisNeighborhoodSearch  *search,
 419                                  const NeighborhoodSearchTestData &data);
 420         void testNearestPoint(gmx::AnalysisNeighborhoodSearch  *search,
 421                               const NeighborhoodSearchTestData &data);
 422         void testPairSearch(gmx::AnalysisNeighborhoodSearch  *search,
 423                             const NeighborhoodSearchTestData &data);
 424         void testPairSearchIndexed(gmx::AnalysisNeighborhood        *nb,
 425                                    const NeighborhoodSearchTestData &data,
 426                                    uint64_t                          seed);
 427         void testPairSearchFull(gmx::AnalysisNeighborhoodSearch          *search,
 428                                 const NeighborhoodSearchTestData         &data,
 429                                 const gmx::AnalysisNeighborhoodPositions &pos,
 430                                 const t_blocka                           *excls,
 431                                 const gmx::ArrayRef<const int>           &refIndices,
 432                                 const gmx::ArrayRef<const int>           &testIndices,
 433                                 bool                                      selfPairs);
 434
 435         gmx::AnalysisNeighborhood        nb_;
 436 };
 437
 438 void NeighborhoodSearchTest::testIsWithin(
 439         gmx::AnalysisNeighborhoodSearch  *search,
 440         const NeighborhoodSearchTestData &data)
 441 {
 442     NeighborhoodSearchTestData::TestPositionList::const_iterator i;
 443     for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
 444     {
 445         const bool bWithin = (i->refMinDist <= data.cutoff_);
 446         EXPECT_EQ(bWithin, search->isWithin(i->x))
 447         << "Distance is " << i->refMinDist;
 448     }
 449 }
 450
 451 void NeighborhoodSearchTest::testMinimumDistance(
 452         gmx::AnalysisNeighborhoodSearch  *search,
 453         const NeighborhoodSearchTestData &data)
 454 {
 455     NeighborhoodSearchTestData::TestPositionList::const_iterator i;
 456
 457     for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
 458     {
 459         const real refDist = i->refMinDist;
 460         EXPECT_REAL_EQ_TOL(refDist, search->minimumDistance(i->x), data.relativeTolerance());
 461     }
 462 }
 463
 464 void NeighborhoodSearchTest::testNearestPoint(
 465         gmx::AnalysisNeighborhoodSearch  *search,
 466         const NeighborhoodSearchTestData &data)
 467 {
 468     NeighborhoodSearchTestData::TestPositionList::const_iterator i;
 469     for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
 470     {
 471         const gmx::AnalysisNeighborhoodPair pair = search->nearestPoint(i->x);
 472         if (pair.isValid())
 473         {
 474             EXPECT_EQ(i->refNearestPoint, pair.refIndex());
 475             EXPECT_EQ(0, pair.testIndex());
 476             EXPECT_REAL_EQ_TOL(i->refMinDist, std::sqrt(pair.distance2()), data.relativeTolerance());
 477         }
 478         else
 479         {
 480             EXPECT_EQ(i->refNearestPoint, -1);
 481         }
 482     }
 483 }
 484
 485 //! Helper function for formatting test failure messages.
 486 std::string formatVector(const rvec x)
 487 {
 488     return gmx::formatString("[%.3f, %.3f, %.3f]", x[XX], x[YY], x[ZZ]);
 489 }
 490
 491 /*! \brief
 492  * Helper function to check that all expected pairs were found.
 493  */
 494 void checkAllPairsFound(const RefPairList &refPairs,
 495                         const std::vector<gmx::RVec> &refPos,
 496                         int testPosIndex, const rvec testPos)
 497 {
 498     // This could be elegantly expressed with Google Mock matchers, but that
 499     // has a significant effect on the runtime of the tests...
 500     int                         count = 0;
 501     RefPairList::const_iterator first;
 502     for (RefPairList::const_iterator i = refPairs.begin(); i != refPairs.end(); ++i)
 503     {
 504         if (!i->bFound)
 505         {
 506             ++count;
 507             first = i;
 508         }
 509     }
 510     if (count > 0)
 511     {
 512         ADD_FAILURE()
 513         << "Some pairs (" << count << "/" << refPairs.size() << ") "
 514         << "within the cutoff were not found. First pair:\n"
 515         << " Ref: " << first->refIndex << " at "
 516         << formatVector(refPos[first->refIndex]) << "\n"
 517         << "Test: " << testPosIndex << " at " << formatVector(testPos) << "\n"
 518         << "Dist: " << first->distance;
 519     }
 520 }
 521
 522 void NeighborhoodSearchTest::testPairSearch(
 523         gmx::AnalysisNeighborhoodSearch  *search,
 524         const NeighborhoodSearchTestData &data)
 525 {
 526     testPairSearchFull(search, data, data.testPositions(), nullptr,
 527                        {}, {}, false);
 528 }
 529
 530 void NeighborhoodSearchTest::testPairSearchIndexed(
 531         gmx::AnalysisNeighborhood        *nb,
 532         const NeighborhoodSearchTestData &data,
 533         uint64_t                          seed)
 534 {
 535     std::vector<int>                refIndices(data.generateIndex(data.refPos_.size(), seed++));
 536     std::vector<int>                testIndices(data.generateIndex(data.testPositions_.size(), seed++));
 537     gmx::AnalysisNeighborhoodSearch search =
 538         nb->initSearch(&data.pbc_,
 539                        data.refPositions().indexed(refIndices));
 540     testPairSearchFull(&search, data, data.testPositions(), nullptr,
 541                        refIndices, testIndices, false);
 542 }
 543
 544 void NeighborhoodSearchTest::testPairSearchFull(
 545         gmx::AnalysisNeighborhoodSearch          *search,
 546         const NeighborhoodSearchTestData         &data,
 547         const gmx::AnalysisNeighborhoodPositions &pos,
 548         const t_blocka                           *excls,
 549         const gmx::ArrayRef<const int>           &refIndices,
 550         const gmx::ArrayRef<const int>           &testIndices,
 551         bool                                      selfPairs)
 552 {
 553     std::map<int, RefPairList> refPairs;
 554     // TODO: Some parts of this code do not work properly if pos does not
 555     // initially contain all the test positions.
 556     if (testIndices.empty())
 557     {
 558         for (size_t i = 0; i < data.testPositions_.size(); ++i)
 559         {
 560             refPairs[i] = data.testPositions_[i].refPairs;
 561         }
 562     }
 563     else
 564     {
 565         for (int index : testIndices)
 566         {
 567             refPairs[index] = data.testPositions_[index].refPairs;
 568         }
 569     }
 570     if (excls != nullptr)
 571     {
 572         GMX_RELEASE_ASSERT(!selfPairs, "Self-pairs testing not implemented with exclusions");
 573         for (auto &entry : refPairs)
 574         {
 575             const int testIndex = entry.first;
 576             ExclusionsHelper::markExcludedPairs(&entry.second, testIndex, excls);
 577         }
 578     }
 579     if (!refIndices.empty())
 580     {
 581         GMX_RELEASE_ASSERT(!selfPairs, "Self-pairs testing not implemented with indexing");
 582         for (auto &entry : refPairs)
 583         {
 584             for (auto &refPair : entry.second)
 585             {
 586                 refPair.bIndexed = false;
 587             }
 588             for (int index : refIndices)
 589             {
 590                 NeighborhoodSearchTestData::RefPair searchPair(index, 0.0);
 591                 auto refPair = std::lower_bound(entry.second.begin(), entry.second.end(), searchPair);
 592                 if (refPair != entry.second.end() && refPair->refIndex == index)
 593                 {
 594                     refPair->bIndexed = true;
 595                 }
 596             }
 597             for (auto &refPair : entry.second)
 598             {
 599                 if (!refPair.bIndexed)
 600                 {
 601                     refPair.bFound = true;
 602                 }
 603             }
 604         }
 605     }
 606
 607     gmx::AnalysisNeighborhoodPositions  posCopy(pos);
 608     if (!testIndices.empty())
 609     {
 610         posCopy.indexed(testIndices);
 611     }
 612     gmx::AnalysisNeighborhoodPairSearch pairSearch
 613         = selfPairs
 614             ? search->startSelfPairSearch()
 615             : search->startPairSearch(posCopy);
 616     gmx::AnalysisNeighborhoodPair       pair;
 617     while (pairSearch.findNextPair(&pair))
 618     {
 619         const int testIndex =
 620             (testIndices.empty() ? pair.testIndex() : testIndices[pair.testIndex()]);
 621         const int refIndex =
 622             (refIndices.empty() ? pair.refIndex() : refIndices[pair.refIndex()]);
 623
 624         if (refPairs.count(testIndex) == 0)
 625         {
 626             ADD_FAILURE()
 627             << "Expected: No pairs are returned for test position " << testIndex << ".\n"
 628             << "  Actual: Pair with ref " << refIndex << " is returned.";
 629             continue;
 630         }
 631         NeighborhoodSearchTestData::RefPair searchPair(refIndex,
 632                                                        std::sqrt(pair.distance2()));
 633         const auto foundRefPair
 634             = std::lower_bound(refPairs[testIndex].begin(), refPairs[testIndex].end(),
 635                                searchPair);
 636         if (foundRefPair == refPairs[testIndex].end() || foundRefPair->refIndex != refIndex)
 637         {
 638             ADD_FAILURE()
 639             << "Expected: Pair (ref: " << refIndex << ", test: " << testIndex
 640             << ") is not within the cutoff.\n"
 641             << "  Actual: It is returned.";
 642         }
 643         else if (foundRefPair->bExcluded)
 644         {
 645             ADD_FAILURE()
 646             << "Expected: Pair (ref: " << refIndex << ", test: " << testIndex
 647             << ") is excluded from the search.\n"
 648             << "  Actual: It is returned.";
 649         }
 650         else if (!foundRefPair->bIndexed)
 651         {
 652             ADD_FAILURE()
 653             << "Expected: Pair (ref: " << refIndex << ", test: " << testIndex
 654             << ") is not part of the indexed set.\n"
 655             << "  Actual: It is returned.";
 656         }
 657         else if (foundRefPair->bFound)
 658         {
 659             ADD_FAILURE()
 660             << "Expected: Pair (ref: " << refIndex << ", test: " << testIndex
 661             << ") is returned only once.\n"
 662             << "  Actual: It is returned multiple times.";
 663             return;
 664         }
 665         else
 666         {
 667             foundRefPair->bFound = true;
 668
 669             EXPECT_REAL_EQ_TOL(foundRefPair->distance, searchPair.distance, data.relativeTolerance())
 670             << "Distance computed by the neighborhood search does not match.";
 671             if (selfPairs)
 672             {
 673                 searchPair = NeighborhoodSearchTestData::RefPair(testIndex, 0.0);
 674                 const auto otherRefPair
 675                     = std::lower_bound(refPairs[refIndex].begin(), refPairs[refIndex].end(),
 676                                        searchPair);
 677                 GMX_RELEASE_ASSERT(otherRefPair != refPairs[refIndex].end(),
 678                                    "Precomputed reference data is not symmetric");
 679                 otherRefPair->bFound = true;
 680             }
 681         }
 682     }
 683
 684     for (auto &entry : refPairs)
 685     {
 686         const int testIndex = entry.first;
 687         checkAllPairsFound(entry.second, data.refPos_, testIndex,
 688                            data.testPositions_[testIndex].x);
 689     }
 690 }
 691
 692 /********************************************************************
 693  * Test data generation
 694  */
 695
 696 class TrivialTestData
 697 {
 698     public:
 699         static const NeighborhoodSearchTestData &get()
 700         {
 701             static TrivialTestData singleton;
 702             return singleton.data_;
 703         }
 704
 705         TrivialTestData() : data_(12345, 1.0)
 706         {
 707             // Make the box so small we are virtually guaranteed to have
 708             // several neighbors for the five test positions
 709             data_.box_[XX][XX] = 3.0;
 710             data_.box_[YY][YY] = 3.0;
 711             data_.box_[ZZ][ZZ] = 3.0;
 712             data_.generateRandomRefPositions(10);
 713             data_.generateRandomTestPositions(5);
 714             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 715             data_.computeReferences(&data_.pbc_);
 716         }
 717
 718     private:
 719         NeighborhoodSearchTestData data_;
 720 };
 721
 722 class TrivialSelfPairsTestData
 723 {
 724     public:
 725         static const NeighborhoodSearchTestData &get()
 726         {
 727             static TrivialSelfPairsTestData singleton;
 728             return singleton.data_;
 729         }
 730
 731         TrivialSelfPairsTestData() : data_(12345, 1.0)
 732         {
 733             data_.box_[XX][XX] = 3.0;
 734             data_.box_[YY][YY] = 3.0;
 735             data_.box_[ZZ][ZZ] = 3.0;
 736             data_.generateRandomRefPositions(20);
 737             data_.useRefPositionsAsTestPositions();
 738             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 739             data_.computeReferences(&data_.pbc_);
 740         }
 741
 742     private:
 743         NeighborhoodSearchTestData data_;
 744 };
 745
 746 class TrivialNoPBCTestData
 747 {
 748     public:
 749         static const NeighborhoodSearchTestData &get()
 750         {
 751             static TrivialNoPBCTestData singleton;
 752             return singleton.data_;
 753         }
 754
 755         TrivialNoPBCTestData() : data_(12345, 1.0)
 756         {
 757             data_.generateRandomRefPositions(10);
 758             data_.generateRandomTestPositions(5);
 759             data_.computeReferences(nullptr);
 760         }
 761
 762     private:
 763         NeighborhoodSearchTestData data_;
 764 };
 765
 766 class RandomBoxFullPBCData
 767 {
 768     public:
 769         static const NeighborhoodSearchTestData &get()
 770         {
 771             static RandomBoxFullPBCData singleton;
 772             return singleton.data_;
 773         }
 774
 775         RandomBoxFullPBCData() : data_(12345, 1.0)
 776         {
 777             data_.box_[XX][XX] = 10.0;
 778             data_.box_[YY][YY] = 5.0;
 779             data_.box_[ZZ][ZZ] = 7.0;
 780             // TODO: Consider whether manually picking some positions would give better
 781             // test coverage.
 782             data_.generateRandomRefPositions(1000);
 783             data_.generateRandomTestPositions(100);
 784             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 785             data_.computeReferences(&data_.pbc_);
 786         }
 787
 788     private:
 789         NeighborhoodSearchTestData data_;
 790 };
 791
 792 class RandomBoxSelfPairsData
 793 {
 794     public:
 795         static const NeighborhoodSearchTestData &get()
 796         {
 797             static RandomBoxSelfPairsData singleton;
 798             return singleton.data_;
 799         }
 800
 801         RandomBoxSelfPairsData() : data_(12345, 1.0)
 802         {
 803             data_.box_[XX][XX] = 10.0;
 804             data_.box_[YY][YY] = 5.0;
 805             data_.box_[ZZ][ZZ] = 7.0;
 806             data_.generateRandomRefPositions(1000);
 807             data_.useRefPositionsAsTestPositions();
 808             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 809             data_.computeReferences(&data_.pbc_);
 810         }
 811
 812     private:
 813         NeighborhoodSearchTestData data_;
 814 };
 815
 816 class RandomBoxXYFullPBCData
 817 {
 818     public:
 819         static const NeighborhoodSearchTestData &get()
 820         {
 821             static RandomBoxXYFullPBCData singleton;
 822             return singleton.data_;
 823         }
 824
 825         RandomBoxXYFullPBCData() : data_(54321, 1.0)
 826         {
 827             data_.box_[XX][XX] = 10.0;
 828             data_.box_[YY][YY] = 5.0;
 829             data_.box_[ZZ][ZZ] = 7.0;
 830             // TODO: Consider whether manually picking some positions would give better
 831             // test coverage.
 832             data_.generateRandomRefPositions(1000);
 833             data_.generateRandomTestPositions(100);
 834             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 835             data_.computeReferencesXY(&data_.pbc_);
 836         }
 837
 838     private:
 839         NeighborhoodSearchTestData data_;
 840 };
 841
 842 class RandomTriclinicFullPBCData
 843 {
 844     public:
 845         static const NeighborhoodSearchTestData &get()
 846         {
 847             static RandomTriclinicFullPBCData singleton;
 848             return singleton.data_;
 849         }
 850
 851         RandomTriclinicFullPBCData() : data_(12345, 1.0)
 852         {
 853             data_.box_[XX][XX] = 5.0;
 854             data_.box_[YY][XX] = 2.5;
 855             data_.box_[YY][YY] = 2.5*std::sqrt(3.0);
 856             data_.box_[ZZ][XX] = 2.5;
 857             data_.box_[ZZ][YY] = 2.5*std::sqrt(1.0/3.0);
 858             data_.box_[ZZ][ZZ] = 5.0*std::sqrt(2.0/3.0);
 859             // TODO: Consider whether manually picking some positions would give better
 860             // test coverage.
 861             data_.generateRandomRefPositions(1000);
 862             data_.generateRandomTestPositions(100);
 863             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 864             data_.computeReferences(&data_.pbc_);
 865         }
 866
 867     private:
 868         NeighborhoodSearchTestData data_;
 869 };
 870
 871 class RandomBox2DPBCData
 872 {
 873     public:
 874         static const NeighborhoodSearchTestData &get()
 875         {
 876             static RandomBox2DPBCData singleton;
 877             return singleton.data_;
 878         }
 879
 880         RandomBox2DPBCData() : data_(12345, 1.0)
 881         {
 882             data_.box_[XX][XX] = 10.0;
 883             data_.box_[YY][YY] = 7.0;
 884             data_.box_[ZZ][ZZ] = 5.0;
 885             // TODO: Consider whether manually picking some positions would give better
 886             // test coverage.
 887             data_.generateRandomRefPositions(1000);
 888             data_.generateRandomTestPositions(100);
 889             set_pbc(&data_.pbc_, epbcXY, data_.box_);
 890             data_.computeReferences(&data_.pbc_);
 891         }
 892
 893     private:
 894         NeighborhoodSearchTestData data_;
 895 };
 896
 897 class RandomBoxNoPBCData
 898 {
 899     public:
 900         static const NeighborhoodSearchTestData &get()
 901         {
 902             static RandomBoxNoPBCData singleton;
 903             return singleton.data_;
 904         }
 905
 906         RandomBoxNoPBCData() : data_(12345, 1.0)
 907         {
 908             data_.box_[XX][XX] = 10.0;
 909             data_.box_[YY][YY] = 5.0;
 910             data_.box_[ZZ][ZZ] = 7.0;
 911             // TODO: Consider whether manually picking some positions would give better
 912             // test coverage.
 913             data_.generateRandomRefPositions(1000);
 914             data_.generateRandomTestPositions(100);
 915             set_pbc(&data_.pbc_, epbcNONE, data_.box_);
 916             data_.computeReferences(nullptr);
 917         }
 918
 919     private:
 920         NeighborhoodSearchTestData data_;
 921 };
 922
 923 /********************************************************************
 924  * Actual tests
 925  */
 926
 927 TEST_F(NeighborhoodSearchTest, SimpleSearch)
 928 {
 929     const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
 930
 931     nb_.setCutoff(data.cutoff_);
 932     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Simple);
 933     gmx::AnalysisNeighborhoodSearch search =
 934         nb_.initSearch(&data.pbc_, data.refPositions());
 935     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
 936
 937     testIsWithin(&search, data);
 938     testMinimumDistance(&search, data);
 939     testNearestPoint(&search, data);
 940     testPairSearch(&search, data);
 941
 942     search.reset();
 943     testPairSearchIndexed(&nb_, data, 123);
 944 }
 945
 946 TEST_F(NeighborhoodSearchTest, SimpleSearchXY)
 947 {
 948     const NeighborhoodSearchTestData &data = RandomBoxXYFullPBCData::get();
 949
 950     nb_.setCutoff(data.cutoff_);
 951     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Simple);
 952     nb_.setXYMode(true);
 953     gmx::AnalysisNeighborhoodSearch search =
 954         nb_.initSearch(&data.pbc_, data.refPositions());
 955     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
 956
 957     testIsWithin(&search, data);
 958     testMinimumDistance(&search, data);
 959     testNearestPoint(&search, data);
 960     testPairSearch(&search, data);
 961 }
 962
 963 TEST_F(NeighborhoodSearchTest, GridSearchBox)
 964 {
 965     const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
 966
 967     nb_.setCutoff(data.cutoff_);
 968     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
 969     gmx::AnalysisNeighborhoodSearch search =
 970         nb_.initSearch(&data.pbc_, data.refPositions());
 971     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
 972
 973     testIsWithin(&search, data);
 974     testMinimumDistance(&search, data);
 975     testNearestPoint(&search, data);
 976     testPairSearch(&search, data);
 977
 978     search.reset();
 979     testPairSearchIndexed(&nb_, data, 456);
 980 }
 981
 982 TEST_F(NeighborhoodSearchTest, GridSearchTriclinic)
 983 {
 984     const NeighborhoodSearchTestData &data = RandomTriclinicFullPBCData::get();
 985
 986     nb_.setCutoff(data.cutoff_);
 987     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
 988     gmx::AnalysisNeighborhoodSearch search =
 989         nb_.initSearch(&data.pbc_, data.refPositions());
 990     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
 991
 992     testPairSearch(&search, data);
 993 }
 994
 995 TEST_F(NeighborhoodSearchTest, GridSearch2DPBC)
 996 {
 997     const NeighborhoodSearchTestData &data = RandomBox2DPBCData::get();
 998
 999     nb_.setCutoff(data.cutoff_);
1000     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
1001     gmx::AnalysisNeighborhoodSearch search =
1002         nb_.initSearch(&data.pbc_, data.refPositions());
1003     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
1004
1005     testIsWithin(&search, data);
1006     testMinimumDistance(&search, data);
1007     testNearestPoint(&search, data);
1008     testPairSearch(&search, data);
1009 }
1010
1011 TEST_F(NeighborhoodSearchTest, GridSearchNoPBC)
1012 {
1013     const NeighborhoodSearchTestData &data = RandomBoxNoPBCData::get();
1014
1015     nb_.setCutoff(data.cutoff_);
1016     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
1017     gmx::AnalysisNeighborhoodSearch search =
1018         nb_.initSearch(&data.pbc_, data.refPositions());
1019     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
1020
1021     testPairSearch(&search, data);
1022 }
1023
1024 TEST_F(NeighborhoodSearchTest, GridSearchXYBox)
1025 {
1026     const NeighborhoodSearchTestData &data = RandomBoxXYFullPBCData::get();
1027
1028     nb_.setCutoff(data.cutoff_);
1029     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
1030     nb_.setXYMode(true);
1031     gmx::AnalysisNeighborhoodSearch search =
1032         nb_.initSearch(&data.pbc_, data.refPositions());
1033     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
1034
1035     testIsWithin(&search, data);
1036     testMinimumDistance(&search, data);
1037     testNearestPoint(&search, data);
1038     testPairSearch(&search, data);
1039 }
1040
1041 TEST_F(NeighborhoodSearchTest, SimpleSelfPairsSearch)
1042 {
1043     const NeighborhoodSearchTestData &data = TrivialSelfPairsTestData::get();
1044
1045     nb_.setCutoff(data.cutoff_);
1046     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Simple);
1047     gmx::AnalysisNeighborhoodSearch search =
1048         nb_.initSearch(&data.pbc_, data.refPositions());
1049     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
1050
1051     testPairSearchFull(&search, data, data.testPositions(), nullptr,
1052                        {}, {}, true);
1053 }
1054
1055 TEST_F(NeighborhoodSearchTest, GridSelfPairsSearch)
1056 {
1057     const NeighborhoodSearchTestData &data = RandomBoxSelfPairsData::get();
1058
1059     nb_.setCutoff(data.cutoff_);
1060     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
1061     gmx::AnalysisNeighborhoodSearch search =
1062         nb_.initSearch(&data.pbc_, data.refPositions());
1063     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
1064
1065     testPairSearchFull(&search, data, data.testPositions(), nullptr,
1066                        {}, {}, true);
1067 }
1068
1069 TEST_F(NeighborhoodSearchTest, HandlesConcurrentSearches)
1070 {
1071     const NeighborhoodSearchTestData &data = TrivialTestData::get();
1072
1073     nb_.setCutoff(data.cutoff_);
1074     gmx::AnalysisNeighborhoodSearch search1 =
1075         nb_.initSearch(&data.pbc_, data.refPositions());
1076     gmx::AnalysisNeighborhoodSearch search2 =
1077         nb_.initSearch(&data.pbc_, data.refPositions());
1078
1079     // These checks are fragile, and unfortunately depend on the random
1080     // engine used to create the test positions. There is no particular reason
1081     // why exactly particles 0 & 2 should have neighbors, but in this case they do.
1082     gmx::AnalysisNeighborhoodPairSearch pairSearch1 =
1083         search1.startPairSearch(data.testPosition(0));
1084     gmx::AnalysisNeighborhoodPairSearch pairSearch2 =
1085         search1.startPairSearch(data.testPosition(2));
1086
1087     testPairSearch(&search2, data);
1088
1089     gmx::AnalysisNeighborhoodPair pair;
1090     ASSERT_TRUE(pairSearch1.findNextPair(&pair))
1091     << "Test data did not contain any pairs for position 0 (problem in the test).";
1092     EXPECT_EQ(0, pair.testIndex());
1093     {
1094         NeighborhoodSearchTestData::RefPair searchPair(pair.refIndex(), std::sqrt(pair.distance2()));
1095         EXPECT_TRUE(data.containsPair(0, searchPair));
1096     }
1097
1098     ASSERT_TRUE(pairSearch2.findNextPair(&pair))
1099     << "Test data did not contain any pairs for position 2 (problem in the test).";
1100     EXPECT_EQ(2, pair.testIndex());
1101     {
1102         NeighborhoodSearchTestData::RefPair searchPair(pair.refIndex(), std::sqrt(pair.distance2()));
1103         EXPECT_TRUE(data.containsPair(2, searchPair));
1104     }
1105 }
1106
1107 TEST_F(NeighborhoodSearchTest, HandlesNoPBC)
1108 {
1109     const NeighborhoodSearchTestData &data = TrivialNoPBCTestData::get();
1110
1111     nb_.setCutoff(data.cutoff_);
1112     gmx::AnalysisNeighborhoodSearch search =
1113         nb_.initSearch(&data.pbc_, data.refPositions());
1114     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
1115
1116     testIsWithin(&search, data);
1117     testMinimumDistance(&search, data);
1118     testNearestPoint(&search, data);
1119     testPairSearch(&search, data);
1120 }
1121
1122 TEST_F(NeighborhoodSearchTest, HandlesNullPBC)
1123 {
1124     const NeighborhoodSearchTestData &data = TrivialNoPBCTestData::get();
1125
1126     nb_.setCutoff(data.cutoff_);
1127     gmx::AnalysisNeighborhoodSearch search =
1128         nb_.initSearch(nullptr, data.refPositions());
1129     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
1130
1131     testIsWithin(&search, data);
1132     testMinimumDistance(&search, data);
1133     testNearestPoint(&search, data);
1134     testPairSearch(&search, data);
1135 }
1136
1137 TEST_F(NeighborhoodSearchTest, HandlesSkippingPairs)
1138 {
1139     const NeighborhoodSearchTestData &data = TrivialTestData::get();
1140
1141     nb_.setCutoff(data.cutoff_);
1142     gmx::AnalysisNeighborhoodSearch     search =
1143         nb_.initSearch(&data.pbc_, data.refPositions());
1144     gmx::AnalysisNeighborhoodPairSearch pairSearch =
1145         search.startPairSearch(data.testPositions());
1146     gmx::AnalysisNeighborhoodPair       pair;
1147     // TODO: This test needs to be adjusted if the grid search gets optimized
1148     // to loop over the test positions in cell order (first, the ordering
1149     // assumption here breaks, and second, it then needs to be tested
1150     // separately for simple and grid searches).
1151     int currentIndex = 0;
1152     while (pairSearch.findNextPair(&pair))
1153     {
1154         while (currentIndex < pair.testIndex())
1155         {
1156             ++currentIndex;
1157         }
1158         EXPECT_EQ(currentIndex, pair.testIndex());
1159         NeighborhoodSearchTestData::RefPair searchPair(pair.refIndex(), std::sqrt(pair.distance2()));
1160         EXPECT_TRUE(data.containsPair(currentIndex, searchPair));
1161         pairSearch.skipRemainingPairsForTestPosition();
1162         ++currentIndex;
1163     }
1164 }
1165
1166 TEST_F(NeighborhoodSearchTest, SimpleSearchExclusions)
1167 {
1168     const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
1169
1170     ExclusionsHelper                  helper(data.refPosCount_, data.testPositions_.size());
1171     helper.generateExclusions();
1172
1173     nb_.setCutoff(data.cutoff_);
1174     nb_.setTopologyExclusions(helper.exclusions());
1175     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Simple);
1176     gmx::AnalysisNeighborhoodSearch search =
1177         nb_.initSearch(&data.pbc_,
1178                        data.refPositions().exclusionIds(helper.refPosIds()));
1179     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
1180
1181     testPairSearchFull(&search, data,
1182                        data.testPositions().exclusionIds(helper.testPosIds()),
1183                        helper.exclusions(), {},
1184                        {}, false);
1185 }
1186
1187 TEST_F(NeighborhoodSearchTest, GridSearchExclusions)
1188 {
1189     const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
1190
1191     ExclusionsHelper                  helper(data.refPosCount_, data.testPositions_.size());
1192     helper.generateExclusions();
1193
1194     nb_.setCutoff(data.cutoff_);
1195     nb_.setTopologyExclusions(helper.exclusions());
1196     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
1197     gmx::AnalysisNeighborhoodSearch search =
1198         nb_.initSearch(&data.pbc_,
1199                        data.refPositions().exclusionIds(helper.refPosIds()));
1200     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
1201
1202     testPairSearchFull(&search, data,
1203                        data.testPositions().exclusionIds(helper.testPosIds()),
1204                        helper.exclusions(), {},
1205                        {}, false);
1206 }
1207
1208 } // namespace