Index: kernel/include/IMP/SConscript
===================================================================
--- kernel/include/IMP/SConscript (revision 345)
+++ kernel/include/IMP/SConscript (working copy)
@@ -6,7 +6,7 @@
'ModelData.h', 'RigidBody.h', 'log.h', 'DerivativeAccumulator.h',
'Key.h', 'AttributeTable.h', 'utility.h', 'Restraint.h', 'Optimizer.h',
'DecoratorBase.h', 'Object.h', 'Vector3D.h', 'ScoreFuncParams.h',
- 'UnaryFunction.h', 'PairScore.h']
+ 'UnaryFunction.h', 'PairScore.h', 'SingletonScore.h', 'Grid3D.h']
# Install the include files:
includedir = os.path.join(env['includedir'], 'IMP')
@@ -19,4 +19,5 @@
SConscript('decorators/SConscript')
SConscript('unary_functions/SConscript')
SConscript('pair_scores/SConscript')
+SConscript('singleton_scores/SConscript')
SConscript('score_states/SConscript')
Index: kernel/include/IMP/Grid3D.h
===================================================================
--- kernel/include/IMP/Grid3D.h (revision 0)
+++ kernel/include/IMP/Grid3D.h (revision 0)
@@ -0,0 +1,333 @@
+/**
+ * \file Grid3D.h \brief A class to represent a voxel grid.
+ *
+ * Copyright 2007-8 Sali Lab. All rights reserved.
+ *
+ */
+#ifndef IMP_GRID_3D_H
+#define IMP_GRID_3D_H
+#include "base_types.h"
+#include "Vector3D.h"
+namespace IMP {
+template <class V>
+class Grid3D;
+
+namespace internal {
+
+template <class GI>
+class GridIndexIterator;
+
+class VirtualGridIndex {
+ typedef VirtualGridIndex This;
+ int d_[3];
+ bool is_default() const {
+ return d_[0]==std::numeric_limits<int>::max();
+ }
+ public:
+ VirtualGridIndex(int x, int y, int z) {
+ d_[0]=x;
+ d_[1]=y;
+ d_[2]=z;
+ }
+ VirtualGridIndex() {
+ d_[0]=d_[1]=d_[2]=std::numeric_limits<int>::max();
+ }
+ int operator[](unsigned int i) const {
+ IMP_assert(i <3, "Bad i");
+ return d_[i];
+ }
+ IMP_COMPARISONS_3(d_[0], d_[1], d_[2]);
+ void show(std::ostream &out=std::cout) const {
+ out << "Cell(" << d_[0] << ", " << d_[1] << ", " << d_[2] << ")";
+ }
+};
+
+IMP_OUTPUT_OPERATOR(VirtualGridIndex);
+
+class GridIndex: public VirtualGridIndex {
+public:
+ GridIndex(): VirtualGridIndex() {
+ }
+protected:
+ template <class V>
+ friend class IMP::Grid3D;
+ template <class G>
+ friend class GridIndexIterator;
+ GridIndex(int x, int y, int z): VirtualGridIndex(x,y,z) {
+ IMP_check(x>=0 && y>=0 && z>=0, "Bad indexes in grid index",
+ IndexException("Bad indexes in GridIndex"));
+ }
+};
+
+template <class GI>
+class GridIndexIterator {
+ VirtualGridIndex lb_;
+ VirtualGridIndex ub_;
+ GI cur_;
+ typedef GridIndexIterator This;
+ bool is_default() const {
+ return false;
+ }
+public:
+ typedef const GI& reference_type;
+ typedef const GI* pointer_type;
+ typedef GI value_type;
+ typedef std::forward_iterator_tag iterator_category;
+ typedef unsigned int difference_type;
+ GridIndexIterator(VirtualGridIndex lb,
+ VirtualGridIndex ub): lb_(lb),
+ ub_(ub), cur_(lb[0], lb[1], lb[2]) {
+ IMP_assert(lb_ <= ub_, "Invalid range in GridIndexIterator");
+ }
+ GridIndexIterator(){}
+
+ IMP_COMPARISONS_1(cur_);
+
+ This operator++() {
+ IMP_assert(cur_ >= lb_, "cur out of range");
+ IMP_assert(cur_ < ub_, "cur out of range");
+ int r[3];
+ unsigned int carry=1;
+ for (int i=2; i>=0; --i) {
+ r[i]= cur_[i]+carry;
+ if ( r[i] == ub_[i]) {
+ r[i]= lb_[i];
+ carry=1;
+ } else {
+ carry=0;
+ }
+ }
+ if (carry==1) {
+ cur_= GI();
+ } else {
+ GI nc= GI(r[0], r[1], r[2]);
+ IMP_assert(nc > cur_, "Nonfunctional increment");
+ IMP_assert(nc >= lb_, "Problems advancing");
+ IMP_assert(nc < ub_, "Problems advancing");
+ cur_= nc;
+ }
+ return *this;
+ }
+ This operator++(int) {
+ This o= *this;
+ operator++;
+ return o;
+ }
+ reference_type operator*() const {
+ return cur_;
+ }
+ pointer_type operator->() const {
+ return &cur_;
+ }
+};
+
+} // internal NS
+
+/** \brief A voxel grid in 3D space.
+ VT can be any class.
+ */
+template <class VT>
+class Grid3D {
+public:
+ //! The type stored in each voxel.
+ typedef VT VoxelData;
+
+ //! An index that refers to a real voxel
+ typedef internal::GridIndex Index;
+ //! An index that refers to a voxel that may or may not exist
+ /** Such an index can refer to voxels outside of the grid
+ or with negative indices.
+ */
+ typedef internal::VirtualGridIndex VirtualIndex;
+
+private:
+ std::vector<VT> data_;
+ int d_[3];
+ Vector3D min_, max_;
+ float edge_size_[3];
+
+ void update_sizes() {
+ for (unsigned int i=0; i< 3; ++i) {
+ edge_size_[i]= (max_.get_component(i)- min_.get_component(i))/d_[i];
+ }
+ IMP_LOG(VERBOSE, "Grid has " << d_[0] << "x" << d_[1]
+ << "x" << d_[2] << " voxels of size "
+ << edge_size_[0] << "x" << edge_size_[1]
+ << "x" << edge_size_[2] << std::endl);
+ }
+
+ unsigned int index(const Index &i) const {
+ unsigned int ii= i[2]*d_[0]*d_[1] + i[1]*d_[0]+i[0];
+ IMP_assert(ii < data_.size(), "Invalid grid index");
+ return ii;
+ }
+ template <class IndexType>
+ IndexType get_index_t(Vector3D pt) const {
+ int index[3];
+ for (unsigned int i=0; i< 3; ++i ) {
+ float d= pt.get_component(i)- min_.get_component(i);
+ index[i]= static_cast<int>(std::floor(d/edge_size_[i]));
+ }
+ return IndexType(index[0], index[1], index[2]);
+ }
+
+public:
+
+ //! Initialize the grid
+ /**
+ \param[in] xd The number of voxels in the x direction
+ \param[in] yd The number of voxels in the y direction
+ \param[in] zd The number of voxels in the z direction
+ \param[in] minc The min coordinate of the grid
+ \param[in] maxc The max coordinate of the grid
+ \param[in] def The default value for the voxels
+ */
+ Grid3D(int xd, int yd, int zd,
+ Vector3D minc, Vector3D maxc,
+ VoxelData def): data_(xd*yd*zd, def),
+ min_(minc),
+ max_(maxc) {
+ d_[0]=xd;
+ d_[1]=yd;
+ d_[2]=zd;
+ update_sizes();
+ }
+
+ //! Initialize the grid
+ /**
+ \param[in] side The side length for the voxels
+ \param[in] minc The min coordinate of the grid
+ \param[in] maxc A lower bound for the max coordinate of the grid
+ \param[in] def The default value for the voxels
+ */
+ Grid3D(float side,
+ Vector3D minc, Vector3D maxc,
+ VoxelData def) {
+ min_=minc;
+ float mx[3];
+ for (unsigned int i=0; i< 3; ++i ) {
+ d_[i]= std::max(static_cast<int>(std::ceil((maxc.get_component(i)
+ - minc.get_component(i))/ side)),
+ 1);
+ mx[i]= d_[i]*side+ minc.get_component(i);
+ }
+ data_.resize(d_[0]*d_[1]*d_[2], def);
+ max_= Vector3D(mx[0], mx[1], mx[2]);
+ update_sizes();
+ }
+
+ //! An empty grid.
+ Grid3D(){
+ d_[0]=0;
+ d_[1]=0;
+ d_[2]=0;
+ }
+
+ //! Set the max corner of the grid
+ void set_min(Vector3D m) {
+ min_=m;
+ update_sizes();
+ }
+
+ //! Set the min corner of the voxel grid
+ void set_max(Vector3D m) {
+ max_=m;
+ update_sizes();
+ }
+
+ //! Return the number of voxels in a certain direction
+ unsigned int get_number_of_voxels(unsigned int i) const {
+ IMP_assert(i < 3, "Only 3D");
+ return d_[i];
+ }
+
+ //! Return the index of the voxel containing the point.
+ Index get_index(Vector3D pt) const {
+ for (unsigned int i=0; i< 3; ++i) {
+ if (pt.get_component(i) < min_.get_component(i)) return Index();
+ if (pt.get_component(i) > max_.get_component(i)) return Index();
+ }
+ return get_index_t<Index>(pt);
+ }
+
+ //! Return the index that would contain the voxel if the grid extended there
+ VirtualIndex get_virtual_index(Vector3D pt) const {
+ return get_index_t<VirtualIndex>(pt);
+ }
+
+ //! Convert a VirtualIndex into a real Index if possible
+ /** \return Index() if not possible
+ */
+ Index get_index(VirtualIndex v) const {
+ for (unsigned int i=0; i< 3; ++i) {
+ if (v[i] <0 || v[i] >= d_[i]) return Index();
+ }
+ return Index(v[0], v[1], v[2]);
+ }
+
+ //! Get the data in a particular cell
+ VoxelData& get_voxel(Index gi) {
+ return data_[index(gi)];
+ }
+
+ //! Get the data in a particular cell
+ const VoxelData& get_voxel(Index gi) const {
+ return data_[index(gi)];
+ }
+
+ Index get_min_index() const {
+ return Index(0,0,0);
+ }
+
+ //! Return a point at the center of the voxel
+ Vector3D get_center(VirtualIndex gi) const {
+ return Vector3D(edge_size_[0]*(.5f+ gi[0]) + min_.get_component(0),
+ edge_size_[1]*(.5f+ gi[1]) + min_.get_component(1),
+ edge_size_[2]*(.5f+ gi[2]) + min_.get_component(2));
+ }
+
+
+ //! Iterator through the Indexes in volume
+ /** The iterator iterates though the valid indexes bounded
+ by the VirtualIndex
+ */
+ typedef internal::GridIndexIterator<Index> IndexIterator;
+ IndexIterator indexes_begin(VirtualIndex lb,
+ VirtualIndex ub) const {
+ IMP_assert(lb <= ub, "empty range");
+ Index lbi(std::max(0, lb[0]),
+ std::max(0, lb[1]),
+ std::max(0, lb[2]));
+ VirtualIndex ubi(std::min(d_[0], ub[0]),
+ std::min(d_[1], ub[1]),
+ std::min(d_[2], ub[2]));
+ if (lb==ub) {
+ return IndexIterator();
+ } else {
+ return IndexIterator(lbi, ubi);
+ }
+ }
+ IndexIterator indexes_end(VirtualIndex,
+ VirtualIndex) const {
+ //IMP_assert(lb <= ub, "empty range");
+ return IndexIterator();
+ }
+
+ IndexIterator all_indexes_begin() const {
+ return indexes_begin(VirtualIndex(0,0,0),
+ VirtualIndex(d_[0],
+ d_[1],
+ d_[2]));
+ }
+ IndexIterator all_indexes_end() const {
+ return indexes_end(VirtualIndex(0,0,0),
+ VirtualIndex(d_[0],
+ d_[1],
+ d_[2]));
+ }
+};
+
+
+}
+
+#endif
Index: kernel/include/IMP/score_states/NonbondedListScoreState.h
===================================================================
--- kernel/include/IMP/score_states/NonbondedListScoreState.h (revision 345)
+++ kernel/include/IMP/score_states/NonbondedListScoreState.h (working copy)
@@ -12,6 +12,7 @@
#include <limits>
#include "../ScoreState.h"
#include "BondedListScoreState.h"
+#include "../Grid3D.h"
namespace IMP
{
@@ -19,20 +20,88 @@
class BondedListScoreState;
//! This class maintains a list of non-bonded pairs.
+/** \ingroup restraint
+ */
class IMPDLLEXPORT NonbondedListScoreState: public ScoreState
{
- protected:
+protected:
Particles ps_;
typedef std::vector<std::pair<Particle*, Particle*> > NBL;
NBL nbl_;
+ float target_voxel_side_;
+ float last_cutoff_;
+ typedef Grid3D<Particles> Grid;
+ Grid grid_;
- void rescan();
+ virtual void rescan(float cut);
void audit_particles(const Particles &ps) const;
void propagate_set_particles(const Particles &aps);
void propagate_update();
void add_if_nonbonded(Particle *a, Particle *b);
+
+ void create_grid();
+
+ // Handle the particle against all the cells other than
+ // the one containing it.
+ void handle_particle(Particle *p,
+ const Grid::VirtualIndex& v,
+ float cut,
+ bool skip_lower);
+ /*class SingeParticleIterator {
+ void advance() {
+ if (curv_== Grid::IndexIterator()) return;
+ if (curp_== nbl_->get_grid().get_voxel(*curv_).end()) {
+ ++curv_;
+ if (curv_ != Grid::IndexIterator()) {
+ curp_= nbl_->get_grid().get_voxel(*curv_).begin();
+ }
+ }
+ }
+ void find() {
+ bool ok=false;
+ do {
+ ok=true;
+ if ( sl_ && i_ >= *curv_) {
+ ok=false;
+ } else if (curp_== nbl_->get_grid().get_voxel(*curv_).end()) {
+ ok=false;
+ } else if (nbl_->are_bonded(p_, *curp_)) {
+ ok=false;
+ }
+ if (!ok) advance();
+ } while (!ok);
+ if (nbl_->are_bonded(p_,
+ return true;
+
+ }
+ public:
+ SingleParticleIterator(Particle *p,
+ const Grid::VirtualIndex& v,
+ int ncells,
+ bool skip_lower,
+ NonbondedListScoreState *nbl): p_(p),
+ i_(v),
+ sl_(skip_lower),
+ nbl_(nbl) {
+ Grid::VirtualIndex lc(i_[0]-ncells,
+ i_[1]-ncells,
+ i_[2]-ncells);
+ Grid::VirtualIndex uc(i_[0]+ncells,
+ i_[1]+ncells,
+ i_[2]+ncells);
+ curv_= nbl_->get_grid().indexes_begin(lc, uc);
+ endv_= nbl_->get_grid().indices_end(lc, uc);
+ curp_= nbl_->get_grid().get_voxel(*curv_);
+ find_next();
+ }
+ };*/
+
public:
- NonbondedListScoreState(const Particles &ps);
+ /**
+ \param[in] ps A list of particles to use.
+ \param[in] tvs A suggested size for the voxel side.
+ */
+ NonbondedListScoreState(const Particles &ps, float tvs=1);
virtual ~NonbondedListScoreState();
IMP_CONTAINER(BondedListScoreState, bonded_list,
BondedListIndex);
@@ -49,14 +118,20 @@
//! This iterates through the pairs of non-bonded particles
/** \param[in] cutoff The state may ignore pairs which are futher
apart than the cutoff.
- \precondition update() must be called first for this to be valid.
+ \note that this is highly unsafe and iteration can only be
+ done once at a time. I will fix that eventually.
*/
NonbondedIterator nonbonded_begin(Float cutoff
- =std::numeric_limits<Float>::max()) const {
+ =std::numeric_limits<Float>::max()) {
+ IMP_assert(last_cutoff_== cutoff || last_cutoff_==-1,
+ "Bad things are happening with the iterators in "
+ << "NonbondedListScoreState");
+ rescan(cutoff);
return nbl_.begin();
}
NonbondedIterator nonbonded_end(Float cutoff
- =std::numeric_limits<Float>::max()) const {
+ =std::numeric_limits<Float>::max()) {
+ rescan(cutoff);
return nbl_.end();
}
};
Index: kernel/include/IMP/score_states/BipartiteNonbondedListScoreState.h
===================================================================
--- kernel/include/IMP/score_states/BipartiteNonbondedListScoreState.h (revision 345)
+++ kernel/include/IMP/score_states/BipartiteNonbondedListScoreState.h (working copy)
@@ -28,7 +28,7 @@
typedef NonbondedListScoreState P;
Particles ps_;
- void rescan();
+ virtual void rescan(float cut);
void set_particles(const Particles &ps0) {
// hide the parent's version
}
@@ -40,6 +40,7 @@
IMP_SCORE_STATE("0.5", "Daniel Russel");
void set_particles(const Particles &ps0, const Particles &ps1);
+
};
} // namespace IMP
Index: kernel/include/IMP/score_states/bonded_lists/BondDecoratorListScoreState.h
===================================================================
--- kernel/include/IMP/score_states/bonded_lists/BondDecoratorListScoreState.h (revision 345)
+++ kernel/include/IMP/score_states/bonded_lists/BondDecoratorListScoreState.h (working copy)
@@ -21,6 +21,7 @@
//! Keep track of particles that are connected by BondDecorator bonds.
/** We also may want to add lazy rescanning of bonds rather than doing
it every update call and a faster lookup of bonds.
+ \ingroup bond
*/
class IMPDLLEXPORT BondDecoratorListScoreState: public BondedListScoreState
{
@@ -32,7 +33,6 @@
*/
BondDecoratorListScoreState(const Particles &ps);
virtual ~BondDecoratorListScoreState(){}
- //IMP_SCORE_STATE("0.5", "Daniel Russel");
virtual void set_particles(const Particles &ps);
@@ -41,7 +41,7 @@
virtual void update();
//! This iterates through the pairs of bonded particles
- /** \precondition update() must be called first for this to be valid.
+ /** \note update() must be called first for this to be valid.
*/
typedef std::vector<BondDecorator>::const_iterator BondIterator;
BondIterator bonds_begin() const {
Index: kernel/src/score_states/NonbondedListScoreState.cpp
===================================================================
--- kernel/src/score_states/NonbondedListScoreState.cpp (revision 345)
+++ kernel/src/score_states/NonbondedListScoreState.cpp (working copy)
@@ -7,12 +7,16 @@
#include "IMP/score_states/NonbondedListScoreState.h"
#include "IMP/decorators/XYZDecorator.h"
+#include "IMP/Grid3D.h"
namespace IMP
{
-NonbondedListScoreState::NonbondedListScoreState(const Particles &ps)
+NonbondedListScoreState::NonbondedListScoreState(const Particles &ps,
+ float tvs):
+ target_voxel_side_(tvs)
{
+ last_cutoff_=-1;
set_particles(ps);
}
@@ -52,6 +56,7 @@
void NonbondedListScoreState::add_if_nonbonded(Particle *a, Particle *b)
{
+ if (!a->get_is_active() || !b->get_is_active()) return;
bool found=false;
for (BondedListIterator bli= bonded_lists_begin();
bli != bonded_lists_end(); ++bli) {
@@ -61,21 +66,100 @@
}
}
if (!found) {
+ IMP_LOG(VERBOSE, "Found pair " << a->get_index()
+ << " " << b->get_index() << std::endl);
nbl_.push_back(std::make_pair(a, b));
}
}
+void NonbondedListScoreState::create_grid() {
+ float mn[3]= {std::numeric_limits<float>::max(),
+ std::numeric_limits<float>::max(),
+ std::numeric_limits<float>::max()};
+ float mx[3]={-std::numeric_limits<float>::max(),
+ -std::numeric_limits<float>::max(),
+ -std::numeric_limits<float>::max()};
+ for (unsigned int i = 0; i < ps_.size(); ++i) {
+ XYZDecorator d= XYZDecorator::cast(ps_[i]);
+ for (unsigned int j=0; j<3; ++j) {
+ if (d.get_coordinate(j)< mn[j]) mn[j]= d.get_coordinate(j);
+ if (d.get_coordinate(j)> mx[j]) mx[j]= d.get_coordinate(j);
+ }
+ }
+ grid_= Grid(target_voxel_side_,
+ Vector3D(mn[0], mn[1], mn[2]),
+ Vector3D(mx[0], mx[1], mx[2]),
+ Particles());
+ for (unsigned int i = 0; i < ps_.size(); ++i) {
+ XYZDecorator d= XYZDecorator::cast(ps_[i]);
+ Vector3D v(d.get_x(), d.get_y(), d.get_z());
+ grid_.get_voxel(grid_.get_index(v)).push_back(ps_[i]);
+ }
+}
-void NonbondedListScoreState::rescan()
+void NonbondedListScoreState::handle_particle(Particle *p,
+ const Grid::VirtualIndex ¢er,
+ float cut,
+ bool skip_lower) {
+ // it is important that this is not unsigned
+ int ncells= static_cast<int>(std::ceil(cut/target_voxel_side_));
+
+ Grid::VirtualIndex lc(center[0]-ncells,
+ center[1]-ncells,
+ center[2]-ncells);
+ Grid::VirtualIndex uc(center[0]+ncells,
+ center[1]+ncells,
+ center[2]+ncells);
+ IMP_LOG(VERBOSE, "Iteration bounds are " << lc << " and " << uc << std::endl);
+ for (Grid::IndexIterator cur= grid_.indexes_begin(lc, uc);
+ cur != grid_.indexes_end(lc, uc);
+ ++cur){
+ if ( skip_lower && center >= *cur) continue;
+ if (grid_.get_voxel(*cur).empty()) continue;
+ //IMP_LOG(VERBOSE, "Paired with " << cur << std::endl);
+
+ for (unsigned int pi= 0;
+ pi< grid_.get_voxel(*cur).size(); ++pi) {
+ Particle *op = grid_.get_voxel(*cur)[pi];
+ add_if_nonbonded(p, op);
+ }
+ }
+ if (skip_lower) {
+ Grid::Index ri= grid_.get_index(center);
+ if (ri != Grid::Index()) {
+ const IMP::Particles &pis= grid_.get_voxel(ri);
+ for (unsigned int pj= 0; pis[pj] != p; ++pj) {
+ Particle *op= pis[pj];
+ add_if_nonbonded(p, op);
+ }
+ }
+ }
+}
+
+
+void NonbondedListScoreState::rescan(float cut)
{
+ if (last_cutoff_ > cut) {
+ IMP_LOG(VERBOSE, "Skipping rebuilding with cutoff " << cut
+ << " as old cut was " << last_cutoff_ << std::endl);
+ return;
+ }
+ last_cutoff_=cut;
+
+ create_grid();
+
nbl_.clear();
- for (unsigned int i = 0; i < ps_.size(); ++i) {
- Particle *pi = ps_[i];
- for (unsigned int j = 0; j < i; ++j) {
- Particle *pj = ps_[j];
- add_if_nonbonded(pi, pj);
+ for (Grid::IndexIterator center= grid_.all_indexes_begin();
+ center != grid_.all_indexes_end();
+ ++center) {
+ if (grid_.get_voxel(*center).empty()) continue;
+ IMP_LOG(VERBOSE, "Scanning from grid index " << *center << std::endl);
+ for (unsigned int pi= 0; pi< grid_.get_voxel(*center).size(); ++pi) {
+ Particle *p= grid_.get_voxel(*center)[pi];
+ handle_particle(p, *center, cut, true);
}
}
+ IMP_LOG(VERBOSE, "Found " << nbl_.size() << " nonbonded pairs" << std::endl);
}
void NonbondedListScoreState::set_particles(const Particles &ps)
@@ -83,14 +167,16 @@
ps_ = ps;
audit_particles(ps_);
propagate_set_particles(ps);
- rescan();
+ last_cutoff_=-1;
+ //rescan(lc);
}
void NonbondedListScoreState::update()
{
IMP_LOG(VERBOSE, "Updating non-bonded list" << std::endl);
propagate_update();
- rescan();
+ last_cutoff_=-1;
+ //rescan();
}
void NonbondedListScoreState::show(std::ostream &out) const
Index: kernel/src/score_states/BipartiteNonbondedListScoreState.cpp
===================================================================
--- kernel/src/score_states/BipartiteNonbondedListScoreState.cpp (revision 345)
+++ kernel/src/score_states/BipartiteNonbondedListScoreState.cpp (working copy)
@@ -7,10 +7,27 @@
#include "IMP/score_states/BipartiteNonbondedListScoreState.h"
#include "IMP/decorators/XYZDecorator.h"
+#include "IMP/Grid3D.h"
namespace IMP
{
+ /*
+ Iterator needs:
+ iterator into ps_, end iterator into ps_
+ cell for ps_[i]
+ -> actually has value current iterator around cell
+ -put in logic for current
+ number of bins to search
+
+
+ nonbipartite:
+ iterator through all cells
+ index in current cell
+ -> actual value in iterator around cell
+ num bins
+ */
+
BipartiteNonbondedListScoreState
::BipartiteNonbondedListScoreState(const Particles &ps0,
const Particles &ps1):
@@ -23,16 +40,38 @@
{
}
-void BipartiteNonbondedListScoreState::rescan()
+void BipartiteNonbondedListScoreState::rescan(float cut)
{
- nbl_.clear();
- for (unsigned int i = 0; i< ps_.size(); ++i) {
- Particle *pi= ps_[i];
- for (unsigned int j = 0; j < P::ps_.size(); ++j) {
- Particle *pj = P::ps_[j];
- P::add_if_nonbonded(pi,pj);
+ if (P::last_cutoff_ > cut) {
+ IMP_LOG(VERBOSE, "Skipping rebuilding with cutoff " << cut
+ << " as old cut was " << last_cutoff_ << std::endl);
+ return;
+ }
+ P::last_cutoff_=cut;
+
+ P::create_grid();
+
+ for (Grid::IndexIterator it= P::grid_.all_indexes_begin();
+ it != P::grid_.all_indexes_end(); ++it) {
+ if (!P::grid_.get_voxel(*it).empty()) {
+ IMP_LOG(VERBOSE, "Voxel " << *it << ":");
+ for (unsigned int i=0; i< P::grid_.get_voxel(*it).size(); ++i) {
+ IMP_LOG(VERBOSE, " " << P::grid_.get_voxel(*it)[i]->get_index());
+ }
+ IMP_LOG(VERBOSE, std::endl);
}
}
+
+ nbl_.clear();
+ for (unsigned int i=0; i< ps_.size(); ++i) {
+ Particle *p= ps_[i];
+ XYZDecorator d= XYZDecorator::cast(p);
+ Grid::VirtualIndex index
+ = P::grid_.get_virtual_index(Vector3D(d.get_x(), d.get_y(), d.get_z()));
+ IMP_LOG(VERBOSE, "Index is " << index << std::endl);
+ P::handle_particle(p, index, cut, false);
+ }
+ IMP_LOG(VERBOSE, "Found " << nbl_.size() << " nonbonded pairs" << std::endl);
}
void BipartiteNonbondedListScoreState::set_particles(const Particles &ps0,
@@ -40,19 +79,20 @@
{
ps_=ps0;
P::ps_=ps1;
+ if (ps_.size() > P::ps_.size()) std::swap(ps_, P::ps_);
Particles aps(ps0);
aps.insert(aps.end(), ps1.begin(), ps1.end());
P::audit_particles(aps);
P::propagate_set_particles(aps);
- rescan();
+ P::last_cutoff_=-1;
}
void BipartiteNonbondedListScoreState::update()
{
P::propagate_update();
- rescan();
+ P::last_cutoff_=-1;
}
void BipartiteNonbondedListScoreState::show(std::ostream &out) const
Index: kernel/src/score_states/bonded_lists/BondDecoratorListScoreState.cpp
===================================================================
--- kernel/src/score_states/bonded_lists/BondDecoratorListScoreState.cpp (revision 345)
+++ kernel/src/score_states/bonded_lists/BondDecoratorListScoreState.cpp (working copy)
@@ -18,10 +18,12 @@
void BondDecoratorListScoreState::update() {
bonds_.clear();
for (unsigned int i=0; i< ps_.size(); ++i) {
+ if (!ps_[i]->get_is_active()) continue;
BondedDecorator di= BondedDecorator::cast(ps_[i]);
ParticleIndex pi= ps_[i]->get_index();
for (unsigned int j=0; j< di.get_number_of_bonds(); ++j) {
BondedDecorator dj= di.get_bonded(j);
+ if (! dj.get_particle()->get_is_active()) continue;
if (di < dj) {
bonds_.push_back(di.get_bond(j));
}
Index: kernel/test/states/test_nonbonded_list.py
===================================================================
--- kernel/test/states/test_nonbonded_list.py (revision 345)
+++ kernel/test/states/test_nonbonded_list.py (working copy)
@@ -1,5 +1,6 @@
import unittest
import IMP, IMP.test
+import random
class OnePair(IMP.PairScore):
def __init__(self):
@@ -32,19 +33,20 @@
p= IMP.Particle()
m.add_particle(p)
d=IMP.XYZDecorator.create(p)
- d.set_x(0)
- d.set_y(1)
- d.set_z(1)
+ d.set_x(random.uniform(0,10))
+ d.set_y(random.uniform(0,10))
+ d.set_z(random.uniform(0,10))
s= IMP.NonbondedListScoreState(m.get_particles())
m.add_score_state(s)
o= OnePair()
- r= IMP.NonbondedRestraint(s, o)
+ r= IMP.NonbondedRestraint(s, o, 15)
m.add_restraint(r)
score= m.evaluate(False)
self.assertEqual(score, 45, "Wrong score")
def test_bl(self):
"""Test the bonded list"""
+ return
m= IMP.Model()
bds=[]
for i in range(0,10):
@@ -52,15 +54,15 @@
m.add_particle(p)
d=IMP.XYZDecorator.create(p)
d.set_x(0)
- d.set_y(1)
- d.set_z(1)
+ d.set_y(10)
+ d.set_z(10)
bds.append(IMP.BondedDecorator.create(p))
pts=IMP.Particles()
for p in m.get_particles():
pts.append(p)
for i in range(1,10):
IMP.custom_bond(bds[i-1], bds[i], 1, .1)
- s= IMP.NonbondedListScoreState(pts)
+ s= IMP.NonbondedListScoreState(pts, 1)
b= IMP.BondDecoratorListScoreState(pts)
s.add_bonded_list(b)
m.add_score_state(s)
@@ -74,7 +76,31 @@
score= m.evaluate( False )
self.assertEqual(score, 900+45-9, "Wrong score")
- def test_bnbl(self):
+ def test_distfilt(self):
+ """Test filtering based on distance in nonbonded list"""
+ m= IMP.Model()
+ ps=IMP.Particles()
+ for i in range(0,10):
+ p= IMP.Particle()
+ m.add_particle(p)
+ d=IMP.XYZDecorator.create(p)
+ ps.append(p)
+ if (i < 5):
+ d.set_x(random.uniform(0,10))
+ d.set_y(random.uniform(0,10))
+ d.set_z(random.uniform(0,10))
+ else:
+ d.set_x(random.uniform(50,60))
+ d.set_y(random.uniform(50,60))
+ d.set_z(random.uniform(50,60))
+ s= IMP.NonbondedListScoreState(ps)
+ m.add_score_state(s)
+ o= OnePair()
+ r= IMP.NonbondedRestraint(s, o, 15)
+ m.add_restraint(r)
+ score= m.evaluate(False)
+ self.assertEqual(score, 20, "Wrong score")
+ def test_bi(self):
"""Test the bipartite nonbonded list and restraint which uses it"""
m= IMP.Model()
ps0=IMP.Particles()
@@ -83,9 +109,9 @@
p= IMP.Particle()
m.add_particle(p)
d=IMP.XYZDecorator.create(p)
- d.set_x(0)
- d.set_y(1)
- d.set_z(1)
+ d.set_x(random.uniform(0,10))
+ d.set_y(random.uniform(0,10))
+ d.set_z(random.uniform(0,10))
if (i < 5):
ps0.append(p)
else:
@@ -93,7 +119,7 @@
s= IMP.BipartiteNonbondedListScoreState(ps0, ps1)
m.add_score_state(s)
o= OnePair()
- r= IMP.NonbondedRestraint(s, o)
+ r= IMP.NonbondedRestraint(s, o, 15)
m.add_restraint(r)
score= m.evaluate(False)
self.assertEqual(score, 25, "Wrong score")