/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2011,2012, by the GROMACS development team, led by
+ * Copyright (c) 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.
* \param[in] selections Thread-local selection collection.
* \param[in] cutoff Cutoff distance for the search
* (<=0 stands for no cutoff).
- * \param[in] posCount Maximum number of reference particles.
*/
- ModuleData(TrajectoryAnalysisModule *module,
+ ModuleData(TrajectoryAnalysisModule *module,
const AnalysisDataParallelOptions &opt,
- const SelectionCollection &selections,
- double cutoff, int posCount)
- : TrajectoryAnalysisModuleData(module, opt, selections),
- nb_(cutoff, posCount)
+ const SelectionCollection &selections,
+ double cutoff)
+ : TrajectoryAnalysisModuleData(module, opt, selections)
{
+ nb_.setCutoff(cutoff);
}
virtual void finish()
}
//! Neighborhood search data for distance calculation.
- NeighborhoodSearch nb_;
+ AnalysisNeighborhood nb_;
};
const SelectionCollection &selections)
{
return TrajectoryAnalysisModuleDataPointer(
- new ModuleData(this, opt, selections, cutoff_, refsel_.posCount()));
+ new ModuleData(this, opt, selections, cutoff_));
}
AnalysisTemplate::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
TrajectoryAnalysisModuleData *pdata)
{
- AnalysisDataHandle dh = pdata->dataHandle(data_);
- NeighborhoodSearch &nb = static_cast<ModuleData *>(pdata)->nb_;
- const Selection &refsel = pdata->parallelSelection(refsel_);
+ AnalysisDataHandle dh = pdata->dataHandle(data_);
+ AnalysisNeighborhood &nb = static_cast<ModuleData *>(pdata)->nb_;
+ const Selection &refsel = pdata->parallelSelection(refsel_);
- nb.init(pbc, refsel.positions());
+ AnalysisNeighborhoodSearch nbsearch =
+ nb.initSearch(pbc, refsel.positions());
dh.startFrame(frnr, fr.time);
for (size_t g = 0; g < sel_.size(); ++g)
{
for (int i = 0; i < nr; ++i)
{
SelectionPosition p = sel.position(i);
- frave += nb.minimumDistance(p.x());
+ frave += nbsearch.minimumDistance(p.x());
}
frave /= nr;
dh.setPoint(g, frave);
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2009,2010,2011,2012, by the GROMACS development team, led by
+ * 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.
#include <math.h>
+#include <algorithm>
+
#include "gromacs/legacyheaders/smalloc.h"
#include "gromacs/legacyheaders/typedefs.h"
#include "gromacs/legacyheaders/pbc.h"
#include "gromacs/selection/nbsearch.h"
#include "gromacs/selection/position.h"
+#include "gromacs/utility/gmxassert.h"
/*! \internal \brief
* Data structure for neighborhood searches.
bool bGrid;
/** Array allocated for storing in-unit-cell reference positions. */
rvec *xref_alloc;
+ /** Allocation count for xref_alloc. */
+ int xref_nalloc;
/** false if the box is rectangular. */
bool bTric;
/** Box vectors of a single grid cell. */
d->exclind = 0;
d->xref_alloc = NULL;
+ d->xref_nalloc = 0;
d->ncells = 0;
d->ncatoms = NULL;
d->catom = NULL;
{
int i;
- if (!d->xref_alloc)
+ if (d->xref_nalloc < n)
{
- snew(d->xref_alloc, d->maxnref);
+ const int allocCount = std::max(d->maxnref, n);
+ srenew(d->xref_alloc, allocCount);
+ d->xref_nalloc = allocCount;
}
d->xref = d->xref_alloc;
grid_clear_cells(d);
gmx_ana_nbsearch_pos_init(gmx_ana_nbsearch_t *d, t_pbc *pbc, const gmx_ana_pos_t *p)
{
gmx_ana_nbsearch_init(d, pbc, p->nr, p->x);
- d->refid = (p->nr < d->maxnref ? p->m.refid : NULL);
+ d->refid = p->m.refid;
}
/*!
d->prevnbi = 0;
d->prevcai = -1;
}
- else
- {
- d->previ = -1;
- }
+ d->previ = -1;
d->exclind = 0;
}
mindist_action(gmx_ana_nbsearch_t *d, int i, real r2)
{
d->cutoff2 = r2;
+ d->previ = i;
return false;
}
*jp = -1;
return false;
}
+
+namespace gmx
+{
+
+namespace internal
+{
+
+/********************************************************************
+ * AnalysisNeighborhoodPairSearchImpl
+ */
+
+class AnalysisNeighborhoodPairSearchImpl
+{
+ public:
+ AnalysisNeighborhoodPairSearchImpl()
+ : nb_(NULL), testIndex_(0)
+ {
+ }
+
+ gmx_ana_nbsearch_t *nb_;
+ int testIndex_;
+};
+
+/********************************************************************
+ * AnalysisNeighborhoodSearchImpl
+ */
+
+class AnalysisNeighborhoodSearchImpl
+{
+ public:
+ typedef AnalysisNeighborhoodPairSearch::ImplPointer
+ PairSearchImplPointer;
+
+ AnalysisNeighborhoodSearchImpl()
+ : nb_(NULL), pairSearch_(new AnalysisNeighborhoodPairSearchImpl)
+ {
+ }
+ ~AnalysisNeighborhoodSearchImpl()
+ {
+ GMX_RELEASE_ASSERT(pairSearch_.unique(),
+ "Dangling AnalysisNeighborhoodPairSearch reference");
+ if (nb_ != NULL)
+ {
+ gmx_ana_nbsearch_free(nb_);
+ }
+ }
+
+ gmx_ana_nbsearch_t *nb_;
+ PairSearchImplPointer pairSearch_;
+};
+
+} // namespace internal
+
+/********************************************************************
+ * AnalysisNeighborhood::Impl
+ */
+
+class AnalysisNeighborhood::Impl
+{
+ public:
+ typedef AnalysisNeighborhoodSearch::ImplPointer SearchImplPointer;
+
+ Impl()
+ : search_(new internal::AnalysisNeighborhoodSearchImpl), bTryGrid_(false)
+ {
+ }
+ ~Impl()
+ {
+ GMX_RELEASE_ASSERT(search_.unique(),
+ "Dangling AnalysisNeighborhoodSearch reference");
+ }
+
+ SearchImplPointer search_;
+ bool bTryGrid_;
+};
+
+/********************************************************************
+ * AnalysisNeighborhood
+ */
+
+AnalysisNeighborhood::AnalysisNeighborhood()
+ : impl_(new Impl)
+{
+}
+
+AnalysisNeighborhood::~AnalysisNeighborhood()
+{
+}
+
+void AnalysisNeighborhood::setCutoff(real cutoff)
+{
+ GMX_RELEASE_ASSERT(impl_->search_->nb_ == NULL,
+ "Multiple calls to setCutoff() not currently supported");
+ impl_->search_->nb_ = gmx_ana_nbsearch_create(cutoff, 0);
+ impl_->search_->pairSearch_->nb_ = impl_->search_->nb_;
+ impl_->bTryGrid_ = impl_->search_->nb_->bTryGrid;
+}
+
+void AnalysisNeighborhood::setMode(SearchMode mode)
+{
+ GMX_RELEASE_ASSERT(impl_->search_->nb_ != NULL,
+ "setParameters() not called");
+ // TODO: Actually implement eSearchMode_Grid.
+ switch (mode)
+ {
+ case eSearchMode_Automatic:
+ impl_->search_->nb_->bTryGrid = impl_->bTryGrid_;
+ break;
+ case eSearchMode_Simple:
+ impl_->search_->nb_->bTryGrid = false;
+ break;
+ case eSearchMode_Grid:
+ impl_->search_->nb_->bTryGrid = impl_->bTryGrid_;
+ break;
+ }
+}
+
+AnalysisNeighborhood::SearchMode AnalysisNeighborhood::mode() const
+{
+ GMX_RELEASE_ASSERT(impl_->search_->nb_ != NULL,
+ "setParameters() not called");
+ if (!impl_->search_->nb_->bTryGrid)
+ {
+ return eSearchMode_Simple;
+ }
+ return eSearchMode_Automatic;
+}
+
+AnalysisNeighborhoodSearch
+AnalysisNeighborhood::initSearch(const t_pbc *pbc, int n, const rvec x[])
+{
+ GMX_RELEASE_ASSERT(impl_->search_->nb_ != NULL,
+ "setParameters() not called");
+ GMX_RELEASE_ASSERT(impl_->search_.unique(),
+ "Multiple concurrent searches not currently supported");
+ gmx_ana_nbsearch_init(impl_->search_->nb_, const_cast<t_pbc *>(pbc), n, x);
+ return AnalysisNeighborhoodSearch(impl_->search_);
+}
+
+AnalysisNeighborhoodSearch
+AnalysisNeighborhood::initSearch(const t_pbc *pbc, const gmx_ana_pos_t *p)
+{
+ GMX_RELEASE_ASSERT(impl_->search_->nb_ != NULL,
+ "setParameters() not called");
+ GMX_RELEASE_ASSERT(impl_->search_.unique(),
+ "Multiple concurrent searches not currently supported");
+ gmx_ana_nbsearch_pos_init(impl_->search_->nb_, const_cast<t_pbc *>(pbc), p);
+ return AnalysisNeighborhoodSearch(impl_->search_);
+}
+
+/********************************************************************
+ * AnalysisNeighborhoodSearch
+ */
+
+AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch()
+{
+}
+
+AnalysisNeighborhoodSearch::AnalysisNeighborhoodSearch(const ImplPointer &impl)
+ : impl_(impl)
+{
+}
+
+void AnalysisNeighborhoodSearch::reset()
+{
+ impl_.reset();
+}
+
+AnalysisNeighborhood::SearchMode AnalysisNeighborhoodSearch::mode() const
+{
+ GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
+ return (impl_->nb_->bGrid
+ ? AnalysisNeighborhood::eSearchMode_Grid
+ : AnalysisNeighborhood::eSearchMode_Simple);
+}
+
+bool AnalysisNeighborhoodSearch::isWithin(const rvec x) const
+{
+ GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
+ GMX_RELEASE_ASSERT(impl_->pairSearch_.unique(),
+ "Calls to other methods concurrently with a pair search "
+ "not currently supported");
+ return gmx_ana_nbsearch_is_within(impl_->nb_, x);
+}
+
+bool AnalysisNeighborhoodSearch::isWithin(const gmx_ana_pos_t *p, int i) const
+{
+ GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
+ GMX_RELEASE_ASSERT(impl_->pairSearch_.unique(),
+ "Calls to other methods concurrently with a pair search "
+ "not currently supported");
+ return gmx_ana_nbsearch_pos_is_within(impl_->nb_, p, i);
+}
+
+real AnalysisNeighborhoodSearch::minimumDistance(const rvec x) const
+{
+ GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
+ GMX_RELEASE_ASSERT(impl_->pairSearch_.unique(),
+ "Calls to other methods concurrently with a pair search "
+ "not currently supported");
+ return gmx_ana_nbsearch_mindist(impl_->nb_, x);
+}
+
+real AnalysisNeighborhoodSearch::minimumDistance(const gmx_ana_pos_t *p, int i) const
+{
+ GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
+ GMX_RELEASE_ASSERT(impl_->pairSearch_.unique(),
+ "Calls to other methods concurrently with a pair search "
+ "not currently supported");
+ return gmx_ana_nbsearch_pos_mindist(impl_->nb_, p, i);
+}
+
+AnalysisNeighborhoodPair
+AnalysisNeighborhoodSearch::nearestPoint(const rvec x) const
+{
+ GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
+ GMX_RELEASE_ASSERT(impl_->pairSearch_.unique(),
+ "Calls to other methods concurrently with a pair search "
+ "not currently supported");
+ (void)gmx_ana_nbsearch_mindist(impl_->nb_, x);
+ return AnalysisNeighborhoodPair(impl_->nb_->previ, 0);
+}
+
+AnalysisNeighborhoodPair
+AnalysisNeighborhoodSearch::nearestPoint(const gmx_ana_pos_t *p, int i) const
+{
+ GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
+ GMX_RELEASE_ASSERT(impl_->pairSearch_.unique(),
+ "Calls to other methods concurrently with a pair search "
+ "not currently supported");
+ (void)gmx_ana_nbsearch_pos_mindist(impl_->nb_, p, i);
+ return AnalysisNeighborhoodPair(impl_->nb_->previ, i);
+}
+
+AnalysisNeighborhoodPairSearch
+AnalysisNeighborhoodSearch::startPairSearch(const rvec x)
+{
+ GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
+ GMX_RELEASE_ASSERT(impl_->pairSearch_.unique(),
+ "Multiple concurrent searches not currently supported");
+ grid_search_start(impl_->nb_, x);
+ impl_->pairSearch_->testIndex_ = 0;
+ return AnalysisNeighborhoodPairSearch(impl_->pairSearch_);
+}
+
+AnalysisNeighborhoodPairSearch
+AnalysisNeighborhoodSearch::startPairSearch(const gmx_ana_pos_t *p, int i)
+{
+ GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
+ GMX_RELEASE_ASSERT(impl_->pairSearch_.unique(),
+ "Multiple concurrent searches not currently supported");
+ grid_search_start(impl_->nb_, p->x[i]);
+ impl_->pairSearch_->testIndex_ = i;
+ return AnalysisNeighborhoodPairSearch(impl_->pairSearch_);
+}
+
+/********************************************************************
+ * AnalysisNeighborhoodPairSearch
+ */
+
+AnalysisNeighborhoodPairSearch::AnalysisNeighborhoodPairSearch(
+ const ImplPointer &impl)
+ : impl_(impl)
+{
+}
+
+bool AnalysisNeighborhoodPairSearch::findNextPair(AnalysisNeighborhoodPair *pair)
+{
+ if (grid_search(impl_->nb_, &within_action))
+ {
+ *pair = AnalysisNeighborhoodPair(impl_->nb_->previ, impl_->testIndex_);
+ return true;
+ }
+ *pair = AnalysisNeighborhoodPair();
+ return false;
+}
+
+} // namespace gmx
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2009,2010,2011,2012, by the GROMACS development team, led by
+ * 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.
* The library also uses the functions internally.
*
* \author Teemu Murtola <teemu.murtola@gmail.com>
+ * \inpublicapi
* \ingroup module_selection
*/
#ifndef GMX_SELECTION_NBSEARCH_H
#define GMX_SELECTION_NBSEARCH_H
+#include <boost/shared_ptr.hpp>
+
#include "../legacyheaders/typedefs.h"
+#include "../utility/common.h"
+#include "../utility/gmxassert.h"
#include "indexutil.h"
namespace gmx
{
-/*
- * C++ wrapper for neighborhood searching.
+namespace internal
+{
+class AnalysisNeighborhoodSearchImpl;
+class AnalysisNeighborhoodPairSearchImpl;
+};
+
+class AnalysisNeighborhoodSearch;
+class AnalysisNeighborhoodPairSearch;
+
+/*! \brief
+ * Neighborhood searching for analysis tools.
+ *
+ * This class implements neighborhood searching routines for analysis tools.
+ * The emphasis is in flexibility and ease of use; one main driver is to have
+ * a common implementation of grid-based searching to avoid replicating this in
+ * multiple tools (and to make more tools take advantage of the significant
+ * performance improvement this allows).
*
+ * To use the search, create an object of this type, call setCutoff() to
+ * initialize it, and then repeatedly call initSearch() to start a search with
+ * different sets of reference positions. For each set of reference positions,
+ * use methods in the returned AnalysisNeighborhoodSearch to find the reference
+ * positions that are within the given cutoff from a provided position.
+ *
+ * \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.
+ *
+ * \inpublicapi
+ * \ingroup module_selection
*/
-class NeighborhoodSearch
+class AnalysisNeighborhood
{
public:
- NeighborhoodSearch(real cutoff, int maxn)
- : d_(gmx_ana_nbsearch_create(cutoff, maxn))
+ //! Searching algorithm to use.
+ enum SearchMode
+ {
+ //! Select algorithm based on heuristic efficiency considerations.
+ eSearchMode_Automatic,
+ //! Use a simple loop over all pairs.
+ eSearchMode_Simple,
+ //! Use grid-based searching whenever possible.
+ eSearchMode_Grid
+ };
+
+ //! Creates an uninitialized neighborhood search.
+ AnalysisNeighborhood();
+ ~AnalysisNeighborhood();
+
+ /*! \brief
+ * Set cutoff distance for the neighborhood searching.
+ *
+ * \param[in] cutoff Cutoff distance for the search
+ * (<=0 stands for no cutoff).
+ *
+ * Currently, needs to be called exactly once, before calling any other
+ * method.
+ */
+ void setCutoff(real cutoff);
+ /*! \brief
+ * Sets the algorithm to use for searching.
+ *
+ * \param[in] mode Search mode to use.
+ *
+ * Note that if \p mode is \ref eSearchMode_Grid, it is still only a
+ * suggestion: grid-based searching may not be possible with the
+ * provided input, in which case a simple search is still used.
+ * This is mainly useful for testing purposes to force a mode.
+ *
+ * Does not throw.
+ */
+ void setMode(SearchMode mode);
+ //! Returns the currently active search mode.
+ SearchMode mode() const;
+
+ /*! \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.
+ * \returns Search object that can be used to find positions from
+ * \p x within the given cutoff.
+ */
+ 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.
+ */
+ AnalysisNeighborhoodSearch
+ initSearch(const t_pbc *pbc, const gmx_ana_pos_t *p);
+
+ private:
+ class Impl;
+
+ PrivateImplPointer<Impl> impl_;
+};
+
+/*! \brief
+ * Value type to represent a pair of positions found in neighborhood searching.
+ *
+ * Methods in this class do not throw.
+ *
+ * \inpublicapi
+ * \ingroup module_selection
+ */
+class AnalysisNeighborhoodPair
+{
+ public:
+ //! Initializes an invalid pair.
+ AnalysisNeighborhoodPair() : refIndex_(-1), testIndex_(0) {}
+ //! Initializes a pair object with the given data.
+ AnalysisNeighborhoodPair(int refIndex, int testIndex)
+ : refIndex_(refIndex), testIndex_(testIndex)
+ {
+ }
+
+ /*! \brief
+ * Whether this pair is valid.
+ *
+ * If isValid() returns false, other methods should not be called.
+ */
+ bool isValid() const { return refIndex_ >= 0; }
+
+ /*! \brief
+ * Returns the index of the reference position in the pair.
+ *
+ * This index is always the index into the position array provided to
+ * AnalysisNeighborhood::initSearch().
+ */
+ int refIndex() const
+ {
+ GMX_ASSERT(isValid(), "Accessing invalid object");
+ return refIndex_;
+ }
+ /*! \brief
+ * Returns the index of the test position in the pair.
+ *
+ * The contents of this index depends on the context (method call) that
+ * produces the pair.
+ * If there was no array in the call, this index is zero.
+ */
+ int testIndex() const
{
+ GMX_ASSERT(isValid(), "Accessing invalid object");
+ return testIndex_;
}
- ~NeighborhoodSearch() { gmx_ana_nbsearch_free(d_); }
- void init(t_pbc *pbc, int n, const rvec x[])
- { gmx_ana_nbsearch_init(d_, pbc, n, x); }
+ private:
+ int refIndex_;
+ int testIndex_;
+};
- void init(t_pbc *pbc, const gmx_ana_pos_t *p)
- { gmx_ana_nbsearch_pos_init(d_, pbc, p); }
+/*! \brief
+ * Initialized neighborhood search with a fixed set of reference positions.
+ *
+ * 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 any method in this class while an
+ * AnalysisNeighborhoodPairSearch object obtained from startPairSearch() of the
+ * same instance exists.
+ *
+ * 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.
+ *
+ * \todo
+ * Make it such that reset() is not necessary to call in code that repeatedly
+ * assigns the result of AnalysisNeighborhood::initSearch() to the same
+ * variable (see sm_distance.cpp).
+ *
+ * \todo
+ * Consider merging nearestPoint() and minimumDistance() by adding the distance
+ * to AnalysisNeighborhoodPair.
+ *
+ * \inpublicapi
+ * \ingroup module_selection
+ */
+class AnalysisNeighborhoodSearch
+{
+ public:
+ /*! \brief
+ * Internal short-hand type for a pointer to the implementation class.
+ *
+ * shared_ptr is used here to automatically keep a reference count to
+ * track whether an implementation class is still used outside the
+ * AnalysisNeighborhood object. Ownership currently always stays with
+ * AnalysisNeighborhood; it always keeps one instance of the pointer.
+ */
+ typedef boost::shared_ptr<internal::AnalysisNeighborhoodSearchImpl>
+ ImplPointer;
- void setExclusions(int nexcl, atom_id *excl)
- { gmx_ana_nbsearch_set_excl(d_, nexcl, excl); }
+ /*! \brief
+ * Initializes an invalid search.
+ *
+ * Such an object cannot be used for searching. It needs to be
+ * assigned a value from AnalysisNeighborhood::initSearch() before it
+ * can be used. Provided to allow declaring a variable to hold the
+ * search before calling AnalysisNeighborhood::initSearch().
+ */
+ AnalysisNeighborhoodSearch();
+ /*! \brief
+ * Internally initialize the search.
+ *
+ * Used to implement AnalysisNeighborhood::initSearch().
+ * Cannot be called from user code.
+ */
+ explicit AnalysisNeighborhoodSearch(const ImplPointer &impl);
+ /*! \brief
+ * Clears this search.
+ *
+ * Equivalent to \c "*this = AnalysisNeighborhoodSearch();".
+ * Currently, this is required if the previous search variable is still
+ * in scope and you want to call AnalysisNeighborhood::initSearch()
+ * again.
+ */
+ void reset();
- bool isWithin(const rvec x)
- { return gmx_ana_nbsearch_is_within(d_, x); }
+ /*! \brief
+ * Returns the searching algorithm that this search is using.
+ *
+ * The return value is never AnalysisNeighborhood::eSearchMode_Automatic.
+ */
+ AnalysisNeighborhood::SearchMode mode() const;
- bool isWithin(const gmx_ana_pos_t *p, int i)
- { return gmx_ana_nbsearch_pos_is_within(d_, p, i); }
+ /*! \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.
+ */
+ real minimumDistance(const rvec x) 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.
+ * \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;
+ /*! \brief
+ * Finds the closest reference point.
+ *
+ * \param[in] p Test positions.
+ * \param[in] i Use the i'th position in \p p for testing.
+ * \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
+ * no reference position is within the cutoff.
+ */
+ AnalysisNeighborhoodPair nearestPoint(const gmx_ana_pos_t *p, int i) const;
- real minimumDistance(const rvec x)
- { return gmx_ana_nbsearch_mindist(d_, x); }
+ /*! \brief
+ * Start a search to find reference positions within a cutoff.
+ *
+ * \param[in] x Test position to search the neighbors for.
+ * \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.
+ */
+ AnalysisNeighborhoodPairSearch startPairSearch(const gmx_ana_pos_t *p, int i);
- real minimumDistance(const gmx_ana_pos_t *p, int i)
- { return gmx_ana_nbsearch_pos_mindist(d_, p, i); }
+ private:
+ ImplPointer impl_;
+};
- bool firstWithin(const rvec x, int *jp)
- { return gmx_ana_nbsearch_first_within(d_, x, jp); }
+/*! \brief
+ * Initialized neighborhood pair search with a fixed set of positions.
+ *
+ * This class is used to loop through pairs of neighbors within the cutoff
+ * provided to AnalysisNeighborhood. The following code demonstrates its use:
+ * \code
+ gmx::AnalysisNeighborhood nb;
+ nb.setCutoff(cutoff);
+ gmx::AnalysisNeighborhoodSearch search = nb.initSearch(pbc, nref, xref);
+ gmx::AnalysisNeighborhoodPairSearch pairSearch = search.startPairSearch(x);
+ gmx::AnalysisNeighborhoodPair pair;
+ while (pairSearch.findNextPair(&pair))
+ {
+ // <do something for each found pair the information in pair>
+ }
+ * \endcode
+ *
+ * It is not possible to use a single search object from multiple threads
+ * concurrently.
+ *
+ * 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.
+ *
+ * \inpublicapi
+ * \ingroup module_selection
+ */
+class AnalysisNeighborhoodPairSearch
+{
+ public:
+ /*! \brief
+ * Internal short-hand type for a pointer to the implementation class.
+ *
+ * See AnalysisNeighborhoodSearch::ImplPointer for rationale of using
+ * shared_ptr and ownership semantics.
+ */
+ typedef boost::shared_ptr<internal::AnalysisNeighborhoodPairSearchImpl>
+ ImplPointer;
- bool firstWithin(const gmx_ana_pos_t *p, int i, int *jp)
- { return gmx_ana_nbsearch_pos_first_within(d_, p, i, jp); }
+ /*! \brief
+ * Internally initialize the search.
+ *
+ * Used to implement AnalysisNeighborhoodSearch::startPairSearch().
+ * Cannot be called from user code.
+ */
+ explicit AnalysisNeighborhoodPairSearch(const ImplPointer &impl);
- bool nextWithin(int *jp)
- { return gmx_ana_nbsearch_next_within(d_, jp); }
+ /*! \brief
+ * Finds the next pair within the cutoff.
+ *
+ * \param[out] pair Information about the found pair.
+ * \returns false if there were no more pairs.
+ *
+ * If the method returns false, \p pair will be invalid.
+ *
+ * \see AnalysisNeighborhoodPair
+ * \see AnalysisNeighborhoodSearch::startPairSearch()
+ */
+ bool findNextPair(AnalysisNeighborhoodPair *pair);
private:
- gmx_ana_nbsearch_t *d_;
+ ImplPointer impl_;
};
} // namespace gmx
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2009,2010,2011,2012, by the GROMACS development team, led by
+ * 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.
*/
#include "gromacs/legacyheaders/macros.h"
#include "gromacs/legacyheaders/pbc.h"
-#include "gromacs/legacyheaders/smalloc.h"
#include "gromacs/legacyheaders/vec.h"
#include "gromacs/selection/nbsearch.h"
*
* The same data structure is used by all the distance-based methods.
*/
-typedef struct
+struct t_methoddata_distance
{
+ t_methoddata_distance() : cutoff(-1.0)
+ {
+ gmx_ana_pos_clear(&p);
+ }
+
/** Cutoff distance. */
- real cutoff;
+ real cutoff;
/** Positions of the reference points. */
- gmx_ana_pos_t p;
+ gmx_ana_pos_t p;
/** Neighborhood search data. */
- gmx_ana_nbsearch_t *nb;
-} t_methoddata_distance;
+ gmx::AnalysisNeighborhood nb;
+ /** Neighborhood search for an invididual frame. */
+ gmx::AnalysisNeighborhoodSearch nbsearch;
+};
/** Allocates data for distance-based selection methods. */
static void *
static void *
init_data_common(int npar, gmx_ana_selparam_t *param)
{
- t_methoddata_distance *data;
-
- snew(data, 1);
- data->cutoff = -1;
+ t_methoddata_distance *data = new t_methoddata_distance();
param[0].val.u.r = &data->cutoff;
param[1].val.u.p = &data->p;
return data;
{
GMX_THROW(gmx::InvalidInputError("Distance cutoff should be > 0"));
}
- d->nb = gmx_ana_nbsearch_create(d->cutoff, d->p.nr);
+ d->nb.setCutoff(d->cutoff);
}
/*!
static void
free_data_common(void *data)
{
- t_methoddata_distance *d = (t_methoddata_distance *)data;
-
- if (d->nb)
- {
- gmx_ana_nbsearch_free(d->nb);
- }
- sfree(d);
+ delete static_cast<t_methoddata_distance *>(data);
}
/*!
{
t_methoddata_distance *d = (t_methoddata_distance *)data;
- gmx_ana_nbsearch_pos_init(d->nb, pbc, &d->p);
+ d->nbsearch.reset();
+ d->nbsearch = d->nb.initSearch(pbc, &d->p);
}
/*!
gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void *data)
{
t_methoddata_distance *d = (t_methoddata_distance *)data;
- int b, i;
- real n;
out->nr = pos->g->isize;
- for (b = 0; b < pos->nr; ++b)
+ for (int b = 0; b < pos->nr; ++b)
{
- n = gmx_ana_nbsearch_pos_mindist(d->nb, pos, b);
- for (i = pos->m.mapb.index[b]; i < pos->m.mapb.index[b+1]; ++i)
+ real dist = d->nbsearch.minimumDistance(pos, b);
+ for (int i = pos->m.mapb.index[b]; i < pos->m.mapb.index[b+1]; ++i)
{
- out->u.r[i] = n;
+ out->u.r[i] = dist;
}
}
}
gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void *data)
{
t_methoddata_distance *d = (t_methoddata_distance *)data;
- int b;
out->u.g->isize = 0;
- for (b = 0; b < pos->nr; ++b)
+ for (int b = 0; b < pos->nr; ++b)
{
- if (gmx_ana_nbsearch_pos_is_within(d->nb, pos, b))
+ if (d->nbsearch.isWithin(pos, b))
{
gmx_ana_pos_append(NULL, out->u.g, pos, b, 0);
}
gmx_add_unit_test(SelectionUnitTests selection-test
indexutil.cpp
+ nbsearch.cpp
selectioncollection.cpp
selectionoption.cpp
toputils.cpp)
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Tests selection neighborhood searching.
+ *
+ * \todo
+ * Increase coverage of these tests for different corner cases: other PBC cases
+ * than full 3D, large cutoffs (larger than half the box size), etc.
+ * At least some of these probably don't work correctly.
+ *
+ * \author Teemu Murtola <teemu.murtola@gmail.com>
+ * \ingroup module_selection
+ */
+#include <gtest/gtest.h>
+
+#include <cmath>
+
+#include <limits>
+#include <set>
+
+#include "gromacs/legacyheaders/gmx_random.h"
+#include "gromacs/legacyheaders/pbc.h"
+#include "gromacs/legacyheaders/smalloc.h"
+#include "gromacs/legacyheaders/vec.h"
+
+#include "gromacs/selection/nbsearch.h"
+
+#include "testutils/testasserts.h"
+
+namespace
+{
+
+/********************************************************************
+ * NeighborhoodSearchTestData
+ */
+
+class NeighborhoodSearchTestData
+{
+ public:
+ struct TestPosition
+ {
+ TestPosition() : refMinDist(0.0), refNearestPoint(-1)
+ {
+ clear_rvec(x);
+ }
+ explicit TestPosition(const rvec x)
+ : refMinDist(0.0), refNearestPoint(-1)
+ {
+ copy_rvec(x, this->x);
+ }
+
+ rvec x;
+ real refMinDist;
+ int refNearestPoint;
+ std::set<int> refPairs;
+ };
+ typedef std::vector<TestPosition> TestPositionList;
+
+ NeighborhoodSearchTestData(int seed, real cutoff);
+ ~NeighborhoodSearchTestData();
+
+ void addTestPosition(const rvec x)
+ {
+ testPositions_.push_back(TestPosition(x));
+ }
+ void generateRandomPosition(rvec x);
+ void generateRandomRefPositions(int count);
+ void generateRandomTestPositions(int count);
+ void computeReferences(t_pbc *pbc);
+
+ gmx_rng_t rng_;
+ real cutoff_;
+ matrix box_;
+ t_pbc pbc_;
+ int refPosCount_;
+ rvec *refPos_;
+ TestPositionList testPositions_;
+};
+
+NeighborhoodSearchTestData::NeighborhoodSearchTestData(int seed, real cutoff)
+ : rng_(NULL), cutoff_(cutoff), refPosCount_(0), refPos_(NULL)
+{
+ // TODO: Handle errors.
+ rng_ = gmx_rng_init(seed);
+ clear_mat(box_);
+ set_pbc(&pbc_, epbcNONE, box_);
+}
+
+NeighborhoodSearchTestData::~NeighborhoodSearchTestData()
+{
+ if (rng_ != NULL)
+ {
+ gmx_rng_destroy(rng_);
+ }
+ sfree(refPos_);
+}
+
+void NeighborhoodSearchTestData::generateRandomPosition(rvec x)
+{
+ rvec fx;
+ fx[XX] = gmx_rng_uniform_real(rng_);
+ fx[YY] = gmx_rng_uniform_real(rng_);
+ fx[ZZ] = gmx_rng_uniform_real(rng_);
+ mvmul(box_, fx, x);
+ // Add a small displacement to allow positions outside the box
+ x[XX] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
+ x[YY] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
+ x[ZZ] += 0.2 * gmx_rng_uniform_real(rng_) - 0.1;
+}
+
+void NeighborhoodSearchTestData::generateRandomRefPositions(int count)
+{
+ refPosCount_ = count;
+ snew(refPos_, refPosCount_);
+ for (int i = 0; i < refPosCount_; ++i)
+ {
+ generateRandomPosition(refPos_[i]);
+ }
+}
+
+void NeighborhoodSearchTestData::generateRandomTestPositions(int count)
+{
+ testPositions_.reserve(count);
+ for (int i = 0; i < count; ++i)
+ {
+ rvec x;
+ generateRandomPosition(x);
+ addTestPosition(x);
+ }
+}
+
+void NeighborhoodSearchTestData::computeReferences(t_pbc *pbc)
+{
+ real cutoff = cutoff_;
+ if (cutoff <= 0)
+ {
+ cutoff = std::numeric_limits<real>::max();
+ }
+ TestPositionList::iterator i;
+ for (i = testPositions_.begin(); i != testPositions_.end(); ++i)
+ {
+ i->refMinDist = cutoff;
+ i->refNearestPoint = -1;
+ i->refPairs.clear();
+ for (int j = 0; j < refPosCount_; ++j)
+ {
+ rvec dx;
+ if (pbc != NULL)
+ {
+ pbc_dx(pbc, i->x, refPos_[j], dx);
+ }
+ else
+ {
+ rvec_sub(i->x, refPos_[j], dx);
+ }
+ const real dist = norm(dx);
+ if (dist < i->refMinDist)
+ {
+ i->refMinDist = dist;
+ i->refNearestPoint = j;
+ }
+ if (dist <= cutoff)
+ {
+ i->refPairs.insert(j);
+ }
+ }
+ }
+}
+
+/********************************************************************
+ * NeighborhoodSearchTest
+ */
+
+class NeighborhoodSearchTest : public ::testing::Test
+{
+ public:
+ void testIsWithin(gmx::AnalysisNeighborhoodSearch *search,
+ const NeighborhoodSearchTestData &data);
+ void testMinimumDistance(gmx::AnalysisNeighborhoodSearch *search,
+ const NeighborhoodSearchTestData &data);
+ void testNearestPoint(gmx::AnalysisNeighborhoodSearch *search,
+ const NeighborhoodSearchTestData &data);
+ void testPairSearch(gmx::AnalysisNeighborhoodSearch *search,
+ const NeighborhoodSearchTestData &data);
+
+ gmx::AnalysisNeighborhood nb_;
+};
+
+void NeighborhoodSearchTest::testIsWithin(
+ gmx::AnalysisNeighborhoodSearch *search,
+ const NeighborhoodSearchTestData &data)
+{
+ NeighborhoodSearchTestData::TestPositionList::const_iterator i;
+ for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
+ {
+ const bool bWithin = (i->refMinDist <= data.cutoff_);
+ EXPECT_EQ(bWithin, search->isWithin(i->x))
+ << "Distance is " << i->refMinDist;
+ }
+}
+
+void NeighborhoodSearchTest::testMinimumDistance(
+ gmx::AnalysisNeighborhoodSearch *search,
+ const NeighborhoodSearchTestData &data)
+{
+ NeighborhoodSearchTestData::TestPositionList::const_iterator i;
+ for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
+ {
+ const real refDist = i->refMinDist;
+ EXPECT_NEAR_REL(refDist, search->minimumDistance(i->x), 20*GMX_REAL_EPS);
+ }
+}
+
+void NeighborhoodSearchTest::testNearestPoint(
+ gmx::AnalysisNeighborhoodSearch *search,
+ const NeighborhoodSearchTestData &data)
+{
+ NeighborhoodSearchTestData::TestPositionList::const_iterator i;
+ for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
+ {
+ const gmx::AnalysisNeighborhoodPair pair = search->nearestPoint(i->x);
+ if (pair.isValid())
+ {
+ EXPECT_EQ(i->refNearestPoint, pair.refIndex());
+ EXPECT_EQ(0, pair.testIndex());
+ }
+ else
+ {
+ EXPECT_EQ(i->refNearestPoint, -1);
+ }
+ }
+}
+
+void NeighborhoodSearchTest::testPairSearch(
+ gmx::AnalysisNeighborhoodSearch *search,
+ const NeighborhoodSearchTestData &data)
+{
+ NeighborhoodSearchTestData::TestPositionList::const_iterator i;
+ for (i = data.testPositions_.begin(); i != data.testPositions_.end(); ++i)
+ {
+ std::set<int> checkSet = i->refPairs;
+ gmx::AnalysisNeighborhoodPairSearch pairSearch =
+ search->startPairSearch(i->x);
+ gmx::AnalysisNeighborhoodPair pair;
+ while (pairSearch.findNextPair(&pair))
+ {
+ EXPECT_EQ(0, pair.testIndex());
+ if (checkSet.erase(pair.refIndex()) == 0)
+ {
+ // TODO: Check whether the same pair was returned more than
+ // once and give a better error message if so.
+ ADD_FAILURE()
+ << "Expected: Position " << pair.refIndex()
+ << " is within cutoff.\n"
+ << " Actual: It is not.";
+ }
+ }
+ EXPECT_TRUE(checkSet.empty()) << "Some positions were not returned by the pair search.";
+ }
+}
+
+/********************************************************************
+ * Test data generation
+ */
+
+class RandomBoxFullPBCData
+{
+ public:
+ static const NeighborhoodSearchTestData &get()
+ {
+ static RandomBoxFullPBCData singleton;
+ return singleton.data_;
+ }
+
+ RandomBoxFullPBCData() : data_(12345, 1.0)
+ {
+ data_.box_[XX][XX] = 10.0;
+ data_.box_[YY][YY] = 5.0;
+ data_.box_[ZZ][ZZ] = 7.0;
+ // TODO: Consider whether manually picking some positions would give better
+ // test coverage.
+ data_.generateRandomRefPositions(1000);
+ data_.generateRandomTestPositions(100);
+ set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
+ data_.computeReferences(&data_.pbc_);
+ }
+
+ private:
+ NeighborhoodSearchTestData data_;
+};
+
+class RandomTriclinicFullPBCData
+{
+ public:
+ static const NeighborhoodSearchTestData &get()
+ {
+ static RandomTriclinicFullPBCData singleton;
+ return singleton.data_;
+ }
+
+ RandomTriclinicFullPBCData() : data_(12345, 1.0)
+ {
+ data_.box_[XX][XX] = 5.0;
+ data_.box_[YY][XX] = 2.5;
+ data_.box_[YY][YY] = 2.5*sqrt(3.0);
+ data_.box_[ZZ][XX] = 2.5;
+ data_.box_[ZZ][YY] = 2.5*sqrt(1.0/3.0);
+ data_.box_[ZZ][ZZ] = 5.0*sqrt(2.0/3.0);
+ // TODO: Consider whether manually picking some positions would give better
+ // test coverage.
+ data_.generateRandomRefPositions(1000);
+ data_.generateRandomTestPositions(100);
+ set_pbc(&data_.pbc_, epbcXYZ, data_.box_);
+ data_.computeReferences(&data_.pbc_);
+ }
+
+ private:
+ NeighborhoodSearchTestData data_;
+};
+
+/********************************************************************
+ * Actual tests
+ */
+
+TEST_F(NeighborhoodSearchTest, SimpleSearch)
+{
+ const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
+
+ nb_.setCutoff(data.cutoff_);
+ nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Simple);
+ gmx::AnalysisNeighborhoodSearch search =
+ nb_.initSearch(&data.pbc_, data.refPosCount_, data.refPos_);
+ ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
+
+ testIsWithin(&search, data);
+ testMinimumDistance(&search, data);
+ testNearestPoint(&search, data);
+ testPairSearch(&search, data);
+}
+
+TEST_F(NeighborhoodSearchTest, GridSearchBox)
+{
+ const NeighborhoodSearchTestData &data = RandomBoxFullPBCData::get();
+
+ nb_.setCutoff(data.cutoff_);
+ nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
+ gmx::AnalysisNeighborhoodSearch search =
+ nb_.initSearch(&data.pbc_, data.refPosCount_, data.refPos_);
+ ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
+
+ testIsWithin(&search, data);
+ testMinimumDistance(&search, data);
+ testNearestPoint(&search, data);
+ testPairSearch(&search, data);
+}
+
+TEST_F(NeighborhoodSearchTest, GridSearchTriclinic)
+{
+ const NeighborhoodSearchTestData &data = RandomTriclinicFullPBCData::get();
+
+ nb_.setCutoff(data.cutoff_);
+ nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
+ gmx::AnalysisNeighborhoodSearch search =
+ nb_.initSearch(&data.pbc_, data.refPosCount_, data.refPos_);
+ ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
+
+ testPairSearch(&search, data);
+}
+
+} // namespace
// Tell the analysis framework that this component exists
registerAnalysisDataset(&data_, "freevolume");
rng_ = NULL;
- nbsearch_ = NULL;
nmol_ = 0;
mtot_ = 0;
cutoff_ = 0;
{
gmx_rng_destroy(rng_);
}
- if (NULL != nbsearch_)
- {
- gmx_ana_nbsearch_free(nbsearch_);
- }
}
rng_ = gmx_rng_init(seed_);
// Initiate the neighborsearching code
- nbsearch_ = gmx_ana_nbsearch_create(cutoff_, atoms->nr);
+ nb_.setCutoff(cutoff_);
}
void
int Ninsert = static_cast<int>(ninsert_*V);
// Use neighborsearching tools!
- gmx_ana_nbsearch_pos_init(nbsearch_, pbc, sel.positions());
+ AnalysisNeighborhoodSearch nbsearch = nb_.initSearch(pbc, sel.positions());
// Then loop over insertions
int NinsTot = 0;
mvmul(fr.box, rand, ins);
// Find the first reference position within the cutoff.
- bool bOverlap = false;
- int jp = NOTSET;
- if (gmx_ana_nbsearch_first_within(nbsearch_, ins, &jp))
+ bool bOverlap = false;
+ AnalysisNeighborhoodPair pair;
+ AnalysisNeighborhoodPairSearch pairSearch = nbsearch.startPairSearch(ins);
+ while (!bOverlap && pairSearch.findNextPair(&pair))
{
- do
- {
- // Compute distance vector to first atom in the neighborlist
- pbc_dx(pbc, ins, sel.position(jp).x(), dx);
+ int jp = pair.refIndex();
+ // Compute distance vector to first atom in the neighborlist
+ pbc_dx(pbc, ins, sel.position(jp).x(), dx);
- // See whether the distance is smaller than allowed
- bOverlap = (norm(dx) <
- probeRadius_+vdw_radius_[sel.position(jp).refId()]);
+ // See whether the distance is smaller than allowed
+ bOverlap = (norm(dx) <
+ probeRadius_+vdw_radius_[sel.position(jp).refId()]);
- // Finds the next reference position within the cutoff.
- }
- while (!bOverlap &&
- (gmx_ana_nbsearch_next_within(nbsearch_, &jp)));
}
if (!bOverlap)
double probeRadius_;
gmx_rng_t rng_;
int seed_, ninsert_;
- gmx_ana_nbsearch_t *nbsearch_;
+ AnalysisNeighborhood nb_;
//! The van der Waals radius per atom
std::vector<double> vdw_radius_;
*/
/*! \libinternal \file
* \brief
- * Improved exception assertions for unit tests.
+ * Extra assertions for unit tests.
*
* This file provides assert macros to replace (ASSERT|EXPECT)(_NO)?_THROW
* from Google Test. They behave otherwise the same as the Google Test ones,
* but also print details of any unexpected exceptions. This makes it much
* easier to see at one glance what went wrong.
*
+ * This file also provides extra floating-point assertions.
+ *
* \if internal
* \todo
* The implementation is somewhat ugly, and accesses some Google Test
#include <gtest/gtest.h>
+#include "gromacs/legacyheaders/maths.h"
+
#include "gromacs/utility/exceptions.h"
+namespace gmx
+{
+namespace test
+{
+
/*! \cond internal */
/*! \internal
* \brief
#define ASSERT_NO_THROW_GMX(statement) \
GMX_TEST_NO_THROW_(statement, GTEST_FATAL_FAILURE_)
+/*! \cond internal */
+/*! \internal \brief
+ * Assertion predicate formatter for comparing two floating-point values.
+ */
+static ::testing::AssertionResult assertWithinRelativeTolerance(
+ const char *expr1, const char *expr2, const char * /*exprTolerance*/,
+ real val1, real val2, real relativeTolerance)
+{
+ if (gmx_within_tol(val1, val2, relativeTolerance))
+ {
+ return ::testing::AssertionSuccess();
+ }
+ return ::testing::AssertionFailure()
+ << "Value of: " << expr2 << " not within tolerance of " << relativeTolerance << "\n"
+ << " Actual: " << val2 << "\n"
+ << "Expected: " << expr1 << "\n"
+ << "Which is: " << val1;
+}
+//! \endcond
+
+/*! \brief
+ * Asserts that two floating-point values are within the given relative error.
+ *
+ * This assert works as EXPECT_NEAR from Google Test, except that it uses a
+ * relative instead of a absolute tolerance.
+ * See gmx_within_tol() for definition of the tolerance.
+ */
+#define EXPECT_NEAR_REL(val1, val2, rel_error) \
+ EXPECT_PRED_FORMAT3(::gmx::test::assertWithinRelativeTolerance, val1, val2, rel_error)
+/*! \brief
+ * Asserts that two floating-point values are within the given relative error.
+ *
+ * This assert works as ASSERT_NEAR from Google Test, except that it uses a
+ * relative instead of a absolute tolerance.
+ * See gmx_within_tol() for definition of the tolerance.
+ */
+#define ASSERT_NEAR_REL(val1, val2, rel_error) \
+ ASSERT_PRED_FORMAT3(::gmx::test::assertWithinRelativeTolerance, val1, val2, rel_error)
+
+} // namespace test
+} // namespace gmx
+
#endif
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff