AnalysisDataHandle dh = pdata->dataHandle(data_);
const Selection &refsel = pdata->parallelSelection(refsel_);
- AnalysisNeighborhoodSearch nbsearch =
- nb_.initSearch(pbc, refsel.positions());
+ AnalysisNeighborhoodSearch nbsearch = nb_.initSearch(pbc, refsel);
dh.startFrame(frnr, fr.time);
for (size_t g = 0; g < sel_.size(); ++g)
{
#include "gromacs/legacyheaders/thread_mpi/mutex.h"
#include "gromacs/selection/position.h"
+#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/gmxassert.h"
namespace gmx
explicit AnalysisNeighborhoodSearchImpl(real cutoff);
~AnalysisNeighborhoodSearchImpl();
- void init(AnalysisNeighborhood::SearchMode mode,
- const t_pbc *pbc, int n, const rvec x[], const int *refid);
+ /*! \brief
+ * Initializes the search with a given box and reference positions.
+ *
+ * \param[in] mode Search mode to use.
+ * \param[in] pbc PBC information.
+ * \param[in] positions Set of reference positions.
+ */
+ void init(AnalysisNeighborhood::SearchMode mode,
+ const t_pbc *pbc,
+ const AnalysisNeighborhoodPositions &positions);
PairSearchImplPointer getPairSearch();
real cutoffSquared() const { return cutoff2_; }
{
public:
explicit AnalysisNeighborhoodPairSearchImpl(const AnalysisNeighborhoodSearchImpl &search)
- : search_(search), testIndex_(0)
+ : search_(search)
{
clear_rvec(xtest_);
clear_ivec(testcell_);
- reset();
+ reset(-1);
}
- //! Clears the loop indices.
- void reset();
//! Initializes a search to find reference positions neighboring \p x.
- void startSearch(const rvec x, int testIndex);
+ void startSearch(const AnalysisNeighborhoodPositions &positions);
//! Searches for the next neighbor.
template <class Action>
bool searchNext(Action action);
//! Initializes a pair representing the pair found by searchNext().
void initFoundPair(AnalysisNeighborhoodPair *pair) const;
+ //! Advances to the next test position, skipping any remaining pairs.
+ void nextTestPosition();
private:
+ //! Clears the loop indices.
+ void reset(int testIndex);
//! Checks whether a reference positiong should be excluded.
bool isExcluded(int j);
+ //! Parent search object.
const AnalysisNeighborhoodSearchImpl &search_;
+ //! Reference to the test positions.
+ ConstArrayRef<rvec> testPositions_;
+ //! Index of the currently active test position in \p testPositions_.
int testIndex_;
//! Stores test position during a pair loop.
rvec xtest_;
}
void AnalysisNeighborhoodSearchImpl::init(
- AnalysisNeighborhood::SearchMode mode,
- const t_pbc *pbc, int n, const rvec x[], const int *refid)
+ AnalysisNeighborhood::SearchMode mode,
+ const t_pbc *pbc,
+ const AnalysisNeighborhoodPositions &positions)
{
+ GMX_RELEASE_ASSERT(positions.index_ == -1,
+ "Individual indexed positions not supported as reference");
pbc_ = const_cast<t_pbc *>(pbc);
- nref_ = n;
+ nref_ = positions.count_;
// TODO: Consider whether it would be possible to support grid searching in
// more cases.
if (mode == AnalysisNeighborhood::eSearchMode_Simple
}
if (bGrid_)
{
- if (xref_nalloc_ < n)
+ if (xref_nalloc_ < nref_)
{
- srenew(xref_alloc_, n);
- xref_nalloc_ = n;
+ srenew(xref_alloc_, nref_);
+ xref_nalloc_ = nref_;
}
xref_ = xref_alloc_;
- for (int i = 0; i < n; ++i)
+ for (int i = 0; i < nref_; ++i)
{
- copy_rvec(x[i], xref_alloc_[i]);
+ copy_rvec(positions.x_[i], xref_alloc_[i]);
}
- put_atoms_in_triclinic_unitcell(ecenterTRIC, pbc_->box, n, xref_alloc_);
- for (int i = 0; i < n; ++i)
+ put_atoms_in_triclinic_unitcell(ecenterTRIC, pbc_->box,
+ nref_, xref_alloc_);
+ for (int i = 0; i < nref_; ++i)
{
ivec refcell;
}
else
{
- xref_ = x;
+ xref_ = positions.x_;
}
- refid_ = refid;
+ // TODO: Once exclusions are supported, this may need to be initialized.
+ refid_ = NULL;
}
#if 0
* AnalysisNeighborhoodPairSearchImpl
*/
-void AnalysisNeighborhoodPairSearchImpl::reset()
+void AnalysisNeighborhoodPairSearchImpl::reset(int testIndex)
{
- previ_ = -1;
- exclind_ = 0;
- prevnbi_ = 0;
- prevcai_ = -1;
+ testIndex_ = testIndex;
+ if (testIndex_ >= 0 && testIndex_ < static_cast<int>(testPositions_.size()))
+ {
+ copy_rvec(testPositions_[testIndex_], xtest_);
+ if (search_.bGrid_)
+ {
+ put_atoms_in_triclinic_unitcell(ecenterTRIC, search_.pbc_->box,
+ 1, &xtest_);
+ search_.mapPointToGridCell(xtest_, testcell_);
+ }
+ }
+ previ_ = -1;
+ exclind_ = 0;
+ prevnbi_ = 0;
+ prevcai_ = -1;
+}
+
+void AnalysisNeighborhoodPairSearchImpl::nextTestPosition()
+{
+ if (testIndex_ < static_cast<int>(testPositions_.size()))
+ {
+ ++testIndex_;
+ reset(testIndex_);
+ }
}
bool AnalysisNeighborhoodPairSearchImpl::isExcluded(int j)
return false;
}
-void AnalysisNeighborhoodPairSearchImpl::startSearch(const rvec x, int testIndex)
+void AnalysisNeighborhoodPairSearchImpl::startSearch(
+ const AnalysisNeighborhoodPositions &positions)
{
- testIndex_ = testIndex;
- copy_rvec(x, xtest_);
- if (search_.bGrid_)
+ if (positions.index_ < 0)
{
- put_atoms_in_triclinic_unitcell(ecenterTRIC, search_.pbc_->box,
- 1, &xtest_);
- search_.mapPointToGridCell(xtest_, testcell_);
+ testPositions_ = ConstArrayRef<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_ = ConstArrayRef<rvec>(positions.x_, positions.index_ + 1);
+ reset(positions.index_);
}
- reset();
}
template <class Action>
bool AnalysisNeighborhoodPairSearchImpl::searchNext(Action action)
{
- if (search_.bGrid_)
+ while (testIndex_ < static_cast<int>(testPositions_.size()))
{
- int nbi = prevnbi_;
- int cai = prevcai_ + 1;
-
- for (; nbi < search_.ngridnb_; ++nbi)
+ if (search_.bGrid_)
{
- ivec cell;
+ int nbi = prevnbi_;
+ int cai = prevcai_ + 1;
+
+ for (; nbi < search_.ngridnb_; ++nbi)
+ {
+ ivec cell;
- ivec_add(testcell_, search_.gnboffs_[nbi], cell);
- cell[XX] = (cell[XX] + search_.ncelldim_[XX]) % search_.ncelldim_[XX];
- cell[YY] = (cell[YY] + search_.ncelldim_[YY]) % search_.ncelldim_[YY];
- cell[ZZ] = (cell[ZZ] + search_.ncelldim_[ZZ]) % search_.ncelldim_[ZZ];
+ ivec_add(testcell_, search_.gnboffs_[nbi], cell);
+ cell[XX] = (cell[XX] + search_.ncelldim_[XX]) % search_.ncelldim_[XX];
+ cell[YY] = (cell[YY] + search_.ncelldim_[YY]) % search_.ncelldim_[YY];
+ cell[ZZ] = (cell[ZZ] + search_.ncelldim_[ZZ]) % search_.ncelldim_[ZZ];
- const int ci = search_.getGridCellIndex(cell);
- const int cellSize = static_cast<int>(search_.cells_[ci].size());
- /* TODO: Calculate the required PBC shift outside the inner loop */
- for (; cai < cellSize; ++cai)
+ const int ci = search_.getGridCellIndex(cell);
+ const int cellSize = static_cast<int>(search_.cells_[ci].size());
+ /* TODO: Calculate the required PBC shift outside the inner loop */
+ for (; cai < cellSize; ++cai)
+ {
+ const int i = search_.cells_[ci][cai];
+ if (isExcluded(i))
+ {
+ continue;
+ }
+ rvec dx;
+ pbc_dx_aiuc(search_.pbc_, xtest_, search_.xref_[i], dx);
+ const real r2 = norm2(dx);
+ if (r2 <= search_.cutoff2_)
+ {
+ if (action(i, r2))
+ {
+ prevnbi_ = nbi;
+ prevcai_ = cai;
+ previ_ = i;
+ return true;
+ }
+ }
+ }
+ exclind_ = 0;
+ cai = 0;
+ }
+ }
+ else
+ {
+ for (int i = previ_ + 1; i < search_.nref_; ++i)
{
- const int i = search_.cells_[ci][cai];
if (isExcluded(i))
{
continue;
}
- rvec dx;
- pbc_dx_aiuc(search_.pbc_, xtest_, search_.xref_[i], dx);
+ rvec dx;
+ if (search_.pbc_)
+ {
+ pbc_dx(search_.pbc_, xtest_, search_.xref_[i], dx);
+ }
+ else
+ {
+ rvec_sub(xtest_, search_.xref_[i], dx);
+ }
const real r2 = norm2(dx);
if (r2 <= search_.cutoff2_)
{
if (action(i, r2))
{
- prevnbi_ = nbi;
- prevcai_ = cai;
- previ_ = i;
+ previ_ = i;
return true;
}
}
}
- exclind_ = 0;
- cai = 0;
}
+ nextTestPosition();
}
- else
- {
- for (int i = previ_ + 1; i < search_.nref_; ++i)
- {
- if (isExcluded(i))
- {
- continue;
- }
- rvec dx;
- if (search_.pbc_)
- {
- pbc_dx(search_.pbc_, xtest_, search_.xref_[i], dx);
- }
- else
- {
- rvec_sub(xtest_, search_.xref_[i], dx);
- }
- const real r2 = norm2(dx);
- if (r2 <= search_.cutoff2_)
- {
- if (action(i, r2))
- {
- previ_ = i;
- return true;
- }
- }
- }
- }
- reset();
return false;
}
}
AnalysisNeighborhoodSearch
-AnalysisNeighborhood::initSearch(const t_pbc *pbc, int n, const rvec x[])
+AnalysisNeighborhood::initSearch(const t_pbc *pbc,
+ const AnalysisNeighborhoodPositions &positions)
{
Impl::SearchImplPointer search(impl_->getSearch());
- search->init(mode(), pbc, n, x, NULL);
- return AnalysisNeighborhoodSearch(search);
-}
-
-AnalysisNeighborhoodSearch
-AnalysisNeighborhood::initSearch(const t_pbc *pbc, const gmx_ana_pos_t *p)
-{
- Impl::SearchImplPointer search(impl_->getSearch());
- search->init(mode(), pbc, p->nr, p->x, p->m.refid);
+ search->init(mode(), pbc, positions);
return AnalysisNeighborhoodSearch(search);
}
: AnalysisNeighborhood::eSearchMode_Simple);
}
-bool AnalysisNeighborhoodSearch::isWithin(const rvec x) const
+bool AnalysisNeighborhoodSearch::isWithin(
+ const AnalysisNeighborhoodPositions &positions) const
{
GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
- pairSearch.startSearch(x, 0);
+ pairSearch.startSearch(positions);
return pairSearch.searchNext(&withinAction);
}
-bool AnalysisNeighborhoodSearch::isWithin(const gmx_ana_pos_t *p, int i) const
-{
- return isWithin(p->x[i]);
-}
-
-real AnalysisNeighborhoodSearch::minimumDistance(const rvec x) const
+real AnalysisNeighborhoodSearch::minimumDistance(
+ const AnalysisNeighborhoodPositions &positions) const
{
GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
- pairSearch.startSearch(x, 0);
+ pairSearch.startSearch(positions);
real minDist2 = impl_->cutoffSquared();
int closestPoint = -1;
MindistAction action(&closestPoint, &minDist2);
return sqrt(minDist2);
}
-real AnalysisNeighborhoodSearch::minimumDistance(const gmx_ana_pos_t *p, int i) const
-{
- return minimumDistance(p->x[i]);
-}
-
AnalysisNeighborhoodPair
-AnalysisNeighborhoodSearch::nearestPoint(const rvec x) const
+AnalysisNeighborhoodSearch::nearestPoint(
+ const AnalysisNeighborhoodPositions &positions) const
{
GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
- pairSearch.startSearch(x, 0);
+ pairSearch.startSearch(positions);
real minDist2 = impl_->cutoffSquared();
int closestPoint = -1;
MindistAction action(&closestPoint, &minDist2);
return AnalysisNeighborhoodPair(closestPoint, 0);
}
-AnalysisNeighborhoodPair
-AnalysisNeighborhoodSearch::nearestPoint(const gmx_ana_pos_t *p, int i) const
-{
- GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
- internal::AnalysisNeighborhoodPairSearchImpl pairSearch(*impl_);
- pairSearch.startSearch(p->x[i], 0);
- real minDist2 = impl_->cutoffSquared();
- int closestPoint = -1;
- MindistAction action(&closestPoint, &minDist2);
- (void)pairSearch.searchNext(action);
- return AnalysisNeighborhoodPair(closestPoint, i);
-}
-
AnalysisNeighborhoodPairSearch
-AnalysisNeighborhoodSearch::startPairSearch(const rvec x)
+AnalysisNeighborhoodSearch::startPairSearch(
+ const AnalysisNeighborhoodPositions &positions) const
{
GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
Impl::PairSearchImplPointer pairSearch(impl_->getPairSearch());
- pairSearch->startSearch(x, 0);
- return AnalysisNeighborhoodPairSearch(pairSearch);
-}
-
-AnalysisNeighborhoodPairSearch
-AnalysisNeighborhoodSearch::startPairSearch(const gmx_ana_pos_t *p, int i)
-{
- GMX_RELEASE_ASSERT(impl_, "Accessing an invalid search object");
- Impl::PairSearchImplPointer pairSearch(impl_->getPairSearch());
- pairSearch->startSearch(p->x[i], i);
+ pairSearch->startSearch(positions);
return AnalysisNeighborhoodPairSearch(pairSearch);
}
return bFound;
}
+void AnalysisNeighborhoodPairSearch::skipRemainingPairsForTestPosition()
+{
+ impl_->nextTestPosition();
+}
+
} // namespace gmx
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)
+ {
+ }
+ /*! \brief
+ * Initializes positions from an array of position vectors.
+ */
+ AnalysisNeighborhoodPositions(const rvec x[], int count)
+ : count_(count), index_(-1), x_(x)
+ {
+ }
+
+ /*! \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_;
+
+ //! 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.
*
/*! \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;
/*! \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.
+ * \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.
*/
- bool isWithin(const gmx_ana_pos_t *p, int i) const;
+ bool isWithin(const AnalysisNeighborhoodPositions &positions) const;
/*! \brief
* Calculates the minimum distance from the reference points.
*
- * \param[in] x Test position.
+ * \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 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;
+ 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.
* \throws std::bad_alloc if out of memory.
- *
- * 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.
- * \throws std::bad_alloc if out of memory.
- *
- * 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);
+ AnalysisNeighborhoodPairSearch
+ startPairSearch(const AnalysisNeighborhoodPositions &positions) const;
private:
typedef internal::AnalysisNeighborhoodSearchImpl Impl;
* \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))
{
* \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_;
*/
#include "selection.h"
+#include "nbsearch.h"
#include "position.h"
#include "selelem.h"
#include "selvalue.h"
namespace internal
{
+/********************************************************************
+ * SelectionData
+ */
+
SelectionData::SelectionData(SelectionTreeElement *elem,
const char *selstr)
: name_(elem->name()), selectionText_(selstr),
} // namespace internal
+/********************************************************************
+ * Selection
+ */
+
+Selection::operator AnalysisNeighborhoodPositions() const
+{
+ return AnalysisNeighborhoodPositions(data().rawPositions_.x,
+ data().rawPositions_.nr);
+}
+
void
Selection::printInfo(FILE *fp) const
fprintf(fp, "\n");
}
+
+/********************************************************************
+ * SelectionPosition
+ */
+
+SelectionPosition::operator AnalysisNeighborhoodPositions() const
+{
+ return AnalysisNeighborhoodPositions(sel_->rawPositions_.x,
+ sel_->rawPositions_.nr)
+ .selectSingleFromArray(i_);
+}
+
} // namespace gmx
class SelectionOptionStorage;
class SelectionTreeElement;
+class AnalysisNeighborhoodPositions;
class Selection;
class SelectionPosition;
return ConstArrayRef<int>(data().rawPositions_.m.mapid, posCount());
}
- //! Deprecated method for direct access to position data.
- const gmx_ana_pos_t *positions() const { return &data().rawPositions_; }
-
//! Returns whether the covered fraction can change between frames.
bool isCoveredFractionDynamic() const { return data().isCoveredFractionDynamic(); }
//! Returns the covered fraction for the current frame.
real coveredFraction() const { return data().coveredFraction_; }
+ /*! \brief
+ * Allows passing a selection directly to neighborhood searching.
+ *
+ * When initialized this way, AnalysisNeighborhoodPair objects return
+ * indices that can be used to index the selection positions with
+ * position().
+ *
+ * Works exactly like if AnalysisNeighborhoodPositions had a
+ * constructor taking a Selection object as a parameter.
+ * See AnalysisNeighborhoodPositions for rationale and additional
+ * discussion.
+ */
+ operator AnalysisNeighborhoodPositions() const;
+
/*! \brief
* Initializes information about covered fractions.
*
return sel_->rawPositions_.m.mapid[i_];
}
+ /*! \brief
+ * Allows passing a selection position directly to neighborhood searching.
+ *
+ * When initialized this way, AnalysisNeighborhoodPair objects return
+ * the index that can be used to access this position using
+ * Selection::position().
+ *
+ * Works exactly like if AnalysisNeighborhoodPositions had a
+ * constructor taking a SelectionPosition object as a parameter.
+ * See AnalysisNeighborhoodPositions for rationale and additional
+ * discussion.
+ */
+ operator AnalysisNeighborhoodPositions() const;
+
private:
const internal::SelectionData *sel_;
int i_;
t_methoddata_distance *d = (t_methoddata_distance *)data;
d->nbsearch.reset();
- d->nbsearch = d->nb.initSearch(pbc, &d->p);
+ gmx::AnalysisNeighborhoodPositions pos(d->p.x, d->p.nr);
+ d->nbsearch = d->nb.initSearch(pbc, pos);
}
/*!
out->nr = pos->g->isize;
for (int b = 0; b < pos->nr; ++b)
{
- real dist = d->nbsearch.minimumDistance(pos, b);
+ real dist = d->nbsearch.minimumDistance(pos->x[b]);
for (int i = pos->m.mapb.index[b]; i < pos->m.mapb.index[b+1]; ++i)
{
out->u.r[i] = dist;
out->u.g->isize = 0;
for (int b = 0; b < pos->nr; ++b)
{
- if (d->nbsearch.isWithin(pos, b))
+ if (d->nbsearch.isWithin(pos->x[b]))
{
gmx_ana_pos_append(NULL, out->u.g, pos, b, 0);
}
NeighborhoodSearchTestData(int seed, real cutoff);
~NeighborhoodSearchTestData();
+ gmx::AnalysisNeighborhoodPositions refPositions() const
+ {
+ return gmx::AnalysisNeighborhoodPositions(refPos_, refPosCount_);
+ }
+ gmx::AnalysisNeighborhoodPositions testPositions() const
+ {
+ if (testPos_ == NULL)
+ {
+ snew(testPos_, testPositions_.size());
+ for (size_t i = 0; i < testPositions_.size(); ++i)
+ {
+ copy_rvec(testPositions_[i].x, testPos_[i]);
+ }
+ }
+ return gmx::AnalysisNeighborhoodPositions(testPos_,
+ testPositions_.size());
+ }
+ gmx::AnalysisNeighborhoodPositions testPosition(int index) const
+ {
+ return testPositions().selectSingleFromArray(index);
+ }
+
void addTestPosition(const rvec x)
{
+ GMX_RELEASE_ASSERT(testPos_ == NULL,
+ "Cannot add positions after testPositions() call");
testPositions_.push_back(TestPosition(x));
}
void generateRandomPosition(rvec x);
int refPosCount_;
rvec *refPos_;
TestPositionList testPositions_;
+
+ private:
+ mutable rvec *testPos_;
};
NeighborhoodSearchTestData::NeighborhoodSearchTestData(int seed, real cutoff)
- : rng_(NULL), cutoff_(cutoff), refPosCount_(0), refPos_(NULL)
+ : rng_(NULL), cutoff_(cutoff), refPosCount_(0), refPos_(NULL), testPos_(NULL)
{
// TODO: Handle errors.
rng_ = gmx_rng_init(seed);
gmx_rng_destroy(rng_);
}
sfree(refPos_);
+ sfree(testPos_);
}
void NeighborhoodSearchTestData::generateRandomPosition(rvec x)
nb_.setCutoff(data.cutoff_);
nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Simple);
gmx::AnalysisNeighborhoodSearch search =
- nb_.initSearch(&data.pbc_, data.refPosCount_, data.refPos_);
+ nb_.initSearch(&data.pbc_, data.refPositions());
ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Simple, search.mode());
testIsWithin(&search, data);
nb_.setCutoff(data.cutoff_);
nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
gmx::AnalysisNeighborhoodSearch search =
- nb_.initSearch(&data.pbc_, data.refPosCount_, data.refPos_);
+ nb_.initSearch(&data.pbc_, data.refPositions());
ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
testIsWithin(&search, data);
nb_.setCutoff(data.cutoff_);
nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
gmx::AnalysisNeighborhoodSearch search =
- nb_.initSearch(&data.pbc_, data.refPosCount_, data.refPos_);
+ nb_.initSearch(&data.pbc_, data.refPositions());
ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
testPairSearch(&search, data);
nb_.setCutoff(data.cutoff_);
nb_.setMode(gmx::AnalysisNeighborhood::eSearchMode_Grid);
gmx::AnalysisNeighborhoodSearch search =
- nb_.initSearch(&data.pbc_, data.refPosCount_, data.refPos_);
+ nb_.initSearch(&data.pbc_, data.refPositions());
// Currently, grid searching not supported with 2D PBC.
//ASSERT_EQ(gmx::AnalysisNeighborhood::eSearchMode_Grid, search.mode());
nb_.setCutoff(data.cutoff_);
gmx::AnalysisNeighborhoodSearch search1 =
- nb_.initSearch(&data.pbc_, data.refPosCount_, data.refPos_);
+ nb_.initSearch(&data.pbc_, data.refPositions());
gmx::AnalysisNeighborhoodSearch search2 =
- nb_.initSearch(&data.pbc_, data.refPosCount_, data.refPos_);
+ nb_.initSearch(&data.pbc_, data.refPositions());
gmx::AnalysisNeighborhoodPairSearch pairSearch1 =
- search1.startPairSearch(data.testPositions_[0].x);
+ search1.startPairSearch(data.testPosition(0));
gmx::AnalysisNeighborhoodPairSearch pairSearch2 =
- search1.startPairSearch(data.testPositions_[1].x);
+ search1.startPairSearch(data.testPosition(1));
testPairSearch(&search2, data);
EXPECT_TRUE(data.testPositions_[0].refPairs.count(pair.refIndex()) == 1);
pairSearch2.findNextPair(&pair);
- EXPECT_EQ(0, pair.testIndex());
+ EXPECT_EQ(1, pair.testIndex());
EXPECT_TRUE(data.testPositions_[1].refPairs.count(pair.refIndex()) == 1);
}
+TEST_F(NeighborhoodSearchTest, HandlesSkippingPairs)
+{
+ const NeighborhoodSearchTestData &data = TrivialTestData::get();
+
+ nb_.setCutoff(data.cutoff_);
+ gmx::AnalysisNeighborhoodSearch search =
+ nb_.initSearch(&data.pbc_, data.refPositions());
+ gmx::AnalysisNeighborhoodPairSearch pairSearch =
+ search.startPairSearch(data.testPositions());
+ gmx::AnalysisNeighborhoodPair pair;
+ // TODO: This test needs to be adjusted if the grid search gets optimized
+ // to loop over the test positions in cell order (first, the ordering
+ // assumption here breaks, and second, it then needs to be tested
+ // separately for simple and grid searches).
+ int currentIndex = 0;
+ while (pairSearch.findNextPair(&pair))
+ {
+ while (currentIndex < pair.testIndex())
+ {
+ ++currentIndex;
+ }
+ EXPECT_EQ(currentIndex, pair.testIndex());
+ EXPECT_TRUE(data.testPositions_[currentIndex].refPairs.count(pair.refIndex()) == 1);
+ pairSearch.skipRemainingPairsForTestPosition();
+ ++currentIndex;
+ }
+}
+
} // namespace
int Ninsert = static_cast<int>(ninsert_*V);
// Use neighborsearching tools!
- AnalysisNeighborhoodSearch nbsearch = nb_.initSearch(pbc, sel.positions());
+ AnalysisNeighborhoodSearch nbsearch = nb_.initSearch(pbc, sel);
// Then loop over insertions
int NinsTot = 0;
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff