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