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