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