src/gromacs/selection/nbsearch.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2017, 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 /*! \file
  36  * \brief API for neighborhood searching for analysis.
  37  *
  38  * The main part of the API is the class gmx::AnalysisNeighborhood.
  39  * See \ref page_analysisnbsearch for an overview.
  40  *
  41  * The classes within this file can be used independently of the other parts
  42  * of the selection module.
  43  *
  44  * \author Teemu Murtola <teemu.murtola@gmail.com>
  45  * \inpublicapi
  46  * \ingroup module_selection
  47  */
  48 #ifndef GMX_SELECTION_NBSEARCH_H
  49 #define GMX_SELECTION_NBSEARCH_H
  50
  51 #include <memory>
  52 #include <vector>
  53
  54 #include "gromacs/math/vec.h"
  55 #include "gromacs/math/vectypes.h"
  56 #include "gromacs/utility/arrayref.h"
  57 #include "gromacs/utility/classhelpers.h"
  58 #include "gromacs/utility/gmxassert.h"
  59 #include "gromacs/utility/real.h"
  60
  61 struct t_blocka;
  62 struct t_pbc;
  63
  64 namespace gmx
  65 {
  66
  67 namespace internal
  68 {
  69 class AnalysisNeighborhoodSearchImpl;
  70 class AnalysisNeighborhoodPairSearchImpl;
  71 };
  72
  73 class AnalysisNeighborhoodSearch;
  74 class AnalysisNeighborhoodPairSearch;
  75
  76 /*! \brief
  77  * Input positions for neighborhood searching.
  78  *
  79  * This class supports uniformly specifying sets of positions for various
  80  * methods in the analysis neighborhood searching classes
  81  * (AnalysisNeighborhood and AnalysisNeighborhoodSearch).
  82  *
  83  * Note that copies are not made: only a reference to the positions passed to
  84  * the constructors are kept.  The caller is responsible to ensure that those
  85  * positions remain in scope as long as the neighborhood search object requires
  86  * access to them.
  87  *
  88  * Also note that in addition to constructors here, Selection and
  89  * SelectionPosition provide conversions operators to this type.  It is done
  90  * this way to not introduce a cyclic dependency between the selection code and
  91  * the neighborhood search code, which in turn allows splitting this search
  92  * code into a separate lower-level module if desired at some point.
  93  *
  94  * Methods in this class do not throw.
  95  *
  96  * \inpublicapi
  97  * \ingroup module_selection
  98  */
  99 class AnalysisNeighborhoodPositions
 100 {
 101     public:
 102         /*! \brief
 103          * Initializes positions from a single position vector.
 104          *
 105          * For positions initialized this way, AnalysisNeighborhoodPair always
 106          * returns zero in the corresponding index.
 107          *
 108          * This constructor is not explicit to allow directly passing an rvec
 109          * to methods that accept positions.
 110          */
 111         AnalysisNeighborhoodPositions(const rvec &x)
 112             : count_(1), index_(-1), x_(&x), exclusionIds_(nullptr), indices_(nullptr)
 113         {
 114         }
 115         /*! \brief
 116          * Initializes positions from an array of position vectors.
 117          */
 118         AnalysisNeighborhoodPositions(const rvec x[], int count)
 119             : count_(count), index_(-1), x_(x), exclusionIds_(nullptr), indices_(nullptr)
 120         {
 121         }
 122         /*! \brief
 123          * Initializes positions from a vector of position vectors.
 124          */
 125         AnalysisNeighborhoodPositions(const std::vector<RVec> &x)
 126             : count_(x.size()), index_(-1), x_(as_rvec_array(x.data())),
 127               exclusionIds_(nullptr), indices_(nullptr)
 128         {
 129         }
 130
 131         /*! \brief
 132          * Sets indices to use for mapping exclusions to these positions.
 133          *
 134          * The exclusion IDs can always be set, but they are ignored unless
 135          * actual exclusions have been set with
 136          * AnalysisNeighborhood::setTopologyExclusions().
 137          */
 138         AnalysisNeighborhoodPositions &
 139         exclusionIds(ConstArrayRef<int> ids)
 140         {
 141             GMX_ASSERT(static_cast<int>(ids.size()) == count_,
 142                        "Exclusion id array should match the number of positions");
 143             exclusionIds_ = ids.data();
 144             return *this;
 145         }
 146         /*! \brief
 147          * Sets indices that select a subset of all positions from the array.
 148          *
 149          * If called, selected positions from the array of positions passed to
 150          * the constructor is used instead of the whole array.
 151          * All returned indices from AnalysisNeighborhoodPair objects are
 152          * indices to the \p indices array passed here.
 153          */
 154         AnalysisNeighborhoodPositions &
 155         indexed(ConstArrayRef<int> indices)
 156         {
 157             count_   = indices.size();
 158             indices_ = indices.data();
 159             return *this;
 160         }
 161
 162         /*! \brief
 163          * Selects a single position to use from an array.
 164          *
 165          * If called, a single position from the array of positions passed to
 166          * the constructor is used instead of the whole array.
 167          * In contrast to the AnalysisNeighborhoodPositions(const rvec &)
 168          * constructor, AnalysisNeighborhoodPair objects return \p index
 169          * instead of zero.
 170          *
 171          * If used together with indexed(), \p index references the index array
 172          * passed to indexed() instead of the position array.
 173          */
 174         AnalysisNeighborhoodPositions &selectSingleFromArray(int index)
 175         {
 176             GMX_ASSERT(index >= 0 && index < count_, "Invalid position index");
 177             index_ = index;
 178             return *this;
 179         }
 180
 181     private:
 182         int                     count_;
 183         int                     index_;
 184         const rvec             *x_;
 185         const int              *exclusionIds_;
 186         const int              *indices_;
 187
 188         //! To access the positions for initialization.
 189         friend class internal::AnalysisNeighborhoodSearchImpl;
 190         //! To access the positions for initialization.
 191         friend class internal::AnalysisNeighborhoodPairSearchImpl;
 192 };
 193
 194 /*! \brief
 195  * Neighborhood searching for analysis tools.
 196  *
 197  * See \ref page_analysisnbsearch for an overview.
 198  *
 199  * To use the search, create an object of this type, call setCutoff() to
 200  * initialize it, and then repeatedly call initSearch() to start a search with
 201  * different sets of reference positions.  For each set of reference positions,
 202  * use methods in the returned AnalysisNeighborhoodSearch to find the reference
 203  * positions that are within the given cutoff from a provided position.
 204  *
 205  * initSearch() is thread-safe and can be called from multiple threads.  Each
 206  * call returns a different instance of the search object that can be used
 207  * independently of the others.  The returned AnalysisNeighborhoodSearch
 208  * objects are also thread-safe, and can be used concurrently from multiple
 209  * threads.  It is also possible to create multiple concurrent searches within
 210  * a single thread.
 211  *
 212  * \todo
 213  * Generalize the exclusion machinery to make it easier to use for other cases
 214  * than atom-atom exclusions from the topology.
 215  *
 216  * \inpublicapi
 217  * \ingroup module_selection
 218  */
 219 class AnalysisNeighborhood
 220 {
 221     public:
 222         //! Searching algorithm to use.
 223         enum SearchMode
 224         {
 225             //! Select algorithm based on heuristic efficiency considerations.
 226             eSearchMode_Automatic,
 227             //! Use a simple loop over all pairs.
 228             eSearchMode_Simple,
 229             //! Use grid-based searching whenever possible.
 230             eSearchMode_Grid
 231         };
 232
 233         //! Creates an uninitialized neighborhood search.
 234         AnalysisNeighborhood();
 235         ~AnalysisNeighborhood();
 236
 237         /*! \brief
 238          * Sets cutoff distance for the neighborhood searching.
 239          *
 240          * \param[in]  cutoff Cutoff distance for the search
 241          *   (<=0 stands for no cutoff).
 242          *
 243          * Currently, can only be called before the first call to initSearch().
 244          * If this method is not called, no cutoff is used in the searches.
 245          *
 246          * Does not throw.
 247          */
 248         void setCutoff(real cutoff);
 249         /*! \brief
 250          * Sets the search to only happen in the XY plane.
 251          *
 252          * Z component of the coordinates is not used in the searching,
 253          * and returned distances are computed in the XY plane.
 254          * Only boxes with the third box vector parallel to the Z axis are
 255          * currently implemented.
 256          *
 257          * Does not throw.
 258          */
 259         void setXYMode(bool bXY);
 260         /*! \brief
 261          * Sets atom exclusions from a topology.
 262          *
 263          * The \p excls structure specifies the exclusions from test positions
 264          * to reference positions, i.e., a block starting at `excls->index[i]`
 265          * specifies the exclusions for test position `i`, and the indices in
 266          * `excls->a` are indices of the reference positions.  If `excls->nr`
 267          * is smaller than a test position id, then such test positions do not
 268          * have any exclusions.
 269          * It is assumed that the indices within a block of indices in
 270          * `excls->a` is ascending.
 271          *
 272          * Does not throw.
 273          *
 274          * \see AnalysisNeighborhoodPositions::exclusionIds()
 275          */
 276         void setTopologyExclusions(const t_blocka *excls);
 277         /*! \brief
 278          * Sets the algorithm to use for searching.
 279          *
 280          * \param[in] mode  Search mode to use.
 281          *
 282          * Note that if \p mode is \ref eSearchMode_Grid, it is still only a
 283          * suggestion: grid-based searching may not be possible with the
 284          * provided input, in which case a simple search is still used.
 285          * This is mainly useful for testing purposes to force a mode.
 286          *
 287          * Does not throw.
 288          */
 289         void setMode(SearchMode mode);
 290         //! Returns the currently active search mode.
 291         SearchMode mode() const;
 292
 293         /*! \brief
 294          * Initializes neighborhood search for a set of positions.
 295          *
 296          * \param[in] pbc        PBC information for the frame.
 297          * \param[in] positions  Set of reference positions to use.
 298          * \returns   Search object that can be used to find positions from
 299          *      \p x within the given cutoff.
 300          * \throws    std::bad_alloc if out of memory.
 301          *
 302          * Currently, the input positions cannot use
 303          * AnalysisNeighborhoodPositions::selectSingleFromArray().
 304          */
 305         AnalysisNeighborhoodSearch
 306         initSearch(const t_pbc                         *pbc,
 307                    const AnalysisNeighborhoodPositions &positions);
 308
 309     private:
 310         class Impl;
 311
 312         PrivateImplPointer<Impl> impl_;
 313 };
 314
 315 /*! \brief
 316  * Value type to represent a pair of positions found in neighborhood searching.
 317  *
 318  * Methods in this class do not throw.
 319  *
 320  * \inpublicapi
 321  * \ingroup module_selection
 322  */
 323 class AnalysisNeighborhoodPair
 324 {
 325     public:
 326         //! Initializes an invalid pair.
 327         AnalysisNeighborhoodPair() : refIndex_(-1), testIndex_(0), distance2_(0.0)
 328         {
 329             clear_rvec(dx_);
 330         }
 331         //! Initializes a pair object with the given data.
 332         AnalysisNeighborhoodPair(int refIndex, int testIndex, real distance2,
 333                                  const rvec dx)
 334             : refIndex_(refIndex), testIndex_(testIndex), distance2_(distance2)
 335         {
 336             copy_rvec(dx, dx_);
 337         }
 338
 339         /*! \brief
 340          * Whether this pair is valid.
 341          *
 342          * If isValid() returns false, other methods should not be called.
 343          */
 344         bool isValid() const { return refIndex_ >= 0; }
 345
 346         /*! \brief
 347          * Returns the index of the reference position in the pair.
 348          *
 349          * This index is always the index into the position array provided to
 350          * AnalysisNeighborhood::initSearch().
 351          */
 352         int refIndex() const
 353         {
 354             GMX_ASSERT(isValid(), "Accessing invalid object");
 355             return refIndex_;
 356         }
 357         /*! \brief
 358          * Returns the index of the test position in the pair.
 359          *
 360          * The contents of this index depends on the context (method call) that
 361          * produces the pair.
 362          * If there was no array in the call, this index is zero.
 363          */
 364         int testIndex() const
 365         {
 366             GMX_ASSERT(isValid(), "Accessing invalid object");
 367             return testIndex_;
 368         }
 369         /*! \brief
 370          * Returns the squared distance between the pair of positions.
 371          */
 372         real distance2() const
 373         {
 374             GMX_ASSERT(isValid(), "Accessing invalid object");
 375             return distance2_;
 376         }
 377         /*! \brief
 378          * Returns the shortest vector between the pair of positions.
 379          *
 380          * The vector is from the test position to the reference position.
 381          */
 382         const rvec &dx() const
 383         {
 384             GMX_ASSERT(isValid(), "Accessing invalid object");
 385             return dx_;
 386         }
 387
 388     private:
 389         int                     refIndex_;
 390         int                     testIndex_;
 391         real                    distance2_;
 392         rvec                    dx_;
 393 };
 394
 395 /*! \brief
 396  * Initialized neighborhood search with a fixed set of reference positions.
 397  *
 398  * An instance of this class is obtained through
 399  * AnalysisNeighborhood::initSearch(), and can be used to do multiple searches
 400  * against the provided set of reference positions.
 401  * It is possible to create concurrent pair searches (including from different
 402  * threads), as well as call other methods in this class while a pair search is
 403  * in progress.
 404  *
 405  * This class works like a pointer: copies of it point to the same search.
 406  * In general, avoid creating copies, and only use the copy/assignment support
 407  * for moving the variable around.  With C++11, this class would best be
 408  * movable.
 409  *
 410  * Methods in this class do not throw unless otherwise indicated.
 411  *
 412  * \todo
 413  * Make it such that reset() is not necessary to call in code that repeatedly
 414  * assigns the result of AnalysisNeighborhood::initSearch() to the same
 415  * variable (see sm_distance.cpp).
 416  *
 417  * \todo
 418  * Consider merging nearestPoint() and minimumDistance() by adding the distance
 419  * to AnalysisNeighborhoodPair.
 420  *
 421  * \inpublicapi
 422  * \ingroup module_selection
 423  */
 424 class AnalysisNeighborhoodSearch
 425 {
 426     public:
 427         /*! \brief
 428          * Internal short-hand type for a pointer to the implementation class.
 429          *
 430          * shared_ptr is used here to automatically keep a reference count to
 431          * track whether an implementation class is still used outside the
 432          * AnalysisNeighborhood object.  Ownership currently always stays with
 433          * AnalysisNeighborhood; it always keeps one instance of the pointer.
 434          */
 435         typedef std::shared_ptr<internal::AnalysisNeighborhoodSearchImpl>
 436             ImplPointer;
 437
 438         /*! \brief
 439          * Initializes an invalid search.
 440          *
 441          * Such an object cannot be used for searching.  It needs to be
 442          * assigned a value from AnalysisNeighborhood::initSearch() before it
 443          * can be used.  Provided to allow declaring a variable to hold the
 444          * search before calling AnalysisNeighborhood::initSearch().
 445          */
 446         AnalysisNeighborhoodSearch();
 447         /*! \brief
 448          * Internally initialize the search.
 449          *
 450          * Used to implement AnalysisNeighborhood::initSearch().
 451          * Cannot be called from user code.
 452          */
 453         explicit AnalysisNeighborhoodSearch(const ImplPointer &impl);
 454
 455         /*! \brief
 456          * Clears this search.
 457          *
 458          * Equivalent to \c "*this = AnalysisNeighborhoodSearch();".
 459          * Currently, this is necessary to avoid unnecessary memory allocation
 460          * if the previous search variable is still in scope when you want to
 461          * call AnalysisNeighborhood::initSearch() again.
 462          */
 463         void reset();
 464
 465         /*! \brief
 466          * Returns the searching algorithm that this search is using.
 467          *
 468          * The return value is never AnalysisNeighborhood::eSearchMode_Automatic.
 469          */
 470         AnalysisNeighborhood::SearchMode mode() const;
 471
 472         /*! \brief
 473          * Check whether a point is within a neighborhood.
 474          *
 475          * \param[in] positions  Set of test positions to use.
 476          * \returns   true if any of the test positions is within the cutoff of
 477          *     any reference position.
 478          */
 479         bool isWithin(const AnalysisNeighborhoodPositions &positions) const;
 480         /*! \brief
 481          * Calculates the minimum distance from the reference points.
 482          *
 483          * \param[in] positions  Set of test positions to use.
 484          * \returns   The distance to the nearest reference position, or the
 485          *     cutoff value if there are no reference positions within the
 486          *     cutoff.
 487          */
 488         real minimumDistance(const AnalysisNeighborhoodPositions &positions) const;
 489         /*! \brief
 490          * Finds the closest reference point.
 491          *
 492          * \param[in] positions  Set of test positions to use.
 493          * \returns   The reference index identifies the reference position
 494          *     that is closest to the test positions.
 495          *     The test index identifies the test position that is closest to
 496          *     the provided test position.  The returned pair is invalid if
 497          *     no reference position is within the cutoff.
 498          */
 499         AnalysisNeighborhoodPair
 500         nearestPoint(const AnalysisNeighborhoodPositions &positions) const;
 501
 502         /*! \brief
 503          * Start a search to find reference positions within a cutoff.
 504          *
 505          * \param[in] positions  Set of test positions to use.
 506          * \returns   Initialized search object to loop through all reference
 507          *     positions within the configured cutoff.
 508          * \throws    std::bad_alloc if out of memory.
 509          */
 510         AnalysisNeighborhoodPairSearch
 511         startPairSearch(const AnalysisNeighborhoodPositions &positions) const;
 512
 513     private:
 514         typedef internal::AnalysisNeighborhoodSearchImpl Impl;
 515
 516         ImplPointer             impl_;
 517 };
 518
 519 /*! \brief
 520  * Initialized neighborhood pair search with a fixed set of positions.
 521  *
 522  * This class is used to loop through pairs of neighbors within the cutoff
 523  * provided to AnalysisNeighborhood.  The following code demonstrates its use:
 524  * \code
 525    gmx::AnalysisNeighborhood       nb;
 526    nb.setCutoff(cutoff);
 527    gmx::AnalysisNeighborhoodPositions refPos(xref, nref);
 528    gmx::AnalysisNeighborhoodSearch search = nb.initSearch(pbc, refPos);
 529    gmx::AnalysisNeighborhoodPairSearch pairSearch = search.startPairSearch(selection);
 530    gmx::AnalysisNeighborhoodPair pair;
 531    while (pairSearch.findNextPair(&pair))
 532    {
 533        // <do something for each found pair the information in pair>
 534    }
 535  * \endcode
 536  *
 537  * It is not possible to use a single search object from multiple threads
 538  * concurrently.
 539  *
 540  * This class works like a pointer: copies of it point to the same search.
 541  * In general, avoid creating copies, and only use the copy/assignment support
 542  * for moving the variable around.  With C++11, this class would best be
 543  * movable.
 544  *
 545  * Methods in this class do not throw.
 546  *
 547  * \inpublicapi
 548  * \ingroup module_selection
 549  */
 550 class AnalysisNeighborhoodPairSearch
 551 {
 552     public:
 553         /*! \brief
 554          * Internal short-hand type for a pointer to the implementation class.
 555          *
 556          * See AnalysisNeighborhoodSearch::ImplPointer for rationale of using
 557          * shared_ptr and ownership semantics.
 558          */
 559         typedef std::shared_ptr<internal::AnalysisNeighborhoodPairSearchImpl>
 560             ImplPointer;
 561
 562         /*! \brief
 563          * Internally initialize the search.
 564          *
 565          * Used to implement AnalysisNeighborhoodSearch::startPairSearch().
 566          * Cannot be called from user code.
 567          */
 568         explicit AnalysisNeighborhoodPairSearch(const ImplPointer &impl);
 569
 570         /*! \brief
 571          * Finds the next pair within the cutoff.
 572          *
 573          * \param[out] pair  Information about the found pair.
 574          * \returns    false if there were no more pairs.
 575          *
 576          * If the method returns false, \p pair will be invalid.
 577          *
 578          * \see AnalysisNeighborhoodPair
 579          * \see AnalysisNeighborhoodSearch::startPairSearch()
 580          */
 581         bool findNextPair(AnalysisNeighborhoodPair *pair);
 582         /*! \brief
 583          * Skip remaining pairs for a test position in the search.
 584          *
 585          * When called after findNextPair(), makes subsequent calls to
 586          * findNextPair() skip any pairs that have the same test position as
 587          * that previously returned.
 588          * This is useful if the caller wants to search whether any reference
 589          * position within the cutoff satisfies some condition.  This method
 590          * can be used to skip remaining pairs after the first such position
 591          * has been found if the remaining pairs would not have an effect on
 592          * the outcome.
 593          */
 594         void skipRemainingPairsForTestPosition();
 595
 596     private:
 597         ImplPointer             impl_;
 598 };
 599
 600 } // namespace gmx
 601
 602 #endif