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 <gtest/gtest.h>
  48
  49 #include <cmath>
  50
  51 #include <limits>
  52 #include <set>
  53 #include <vector>
  54
  55 #include "gromacs/legacyheaders/pbc.h"
  56 #include "gromacs/legacyheaders/vec.h"
  57
  58 #include "gromacs/selection/nbsearch.h"
  59 #include "gromacs/random/random.h"
  60 #include "gromacs/utility/smalloc.h"
  61
  62 #include "testutils/testasserts.h"
  63
  64 namespace
  65 {
  66
  67 /********************************************************************
  68  * NeighborhoodSearchTestData
  69  */
  70
  71 class NeighborhoodSearchTestData
  72 {
  73     public:
  74         struct TestPosition
  75         {
  76             explicit TestPosition(const rvec x)
  77                 : refMinDist(0.0), refNearestPoint(-1)
  78             {
  79                 copy_rvec(x, this->x);
  80             }
  81
  82             rvec                x;
  83             real                refMinDist;
  84             int                 refNearestPoint;
  85             std::set<int>       refPairs;
  86         };
  87         typedef std::vector<TestPosition> TestPositionList;
  88
  89         NeighborhoodSearchTestData(int seed, real cutoff);
  90         ~NeighborhoodSearchTestData();
  91
  92         gmx::AnalysisNeighborhoodPositions refPositions() const
  93         {
  94             return gmx::AnalysisNeighborhoodPositions(refPos_, refPosCount_);
  95         }
  96         gmx::AnalysisNeighborhoodPositions testPositions() const
  97         {
  98             if (testPos_ == NULL)
  99             {
 100                 snew(testPos_, testPositions_.size());
 101                 for (size_t i = 0; i < testPositions_.size(); ++i)
 102                 {
 103                     copy_rvec(testPositions_[i].x, testPos_[i]);
 104                 }
 105             }
 106             return gmx::AnalysisNeighborhoodPositions(testPos_,
 107                                                       testPositions_.size());
 108         }
 109         gmx::AnalysisNeighborhoodPositions testPosition(int index) const
 110         {
 111             return testPositions().selectSingleFromArray(index);
 112         }
 113
 114         void addTestPosition(const rvec x)
 115         {
 116             GMX_RELEASE_ASSERT(testPos_ == NULL,
 117                                "Cannot add positions after testPositions() call");
 118             testPositions_.push_back(TestPosition(x));
 119         }
 120         void generateRandomPosition(rvec x);
 121         void generateRandomRefPositions(int count);
 122         void generateRandomTestPositions(int count);
 123         void computeReferences(t_pbc *pbc);
 124
 125         gmx_rng_t                        rng_;
 126         real                             cutoff_;
 127         matrix                           box_;
 128         t_pbc                            pbc_;
 129         int                              refPosCount_;
 130         rvec                            *refPos_;
 131         TestPositionList                 testPositions_;
 132
 133     private:
 134         mutable rvec                    *testPos_;
 135 };
 136
 137 NeighborhoodSearchTestData::NeighborhoodSearchTestData(int seed, real cutoff)
 138     : rng_(NULL), cutoff_(cutoff), refPosCount_(0), refPos_(NULL), testPos_(NULL)
 139 {
 140     // TODO: Handle errors.
 141     rng_ = gmx_rng_init(seed);
 142     clear_mat(box_);
 143     set_pbc(&pbc_, epbcNONE, box_);
 144 }
 145
 146 NeighborhoodSearchTestData::~NeighborhoodSearchTestData()
 147 {
 148     if (rng_ != NULL)
 149     {
 150         gmx_rng_destroy(rng_);
 151     }
 152     sfree(refPos_);
 153     sfree(testPos_);
 154 }
 155
 156 void NeighborhoodSearchTestData::generateRandomPosition(rvec x)
 157 {
 158     rvec fx;
 159     fx[XX] = gmx_rng_uniform_real(rng_);
 160     fx[YY] = gmx_rng_uniform_real(rng_);
 161     fx[ZZ] = gmx_rng_uniform_real(rng_);
 162     mvmul(box_, fx, x);
 163     // Add a small displacement to allow positions outside the box
 164     x[XX] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
 165     x[YY] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
 166     x[ZZ] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
 167 }
 168
 169 void NeighborhoodSearchTestData::generateRandomRefPositions(int count)
 170 {
 171     refPosCount_ = count;
 172     snew(refPos_, refPosCount_);
 173     for (int i = 0; i < refPosCount_; ++i)
 174     {
 175         generateRandomPosition(refPos_[i]);
 176     }
 177 }
 178
 179 void NeighborhoodSearchTestData::generateRandomTestPositions(int count)
 180 {
 181     testPositions_.reserve(count);
 182     for (int i = 0; i < count; ++i)
 183     {
 184         rvec x;
 185         generateRandomPosition(x);
 186         addTestPosition(x);
 187     }
 188 }
 189
 190 void NeighborhoodSearchTestData::computeReferences(t_pbc *pbc)
 191 {
 192     real cutoff = cutoff_;
 193     if (cutoff <= 0)
 194     {
 195         cutoff = std::numeric_limits<real>::max();
 196     }
 197     TestPositionList::iterator i;
 198     for (i = testPositions_.begin(); i != testPositions_.end(); ++i)
 199     {
 200         i->refMinDist      = cutoff;
 201         i->refNearestPoint = -1;
 202         i->refPairs.clear();
 203         for (int j = 0; j < refPosCount_; ++j)
 204         {
 205             rvec dx;
 206             if (pbc != NULL)
 207             {
 208                 pbc_dx(pbc, i->x, refPos_[j], dx);
 209             }
 210             else
 211             {
 212                 rvec_sub(i->x, refPos_[j], dx);
 213             }
 214             const real dist = norm(dx);
 215             if (dist < i->refMinDist)
 216             {
 217                 i->refMinDist      = dist;
 218                 i->refNearestPoint = j;
 219             }
 220             if (dist <= cutoff)
 221             {
 222                 i->refPairs.insert(j);
 223             }
 224         }
 225     }
 226 }
 227
 228 /********************************************************************
 229  * NeighborhoodSearchTest
 230  */
 231
 232 class NeighborhoodSearchTest : public ::testing::Test
 233 {
 234     public:
 235         void testIsWithin(gmx::AnalysisNeighborhoodSearch  *search,
 236                           const NeighborhoodSearchTestData &data);
 237         void testMinimumDistance(gmx::AnalysisNeighborhoodSearch  *search,
 238                                  const NeighborhoodSearchTestData &data);
 239         void testNearestPoint(gmx::AnalysisNeighborhoodSearch  *search,
 240                               const NeighborhoodSearchTestData &data);
 241         void testPairSearch(gmx::AnalysisNeighborhoodSearch  *search,
 242                             const NeighborhoodSearchTestData &data);
 243
 244         gmx::AnalysisNeighborhood        nb_;
 245 };
 246
 247 void NeighborhoodSearchTest::testIsWithin(
 248         gmx::AnalysisNeighborhoodSearch  *search,
 249         const NeighborhoodSearchTestData &data)
 250 {
 251     NeighborhoodSearchTestData::TestPositionList::const_iterator i;
 252     for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
 253     {
 254         const bool bWithin = (i->refMinDist <= data.cutoff_);
 255         EXPECT_EQ(bWithin, search->isWithin(i->x))
 256         << "Distance is " << i->refMinDist;
 257     }
 258 }
 259
 260 void NeighborhoodSearchTest::testMinimumDistance(
 261         gmx::AnalysisNeighborhoodSearch  *search,
 262         const NeighborhoodSearchTestData &data)
 263 {
 264     NeighborhoodSearchTestData::TestPositionList::const_iterator i;
 265     for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
 266     {
 267         const real refDist = i->refMinDist;
 268         EXPECT_REAL_EQ_TOL(refDist, search->minimumDistance(i->x),
 269                            gmx::test::ulpTolerance(20));
 270     }
 271 }
 272
 273 void NeighborhoodSearchTest::testNearestPoint(
 274         gmx::AnalysisNeighborhoodSearch  *search,
 275         const NeighborhoodSearchTestData &data)
 276 {
 277     NeighborhoodSearchTestData::TestPositionList::const_iterator i;
 278     for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
 279     {
 280         const gmx::AnalysisNeighborhoodPair pair = search->nearestPoint(i->x);
 281         if (pair.isValid())
 282         {
 283             EXPECT_EQ(i->refNearestPoint, pair.refIndex());
 284             EXPECT_EQ(0, pair.testIndex());
 285         }
 286         else
 287         {
 288             EXPECT_EQ(i->refNearestPoint, -1);
 289         }
 290     }
 291 }
 292
 293 void NeighborhoodSearchTest::testPairSearch(
 294         gmx::AnalysisNeighborhoodSearch  *search,
 295         const NeighborhoodSearchTestData &data)
 296 {
 297     NeighborhoodSearchTestData::TestPositionList::const_iterator i;
 298     for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
 299     {
 300         std::set<int> checkSet                         = i->refPairs;
 301         gmx::AnalysisNeighborhoodPairSearch pairSearch =
 302             search->startPairSearch(i->x);
 303         gmx::AnalysisNeighborhoodPair       pair;
 304         while (pairSearch.findNextPair(&pair))
 305         {
 306             EXPECT_EQ(0, pair.testIndex());
 307             if (checkSet.erase(pair.refIndex()) == 0)
 308             {
 309                 // TODO: Check whether the same pair was returned more than
 310                 // once and give a better error message if so.
 311                 ADD_FAILURE()
 312                 << "Expected: Position " << pair.refIndex()
 313                 << " is within cutoff.\n"
 314                 << "  Actual: It is not.";
 315             }
 316         }
 317         EXPECT_TRUE(checkSet.empty()) << "Some positions were not returned by the pair search.";
 318     }
 319 }
 320
 321 /********************************************************************
 322  * Test data generation
 323  */
 324
 325 class TrivialTestData
 326 {
 327     public:
 328         static const NeighborhoodSearchTestData &get()
 329         {
 330             static TrivialTestData singleton;
 331             return singleton.data_;
 332         }
 333
 334         TrivialTestData() : data_(12345, 1.0)
 335         {
 336             data_.box_[XX][XX] = 5.0;
 337             data_.box_[YY][YY] = 5.0;
 338             data_.box_[ZZ][ZZ] = 5.0;
 339             data_.generateRandomRefPositions(10);
 340             data_.generateRandomTestPositions(5);
 341             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 342             data_.computeReferences(&data_.pbc_);
 343         }
 344
 345     private:
 346         NeighborhoodSearchTestData data_;
 347 };
 348
 349 class RandomBoxFullPBCData
 350 {
 351     public:
 352         static const NeighborhoodSearchTestData &get()
 353         {
 354             static RandomBoxFullPBCData singleton;
 355             return singleton.data_;
 356         }
 357
 358         RandomBoxFullPBCData() : data_(12345, 1.0)
 359         {
 360             data_.box_[XX][XX] = 10.0;
 361             data_.box_[YY][YY] = 5.0;
 362             data_.box_[ZZ][ZZ] = 7.0;
 363             // TODO: Consider whether manually picking some positions would give better
 364             // test coverage.
 365             data_.generateRandomRefPositions(1000);
 366             data_.generateRandomTestPositions(100);
 367             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 368             data_.computeReferences(&data_.pbc_);
 369         }
 370
 371     private:
 372         NeighborhoodSearchTestData data_;
 373 };
 374
 375 class RandomTriclinicFullPBCData
 376 {
 377     public:
 378         static const NeighborhoodSearchTestData &get()
 379         {
 380             static RandomTriclinicFullPBCData singleton;
 381             return singleton.data_;
 382         }
 383
 384         RandomTriclinicFullPBCData() : data_(12345, 1.0)
 385         {
 386             data_.box_[XX][XX] = 5.0;
 387             data_.box_[YY][XX] = 2.5;
 388             data_.box_[YY][YY] = 2.5*sqrt(3.0);
 389             data_.box_[ZZ][XX] = 2.5;
 390             data_.box_[ZZ][YY] = 2.5*sqrt(1.0/3.0);
 391             data_.box_[ZZ][ZZ] = 5.0*sqrt(2.0/3.0);
 392             // TODO: Consider whether manually picking some positions would give better
 393             // test coverage.
 394             data_.generateRandomRefPositions(1000);
 395             data_.generateRandomTestPositions(100);
 396             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 397             data_.computeReferences(&data_.pbc_);
 398         }
 399
 400     private:
 401         NeighborhoodSearchTestData data_;
 402 };
 403
 404 class RandomBox2DPBCData
 405 {
 406     public:
 407         static const NeighborhoodSearchTestData &get()
 408         {
 409             static RandomBox2DPBCData singleton;
 410             return singleton.data_;
 411         }
 412
 413         RandomBox2DPBCData() : data_(12345, 1.0)
 414         {
 415             data_.box_[XX][XX] = 10.0;
 416             data_.box_[YY][YY] = 7.0;
 417             data_.box_[ZZ][ZZ] = 5.0;
 418             // TODO: Consider whether manually picking some positions would give better
 419             // test coverage.
 420             data_.generateRandomRefPositions(1000);
 421             data_.generateRandomTestPositions(100);
 422             set_pbc(&data_.pbc_, epbcXY, data_.box_);
 423             data_.computeReferences(&data_.pbc_);
 424         }
 425
 426     private:
 427         NeighborhoodSearchTestData data_;
 428 };
 429
 430 /********************************************************************
 431  * Actual tests
 432  */
 433
 434 TEST_F(NeighborhoodSearchTest, SimpleSearch)
 435 {
 436     const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
 437
 438     nb_.setCutoff(data.cutoff_);
 439     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Simple);
 440     gmx::AnalysisNeighborhoodSearch search =
 441         nb_.initSearch(&data.pbc_, data.refPositions());
 442     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
 443
 444     testIsWithin(&search, data);
 445     testMinimumDistance(&search, data);
 446     testNearestPoint(&search, data);
 447     testPairSearch(&search, data);
 448 }
 449
 450 TEST_F(NeighborhoodSearchTest, GridSearchBox)
 451 {
 452     const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
 453
 454     nb_.setCutoff(data.cutoff_);
 455     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
 456     gmx::AnalysisNeighborhoodSearch search =
 457         nb_.initSearch(&data.pbc_, data.refPositions());
 458     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
 459
 460     testIsWithin(&search, data);
 461     testMinimumDistance(&search, data);
 462     testNearestPoint(&search, data);
 463     testPairSearch(&search, data);
 464 }
 465
 466 TEST_F(NeighborhoodSearchTest, GridSearchTriclinic)
 467 {
 468     const NeighborhoodSearchTestData &data = RandomTriclinicFullPBCData::get();
 469
 470     nb_.setCutoff(data.cutoff_);
 471     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
 472     gmx::AnalysisNeighborhoodSearch search =
 473         nb_.initSearch(&data.pbc_, data.refPositions());
 474     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
 475
 476     testPairSearch(&search, data);
 477 }
 478
 479 TEST_F(NeighborhoodSearchTest, GridSearch2DPBC)
 480 {
 481     const NeighborhoodSearchTestData &data = RandomBox2DPBCData::get();
 482
 483     nb_.setCutoff(data.cutoff_);
 484     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
 485     gmx::AnalysisNeighborhoodSearch search =
 486         nb_.initSearch(&data.pbc_, data.refPositions());
 487     // Currently, grid searching not supported with 2D PBC.
 488     //ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
 489
 490     testIsWithin(&search, data);
 491     testMinimumDistance(&search, data);
 492     testNearestPoint(&search, data);
 493     testPairSearch(&search, data);
 494 }
 495
 496 TEST_F(NeighborhoodSearchTest, HandlesConcurrentSearches)
 497 {
 498     const NeighborhoodSearchTestData &data = TrivialTestData::get();
 499
 500     nb_.setCutoff(data.cutoff_);
 501     gmx::AnalysisNeighborhoodSearch search1 =
 502         nb_.initSearch(&data.pbc_, data.refPositions());
 503     gmx::AnalysisNeighborhoodSearch search2 =
 504         nb_.initSearch(&data.pbc_, data.refPositions());
 505
 506     gmx::AnalysisNeighborhoodPairSearch pairSearch1 =
 507         search1.startPairSearch(data.testPosition(0));
 508     gmx::AnalysisNeighborhoodPairSearch pairSearch2 =
 509         search1.startPairSearch(data.testPosition(1));
 510
 511     testPairSearch(&search2, data);
 512
 513     gmx::AnalysisNeighborhoodPair pair;
 514     pairSearch1.findNextPair(&pair);
 515     EXPECT_EQ(0, pair.testIndex());
 516     EXPECT_TRUE(data.testPositions_[0].refPairs.count(pair.refIndex()) == 1);
 517
 518     pairSearch2.findNextPair(&pair);
 519     EXPECT_EQ(1, pair.testIndex());
 520     EXPECT_TRUE(data.testPositions_[1].refPairs.count(pair.refIndex()) == 1);
 521 }
 522
 523 TEST_F(NeighborhoodSearchTest, HandlesSkippingPairs)
 524 {
 525     const NeighborhoodSearchTestData &data = TrivialTestData::get();
 526
 527     nb_.setCutoff(data.cutoff_);
 528     gmx::AnalysisNeighborhoodSearch     search =
 529         nb_.initSearch(&data.pbc_, data.refPositions());
 530     gmx::AnalysisNeighborhoodPairSearch pairSearch =
 531         search.startPairSearch(data.testPositions());
 532     gmx::AnalysisNeighborhoodPair       pair;
 533     // TODO: This test needs to be adjusted if the grid search gets optimized
 534     // to loop over the test positions in cell order (first, the ordering
 535     // assumption here breaks, and second, it then needs to be tested
 536     // separately for simple and grid searches).
 537     int currentIndex = 0;
 538     while (pairSearch.findNextPair(&pair))
 539     {
 540         while (currentIndex < pair.testIndex())
 541         {
 542             ++currentIndex;
 543         }
 544         EXPECT_EQ(currentIndex, pair.testIndex());
 545         EXPECT_TRUE(data.testPositions_[currentIndex].refPairs.count(pair.refIndex()) == 1);
 546         pairSearch.skipRemainingPairsForTestPosition();
 547         ++currentIndex;
 548     }
 549 }
 550
 551 } // namespace