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