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