Index: kernel/test/states/test_nonbonded_list.py =================================================================== --- kernel/test/states/test_nonbonded_list.py (revision 354) +++ 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,20 +54,20 @@ 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) o= OnePair() - r= IMP.NonbondedRestraint(s, o) + r= IMP.NonbondedRestraint(s, o, 10) os=OneScore() print os(6) br= IMP.BondDecoratorRestraint(b, os) @@ -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,17 +109,17 @@ 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: ps1.append(p) - s= IMP.BipartiteNonbondedListScoreState(ps0, ps1) + s= IMP.BipartiteNonbondedListScoreState(ps0, ps1, 1000) 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") Index: kernel/include/IMP/score_states/NonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/NonbondedListScoreState.h (revision 354) +++ kernel/include/IMP/score_states/NonbondedListScoreState.h (working copy) @@ -10,29 +10,102 @@ #include #include +#include "BondedListScoreState.h" +#include "MaxChangeScoreState.h" #include "../ScoreState.h" -#include "BondedListScoreState.h" +#include "../Grid3D.h" namespace IMP { class BondedListScoreState; +class MaxChangeScoreState; //! This class maintains a list of non-bonded pairs. +/** \ingroup restraint + */ class IMPDLLEXPORT NonbondedListScoreState: public ScoreState { - protected: - Particles ps_; +protected: + float target_voxel_side_; typedef std::vector > NBL; NBL nbl_; + float last_cutoff_; + bool grid_valid_; + typedef Grid3D Grid; + Grid grid_; + std::auto_ptr mc_; - 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 invalidate_nbl() { + last_cutoff_=-1; + } + 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,16 +122,26 @@ //! 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::max()) const { + =std::numeric_limits::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::max()) const { + =std::numeric_limits::max()) { + rescan(cutoff); return nbl_.end(); } + + const Particles &get_particles() const { + return mc_->get_particles(); + } }; } // namespace IMP Index: kernel/include/IMP/score_states/BipartiteNonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/BipartiteNonbondedListScoreState.h (revision 354) +++ kernel/include/IMP/score_states/BipartiteNonbondedListScoreState.h (working copy) @@ -26,20 +26,24 @@ public NonbondedListScoreState { typedef NonbondedListScoreState P; - Particles ps_; + std::auto_ptr mc_; - void rescan(); + virtual void rescan(float cut); void set_particles(const Particles &ps0) { // hide the parent's version } public: BipartiteNonbondedListScoreState(const Particles &ps0, - const Particles &ps1); + const Particles &ps1, + float target_side); virtual ~BipartiteNonbondedListScoreState(); IMP_SCORE_STATE("0.5", "Daniel Russel"); void set_particles(const Particles &ps0, const Particles &ps1); + const Particles &get_particles() const { + return mc_->get_particles(); + } }; } // namespace IMP Index: kernel/include/IMP/score_states/bonded_lists/BondDecoratorListScoreState.h =================================================================== --- kernel/include/IMP/score_states/bonded_lists/BondDecoratorListScoreState.h (revision 354) +++ 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::const_iterator BondIterator; BondIterator bonds_begin() const { Index: kernel/include/IMP/score_states/MaxChangeScoreState.h =================================================================== --- kernel/include/IMP/score_states/MaxChangeScoreState.h (revision 0) +++ kernel/include/IMP/score_states/MaxChangeScoreState.h (revision 0) @@ -0,0 +1,70 @@ +/** + * \file MaxChangeScoreState.h + * \brief Keep track of the maximum change of a set of attributes. + * + * Copyright 2007-8 Sali Lab. All rights reserved. + */ + +#ifndef __IMP_MAX_CHANGE_SCORE_STATE_H +#define __IMP_MAX_CHANGE_SCORE_STATE_H + +#include "../ScoreState.h" +#include "../Index.h" +#include "../Particle.h" + +#include + +namespace IMP +{ + +//! This abstract class maintains a list of bonded pairs. +class IMPDLLEXPORT MaxChangeScoreState: public ScoreState +{ + FloatKeys keys_; + std::vector > orig_; + Particles ps_; + float max_change_; + public: + MaxChangeScoreState(const FloatKeys &keys):keys_(keys){} + virtual ~MaxChangeScoreState(){} + + virtual void update() { + max_change_=0; + for (unsigned int i=0; i < orig_.size(); ++i) { + for (unsigned int j=0; j < keys_.size(); ++j) { + max_change_= std::max(max_change_, + std::abs(ps_[i]->get_value(keys_[j]) + - orig_[i][j])); + } + } + } + + void reset() { + orig_.resize(ps_.size(), std::vector(keys_.size(), 0)); + for (unsigned int i=0; i < orig_.size(); ++i) { + for (unsigned int j=0; j < keys_.size(); ++j) { + orig_[i][j]=-ps_[i]->get_value(keys_[j]); + } + } + max_change_=0; + } + + float get_max() const { + return max_change_; + } + + //! Set the set of particles used + void set_particles(const Particles &ps) { + ps_=ps; + reset(); + } + + const Particles &get_particles() const { + return ps_; + } + +}; + +} // namespace IMP + +#endif /* __IMP_BONDED_LIST_SCORE_STATE_H */ Index: kernel/include/IMP/restraints/NonbondedRestraint.h =================================================================== --- kernel/include/IMP/restraints/NonbondedRestraint.h (revision 354) +++ kernel/include/IMP/restraints/NonbondedRestraint.h (working copy) @@ -21,6 +21,9 @@ class PairScore; //! Restrain all pairs of non-bonded particles +/** + \ingroup restraint + */ class IMPDLLEXPORT NonbondedRestraint : public Restraint { public: Index: kernel/include/IMP/restraints/BondDecoratorRestraint.h =================================================================== --- kernel/include/IMP/restraints/BondDecoratorRestraint.h (revision 354) +++ kernel/include/IMP/restraints/BondDecoratorRestraint.h (working copy) @@ -27,6 +27,8 @@ The bond stiffness is assumed to be 1 if it is not specified in the bond. This can become a parameter eventually. + + \ingroup bond */ class IMPDLLEXPORT BondDecoratorRestraint : public Restraint { @@ -36,7 +38,7 @@ of bonds. \param[in] f The UnaryFunction to apply to the particles. */ - BondDecoratorRestraint(BondDecoratorListScoreState *nbl, UnaryFunction *f); + BondDecoratorRestraint(BondDecoratorListScoreState *bl, UnaryFunction *f); virtual ~BondDecoratorRestraint(){} IMP_RESTRAINT("0.5", "Daniel Russel"); Index: kernel/include/IMP/Grid3D.h =================================================================== --- kernel/include/IMP/Grid3D.h (revision 0) +++ kernel/include/IMP/Grid3D.h (revision 0) @@ -0,0 +1,371 @@ +/** + * \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 Grid3D; + +namespace internal { + +template +class GridIndexIterator; + +class VirtualGridIndex { + typedef VirtualGridIndex This; + int d_[3]; + bool is_default() const { + return d_[0]==std::numeric_limits::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::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] << ")"; + } + bool strictly_larger_than(const VirtualGridIndex &o) const { + return d_[0] > o.d_[0] && d_[1] > o.d_[1] && d_[2] > o.d_[2]; + } +}; + +IMP_OUTPUT_OPERATOR(VirtualGridIndex); + +class GridIndex: public VirtualGridIndex { +public: + GridIndex(): VirtualGridIndex() { + } +protected: + template + friend class IMP::Grid3D; + template + 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 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(ub_.strictly_larger_than(lb_), + "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(ub_.strictly_larger_than(nc), "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 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 data_; + int d_[3]; + Vector3D min_, max_; + float edge_size_[3]; + + void update_sizes() { + for (unsigned int i=0; i< 3; ++i) { + // hack to try to handle roundoff errors + // I would like to find something more reliable + edge_size_[i]= 1.05*(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 + 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(std::floor(d/edge_size_[i])); + } + return IndexType(index[0], index[1], index[2]); + } + + int snap(int dim, int v) const { + if (v < 0) return 0; + else if (v > d_[dim]) return d_[dim]; + else return v; + } + + Index snap(const VirtualIndex &v) const { + return Index(snap(0, v[0]), + snap(1, v[1]), + snap(2, v[2])); + } + std::pair empty_range() const { + return std::make_pair(Index(0,0,0), VirtualIndex(0,0,0)); + } + + std::pair intersect(VirtualIndex l, + VirtualIndex u) const { + Index rlb; + VirtualIndex rub; + for (unsigned int i=0; i< 3; ++i) { + if (u[i] <= 0) return empty_range(); + if (l[i] >= d_[i]) return empty_range(); + } + return std::make_pair(snap(l), snap(u)); + } + + +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(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(); + } + + //! Get the min corner + const Vector3D &get_min() const { + return min_; + } + + //! Get the max corner + const Vector3D &get_max() const { + return max_; + } + + //! 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(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(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 IndexIterator; + IndexIterator indexes_begin(VirtualIndex lb, + VirtualIndex ub) const { + std::pair bp= intersect(lb,ub); + if (bp.first== bp.second) { + return IndexIterator(); + } else { + IMP_assert(bp.second.strictly_larger_than(bp.first), "empty range"); + return IndexIterator(bp.first, bp.second); + } + } + 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/src/restraints/NonbondedRestraint.cpp =================================================================== --- kernel/src/restraints/NonbondedRestraint.cpp (revision 354) +++ kernel/src/restraints/NonbondedRestraint.cpp (working copy) @@ -30,11 +30,20 @@ IMP_CHECK_OBJECT(sf_.get()); IMP_CHECK_OBJECT(nbl_); Float score=0; - + IMP_LOG(VERBOSE, "Nonbonded restraint on " + << std::distance(nbl_->nonbonded_begin(max_dist_), + nbl_->nonbonded_end(max_dist_)) + << " pairs" << std::endl); for (NonbondedListScoreState::NonbondedIterator it = nbl_->nonbonded_begin(max_dist_); it != nbl_->nonbonded_end(max_dist_); ++it) { - score += sf_->evaluate(it->first, it->second, accum); + float thisscore = sf_->evaluate(it->first, it->second, accum); + if (thisscore != 0) { + IMP_LOG(VERBOSE, "Pair " << it->first->get_index() + << " and " << it->second->get_index() << " have score " + << thisscore << std::endl); + } + score+= thisscore; } return score; Index: kernel/src/restraints/BondDecoratorRestraint.cpp =================================================================== --- kernel/src/restraints/BondDecoratorRestraint.cpp (revision 354) +++ kernel/src/restraints/BondDecoratorRestraint.cpp (working copy) @@ -33,6 +33,10 @@ } if (s <0) s=1; try { + /*IMP_LOG(VERBOSE, "Bonded pair " + << bd.get_bonded(0).get_particle()->get_index() + << " " << bd.get_bonded(1).get_particle()->get_index() + << " with length " << l << " and stiffness " << s << std::endl);*/ sum+= internal::evaluate_distance_pair_score(bd.get_bonded(0).get_particle(), bd.get_bonded(1).get_particle(), Index: kernel/src/score_states/NonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/NonbondedListScoreState.cpp (revision 354) +++ kernel/src/score_states/NonbondedListScoreState.cpp (working copy) @@ -7,12 +7,22 @@ #include "IMP/score_states/NonbondedListScoreState.h" #include "IMP/decorators/XYZDecorator.h" +#include "IMP/Grid3D.h" +#include "IMP/score_states/MaxChangeScoreState.h" namespace IMP { -NonbondedListScoreState::NonbondedListScoreState(const Particles &ps) +NonbondedListScoreState::NonbondedListScoreState(const Particles &ps, + float tvs): + target_voxel_side_(tvs), grid_valid_(false) { + last_cutoff_=-1; + FloatKeys fks; + fks.push_back(FloatKey("x")); + fks.push_back(FloatKey("y")); + fks.push_back(FloatKey("z")); + mc_= std::auto_ptr(new MaxChangeScoreState(fks)); set_particles(ps); } @@ -34,24 +44,9 @@ } } -void NonbondedListScoreState::propagate_set_particles(const Particles &aps) -{ - for (BondedListScoreStateIterator bli= bonded_lists_begin(); - bli != bonded_lists_end(); ++bli) { - (*bli)->set_particles(aps); - } -} - -void NonbondedListScoreState::propagate_update() -{ - for (BondedListScoreStateIterator bli= bonded_lists_begin(); - bli != bonded_lists_end(); ++bli) { - (*bli)->update(); - } -} - 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,36 +56,145 @@ } } 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() { + IMP_LOG(VERBOSE, "Creating nonbonded grid..." << std::flush); + float mn[3]= {std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max()}; + float mx[3]={-std::numeric_limits::max(), + -std::numeric_limits::max(), + -std::numeric_limits::max()}; + for (unsigned int i = 0; i < get_particles().size(); ++i) { + XYZDecorator d= XYZDecorator::cast(get_particles()[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 < get_particles().size(); ++i) { + XYZDecorator d= XYZDecorator::cast(get_particles()[i]); + Vector3D v(d.get_x(), d.get_y(), d.get_z()); + grid_.get_voxel(grid_.get_index(v)).push_back(get_particles()[i]); + } + grid_valid_=true; + mc_->reset(); + invalidate_nbl(); + IMP_LOG(VERBOSE, "done." << std::endl); +} -void NonbondedListScoreState::rescan() +void NonbondedListScoreState::handle_particle(Particle *p, + const Grid::VirtualIndex ¢er, + float cut, + bool skip_lower) { + Grid::VirtualIndex lc, uc; + + if ( cut > target_voxel_side_*1000 ) { + // This is needed to handle overflow + lc=Grid::VirtualIndex(0,0,0); + uc=Grid::VirtualIndex(grid_.get_number_of_voxels(0), + grid_.get_number_of_voxels(1), + grid_.get_number_of_voxels(2)); + } else { + // it is important that this is not unsigned + int ncells= static_cast(std::ceil(cut/target_voxel_side_)); + // to allow laziness in rebuilding + ++ncells; + lc=Grid::VirtualIndex(center[0]-ncells, + center[1]-ncells, + center[2]-ncells); + uc=Grid::VirtualIndex(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) { + return; + } + IMP_LOG(VERBOSE, "Rescanning " << cut + << "(" << last_cutoff_ << ")" << std::endl); + + last_cutoff_=cut; + 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); + IMP_LOG(VERBOSE, "Scanning for nonbonded pairs..." << std::flush); + 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, "done" << std::endl); + //IMP_LOG(VERBOSE, "Found " << nbl_.size() << " nonbonded pairs" + //<< std::endl); + IMP_LOG(VERBOSE, "Found " << nbl_.size() << " nonbonded pairs" + << " for " << get_particles().size() + << " particles" << std::endl); } void NonbondedListScoreState::set_particles(const Particles &ps) { - ps_ = ps; - audit_particles(ps_); - propagate_set_particles(ps); - rescan(); + audit_particles(get_particles()); + for (BondedListScoreStateIterator bli= bonded_lists_begin(); + bli != bonded_lists_end(); ++bli) { + (*bli)->set_particles(ps); + } + mc_->set_particles(ps); + grid_valid_=false; } void NonbondedListScoreState::update() { IMP_LOG(VERBOSE, "Updating non-bonded list" << std::endl); - propagate_update(); - rescan(); + for (BondedListScoreStateIterator bli= bonded_lists_begin(); + bli != bonded_lists_end(); ++bli) { + (*bli)->update(); + } + if (!grid_valid_ || mc_->get_max() > target_voxel_side_) { + create_grid(); + } } void NonbondedListScoreState::show(std::ostream &out) const Index: kernel/src/score_states/BipartiteNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/BipartiteNonbondedListScoreState.cpp (revision 354) +++ kernel/src/score_states/BipartiteNonbondedListScoreState.cpp (working copy) @@ -7,15 +7,37 @@ #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): - NonbondedListScoreState(ps0) + const Particles &ps1, float target): + NonbondedListScoreState(ps0, target) { + FloatKeys fks; + fks.push_back(FloatKey("x")); + fks.push_back(FloatKey("y")); + fks.push_back(FloatKey("z")); + mc_= std::auto_ptr(new MaxChangeScoreState(fks)); set_particles(ps0, ps1); } @@ -23,36 +45,58 @@ { } -void BipartiteNonbondedListScoreState::rescan() +void BipartiteNonbondedListScoreState::rescan(float cut) { + if (P::last_cutoff_ >= cut) { + return; + } + P::last_cutoff_=cut; + + /*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 *pi= ps_[i]; - for (unsigned int j = 0; j < P::ps_.size(); ++j) { - Particle *pj = P::ps_[j]; - P::add_if_nonbonded(pi,pj); - } + for (unsigned int i=0; i< get_particles().size(); ++i) { + Particle *p= get_particles()[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() << " bipartite nonbonded pairs" + << " for " << P::get_particles().size() << " and " + << get_particles().size() << " sized sets" << std::endl); } void BipartiteNonbondedListScoreState::set_particles(const Particles &ps0, const Particles &ps1) { - ps_=ps0; - P::ps_=ps1; - Particles aps(ps0); - aps.insert(aps.end(), ps1.begin(), ps1.end()); - - P::audit_particles(aps); - P::propagate_set_particles(aps); - - rescan(); + IMP_LOG(VERBOSE, "Setting bipartite particles " << ps0.size() + << " and " << ps1.size() << std::endl); + mc_->set_particles(ps1); + P::set_particles(ps0); + P::audit_particles(ps1); + IMP_assert(ps0.size() == P::get_particles().size(), + "Where did they go?"); + IMP_assert(ps1.size() == get_particles().size(), + "Where did we go?"); } void BipartiteNonbondedListScoreState::update() { - P::propagate_update(); - rescan(); + P::update(); + if (mc_->get_max() > P::target_voxel_side_) { + P::invalidate_nbl(); + } } 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 354) +++ 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)); }