src/gromacs/selection/tests/nbsearch.cpp

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