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             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_REAL_EQ_TOL(refDist, search->minimumDistance(i->x),
 273                            gmx::test::ulpTolerance(20));
 274     }
 275 }
 276
 277 void NeighborhoodSearchTest::testNearestPoint(
 278         gmx::AnalysisNeighborhoodSearch  *search,
 279         const NeighborhoodSearchTestData &data)
 280 {
 281     NeighborhoodSearchTestData::TestPositionList::const_iterator i;
 282     for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
 283     {
 284         const gmx::AnalysisNeighborhoodPair pair = search->nearestPoint(i->x);
 285         if (pair.isValid())
 286         {
 287             EXPECT_EQ(i->refNearestPoint, pair.refIndex());
 288             EXPECT_EQ(0, pair.testIndex());
 289         }
 290         else
 291         {
 292             EXPECT_EQ(i->refNearestPoint, -1);
 293         }
 294     }
 295 }
 296
 297 void NeighborhoodSearchTest::testPairSearch(
 298         gmx::AnalysisNeighborhoodSearch  *search,
 299         const NeighborhoodSearchTestData &data)
 300 {
 301     NeighborhoodSearchTestData::TestPositionList::const_iterator i;
 302     for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
 303     {
 304         std::set<int> checkSet                         = i->refPairs;
 305         gmx::AnalysisNeighborhoodPairSearch pairSearch =
 306             search->startPairSearch(i->x);
 307         gmx::AnalysisNeighborhoodPair       pair;
 308         while (pairSearch.findNextPair(&pair))
 309         {
 310             EXPECT_EQ(0, pair.testIndex());
 311             if (checkSet.erase(pair.refIndex()) == 0)
 312             {
 313                 // TODO: Check whether the same pair was returned more than
 314                 // once and give a better error message if so.
 315                 ADD_FAILURE()
 316                 << "Expected: Position " << pair.refIndex()
 317                 << " is within cutoff.\n"
 318                 << "  Actual: It is not.";
 319             }
 320         }
 321         EXPECT_TRUE(checkSet.empty()) << "Some positions were not returned by the pair search.";
 322     }
 323 }
 324
 325 /********************************************************************
 326  * Test data generation
 327  */
 328
 329 class TrivialTestData
 330 {
 331     public:
 332         static const NeighborhoodSearchTestData &get()
 333         {
 334             static TrivialTestData singleton;
 335             return singleton.data_;
 336         }
 337
 338         TrivialTestData() : data_(12345, 1.0)
 339         {
 340             data_.box_[XX][XX] = 5.0;
 341             data_.box_[YY][YY] = 5.0;
 342             data_.box_[ZZ][ZZ] = 5.0;
 343             data_.generateRandomRefPositions(10);
 344             data_.generateRandomTestPositions(5);
 345             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 346             data_.computeReferences(&data_.pbc_);
 347         }
 348
 349     private:
 350         NeighborhoodSearchTestData data_;
 351 };
 352
 353 class RandomBoxFullPBCData
 354 {
 355     public:
 356         static const NeighborhoodSearchTestData &get()
 357         {
 358             static RandomBoxFullPBCData singleton;
 359             return singleton.data_;
 360         }
 361
 362         RandomBoxFullPBCData() : data_(12345, 1.0)
 363         {
 364             data_.box_[XX][XX] = 10.0;
 365             data_.box_[YY][YY] = 5.0;
 366             data_.box_[ZZ][ZZ] = 7.0;
 367             // TODO: Consider whether manually picking some positions would give better
 368             // test coverage.
 369             data_.generateRandomRefPositions(1000);
 370             data_.generateRandomTestPositions(100);
 371             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 372             data_.computeReferences(&data_.pbc_);
 373         }
 374
 375     private:
 376         NeighborhoodSearchTestData data_;
 377 };
 378
 379 class RandomTriclinicFullPBCData
 380 {
 381     public:
 382         static const NeighborhoodSearchTestData &get()
 383         {
 384             static RandomTriclinicFullPBCData singleton;
 385             return singleton.data_;
 386         }
 387
 388         RandomTriclinicFullPBCData() : data_(12345, 1.0)
 389         {
 390             data_.box_[XX][XX] = 5.0;
 391             data_.box_[YY][XX] = 2.5;
 392             data_.box_[YY][YY] = 2.5*sqrt(3.0);
 393             data_.box_[ZZ][XX] = 2.5;
 394             data_.box_[ZZ][YY] = 2.5*sqrt(1.0/3.0);
 395             data_.box_[ZZ][ZZ] = 5.0*sqrt(2.0/3.0);
 396             // TODO: Consider whether manually picking some positions would give better
 397             // test coverage.
 398             data_.generateRandomRefPositions(1000);
 399             data_.generateRandomTestPositions(100);
 400             set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
 401             data_.computeReferences(&data_.pbc_);
 402         }
 403
 404     private:
 405         NeighborhoodSearchTestData data_;
 406 };
 407
 408 class RandomBox2DPBCData
 409 {
 410     public:
 411         static const NeighborhoodSearchTestData &get()
 412         {
 413             static RandomBox2DPBCData singleton;
 414             return singleton.data_;
 415         }
 416
 417         RandomBox2DPBCData() : data_(12345, 1.0)
 418         {
 419             data_.box_[XX][XX] = 10.0;
 420             data_.box_[YY][YY] = 7.0;
 421             data_.box_[ZZ][ZZ] = 5.0;
 422             // TODO: Consider whether manually picking some positions would give better
 423             // test coverage.
 424             data_.generateRandomRefPositions(1000);
 425             data_.generateRandomTestPositions(100);
 426             set_pbc(&data_.pbc_, epbcXY, data_.box_);
 427             data_.computeReferences(&data_.pbc_);
 428         }
 429
 430     private:
 431         NeighborhoodSearchTestData data_;
 432 };
 433
 434 /********************************************************************
 435  * Actual tests
 436  */
 437
 438 TEST_F(NeighborhoodSearchTest, SimpleSearch)
 439 {
 440     const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
 441
 442     nb_.setCutoff(data.cutoff_);
 443     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Simple);
 444     gmx::AnalysisNeighborhoodSearch search =
 445         nb_.initSearch(&data.pbc_, data.refPositions());
 446     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
 447
 448     testIsWithin(&search, data);
 449     testMinimumDistance(&search, data);
 450     testNearestPoint(&search, data);
 451     testPairSearch(&search, data);
 452 }
 453
 454 TEST_F(NeighborhoodSearchTest, GridSearchBox)
 455 {
 456     const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
 457
 458     nb_.setCutoff(data.cutoff_);
 459     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
 460     gmx::AnalysisNeighborhoodSearch search =
 461         nb_.initSearch(&data.pbc_, data.refPositions());
 462     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
 463
 464     testIsWithin(&search, data);
 465     testMinimumDistance(&search, data);
 466     testNearestPoint(&search, data);
 467     testPairSearch(&search, data);
 468 }
 469
 470 TEST_F(NeighborhoodSearchTest, GridSearchTriclinic)
 471 {
 472     const NeighborhoodSearchTestData &data = RandomTriclinicFullPBCData::get();
 473
 474     nb_.setCutoff(data.cutoff_);
 475     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
 476     gmx::AnalysisNeighborhoodSearch search =
 477         nb_.initSearch(&data.pbc_, data.refPositions());
 478     ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
 479
 480     testPairSearch(&search, data);
 481 }
 482
 483 TEST_F(NeighborhoodSearchTest, GridSearch2DPBC)
 484 {
 485     const NeighborhoodSearchTestData &data = RandomBox2DPBCData::get();
 486
 487     nb_.setCutoff(data.cutoff_);
 488     nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
 489     gmx::AnalysisNeighborhoodSearch search =
 490         nb_.initSearch(&data.pbc_, data.refPositions());
 491     // Currently, grid searching not supported with 2D PBC.
 492     //ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
 493
 494     testIsWithin(&search, data);
 495     testMinimumDistance(&search, data);
 496     testNearestPoint(&search, data);
 497     testPairSearch(&search, data);
 498 }
 499
 500 TEST_F(NeighborhoodSearchTest, HandlesConcurrentSearches)
 501 {
 502     const NeighborhoodSearchTestData &data = TrivialTestData::get();
 503
 504     nb_.setCutoff(data.cutoff_);
 505     gmx::AnalysisNeighborhoodSearch search1 =
 506         nb_.initSearch(&data.pbc_, data.refPositions());
 507     gmx::AnalysisNeighborhoodSearch search2 =
 508         nb_.initSearch(&data.pbc_, data.refPositions());
 509
 510     gmx::AnalysisNeighborhoodPairSearch pairSearch1 =
 511         search1.startPairSearch(data.testPosition(0));
 512     gmx::AnalysisNeighborhoodPairSearch pairSearch2 =
 513         search1.startPairSearch(data.testPosition(1));
 514
 515     testPairSearch(&search2, data);
 516
 517     gmx::AnalysisNeighborhoodPair pair;
 518     pairSearch1.findNextPair(&pair);
 519     EXPECT_EQ(0, pair.testIndex());
 520     EXPECT_TRUE(data.testPositions_[0].refPairs.count(pair.refIndex()) == 1);
 521
 522     pairSearch2.findNextPair(&pair);
 523     EXPECT_EQ(1, pair.testIndex());
 524     EXPECT_TRUE(data.testPositions_[1].refPairs.count(pair.refIndex()) == 1);
 525 }
 526
 527 TEST_F(NeighborhoodSearchTest, HandlesSkippingPairs)
 528 {
 529     const NeighborhoodSearchTestData &data = TrivialTestData::get();
 530
 531     nb_.setCutoff(data.cutoff_);
 532     gmx::AnalysisNeighborhoodSearch     search =
 533         nb_.initSearch(&data.pbc_, data.refPositions());
 534     gmx::AnalysisNeighborhoodPairSearch pairSearch =
 535         search.startPairSearch(data.testPositions());
 536     gmx::AnalysisNeighborhoodPair       pair;
 537     // TODO: This test needs to be adjusted if the grid search gets optimized
 538     // to loop over the test positions in cell order (first, the ordering
 539     // assumption here breaks, and second, it then needs to be tested
 540     // separately for simple and grid searches).
 541     int currentIndex = 0;
 542     while (pairSearch.findNextPair(&pair))
 543     {
 544         while (currentIndex < pair.testIndex())
 545         {
 546             ++currentIndex;
 547         }
 548         EXPECT_EQ(currentIndex, pair.testIndex());
 549         EXPECT_TRUE(data.testPositions_[currentIndex].refPairs.count(pair.refIndex()) == 1);
 550         pairSearch.skipRemainingPairsForTestPosition();
 551         ++currentIndex;
 552     }
 553 }
 554
 555 } // namespace