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