Support subset of rvec array as input for analysis nbsearch
authorTeemu Murtola <teemu.murtola@gmail.com>
Sat, 20 Dec 2014 10:29:25 +0000 (12:29 +0200)
committerMark Abraham <mark.j.abraham@gmail.com>
Mon, 5 Jan 2015 23:08:32 +0000 (00:08 +0100)
Make it possible to perform a neighbor search such that instead of using
a continuous rvec array as the reference/test positions, it is possible
considers a subset of such an array, specified by a separate index
array.

Also, add support for accessing the components of the shortest distance
vector between found pairs through AnalysisNeighborhoodPair.  This frees
the caller from a separate pbc_dx() call in case that distance vector is
of interest.

Tests could still be extended and improved, but left for a later change.

Change-Id: I253ddda483122367d9c41a5ed8bb172e476a90bb

src/gromacs/selection/nbsearch.cpp
src/gromacs/selection/nbsearch.h
src/gromacs/selection/tests/nbsearch.cpp
src/gromacs/utility/arrayref.h

index a8bd1addda128714a7da8ec9c851da8433200ed7..2fe48376b843f09b9f3678f77d58fb1b8a420abd 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2009,2010,2011,2012,2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2009,2010,2011,2012,2013,2014,2015, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -272,6 +272,8 @@ class AnalysisNeighborhoodSearchImpl
         const rvec             *xref_;
         //! Reference position exclusion IDs.
         const int              *refExclusionIds_;
+        //! Reference position indices (NULL if no indices).
+        const int              *refIndices_;
         //! Exclusions.
         const t_blocka         *excls_;
         //! PBC data.
@@ -325,7 +327,10 @@ class AnalysisNeighborhoodPairSearchImpl
         explicit AnalysisNeighborhoodPairSearchImpl(const AnalysisNeighborhoodSearchImpl &search)
             : search_(search)
         {
+            testPosCount_     = 0;
+            testPositions_    = NULL;
             testExclusionIds_ = NULL;
+            testIndices_      = NULL;
             nexcl_            = 0;
             excl_             = NULL;
             clear_rvec(xtest_);
@@ -353,10 +358,14 @@ class AnalysisNeighborhoodPairSearchImpl
 
         //! Parent search object.
         const AnalysisNeighborhoodSearchImpl   &search_;
+        //! Number of test positions.
+        int                                     testPosCount_;
         //! Reference to the test positions.
-        ConstArrayRef<rvec>                     testPositions_;
+        const rvec                             *testPositions_;
         //! Reference to the test exclusion indices.
         const int                              *testExclusionIds_;
+        //! Reference to the test position indices.
+        const int                              *testIndices_;
         //! Number of excluded reference positions for current test particle.
         int                                     nexcl_;
         //! Exclusions for current test particle.
@@ -369,6 +378,8 @@ class AnalysisNeighborhoodPairSearchImpl
         int                                     previ_;
         //! Stores the pair distance corresponding to previ_;
         real                                    prevr2_;
+        //! Stores the shortest distance vector corresponding to previ_;
+        rvec                                    prevdx_;
         //! Stores the current exclusion index during loops.
         int                                     exclind_;
         //! Stores the fractional test particle cell location during loops.
@@ -404,6 +415,7 @@ AnalysisNeighborhoodSearchImpl::AnalysisNeighborhoodSearchImpl(real cutoff)
     nref_            = 0;
     xref_            = NULL;
     refExclusionIds_ = NULL;
+    refIndices_      = NULL;
     std::memset(&pbc_, 0, sizeof(pbc_));
 
     bGrid_          = false;
@@ -916,6 +928,7 @@ void AnalysisNeighborhoodSearchImpl::init(
         bGrid_ = initGrid(pbc_, positions.count_, positions.x_, bUseBoundingBox,
                           mode == AnalysisNeighborhood::eSearchMode_Grid);
     }
+    refIndices_ = positions.indices_;
     if (bGrid_)
     {
         xrefAlloc_.resize(nref_);
@@ -923,11 +936,21 @@ void AnalysisNeighborhoodSearchImpl::init(
 
         for (int i = 0; i < nref_; ++i)
         {
-            rvec refcell;
-            mapPointToGridCell(positions.x_[i], refcell, xrefAlloc_[i]);
+            const int ii = (refIndices_ != NULL) ? refIndices_[i] : i;
+            rvec      refcell;
+            mapPointToGridCell(positions.x_[ii], refcell, xrefAlloc_[i]);
             addToGridCell(refcell, i);
         }
     }
+    else if (refIndices_ != NULL)
+    {
+        xrefAlloc_.resize(nref_);
+        xref_ = as_rvec_array(&xrefAlloc_[0]);
+        for (int i = 0; i < nref_; ++i)
+        {
+            copy_rvec(positions.x_[refIndices_[i]], xrefAlloc_[i]);
+        }
+    }
     else
     {
         xref_ = positions.x_;
@@ -951,22 +974,24 @@ void AnalysisNeighborhoodSearchImpl::init(
 void AnalysisNeighborhoodPairSearchImpl::reset(int testIndex)
 {
     testIndex_ = testIndex;
-    if (testIndex_ >= 0 && testIndex_ < static_cast<int>(testPositions_.size()))
+    if (testIndex_ >= 0 && testIndex_ < testPosCount_)
     {
+        const int index =
+            (testIndices_ != NULL ? testIndices_[testIndex] : testIndex);
         if (search_.bGrid_)
         {
-            search_.mapPointToGridCell(testPositions_[testIndex], testcell_, xtest_);
+            search_.mapPointToGridCell(testPositions_[index], testcell_, xtest_);
             search_.initCellRange(testcell_, currCell_, cellBound_, ZZ);
             search_.initCellRange(testcell_, currCell_, cellBound_, YY);
             search_.initCellRange(testcell_, currCell_, cellBound_, XX);
         }
         else
         {
-            copy_rvec(testPositions_[testIndex_], xtest_);
+            copy_rvec(testPositions_[index], xtest_);
         }
         if (search_.excls_ != NULL)
         {
-            const int exclIndex  = testExclusionIds_[testIndex];
+            const int exclIndex  = testExclusionIds_[index];
             if (exclIndex < search_.excls_->nr)
             {
                 const int startIndex = search_.excls_->index[exclIndex];
@@ -982,13 +1007,14 @@ void AnalysisNeighborhoodPairSearchImpl::reset(int testIndex)
     }
     previ_     = -1;
     prevr2_    = 0.0;
+    clear_rvec(prevdx_);
     exclind_   = 0;
     prevcai_   = -1;
 }
 
 void AnalysisNeighborhoodPairSearchImpl::nextTestPosition()
 {
-    if (testIndex_ < static_cast<int>(testPositions_.size()))
+    if (testIndex_ < testPosCount_)
     {
         ++testIndex_;
         reset(testIndex_);
@@ -999,7 +1025,9 @@ bool AnalysisNeighborhoodPairSearchImpl::isExcluded(int j)
 {
     if (exclind_ < nexcl_)
     {
-        const int refId = search_.refExclusionIds_[j];
+        const int index =
+            (search_.refIndices_ != NULL ? search_.refIndices_[j] : j);
+        const int refId = search_.refExclusionIds_[index];
         while (exclind_ < nexcl_ && excl_[exclind_] < refId)
         {
             ++exclind_;
@@ -1016,19 +1044,21 @@ bool AnalysisNeighborhoodPairSearchImpl::isExcluded(int j)
 void AnalysisNeighborhoodPairSearchImpl::startSearch(
         const AnalysisNeighborhoodPositions &positions)
 {
+    testPosCount_     = positions.count_;
+    testPositions_    = positions.x_;
     testExclusionIds_ = positions.exclusionIds_;
+    testIndices_      = positions.indices_;
     GMX_RELEASE_ASSERT(search_.excls_ == NULL || testExclusionIds_ != NULL,
                        "Exclusion IDs must be set when exclusions are enabled");
     if (positions.index_ < 0)
     {
-        testPositions_ = constArrayRefFromArray<rvec>(positions.x_, positions.count_);
         reset(0);
     }
     else
     {
         // Somewhat of a hack: setup the array such that only the last position
         // will be used.
-        testPositions_ = constArrayRefFromArray<rvec>(positions.x_, positions.index_ + 1);
+        testPosCount_ = positions.index_ + 1;
         reset(positions.index_);
     }
 }
@@ -1036,7 +1066,7 @@ void AnalysisNeighborhoodPairSearchImpl::startSearch(
 template <class Action>
 bool AnalysisNeighborhoodPairSearchImpl::searchNext(Action action)
 {
-    while (testIndex_ < static_cast<int>(testPositions_.size()))
+    while (testIndex_ < testPosCount_)
     {
         if (search_.bGrid_)
         {
@@ -1055,19 +1085,20 @@ bool AnalysisNeighborhoodPairSearchImpl::searchNext(Action action)
                         continue;
                     }
                     rvec       dx;
-                    rvec_sub(xtest_, search_.xref_[i], dx);
-                    rvec_add(dx, shift, dx);
+                    rvec_sub(search_.xref_[i], xtest_, dx);
+                    rvec_sub(dx, shift, dx);
                     const real r2
                         = search_.bXY_
                             ? dx[XX]*dx[XX] + dx[YY]*dx[YY]
                             : norm2(dx);
                     if (r2 <= search_.cutoff2_)
                     {
-                        if (action(i, r2))
+                        if (action(i, r2, dx))
                         {
                             prevcai_ = cai;
                             previ_   = i;
                             prevr2_  = r2;
+                            copy_rvec(dx, prevdx_);
                             return true;
                         }
                     }
@@ -1088,11 +1119,11 @@ bool AnalysisNeighborhoodPairSearchImpl::searchNext(Action action)
                 rvec dx;
                 if (search_.pbc_.ePBC != epbcNONE)
                 {
-                    pbc_dx(&search_.pbc_, xtest_, search_.xref_[i], dx);
+                    pbc_dx(&search_.pbc_, search_.xref_[i], xtest_, dx);
                 }
                 else
                 {
-                    rvec_sub(xtest_, search_.xref_[i], dx);
+                    rvec_sub(search_.xref_[i], xtest_, dx);
                 }
                 const real r2
                     = search_.bXY_
@@ -1100,10 +1131,11 @@ bool AnalysisNeighborhoodPairSearchImpl::searchNext(Action action)
                         : norm2(dx);
                 if (r2 <= search_.cutoff2_)
                 {
-                    if (action(i, r2))
+                    if (action(i, r2, dx))
                     {
                         previ_  = i;
                         prevr2_ = r2;
+                        copy_rvec(dx, prevdx_);
                         return true;
                     }
                 }
@@ -1123,7 +1155,7 @@ void AnalysisNeighborhoodPairSearchImpl::initFoundPair(
     }
     else
     {
-        *pair = AnalysisNeighborhoodPair(previ_, testIndex_, prevr2_);
+        *pair = AnalysisNeighborhoodPair(previ_, testIndex_, prevr2_, prevdx_);
     }
 }
 
@@ -1140,7 +1172,7 @@ namespace
  *
  * Simply breaks the loop on the first found neighbor.
  */
-bool withinAction(int /*i*/, real /*r2*/)
+bool withinAction(int /*i*/, real /*r2*/, const rvec /* dx */)
 {
     return true;
 }
@@ -1164,23 +1196,25 @@ class MindistAction
          *
          * \param[out] closestPoint Index of the closest reference location.
          * \param[out] minDist2     Minimum distance squared.
+         * \param[out] dx           Shortest distance vector.
          *
          * The constructor call does not modify the pointed values, but only
          * stores the pointers for later use.
          * See the class description for additional semantics.
          */
-        MindistAction(int *closestPoint, real *minDist2)
-            : closestPoint_(*closestPoint), minDist2_(*minDist2)
+        MindistAction(int *closestPoint, real *minDist2, rvec *dx)
+            : closestPoint_(*closestPoint), minDist2_(*minDist2), dx_(*dx)
         {
         }
 
         //! Processes a neighbor to find the nearest point.
-        bool operator()(int i, real r2)
+        bool operator()(int i, real r2, const rvec dx)
         {
             if (r2 < minDist2_)
             {
                 closestPoint_ = i;
                 minDist2_     = r2;
+                copy_rvec(dx, dx_);
             }
             return false;
         }
@@ -1188,6 +1222,7 @@ class MindistAction
     private:
         int     &closestPoint_;
         real    &minDist2_;
+        rvec    &dx_;
 
         GMX_DISALLOW_ASSIGN(MindistAction);
 };
@@ -1349,7 +1384,8 @@ real AnalysisNeighborhoodSearch::minimumDistance(
     pairSearch.startSearch(positions);
     real          minDist2     = impl_->cutoffSquared();
     int           closestPoint = -1;
-    MindistAction action(&closestPoint, &minDist2);
+    rvec          dx           = {0.0, 0.0, 0.0};
+    MindistAction action(&closestPoint, &minDist2, &dx);
     (void)pairSearch.searchNext(action);
     return sqrt(minDist2);
 }
@@ -1363,9 +1399,10 @@ AnalysisNeighborhoodSearch::nearestPoint(
     pairSearch.startSearch(positions);
     real          minDist2     = impl_->cutoffSquared();
     int           closestPoint = -1;
-    MindistAction action(&closestPoint, &minDist2);
+    rvec          dx           = {0.0, 0.0, 0.0};
+    MindistAction action(&closestPoint, &minDist2, &dx);
     (void)pairSearch.searchNext(action);
-    return AnalysisNeighborhoodPair(closestPoint, 0, minDist2);
+    return AnalysisNeighborhoodPair(closestPoint, 0, minDist2, dx);
 }
 
 AnalysisNeighborhoodPairSearch
index 4419c2b27d07c6a151e52ed521b1d9b99835a358..bc5baba8d156311576a793b962ea429b0a2b30d8 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2009,2010,2011,2012,2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2009,2010,2011,2012,2013,2014,2015, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -53,6 +53,7 @@
 
 #include <boost/shared_ptr.hpp>
 
+#include "gromacs/math/vec.h"
 #include "gromacs/math/vectypes.h"
 #include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/classhelpers.h"
@@ -110,14 +111,14 @@ class AnalysisNeighborhoodPositions
          * to methods that accept positions.
          */
         AnalysisNeighborhoodPositions(const rvec &x)
-            : count_(1), index_(-1), x_(&x), exclusionIds_(NULL)
+            : count_(1), index_(-1), x_(&x), exclusionIds_(NULL), indices_(NULL)
         {
         }
         /*! \brief
          * Initializes positions from an array of position vectors.
          */
         AnalysisNeighborhoodPositions(const rvec x[], int count)
-            : count_(count), index_(-1), x_(x), exclusionIds_(NULL)
+            : count_(count), index_(-1), x_(x), exclusionIds_(NULL), indices_(NULL)
         {
         }
         /*! \brief
@@ -125,7 +126,7 @@ class AnalysisNeighborhoodPositions
          */
         AnalysisNeighborhoodPositions(const std::vector<RVec> &x)
             : count_(x.size()), index_(-1), x_(as_rvec_array(&x[0])),
-              exclusionIds_(NULL)
+              exclusionIds_(NULL), indices_(NULL)
         {
         }
 
@@ -144,6 +145,21 @@ class AnalysisNeighborhoodPositions
             exclusionIds_ = ids.data();
             return *this;
         }
+        /*! \brief
+         * Sets indices that select a subset of all positions from the array.
+         *
+         * If called, selected positions from the array of positions passed to
+         * the constructor is used instead of the whole array.
+         * All returned indices from AnalysisNeighborhoodPair objects are
+         * indices to the \p indices array passed here.
+         */
+        AnalysisNeighborhoodPositions &
+        indexed(ConstArrayRef<int> indices)
+        {
+            count_   = indices.size();
+            indices_ = indices.data();
+            return *this;
+        }
 
         /*! \brief
          * Selects a single position to use from an array.
@@ -153,6 +169,9 @@ class AnalysisNeighborhoodPositions
          * In contrast to the AnalysisNeighborhoodPositions(const rvec &)
          * constructor, AnalysisNeighborhoodPair objects return \p index
          * instead of zero.
+         *
+         * If used together with indexed(), \p index references the index array
+         * passed to indexed() instead of the position array.
          */
         AnalysisNeighborhoodPositions &selectSingleFromArray(int index)
         {
@@ -166,6 +185,7 @@ class AnalysisNeighborhoodPositions
         int                     index_;
         const rvec             *x_;
         const int              *exclusionIds_;
+        const int              *indices_;
 
         //! To access the positions for initialization.
         friend class internal::AnalysisNeighborhoodSearchImpl;
@@ -332,11 +352,16 @@ class AnalysisNeighborhoodPair
 {
     public:
         //! Initializes an invalid pair.
-        AnalysisNeighborhoodPair() : refIndex_(-1), testIndex_(0), distance2_(0.0) {}
+        AnalysisNeighborhoodPair() : refIndex_(-1), testIndex_(0), distance2_(0.0)
+        {
+            clear_rvec(dx_);
+        }
         //! Initializes a pair object with the given data.
-        AnalysisNeighborhoodPair(int refIndex, int testIndex, real distance2)
+        AnalysisNeighborhoodPair(int refIndex, int testIndex, real distance2,
+                                 const rvec dx)
             : refIndex_(refIndex), testIndex_(testIndex), distance2_(distance2)
         {
+            copy_rvec(dx, dx_);
         }
 
         /*! \brief
@@ -377,11 +402,22 @@ class AnalysisNeighborhoodPair
             GMX_ASSERT(isValid(), "Accessing invalid object");
             return distance2_;
         }
+        /*! \brief
+         * Returns the shortest vector between the pair of positions.
+         *
+         * The vector is from the test position to the reference position.
+         */
+        const rvec &dx() const
+        {
+            GMX_ASSERT(isValid(), "Accessing invalid object");
+            return dx_;
+        }
 
     private:
         int                     refIndex_;
         int                     testIndex_;
         real                    distance2_;
+        rvec                    dx_;
 };
 
 /*! \brief
index 5d876aac3e4018c46a7211f50c5731649bfbab75..471dae9ab8a3c758ad16f70b1cdd22f8ba3da15b 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -80,7 +80,7 @@ class NeighborhoodSearchTestData
         {
             RefPair(int refIndex, real distance)
                 : refIndex(refIndex), distance(distance), bFound(false),
-                  bExcluded(false)
+                  bExcluded(false), bIndexed(true)
             {
             }
 
@@ -97,6 +97,7 @@ class NeighborhoodSearchTestData
             // Simpler to have just a single structure for both purposes.
             bool                bFound;
             bool                bExcluded;
+            bool                bIndexed;
         };
 
         struct TestPosition
@@ -145,7 +146,8 @@ class NeighborhoodSearchTestData
                                "Cannot add positions after testPositions() call");
             testPositions_.push_back(TestPosition(x));
         }
-        gmx::RVec generateRandomPosition() const;
+        gmx::RVec generateRandomPosition();
+        std::vector<int> generateIndex(int count) const;
         void generateRandomRefPositions(int count);
         void generateRandomTestPositions(int count);
         void computeReferences(t_pbc *pbc)
@@ -203,7 +205,7 @@ NeighborhoodSearchTestData::~NeighborhoodSearchTestData()
     }
 }
 
-gmx::RVec NeighborhoodSearchTestData::generateRandomPosition() const
+gmx::RVec NeighborhoodSearchTestData::generateRandomPosition()
 {
     rvec fx, x;
     fx[XX] = gmx_rng_uniform_real(rng_);
@@ -217,6 +219,19 @@ gmx::RVec NeighborhoodSearchTestData::generateRandomPosition() const
     return x;
 }
 
+std::vector<int> NeighborhoodSearchTestData::generateIndex(int count) const
+{
+    std::vector<int> result;
+    for (int i = 0; i < count; ++i)
+    {
+        if (gmx_rng_uniform_real(rng_) > 0.5)
+        {
+            result.push_back(i);
+        }
+    }
+    return result;
+}
+
 void NeighborhoodSearchTestData::generateRandomRefPositions(int count)
 {
     refPosCount_ = count;
@@ -395,10 +410,14 @@ class NeighborhoodSearchTest : public ::testing::Test
                               const NeighborhoodSearchTestData &data);
         void testPairSearch(gmx::AnalysisNeighborhoodSearch  *search,
                             const NeighborhoodSearchTestData &data);
+        void testPairSearchIndexed(gmx::AnalysisNeighborhood        *nb,
+                                   const NeighborhoodSearchTestData &data);
         void testPairSearchFull(gmx::AnalysisNeighborhoodSearch          *search,
                                 const NeighborhoodSearchTestData         &data,
                                 const gmx::AnalysisNeighborhoodPositions &pos,
-                                const t_blocka                           *excls);
+                                const t_blocka                           *excls,
+                                const gmx::ConstArrayRef<int>            &refIndices,
+                                const gmx::ConstArrayRef<int>            &testIndices);
 
         gmx::AnalysisNeighborhood        nb_;
 };
@@ -492,39 +511,69 @@ void NeighborhoodSearchTest::testPairSearch(
         gmx::AnalysisNeighborhoodSearch  *search,
         const NeighborhoodSearchTestData &data)
 {
-    testPairSearchFull(search, data, data.testPositions(), NULL);
+    testPairSearchFull(search, data, data.testPositions(), NULL,
+                       gmx::EmptyArrayRef(), gmx::EmptyArrayRef());
+}
+
+void NeighborhoodSearchTest::testPairSearchIndexed(
+        gmx::AnalysisNeighborhood        *nb,
+        const NeighborhoodSearchTestData &data)
+{
+    std::vector<int>                refIndices(data.generateIndex(data.refPos_.size()));
+    std::vector<int>                testIndices(data.generateIndex(data.testPositions_.size()));
+    gmx::AnalysisNeighborhoodSearch search =
+        nb->initSearch(&data.pbc_,
+                       data.refPositions().indexed(refIndices));
+    testPairSearchFull(&search, data, data.testPositions(), NULL,
+                       refIndices, testIndices);
 }
 
 void NeighborhoodSearchTest::testPairSearchFull(
         gmx::AnalysisNeighborhoodSearch          *search,
         const NeighborhoodSearchTestData         &data,
         const gmx::AnalysisNeighborhoodPositions &pos,
-        const t_blocka                           *excls)
+        const t_blocka                           *excls,
+        const gmx::ConstArrayRef<int>            &refIndices,
+        const gmx::ConstArrayRef<int>            &testIndices)
 {
     // TODO: Some parts of this code do not work properly if pos does not
-    // contain all the test positions.
+    // initially contain all the test positions.
     std::set<int> remainingTestPositions;
-    for (size_t i = 0; i < data.testPositions_.size(); ++i)
+    gmx::AnalysisNeighborhoodPositions  posCopy(pos);
+    if (testIndices.empty())
     {
-        remainingTestPositions.insert(i);
+        for (size_t i = 0; i < data.testPositions_.size(); ++i)
+        {
+            remainingTestPositions.insert(i);
+        }
+    }
+    else
+    {
+        remainingTestPositions.insert(testIndices.begin(), testIndices.end());
+        posCopy.indexed(testIndices);
     }
     gmx::AnalysisNeighborhoodPairSearch pairSearch
-        = search->startPairSearch(pos);
+        = search->startPairSearch(posCopy);
     gmx::AnalysisNeighborhoodPair       pair;
-    // TODO: There is an ordering assumption here that may break in the future:
-    // all pairs for a test position are assumed to be returned consencutively.
+    // There is an ordering assumption here that all pairs for a test position
+    // are returned consencutively; with the current optimizations in the
+    // search code, this is reasoable, as the set of grid cell pairs searched
+    // depends on the test position.
     RefPairList refPairs;
     int         prevTestPos = -1;
     while (pairSearch.findNextPair(&pair))
     {
-        if (pair.testIndex() != prevTestPos)
+        const int testIndex =
+            (testIndices.empty() ? pair.testIndex() : testIndices[pair.testIndex()]);
+        const int refIndex =
+            (refIndices.empty() ? pair.refIndex() : refIndices[pair.refIndex()]);
+        if (testIndex != prevTestPos)
         {
             if (prevTestPos != -1)
             {
                 checkAllPairsFound(refPairs, data.refPos_, prevTestPos,
                                    data.testPositions_[prevTestPos].x);
             }
-            const int testIndex = pair.testIndex();
             if (remainingTestPositions.count(testIndex) == 0)
             {
                 ADD_FAILURE()
@@ -537,32 +586,63 @@ void NeighborhoodSearchTest::testPairSearchFull(
             {
                 ExclusionsHelper::markExcludedPairs(&refPairs, testIndex, excls);
             }
+            if (!refIndices.empty())
+            {
+                RefPairList::iterator refPair;
+                for (refPair = refPairs.begin(); refPair != refPairs.end(); ++refPair)
+                {
+                    refPair->bIndexed = false;
+                }
+                for (size_t i = 0; i < refIndices.size(); ++i)
+                {
+                    NeighborhoodSearchTestData::RefPair searchPair(refIndices[i], 0.0);
+                    refPair = std::lower_bound(refPairs.begin(), refPairs.end(), searchPair);
+                    if (refPair != refPairs.end() && refPair->refIndex == refIndices[i])
+                    {
+                        refPair->bIndexed = true;
+                    }
+                }
+                for (refPair = refPairs.begin(); refPair != refPairs.end(); ++refPair)
+                {
+                    if (!refPair->bIndexed)
+                    {
+                        refPair->bFound = true;
+                    }
+                }
+            }
             prevTestPos = testIndex;
         }
 
-        NeighborhoodSearchTestData::RefPair searchPair(pair.refIndex(),
+        NeighborhoodSearchTestData::RefPair searchPair(refIndex,
                                                        sqrt(pair.distance2()));
         RefPairList::iterator               foundRefPair
             = std::lower_bound(refPairs.begin(), refPairs.end(), searchPair);
-        if (foundRefPair == refPairs.end() || foundRefPair->refIndex != pair.refIndex())
+        if (foundRefPair == refPairs.end() || foundRefPair->refIndex != refIndex)
         {
             ADD_FAILURE()
-            << "Expected: Pair (ref: " << pair.refIndex() << ", test: "
-            << pair.testIndex() << ") is not within the cutoff.\n"
+            << "Expected: Pair (ref: " << refIndex << ", test: " << testIndex
+            << ") is not within the cutoff.\n"
             << "  Actual: It is returned.";
         }
         else if (foundRefPair->bExcluded)
         {
             ADD_FAILURE()
-            << "Expected: Pair (ref: " << pair.refIndex() << ", test: "
-            << pair.testIndex() << ") is excluded from the search.\n"
+            << "Expected: Pair (ref: " << refIndex << ", test: " << testIndex
+            << ") is excluded from the search.\n"
+            << "  Actual: It is returned.";
+        }
+        else if (!foundRefPair->bIndexed)
+        {
+            ADD_FAILURE()
+            << "Expected: Pair (ref: " << refIndex << ", test: " << testIndex
+            << ") is not part of the indexed set.\n"
             << "  Actual: It is returned.";
         }
         else if (foundRefPair->bFound)
         {
             ADD_FAILURE()
-            << "Expected: Pair (ref: " << pair.refIndex() << ", test: "
-            << pair.testIndex() << ") is returned only once.\n"
+            << "Expected: Pair (ref: " << refIndex << ", test: " << testIndex
+            << ") is returned only once.\n"
             << "  Actual: It is returned multiple times.";
             return;
         }
@@ -788,6 +868,9 @@ TEST_F(NeighborhoodSearchTest, SimpleSearch)
     testMinimumDistance(&search, data);
     testNearestPoint(&search, data);
     testPairSearch(&search, data);
+
+    search.reset();
+    testPairSearchIndexed(&nb_, data);
 }
 
 TEST_F(NeighborhoodSearchTest, SimpleSearchXY)
@@ -821,6 +904,9 @@ TEST_F(NeighborhoodSearchTest, GridSearchBox)
     testMinimumDistance(&search, data);
     testNearestPoint(&search, data);
     testPairSearch(&search, data);
+
+    search.reset();
+    testPairSearchIndexed(&nb_, data);
 }
 
 TEST_F(NeighborhoodSearchTest, GridSearchTriclinic)
@@ -993,7 +1079,7 @@ TEST_F(NeighborhoodSearchTest, SimpleSearchExclusions)
 
     testPairSearchFull(&search, data,
                        data.testPositions().exclusionIds(helper.testPosIds()),
-                       helper.exclusions());
+                       helper.exclusions(), gmx::EmptyArrayRef(), gmx::EmptyArrayRef());
 }
 
 TEST_F(NeighborhoodSearchTest, GridSearchExclusions)
@@ -1013,7 +1099,7 @@ TEST_F(NeighborhoodSearchTest, GridSearchExclusions)
 
     testPairSearchFull(&search, data,
                        data.testPositions().exclusionIds(helper.testPosIds()),
-                       helper.exclusions());
+                       helper.exclusions(), gmx::EmptyArrayRef(), gmx::EmptyArrayRef());
 }
 
 } // namespace
index fb4778d24a7573b6cc6d3cad83538c712c4bf7d1..d59a3548577305800b7836d5f6a7d1acc1ca58f5 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2012,2013,2014,2015, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -153,6 +153,22 @@ class ArrayRef
         {
             GMX_ASSERT(end >= begin, "Invalid range");
         }
+        /*! \brief
+         * Constructs a reference to a whole vector.
+         *
+         * \param[in] v  Vector to reference.
+         *
+         * Passed vector must remain valid and not be reallocated for the
+         * lifetime of this object.
+         *
+         * This constructor is not explicit to allow directly passing
+         * std::vector to a method that takes ArrayRef.
+         */
+        ArrayRef(std::vector<T> &v)
+            : begin_((!v.empty()) ? &v[0] : NULL),
+              end_((!v.empty()) ? &v[0] + v.size() : NULL)
+        {
+        }
         //! \cond
         // Doxygen 1.8.5 doesn't parse the declaration correctly...
         /*! \brief
@@ -393,6 +409,22 @@ class ConstArrayRef
         {
             GMX_ASSERT(end >= begin, "Invalid range");
         }
+        /*! \brief
+         * Constructs a reference to a whole vector.
+         *
+         * \param[in] v  Vector to reference.
+         *
+         * Passed vector must remain valid and not be reallocated for the
+         * lifetime of this object.
+         *
+         * This constructor is not explicit to allow directly passing
+         * std::vector to a method that takes ConstArrayRef.
+         */
+        ConstArrayRef(const std::vector<T> &v)
+            : begin_((!v.empty()) ? &v[0] : NULL),
+              end_((!v.empty()) ? &v[0] + v.size() : NULL)
+        {
+        }
         //! \cond
         // Doxygen 1.8.5 doesn't parse the declaration correctly...
         /*! \brief