Here are some changes to all spheres nonbonded list to fix a bug that Frido found as well as make things faster. The nonbonded lists now have checks built into them so asserts should fail if anything goes wrong. Note that these checks are very expensive for large numbers of particles (quadratic) and so any real runs with IMP should use builds with NDEBUG defined. We may want to define macros like IMP_CHECK_EXPENSIVE to provide finer grain control over debugging checks.

Index: kernel/include/IMP/score_states/NonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/NonbondedListScoreState.h (revision 456) +++ kernel/include/IMP/score_states/NonbondedListScoreState.h (working copy) @@ -27,11 +27,13 @@ class IMPDLLEXPORT NonbondedListScoreState: public ScoreState { private: +protected: + // made protected for debuggin code, do not use otherwise typedef std::vector<std::pair<Particle*, Particle*> > NBL; NBL nbl_; float last_cutoff_; -protected:

+ unsigned int size_nbl() const {return nbl_.size();}

//! rebuild the nonbonded list @@ -60,20 +62,23 @@ } };

+ bool are_bonded(Particle *a, Particle *b) const { + for (BondedListConstIterator bli= bonded_lists_begin(); + bli != bonded_lists_end(); ++bli) { + if ((*bli)->are_bonded(a, b)) { + return true; + } + } + return false; + } + //! tell the bonded lists what particles to pay attention to void propagate_particles(const Particles &ps);

void add_to_nbl(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) { - if ((*bli)->are_bonded(a, b)) { - found = true; - break; - } - } - if (!found) { + + if (!are_bonded(a,b)) { IMP_LOG(VERBOSE, "Found pair " << a->get_index() << " " << b->get_index() << std::endl); nbl_.push_back(std::make_pair(a, b)); @@ -90,6 +95,7 @@ // kind of evil hack to make the names better // perhaps the macro should be made more flexible typedef BondedListScoreStateIterator BondedListIterator; + typedef BondedListScoreStateConstIterator BondedListConstIterator;

IMP_SCORE_STATE(internal::kernel_version_info)

Index: kernel/test/states/test_nonbonded_list.py =================================================================== --- kernel/test/states/test_nonbonded_list.py (revision 456) +++ kernel/test/states/test_nonbonded_list.py (working copy) @@ -27,7 +27,7 @@ def test_it(self): """Test the nonbonded list and restraint which uses it""" m= IMP.Model() - for i in range(0,10): + for i in range(0,100): p= IMP.Particle() m.add_particle(p) d=IMP.XYZDecorator.create(p) @@ -39,13 +39,13 @@ r= IMP.NonbondedRestraint(o, s, 15) m.add_restraint(r) score= m.evaluate(False) - self.assertEqual(score, 45, "Wrong score") + self.assertEqual(score, 4950, "Wrong score")

def test_bl(self): """Test the bonded list""" m= IMP.Model() bds=[] - for i in range(0,10): + for i in range(0,100): p= IMP.Particle() m.add_particle(p) d=IMP.XYZDecorator.create(p) @@ -56,7 +56,7 @@ pts=IMP.Particles() for p in m.get_particles(): pts.append(p) - for i in range(1,10): + for i in range(1,100): IMP.custom_bond(bds[i-1], bds[i], 1, .1) s= IMP.AllNonbondedListScoreState(pts, 1) b= IMP.BondDecoratorListScoreState(pts) @@ -70,7 +70,8 @@ m.add_restraint(r) m.add_restraint(br) score= m.evaluate( False ) - self.assertEqual(score, 900+45-9, "Wrong score") + print "Score with bonds is " + str(score) + self.assertEqual(score, 9900+4950-99, "Wrong score")

def test_distfilt(self): """Test filtering based on distance in nonbonded list""" @@ -120,7 +121,7 @@ """Test the nonbonded list of spheres (num pairs)""" m= IMP.Model() rk= IMP.FloatKey("radius") - for i in range(0,10): + for i in range(0,100): p= IMP.Particle() m.add_particle(p) d=IMP.XYZDecorator.create(p) @@ -134,18 +135,149 @@ r= IMP.NonbondedRestraint(o, s, 15) m.add_restraint(r) score= m.evaluate(False) - self.assertEqual(score, 45, "Wrong score") + self.assertEqual(score, 4950, "Wrong score") + def test_frido_spheres(self): + """Test the nonbonded list with frido's spheres""" + m= IMP.Model() + rk= IMP.FloatKey("radius") + print "Frido begin" + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(87.0621490479) + d.set_y(140.299957275) + d.set_z(76.7119979858) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 20.1244087219, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(80.805770874) + d.set_y(99.9667434692) + d.set_z(66.9167098999) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 19.8160457611, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(112.992515564) + d.set_y(119.602111816) + d.set_z(53.6557235718) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 20.3428726196, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(118.784294128) + d.set_y(86.842376709) + d.set_z(65.2264938354) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 20.1794681549, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(132.930664062) + d.set_y(85.9250793457) + d.set_z(62.0034713745) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 19.730260849, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(97.6803894043) + d.set_y(77.0734939575) + d.set_z(45.9188423157) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 20.0133743286, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(74.9787445068) + d.set_y(103.48789978) + d.set_z(65.8464813232) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 20.1519756317, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(74.5077972412) + d.set_y(114.997116089) + d.set_z(103.185188293) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 18.6368904114, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(110.806472778) + d.set_y(100.937240601) + d.set_z(91.8201828003) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 19.2294254303, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(125.188072205) + d.set_y(135.360610962) + d.set_z(84.0617752075) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 18.9221935272, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(122.894897461) + d.set_y(114.53062439) + d.set_z(92.5285186768) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 18.9533672333, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(124.656105042) + d.set_y(88.4534988403) + d.set_z(99.5341949463) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 18.8280506134, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(109.428604126) + d.set_y(60.1015129089) + d.set_z(79.1919250488) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 19.1991424561, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(79.7605895996) + d.set_y(78.1874771118) + d.set_z(95.7365951538) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 19.3197059631, False) + for p in m.get_particles(): + d= IMP.XYZDecorator.cast(p) + print ".sphere "+str(d.get_x())+ " " + str(d.get_y())\ + + " " + str(d.get_z()) + " " + str(p.get_value(rk)) + + s= IMP.AllSphereNonbondedListScoreState(m.get_particles(), rk) + m.add_score_state(s) + sd= IMP.SphereDistancePairScore(IMP.HarmonicLowerBound(0,1), + rk) + r= IMP.NonbondedRestraint(sd, s, 1) + m.add_restraint(r) + score= m.evaluate(False) + print "Frido score is" + print 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): + for i in range(0,100): p= IMP.Particle() m.add_particle(p) d=IMP.XYZDecorator.create(p) d.randomize_in_box(IMP.Vector3D(0,0,0), - IMP.Vector3D(10,10,10)); - p.add_attribute(rk, random.uniform(0,1000), False) + IMP.Vector3D(20,20,20)); + p.add_attribute(rk, random.uniform(0,100), False) d.set_coordinates_are_optimized(True) s= IMP.AllSphereNonbondedListScoreState(m.get_particles(), rk) m.add_score_state(s) @@ -156,6 +288,7 @@ score= m.evaluate(False) opt= IMP.ConjugateGradients() opt.set_model(m) + IMP.set_log_level(IMP.TERSE) score =opt.optimize(10000) print score for p in m.get_particles(): Index: kernel/include/IMP/internal/ParticleGrid.h =================================================================== --- kernel/include/IMP/internal/ParticleGrid.h (revision 456) +++ kernel/include/IMP/internal/ParticleGrid.h (working copy) @@ -42,6 +42,8 @@ const Particles& get_particles() const {return mc_->get_particles();} bool update();

+ void show(std::ostream &out) const; + typedef Grid::VirtualIndex VirtualIndex; typedef Grid::Index Index; Grid::VirtualIndex get_virtual_index(Vector3D pt) const { @@ -148,6 +150,8 @@ if (curp_== grid_->get_voxel(*cvoxel_).end()) { ++cvoxel_; find_voxel(); + } else { + temp_= std::make_pair(*curp_, *cvoxel_); } }

@@ -178,6 +182,8 @@ } };

+IMP_OUTPUT_OPERATOR(ParticleGrid); + } // namespace internal

} // namespace IMP Index: kernel/src/score_states/AllSphereNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/AllSphereNonbondedListScoreState.cpp (revision 456) +++ kernel/src/score_states/AllSphereNonbondedListScoreState.cpp (working copy) @@ -13,6 +13,8 @@ namespace IMP {

+static unsigned int min_grid_size=20; + AllSphereNonbondedListScoreState ::AllSphereNonbondedListScoreState(const Particles &ps, FloatKey radius): rk_(radius) @@ -59,29 +61,40 @@ if ( r > maxr) maxr=r; if ( r > 0 && r < minr) minr=r; } - float curr=minr; + float curr=minr*2; Floats cuts; do { cuts.push_back(curr); curr *= 2; } while (curr < maxr); - cuts.push_back(maxr); + cuts.push_back(2*maxr);

std::vector<Particles> 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) { + IMP_assert(j< cuts.size(), "Internal error in ASNBLSS"); + if (cuts[j] > r) { ops[j].push_back(ps[i]); break; } } } - + // consolidate + for (unsigned int i=1; i< ops.size(); ++i) { + if (ops[i-1].size() + ops[i].size() < min_grid_size) { + ops[i].insert(ops[i].end(), ops[i-1].begin(), ops[i-1].end()); + ops[i-1].clear(); + } + } for (unsigned int i=0; i< cuts.size(); ++i) { if (ops[i].empty()) continue; out.push_back(Bin()); - out.back().rmax= cuts[i]; + float rmax=0; + for (unsigned int j=0; j< ops[i].size(); ++j) { + rmax= std::max(rmax, ops[i][j]->get_value(rk_)); + } + out.back().rmax= rmax; internal::ParticleGrid *pg = new internal::ParticleGrid(side_from_r(out.back().rmax)); out.back().grid= pg; @@ -90,7 +103,7 @@ 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); + << ": " << *out[i].grid << std::endl); }

#ifndef NDEBUG @@ -121,6 +134,7 @@ bool bad=false; for (unsigned int i=0; i< bins_.size(); ++i) { if (bins_[i].grid->update()) bad=true; + IMP_LOG(VERBOSE, bins_[i].rmax << std::endl << *bins_[i].grid << std::endl); } if (bad) { IMP_LOG(VERBOSE, "Destroying nbl in Sphere list"<< std::endl); @@ -186,6 +200,37 @@ last_index=it->second; } } + +#ifndef NDEBUG + Particles ps; + for (unsigned int i=0; i< bins_.size(); ++i) { + if (bins_[i].grid) { + ps.insert(ps.end(), bins_[i].grid->get_particles().begin(), + bins_[i].grid->get_particles().end()); + } + } + for (unsigned int i=0; i< ps.size(); ++i) { + XYZDecorator di= XYZDecorator::cast(ps[i]); + for (unsigned int j=0; j< i; ++j) { + XYZDecorator dj= XYZDecorator::cast(ps[j]); + if (distance(di, dj) - ps[i]->get_value(rk_) - ps[j]->get_value(rk_) + <= cut && !are_bonded(ps[i], ps[j])) { + bool found=false; + for (NonbondedIterator nit= nbl_.begin(); + nit != nbl_.end(); ++nit) { + if (nit->first == ps[i] && nit->second == ps[j] + || nit->first == ps[j] && nit->second == ps[i]) { + IMP_assert(!found, "Entry is in list twice"); + found=true; + } + } + IMP_assert(found, "Nonbonded list is missing " + << ps[i]->get_index() << " " << di + << " and " << ps[j]->get_index() << " " << dj << std::endl); + } + } + } +#endif }

Index: kernel/src/score_states/AllNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/AllNonbondedListScoreState.cpp (revision 456) +++ kernel/src/score_states/AllNonbondedListScoreState.cpp (working copy) @@ -41,6 +41,30 @@ IMP_LOG(VERBOSE, "Found " << NonbondedListScoreState::size_nbl() << " nonbonded pairs" << std::endl); + +#ifndef NDEBUG + Particles ps= grid_.get_particles(); + for (unsigned int i=0; i< ps.size(); ++i) { + XYZDecorator di= XYZDecorator::cast(ps[i]); + for (unsigned int j=0; j< i; ++j) { + XYZDecorator dj= XYZDecorator::cast(ps[j]); + if (distance(di, dj) <= cut && !are_bonded(ps[i], ps[j]) ) { + bool found=false; + for (NonbondedIterator nit= nbl_.begin(); + nit != nbl_.end(); ++nit) { + if (nit->first == ps[i] && nit->second == ps[j] + || nit->first == ps[j] && nit->second == ps[i]) { + IMP_assert(!found, "Entry is in list twice"); + found=true; + } + } + IMP_assert(found, "Nonbonded list is missing " + << ps[i]->get_index() << " " << di + << " and " << ps[j]->get_index() << " " << dj << std::endl); + } + } + } +#endif }

void AllNonbondedListScoreState::set_particles(const Particles &ps) @@ -56,6 +80,7 @@ if (grid_.update()) { NonbondedListScoreState::clear_nbl(); } + IMP_LOG(VERBOSE, grid_); }

void AllNonbondedListScoreState::show(std::ostream &out) const Index: kernel/src/internal/ParticleGrid.cpp =================================================================== --- kernel/src/internal/ParticleGrid.cpp (revision 456) +++ kernel/src/internal/ParticleGrid.cpp (working copy) @@ -108,6 +108,19 @@ } }

+ +void ParticleGrid::show(std::ostream &out) const { + for (Grid::IndexIterator it= grid_.all_indexes_begin(); + it != grid_.all_indexes_end(); ++it) { + out << *it << ": "; + //Grid::Index + for (unsigned int i=0; i< grid_.get_voxel(*it).size(); ++i) { + out << grid_.get_voxel(*it)[i]->get_index() << " "; + } + out << std::endl; + } +} + } // namespace internal

} // namespace IMP