src/gromacs/selection/tests/nbsearch.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013,2014, 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 "gromacs/selection/nbsearch.h"
  48
  49 #include <gtest/gtest.h>
  50
  51 #include <cmath>
  52
  53 #include <algorithm>
  54 #include <limits>
  55 #include <numeric>
  56 #include <vector>
  57
  58 #include "gromacs/math/vec.h"
  59 #include "gromacs/pbcutil/pbc.h"
  60 #include "gromacs/random/random.h"
  61 #include "gromacs/topology/block.h"
  62 #include "gromacs/utility/smalloc.h"
  63
  64 #include "testutils/testasserts.h"
  65
  66 namespace
  67 {
  68
  69 /********************************************************************
  70  * NeighborhoodSearchTestData
  71  */
  72
  73 class NeighborhoodSearchTestData
  74 {
  75     public:
  76         struct RefPair
  77         {
  78             RefPair(int refIndex, real distance)
  79                 : refIndex(refIndex), distance(distance), bFound(false),
  80                   bExcluded(false)
  81             {
  82             }
  83
  84             bool operator<(const RefPair &other) const
  85             {
  86                 return refIndex < other.refIndex;
  87             }
  88
  89             int                 refIndex;
  90             real                distance;
  91             // The variables below are state variables that are only used
  92             // during the actual testing after creating a copy of the reference
  93             // pair list, not as part of the reference data.
  94             // Simpler to have just a single structure for both purposes.
  95             bool                bFound;
  96             bool                bExcluded;
  97         };
  98
  99         struct TestPosition
 100         {
 101             TestPosition() : refMinDist(0.0), refNearestPoint(-1)
 102             {
 103                 clear_rvec(x);
 104             }
 105             explicit TestPosition(const rvec x)
 106                 : refMinDist(0.0), refNearestPoint(-1)
 107             {
 108                 copy_rvec(x, this->x);
 109             }
 110
 111             rvec                 x;
 112             real                 refMinDist;
 113             int                  refNearestPoint;
 114             std::vector<RefPair> refPairs;
 115         };
 116
 117         typedef std::vector<TestPosition> TestPositionList;
 118
 119         NeighborhoodSearchTestData(int seed, real cutoff);
 120         ~NeighborhoodSearchTestData();
 121
 122         gmx::AnalysisNeighborhoodPositions refPositions() const
 123         {
 124             return gmx::AnalysisNeighborhoodPositions(refPos_, refPosCount_);
 125         }
 126         gmx::AnalysisNeighborhoodPositions testPositions() const
 127         {
 128             if (testPos_ == NULL)
 129             {
 130                 snew(testPos_, testPositions_.size());
 131                 for (size_t i = 0; i < testPositions_.size(); ++i)
 132                 {
 133                     copy_rvec(testPositions_[i].x, testPos_[i]);
 134                 }
 135             }
 136             return gmx::AnalysisNeighborhoodPositions(testPos_,
 137                                                       testPositions_.size());
 138         }
 139         gmx::AnalysisNeighborhoodPositions testPosition(int index) const
 140         {
 141             return testPositions().selectSingleFromArray(index);
 142         }
 143
 144         void addTestPosition(const rvec x)
 145         {
 146             GMX_RELEASE_ASSERT(testPos_ == NULL,
 147                                "Cannot add positions after testPositions() call");
 148             testPositions_.push_back(TestPosition(x));
 149         }
 150         void generateRandomPosition(rvec x);
 151         void generateRandomRefPositions(int count);
 152         void generateRandomTestPositions(int count);
 153         void computeReferences(t_pbc *pbc)
 154         {
 155             computeReferencesInternal(pbc, false);
 156         }
 157         void computeReferencesXY(t_pbc *pbc)
 158         {
 159             computeReferencesInternal(pbc, true);
 160         }
 161
 162         bool containsPair(int testIndex, const RefPair &pair) const
 163         {
 164             const std::vector<RefPair>          &refPairs = testPositions_[testIndex].refPairs;
 165             std::vector<RefPair>::const_iterator foundRefPair
 166                 = std::lower_bound(refPairs.begin(), refPairs.end(), pair);
 167             if (foundRefPair == refPairs.end() || foundRefPair->refIndex != pair.refIndex)
 168             {
 169                 return false;
 170             }
 171             return true;
 172         }
 173
 174         gmx_rng_t                        rng_;
 175         real                             cutoff_;
 176         matrix                           box_;
 177         t_pbc                            pbc_;
 178         int                              refPosCount_;
 179         rvec                            *refPos_;
 180         TestPositionList                 testPositions_;
 181
 182     private:
 183         void computeReferencesInternal(t_pbc *pbc, bool bXY);
 184
 185         mutable rvec                    *testPos_;
 186 };
 187
 188 //! Shorthand for a collection of reference pairs.
 189 typedef std::vector<NeighborhoodSearchTestData::RefPair> RefPairList;
 190
 191 NeighborhoodSearchTestData::NeighborhoodSearchTestData(int seed, real cutoff)
 192     : rng_(NULL), cutoff_(cutoff), refPosCount_(0), refPos_(NULL), testPos_(NULL)
 193 {
 194     // TODO: Handle errors.
 195     rng_ = gmx_rng_init(seed);
 196     clear_mat(box_);
 197     set_pbc(&pbc_, epbcNONE, box_);
 198 }
 199
 200 NeighborhoodSearchTestData::~NeighborhoodSearchTestData()
 201 {
 202     if (rng_ != NULL)
 203     {
 204         gmx_rng_destroy(rng_);
 205     }
 206     sfree(refPos_);
 207     sfree(testPos_);
 208 }
 209
 210 void NeighborhoodSearchTestData::generateRandomPosition(rvec x)
 211 {
 212     rvec fx;
 213     fx[XX] = gmx_rng_uniform_real(rng_);
 214     fx[YY] = gmx_rng_uniform_real(rng_);
 215     fx[ZZ] = gmx_rng_uniform_real(rng_);
 216     mvmul(box_, fx, x);
 217     // Add a small displacement to allow positions outside the box
 218     x[XX] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
 219     x[YY] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
 220     x[ZZ] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
 221 }
 222
 223 void NeighborhoodSearchTestData::generateRandomRefPositions(int count)
 224 {
 225     refPosCount_ = count;
 226     snew(refPos_, refPosCount_);
 227     for (int i = 0; i < refPosCount_; ++i)
 228     {
 229         generateRandomPosition(refPos_[i]);
 230     }
 231 }
 232
 233 void NeighborhoodSearchTestData::generateRandomTestPositions(int count)
 234 {
 235     testPositions_.reserve(count);
 236     for (int i = 0; i < count; ++i)
 237     {
 238         rvec x;
 239         generateRandomPosition(x);
 240         addTestPosition(x);
 241     }
 242 }
 243
 244 void NeighborhoodSearchTestData::computeReferencesInternal(t_pbc *pbc, bool bXY)
 245 {
 246     real cutoff = cutoff_;
 247     if (cutoff <= 0)
 248     {
 249         cutoff = std::numeric_limits<real>::max();
 250     }
 251     TestPositionList::iterator i;
 252     for (i = testPositions_.begin(); i != testPositions_.end(); ++i)
 253     {
 254         i->refMinDist      = cutoff;
 255         i->refNearestPoint = -1;
 256         i->refPairs.clear();
 257         for (int j = 0; j < refPosCount_; ++j)
 258         {
 259             rvec dx;
 260             if (pbc != NULL)
 261             {
 262                 pbc_dx(pbc, i->x, refPos_[j], dx);
 263             }
 264             else
 265             {
 266                 rvec_sub(i->x, refPos_[j], dx);
 267             }
 268             // TODO: This may not work intuitively for 2D with the third box
 269             // vector not parallel to the Z axis, but neither does the actual
 270             // neighborhood search.
 271             const real dist =
 272                 !bXY ? norm(dx) : sqrt(sqr(dx[XX]) + sqr(dx[YY]));
 273             if (dist < i->refMinDist)
 274             {
 275                 i->refMinDist      = dist;
 276                 i->refNearestPoint = j;
 277             }
 278             if (dist <= cutoff)
 279             {
 280                 RefPair pair(j, dist);
 281                 GMX_RELEASE_ASSERT(i->refPairs.empty() || i->refPairs.back() < pair,
 282                                    "Reference pairs should be generated in sorted order");
 283                 i->refPairs.push_back(pair);
 284             }
 285         }
 286     }
 287 }
 288
 289 /********************************************************************
 290  * ExclusionsHelper
 291  */
 292
 293 class ExclusionsHelper
 294 {
 295     public:
 296         static void markExcludedPairs(RefPairList *refPairs, int testIndex,
 297                                       const t_blocka *excls);
 298
 299         ExclusionsHelper(int refPosCount, int testPosCount);
 300
 301         void generateExclusions();
 302
 303         const t_blocka *exclusions() const { return &excls_; }
 304
 305         gmx::ConstArrayRef<int> refPosIds() const
 306         {
 307             return gmx::constArrayRefFromVector<int>(exclusionIds_.begin(),
 308                                                      exclusionIds_.begin() + refPosCount_);
 309         }
 310         gmx::ConstArrayRef<int> testPosIds() const
 311         {
 312             return gmx::constArrayRefFromVector<int>(exclusionIds_.begin(),
 313                                                      exclusionIds_.begin() + testPosCount_);
 314         }
 315
 316     private:
 317         int              refPosCount_;
 318         int              testPosCount_;
 319         std::vector<int> exclusionIds_;
 320         std::vector<int> exclsIndex_;
 321         std::vector<int> exclsAtoms_;
 322         t_blocka         excls_;
 323 };
 324
 325 // static
 326 void ExclusionsHelper::markExcludedPairs(RefPairList *refPairs, int testIndex,
 327                                          const t_blocka *excls)
 328 {
 329     int count = 0;
 330     for (int i = excls->index[testIndex]; i < excls->index[testIndex + 1]; ++i)
 331     {
 332         const int                           excludedIndex = excls->a[i];
 333         NeighborhoodSearchTestData::RefPair searchPair(excludedIndex, 0.0);
 334         RefPairList::iterator               excludedRefPair
 335             = std::lower_bound(refPairs->begin(), refPairs->end(), searchPair);
 336         if (excludedRefPair != refPairs->end()
 337             && excludedRefPair->refIndex == excludedIndex)
 338         {
 339             excludedRefPair->bFound    = true;
 340             excludedRefPair->bExcluded = true;
 341             ++count;
 342         }
 343     }
 344 }
 345
 346 ExclusionsHelper::ExclusionsHelper(int refPosCount, int testPosCount)
 347     : refPosCount_(refPosCount), testPosCount_(testPosCount)
 348 {
 349     // Generate an array of 0, 1, 2, ...
 350     // TODO: Make the tests work also with non-trivial exclusion IDs,
 351     // and test that.
 352     exclusionIds_.resize(std::max(refPosCount, testPosCount), 1);
 353     exclusionIds_[0] = 0;
 354     std::partial_sum(exclusionIds_.begin(), exclusionIds_.end(),
 355                      exclusionIds_.begin());
 356
 357     excls_.nr           = 0;
 358     excls_.index        = NULL;
 359     excls_.nra          = 0;
 360     excls_.a            = NULL;
 361     excls_.nalloc_index = 0;
 362     excls_.nalloc_a     = 0;
 363 }
 364
 365 void ExclusionsHelper::generateExclusions()
 366 {
 367     // TODO: Consider a better set of test data, where the density of the
 368     // particles would be higher, or where the exclusions would not be random,
 369     // to make a higher percentage of the exclusions to actually be within the
 370     // cutoff.
 371     exclsIndex_.reserve(testPosCount_ + 1);
 372     exclsAtoms_.reserve(testPosCount_ * 20);
 373     exclsIndex_.push_back(0);
 374     for (int i = 0; i < testPosCount_; ++i)
 375     {
 376         for (int j = 0; j < 20; ++j)
 377         {
 378             exclsAtoms_.push_back(i + j*3);
 379         }
 380         exclsIndex_.push_back(exclsAtoms_.size());
 381     }
 382     excls_.nr    = exclsIndex_.size();
 383     excls_.index = &exclsIndex_[0];
 384     excls_.nra   = exclsAtoms_.size();
 385     excls_.a     = &exclsAtoms_[0];
 386 }
 387
 388 /********************************************************************
 389  * NeighborhoodSearchTest
 390  */
 391
 392 class NeighborhoodSearchTest : public ::testing::Test
 393 {
 394     public:
 395         void testIsWithin(gmx::AnalysisNeighborhoodSearch  *search,
 396                           const NeighborhoodSearchTestData &data);
 397         void testMinimumDistance(gmx::AnalysisNeighborhoodSearch  *search,
 398                                  const NeighborhoodSearchTestData &data);
 399         void testNearestPoint(gmx::AnalysisNeighborhoodSearch  *search,
 400                               const NeighborhoodSearchTestData &data);
 401         void testPairSearch(gmx::AnalysisNeighborhoodSearch  *search,
 402                             const NeighborhoodSearchTestData &data);
 403         void testPairSearchFull(gmx::AnalysisNeighborhoodSearch          *search,
 404                                 const NeighborhoodSearchTestData         &data,
 405                                 const gmx::AnalysisNeighborhoodPositions &pos,
 406                                 const t_blocka                           *excls);
 407
 408         gmx::AnalysisNeighborhood        nb_;
 409 };
 410
 411 void NeighborhoodSearchTest::testIsWithin(
 412         gmx::AnalysisNeighborhoodSearch  *search,
 413         const NeighborhoodSearchTestData &data)
 414 {
 415     NeighborhoodSearchTestData::TestPositionList::const_iterator i;
 416     for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
 417     {
 418         const bool bWithin = (i->refMinDist <= data.cutoff_);
 419         EXPECT_EQ(bWithin, search->isWithin(i->x))
 420         << "Distance is " << i->refMinDist;
 421     }
 422 }
 423
 424 void NeighborhoodSearchTest::testMinimumDistance(
 425         gmx::AnalysisNeighborhoodSearch  *search,
 426         const NeighborhoodSearchTestData &data)
 427 {
 428     NeighborhoodSearchTestData::TestPositionList::const_iterator i;
 429     for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
 430     {
 431         const real refDist = i->refMinDist;
 432         EXPECT_REAL_EQ_TOL(refDist, search->minimumDistance(i->x),
 433                            gmx::test::ulpTolerance(20));
 434     }
 435 }
 436
 437 void NeighborhoodSearchTest::testNearestPoint(
 438         gmx::AnalysisNeighborhoodSearch  *search,
 439         const NeighborhoodSearchTestData &data)
 440 {
 441     NeighborhoodSearchTestData::TestPositionList::const_iterator i;
 442     for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
 443     {
 444         const gmx::AnalysisNeighborhoodPair pair = search->nearestPoint(i->x);
 445         if (pair.isValid())
 446         {
 447             EXPECT_EQ(i->refNearestPoint, pair.refIndex());
 448             EXPECT_EQ(0, pair.testIndex());
 449             EXPECT_REAL_EQ_TOL(i->refMinDist, sqrt(pair.distance2()),
 450                                gmx::test::ulpTolerance(64));
 451         }
 452         else
 453         {
 454             EXPECT_EQ(i->refNearestPoint, -1);
 455         }
 456     }
 457 }
 458
 459 /*! \brief
 460  * Helper function to check that all expected pairs were found.
 461  */
 462 static void checkAllPairsFound(const RefPairList &refPairs)
 463 {
 464     // This could be elegantly expressed with Google Mock matchers, but that
 465     // has a significant effect on the runtime of the tests...
 466     for (RefPairList::const_iterator i = refPairs.begin(); i != refPairs.end(); ++i)
 467     {
 468         if (!i->bFound)
 469         {
 470             ADD_FAILURE()
 471             << "Some pairs within the cutoff were not found.";
 472             break;
 473         }
 474     }
 475 }
 476
 477 void NeighborhoodSearchTest::testPairSearch(
 478         gmx::AnalysisNeighborhoodSearch  *search,
 479         const NeighborhoodSearchTestData &data)
 480 {
 481     testPairSearchFull(search, data, data.testPositions(), NULL);
 482 }
 483
 484 void NeighborhoodSearchTest::testPairSearchFull(
 485         gmx::AnalysisNeighborhoodSearch          *search,
 486         const NeighborhoodSearchTestData         &data,
 487         const gmx::AnalysisNeighborhoodPositions &pos,
 488         const t_blocka                           *excls)
 489 {
 490     // TODO: Some parts of this code do not work properly if pos does not
 491     // contain all the test positions.
 492     std::set<int> remainingTestPositions;
 493     for (size_t i = 0; i < data.testPositions_.size(); ++i)
 494     {
 495         remainingTestPositions.insert(i);
 496     }
 497     gmx::AnalysisNeighborhoodPairSearch pairSearch
 498         = search->startPairSearch(pos);
 499     gmx::AnalysisNeighborhoodPair       pair;
 500     // TODO: There is an ordering assumption here that may break in the future:
 501     // all pairs for a test position are assumed to be returned consencutively.
 502     RefPairList refPairs;
 503     int         prevTestPos = -1;
 504     while (pairSearch.findNextPair(&pair))
 505     {
 506         if (pair.testIndex() != prevTestPos)
 507         {
 508             if (prevTestPos != -1)
 509             {
 510                 checkAllPairsFound(refPairs);
 511             }
 512             const int testIndex = pair.testIndex();
 513             if (remainingTestPositions.count(testIndex) == 0)
 514             {
 515                 ADD_FAILURE()
 516                 << "Pairs for test position " << testIndex
 517                 << " are returned more than once.";
 518             }
 519             remainingTestPositions.erase(testIndex);
 520             refPairs = data.testPositions_[testIndex].refPairs;
 521             if (excls != NULL)
 522             {
 523                 ExclusionsHelper::markExcludedPairs(&refPairs, testIndex, excls);
 524             }
 525             prevTestPos = testIndex;
 526         }
 527
 528         NeighborhoodSearchTestData::RefPair searchPair(pair.refIndex(),
 529                                                        sqrt(pair.distance2()));
 530         RefPairList::iterator               foundRefPair
 531             = std::lower_bound(refPairs.begin(), refPairs.end(), searchPair);
 532         if (foundRefPair == refPairs.end() || foundRefPair->refIndex != pair.refIndex())
 533         {
 534             ADD_FAILURE()
 535             << "Expected: Pair (ref: " << pair.refIndex() << ", test: "
 536             << pair.testIndex() << ") is not within the cutoff.\n"
 537             << "  Actual: It is returned.";
 538         }
 539         else if (foundRefPair->bExcluded)
 540         {
 541             ADD_FAILURE()
 542             << "Expected: Pair (ref: " << pair.refIndex() << ", test: "
 543             << pair.testIndex() << ") is excluded from the search.\n"
 544             << "  Actual: It is returned.";
 545         }
 546         else if (foundRefPair->bFound)
 547         {
 548             ADD_FAILURE()
 549             << "Expected: Pair (ref: " << pair.refIndex() << ", test: "
 550             << pair.testIndex() << ") is returned only once.\n"
 551             << "  Actual: It is returned multiple times.";
 552         }
 553         else
 554         {
 555             foundRefPair->bFound = true;
 556             EXPECT_REAL_EQ_TOL(foundRefPair->distance, searchPair.distance,
 557                                gmx::test::ulpTolerance(64))
 558             << "Distance computed by the neighborhood search does not match.";
 559         }
 560     }
 561     checkAllPairsFound(refPairs);
 562     for (std::set<int>::const_iterator i = remainingTestPositions.begin();
 563          i != remainingTestPositions.end(); ++i)
 564     {
 565         if (!data.testPositions_[*i].refPairs.empty())
 566         {
 567             ADD_FAILURE()
 568             << "Expected: Pairs would be returned for test position " << *i << ".\n"
 569             << "  Actual: None were returned.";
 570             break;
 571         }
 572     }
 573 }
 574
 575 /********************************************************************
 576  * Test data generation
 577  */
 578
 579 class TrivialTestData
 580 {
 581     public:
 582         static const NeighborhoodSearchTestData &get()
 583         {
 584             static TrivialTestData singleton;
 585             return singleton.data_;
 586         }
 587
 588         TrivialTestData() : data_(12345, 1.0)
 589         {
 590             data_.box_[XX][XX] = 5.0;
 591             data_.box_[YY][YY] = 5.0;
 592             data_.box_[ZZ][ZZ] = 5.0;
 593             data_.generateRandomRefPositions(10);
 594             data_.generateRandomTestPositions(5);
 595             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 596             data_.computeReferences(&data_.pbc_);
 597         }
 598
 599     private:
 600         NeighborhoodSearchTestData data_;
 601 };
 602
 603 class RandomBoxFullPBCData
 604 {
 605     public:
 606         static const NeighborhoodSearchTestData &get()
 607         {
 608             static RandomBoxFullPBCData singleton;
 609             return singleton.data_;
 610         }
 611
 612         RandomBoxFullPBCData() : data_(12345, 1.0)
 613         {
 614             data_.box_[XX][XX] = 10.0;
 615             data_.box_[YY][YY] = 5.0;
 616             data_.box_[ZZ][ZZ] = 7.0;
 617             // TODO: Consider whether manually picking some positions would give better
 618             // test coverage.
 619             data_.generateRandomRefPositions(1000);
 620             data_.generateRandomTestPositions(100);
 621             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 622             data_.computeReferences(&data_.pbc_);
 623         }
 624
 625     private:
 626         NeighborhoodSearchTestData data_;
 627 };
 628
 629 class RandomBoxXYFullPBCData
 630 {
 631     public:
 632         static const NeighborhoodSearchTestData &get()
 633         {
 634             static RandomBoxXYFullPBCData singleton;
 635             return singleton.data_;
 636         }
 637
 638         RandomBoxXYFullPBCData() : data_(54321, 1.0)
 639         {
 640             data_.box_[XX][XX] = 10.0;
 641             data_.box_[YY][YY] = 5.0;
 642             data_.box_[ZZ][ZZ] = 7.0;
 643             // TODO: Consider whether manually picking some positions would give better
 644             // test coverage.
 645             data_.generateRandomRefPositions(1000);
 646             data_.generateRandomTestPositions(100);
 647             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 648             data_.computeReferencesXY(&data_.pbc_);
 649         }
 650
 651     private:
 652         NeighborhoodSearchTestData data_;
 653 };
 654
 655 class RandomTriclinicFullPBCData
 656 {
 657     public:
 658         static const NeighborhoodSearchTestData &get()
 659         {
 660             static RandomTriclinicFullPBCData singleton;
 661             return singleton.data_;
 662         }
 663
 664         RandomTriclinicFullPBCData() : data_(12345, 1.0)
 665         {
 666             data_.box_[XX][XX] = 5.0;
 667             data_.box_[YY][XX] = 2.5;
 668             data_.box_[YY][YY] = 2.5*sqrt(3.0);
 669             data_.box_[ZZ][XX] = 2.5;
 670             data_.box_[ZZ][YY] = 2.5*sqrt(1.0/3.0);
 671             data_.box_[ZZ][ZZ] = 5.0*sqrt(2.0/3.0);
 672             // TODO: Consider whether manually picking some positions would give better
 673             // test coverage.
 674             data_.generateRandomRefPositions(1000);
 675             data_.generateRandomTestPositions(100);
 676             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 677             data_.computeReferences(&data_.pbc_);
 678         }
 679
 680     private:
 681         NeighborhoodSearchTestData data_;
 682 };
 683
 684 class RandomBox2DPBCData
 685 {
 686     public:
 687         static const NeighborhoodSearchTestData &get()
 688         {
 689             static RandomBox2DPBCData singleton;
 690             return singleton.data_;
 691         }
 692
 693         RandomBox2DPBCData() : data_(12345, 1.0)
 694         {
 695             data_.box_[XX][XX] = 10.0;
 696             data_.box_[YY][YY] = 7.0;
 697             data_.box_[ZZ][ZZ] = 5.0;
 698             // TODO: Consider whether manually picking some positions would give better
 699             // test coverage.
 700             data_.generateRandomRefPositions(1000);
 701             data_.generateRandomTestPositions(100);
 702             set_pbc(&data_.pbc_, epbcXY, data_.box_);
 703             data_.computeReferences(&data_.pbc_);
 704         }
 705
 706     private:
 707         NeighborhoodSearchTestData data_;
 708 };
 709
 710 /********************************************************************
 711  * Actual tests
 712  */
 713
 714 TEST_F(NeighborhoodSearchTest, SimpleSearch)
 715 {
 716     const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
 717
 718     nb_.setCutoff(data.cutoff_);
 719     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Simple);
 720     gmx::AnalysisNeighborhoodSearch search =
 721         nb_.initSearch(&data.pbc_, data.refPositions());
 722     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
 723
 724     testIsWithin(&search, data);
 725     testMinimumDistance(&search, data);
 726     testNearestPoint(&search, data);
 727     testPairSearch(&search, data);
 728 }
 729
 730 TEST_F(NeighborhoodSearchTest, GridSearchBox)
 731 {
 732     const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
 733
 734     nb_.setCutoff(data.cutoff_);
 735     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
 736     gmx::AnalysisNeighborhoodSearch search =
 737         nb_.initSearch(&data.pbc_, data.refPositions());
 738     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
 739
 740     testIsWithin(&search, data);
 741     testMinimumDistance(&search, data);
 742     testNearestPoint(&search, data);
 743     testPairSearch(&search, data);
 744 }
 745
 746 TEST_F(NeighborhoodSearchTest, GridSearchTriclinic)
 747 {
 748     const NeighborhoodSearchTestData &data = RandomTriclinicFullPBCData::get();
 749
 750     nb_.setCutoff(data.cutoff_);
 751     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
 752     gmx::AnalysisNeighborhoodSearch search =
 753         nb_.initSearch(&data.pbc_, data.refPositions());
 754     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
 755
 756     testPairSearch(&search, data);
 757 }
 758
 759 TEST_F(NeighborhoodSearchTest, GridSearch2DPBC)
 760 {
 761     const NeighborhoodSearchTestData &data = RandomBox2DPBCData::get();
 762
 763     nb_.setCutoff(data.cutoff_);
 764     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
 765     gmx::AnalysisNeighborhoodSearch search =
 766         nb_.initSearch(&data.pbc_, data.refPositions());
 767     // Currently, grid searching not supported with 2D PBC.
 768     //ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
 769
 770     testIsWithin(&search, data);
 771     testMinimumDistance(&search, data);
 772     testNearestPoint(&search, data);
 773     testPairSearch(&search, data);
 774 }
 775
 776 TEST_F(NeighborhoodSearchTest, GridSearchXYBox)
 777 {
 778     const NeighborhoodSearchTestData &data = RandomBoxXYFullPBCData::get();
 779
 780     nb_.setCutoff(data.cutoff_);
 781     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
 782     nb_.setXYMode(true);
 783     gmx::AnalysisNeighborhoodSearch search =
 784         nb_.initSearch(&data.pbc_, data.refPositions());
 785     // Currently, grid searching not supported with XY.
 786     //ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
 787
 788     testIsWithin(&search, data);
 789     testMinimumDistance(&search, data);
 790     testNearestPoint(&search, data);
 791     testPairSearch(&search, data);
 792 }
 793
 794 TEST_F(NeighborhoodSearchTest, HandlesConcurrentSearches)
 795 {
 796     const NeighborhoodSearchTestData &data = TrivialTestData::get();
 797
 798     nb_.setCutoff(data.cutoff_);
 799     gmx::AnalysisNeighborhoodSearch search1 =
 800         nb_.initSearch(&data.pbc_, data.refPositions());
 801     gmx::AnalysisNeighborhoodSearch search2 =
 802         nb_.initSearch(&data.pbc_, data.refPositions());
 803
 804     gmx::AnalysisNeighborhoodPairSearch pairSearch1 =
 805         search1.startPairSearch(data.testPosition(0));
 806     gmx::AnalysisNeighborhoodPairSearch pairSearch2 =
 807         search1.startPairSearch(data.testPosition(1));
 808
 809     testPairSearch(&search2, data);
 810
 811     gmx::AnalysisNeighborhoodPair pair;
 812     ASSERT_TRUE(pairSearch1.findNextPair(&pair))
 813     << "Test data did not contain any pairs for position 0 (problem in the test).";
 814     EXPECT_EQ(0, pair.testIndex());
 815     {
 816         NeighborhoodSearchTestData::RefPair searchPair(pair.refIndex(), sqrt(pair.distance2()));
 817         EXPECT_TRUE(data.containsPair(0, searchPair));
 818     }
 819
 820     ASSERT_TRUE(pairSearch2.findNextPair(&pair))
 821     << "Test data did not contain any pairs for position 1 (problem in the test).";
 822     EXPECT_EQ(1, pair.testIndex());
 823     {
 824         NeighborhoodSearchTestData::RefPair searchPair(pair.refIndex(), sqrt(pair.distance2()));
 825         EXPECT_TRUE(data.containsPair(1, searchPair));
 826     }
 827 }
 828
 829 TEST_F(NeighborhoodSearchTest, HandlesSkippingPairs)
 830 {
 831     const NeighborhoodSearchTestData &data = TrivialTestData::get();
 832
 833     nb_.setCutoff(data.cutoff_);
 834     gmx::AnalysisNeighborhoodSearch     search =
 835         nb_.initSearch(&data.pbc_, data.refPositions());
 836     gmx::AnalysisNeighborhoodPairSearch pairSearch =
 837         search.startPairSearch(data.testPositions());
 838     gmx::AnalysisNeighborhoodPair       pair;
 839     // TODO: This test needs to be adjusted if the grid search gets optimized
 840     // to loop over the test positions in cell order (first, the ordering
 841     // assumption here breaks, and second, it then needs to be tested
 842     // separately for simple and grid searches).
 843     int currentIndex = 0;
 844     while (pairSearch.findNextPair(&pair))
 845     {
 846         while (currentIndex < pair.testIndex())
 847         {
 848             ++currentIndex;
 849         }
 850         EXPECT_EQ(currentIndex, pair.testIndex());
 851         NeighborhoodSearchTestData::RefPair searchPair(pair.refIndex(), sqrt(pair.distance2()));
 852         EXPECT_TRUE(data.containsPair(currentIndex, searchPair));
 853         pairSearch.skipRemainingPairsForTestPosition();
 854         ++currentIndex;
 855     }
 856 }
 857
 858 TEST_F(NeighborhoodSearchTest, SimpleSearchExclusions)
 859 {
 860     const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
 861
 862     ExclusionsHelper                  helper(data.refPosCount_, data.testPositions_.size());
 863     helper.generateExclusions();
 864
 865     nb_.setCutoff(data.cutoff_);
 866     nb_.setTopologyExclusions(helper.exclusions());
 867     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Simple);
 868     gmx::AnalysisNeighborhoodSearch search =
 869         nb_.initSearch(&data.pbc_,
 870                        data.refPositions().exclusionIds(helper.refPosIds()));
 871     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
 872
 873     testPairSearchFull(&search, data,
 874                        data.testPositions().exclusionIds(helper.testPosIds()),
 875                        helper.exclusions());
 876 }
 877
 878 TEST_F(NeighborhoodSearchTest, GridSearchExclusions)
 879 {
 880     const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
 881
 882     ExclusionsHelper                  helper(data.refPosCount_, data.testPositions_.size());
 883     helper.generateExclusions();
 884
 885     nb_.setCutoff(data.cutoff_);
 886     nb_.setTopologyExclusions(helper.exclusions());
 887     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
 888     gmx::AnalysisNeighborhoodSearch search =
 889         nb_.initSearch(&data.pbc_,
 890                        data.refPositions().exclusionIds(helper.refPosIds()));
 891     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
 892
 893     testPairSearchFull(&search, data,
 894                        data.testPositions().exclusionIds(helper.testPosIds()),
 895                        helper.exclusions());
 896 }
 897
 898 } // namespace