Sort all includes in src/gromacs
[alexxy/gromacs.git] / src / gromacs / selection / nbsearch.h
index a6127f76ba9ace9b64ebb412fdbaf743fc56f52b..fd7b41b129b2ea892f6f6fdc5ae988b5989de1bb 100644 (file)
@@ -1,10 +1,10 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2009,2010,2011,2012,2013, by the GROMACS development team, led by
- * David van der Spoel, Berk Hess, Erik Lindahl, and including many
- * others, as listed in the AUTHORS file in the top-level source
- * directory and at http://www.gromacs.org.
+ * Copyright (c) 2009,2010,2011,2012,2013,2014, 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.
  *
  * GROMACS is free software; you can redistribute it and/or
  * modify it under the terms of the GNU Lesser General Public License
 
 #include <boost/shared_ptr.hpp>
 
-#include "../legacyheaders/typedefs.h"
-#include "../utility/common.h"
-#include "../utility/gmxassert.h"
+#include "gromacs/math/vectypes.h"
+#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/common.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/real.h"
 
-#include "indexutil.h"
-
-struct gmx_ana_pos_t;
+struct t_blocka;
+struct t_pbc;
 
 namespace gmx
 {
@@ -71,6 +72,97 @@ class AnalysisNeighborhoodPairSearchImpl;
 class AnalysisNeighborhoodSearch;
 class AnalysisNeighborhoodPairSearch;
 
+/*! \brief
+ * Input positions for neighborhood searching.
+ *
+ * This class supports uniformly specifying sets of positions for various
+ * methods in the analysis neighborhood searching classes
+ * (AnalysisNeighborhood and AnalysisNeighborhoodSearch).
+ *
+ * Note that copies are not made: only a reference to the positions passed to
+ * the constructors are kept.  The caller is responsible to ensure that those
+ * positions remain in scope as long as the neighborhood search object requires
+ * access to them.
+ *
+ * Also note that in addition to constructors here, Selection and
+ * SelectionPosition provide conversions operators to this type.  It is done
+ * this way to not introduce a cyclic dependency between the selection code and
+ * the neighborhood search code, which in turn allows splitting this search
+ * code into a separate lower-level module if desired at some point.
+ *
+ * Methods in this class do not throw.
+ *
+ * \inpublicapi
+ * \ingroup module_selection
+ */
+class AnalysisNeighborhoodPositions
+{
+    public:
+        /*! \brief
+         * Initializes positions from a single position vector.
+         *
+         * For positions initialized this way, AnalysisNeighborhoodPair always
+         * returns zero in the corresponding index.
+         *
+         * This constructor is not explicit to allow directly passing an rvec
+         * to methods that accept positions.
+         */
+        AnalysisNeighborhoodPositions(const rvec &x)
+            : count_(1), index_(-1), x_(&x), exclusionIds_(NULL)
+        {
+        }
+        /*! \brief
+         * Initializes positions from an array of position vectors.
+         */
+        AnalysisNeighborhoodPositions(const rvec x[], int count)
+            : count_(count), index_(-1), x_(x), exclusionIds_(NULL)
+        {
+        }
+
+        /*! \brief
+         * Sets indices to use for mapping exclusions to these positions.
+         *
+         * The exclusion IDs can always be set, but they are ignored unless
+         * actual exclusions have been set with
+         * AnalysisNeighborhood::setTopologyExclusions().
+         */
+        AnalysisNeighborhoodPositions &
+        exclusionIds(ConstArrayRef<int> ids)
+        {
+            GMX_ASSERT(static_cast<int>(ids.size()) == count_,
+                       "Exclusion id array should match the number of positions");
+            exclusionIds_ = ids.data();
+            return *this;
+        }
+
+        /*! \brief
+         * Selects a single position to use from an array.
+         *
+         * If called, a single position from the array of positions passed to
+         * the constructor is used instead of the whole array.
+         * In contrast to the AnalysisNeighborhoodPositions(const rvec &)
+         * constructor, AnalysisNeighborhoodPair objects return \p index
+         * instead of zero.
+         */
+        AnalysisNeighborhoodPositions &selectSingleFromArray(int index)
+        {
+            GMX_ASSERT(index >= 0 && index < count_, "Invalid position index");
+            index_ = index;
+            return *this;
+        }
+
+    private:
+        int                     count_;
+        int                     index_;
+        const rvec             *x_;
+        const int              *exclusionIds_;
+
+        //! To access the positions for initialization.
+        friend class internal::AnalysisNeighborhoodSearchImpl;
+        //! To access the positions for initialization.
+        friend class internal::AnalysisNeighborhoodPairSearchImpl;
+};
+
 /*! \brief
  * Neighborhood searching for analysis tools.
  *
@@ -88,17 +180,14 @@ class AnalysisNeighborhoodPairSearch;
  *
  * initSearch() is thread-safe and can be called from multiple threads.  Each
  * call returns a different instance of the search object that can be used
- * independently of the others.  However, the returned search objects can only
- * be used within a single thread.  It is also possible to create multiple
- * concurrent searches within a single thread.
+ * independently of the others.  The returned AnalysisNeighborhoodSearch
+ * objects are also thread-safe, and can be used concurrently from multiple
+ * threads.  It is also possible to create multiple concurrent searches within
+ * a single thread.
  *
  * \todo
- * Support for exclusions.
- * The 4.5/4.6 C API had very low-level support for exclusions, which was not
- * very convenient to use, and hadn't been tested much.  The internal code that
- * it used to do the exclusion during the search itself is still there, but it
- * needs more thought on what would be a convenient way to initialize it.
- * Can be implemented once there is need for it in some calling code.
+ * Generalize the exclusion machinery to make it easier to use for other cases
+ * than atom-atom exclusions from the topology.
  *
  * \inpublicapi
  * \ingroup module_selection
@@ -122,7 +211,7 @@ class AnalysisNeighborhood
         ~AnalysisNeighborhood();
 
         /*! \brief
-         * Set cutoff distance for the neighborhood searching.
+         * Sets cutoff distance for the neighborhood searching.
          *
          * \param[in]  cutoff Cutoff distance for the search
          *   (<=0 stands for no cutoff).
@@ -133,6 +222,34 @@ class AnalysisNeighborhood
          * Does not throw.
          */
         void setCutoff(real cutoff);
+        /*! \brief
+         * Sets the search to only happen in the XY plane.
+         *
+         * Z component of the coordinates is not used in the searching,
+         * and returned distances are computed in the XY plane.
+         * Only boxes with the third box vector parallel to the Z axis are
+         * currently implemented.
+         *
+         * Does not throw.
+         */
+        void setXYMode(bool bXY);
+        /*! \brief
+         * Sets atom exclusions from a topology.
+         *
+         * The \p excls structure specifies the exclusions from test positions
+         * to reference positions, i.e., a block starting at `excls->index[i]`
+         * specifies the exclusions for test position `i`, and the indices in
+         * `excls->a` are indices of the reference positions.  If `excls->nr`
+         * is smaller than a test position id, then such test positions do not
+         * have any exclusions.
+         * It is assumed that the indices within a block of indices in
+         * `excls->a` is ascending.
+         *
+         * Does not throw.
+         *
+         * \see AnalysisNeighborhoodPositions::exclusionIds()
+         */
+        void setTopologyExclusions(const t_blocka *excls);
         /*! \brief
          * Sets the algorithm to use for searching.
          *
@@ -152,26 +269,18 @@ class AnalysisNeighborhood
         /*! \brief
          * Initializes neighborhood search for a set of positions.
          *
-         * \param[in] pbc PBC information for the frame.
-         * \param[in] n   Number of reference positions for the frame.
-         * \param[in] x   \p n reference positions for the frame.
+         * \param[in] pbc        PBC information for the frame.
+         * \param[in] positions  Set of reference positions to use.
          * \returns   Search object that can be used to find positions from
          *      \p x within the given cutoff.
          * \throws    std::bad_alloc if out of memory.
-         */
-        AnalysisNeighborhoodSearch
-        initSearch(const t_pbc *pbc, int n, const rvec x[]);
-        /*! \brief
-         * Initializes neighborhood search for a set of positions.
          *
-         * \param[in] pbc PBC information for the frame.
-         * \param[in] p   Reference positions for the frame.
-         * \returns   Search object that can be used to find positions from
-         *      \p p within the given cutoff.
-         * \throws    std::bad_alloc if out of memory.
+         * Currently, the input positions cannot use
+         * AnalysisNeighborhoodPositions::selectSingleFromArray().
          */
         AnalysisNeighborhoodSearch
-        initSearch(const t_pbc *pbc, const gmx_ana_pos_t *p);
+        initSearch(const t_pbc                         *pbc,
+                   const AnalysisNeighborhoodPositions &positions);
 
     private:
         class Impl;
@@ -191,10 +300,10 @@ class AnalysisNeighborhoodPair
 {
     public:
         //! Initializes an invalid pair.
-        AnalysisNeighborhoodPair() : refIndex_(-1), testIndex_(0) {}
+        AnalysisNeighborhoodPair() : refIndex_(-1), testIndex_(0), distance2_(0.0) {}
         //! Initializes a pair object with the given data.
-        AnalysisNeighborhoodPair(int refIndex, int testIndex)
-            : refIndex_(refIndex), testIndex_(testIndex)
+        AnalysisNeighborhoodPair(int refIndex, int testIndex, real distance2)
+            : refIndex_(refIndex), testIndex_(testIndex), distance2_(distance2)
         {
         }
 
@@ -228,10 +337,19 @@ class AnalysisNeighborhoodPair
             GMX_ASSERT(isValid(), "Accessing invalid object");
             return testIndex_;
         }
+        /*! \brief
+         * Returns the squared distance between the pair of positions.
+         */
+        real distance2() const
+        {
+            GMX_ASSERT(isValid(), "Accessing invalid object");
+            return distance2_;
+        }
 
     private:
         int                     refIndex_;
         int                     testIndex_;
+        real                    distance2_;
 };
 
 /*! \brief
@@ -240,16 +358,16 @@ class AnalysisNeighborhoodPair
  * An instance of this class is obtained through
  * AnalysisNeighborhood::initSearch(), and can be used to do multiple searches
  * against the provided set of reference positions.
- * Currently, it is not possible to call startPairSearch() while a previous
- * AnalysisNeighborhoodPairSearch object obtained from startPairSearch() of the
- * same instance still exists.
+ * It is possible to create concurrent pair searches (including from different
+ * threads), as well as call other methods in this class while a pair search is
+ * in progress.
  *
  * This class works like a pointer: copies of it point to the same search.
  * In general, avoid creating copies, and only use the copy/assignment support
  * for moving the variable around.  With C++11, this class would best be
  * movable.
  *
- * Methods in this class do not throw.
+ * Methods in this class do not throw unless otherwise indicated.
  *
  * \todo
  * Make it such that reset() is not necessary to call in code that repeatedly
@@ -314,86 +432,47 @@ class AnalysisNeighborhoodSearch
         /*! \brief
          * Check whether a point is within a neighborhood.
          *
-         * \param[in] x  Test position.
-         * \returns   true if the test position is within the cutoff of any
-         *     reference position.
-         */
-        bool isWithin(const rvec x) const;
-        /*! \brief
-         * Check whether a point is within a neighborhood.
-         *
-         * \param[in] p  Test positions.
-         * \param[in] i  Use the i'th position in \p p for testing.
-         * \returns   true if the test position is within the cutoff of any
-         *     reference position.
-         */
-        bool isWithin(const gmx_ana_pos_t *p, int i) const;
-        /*! \brief
-         * Calculates the minimum distance from the reference points.
-         *
-         * \param[in] x  Test position.
-         * \returns   The distance to the nearest reference position, or the
-         *     cutoff value if there are no reference positions within the
-         *     cutoff.
+         * \param[in] positions  Set of test positions to use.
+         * \returns   true if any of the test positions is within the cutoff of
+         *     any reference position.
          */
-        real minimumDistance(const rvec x) const;
+        bool isWithin(const AnalysisNeighborhoodPositions &positions) const;
         /*! \brief
          * Calculates the minimum distance from the reference points.
          *
-         * \param[in] p  Test positions.
-         * \param[in] i  Use the i'th position in \p p for testing.
+         * \param[in] positions  Set of test positions to use.
          * \returns   The distance to the nearest reference position, or the
          *     cutoff value if there are no reference positions within the
          *     cutoff.
          */
-        real minimumDistance(const gmx_ana_pos_t *p, int i) const;
-        /*! \brief
-         * Finds the closest reference point.
-         *
-         * \param[in] x  Test position.
-         * \returns   The reference index identifies the reference position
-         *     that is closest to the test position.
-         *     The test index is always zero.  The returned pair is invalid if
-         *     no reference position is within the cutoff.
-         */
-        AnalysisNeighborhoodPair nearestPoint(const rvec x) const;
+        real minimumDistance(const AnalysisNeighborhoodPositions &positions) const;
         /*! \brief
          * Finds the closest reference point.
          *
-         * \param[in] p  Test positions.
-         * \param[in] i  Use the i'th position in \p p for testing.
+         * \param[in] positions  Set of test positions to use.
          * \returns   The reference index identifies the reference position
-         *     that is closest to the test position.
-         *     The test index is always \p i.  The returned pair is invalid if
+         *     that is closest to the test positions.
+         *     The test index identifies the test position that is closest to
+         *     the provided test position.  The returned pair is invalid if
          *     no reference position is within the cutoff.
          */
-        AnalysisNeighborhoodPair nearestPoint(const gmx_ana_pos_t *p, int i) const;
+        AnalysisNeighborhoodPair
+        nearestPoint(const AnalysisNeighborhoodPositions &positions) const;
 
         /*! \brief
          * Start a search to find reference positions within a cutoff.
          *
-         * \param[in] x  Test position to search the neighbors for.
+         * \param[in] positions  Set of test positions to use.
          * \returns   Initialized search object to loop through all reference
          *     positions within the configured cutoff.
-         *
-         * In the AnalysisNeighborhoodPair objects returned by the search, the
-         * test index is always zero.
-         */
-        AnalysisNeighborhoodPairSearch startPairSearch(const rvec x);
-        /*! \brief
-         * Start a search to find reference positions within a cutoff.
-         *
-         * \param[in] p  Test positions.
-         * \param[in] i  Use the i'th position in \p p for testing.
-         * \returns   Initialized search object to loop through all reference
-         *     positions within the configured cutoff.
-         *
-         * In the AnalysisNeighborhoodPair objects returned by the search, the
-         * test index is always \p i.
+         * \throws    std::bad_alloc if out of memory.
          */
-        AnalysisNeighborhoodPairSearch startPairSearch(const gmx_ana_pos_t *p, int i);
+        AnalysisNeighborhoodPairSearch
+        startPairSearch(const AnalysisNeighborhoodPositions &positions) const;
 
     private:
+        typedef internal::AnalysisNeighborhoodSearchImpl Impl;
+
         ImplPointer             impl_;
 };
 
@@ -405,8 +484,9 @@ class AnalysisNeighborhoodSearch
  * \code
    gmx::AnalysisNeighborhood       nb;
    nb.setCutoff(cutoff);
-   gmx::AnalysisNeighborhoodSearch search = nb.initSearch(pbc, nref, xref);
-   gmx::AnalysisNeighborhoodPairSearch pairSearch = search.startPairSearch(x);
+   gmx::AnalysisNeighborhoodPositions refPos(xref, nref);
+   gmx::AnalysisNeighborhoodSearch search = nb.initSearch(pbc, refPos);
+   gmx::AnalysisNeighborhoodPairSearch pairSearch = search.startPairSearch(selection);
    gmx::AnalysisNeighborhoodPair pair;
    while (pairSearch.findNextPair(&pair))
    {
@@ -459,6 +539,19 @@ class AnalysisNeighborhoodPairSearch
          * \see AnalysisNeighborhoodSearch::startPairSearch()
          */
         bool findNextPair(AnalysisNeighborhoodPair *pair);
+        /*! \brief
+         * Skip remaining pairs for a test position in the search.
+         *
+         * When called after findNextPair(), makes subsequent calls to
+         * findNextPair() skip any pairs that have the same test position as
+         * that previously returned.
+         * This is useful if the caller wants to search whether any reference
+         * position within the cutoff satisfies some condition.  This method
+         * can be used to skip remaining pairs after the first such position
+         * has been found if the remaining pairs would not have an effect on
+         * the outcome.
+         */
+        void skipRemainingPairsForTestPosition();
 
     private:
         ImplPointer             impl_;