Index: kernel/pyext/IMP.i =================================================================== --- kernel/pyext/IMP.i (revision 393) +++ kernel/pyext/IMP.i (working copy) @@ -147,6 +147,7 @@ %include "IMP/score_states/MaxChangeScoreState.h" %include "IMP/score_states/NonbondedListScoreState.h" %include "IMP/score_states/AllNonbondedListScoreState.h" +%include "IMP/score_states/AllSphereNonbondedListScoreState.h" %include "IMP/score_states/BipartiteNonbondedListScoreState.h" %include "IMP/score_states/BondDecoratorListScoreState.h" Index: kernel/src/internal/ParticleGrid.cpp =================================================================== --- kernel/src/internal/ParticleGrid.cpp (revision 393) +++ kernel/src/internal/ParticleGrid.cpp (working copy) @@ -24,9 +24,10 @@ fks.push_back(FloatKey("y")); fks.push_back(FloatKey("z")); mc_= std::auto_ptr(new MaxChangeScoreState(fks)); - } +ParticleGrid::ParticleGrid(): target_voxel_side_(0), grid_valid_(0){} + void ParticleGrid::build_grid() { IMP_LOG(VERBOSE, "Creating nonbonded grid..." << std::flush); Index: kernel/src/score_states/AllSphereNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/AllSphereNonbondedListScoreState.cpp (revision 0) +++ kernel/src/score_states/AllSphereNonbondedListScoreState.cpp (revision 0) @@ -0,0 +1,178 @@ +#include "IMP/score_states/AllSphereNonbondedListScoreState.h" +#include "IMP/decorators/XYZDecorator.h" + +namespace IMP +{ + +AllSphereNonbondedListScoreState +::AllSphereNonbondedListScoreState(const Particles &ps, FloatKey radius): + rk_(radius){ + set_particles(ps); +} + +AllSphereNonbondedListScoreState +::~AllSphereNonbondedListScoreState() { + cleanup(bins_); + cleanup(fixed_bins_); +} + +void AllSphereNonbondedListScoreState::set_particles(const Particles &ps) { + NonbondedListScoreState::clear_nbl(); + // split into moving and fixed and call repartition twice + Particles moving, fixed; + for (unsigned int i=0; i< ps.size(); ++i) { + XYZDecorator d= XYZDecorator::cast(ps[i]); + if (d.get_coordinates_are_optimized()) { + moving.push_back(ps[i]); + } else { + fixed.push_back(ps[i]); + } + } + cleanup(bins_); + cleanup(fixed_bins_); + repartition_points(moving, bins_); + repartition_points(fixed, fixed_bins_); + for (unsigned int i=0; i< fixed_bins_.size(); ++i) { + fixed_bins_[i].grid->update(); + } +} + +void AllSphereNonbondedListScoreState::repartition_points(const Particles &ps, + std::vector &out){ + cleanup(out); + if (ps.empty()) return; + float minr=std::numeric_limits::max(), maxr=0; + for (unsigned int i=0; i< ps.size(); ++i) { + float r= ps[i]->get_value(rk_); + if ( r > maxr) maxr=r; + if ( r > 0 && r < minr) minr=r; + } + float curr=minr; + Floats cuts; + do { + cuts.push_back(curr); + curr *= 2; + } while (curr < maxr); + cuts.push_back(maxr); + + std::vector ops(cuts.size()); + for (unsigned int i=0; i< ps.size(); ++i) { + float r= ps[i]->get_value(rk_); + for (unsigned int j=0; ; ++j) { + if (cuts[j] >= r) { + ops[j].push_back(ps[i]); + break; + } + } + } + + for (unsigned int i=0; i< cuts.size(); ++i) { + if (ops[i].empty()) continue; + out.push_back(Bin()); + out.back().rmax= cuts[i]; + internal::ParticleGrid *pg + = new internal::ParticleGrid(side_from_r(out.back().rmax)); + out.back().grid= pg; + out.back().grid->add_particles(ops[i]); + } + IMP_LOG(VERBOSE, "Created " << out.size() << " grids" << std::endl); + for (unsigned int i=0; i< out.size(); ++i) { + IMP_LOG(VERBOSE, out[i].rmax + << ": " << out[i].grid->get_particles().size() << std::endl); + } + +#ifndef NDEBUG + Particles all; + for (unsigned int i=0; i< out.size(); ++i) { + all.insert(all.end(), out[i].grid->get_particles().begin(), + out[i].grid->get_particles().end()); + } + std::sort(all.begin(), all.end()); + Particles::iterator it = std::unique(all.begin(), all.end()); + IMP_assert(it == all.end(), "Duplicate particles " << all.size()); + IMP_assert(all.size() == ps.size(), "Wrong number of particles at end " + << all.size() << " " << ps.size()); +#endif +} + +void AllSphereNonbondedListScoreState::cleanup(std::vector &bins) { + for (unsigned int i=0; i< bins.size(); ++i) { + delete bins[i].grid; + } + bins.clear(); +} + +void AllSphereNonbondedListScoreState::update() { + bool bad=false; + for (unsigned int i=0; i< bins_.size(); ++i) { + if (bins_[i].grid->update()) bad=true; + } + if (bad) { + IMP_LOG(VERBOSE, "Destroying nbl in Sphere list"<< std::endl); + NonbondedListScoreState::clear_nbl(); + } +} + +void AllSphereNonbondedListScoreState::generate_nbl(const Bin &particle_bin, + const Bin &grid_bin, + float cut) { + IMP_CHECK_OBJECT(particle_bin.grid); + IMP_CHECK_OBJECT(grid_bin.grid); + IMP_LOG(VERBOSE, "Generate nbl for " << particle_bin.rmax + << " and " << grid_bin.rmax << std::endl); + for (unsigned int k=0; k< particle_bin.grid->get_particles().size(); ++k) { + Particle *p= particle_bin.grid->get_particles()[k]; + NonbondedListScoreState::AddToNBL f(this, p); + XYZDecorator d= XYZDecorator::cast(p); + internal::ParticleGrid::VirtualIndex index + = grid_bin.grid->get_virtual_index(Vector3D(d.get_x(), + d.get_y(), + d.get_z())); + IMP_LOG(VERBOSE, "Searching for " << p->get_index() + << " from " << index + << " of bin " << grid_bin.rmax << std::endl); + grid_bin.grid->apply_to_nearby(f, index, + cut+2*particle_bin.rmax + grid_bin.rmax, + false); + } + +} + +void AllSphereNonbondedListScoreState::rebuild_nbl(Float cut) { + IMP_LOG(VERBOSE, "Rebuilding NBL with " << bins_.size() << " dynamic and " + << fixed_bins_.size() << " fixed" << std::endl); + for (unsigned int i=0; i< bins_.size(); ++i) { + for (unsigned int j=i+1; j< bins_.size(); ++j) { + generate_nbl(bins_[i], bins_[j], cut); + } + + for (unsigned int j=0; j< fixed_bins_.size(); ++j) { + generate_nbl(bins_[i], fixed_bins_[j], cut); + } + + // same code as AllNonbonded. Would be nice to consolidate + internal::ParticleGrid::Index last_index; + for (internal::ParticleGrid::ParticleVoxelIterator it + = bins_[i].grid->particle_voxels_begin(); + it != bins_[i].grid->particle_voxels_end(); ++it) { + IMP_LOG(VERBOSE, "Searching with particle " << it->first->get_index() + << std::endl); + NonbondedListScoreState::AddToNBL f(this, it->first); + bins_[i].grid->apply_to_nearby(f, it->second, + cut+3*bins_[i].rmax, + true); + if (it->second != last_index) { + IMP_LOG(VERBOSE, "Searching in " << it->second + << std::endl); + bins_[i].grid->apply_to_cell_pairs(f, it->second); + } + last_index=it->second; + } + } +} + + +void AllSphereNonbondedListScoreState::show(std::ostream &out) const { + out << "AllSphereNonbondedListScoreState" << std::endl; +} +} Index: kernel/src/score_states/BipartiteNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/BipartiteNonbondedListScoreState.cpp (revision 393) +++ kernel/src/score_states/BipartiteNonbondedListScoreState.cpp (working copy) @@ -49,7 +49,8 @@ XYZDecorator d= XYZDecorator::cast(p); internal::ParticleGrid::VirtualIndex index = grid_.get_virtual_index(Vector3D(d.get_x(), d.get_y(), d.get_z())); - //IMP_LOG(VERBOSE, "Index is " << index << std::endl); + IMP_LOG(VERBOSE, "Searching for " << p->get_index() + << " from cell " << index << std::endl); grid_.apply_to_nearby(f, index, cut, false); } IMP_LOG(VERBOSE, "Found " << P::size_nbl() Index: kernel/include/IMP/internal/Grid3D.h =================================================================== --- kernel/include/IMP/internal/Grid3D.h (revision 393) +++ kernel/include/IMP/internal/Grid3D.h (working copy) @@ -180,32 +180,28 @@ private: std::vector data_; int d_[3]; - Vector3D min_, max_; + Vector3D min_; 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"); + IMP_assert(ii < data_.size(), "Invalid grid index " + << i[0] << " " << i[1] << " " << i[2] + << ": " << d_[0] << " " << d_[1] << " " << d_[2]); return ii; } template IndexType get_index_t(Vector3D pt) const { int index[3]; for (unsigned int i=0; i< 3; ++i ) { + IMP_assert(d_[i] != 0, "Invalid grid in Index"); float d= pt.get_component(i)- min_.get_component(i); index[i]= static_cast(std::floor(d/edge_size_[i])); + IMP_assert(std::abs(index[i]) < 200000000, + "Something is probably wrong " << d + << " " << pt.get_component(i) + << " " << min_.get_component(i) + << " " << edge_size_[i]); } return IndexType(index[0], index[1], index[2]); } @@ -249,14 +245,30 @@ \param[in] def The default value for the voxels */ Grid3D(int xd, int yd, int zd, - Vector3D minc, Vector3D maxc, + Vector3D minc, Vector3D max, VoxelData def): data_(xd*yd*zd, def), - min_(minc), - max_(maxc) { + min_(minc) { d_[0]=xd; d_[1]=yd; d_[2]=zd; - update_sizes(); + /** + \todo Grid3D::update_sizes is a mess of hacks. Figure out a clean + implementation. + */ + for (unsigned int i=0; i< 3; ++i) { + // hack to try to handle roundoff errors + // I would like to find something more reliable + if (d_[i]==1) { + // make sure that the total grid size is not vanishing + // this is probably not hte right thing to do + edge_size_[i]= std::max(1.05*(max.get_component(i) + - min_.get_component(i))/d_[i], + 1.0); + } else { + edge_size_[i]= 1.05*(max.get_component(i) + - min_.get_component(i))/d_[i]; + } + } } //! Initialize the grid @@ -269,16 +281,22 @@ Vector3D minc, Vector3D maxc, VoxelData def) { min_=minc; - float mx[3]; for (unsigned int i=0; i< 3; ++i ) { + IMP_assert(minc.get_component(i) <= maxc.get_component(i), + "Min must not be larger than max"); d_[i]= std::max(static_cast(std::ceil((maxc.get_component(i) - - minc.get_component(i))/ side)), + - 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(); + for (unsigned int i=0; i< 3; ++i) { + edge_size_[i]=side; + // hack to make sure the grid is big enough. + while (edge_size_[i]*d_[i] < maxc.get_component(i)) { + edge_size_[i]*= 1.05; + } + } } //! An empty grid. @@ -288,26 +306,16 @@ 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_; + Vector3D get_max() const { + return Vector3D(d_[0]*edge_size_[0], + d_[1]*edge_size_[1], + d_[2]*edge_size_[2]); } //! Return the number of voxels in a certain direction @@ -318,11 +326,12 @@ //! Return the index of the voxel containing the point. Index get_index(Vector3D pt) const { + VirtualIndex v= get_virtual_index(pt); 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(); + if (v[i] < 0) return Index(); + else if (v[i] >= d_[i]) return Index(); } - return get_index_t(pt); + return Index(v[0], v[1], v[2]); } //! Return the index that would contain the voxel if the grid extended there Index: kernel/include/IMP/internal/ParticleGrid.h =================================================================== --- kernel/include/IMP/internal/ParticleGrid.h (revision 393) +++ kernel/include/IMP/internal/ParticleGrid.h (working copy) @@ -9,18 +9,17 @@ #define __IMP_PARTICLE_GRID_H #include "Grid3D.h" +#include "../score_states/MaxChangeScoreState.h" #include "../base_types.h" namespace IMP { -class MaxChangeScoreState; - namespace internal { /** \internal */ -class ParticleGrid + class ParticleGrid: public Object { typedef internal::Grid3D Grid; Grid grid_; @@ -31,6 +30,7 @@ void build_grid(); void audit_particles(const Particles &ps) const; public: + ParticleGrid(); //! suggested grid edge size. ParticleGrid(Float sz); @@ -38,12 +38,13 @@ void add_particles(const Particles &ps); void clear_particles(); + const Particles& get_particles() const {return mc_->get_particles();} bool update(); typedef Grid::VirtualIndex VirtualIndex; typedef Grid::Index Index; Grid::VirtualIndex get_virtual_index(Vector3D pt) const { - return grid_.get_index(pt); + return grid_.get_virtual_index(pt); } //! Apply the function F to all particles near to the center @@ -59,7 +60,7 @@ template void apply_to_nearby(F f, const Grid::VirtualIndex ¢er, float cut, - bool skip_lower) { + bool skip_lower) const { Grid::VirtualIndex lc, uc; if ( cut > target_voxel_side_*1000 ) { @@ -99,7 +100,7 @@ } template - void apply_to_cell_pairs(F f, const Grid::Index ¢er) { + void apply_to_cell_pairs(F f, const Grid::Index ¢er) const { const Particles &ps= grid_.get_voxel(center); for (unsigned int i=0; i< ps.size(); ++i) { for (unsigned int j=0; j< i; ++j) { Index: kernel/include/IMP/score_states/AllSphereNonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/AllSphereNonbondedListScoreState.h (revision 0) +++ kernel/include/IMP/score_states/AllSphereNonbondedListScoreState.h (revision 0) @@ -0,0 +1,74 @@ +/** + * \file AllSphereNonbondedListScoreState.h + * \brief Allow iteration through pairs of a set of spheres. + * + * Copyright 2007-8 Sali Lab. All rights reserved. + */ + +#ifndef __IMP_NONBONDED_SPHERE_LIST_SCORE_STATE_H +#define __IMP_NONBONDED_SPHERE_LIST_SCORE_STATE_H + +#include "NonbondedListScoreState.h" +#include "../internal/ParticleGrid.h" +#include "MaxChangeScoreState.h" + +#include +#include + +namespace IMP +{ + +//! This class maintains a list of non-bonded pairs of spheres +/** \note The radius is currently assumed not to change. This could + be fixed later. + \note Points whose coordinates are not optimized are assumed to + stay that way (and vice + versa, although that direction doesn't matter so much). + \note The structure is slightly dumb about rebuilding and will + rebuild the whole list of any of the grids become invalidated. + This could be improved as each piece is computed separately (so + they could be cached). + + \ingroup restraint + */ +class IMPDLLEXPORT AllSphereNonbondedListScoreState: + public NonbondedListScoreState +{ + protected: + //! \internal + struct Bin { + internal::ParticleGrid *grid; + Float rmax; + Bin(): grid(NULL), rmax(-1){} + Bin(const Bin &o): grid(o.grid), rmax(o.rmax){} + }; + std::vector bins_; + std::vector fixed_bins_; + FloatKey rk_; + + //! \internal + void rebuild_nbl(Float cut); + + void repartition_points(const Particles &ps, std::vector &out); + + float side_from_r(float r) const {return r*1.6;} + + void generate_nbl(const Bin &particle_bin, const Bin &grid_bin, float cut); + + void cleanup(std::vector &bins); + +public: + /** + \param[in] ps A list of particles to use. + \param[in] radius The key to use to get the radius + */ + AllSphereNonbondedListScoreState(const Particles &ps, FloatKey radius); + ~AllSphereNonbondedListScoreState(); + IMP_SCORE_STATE("0.5", "Daniel Russel"); + + void set_particles(const Particles &ps); +}; + +} // namespace IMP + +#endif /* __IMP_NONBONDED_LIST_SCORE_STATE_H */ Index: kernel/include/IMP/score_states/NonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/NonbondedListScoreState.h (revision 393) +++ kernel/include/IMP/score_states/NonbondedListScoreState.h (working copy) @@ -73,8 +73,8 @@ } } if (!found) { - /*IMP_LOG(VERBOSE, "Found pair " << a->get_index() - << " " << b->get_index() << std::endl);*/ + IMP_LOG(VERBOSE, "Found pair " << a->get_index() + << " " << b->get_index() << std::endl); nbl_.push_back(std::make_pair(a, b)); } } Index: kernel/test/states/test_nonbonded_list.py =================================================================== --- kernel/test/states/test_nonbonded_list.py (revision 393) +++ kernel/test/states/test_nonbonded_list.py (working copy) @@ -46,7 +46,6 @@ def test_bl(self): """Test the bonded list""" - return m= IMP.Model() bds=[] for i in range(0,10): @@ -123,6 +122,69 @@ m.add_restraint(r) score= m.evaluate(False) self.assertEqual(score, 25, "Wrong score") + def test_spheres(self): + """Test the nonbonded list of spheres (num pairs)""" + m= IMP.Model() + rk= IMP.FloatKey("radius") + for i in range(0,10): + p= IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(random.uniform(0,10)) + d.set_y(random.uniform(0,10)) + d.set_z(random.uniform(0,10)) + p.add_attribute(rk, i, False) + d.set_coordinates_are_optimized(True) + s= IMP.AllSphereNonbondedListScoreState(m.get_particles(), rk) + m.add_score_state(s) + o= OnePair() + r= IMP.NonbondedRestraint(s, o, 15) + m.add_restraint(r) + score= m.evaluate(False) + self.assertEqual(score, 45, "Wrong score") + def test_spheres(self): + """Test the nonbonded list of spheres (collision detection)""" + m= IMP.Model() + rk= IMP.FloatKey("radius") + for i in range(0,10): + p= IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(random.uniform(0,10)) + d.set_y(random.uniform(0,10)) + d.set_z(random.uniform(0,10)) + p.add_attribute(rk, random.uniform(0,1000), False) + d.set_coordinates_are_optimized(True) + s= IMP.AllSphereNonbondedListScoreState(m.get_particles(), rk) + m.add_score_state(s) + sd= IMP.SphereDistancePairScore(IMP.HarmonicLowerBound(0,1), + rk) + r= IMP.NonbondedRestraint(s, sd, 1) + m.add_restraint(r) + score= m.evaluate(False) + opt= IMP.ConjugateGradients() + opt.set_model(m) + score =opt.optimize(10000) + print score + for p in m.get_particles(): + dp= IMP.XYZDecorator.cast(p) + print ".sphere "+str(dp.get_x()) + " " + str(dp.get_y())\ + + " " + str(dp.get_z()) + " " +str( p.get_value(rk)) + for p in m.get_particles(): + p.show() + dp= IMP.XYZDecorator.cast(p) + for q in m.get_particles(): + dq= IMP.XYZDecorator.cast(q) + if (p.get_index() != q.get_index()): + d = IMP.distance(dp,dq) + rd= p.get_value(rk) + q.get_value(rk) + if (rd > d): + p.show() + q.show() + print d + print rd + self.assert_(rd <= d, "Some sphere are not repelled") + if __name__ == '__main__': unittest.main()