Index: kernel/test/states/test_score_state.py =================================================================== --- kernel/test/states/test_score_state.py (revision 541) +++ kernel/test/states/test_score_state.py (working copy) @@ -27,10 +27,10 @@ IMP.ScoreState.__init__(self) self.log = log - def update(self): + def do_before_evaluate(self): self.log.append('update') - def after_evaluate(self, accum): + def do_after_evaluate(self, accum): self.log.append('after_evaluate') class TestScoreState(IMP.test.TestCase): Index: kernel/test/states/test_vrml_log.py =================================================================== --- kernel/test/states/test_vrml_log.py (revision 541) +++ kernel/test/states/test_vrml_log.py (working copy) @@ -6,7 +6,7 @@ def setUp(self): IMP.set_log_level(IMP.TERSE) - def _testit(self, rk, r,g,b, pref): + def _testit(self, rk,pref): """Test logging to a VRML file""" m= IMP.Model() o= IMP.SteepestDescent() @@ -23,15 +23,11 @@ p1= IMP.Particle() m.add_particle(p1) d1= IMP.XYZDecorator.create(p1) - p1.add_attribute(r, 1.0, False) - p1.add_attribute(g, 0.0, False) - p1.add_attribute(b, 0.0, False) d1.set_x(1) d1.set_y(1) d1.set_z(1) a= IMP.VRMLLogOptimizerState(nm, IMP.Particles([p0,p1])) - a.set_radius(rk) - a.set_color(r, g, b) + a.set_radius_key(rk) o.add_optimizer_state(a) a.update() @@ -40,17 +36,12 @@ def test_1(self): """Testing the VRML log""" self._testit(IMP.FloatKey("radius"), - IMP.FloatKey("red"), - IMP.FloatKey("green"), - IMP.FloatKey("blue"), "test1") + "test1") def test_2(self): """Testing the VRML log with new attribute names""" self._testit(IMP.FloatKey("another_radius"), - IMP.FloatKey("red5"), - IMP.FloatKey("green5"), - IMP.FloatKey("blue5"), - "test1") + "test1") def test_skip(self): """Test skipping steps in the VRML log""" m= IMP.Model() Index: kernel/test/states/test_nonbonded_list.py =================================================================== --- kernel/test/states/test_nonbonded_list.py (revision 541) +++ kernel/test/states/test_nonbonded_list.py (working copy) @@ -6,6 +6,14 @@ def __init__(self): IMP.PairScore.__init__(self) def evaluate(self, pa, pb, da): + if (not pa.get_is_active()): + print "Inactive particle" + pa.show() + pa.get_value(IMP.FloatKey("x")) + if (not pb.get_is_active()): + print "Inactive particle" + pb.show() + pb.get_value(IMP.FloatKey("x")) return 1 def get_version_info(self): return IMP.VersionInfo("Me", "0.5") @@ -24,58 +32,118 @@ class TestNBL(IMP.test.TestCase): def setUp(self): - IMP.set_log_level(IMP.VERBOSE) + IMP.set_log_level(IMP.TERSE) + self.rk= IMP.FloatKey("radius") - def test_it(self): - """Test the nonbonded list and restraint which uses it""" - m= IMP.Model() - for i in range(0,100): + def test_quadratic(self): + """Test quadratic NBL""" + ss='IMP.AllNonbondedListScoreState(md, self.rk, IMP.AllNonbondedListScoreState.QUADRATIC)' + #eval(ss) + self.do_test_all(ss) + + + def test_grid(self): + """Test grid NBL""" + #IMP.set_log_level(IMP.VERBOSE) + ss='IMP.AllNonbondedListScoreState(md, self.rk, IMP.AllNonbondedListScoreState.GRID)' + self.do_test_all(ss) + + def test_bbox(self): + """Test bbox NBL""" + try: + q= IMP.AllNonbondedListScoreState.BBOX + except: + print "No CGAL support" + else: + ss= 'IMP.AllNonbondedListScoreState(md, self.rk, IMP.AllNonbondedListScoreState.BBOX)' + self.do_test_all(ss) + + def test_default(self): + """Test default NBL""" + ss= 'IMP.AllNonbondedListScoreState(md, self.rk)' + self.do_test_all(ss) + + def test_bipartite_quadratic(self): + """Test bipartite quadratic NBL""" + ss='IMP.BipartiteNonbondedListScoreState(md, self.rk, IMP.BipartiteNonbondedListScoreState.QUADRATIC)' + self.do_test_all_bi(ss) + + def test_bipartite_bbox(self): + """Test bipartite bbox NBL""" + try: + q= IMP.BipartiteNonbondedListScoreState.BBOX + except: + print "No CGAL support" + else: + ss= 'IMP.BipartiteNonbondedListScoreState(md, self.rk, IMP.BipartiteNonbondedListScoreState.BBOX)' + self.do_test_all_bi(ss) + + + + def do_test_all(self, ss): + self.do_test_bl(ss) + self.do_test_distfilt(ss) + self.do_test_spheres(ss) + + def do_test_all_bi(self, ss): + self.do_test_bi(ss) + self.do_test_bi_update(ss) + + def make_spheres(self, m, num, lbv, ubv, minr, maxr): + ps=self.make_particles(m, num, lbv, ubv) + for p in ps: + p.add_attribute(self.rk, random.uniform(minr, maxr), False) + return ps + + def make_particles(self, m, num, lbv, ubv): + ps=IMP.Particles() + lb=IMP.Vector3D(lbv[0],lbv[1],lbv[2]) + ub=IMP.Vector3D(ubv[0],ubv[1],ubv[2]) + for i in range(0, num): p= IMP.Particle() m.add_particle(p) d=IMP.XYZDecorator.create(p) d.set_coordinates_are_optimized(True) - d.randomize_in_box(IMP.Vector3D(0,0,0), - IMP.Vector3D(10,10,10)); - s= IMP.AllNonbondedListScoreState(IMP.FloatKey(), - m.get_particles()) + d.randomize_in_box(lb, ub) + ps.append(p) + return ps + + def do_test_update(self, ss): + m= IMP.Model() + self.make_spheres(m, 20, [0,0,0], [10,10,10], .1, 1) + ds=[] + for p in m.get_particles(): + ds.append(IMP.XYZDecorator.cast(p)) + md=5 + s=eval(ss) + s.set_particles(m.get_particles()) m.add_score_state(s) o= OnePair() - r= IMP.NonbondedRestraint(o, s, 15) + r= IMP.NonbondedRestraint(o, s) m.add_restraint(r) - score= m.evaluate(False) - print score - self.assertEqual(score, 4950, "Wrong score") - p = m.get_particle(IMP.ParticleIndex(3)) - p.set_is_active(False) - score= m.evaluate(False) - print score - self.assertEqual(score, 4851, "Wrong score with removal") + # use the internal checks + for i in range(0,10): + score= m.evaluate(False) + for d in ds: + d.randomize_in_sphere(d.get_vector(), 2.0) - def test_bl(self): - """Test the bonded list""" + def do_test_bl(self, ss): + """Test the bond decorator list""" m= IMP.Model() bds=[] - for i in range(0,100): - p= IMP.Particle() - m.add_particle(p) - d=IMP.XYZDecorator.create(p) - d.set_coordinates_are_optimized(True) - d.set_x(0) - d.set_y(10) - d.set_z(10) + pts=self.make_spheres(m, 20, [0,10,10], [10,20,20], .1, 1) + for p in pts: bds.append(IMP.BondedDecorator.create(p)) - pts=IMP.Particles() - for p in m.get_particles(): - pts.append(p) - for i in range(1,100): + for i in range(1,len(pts)): IMP.custom_bond(bds[i-1], bds[i], 1, .1) - s= IMP.AllNonbondedListScoreState(IMP.FloatKey(), - pts) + md= 20 + s=eval(ss) + s.set_particles(pts) b= IMP.BondDecoratorListScoreState(pts) s.add_bonded_list(b) m.add_score_state(s) o= OnePair() - r= IMP.NonbondedRestraint(o, s, 10) + r= IMP.NonbondedRestraint(o, s) os=OneScore() print os.evaluate(6) br= IMP.BondDecoratorRestraint(os, b) @@ -83,39 +151,37 @@ m.add_restraint(br) score= m.evaluate( False ) print "Score with bonds is " + str(score) - self.assertEqual(score, 9900+4950-99, "Wrong score") + self.assertEqual(score, 190+1900-19, "Wrong score") - def test_distfilt(self): + def do_test_distfilt(self, ss): """Test filtering based on distance in nonbonded list""" - IMP.set_log_level(IMP.TERSE) + #IMP.set_log_level(IMP.TERSE) m= IMP.Model() ps=IMP.Particles() - for i in range(0,200): - p= IMP.Particle() - m.add_particle(p) - d=IMP.XYZDecorator.create(p) + ps= self.make_particles(m, 20, [0,0,0], [10,10,10]) + pts= self.make_particles(m, 20, [160,160,160], [170,170,170]) + for p in pts: ps.append(p) - d.set_coordinates_are_optimized(True) - if (i < 100): - d.randomize_in_box(IMP.Vector3D(0,0,0), - IMP.Vector3D(10,10,10)); - else: - d.randomize_in_box(IMP.Vector3D(160,160,160), - IMP.Vector3D(170,170,170)); - s= IMP.AllNonbondedListScoreState(IMP.FloatKey(), ps) + md=15 + s=eval(ss) + s.set_particles(ps) m.add_score_state(s) o= OnePair() - r= IMP.NonbondedRestraint(o, s, 15) + r= IMP.NonbondedRestraint(o, s) m.add_restraint(r) score= m.evaluate(False) - self.assertEqual(score, 9900, "Wrong score") + self.assertEqual(score, 2*190, "Wrong score") - ps[3].set_is_active(False) + m.remove_particle(IMP.ParticleIndex(3)) + self.assert_(not ps[3].get_is_active(), "Particle not inactive") ps=None score= m.evaluate(False) print score - self.assertEqual(score, 9801, "Wrong score with removal") + self.assertEqual(score, 171+190, "Wrong score with removal") + #for p in s.get_particles(): + # self.assert_(p.get_is_active(), "Inactive particle not removed") + p= IMP.Particle() m.add_particle(p) print "Index is " +str(p.get_index().get_index()) @@ -127,246 +193,83 @@ s.add_particles(nps) score= m.evaluate(False) print score - self.assertEqual(score, 9900, "Wrong score after insert") + self.assertEqual(score, 2*190, "Wrong score after insert") - def test_bi(self): + def do_test_bi(self, ss): """Test the bipartite nonbonded list and restraint which uses it""" m= IMP.Model() - ps0=IMP.Particles() - ps1=IMP.Particles() - for i in range(0,10): - p= IMP.Particle() - m.add_particle(p) - d=IMP.XYZDecorator.create(p) - d.set_coordinates_are_optimized(True) - d.randomize_in_box(IMP.Vector3D(0,0,0), - IMP.Vector3D(10,10,10)); - if (i < 5): - ps0.append(p) - else: - ps1.append(p) - s= IMP.QuadraticBipartiteNonbondedListScoreState(IMP.FloatKey(), - ps0, ps1) + ps0=self.make_spheres(m, 5, [0,0,0], [10,10,10], .1, 1) + ps1=self.make_spheres(m, 5, [0,0,0], [10,10,10], .1, 1) + md=15 + s=eval(ss) + #IMP.QuadraticBipartiteNonbondedListScoreState(IMP.FloatKey(), + # ps0, ps1) + s.set_particles(ps0, ps1) m.add_score_state(s) o= OnePair() - r= IMP.NonbondedRestraint(o, s, 15) + r= IMP.NonbondedRestraint(o, s) m.add_restraint(r) score= m.evaluate(False) + print score self.assertEqual(score, 25, "Wrong score") - def test_spheres2(self): - """Test the (quadratic and optimized) nonbonded lists of spheres (num pairs)""" + + def do_test_bi_update(self, ss): + """Test the bipartite nonbonded list and restraint which uses it""" m= IMP.Model() - rk= IMP.FloatKey("radius") - 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, i, False) - d.set_coordinates_are_optimized(True) - s= IMP.AllNonbondedListScoreState(rk, - m.get_particles()) - s2= IMP.QuadraticAllNonbondedListScoreState(rk, m.get_particles()) + ps0=self.make_spheres(m, 5, [0,0,0], [10,10,10], .1, 1) + ps1=self.make_spheres(m, 5, [0,0,0], [10,10,10], .1, 1) + md=5 + s=eval(ss) + #IMP.QuadraticBipartiteNonbondedListScoreState(IMP.FloatKey(), + # ps0, ps1) + s.set_particles(ps0, ps1) m.add_score_state(s) - m.add_score_state(s2) o= OnePair() - o2= OnePair() - r= IMP.NonbondedRestraint(o, s, 15) - r2= IMP.NonbondedRestraint(o2, s2, 15) + r= IMP.NonbondedRestraint(o, s) m.add_restraint(r) - m.add_restraint(r2) - score= m.evaluate(False) - self.assertEqual(score, 2*4950, "Wrong score") - def test_spheres3(self): - """Test the quadratic nonbonded list of spheres (num pairs)""" - m= IMP.Model() - rk= IMP.FloatKey("radius") - 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, i, False) - d.set_coordinates_are_optimized(True) - s= IMP.QuadraticAllNonbondedListScoreState(rk, m.get_particles()) - m.add_score_state(s) - o= OnePair() - r= IMP.NonbondedRestraint(o, s, 15) - m.add_restraint(r) - score= m.evaluate(False) - self.assertEqual(score, 4950, "Wrong score") - def test_frido_spheres(self): - """Test the nonbonded list with frido's spheres""" - # This uses the internal checks of the nonbonded list to - # verify correctness - 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)) + for i in range(0,20): + score= m.evaluate(False) + for p in ps0: + d= IMP.XYZDecorator.cast(p) + d.randomize_in_sphere(d.get_vector(), 1) + for p in ps1: + d= IMP.XYZDecorator.cast(p) + d.randomize_in_sphere(d.get_vector(), 1) - s= IMP.AllNonbondedListScoreState(rk, m.get_particles()) - 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): + + def do_test_spheres(self, ss): """Test the nonbonded list of spheres (collision detection)""" m= IMP.Model() - rk= IMP.FloatKey("radius") - 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(20,20,20)); - p.add_attribute(rk, random.uniform(0,100), False) - d.set_coordinates_are_optimized(True) - s= IMP.AllNonbondedListScoreState(rk, m.get_particles()) + ps= self.make_spheres(m, 20, IMP.Vector3D(0,0,0), + IMP.Vector3D(10,10,10), 0, 100) + md= 1 + s=eval(ss) + s.set_particles(m.get_particles()) + #IMP.AllNonbondedListScoreState(rk, m.get_particles()) m.add_score_state(s) sd= IMP.SphereDistancePairScore(IMP.HarmonicLowerBound(0,1), - rk) - r= IMP.NonbondedRestraint(sd, s, 1) + self.rk) + r= IMP.NonbondedRestraint(sd, s) m.add_restraint(r) score= m.evaluate(False) opt= IMP.ConjugateGradients() opt.set_model(m) - IMP.set_log_level(IMP.TERSE) + #IMP.set_log_level(IMP.TERSE) 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(self.rk)) for p in m.get_particles(): + #p.show() 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) + rd= p.get_value(self.rk) + q.get_value(self.rk) if (rd > d): p.show() q.show() @@ -375,6 +278,7 @@ # Allow a little extra, for imperfect optimization: self.assert_(rd <= d + 1e-3, "Some spheres are not repelled") + print "Memory usage is " + str(m.get_memory_usage()) if __name__ == '__main__': Index: kernel/test/states/test_gravity_center.py =================================================================== --- kernel/test/states/test_gravity_center.py (revision 541) +++ kernel/test/states/test_gravity_center.py (working copy) @@ -42,7 +42,7 @@ for weighted in (False, True): (ps, center, gc) = self._make_gravity_center(model, 10.0, 20.0, 30.0, weighted) - gc.update() + gc.update_position() # The gravity center's xyz should not be optimizable: self.assertEqual(center.get_is_optimized(xkey), False) self.assertEqual(center.get_is_optimized(ykey), False) Index: kernel/test/particles/test_clone.py =================================================================== --- kernel/test/particles/test_clone.py (revision 0) +++ kernel/test/particles/test_clone.py (revision 0) @@ -0,0 +1,42 @@ +import unittest +import IMP +import IMP.utils +import IMP.test + +class ParticleTests(IMP.test.TestCase): + """Test clone particles""" + + def test_clone_one(self): + """Check that cloned particle has same attributes""" + p= IMP.Particle() + m= IMP.Model() + m.add_particle(p) + + ik= IMP.IntKey("ik") + sk= IMP.StringKey("sk") + fk= IMP.FloatKey("fk") + pk= IMP.ParticleKey("pk") + + p.add_attribute(ik, 3) + p.add_attribute(fk, 4.0) + p.add_attribute(sk, "hello world") + p.add_attribute(pk, p) + + pc= p.clone() + self.assert_(pc.get_index() != p.get_index(), "Indexes should be different") + self.assertEqual(pc.get_float_attributes().size(), 1, + "Should only be one float key in particle") + self.assertEqual(pc.get_string_attributes().size(), 1, + "Should only be one string key in particle") + self.assertEqual(pc.get_particle_attributes().size(), 1, + "Should only be one particle key in particle") + self.assertEqual(pc.get_int_attributes().size(), 1, + "Should only be one int key in particle") + self.assertEqual(pc.get_value(ik), 3, "Int value not copied") + self.assertEqual(pc.get_value(fk), 4.0, "Float value not copied") + self.assertEqual(pc.get_value(sk), "hello world", + "String value not copied") + self.assertEqual(pc.get_value(pk).get_index(), p.get_index(), "Particle value not copied") + +if __name__ == '__main__': + unittest.main() Index: kernel/test/unary_functions/test_splines.py =================================================================== --- kernel/test/unary_functions/test_splines.py (revision 541) +++ kernel/test/unary_functions/test_splines.py (working copy) @@ -31,6 +31,8 @@ self.assertEqual(closed_spline.evaluate(25.0), 0.0) self.assertRaises(ValueError, open_spline.evaluate, 9.9) self.assertRaises(ValueError, open_spline.evaluate, 25.1) + self.check_unary_function_deriv(open_spline, 10, 20, .1); + self.check_unary_function_deriv(closed_spline, 10.0, 25.0, .1); def test_interpolate(self): """Test that spline-interpolated values are correct""" @@ -48,7 +50,7 @@ for i in range(10): floats.append(test_func(minrange + spline_spacing * i)[0]) spline = spline_func(floats, minrange, spline_spacing) - + self.check_unary_function_deriv(spline, minrange+.1, minrange+spline_spacing*9-.1, .1) # Now test the spline against the test function for intermediate points for i in range(30): val = minrange + test_spacing * i Index: kernel/test/unary_functions/test_wlc.py =================================================================== --- kernel/test/unary_functions/test_wlc.py (revision 0) +++ kernel/test/unary_functions/test_wlc.py (revision 0) @@ -0,0 +1,19 @@ +import unittest +import IMP +import IMP.test + + +class WLCTests(IMP.test.TestCase): + """Tests for WLC unary function""" + + def test_wlc(self): + """Test that the wlc values are sane""" + wlc= IMP.WormLikeChain(200, 3.4) + self.check_unary_function_deriv(wlc, 180, 250, .5) + self.check_unary_function_min(wlc, 0, 250, .5, 0) + self.check_unary_function_deriv(wlc, 0, 180, .5) + + self.assert_(wlc.evaluate_deriv(180)[1] >4.2, "Bad force values") + +if __name__ == '__main__': + unittest.main() Index: kernel/test/restraints/test_list.py =================================================================== --- kernel/test/restraints/test_list.py (revision 541) +++ kernel/test/restraints/test_list.py (working copy) @@ -55,6 +55,30 @@ e= m.evaluate(False) self.assertEqual(e,5, "Wrong distance in score") + def test_ss2(self): + """Test the enclosing sphere """ + m= IMP.Model() + p= IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(100) + d.set_y(1) + d.set_z(1) + d.set_coordinates_are_optimized(True) + v= IMP.Vector3D(5,5,5) + h= IMP.HarmonicUpperBound(10, 10) + s= IMP.DistanceToSingletonScore(h, v) + r= IMP.SingletonListRestraint(s, m.get_particles()) + m.add_restraint(r) + e= m.evaluate(False) + o= IMP.ConjugateGradients() + o.set_model(m) + o.optimize(100) + d= IMP.XYZDecorator.cast(m.get_particle(IMP.ParticleIndex(0))) + dist2 = (d.get_x()-5)**2+(d.get_y()-5)**2+(d.get_y()-5)**2 + print "Final" + d.show() + self.assert_(dist2 < 100, "Enclosing sphere not enclosing") if __name__ == '__main__': unittest.main() Index: kernel/include/IMP/singleton_scores/DistanceToSingletonScore.h =================================================================== --- kernel/include/IMP/singleton_scores/DistanceToSingletonScore.h (revision 541) +++ kernel/include/IMP/singleton_scores/DistanceToSingletonScore.h (working copy) @@ -31,6 +31,23 @@ virtual void show(std::ostream &out=std::cout) const; }; +// for some reason swig gets linking errors on linux with this function +#ifndef SWIG +namespace internal +{ + +//! An internal helper function for evaluating distance potentials +/** The function applies f to scale*(distance-offset). + */ +Float evaluate_distance_to_singleton_score(const Vector3D &v, + Particle *a, + DerivativeAccumulator *da, + UnaryFunction *f, + Float offset, + Float scale); +} +#endif + } // namespace IMP #endif /* __IMP_DISTANCE_TO_SINGLETON_SCORE_H */ Index: kernel/include/IMP/singleton_scores/SConscript =================================================================== --- kernel/include/IMP/singleton_scores/SConscript (revision 541) +++ kernel/include/IMP/singleton_scores/SConscript (working copy) @@ -1,9 +1,9 @@ +Import('env') import os.path -Import('env') - -files = ['DistanceToSingletonScore.h', 'AttributeSingletonScore.h'] - -# Install the include files: -includedir = os.path.join(env['includedir'], 'IMP', 'singleton_scores') +files=[ + 'AttributeSingletonScore.h', + 'DistanceToSingletonScore.h', + ] +includedir = os.path.join(env['includedir'], 'IMP', 'singleton_scores' ) inst = env.Install(includedir, files) env.Alias('install', inst) Index: kernel/include/IMP/Model.h =================================================================== --- kernel/include/IMP/Model.h (revision 541) +++ kernel/include/IMP/Model.h (working copy) @@ -33,6 +33,8 @@ class IMPDLLEXPORT Model: public internal::Object { friend class Restraint; + unsigned int iteration_; + public: Model(); ~Model(); @@ -58,6 +60,8 @@ VersionInfo get_version_info() const { return internal::kernel_version_info; } + + unsigned int get_memory_usage() const; }; Index: kernel/include/IMP/score_states/QuadraticNonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/QuadraticNonbondedListScoreState.h (revision 541) +++ kernel/include/IMP/score_states/QuadraticNonbondedListScoreState.h (working copy) @@ -1,51 +0,0 @@ -/** - * \file QuadraticNonbondedListScoreState.h - * \brief Allow iteration through pairs of a set of spheres. - * - * Copyright 2007-8 Sali Lab. All rights reserved. - */ - -#ifndef __IMP_QUADRATIC_NONBONDED_LIST_SCORE_STATE_H -#define __IMP_QUADRATIC_NONBONDED_LIST_SCORE_STATE_H - -#include "NonbondedListScoreState.h" -#include "../internal/kernel_version_info.h" -#include "MaxChangeScoreState.h" - -#include -#include - -namespace IMP -{ - -//! A base class for nonbonded lists which test all pairs of particles -/** This should not be used by end users. But since it needs to be in the - inheritance hierarchy, it should be in the IMP namespace. - */ -class IMPDLLEXPORT QuadraticNonbondedListScoreState: - public NonbondedListScoreState -{ - typedef NonbondedListScoreState P; - internal::ObjectPointer mc_; - float slack_; - protected: - void handle_nbl_pair( Particle *a, Particle *b, float d); - const Particles &get_particles() const { - return mc_->get_particles(); - } - void set_particles(const Particles &ps) { - P::clear_nbl(); - mc_->clear_particles(); - mc_->add_particles(ps); - } - - QuadraticNonbondedListScoreState(FloatKey radius); - ~QuadraticNonbondedListScoreState(); - - public: - void update(); -}; - -} // namespace IMP - -#endif /* __IMP_QUADRATIC_NONBONDED_LIST_SCORE_STATE_H */ Index: kernel/include/IMP/score_states/BondDecoratorListScoreState.h =================================================================== --- kernel/include/IMP/score_states/BondDecoratorListScoreState.h (revision 541) +++ kernel/include/IMP/score_states/BondDecoratorListScoreState.h (working copy) @@ -8,10 +8,11 @@ #ifndef __IMP_BOND_DECORATOR_LIST_SCORE_STATE_H #define __IMP_BOND_DECORATOR_LIST_SCORE_STATE_H -#include #include "BondedListScoreState.h" #include "../decorators/bond_decorators.h" +#include + namespace IMP { @@ -25,6 +26,7 @@ */ class IMPDLLEXPORT BondDecoratorListScoreState: public BondedListScoreState { + // This is only used to provide the iterators. Probably should be lazy. std::vector bonds_; Particles ps_; public: @@ -38,8 +40,6 @@ virtual bool are_bonded(Particle *a, Particle *b) const; - virtual void update(); - //! This iterates through the pairs of bonded particles /** \note update() must be called first for this to be valid. */ @@ -54,6 +54,15 @@ unsigned int number_of_bonds() const { return bonds_.size(); } + + unsigned int get_memory_usage() const { + return BondedListScoreState::get_memory_usage() + + bonds_.capacity() * sizeof(BondDecorator) + + sizeof(Particle*)* ps_.capacity(); + } + + protected: + virtual void do_before_evaluate(); }; } // namespace IMP Index: kernel/include/IMP/score_states/NonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/NonbondedListScoreState.h (revision 541) +++ kernel/include/IMP/score_states/NonbondedListScoreState.h (working copy) @@ -9,44 +9,65 @@ #define __IMP_NONBONDED_LIST_SCORE_STATE_H #include "../ScoreState.h" -#include "../internal/ParticleGrid.h" #include "../internal/kernel_version_info.h" +#include "../decorators/XYZDecorator.h" #include "BondedListScoreState.h" +#include + #include #include namespace IMP { +namespace internal +{ +class NBLAddPairIfNonbonded; +class NBLAddIfNonbonded; +} + typedef std::vector BondedListScoreStates; -//! A base class for classes that maintain a list of non-bonded pairs. -/** +//! An abstract base class for classes that maintain a list of non-bonded pairs. +/** \note If no value for the radius key is specified, all radii are + considered to be zero. */ class IMPDLLEXPORT NonbondedListScoreState: public ScoreState { private: FloatKey rk_; - -protected: - // made protected for debugging code, do not use otherwise + // the distance cutoff to count a pair as adjacent + float cutoff_; + /* How much to add to the size of particles to allow particles to move + without rebuilding the list*/ + Float slack_; + /* An estimate of what the slack should be next time the list is recomputed*/ + Float next_slack_; + int num_steps_; + bool nbl_is_valid_; + int number_of_rebuilds_; + int number_of_updates_; typedef std::vector > NBL; NBL nbl_; - float last_cutoff_; +protected: + + Float get_slack() const {return slack_;} + Float get_cutoff() const {return cutoff_;} + unsigned int size_nbl() const {return nbl_.size();} //! rebuild the nonbonded list - /** \internal + /** This is used by classes which inherit from this to rebuild the NBL + + \internal */ - virtual void rebuild_nbl(float cut)=0; + virtual void rebuild_nbl()=0; - void clear_nbl() { - last_cutoff_=-1; - nbl_.clear(); - } + + // Return true if two particles are bonded bool are_bonded(Particle *a, Particle *b) const { IMP_assert(a->get_is_active() && b->get_is_active(), "Inactive particles should have been removed" @@ -59,28 +80,13 @@ } return false; } + friend struct internal::NBLAddPairIfNonbonded; + friend struct internal::NBLAddIfNonbonded; - struct AddToNBL; - friend struct AddToNBL; - // these should be two separate classes, but they are not - struct AddToNBL{ - NonbondedListScoreState *state_; - Particle *p_; - AddToNBL(NonbondedListScoreState *s, Particle *p): state_(s), - p_(p){} - void operator()(Particle *p) { - operator()(p_, p); - } - void operator()(Particle *a, Particle *b) { - state_->add_to_nbl(a,b); - } - }; - //! tell the bonded lists what particles to pay attention to - void propagate_particles(const Particles &ps); - - void add_to_nbl(Particle *a, Particle *b) { + // Add the pair to the nbl if the particles are not bonded + void add_if_nonbonded(Particle *a, Particle *b) { IMP_assert(a->get_is_active() && b->get_is_active(), "Inactive particles should have been stripped"); @@ -95,71 +101,162 @@ } } - Float get_radius(Particle *a) const { - if (rk_ != FloatKey() && a->has_attribute(rk_)) { - return a->get_value(rk_); - } else { - return 0; + struct GetRadius{ + FloatKey rk_; + GetRadius(FloatKey rk): rk_(rk){ } + GetRadius(){} + Float operator() (Particle *a) const { + if (rk_ != FloatKey() && a->has_attribute(rk_)) { + return a->get_value(rk_); + } else { + return 0; + } + } + }; + + // Check if the bounding boxes of the spheres overlap + void add_if_box_overlap(Particle *a, Particle *b) { + BoxesOverlap bo= boxes_overlap_object(slack_+ cutoff_); + if (!bo(a, b)) { + IMP_LOG(VERBOSE, "Pair " << a->get_index() + << " and " << b->get_index() << " rejected on coordinate " + << std::endl); + } + IMP_LOG(VERBOSE, "Adding pair " << a->get_index() + << " and " << b->get_index() << std::endl); + add_if_nonbonded(a, b); } + GetRadius get_radius_object() const { + return GetRadius(rk_); + } + + // return true if the nbl was invalidated by a move of mc + // rebuild_cost is an estimate of the rebuilding cost used + // to try to get the slack right. The units are that of + // evaluating a distance + bool update(Float mc, Float rebuild_cost); + + // Return a list of all the particles in in with a radius field + Particles particles_with_radius(const Particles &in) const; + + + // Return true if the nbl is valid + bool get_nbl_is_valid() const {return nbl_is_valid_;} + + void set_nbl_is_valid(bool tf); + + struct BoxesOverlap { + GetRadius gr_; + Float slack_; + BoxesOverlap(GetRadius rk, Float sl): gr_(rk), slack_(sl){} + BoxesOverlap(){} + bool operator()(Particle *a, Particle *b) const { + XYZDecorator da(a); + XYZDecorator db(b); + float ra= gr_(a); + float rb= gr_(b); + for (unsigned int i=0; i< 3; ++i) { + float delta=std::abs(da.get_coordinate(i) - db.get_coordinate(i)); + if (delta -ra -rb > slack_) { + return false; + } + } + return true; + } + + bool operator()(std::pair p) const { + return operator()(p.first, p.second); + } + }; + + BoxesOverlap boxes_overlap_object(float cut) const { + return BoxesOverlap(get_radius_object(), cut); + } public: /** */ - NonbondedListScoreState(FloatKey rk); + NonbondedListScoreState(Float cutoff, FloatKey rk); + ~NonbondedListScoreState(); FloatKey get_radius_key() const {return rk_;} void set_radius_key(FloatKey rk) {rk_=rk;} IMP_CONTAINER(BondedListScoreState, bonded_list, BondedListIndex); + // 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) + virtual void show(std::ostream &out=std::cout) const; + virtual IMP::VersionInfo get_version_info() const { + return internal::kernel_version_info; + } //! An iterator through nonbonded particles - /** The value type is an std::pair + /** The value type is an ParticlePair. */ - typedef NBL::const_iterator NonbondedIterator; + typedef boost::filter_iterator NonbondedIterator; //! This iterates through the pairs of non-bonded particles - /** \param[in] cutoff The state may ignore pairs which are futher - apart than the cutoff. - \note that this is highly unsafe and iteration can only be - done once at a time. I will fix that eventually. + NonbondedIterator nonbonded_begin() const { + IMP_check(get_nbl_is_valid(), "Must call update first", + ValueException("Must call update first")); + return NonbondedIterator(boxes_overlap_object(cutoff_), + nbl_.begin(), nbl_.end()); + } + NonbondedIterator nonbonded_end() const { + IMP_check(get_nbl_is_valid(), "Must call update first", + ValueException("Must call update first")); + return NonbondedIterator(boxes_overlap_object(cutoff_), + nbl_.end(), nbl_.end()); + } - \note The distance cutoff is the l2 norm between the 3D coordinates - of the Particles. It ignore any size that may be associated with - the particles. - */ - NonbondedIterator nonbonded_begin(Float cutoff - =std::numeric_limits::max()) { - IMP_assert(last_cutoff_== cutoff || last_cutoff_==-1, - "Bad things are happening with the iterators in " - << "NonbondedListScoreState"); - if (last_cutoff_ < cutoff) { - IMP_LOG(VERBOSE, "Rebuilding NBL cutoff " << cutoff << std::endl); - clear_nbl(); - rebuild_nbl(cutoff); - last_cutoff_=cutoff; - } - return nbl_.begin(); + + unsigned int number_of_nonbonded() const { + IMP_check(get_nbl_is_valid(), "Must call update first", + ValueException("Must call update first")); + return nbl_.size(); } - NonbondedIterator nonbonded_end(Float cutoff - =std::numeric_limits::max()) { - if (last_cutoff_ < cutoff) { - IMP_LOG(VERBOSE, "Rebuilding NBL cutoff " << cutoff << std::endl); - clear_nbl(); - rebuild_nbl(cutoff); - last_cutoff_=cutoff; - } - return nbl_.end(); + + + unsigned int get_memory_usage() const { + return ScoreState::get_memory_usage() + + sizeof(FloatKey) + 3*sizeof(Float) + 3*sizeof(int) + + nbl_.capacity() * sizeof(std::pair) + + get_bonded_lists_memory_usage(); } + + void show_statistics(std::ostream &out=std::cout) const; + }; + +namespace internal +{ + struct NBLAddPairIfNonbonded{ + NonbondedListScoreState *state_; + NBLAddPairIfNonbonded(NonbondedListScoreState *s): state_(s){} + void operator()(Particle *a, Particle *b) { + state_->add_if_nonbonded(a,b); + } + }; + + struct NBLAddIfNonbonded{ + NonbondedListScoreState *state_; + Particle *p_; + NBLAddIfNonbonded(NonbondedListScoreState *s, Particle *p): state_(s), + p_(p){} + void operator()(Particle *p) { + state_->add_if_nonbonded(p_, p); + } + }; +} + } // namespace IMP #endif /* __IMP_NONBONDED_LIST_SCORE_STATE_H */ Index: kernel/include/IMP/score_states/SConscript =================================================================== --- kernel/include/IMP/score_states/SConscript (revision 541) +++ kernel/include/IMP/score_states/SConscript (working copy) @@ -1,15 +1,15 @@ +Import('env') import os.path -Import('env') - -files = ['BondedListScoreState.h', 'NonbondedListScoreState.h', - 'BipartiteNonbondedListScoreState.h', 'MaxChangeScoreState.h', - 'BondDecoratorListScoreState.h', 'AllNonbondedListScoreState.h', - 'QuadraticBipartiteNonbondedListScoreState.h', - 'QuadraticAllNonbondedListScoreState.h', - 'QuadraticNonbondedListScoreState.h', 'GravityCenterScoreState.h', - 'ChainBondedListScoreState.h'] - -# Install the include files: -includedir = os.path.join(env['includedir'], 'IMP', 'score_states') +files=[ + 'AllNonbondedListScoreState.h', + 'BipartiteNonbondedListScoreState.h', + 'BondDecoratorListScoreState.h', + 'BondedListScoreState.h', + 'CoverBondsScoreState.h', + 'GravityCenterScoreState.h', + 'MaxChangeScoreState.h', + 'NonbondedListScoreState.h', + ] +includedir = os.path.join(env['includedir'], 'IMP', 'score_states' ) inst = env.Install(includedir, files) env.Alias('install', inst) Index: kernel/include/IMP/score_states/QuadraticBipartiteNonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/QuadraticBipartiteNonbondedListScoreState.h (revision 541) +++ kernel/include/IMP/score_states/QuadraticBipartiteNonbondedListScoreState.h (working copy) @@ -1,49 +0,0 @@ -/** - * \file QuadraticBipartiteNonbondedListScoreState.h - * \brief Allow iteration through pairs of a set of atoms. - * - * Copyright 2007-8 Sali Lab. All rights reserved. - */ - -#ifndef __IMP_QUADRATIC_BIPARTITE_NONBONDED_LIST_SCORE_STATE_H -#define __IMP_QUADRATIC_BIPARTITE_NONBONDED_LIST_SCORE_STATE_H - -#include "QuadraticNonbondedListScoreState.h" -#include "../internal/ParticleGrid.h" -#include "../internal/kernel_version_info.h" - -#include -#include - -namespace IMP -{ - -//! This class maintains a list of non-bonded pairs between two sets. -/** The class works roughly like the NonbondedListScoreState except - only pairs where one particle is taken from each set are returned. - - \note QuadraticBipartiteNonbondedListScoreState is basically an - implementation detail for performance analysis and should not - be used by end users. - */ -class IMPDLLEXPORT QuadraticBipartiteNonbondedListScoreState: - public QuadraticNonbondedListScoreState -{ - typedef QuadraticNonbondedListScoreState P; - unsigned int na_; - - virtual void rebuild_nbl(float cut); -public: - QuadraticBipartiteNonbondedListScoreState(FloatKey rk, - const Particles &ps0, - const Particles &ps1); - QuadraticBipartiteNonbondedListScoreState(FloatKey rk); - - IMP_SCORE_STATE(internal::kernel_version_info) - - void set_particles(const Particles &ps0, const Particles &ps1); -}; - -} // namespace IMP - -#endif /* __IMP_QUADRATIC_BIPARTITE_NONBONDED_LIST_SCORE_STATE_H */ Index: kernel/include/IMP/score_states/CoverBondsScoreState.h =================================================================== --- kernel/include/IMP/score_states/CoverBondsScoreState.h (revision 0) +++ kernel/include/IMP/score_states/CoverBondsScoreState.h (revision 0) @@ -0,0 +1,55 @@ +/** + * \file CoverBondsScoreState.h + * \brief Covers a set of bonds with spheres. + * + * Copyright 2007-8 Sali Lab. All rights reserved. + */ + +#ifndef __IMP_COVER_BONDS_SCORE_STATE_H +#define __IMP_COVER_BONDS_SCORE_STATE_H + +#include "../ScoreState.h" +#include "../internal/kernel_version_info.h" + +namespace IMP +{ + class BondDecoratorListScoreState; + +//! This class sets the position and radius of each bond to cover the endpoints. +/** For each bond in the list of particles set the x,y,z to be the + average of that of the bonds endpoints and the radius to be the + sqrt of half the distance. + + Note that the ends of the bond are treated as points so any radius + they might have is ignored. This would be easy to change if + desired. + + This is designed to be used in conjunction with a bonded list. Make sure + it is updated first. + */ +class IMPDLLEXPORT CoverBondsScoreState: public ScoreState +{ + internal::ObjectPointer bl_; + FloatKey rk_; +public: + /** Get the list of bonds from the BondDecoratorListScoreState. This list is + not owned and update is not called on this list automatically. + */ + CoverBondsScoreState(BondDecoratorListScoreState *bl, + FloatKey rk=FloatKey("radius")); + ~CoverBondsScoreState(); + + IMP_SCORE_STATE(internal::kernel_version_info); + + //! apply the derivatives to the endpoints + void after_evaluate(DerivativeAccumulator *accpt); + + unsigned int get_memory_usage() const { + return ScoreState::get_memory_usage() + + sizeof(rk_) + sizeof(bl_); + } +}; + +} // namespace IMP + +#endif /* __IMP_NONBONDED_LIST_SCORE_STATE_H */ Index: kernel/include/IMP/score_states/BipartiteNonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/BipartiteNonbondedListScoreState.h (revision 541) +++ kernel/include/IMP/score_states/BipartiteNonbondedListScoreState.h (working copy) @@ -1,32 +1,91 @@ /** * \file BipartiteNonbondedListScoreState.h - * \brief Allow iteration through pairs of a set of atoms. + * \brief Bipartiteow iteration through pairs of a set of s. * - * Copyright 2007-8 Sali Lab. All rights reserved. + * Copyright 2007-8 Sali Lab. Bipartite rights reserved. */ #ifndef __IMP_BIPARTITE_NONBONDED_LIST_SCORE_STATE_H #define __IMP_BIPARTITE_NONBONDED_LIST_SCORE_STATE_H -#include "QuadraticBipartiteNonbondedListScoreState.h" +#include "NonbondedListScoreState.h" +#include "../internal/kernel_version_info.h" +#include "MaxChangeScoreState.h" + namespace IMP { -//! Maintain a nonbonded list between two disjoint sets. -/** - \note If no value for the radius key is specified, all radii are - considered to be zero. - \ingroup restraint +//! This class maintains a list of non-bonded pairs of spheres between two sets +/** To iterate through the list of pairs use the NonbondedListScoreState::begin, + NonbondedListScoreState::end functions. + + The radius key can be the default key. + + Changes in coordinates and radii are handled properly. + + \ingroup restraint */ -class BipartiteNonbondedListScoreState: - public QuadraticBipartiteNonbondedListScoreState { - typedef QuadraticBipartiteNonbondedListScoreState P; +class IMPDLLEXPORT BipartiteNonbondedListScoreState: + public NonbondedListScoreState +{ + typedef NonbondedListScoreState P; public: - BipartiteNonbondedListScoreState(FloatKey rk, + //! What algorithm to use to perform the computations + enum Algorithm { + //! Check all pairs of particles to see if they are close enough + QUADRATIC, + //! Sweep space looking for intersecting bounding boxes. + BBOX, + //! Choose whatever is best that is available + DEFAULT + }; +protected: + Algorithm a_; + internal::ObjectPointer mc0_, mc1_, mcr_; + + void process_sets(const Particles &p0, + const Particles &p1); + + //! \internal + void rebuild_nbl(); + void check_nbl() const; +public: + + /**\param[in] cutoff The distance cutoff to use. + \param[in] radius The key to use to get the radius + \param[in] a Which algorithm to use. The default is the best. + */ + BipartiteNonbondedListScoreState(Float cutoff, FloatKey radius, + Algorithm a=DEFAULT); + /**\param[in] cutoff The distance cutoff to use. + \param[in] radius The key to use to get the radius + \param[in] ps0 The first set. + \param[in] ps1 The second set. + \param[in] a Which algorithm to use. The default is the best. + */ + BipartiteNonbondedListScoreState(Float cutoff, + FloatKey radius, const Particles &ps0, - const Particles &ps1): P(rk, ps0, ps1){} - BipartiteNonbondedListScoreState(FloatKey rk): P(rk){} + const Particles &ps1, + Algorithm a=DEFAULT); + + ~BipartiteNonbondedListScoreState(); + IMP_SCORE_STATE(internal::kernel_version_info); + + //! Add the particles to the first set + void add_particles_0(const Particles &ps); + //! Add the particles to the second set + void add_particles_1(const Particles &ps); + //! Remove all particles + void clear_particles(); + //! Replace the set of particles + void set_particles(const Particles &ps0, const Particles &ps1); + + /** If there is CGAL support, a more efficient algorithm BBOX can be used*/ + void set_algorithm(Algorithm a); + + unsigned int get_memory_usage() const; }; } // namespace IMP Index: kernel/include/IMP/score_states/GravityCenterScoreState.h =================================================================== --- kernel/include/IMP/score_states/GravityCenterScoreState.h (revision 541) +++ kernel/include/IMP/score_states/GravityCenterScoreState.h (working copy) @@ -36,22 +36,21 @@ */ GravityCenterScoreState(Particle *center, FloatKey weightkey=FloatKey(), const Particles &ps= Particles()); - virtual ~GravityCenterScoreState() {} + virtual ~GravityCenterScoreState(); //! Set the position of the center particle from the original points. - void set_position(); + void update_position(); - void update() { - set_position(); - } + IMP_SCORE_STATE(internal::kernel_version_info); - void after_evaluate(DerivativeAccumulator *accpt) { +protected: + + void do_after_evaluate(DerivativeAccumulator *accpt) { if (accpt) { transform_derivatives(accpt); } } -protected: //! Back-transform any forces on the center particle to the original points. void transform_derivatives(DerivativeAccumulator *accpt); Index: kernel/include/IMP/score_states/BondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/BondedListScoreState.h (revision 541) +++ kernel/include/IMP/score_states/BondedListScoreState.h (working copy) @@ -25,9 +25,6 @@ BondedListScoreState(){} virtual ~BondedListScoreState(){} - //! Set the set of particles used - virtual void set_particles(const Particles &ps)=0; - //! Return true if the two particles are bonded virtual bool are_bonded(Particle *a, Particle *b) const =0; }; Index: kernel/include/IMP/score_states/QuadraticAllNonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/QuadraticAllNonbondedListScoreState.h (revision 541) +++ kernel/include/IMP/score_states/QuadraticAllNonbondedListScoreState.h (working copy) @@ -1,49 +0,0 @@ -/** - * \file QuadraticAllNonbondedListScoreState.h - * \brief Allow iteration through pairs of a set of spheres. - * - * Copyright 2007-8 Sali Lab. All rights reserved. - */ - -#ifndef __IMP_QUADRATIC_ALL_NONBONDED_LIST_SCORE_STATE_H -#define __IMP_QUADRATIC_ALL_NONBONDED_LIST_SCORE_STATE_H - -#include "NonbondedListScoreState.h" -#include "../internal/kernel_version_info.h" -#include "QuadraticNonbondedListScoreState.h" - -#include -#include - -namespace IMP -{ - -//! This class maintains a list of non-bonded pairs of spheres -/** \note QuadraticBipartiteNonbondedListScoreState is basically an - implementation detail for performance analysis and should not - be used by end users. - */ -class IMPDLLEXPORT QuadraticAllNonbondedListScoreState: - public QuadraticNonbondedListScoreState -{ - typedef QuadraticNonbondedListScoreState P; - - //! \internal - void rebuild_nbl(Float cut); - -public: - /** - \param[in] ps A list of particles to use. - \param[in] radius The key to use to get the radius - */ - QuadraticAllNonbondedListScoreState(FloatKey radius, - const Particles &ps= Particles()); - ~QuadraticAllNonbondedListScoreState(); - IMP_SCORE_STATE(internal::kernel_version_info) - - void set_particles(const Particles &ps); -}; - -} // namespace IMP - -#endif /* __IMP_QUADRATIC_ALL_NONBONDED_LIST_SCORE_STATE_H */ Index: kernel/include/IMP/score_states/ChainBondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/ChainBondedListScoreState.h (revision 541) +++ kernel/include/IMP/score_states/ChainBondedListScoreState.h (working copy) @@ -1,44 +0,0 @@ -/** - * \file ChainBondedListScoreState.h - * \brief Allow iteration through pairs of a set of atoms. - * - * Copyright 2007-8 Sali Lab. All rights reserved. - */ - -#ifndef __IMP_CHAIN_BONDED_LIST_SCORE_STATE_H -#define __IMP_CHAIN_BONDED_LIST_SCORE_STATE_H - -#include "BondedListScoreState.h" -#include "../decorators/bond_decorators.h" -#include "../internal/Vector.h" - -namespace IMP -{ - -//! Exclude consecutive particles in a chain. -/** \ingroup bond - */ -class IMPDLLEXPORT ChainBondedListScoreState: public BondedListScoreState -{ - std::vector > chains_; - IntKey cik_; - unsigned int next_index_; -public: - //! Find bonds amongst the following points. - /** \param [in] ps The set of particles to use. - */ - ChainBondedListScoreState(); - virtual ~ChainBondedListScoreState() {} - - void add_chain(const Particles &ps); - - void clear_chains(); - - virtual bool are_bonded(Particle *a, Particle *b) const; - - virtual void update(); -}; - -} // namespace IMP - -#endif /* __IMP_CHAIN_BONDED_LIST_SCORE_STATE_H */ Index: kernel/include/IMP/score_states/MaxChangeScoreState.h =================================================================== --- kernel/include/IMP/score_states/MaxChangeScoreState.h (revision 541) +++ kernel/include/IMP/score_states/MaxChangeScoreState.h (working copy) @@ -27,7 +27,6 @@ { FloatKeys keys_; FloatKeys origkeys_; - Particles ps_; float max_change_; public: //! Track the changes with the specified keys. @@ -36,7 +35,7 @@ virtual ~MaxChangeScoreState(){} - virtual void update(); + IMP_SCORE_STATE(internal::kernel_version_info); //! Measure differences from the current value. void reset(); @@ -46,6 +45,13 @@ return max_change_; } + unsigned int get_memory_usage() const { + return keys_.capacity()*sizeof(FloatKey) + + origkeys_.capacity()*sizeof(FloatKey) + + get_particles_memory_usage() + sizeof(float) + + ScoreState::get_memory_usage(); + } + IMP_LIST(public, Particle, particle, Particle*); }; Index: kernel/include/IMP/score_states/AllNonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/AllNonbondedListScoreState.h (revision 541) +++ kernel/include/IMP/score_states/AllNonbondedListScoreState.h (working copy) @@ -9,76 +9,98 @@ #define __IMP_ALL_NONBONDED_LIST_SCORE_STATE_H #include "NonbondedListScoreState.h" -#include "../internal/ParticleGrid.h" #include "../internal/kernel_version_info.h" -#include "MaxChangeScoreState.h" +#include "../internal/Vector.h" -#include -#include namespace IMP { -//! This class maintains a list of non-bonded pairs of particles -/** \note If no value for the radius key is specified, all radii are - considered to be zero. +class MaxChangeScoreState; +namespace internal +{ +class ParticleGrid; +} - \note The radius is currently assumed not to change. This could - be fixed later. +//! This class maintains a list of non-bonded pairs of spheres +/** To iterate through the list of pairs use the NonbondedListScoreState::begin, + NonbondedListScoreState::end functions. - \todo 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). + The radius key can be the default key. + Changes in coordinates and radii are handled properly unless grid is used, + then changes in radii may not be handled properly. + \ingroup restraint */ class IMPDLLEXPORT AllNonbondedListScoreState: public NonbondedListScoreState { typedef NonbondedListScoreState P; + public: + //! What algorithm to use to perform the computations + enum Algorithm { + //! Check all pairs of particles to see if they are close enough + QUADRATIC, + GRID, + //! Sweep space looking for intersecting bounding boxes. + BBOX, + //! Choose the best algorithm available + DEFAULT + }; 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_; + Algorithm a_; + internal::ObjectPointer mc_, mcr_; + void check_nbl() const; + //! \internal - void rebuild_nbl(Float cut); + void rebuild_nbl(); - void repartition_points(const Particles &ps, std::vector &out); - float side_from_r(float r) const; + // methods for grid + void grid_rebuild_nbl(); - void generate_nbl(const Bin &particle_bin, const Bin &grid_bin, float cut); + void grid_partition_points(internal::Vector &bins); - void cleanup(std::vector &bins); + void grid_generate_nbl(const internal::ParticleGrid *particle_bin, + const internal::ParticleGrid *grid_bin); + float grid_side_from_r(float r) const; public: - /** - \param[in] ps A list of particles to use. + + + /**\param[in] cutoff The distance cutoff to use. \param[in] radius The key to use to get the radius + \param[in] a Which algorithm to use. The default is the best. */ - AllNonbondedListScoreState(FloatKey radius, - const Particles &ps=Particles()); + AllNonbondedListScoreState(Float cutoff, + FloatKey radius, + Algorithm a= DEFAULT); + /**\param[in] cutoff The distance cutoff to use. + \param[in] radius The key to use to get the radius + \param[in] ps A list of particles to use. + \param[in] a Which algorithm to use. The default is the best. + */ + AllNonbondedListScoreState(Float cutoff, + FloatKey radius, + const Particles &ps, + Algorithm a= DEFAULT); ~AllNonbondedListScoreState(); - IMP_SCORE_STATE(internal::kernel_version_info) - void set_particles(const Particles &ps); + IMP_SCORE_STATE(internal::kernel_version_info); - //! Add a few particles to the nonbonded list - /** Note that this invalidates the nonbonded list. - \todo We could just add newly created pairs to the nonbonded list. - */ + //! Add the particles and add them to the NBL void add_particles(const Particles &ps); + //! Remove all particles + void clear_particles(); + //! Replace the set of particles + void set_particles(const Particles &ps); - //! Return a list of all the particles used - Particles get_particles() const; + /** If there is CGAL support, a more efficient algorithm BBOX can be used*/ + void set_algorithm(Algorithm a); + + unsigned int get_memory_usage() const; }; } // namespace IMP Index: kernel/include/IMP/restraints/SingletonListRestraint.h =================================================================== --- kernel/include/IMP/restraints/SingletonListRestraint.h (revision 541) +++ kernel/include/IMP/restraints/SingletonListRestraint.h (working copy) @@ -39,6 +39,10 @@ using Restraint::clear_particles; using Restraint::set_particles; + unsigned int get_memory_usage() const { + return sizeof(ss_) + Restraint::get_memory_usage(); + } + protected: internal::ObjectPointer ss_; }; Index: kernel/include/IMP/restraints/ConnectivityRestraint.h =================================================================== --- kernel/include/IMP/restraints/ConnectivityRestraint.h (revision 541) +++ kernel/include/IMP/restraints/ConnectivityRestraint.h (working copy) @@ -50,6 +50,10 @@ IMP_RESTRAINT(internal::kernel_version_info) + unsigned int get_memory_usage() const { + return Restraint::get_memory_usage() + sizeof(ps_) + + set_offsets_.capacity()*sizeof(unsigned int); + } protected: internal::ObjectPointer ps_; Index: kernel/include/IMP/restraints/PairChainRestraint.h =================================================================== --- kernel/include/IMP/restraints/PairChainRestraint.h (revision 541) +++ kernel/include/IMP/restraints/PairChainRestraint.h (working copy) @@ -42,6 +42,11 @@ //! Clear all the stored chains void clear_chains(); + unsigned int get_memory_usage() const { + return sizeof(ts_) + sizeof(unsigned int)*chain_splits_.capacity() + + Restraint::get_memory_usage(); + } + protected: internal::ObjectPointer ts_; std::vector chain_splits_; Index: kernel/include/IMP/restraints/SConscript =================================================================== --- kernel/include/IMP/restraints/SConscript (revision 541) +++ kernel/include/IMP/restraints/SConscript (working copy) @@ -1,14 +1,18 @@ +Import('env') import os.path -Import('env') - -files = ['ConnectivityRestraint.h', - 'DistanceRestraint.h', 'AngleRestraint.h', - 'DihedralRestraint.h', 'RestraintSet.h', - 'NonbondedRestraint.h', 'BondDecoratorRestraint.h', - 'SingletonListRestraint.h', 'PairListRestraint.h', - 'TripletChainRestraint.h', 'PairChainRestraint.h'] - -# Install the include files: -includedir = os.path.join(env['includedir'], 'IMP', 'restraints') +files=[ + 'AngleRestraint.h', + 'BondDecoratorRestraint.h', + 'ConnectivityRestraint.h', + 'DihedralRestraint.h', + 'DistanceRestraint.h', + 'NonbondedRestraint.h', + 'PairChainRestraint.h', + 'PairListRestraint.h', + 'RestraintSet.h', + 'SingletonListRestraint.h', + 'TripletChainRestraint.h', + ] +includedir = os.path.join(env['includedir'], 'IMP', 'restraints' ) inst = env.Install(includedir, files) env.Alias('install', inst) Index: kernel/include/IMP/restraints/BondDecoratorRestraint.h =================================================================== --- kernel/include/IMP/restraints/BondDecoratorRestraint.h (revision 541) +++ kernel/include/IMP/restraints/BondDecoratorRestraint.h (working copy) @@ -45,6 +45,13 @@ IMP_RESTRAINT(internal::kernel_version_info) + void set_function(UnaryFunction *f) {f_=f;} + + unsigned int get_memory_usage() const { + return Restraint::get_memory_usage() + sizeof(bl_) + + sizeof(f_); + } + protected: BondDecoratorListScoreState *bl_; internal::ObjectPointer f_; Index: kernel/include/IMP/restraints/NonbondedRestraint.h =================================================================== --- kernel/include/IMP/restraints/NonbondedRestraint.h (revision 541) +++ kernel/include/IMP/restraints/NonbondedRestraint.h (working copy) @@ -34,18 +34,17 @@ object is deleted upon destruction. \param[in] nbl The non-bonded list to use to get the list of particles. - \param[in] max_dist Pairs beyond this distance may be dropped. */ - NonbondedRestraint(PairScore *ps, NonbondedListScoreState *nbl, - Float max_dist= std::numeric_limits::max()); + NonbondedRestraint(PairScore *ps, NonbondedListScoreState *nbl); virtual ~NonbondedRestraint(){} IMP_RESTRAINT(internal::kernel_version_info) + unsigned int get_memory_usage() const; + protected: NonbondedListScoreState *nbl_; internal::ObjectPointer sf_; - Float max_dist_; }; } // namespace IMP Index: kernel/include/IMP/restraints/AngleRestraint.h =================================================================== --- kernel/include/IMP/restraints/AngleRestraint.h (revision 541) +++ kernel/include/IMP/restraints/AngleRestraint.h (working copy) @@ -16,6 +16,7 @@ namespace IMP { + class UnaryFunction; //! Angle restraint between three particles @@ -34,6 +35,8 @@ IMP_RESTRAINT(internal::kernel_version_info) + unsigned int get_memory_usage() const ; + protected: internal::ObjectPointer sf_; }; Index: kernel/include/IMP/restraints/TripletChainRestraint.h =================================================================== --- kernel/include/IMP/restraints/TripletChainRestraint.h (revision 541) +++ kernel/include/IMP/restraints/TripletChainRestraint.h (working copy) @@ -42,6 +42,10 @@ //! Clear all the stored chains void clear_chains(); + unsigned int get_memory_usage() const { + return sizeof(ts_) + sizeof(unsigned int)*chain_splits_.capacity() + + Restraint::get_memory_usage(); + } protected: internal::ObjectPointer ts_; std::vector chain_splits_; Index: kernel/include/IMP/restraints/PairListRestraint.h =================================================================== --- kernel/include/IMP/restraints/PairListRestraint.h (revision 541) +++ kernel/include/IMP/restraints/PairListRestraint.h (working copy) @@ -39,7 +39,15 @@ void add_particle_pair(ParticlePair p); void clear_particle_pairs(); void add_particle_pairs(const ParticlePairs &ps); + void set_particle_pairs(const ParticlePairs &ps) { + clear_particle_pairs(); + add_particle_pairs(ps); + } + unsigned int get_memory_usage() const { + return sizeof(ss_) + Restraint::get_memory_usage(); + } + protected: internal::ObjectPointer ss_; }; Index: kernel/include/IMP/restraints/DistanceRestraint.h =================================================================== --- kernel/include/IMP/restraints/DistanceRestraint.h (revision 541) +++ kernel/include/IMP/restraints/DistanceRestraint.h (working copy) @@ -39,6 +39,10 @@ IMP_RESTRAINT(internal::kernel_version_info) + unsigned int get_memory_usage() const { + return sizeof(dp_) + Restraint::get_memory_usage(); + } + protected: //! scoring function for this restraint DistancePairScore dp_; Index: kernel/include/IMP/restraints/RestraintSet.h =================================================================== --- kernel/include/IMP/restraints/RestraintSet.h (revision 541) +++ kernel/include/IMP/restraints/RestraintSet.h (working copy) @@ -42,6 +42,11 @@ //! Get weight for all restraints contained by this set. Float get_weight() const { return weight_; } + unsigned int get_memory_usage() const { + return sizeof(weight_) + Restraint::get_memory_usage() + + sizeof(char)* name_.size() + get_restraints_memory_usage(); + } + protected: //! Weight for all restraints. Index: kernel/include/IMP/restraints/DihedralRestraint.h =================================================================== --- kernel/include/IMP/restraints/DihedralRestraint.h (revision 541) +++ kernel/include/IMP/restraints/DihedralRestraint.h (working copy) @@ -34,6 +34,10 @@ IMP_RESTRAINT(internal::kernel_version_info) + unsigned int get_memory_usage() const { + return sizeof(score_func_) + Restraint::get_memory_usage(); + } + protected: //! scoring function for this restraint UnaryFunction* score_func_; Index: kernel/include/IMP/Restraint.h =================================================================== --- kernel/include/IMP/Restraint.h (revision 541) +++ kernel/include/IMP/Restraint.h (working copy) @@ -94,6 +94,10 @@ return model_.get(); } + unsigned int get_memory_usage() const { + return sizeof(model_)+ sizeof(bool) + get_particles_memory_usage(); + } + IMP_LIST(protected, Particle, particle, Particle*) private: Index: kernel/include/IMP/decorators/XYZDecorator.h =================================================================== --- kernel/include/IMP/decorators/XYZDecorator.h (revision 541) +++ kernel/include/IMP/decorators/XYZDecorator.h (working copy) @@ -51,6 +51,12 @@ void set_coordinate(unsigned int i, Float v) { get_particle()->set_value(get_coordinate_key(i), v); } + //! set the ith coordinate + void set_coordinates(const Vector3D &v) { + set_x(v[0]); + set_y(v[1]); + set_z(v[2]); + } //! Get the ith coordinate Float get_coordinate(int i) const { return get_particle()->get_value(get_coordinate_key(i)); @@ -117,11 +123,18 @@ IMP_OUTPUT_OPERATOR(XYZDecorator); -//! Compute the distance between a pair of particles +//! Compute the distance between a pair of points /** \ingroup helper */ IMPDLLEXPORT Float distance(XYZDecorator a, XYZDecorator b); + +//! Compute the distance between a pair of spheres +/** \ingroup helper + */ +IMPDLLEXPORT Float sphere_distance(XYZDecorator a, XYZDecorator b, FloatKey rk); + + } // namespace IMP #endif /* __IMP_XYZ_DECORATOR_H */ Index: kernel/include/IMP/decorators/SConscript =================================================================== --- kernel/include/IMP/decorators/SConscript (revision 541) +++ kernel/include/IMP/decorators/SConscript (working copy) @@ -1,11 +1,16 @@ +Import('env') import os.path -Import('env') - -files = ['HierarchyDecorator.h', 'NameDecorator.h', 'utility.h', - 'XYZDecorator.h', 'ResidueDecorator.h', 'AtomDecorator.h', - 'MolecularHierarchyDecorator.h', 'bond_decorators.h', 'macros.h'] - -# Install the include files: -includedir = os.path.join(env['includedir'], 'IMP', 'decorators') +files=[ + 'AtomDecorator.h', + 'bond_decorators.h', + 'HierarchyDecorator.h', + 'macros.h', + 'MolecularHierarchyDecorator.h', + 'NameDecorator.h', + 'ResidueDecorator.h', + 'utility.h', + 'XYZDecorator.h', + ] +includedir = os.path.join(env['includedir'], 'IMP', 'decorators' ) inst = env.Install(includedir, files) env.Alias('install', inst) Index: kernel/include/IMP/optimizers/states/SConscript =================================================================== --- kernel/include/IMP/optimizers/states/SConscript (revision 541) +++ kernel/include/IMP/optimizers/states/SConscript (working copy) @@ -1,10 +1,10 @@ +Import('env') import os.path -Import('env') - -files = ['VRMLLogOptimizerState.h', 'CMMLogOptimizerState.h', - 'VelocityScalingOptimizerState.h'] - -# Install the include files: -includedir = os.path.join(env['includedir'], 'IMP', 'optimizers', 'states') +files=[ + 'CMMLogOptimizerState.h', + 'VelocityScalingOptimizerState.h', + 'VRMLLogOptimizerState.h', + ] +includedir = os.path.join(env['includedir'], 'IMP', 'optimizers', 'states' ) inst = env.Install(includedir, files) env.Alias('install', inst) Index: kernel/include/IMP/optimizers/states/VRMLLogOptimizerState.h =================================================================== --- kernel/include/IMP/optimizers/states/VRMLLogOptimizerState.h (revision 541) +++ kernel/include/IMP/optimizers/states/VRMLLogOptimizerState.h (working copy) @@ -16,6 +16,9 @@ #include "../../OptimizerState.h" #include "../../internal/kernel_version_info.h" +#include +#include + namespace IMP { @@ -50,39 +53,36 @@ //! The float key to use for the radius /** Particles without such an attribute are drawn as fixed sized markers. */ - void set_radius(FloatKey k) { + void set_radius_key(FloatKey k) { radius_=k; } //! The three color components /** Color values should be between 0 and 1. They will be snapped if needed. */ - void set_color(FloatKey r, FloatKey g, FloatKey b) { - r_=r; g_=g; b_=b; + void set_color_key(IntKey k) { + color_=k; } - //! Set the particles to use. - void set_particles(const Particles &pis) { - pis_=pis; - } + //! Add a value to the color map + void set_color(int c, Float r, Float g, Float b); + + IMP_LIST(public, Particle, particle, Particle*); + //! Force it to write the next file void write_next_file(); void write(std::string name) const; +protected: //! A helper function to just write a list of particles to a file - static void write(const Particles &pis, - FloatKey radius, FloatKey r, - FloatKey g, FloatKey b, - std::ostream &out); - -protected: - Particles pis_; + void write(std::ostream &out) const; std::string filename_; int file_number_; int call_number_; int skip_steps_; FloatKey radius_; - FloatKey r_, g_, b_; + IntKey color_; + std::map > colors_; }; } // namespace IMP Index: kernel/include/IMP/optimizers/SConscript =================================================================== --- kernel/include/IMP/optimizers/SConscript (revision 541) +++ kernel/include/IMP/optimizers/SConscript (working copy) @@ -1,12 +1,16 @@ +Import('env') import os.path -Import('env') - -files = ['ConjugateGradients.h', 'SteepestDescent.h', 'MolecularDynamics.h', - 'MoverBase.h', 'MonteCarlo.h', 'Mover.h'] - -# Install the include files: -includedir = os.path.join(env['includedir'], 'IMP', 'optimizers') +files=[ + 'BrownianDynamics.h', + 'ConjugateGradients.h', + 'MolecularDynamics.h', + 'MonteCarlo.h', + 'MoverBase.h', + 'Mover.h', + 'SteepestDescent.h', + ] +includedir = os.path.join(env['includedir'], 'IMP', 'optimizers' ) inst = env.Install(includedir, files) env.Alias('install', inst) -SConscript('states/SConscript') -SConscript('movers/SConscript') +SConscript( 'movers/SConscript' ) +SConscript( 'states/SConscript' ) Index: kernel/include/IMP/optimizers/movers/SConscript =================================================================== --- kernel/include/IMP/optimizers/movers/SConscript (revision 541) +++ kernel/include/IMP/optimizers/movers/SConscript (working copy) @@ -1,8 +1,9 @@ +Import('env') import os.path -Import('env') - -files = ['BallMover.h', 'NormalMover.h'] - +files=[ + 'BallMover.h', + 'NormalMover.h', + ] includedir = os.path.join(env['includedir'], 'IMP', 'optimizers', 'movers' ) inst = env.Install(includedir, files) env.Alias('install', inst) Index: kernel/include/IMP/optimizers/BrownianDynamics.h =================================================================== --- kernel/include/IMP/optimizers/BrownianDynamics.h (revision 0) +++ kernel/include/IMP/optimizers/BrownianDynamics.h (revision 0) @@ -0,0 +1,182 @@ +/** + * \file BrownianDynamics.h \brief Simple molecular dynamics optimizer. + * + * Copyright 2007-8 Sali Lab. All rights reserved. + * + */ + +#ifndef __IMP_BROWNIAN_DYNAMICS_H +#define __IMP_BROWNIAN_DYNAMICS_H + +#include "../IMP_config.h" +#include "../Particle.h" +#include "../Optimizer.h" +#include "../internal/units.h" +#include "../internal/kernel_version_info.h" +#include "../Vector3D.h" + +namespace IMP +{ + +//! Simple brownian dynamics optimizer. +/** The particles to be optimized must have optimizable x,y,z attributes + and a non-optimizable "stokes radius"; this optimizer assumes + the score to be energy in kcal/mol, the xyz coordinates to be in + angstroms and the diffusion coefficent be in cm^2/s + + Particles without optimized x,y,z and nonoptimized D are skipped. + + \todo Currently dynamic time steps are handled by scaling back dt + if forces are too high at the current position. It would be better + to scale back the time step if the forces will be too high at the + next position. + + \ingroup optimizer + */ +class IMPDLLEXPORT BrownianDynamics : public Optimizer +{ +public: + //! \param[in] radkey The key for the stokes radius + BrownianDynamics(FloatKey dkey= FloatKey("D")); + virtual ~BrownianDynamics(); + + IMP_OPTIMIZER(internal::kernel_version_info) + + //! Set the max step that should be allowed in angstroms + /** The timestep is shrunk if any particles are moved more than this */ + void set_max_step(Float a) { + set_max_step(internal::Angstrom(a)); + } + + //! Set time step in fs + void set_time_step(Float t) { + set_time_step(internal::FemtoSecond(t)); + } + + //! Time step in fs + Float get_time_step() const {return get_time_step_units().get_value();} + + //! The last time step used in fs + Float get_last_time_step() const { + return get_last_time_step_units().get_value(); + } + + + void set_temperature(Float t) { set_temperature(internal::Kelvin(t));} + Float get_temperature() const { return get_temperature_units().get_value();} + + FloatKey get_d_key() const { return dkey_;} + + //! Estimate the diffusion coeffient from the mass in kD + /** This method does not do anything fancy, but is OK for large + globular proteins. + + Note that it depends on temperature and so can't be static. + */ + Float estimate_radius_from_mass(Float mass_in_kd) const { + return estimate_radius_from_mass_units(mass_in_kd).get_value(); + } + + //! Return the expected distance moved for a particle with a given D + /** The units on D are cm^2/sec and the return has units of Angstroms. + */ + Float compute_sigma_from_D(Float D) const { + return compute_sigma_from_D(internal::Cm2PerSecond(D)).get_value(); + } + + //! Estimate the diffusion coefficient from the radius in Angstroms + /** This depends on the temperature. + */ + Float estimate_D_from_radius(Float radius_in_angstroms) const { + internal::Angstrom r(radius_in_angstroms); + return estimate_D_from_radius(r).get_value(); + } + + //! Returns a force value which would move the particle by sigma + /** This value is in KCal/A/mol + */ + Float get_force_scale_from_D(Float D) const { + return get_force_scale_from_D(internal::Cm2PerSecond(D)).get_value(); + } + + //! Get the current time in femptoseconds + Float get_current_time() const { + return get_current_time_units().get_value(); + } + + //! Set the current time in femptoseconds + void set_current_time(Float fs) { + set_current_time( internal::FemtoSecond(fs)); + } + + //! Return a histogram of timesteps + /** + The 0 value is the bin from get_time_step() to get_time_step()/2 + and they go as factors of two from there. + */ + const std::vector get_time_step_histogram() const { + return time_steps_; + } + + IMP_LIST(private, Particle, particle, Particle*); + +protected: + //! Perform a single dynamics step. + // return true if the initial step size was OK + bool step(); + + // Propose a single step, this will be accepted if all moves are + // small enough. Return true if it should be accepted. + bool propose_step(std::vector &proposal); + + internal::FemtoJoule kt() const; + + void setup_particles(); + + void set_max_step(internal::Angstrom a) { + max_change_= a; + } + + void set_time_step(internal::FemtoSecond t); + + internal::FemtoSecond get_time_step_units() const {return max_dt_;} + + internal::FemtoSecond get_last_time_step_units() const {return cur_dt_;} + + void set_temperature(internal::Kelvin t) { T_=t;} + + internal::Kelvin get_temperature_units() const { return T_;} + + internal::Angstrom estimate_radius_from_mass_units(Float mass_in_kd) const; + + internal::Angstrom compute_sigma_from_D(internal::Cm2PerSecond D) const; + + internal::Cm2PerSecond + estimate_D_from_radius(internal::Angstrom radius) const; + + internal::KCalPerAMol get_force_scale_from_D(internal::Cm2PerSecond D) const; + + internal::FemtoSecond get_current_time_units() const { + return cur_time_; + } + + //! Set the current time in femptoseconds + void set_current_time(internal::FemtoSecond fs) { + cur_time_= fs; + } + + + + internal::Angstrom max_change_; + internal::FemtoSecond max_dt_, cur_dt_, cur_time_; + internal::Kelvin T_; + + + FloatKey dkey_; + + std::vector time_steps_; +}; + +} // namespace IMP + +#endif /* __IMP_BROWNIAN_DYNAMICS_H */ Index: kernel/include/IMP/pair_scores/SConscript =================================================================== --- kernel/include/IMP/pair_scores/SConscript (revision 541) +++ kernel/include/IMP/pair_scores/SConscript (working copy) @@ -1,9 +1,9 @@ +Import('env') import os.path -Import('env') - -files = ['DistancePairScore.h', 'SphereDistancePairScore.h'] - -# Install the include files: -includedir = os.path.join(env['includedir'], 'IMP', 'pair_scores') +files=[ + 'DistancePairScore.h', + 'SphereDistancePairScore.h', + ] +includedir = os.path.join(env['includedir'], 'IMP', 'pair_scores' ) inst = env.Install(includedir, files) env.Alias('install', inst) Index: kernel/include/IMP/internal/SConscript =================================================================== --- kernel/include/IMP/internal/SConscript (revision 541) +++ kernel/include/IMP/internal/SConscript (working copy) @@ -1,12 +1,21 @@ +Import('env') import os.path -Import('env') - -files = ['AttributeTable.h', 'graph_base.h', 'Grid3D.h', 'Vector.h', - 'Object.h', 'ObjectPointer.h', 'ObjectContainer.h', 'ParticleGrid.h', - 'kernel_version_info.h', 'constants.h', 'units.h', - 'RefCountedObject.h', 'utility.h'] - -# Install the include files: -includedir = os.path.join(env['includedir'], 'IMP', 'internal') +files=[ + 'AttributeTable.h', + 'bbox_nbl_helpers.h', + 'constants.h', + 'graph_base.h', + 'Grid3D.h', + 'kernel_version_info.h', + 'ObjectContainer.h', + 'Object.h', + 'ObjectPointer.h', + 'ParticleGrid.h', + 'RefCountedObject.h', + 'units.h', + 'utility.h', + 'Vector.h', + ] +includedir = os.path.join(env['includedir'], 'IMP', 'internal' ) inst = env.Install(includedir, files) env.Alias('install', inst) Index: kernel/include/IMP/internal/Vector.h =================================================================== --- kernel/include/IMP/internal/Vector.h (revision 541) +++ kernel/include/IMP/internal/Vector.h (working copy) @@ -30,33 +30,41 @@ template class Vector: public std::vector { + typedef Vector This; public: typedef std::vector P; Vector(){} + template Vector(It b, It e) { insert(P::end(), b, e); } + Vector(const P &p) { insert(P::end(), p.begin(), p.end()); } + ~Vector(){clear();} + const D& operator[](unsigned int i) const { IMP_check(i < P::size(), "Index " << i << " out of range", IndexException("")); return P::operator[](i); } + D& operator[](unsigned int i) { IMP_check(i < P::size(), "Index " << i << " out of range", IndexException("")); return P::operator[](i); } + void erase(typename P::iterator it) { unref(*it); P::erase(it); } + void erase(typename P::iterator b, typename P::iterator e) { for (typename P::iterator c= b; c != e; ++c) { @@ -64,21 +72,25 @@ } P::erase(b,e); } + unsigned int push_back(D d) { ref(d); P::push_back(d); return P::size()-1; } + void pop_back() { unref(P::back()); P::pop_back(); } + void clear() { for (typename P::iterator it= P::begin(); it != P::end(); ++it) { unref(*it); } P::clear(); } + template void insert(typename P::iterator it, It b, It e) { for (It c= b; c != e; ++c) { @@ -87,6 +99,11 @@ } P::insert(it, b, e); } + + + unsigned int get_memory_usage() const { + return sizeof(This) + sizeof(D) * P::capacity(); + } }; Index: kernel/include/IMP/internal/units.h =================================================================== --- kernel/include/IMP/internal/units.h (revision 541) +++ kernel/include/IMP/internal/units.h (working copy) @@ -96,6 +96,10 @@ return ExponentialNumber(v_*o.v_); } + This operator-() const { + return This(-v_); + } + void show(std::ostream &out) const { std::ios::fmtflags of = out.flags(); out << setiosflags(std::ios::fixed) << v_; @@ -180,6 +184,17 @@ bool operator<(MKSUnit o) const { return v_ < o.v_; } + + //! MKS units can be compared if their types match + template + bool operator<=(MKSUnit o) const { + return v_ <= o.v_; + } + //! MKS units can be compared if their types match + template + bool operator>=(MKSUnit o) const { + return v_ >= o.v_; + } //! only works with matching types template bool operator>(MKSUnit o) const { @@ -228,6 +243,10 @@ return This(v_-o.v_); } + This operator-() const { + return This(-v_); + } + //! Get a string version of the name of the units std::string get_name() const { int m=M; @@ -275,6 +294,22 @@ return MKSUnit(IMP::square(o.get_value())); } +//! They have to be the same EXP to make the return type make sense +template +MKSUnit +min(MKSUnit a, MKSUnit b) { + if (a +MKSUnit +max(MKSUnit a, MKSUnit b) { + if (a>b) return a; + else return b; +} + typedef MKSUnit<-10, 1, 0, 0, 0> Angstrom; typedef MKSUnit<-20, 2, 0, 0, 0> SquaredAngstrom; typedef MKSUnit<0, 2, 1, 0, -2> Joule; @@ -334,6 +369,9 @@ //! Get the value ignoring the exponent double get_value() const {return v_.get_value();} + This operator-() const { + return This(-v_); + } //! Get the value ignoring the exponent V get_exponential_value() const {return v_;} Index: kernel/include/IMP/internal/utility.h =================================================================== --- kernel/include/IMP/internal/utility.h (revision 541) +++ kernel/include/IMP/internal/utility.h (working copy) @@ -26,17 +26,24 @@ } }; -//! Remove all the inactive particles from a vector of particles -inline void remove_inactive_particles(Particles &ps) +//! Remove all particles from a range in a container +inline void remove_inactive_particles(Particles &ps, + Particles::iterator b, + Particles::iterator e) { - ps.erase(std::remove_if(ps.begin(), ps.end(), - internal::IsInactiveParticle()), ps.end()); - for (Particles::iterator c = ps.begin(); c != ps.end(); ++c) { + ps.erase(std::remove_if(b, e, internal::IsInactiveParticle()), e); + for (Particles::iterator c= b; c!= e; ++c) { (*c)->assert_is_valid(); IMP_assert((*c)->get_is_active(), "Did not remove inactive particle"); } } +//! Remove all the inactive particles from a vector of particles +inline void remove_inactive_particles(Particles &ps) +{ + remove_inactive_particles(ps, ps.begin(), ps.end()); +} + } // namespace internal } // namespace IMP Index: kernel/include/IMP/internal/Object.h =================================================================== --- kernel/include/IMP/internal/Object.h (revision 541) +++ kernel/include/IMP/internal/Object.h (working copy) @@ -35,7 +35,7 @@ { protected: Object(); - ~Object(); + virtual ~Object(); public: //! Throw an assertion if the object has been freed @@ -54,6 +54,10 @@ IMP_assert(count_ ==1, "Too many unrefs on object"); --count_; } + + virtual unsigned int get_memory_usage() const { + return sizeof(Object); + } protected: int count_; private: Index: kernel/include/IMP/internal/bbox_nbl_helpers.h =================================================================== --- kernel/include/IMP/internal/bbox_nbl_helpers.h (revision 0) +++ kernel/include/IMP/internal/bbox_nbl_helpers.h (revision 0) @@ -0,0 +1,40 @@ +/** + * \file bbox_nbl_helpers.h + * \brief Helpers for the bbox NBL. + * + * Copyright 2007-8 Sali Lab. All rights reserved. + * + */ + +#ifndef __IMP_BBOX_NBL_HELPERS_H +#define __IMP_BBOX_NBL_HELPERS_H + +#include "../base_types.h" +#include "../Particle.h" + +namespace IMP +{ + +namespace internal +{ + + +// defined in NonbondedListScoreState.h +class NBLAddPairIfNonbonded; + +/* scan for overlapping bounding boxes within the set of particles with + radii. The radii are assume to be 0 if missing. Slack and cutoff/2 are + added to the radii. ap(Paritcle*, Particle*) is called for each + intersecting pair*/ +void bbox_scan(const Particles &ps, FloatKey rk, Float slack, Float cutoff, + const NBLAddPairIfNonbonded &ap); +/* bipartite version of above */ +void bipartite_bbox_scan(const Particles &ps0, const Particles &ps1, + FloatKey rk, Float slack, Float cutoff, + const NBLAddPairIfNonbonded &ap); + +} + +} + +#endif Index: kernel/include/IMP/internal/ObjectContainer.h =================================================================== --- kernel/include/IMP/internal/ObjectContainer.h (revision 541) +++ kernel/include/IMP/internal/ObjectContainer.h (working copy) @@ -8,6 +8,8 @@ #ifndef __IMP_OBJECT_CONTAINER_H #define __IMP_OBJECT_CONTAINER_H +#include "Object.h" + #include #include @@ -31,7 +33,11 @@ template class ObjectContainer: public std::vector { + typedef ObjectContainer This; + std::vector free_; + + struct OK { bool operator()(const O*a) const { return a != NULL; @@ -124,7 +130,7 @@ } template void insert(iterator c, It b, It e) { -#ifndef NDEBUG + IMP_BEGIN_CHECK(CHEAP); for (It cc= b; cc != e; ++cc) { IMP_CHECK_OBJECT(*cc); for (typename Vector::const_iterator it= Vector::begin(); @@ -132,8 +138,8 @@ IMP_assert(*it != *cc, "IMP Containers can only have one copy of " << " each object"); } - } -#endif + } + IMP_END_CHECK; for (It cc= b; cc != e; ++cc) { ref(*cc); } @@ -144,6 +150,15 @@ Vector::insert(c.base(), b, e); } + unsigned int get_memory_usage() const { + unsigned int owned=0; + for (const_iterator it= begin(); it != end(); ++it) { + owned += static_cast(*it)->get_memory_usage(); + } + return owned + sizeof(This) + sizeof(void*) * Vector::capacity() + + free_.capacity() *sizeof(unsigned int); + } + }; } // namespace internal Index: kernel/include/IMP/internal/ParticleGrid.h =================================================================== --- kernel/include/IMP/internal/ParticleGrid.h (revision 541) +++ kernel/include/IMP/internal/ParticleGrid.h (working copy) @@ -25,26 +25,18 @@ // don't need ref counting since mc_ has the same set of points typedef internal::Grid3D Grid; Grid grid_; - internal::ObjectPointer mc_; Float target_voxel_side_; - bool grid_valid_; - void build_grid(); + void build_grid(const Particles &ps); void audit_particles(const Particles &ps) const; void add_particle_to_grid(Particle *p); public: ParticleGrid(); //! suggested grid edge size. - ParticleGrid(Float sz); + ParticleGrid(Float sz, const Particles &ps); Float get_voxel_size() const {return target_voxel_side_;} - void add_particles(const Particles &ps); - void add_particle(Particle *p); - void clear_particles(); - const Particles& get_particles() const {return mc_->get_particles();} - bool update(); - void show(std::ostream &out) const; typedef Grid::VirtualIndex VirtualIndex; Index: kernel/include/IMP/internal/constants.h =================================================================== --- kernel/include/IMP/internal/constants.h (revision 541) +++ kernel/include/IMP/internal/constants.h (working copy) @@ -28,6 +28,8 @@ //! the default temperature extern IMPDLLEXPORT const Kelvin T; +static const double PI= 3.1415926535897931; + } // namespace internal } // namespace IMP Index: kernel/include/IMP/Key.h =================================================================== --- kernel/include/IMP/Key.h (revision 541) +++ kernel/include/IMP/Key.h (working copy) @@ -14,7 +14,6 @@ #include #include - /** \internal \page keys How Keys work in IMP @@ -91,7 +90,6 @@ { typedef std::map Map; typedef std::vector RMap; - void show(std::ostream &out= std::cout) const; KeyData(); void assert_is_initialized() const; @@ -119,6 +117,30 @@ template struct KeyDataHolder {}; +/** + Define a new key type. There must be an accompanying IMP_DEFINE_KEY_TYPE. + \note The traits class includes a hack to try to check that the data + is initialized before the tables are used. + + \param[in] Name The name for the new type. + \param[in] Tag A class which is unique for this type. For attributes + use the type of the attributes. For other Keys, declare an empty + class with a unique name and use it. + */ +#define IMP_DECLARE_KEY_TYPE(Name, Tag) \ + IMP_KEY_DATA_HOLDER(Name, Tag); \ + typedef Key Name; \ + typedef std::vector Name##s + + +/** This must occur in exactly one .o in the internal namespace. Should + be used in the IMP namespace.*/ +#define IMP_DEFINE_KEY_TYPE(Name, Tag) \ + namespace internal { \ + KeyData KeyDataHolder::data; \ + } + + } // namespace internal Index: kernel/include/IMP/exception.h =================================================================== --- kernel/include/IMP/exception.h (revision 541) +++ kernel/include/IMP/exception.h (working copy) @@ -18,6 +18,9 @@ namespace IMP { + +#ifndef SWIG + //! The general base class for IMP exceptions /** This way we can catch IMP exceptions without getting memory allocation errors and everything. And it enforces having a description. @@ -45,11 +48,13 @@ destroy(); } Exception(const Exception &o) {copy(o);} + Exception &operator=(const Exception &o) { destroy(); copy(o); return *this; } + private: void destroy() { if (str_ != NULL) { @@ -109,6 +114,34 @@ }; +#endif // SWIG + + + +//! Determine the level of runtime checks performed +/** NONE means that minimial checks are used. CHEAP + means that only constant time checks are performed + and with EXPENSIVE non-linear time checks will be run. + */ +enum CheckLevel {NONE=0, CHEAP=1, EXPENSIVE=2}; + + +//! Control runtime checks in the code +/** The default level of checks is CHEAP. + */ +IMPDLLEXPORT void set_check_level(CheckLevel tf); + +//! Get the current audit mode +IMPDLLEXPORT CheckLevel get_check_level(); + + + +#define IMP_BEGIN_CHECK(level)\ + if (level <= ::IMP::get_check_level()) { do {} while(false) + +#define IMP_END_CHECK } do {} while(false) + + namespace internal { @@ -129,7 +162,7 @@ } // namespace IMP -#ifndef NDEBUG +//#ifndef NDEBUG //! An assertion for IMP. An IMP::ErrorException will be thrown. /** Since it is a debug-only check and no attempt should be made to @@ -139,16 +172,16 @@ \param[in] message Write this message if the assertion fails. \ingroup assert */ -#define IMP_assert(expr, message) \ - do { \ - if (!(expr)) { \ - IMP_ERROR(message); \ - IMP::internal::assert_fail(); \ - } \ +#define IMP_assert(expr, message) \ + do { \ + if (IMP::get_check_level() >= IMP::EXPENSIVE && !(expr)) { \ + IMP_ERROR(message); \ + IMP::internal::assert_fail(); \ + } \ } while(false) -#else +/*#else #define IMP_assert(expr, message) -#endif +#endif*/ //! A runtime check for IMP. /** \param[in] expr The assertion expression. @@ -156,12 +189,12 @@ \param[in] exception Throw the object constructed by this expression. \ingroup assert */ -#define IMP_check(expr, message, exception) \ - do { \ - if (!(expr)) { \ - IMP_ERROR(message); \ - throw exception; \ - } \ +#define IMP_check(expr, message, exception) \ + do { \ + if (IMP::get_check_level() >= IMP::CHEAP && !(expr)) { \ + IMP_ERROR(message); \ + throw exception; \ + } \ } while (false) //! A runtime failure for IMP. Index: kernel/include/IMP/unary_functions/WormLikeChain.h =================================================================== --- kernel/include/IMP/unary_functions/WormLikeChain.h (revision 0) +++ kernel/include/IMP/unary_functions/WormLikeChain.h (revision 0) @@ -0,0 +1,115 @@ +/** + * \file WormLikeChainEnergy.h \brief WormLikeChainEnergy function. + * + * Copyright 2007-8 Sali Lab. All rights reserved. + */ + +#ifndef __IMP_WORMLIKECHAIN_H +#define __IMP_WORMLIKECHAIN_H + +#include "../UnaryFunction.h" +#include "../internal/constants.h" +#include "../internal/units.h" + +namespace IMP +{ + +//! Worm-like-chain energy for polymer chains +/** This function implements one polymer force/extension curve. It + is a physical energy and all the inputs are in Angstroms + and the outputs in kCal/mol (/angstrom). + + \note The actual worm like chain force blows up as the extension approaches + the maximim, since that makes optimization problematic, we + just linearly extend the force after 99% of maximum. + \ingroup unaryf + */ +class WormLikeChain : public UnaryFunction +{ +public: + //! Define the energy term + /** \param[in] l_max Is the maximum length of the chain in angstroms + \param[in] lp is the persistence length also in angstroms + */ + WormLikeChain(Float l_max, Float lp) : lmax_(l_max), lp_(lp) {} + + virtual ~WormLikeChain() {} + + //! Calculate the WormLikeChain energy given the length + /** \param[in] l Current length in Angstroms + \return Energy in kCal/Mol + */ + virtual Float evaluate(Float lf) { + static const internal::PicoJoule zero=eval(internal::Angstrom(0)); + internal::Angstrom l(lf); + if (l < internal::Angstrom(0)) l=internal::Angstrom(0); + internal::PicoJoule ret; + if (l < cutoff()) { + ret= (eval(l) - zero); + } else { + internal::PicoJoule springterm=(l-cutoff())*cderiv(cutoff()); + ret= (eval(cutoff())+ springterm -zero); + } + /*std::cout << "Return is " << ret <<" " << l << " " << lp_ << " " + << lmax_ << std::endl;*/ + return internal::KCalPerMol(convert_to_kcal(ret)).get_value(); + } + + //! Calculate the WormLikeChain energy given the length + /** \param[in] l Current length in Angstroms + \param[out] deriv The force in kCal/Angstrom Mol + \return Score + */ + virtual Float evaluate_deriv(Float fl, Float& deriv) { + internal::Angstrom l(fl); + if (l < internal::Angstrom(0)) l=internal::Angstrom(0); + internal::PicoNewton doubled; + if (l < cutoff()) { + doubled= cderiv(l); + } else { + doubled= cderiv(cutoff()); + IMP_LOG(VERBOSE, "Overstretched " << cderiv(cutoff()) << " " << doubled + << " " << l << " " << lmax_ << " " << cutoff() + << std::endl); + } + //std::cout << "Force is " << doubled << " for length " << l << std::endl; + // convert from picoNewton + deriv = internal::KCalPerAMol(internal::convert_to_kcal(doubled)) + .get_value(); + //std::cout << "Which converts to " << d << std::endl; + return evaluate(fl); + } + + void show(std::ostream &out=std::cout) const { + out << "WormLikeChain " << lmax_ << " " << lp_ << std::endl; + } + +protected: + // have to name it this since the SWIG typemap hack takes the string deriv + internal::PicoNewton cderiv(internal::Angstrom l) const { + internal::PicoNewton pn= internal::KB*internal::T + /lp_*(internal::Scalar(.25)/ square(internal::Scalar(1)-l/lmax_) + -internal::Scalar(.25)+l/lmax_); + return pn; + } + + internal::PicoJoule eval(internal::Angstrom m) const { + internal::PicoJoule J + = internal::KB*internal::T/lp_*(internal::Scalar(.25)*square(lmax_) + /(lmax_-m) + -m*internal::Scalar(.25) + +internal::Scalar(.5)*square(m) + /lmax_); + return J; + } + + internal::Angstrom cutoff() const { + return internal::Scalar(.99)*lmax_; + } + + internal::Angstrom lmax_, lp_; +}; + +} // namespace IMP + +#endif /* __IMP_HARMONIC_H */ Index: kernel/include/IMP/unary_functions/Harmonic.h =================================================================== --- kernel/include/IMP/unary_functions/Harmonic.h (revision 541) +++ kernel/include/IMP/unary_functions/Harmonic.h (working copy) @@ -8,6 +8,7 @@ #define __IMP_HARMONIC_H #include "../UnaryFunction.h" +#include "../utility.h" namespace IMP { Index: kernel/include/IMP/unary_functions/SConscript =================================================================== --- kernel/include/IMP/unary_functions/SConscript (revision 541) +++ kernel/include/IMP/unary_functions/SConscript (working copy) @@ -1,10 +1,15 @@ +Import('env') import os.path -Import('env') - -files = ['Harmonic.h', 'HarmonicLowerBound.h', 'HarmonicUpperBound.h', - 'OpenCubicSpline.h', 'ClosedCubicSpline.h', 'Cosine.h', 'Linear.h'] - -# Install the include files: -includedir = os.path.join(env['includedir'], 'IMP', 'unary_functions') +files=[ + 'ClosedCubicSpline.h', + 'Cosine.h', + 'Harmonic.h', + 'HarmonicLowerBound.h', + 'HarmonicUpperBound.h', + 'Linear.h', + 'OpenCubicSpline.h', + 'WormLikeChain.h', + ] +includedir = os.path.join(env['includedir'], 'IMP', 'unary_functions' ) inst = env.Install(includedir, files) env.Alias('install', inst) Index: kernel/include/IMP/SConscript =================================================================== --- kernel/include/IMP/SConscript (revision 541) +++ kernel/include/IMP/SConscript (working copy) @@ -1,26 +1,39 @@ +Import('env') import os.path -Import('env') - -files = ['base_types.h', 'random.h', 'Index.h', 'Model.h', - 'Particle.h', 'ScoreState.h', 'OptimizerState.h', 'IMP_config.h', - 'log.h', 'DerivativeAccumulator.h', - 'Key.h', 'utility.h', 'Restraint.h', 'Optimizer.h', - 'DecoratorBase.h', 'Vector3D.h', - 'UnaryFunction.h', 'PairScore.h', 'SingletonScore.h', 'macros.h', - 'TripletScore.h', 'exception.h', 'VersionInfo.h'] - -# Install the include files: -includedir = os.path.join(env['includedir'], 'IMP') +files=[ + 'base_types.h', + 'DecoratorBase.h', + 'DerivativeAccumulator.h', + 'exception.h', + 'IMP_config.h', + 'Index.h', + 'Key.h', + 'log.h', + 'macros.h', + 'Model.h', + 'Optimizer.h', + 'OptimizerState.h', + 'PairScore.h', + 'Particle.h', + 'random.h', + 'Restraint.h', + 'ScoreState.h', + 'SingletonScore.h', + 'TripletScore.h', + 'UnaryFunction.h', + 'utility.h', + 'Vector3D.h', + 'VersionInfo.h', + ] +includedir = os.path.join(env['includedir'], 'IMP' ) inst = env.Install(includedir, files) env.Alias('install', inst) - -# Subdirectories -SConscript('restraints/SConscript') -SConscript('optimizers/SConscript') -SConscript('decorators/SConscript') -SConscript('internal/SConscript') -SConscript('unary_functions/SConscript') -SConscript('pair_scores/SConscript') -SConscript('singleton_scores/SConscript') -SConscript('triplet_scores/SConscript') -SConscript('score_states/SConscript') +SConscript( 'decorators/SConscript' ) +SConscript( 'internal/SConscript' ) +SConscript( 'optimizers/SConscript' ) +SConscript( 'pair_scores/SConscript' ) +SConscript( 'restraints/SConscript' ) +SConscript( 'score_states/SConscript' ) +SConscript( 'singleton_scores/SConscript' ) +SConscript( 'triplet_scores/SConscript' ) +SConscript( 'unary_functions/SConscript' ) Index: kernel/include/IMP/Particle.h =================================================================== --- kernel/include/IMP/Particle.h (revision 541) +++ kernel/include/IMP/Particle.h (working copy) @@ -321,7 +321,17 @@ return particle_indexes_.attribute_keys_end(); } + //! Create another particle in the same model with all the same attributes + Particle *clone() const; + + unsigned int get_memory_usage() const { + return sizeof(Particle)+float_indexes_.get_heap_memory_usage() + + int_indexes_.get_heap_memory_usage() + + string_indexes_.get_heap_memory_usage() + + particle_indexes_.get_heap_memory_usage(); + } + protected: void zero_derivatives(); Index: kernel/include/IMP/triplet_scores/SConscript =================================================================== --- kernel/include/IMP/triplet_scores/SConscript (revision 541) +++ kernel/include/IMP/triplet_scores/SConscript (working copy) @@ -1,9 +1,8 @@ -import os.path Import('env') - -files = ['AngleTripletScore.h'] - -# Install the include files: -includedir = os.path.join(env['includedir'], 'IMP', 'triplet_scores') +import os.path +files=[ + 'AngleTripletScore.h', + ] +includedir = os.path.join(env['includedir'], 'IMP', 'triplet_scores' ) inst = env.Install(includedir, files) env.Alias('install', inst) Index: kernel/include/IMP/ScoreState.h =================================================================== --- kernel/include/IMP/ScoreState.h (revision 541) +++ kernel/include/IMP/ScoreState.h (working copy) @@ -46,17 +46,19 @@ ScoreState(std::string name=std::string()); virtual ~ScoreState(); - //! Update the state given the current state of the model. - /** This is also called prior to every calculation of the model score. + //! Update if needed + /** The protected update() method will be called if the iteration + count has not yet been seen. */ - virtual void update() = 0; + void before_evaluate(unsigned int iteration); - //! Do any necessary updates after the model score is calculated. - /** \param[in] accpt The object used to scale derivatives in the score - calculation, or NULL if derivatives were not requested. + //! Do post evaluation work if needed + /** The protected after_evaluate method will be called if needed. */ - virtual void after_evaluate(DerivativeAccumulator *accpt) {} + void after_evaluate(unsigned int iteration, + DerivativeAccumulator *accpt); + //! Show the ScoreState /** The output of show may take several lines and should end in a newline. */ @@ -82,10 +84,44 @@ "Must call set_model before get_model on state"); return model_.get(); } + protected: + //! Update the state given the current state of the model. + /** This is also called prior to every calculation of the model score. + It should be implemented by ScoreStates in order to provide functionality. + + \note This can't have the same name as the public function due to the + way C++ handles overloading and name lookups--if only one is implemented + in the child class it will only find that one. + */ + virtual void do_before_evaluate() = 0; + + //! Do any necessary updates after the model score is calculated. + /** \param[in] accpt The object used to scale derivatives in the score + calculation, or NULL if derivatives were not requested. + */ + virtual void do_after_evaluate(DerivativeAccumulator *accpt) {} + + //! Get the current iteration count value. + /** The value is updated before update() is called + */ + unsigned int get_before_evaluate_iteration() const { + return update_iteration_; + } + + //! Get the current after_evaluate iteration count value. + /** The value is updated before after_evaluate() is called + */ + unsigned int get_after_evaluate_iteration() const { + return after_iteration_; + } + + private: + + unsigned int update_iteration_; + unsigned int after_iteration_; // all of the particle data internal::ObjectPointer model_; - std::string name_; }; Index: kernel/include/IMP/macros.h =================================================================== --- kernel/include/IMP/macros.h (revision 541) +++ kernel/include/IMP/macros.h (working copy) @@ -188,10 +188,26 @@ /** \return version and authorship information */ \ virtual IMP::VersionInfo get_version_info() const { return version_info; } -//! See IMP_OPTIMIZER_STATE -#define IMP_SCORE_STATE(version_info)\ - IMP_OPTIMIZER_STATE(version_info) +//! Define the basics needed for an ScoreState +/** This macro declares the required functions + - void do_before_update() + - void show(std::ostream &out) const + and defines the function + - get_version_info + \param[in] version_info The version info object to return. +*/ +#define IMP_SCORE_STATE(version_info) \ +protected: \ + /** update the state*/ \ + virtual void do_before_evaluate(); \ +public: \ + /** write information about the state to the stream*/ \ + virtual void show(std::ostream &out=std::cout) const; \ + /** \return version and authorship information */ \ + virtual IMP::VersionInfo get_version_info() const { return version_info; } + + //! Use the swap_with member function to swap two objects #define IMP_SWAP(name) \ inline void swap(name &a, name &b) { \ @@ -251,6 +267,7 @@ private: \ /** \internal */ \ Container lcname##_vector_; \ + unsigned int get_##lcname##s_memory_usage() const ; \ protection: /** \internal @@ -272,7 +289,7 @@ Data obj= lcname##_vector_[osz+i]; \ IndexType index(osz+i); \ Init_obj; \ - if (false) std::cout << *obj << index; \ + if (false) std::cout << *obj << index; \ } \ Onchanged; \ } \ @@ -281,6 +298,9 @@ lcname##_vector_.clear(); \ Onchanged; \ } \ + unsigned int Class::get_##lcname##s_memory_usage() const { \ + return lcname##_vector_.get_memory_usage(); \ + } \ //! Use this to add a container of IMP objects @@ -291,6 +311,8 @@ \param[in] Ucname The name of the type in uppercase \param[in] lcname The name of the type in lower case \param[in] Data The type of the data to store. + \param[in] Onchanged Code to get executed every time the container + changes \note the type Ucnames must be declared and be a vector of Data. @@ -376,7 +398,6 @@ IMP_CONTAINER. \param[in] init Code to modify the passed in object. The object is obj its index index. - \param[in] onchanged Code to execute when the container is changed. \param[in] onremove Code to execute when an object is removed. The object being removed is obj. */ Index: kernel/include/SConscript =================================================================== --- kernel/include/SConscript (revision 541) +++ kernel/include/SConscript (working copy) @@ -1,11 +1,9 @@ Import('env') - -files = ['IMP.h'] - -# Install the include files: +import os.path +files=[ + 'IMP.h', + ] includedir = env['includedir'] inst = env.Install(includedir, files) env.Alias('install', inst) - -# Subdirectories -SConscript('IMP/SConscript') +SConscript( 'IMP/SConscript' ) Index: kernel/include/IMP.h =================================================================== --- kernel/include/IMP.h (revision 541) +++ kernel/include/IMP.h (working copy) @@ -1,84 +1,99 @@ /** - * \file IMP.h \brief IMP, an Integrative Modeling Platform. - * - * Copyright 2007-8 Sali Lab. All rights reserved. - * - */ - +* \file IMP.h \brief IMP, an Integrative Modeling Platform. +* +* Copyright 2007-8 Sali Lab. All rights reserved. +* +*/ #ifndef __IMP_H #define __IMP_H - -#include "IMP/IMP_config.h" -#include "IMP/log.h" -#include "IMP/random.h" -#include "IMP/base_types.h" -#include "IMP/Particle.h" -#include "IMP/Optimizer.h" -#include "IMP/Restraint.h" -#include "IMP/exception.h" -#include "IMP/UnaryFunction.h" -#include "IMP/unary_functions/Harmonic.h" -#include "IMP/unary_functions/HarmonicLowerBound.h" -#include "IMP/unary_functions/HarmonicUpperBound.h" -#include "IMP/unary_functions/OpenCubicSpline.h" -#include "IMP/unary_functions/ClosedCubicSpline.h" -#include "IMP/unary_functions/Cosine.h" -#include "IMP/unary_functions/Linear.h" -#include "IMP/Model.h" -#include "IMP/PairScore.h" -#include "IMP/SingletonScore.h" -#include "IMP/TripletScore.h" -#include "IMP/Vector3D.h" -#include "IMP/VersionInfo.h" -#include "IMP/decorators/HierarchyDecorator.h" -#include "IMP/decorators/MolecularHierarchyDecorator.h" -#include "IMP/decorators/NameDecorator.h" -#include "IMP/decorators/AtomDecorator.h" -#include "IMP/decorators/ResidueDecorator.h" -#include "IMP/decorators/XYZDecorator.h" -#include "IMP/decorators/bond_decorators.h" -#include "IMP/optimizers/SteepestDescent.h" -#include "IMP/optimizers/ConjugateGradients.h" -#include "IMP/optimizers/MolecularDynamics.h" -#include "IMP/optimizers/MonteCarlo.h" -#include "IMP/optimizers/Mover.h" -#include "IMP/optimizers/MoverBase.h" -#include "IMP/optimizers/movers/BallMover.h" -#include "IMP/optimizers/movers/NormalMover.h" -#include "IMP/optimizers/states/VRMLLogOptimizerState.h" -#include "IMP/optimizers/states/CMMLogOptimizerState.h" -#include "IMP/optimizers/states/VelocityScalingOptimizerState.h" -#include "IMP/pair_scores/DistancePairScore.h" -#include "IMP/pair_scores/SphereDistancePairScore.h" -#include "IMP/singleton_scores/DistanceToSingletonScore.h" -#include "IMP/singleton_scores/AttributeSingletonScore.h" -#include "IMP/triplet_scores/AngleTripletScore.h" -#include "IMP/restraints/RestraintSet.h" -#include "IMP/restraints/DistanceRestraint.h" -#include "IMP/restraints/AngleRestraint.h" -#include "IMP/restraints/DihedralRestraint.h" -#include "IMP/restraints/ConnectivityRestraint.h" -#include "IMP/restraints/NonbondedRestraint.h" -#include "IMP/restraints/BondDecoratorRestraint.h" -#include "IMP/restraints/SingletonListRestraint.h" -#include "IMP/restraints/PairListRestraint.h" -#include "IMP/restraints/TripletChainRestraint.h" -#include "IMP/restraints/PairChainRestraint.h" -#include "IMP/score_states/BipartiteNonbondedListScoreState.h" -#include "IMP/score_states/MaxChangeScoreState.h" -#include "IMP/score_states/NonbondedListScoreState.h" -#include "IMP/score_states/BondedListScoreState.h" -#include "IMP/score_states/BondDecoratorListScoreState.h" -#include "IMP/score_states/AllNonbondedListScoreState.h" -#include "IMP/score_states/QuadraticBipartiteNonbondedListScoreState.h" -#include "IMP/score_states/QuadraticAllNonbondedListScoreState.h" -#include "IMP/score_states/QuadraticNonbondedListScoreState.h" -#include "IMP/score_states/GravityCenterScoreState.h" -#include "IMP/score_states/ChainBondedListScoreState.h" - - -/** - \namespace IMP The IMP namespace. - */ - +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #endif /* __IMP_H */ Index: kernel/doc/internal/doxygen.conf =================================================================== --- kernel/doc/internal/doxygen.conf (revision 541) +++ kernel/doc/internal/doxygen.conf (working copy) @@ -490,7 +490,7 @@ # excluded from the INPUT source files. This way you can easily exclude a # subdirectory from a directory tree whose root is specified with the INPUT tag. -EXCLUDE = +EXCLUDE = ../../include/IMP/internal/ ../../src/internal # The EXCLUDE_SYMLINKS tag can be used select whether or not files or # directories that are symbolic links (a Unix filesystem feature) are excluded @@ -677,7 +677,7 @@ # each generated HTML page. If it is left blank doxygen will generate a # standard footer. -HTML_FOOTER = +HTML_FOOTER = footer.html # The HTML_STYLESHEET tag can be used to specify a user-defined cascading # style sheet that is used by each HTML page. It can be used to Index: kernel/doc/examples/chain.py =================================================================== --- kernel/doc/examples/chain.py (revision 541) +++ kernel/doc/examples/chain.py (working copy) @@ -7,6 +7,10 @@ # another. #IMP.set_log_level(IMP.VERBOSE) + +# turn on all debugging checks +IMP.set_check_level(IMP.EXPENSIVE) + radius =1.0 rk= IMP.FloatKey("radius") m= IMP.Model() @@ -34,7 +38,7 @@ p.show() # Set up the nonbonded list -nbl= IMP.AllNonbondedListScoreState(rk, chain) +nbl= IMP.AllNonbondedListScoreState(1, rk, chain) nbli= m.add_score_state(nbl) # This ScoreState uses the bonds constructed above to restrain bl= IMP.BondDecoratorListScoreState(chain) @@ -43,7 +47,7 @@ # Set up excluded volume ps= IMP.SphereDistancePairScore(IMP.HarmonicLowerBound(0,1), rk) -evr= IMP.NonbondedRestraint(ps, nbl, 1) +evr= IMP.NonbondedRestraint(ps, nbl) evri= m.add_restraint(evr) # Restraint for bonds @@ -72,7 +76,7 @@ # Write the progression of states as the system is optimized to # the files state.000.vrml, state.001.vrml etc. vrml= IMP.VRMLLogOptimizerState("state.%03d.vrml", chain) -vrml.set_radius(rk) +vrml.set_radius_key(rk) vrml.update() vrml.set_skip_steps(100) o.add_optimizer_state(vrml) Index: kernel/src/exception.cpp =================================================================== --- kernel/src/exception.cpp (revision 0) +++ kernel/src/exception.cpp (revision 0) @@ -0,0 +1,27 @@ +/** + * \file exception.cpp \brief Check handling. + * + * Copyright 2007-8 Sali Lab. All rights reserved. + * + */ + +#include "IMP/exception.h" + +namespace IMP +{ +static CheckLevel check_mode= +#ifdef NDEBUG + CHEAP; +#else + EXPENSIVE; +#endif + +void set_check_level(CheckLevel cm) { + check_mode= cm; +} + +CheckLevel get_check_level() { + return check_mode; +} + +} // namespace IMP Index: kernel/src/singleton_scores/SConscript =================================================================== --- kernel/src/singleton_scores/SConscript (revision 541) +++ kernel/src/singleton_scores/SConscript (working copy) @@ -1,6 +1,5 @@ Import('env') - -files = ['DistanceToSingletonScore.cpp', 'AttributeSingletonScore.cpp'] - -files = [File(x) for x in files] +files=[] +files.append(File( 'AttributeSingletonScore.cpp' )) +files.append(File( 'DistanceToSingletonScore.cpp' )) Return('files') Index: kernel/src/singleton_scores/DistanceToSingletonScore.cpp =================================================================== --- kernel/src/singleton_scores/DistanceToSingletonScore.cpp (revision 541) +++ kernel/src/singleton_scores/DistanceToSingletonScore.cpp (working copy) @@ -6,8 +6,9 @@ */ #include "IMP/singleton_scores/DistanceToSingletonScore.h" +#include "IMP/UnaryFunction.h" +#include "IMP/Particle.h" #include "IMP/decorators/XYZDecorator.h" -#include "IMP/UnaryFunction.h" namespace IMP { @@ -16,6 +17,7 @@ { static const Float MIN_DISTANCE = .00001; + Float evaluate_distance_to_singleton_score(const Vector3D &v, Particle *b, DerivativeAccumulator *da, @@ -75,7 +77,7 @@ void DistanceToSingletonScore::show(std::ostream &out) const { - out << "DistanceToSingletScore using "; + out << "DistanceToSingletonScore using "; f_->show(out); } Index: kernel/src/Model.cpp =================================================================== --- kernel/src/Model.cpp (revision 541) +++ kernel/src/Model.cpp (working copy) @@ -19,6 +19,7 @@ //! Constructor Model::Model() { + iteration_=0; } @@ -59,7 +60,7 @@ for (ScoreStateIterator it = score_states_begin(); it != score_states_end(); ++it) { IMP_CHECK_OBJECT(*it); - (*it)->update(); + (*it)->before_evaluate(iteration_); IMP_LOG(VERBOSE, "." << std::flush); } IMP_LOG(TERSE, "End update ScoreStates." << std::endl); @@ -74,6 +75,7 @@ "Begin evaluate restraints " << (calc_derivs?"with derivatives":"without derivatives") << std::endl); + for (RestraintIterator it = restraints_begin(); it != restraints_end(); ++it) { IMP_CHECK_OBJECT(*it); @@ -86,19 +88,38 @@ IMP_LOG(TERSE, "Restraint score is " << tscore << std::endl); score+= tscore; } + IMP_LOG(TERSE, "End evaluate restraints." << std::endl); IMP_LOG(TERSE, "Begin after_evaluate of ScoreStates " << std::endl); + for (ScoreStateIterator it = score_states_begin(); it != score_states_end(); ++it) { IMP_CHECK_OBJECT(*it); - (*it)->after_evaluate(accpt); + (*it)->after_evaluate(iteration_, accpt); IMP_LOG(VERBOSE, "." << std::flush); } + IMP_LOG(TERSE, "End after_evaluate of ScoreStates." << std::endl); IMP_LOG(TERSE, "End Model::evaluate. Final score: " << score << std::endl); + + for (ParticleIterator pit= particles_begin(); + pit != particles_end(); ++pit) { + for (Particle::FloatKeyIterator fkit = (*pit)->float_keys_begin(); + fkit != (*pit)->float_keys_end(); ++fkit) { + Float v= (*pit)->get_value(*fkit); + IMP_check(v==v, "NaN found in particle " << (*pit)->get_index() + << " attribute " << *fkit, + ValueException("NaN found")); + Float d= (*pit)->get_derivative(*fkit); + IMP_check(d==d, "NaN found in particle derivative " << (*pit)->get_index() + << " attribute " << *fkit, + ValueException("NaN found")); + } + } + ++iteration_; return score; } @@ -117,5 +138,11 @@ } +unsigned int Model::get_memory_usage() const { + return sizeof(Model) + get_particles_memory_usage() + + get_score_states_memory_usage() + + get_restraints_memory_usage(); +} + } // namespace IMP Index: kernel/src/Particle.cpp =================================================================== --- kernel/src/Particle.cpp (revision 541) +++ kernel/src/Particle.cpp (working copy) @@ -45,6 +45,16 @@ } +Particle* Particle::clone() const { + Particle *p= new Particle(); + get_model()->add_particle(p); + p->float_indexes_= float_indexes_; + p->int_indexes_= int_indexes_; + p->string_indexes_= string_indexes_; + p->particle_indexes_= particle_indexes_; + return p; +} + void Particle::show(std::ostream& out) const { char* inset = " "; Index: kernel/src/unary_functions/SConscript =================================================================== --- kernel/src/unary_functions/SConscript (revision 541) +++ kernel/src/unary_functions/SConscript (working copy) @@ -1,6 +1,6 @@ Import('env') - -files = ['OpenCubicSpline.cpp', 'ClosedCubicSpline.cpp', 'Cosine.cpp'] - -files = [File(x) for x in files] +files=[] +files.append(File( 'ClosedCubicSpline.cpp' )) +files.append(File( 'Cosine.cpp' )) +files.append(File( 'OpenCubicSpline.cpp' )) Return('files') Index: kernel/src/SConscript =================================================================== --- kernel/src/SConscript (revision 541) +++ kernel/src/SConscript (working copy) @@ -1,35 +1,29 @@ Import('env', 'get_sharedlib_environment') - -# Get an environment suitable for building a shared library: env = get_sharedlib_environment(env, 'IMP_EXPORTS', cplusplus=True) -env.Append(CPPPATH=["#/kernel/include", env['BOOST_CPPPATH']]) - -# Subdirectories: -restraints_files = SConscript('restraints/SConscript') -optimizers_files = SConscript('optimizers/SConscript') -decorators_files = SConscript('decorators/SConscript') -unary_functions_files = SConscript('unary_functions/SConscript') -score_states_files = SConscript('score_states/SConscript') -pair_scores_files = SConscript('pair_scores/SConscript') -singleton_scores_files = SConscript('singleton_scores/SConscript') -triplet_scores_files = SConscript('triplet_scores/SConscript') -score_states_files = SConscript('score_states/SConscript') -internal_files = SConscript('internal/SConscript') - -# Source files -files = ['base_types.cpp', 'Model.cpp', - 'Particle.cpp', 'ScoreState.cpp', - 'OptimizerState.cpp', 'Log.cpp', 'Restraint.cpp', 'Optimizer.cpp', - 'random.cpp', 'Key.cpp' - ] + decorators_files + restraints_files + optimizers_files \ - + unary_functions_files + pair_scores_files + singleton_scores_files \ - + triplet_scores_files + score_states_files + internal_files - -# Build the shared library: -lib = env.SharedLibrary('imp', files) - -# Install the library: +env.Prepend(CPPPATH=["#/kernel/include", env["BOOST_CPPPATH"]]) +env.Append(LIBS= '' ) +files=[] +files.append(File( 'base_types.cpp' )) +files.append(File( 'exception.cpp' )) +files.append(File( 'Key.cpp' )) +files.append(File( 'Log.cpp' )) +files.append(File( 'Model.cpp' )) +files.append(File( 'Optimizer.cpp' )) +files.append(File( 'OptimizerState.cpp' )) +files.append(File( 'Particle.cpp' )) +files.append(File( 'random.cpp' )) +files.append(File( 'Restraint.cpp' )) +files.append(File( 'ScoreState.cpp' )) +files.append(SConscript( 'decorators/SConscript' )) +files.append(SConscript( 'internal/SConscript' )) +files.append(SConscript( 'optimizers/SConscript' )) +files.append(SConscript( 'pair_scores/SConscript' )) +files.append(SConscript( 'restraints/SConscript' )) +files.append(SConscript( 'score_states/SConscript' )) +files.append(SConscript( 'singleton_scores/SConscript' )) +files.append(SConscript( 'triplet_scores/SConscript' )) +files.append(SConscript( 'unary_functions/SConscript' )) +lib = env.SharedLibrary( 'imp' , files) libinst = env.Install(env['libdir'], lib) env.Alias('install', [libinst]) - Return('lib') Index: kernel/src/ScoreState.cpp =================================================================== --- kernel/src/ScoreState.cpp (revision 541) +++ kernel/src/ScoreState.cpp (working copy) @@ -9,6 +9,7 @@ #include "IMP/ScoreState.h" #include +#include namespace IMP { @@ -16,11 +17,29 @@ //! Constructor ScoreState::ScoreState(std::string name) : name_(name) { + update_iteration_= std::numeric_limits::max(); + after_iteration_= std::numeric_limits::max(); model_ = NULL; IMP_LOG(VERBOSE, "ScoreState constructed " << name << std::endl); } +void ScoreState::before_evaluate(unsigned int iter) { + if (update_iteration_ != iter) { + update_iteration_= iter; + do_before_evaluate(); + } +} + + + void ScoreState::after_evaluate(unsigned int iter, + DerivativeAccumulator *da) { + if (after_iteration_ != iter) { + after_iteration_= iter; + do_after_evaluate(da); + } +} + //! Destructor ScoreState::~ScoreState() { Index: kernel/src/score_states/BondDecoratorListScoreState.cpp =================================================================== --- kernel/src/score_states/BondDecoratorListScoreState.cpp (revision 541) +++ kernel/src/score_states/BondDecoratorListScoreState.cpp (working copy) @@ -17,14 +17,14 @@ set_particles(ps); } -void BondDecoratorListScoreState::update() +void BondDecoratorListScoreState::do_before_evaluate() { IMP_LOG(TERSE, "Updating BondDecoratorList for " << ps_.size() << " particles" << std::endl); bonds_.clear(); for (unsigned int i=0; i< ps_.size(); ++i) { if (!ps_[i]->get_is_active()) continue; - BondedDecorator di= BondedDecorator::cast(ps_[i]); + BondedDecorator di(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); Index: kernel/src/score_states/NonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/NonbondedListScoreState.cpp (revision 541) +++ kernel/src/score_states/NonbondedListScoreState.cpp (working copy) @@ -15,24 +15,7 @@ namespace IMP { -NonbondedListScoreState::NonbondedListScoreState(FloatKey rk): rk_(rk) -{ - last_cutoff_=-1; -} - - - - -void NonbondedListScoreState::propagate_particles(const Particles&ps) -{ - clear_nbl(); - for (BondedListScoreStateIterator bli= bonded_lists_begin(); - bli != bonded_lists_end(); ++bli) { - (*bli)->set_particles(ps); - } -} - namespace internal { @@ -45,17 +28,122 @@ } // namespace internal -void NonbondedListScoreState::update() + + + +NonbondedListScoreState +::NonbondedListScoreState(Float cut, + FloatKey rk): rk_(rk), + cutoff_(cut), + nbl_is_valid_(false) { + slack_=20; + next_slack_=slack_; + num_steps_=1; + number_of_updates_=1; + number_of_rebuilds_=0; +} + + + + +NonbondedListScoreState::~NonbondedListScoreState(){ +} + +void NonbondedListScoreState::show_statistics(std::ostream &out) const { + out << "Nonbonded list averaged " + << static_cast(number_of_updates_) + / number_of_rebuilds_ << " steps between rebuilds" << std::endl; +} + + + +Particles NonbondedListScoreState +::particles_with_radius(const Particles &in) const { + Particles ret; + if (rk_== FloatKey()) return ret; + ret.reserve(in.size()); + for (unsigned int i=0; i< in.size(); ++i) { + if (in[i]->has_attribute(rk_)) { + ret.push_back(in[i]); + } + } + return ret; +} + + +/* + Want to minimize rebuild_cost*rebuild_frequence+ 2*nbl_size. + + To do this assume 0 slack has nbl of zero. So the nbl_size + goes as (rate*time_steps)^3= current_size, so + rate= current_size^(1/3)/time_steps + + rebuild_cost/time_steps + 2*(rate*timesteps)^3 + + Min is at 3^(-1/3) (cost)^(1/4)*r^(-3/4); + + Then the slack should be something or another XXXXX + + */ + + +bool NonbondedListScoreState::update(Float mc, Float rebuild_cost) +{ + ++number_of_updates_; IMP_LOG(VERBOSE, "Updating non-bonded list" << std::endl); for (BondedListScoreStateIterator bli= bonded_lists_begin(); bli != bonded_lists_end(); ++bli) { - (*bli)->update(); + (*bli)->before_evaluate(ScoreState::get_before_evaluate_iteration()); } - // if the list is not deleted, we need to scan for inactive particles - nbl_.erase(std::remove_if(nbl_.begin(), nbl_.end(), internal::HasInactive()), - nbl_.end()); + if (nbl_is_valid_) { + /*std::cout << "Rate is " << rate << " target is " << target_steps + << " so slack is " << target_slack << " mc " << mc + << " nbl " << nbl_.size() << " cost " + << rebuild_cost << std::endl;*/ + if (mc > slack_) { + /* Float rate= std::pow(static_cast(nbl_.size()), + .333f)/ num_steps_; + Float target_steps= .6*std::pow(rebuild_cost, .25f) + *std::pow(rate, -.75f); + Float target_slack= (target_steps+1)*mc/num_steps_; + next_slack_= target_slack*.5 + .5*next_slack_; + */ + + /*std::cout << "Killing nbl because " << mc << " + << slack_ << " " << next_slack_ + << " " << num_steps_ << std::endl;*/ + if (num_steps_ < 50) { + //slack_= next_slack_; + } + num_steps_=1; + //next_slack_= std::max(2.0*mc, 2.0*slack_); + set_nbl_is_valid(false); + ++number_of_rebuilds_; + slack_=next_slack_; + /*} else if (num_steps_ > 100) { + //next_slack_=slack_/2.0; + slack_=next_slack_; + num_steps_=1; + set_nbl_is_valid(false); + ++number_of_rebuilds_;*/ + } else { + ++num_steps_; + //next_slack_= next_slack_*.98; + } + } + + bool rebuilt=false; + if (!get_nbl_is_valid()) { + rebuild_nbl(); + rebuilt=true; + } else { + nbl_.erase(std::remove_if(nbl_.begin(), nbl_.end(), + internal::HasInactive()), + nbl_.end()); + } + return rebuilt; } void NonbondedListScoreState::show(std::ostream &out) const @@ -63,6 +151,15 @@ out << "NonbondedList" << std::endl; } +void NonbondedListScoreState::set_nbl_is_valid(bool tf) { + nbl_is_valid_= tf; + if (!nbl_is_valid_) { + NBL empty; + // free memory to make sure it shrinks + std::swap(empty, nbl_); + } +} + IMP_CONTAINER_IMPL(NonbondedListScoreState, BondedListScoreState, bonded_list, BondedListIndex, { if (0) std::cout <<*obj; Index: kernel/src/score_states/QuadraticBipartiteNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/QuadraticBipartiteNonbondedListScoreState.cpp (revision 541) +++ kernel/src/score_states/QuadraticBipartiteNonbondedListScoreState.cpp (working copy) @@ -1,66 +0,0 @@ -/** - * \file QuadraticBipartiteNonbondedListScoreState.cpp - * \brief Allow iteration through pairs of a set of atoms. - * - * Copyright 2007-8 Sali Lab. All rights reserved. - */ - -#include "IMP/score_states/QuadraticBipartiteNonbondedListScoreState.h" -#include "IMP/decorators/XYZDecorator.h" -#include "IMP/score_states/MaxChangeScoreState.h" - -namespace IMP -{ - -QuadraticBipartiteNonbondedListScoreState -::QuadraticBipartiteNonbondedListScoreState(FloatKey rk, - const Particles &ps0, - const Particles &ps1): P(rk) -{ - set_particles(ps0, ps1); -} - - -QuadraticBipartiteNonbondedListScoreState -::QuadraticBipartiteNonbondedListScoreState(FloatKey rk): P(rk) -{ -} - -void QuadraticBipartiteNonbondedListScoreState::rebuild_nbl(float cut) -{ - IMP_LOG(TERSE, "Rebuilding QBNBL with lists of size " << na_ - << " and " << P::get_particles().size() -na_ - << " and cutoff " << cut << std::endl); - for (unsigned int i=0; i< na_; ++i) { - for (unsigned int j=na_; j < P::get_particles().size(); ++j) { - P::handle_nbl_pair(P::get_particles()[i], - P::get_particles()[j], - cut); - } - } - IMP_LOG(TERSE, "NBL has " << P::nbl_.size() << " pairs" << std::endl); -} - -void QuadraticBipartiteNonbondedListScoreState -::set_particles(const Particles &ps0, - const Particles &ps1) -{ - IMP_LOG(TERSE, "Setting QBNBLSS particles " << ps0.size() - << " and " << ps1.size() << std::endl); - Particles all(ps0); - all.insert(all.end(), ps1.begin(), ps1.end()); - P::set_particles(all); - na_= ps0.size(); -} - -void QuadraticBipartiteNonbondedListScoreState::update() -{ - P::update(); -} - -void QuadraticBipartiteNonbondedListScoreState::show(std::ostream &out) const -{ - out << "QuadraticBipartiteNonbondedList" << std::endl; -} - -} // namespace IMP Index: kernel/src/score_states/CoverBondsScoreState.cpp =================================================================== --- kernel/src/score_states/CoverBondsScoreState.cpp (revision 0) +++ kernel/src/score_states/CoverBondsScoreState.cpp (revision 0) @@ -0,0 +1,78 @@ +#include "IMP/score_states/CoverBondsScoreState.h" +#include "IMP/score_states/BondDecoratorListScoreState.h" +#include "IMP/decorators/bond_decorators.h" +#include "IMP/decorators/XYZDecorator.h" + +namespace IMP +{ + + +CoverBondsScoreState::CoverBondsScoreState(BondDecoratorListScoreState *bl, + FloatKey rk): bl_(bl), rk_(rk) { +} + +CoverBondsScoreState::~CoverBondsScoreState() { +} + +void CoverBondsScoreState::do_before_evaluate() { + for (BondDecoratorListScoreState::BondIterator it= bl_->bonds_begin(); + it != bl_->bonds_end(); ++it) { + BondDecorator bd= *it; + BondedDecorator pa= bd.get_bonded(0); + BondedDecorator pb= bd.get_bonded(1); + IMP_LOG(VERBOSE, "Processing bond between " + << pa.get_particle()->get_index() + << " and " << pb.get_particle()->get_index() << std::endl); + XYZDecorator da(pa.get_particle()); + XYZDecorator db(pb.get_particle()); + IMP_LOG(VERBOSE, "Endpoints are " << da << " and " + << db << std::endl); + Vector3D diff= da.get_vector_to(db); + float len= diff.get_magnitude(); + Vector3D cv= da.get_vector() + .5*diff; + + XYZDecorator dxyz; + // slightly evil + try { + dxyz= XYZDecorator::cast(bd.get_particle()); + } catch (const InvalidStateException &ve) { + dxyz= XYZDecorator::create(bd.get_particle()); + } + dxyz.set_coordinates(cv); + + if (dxyz.get_particle()->has_attribute(rk_)) { + dxyz.get_particle()->set_value(rk_, len/2.0); + } else { + dxyz.get_particle()->add_attribute(rk_, len/2.0, false); + } + IMP_LOG(VERBOSE, "Bond cover is " << cv << ": " << len/2.0 << std::endl); + } +} + +void CoverBondsScoreState::after_evaluate(DerivativeAccumulator *dva) { + if (dva) { + for (BondDecoratorListScoreState::BondIterator it= bl_->bonds_begin(); + it != bl_->bonds_end(); ++it) { + XYZDecorator d(it->get_particle()); + Vector3D deriv; + // divide derivatives equally between endpoints + for (int i = 0; i < 3; ++i) { + deriv[i] = d.get_coordinate_derivative(i) / 2.0; + } + + for (int i=0; i< 2; ++i) { + BondedDecorator bd= it->get_bonded(i); + XYZDecorator d(bd.get_particle()); + for (int i = 0; i < 3; ++i) { + d.add_to_coordinate_derivative(i, deriv[i], *dva); + } + } + } + } +} + +void CoverBondsScoreState::show(std::ostream &out) const { + out << "CoverBondsScoreState on " << bl_->number_of_bonds() + << " bonds " << std::endl; +} +} Index: kernel/src/score_states/GravityCenterScoreState.cpp =================================================================== --- kernel/src/score_states/GravityCenterScoreState.cpp (revision 541) +++ kernel/src/score_states/GravityCenterScoreState.cpp (working copy) @@ -21,11 +21,19 @@ add_particles(ps); } + +GravityCenterScoreState::~GravityCenterScoreState(){} + + +void GravityCenterScoreState::do_before_evaluate() { + update_position(); +} + // check that the particle is an xyz particle IMP_LIST_IMPL(GravityCenterScoreState, Particle, particle, Particle*, - XYZDecorator::cast(obj), set_position()); + XYZDecorator::cast(obj), update_position()); -void GravityCenterScoreState::set_position() +void GravityCenterScoreState::update_position() { Vector3D cvect(0.0, 0.0, 0.0); bool do_weighting = (weightkey_ != FloatKey()); @@ -78,4 +86,10 @@ } + +void GravityCenterScoreState::show(std::ostream &out) const { + out << "GravityCenter" < + +namespace IMP +{ + + + +// Turn the default into an actual algorithm and work around missing algorithms +/* This cannot be shared with AllNBL because that one has grid and this does + not. + */ +static BipartiteNonbondedListScoreState::Algorithm +translate_algorithm(BipartiteNonbondedListScoreState::Algorithm a){ +#ifdef IMP_USE_CGAL + switch (a) { + case BipartiteNonbondedListScoreState::DEFAULT: + return BipartiteNonbondedListScoreState::BBOX; + default: + return a; + } +#else + switch (a) { + case BipartiteNonbondedListScoreState::BBOX: + IMP_WARN("BipartiteNonbondedListScoreState::BBOX requires CGAL support. " + << "GRID used instead." << std::endl); + case BipartiteNonbondedListScoreState::DEFAULT: + return BipartiteNonbondedListScoreState::QUADRATIC; + default: + return a; + } +#endif +} + + + +BipartiteNonbondedListScoreState +::BipartiteNonbondedListScoreState(Float cut, + FloatKey rk, + Algorithm a): + P(cut, rk), a_(translate_algorithm(a)) +{ + FloatKeys ks; + ks.push_back(rk); + mc0_= new MaxChangeScoreState(XYZDecorator::get_xyz_keys()); + mc1_= new MaxChangeScoreState(XYZDecorator::get_xyz_keys()); + mcr_= new MaxChangeScoreState(ks); +} + +BipartiteNonbondedListScoreState +::BipartiteNonbondedListScoreState(Float cut, FloatKey rk, + const Particles &ps0, + const Particles &ps1, + Algorithm a): + P(cut, rk), a_(translate_algorithm(a)) +{ + FloatKeys ks; + ks.push_back(rk); + mc0_= new MaxChangeScoreState(XYZDecorator::get_xyz_keys()); + mc1_= new MaxChangeScoreState(XYZDecorator::get_xyz_keys()); + mcr_= new MaxChangeScoreState(ks); + set_particles(ps0, ps1); +} + + +BipartiteNonbondedListScoreState::~BipartiteNonbondedListScoreState(){} + +void BipartiteNonbondedListScoreState::do_before_evaluate() { + mc0_->before_evaluate(ScoreState::get_before_evaluate_iteration()); + mc0_->before_evaluate(ScoreState::get_before_evaluate_iteration()); + mcr_->before_evaluate(ScoreState::get_before_evaluate_iteration()); + Float mc= std::max(mc0_->get_max(), mc1_->get_max()); + + Float cost; + switch (a_){ + case QUADRATIC: + cost= 10*mc0_->number_of_particles()* mc1_->number_of_particles(); + break; + case BBOX: + cost= 1000*(mcr_->number_of_particles()); + break; + default: + IMP_assert(0, "Bad algorithm"); + cost= 1000*(mcr_->number_of_particles()); + } + + if (P::update(mc+ mcr_->get_max(), cost)) { + mc0_->reset(); + mc1_->reset(); + mcr_->reset(); + } + IMP_BEGIN_CHECK(EXPENSIVE); + check_nbl(); + IMP_END_CHECK; +} + +void BipartiteNonbondedListScoreState::process_sets(const Particles &p0, + const Particles &p1) { + switch (a_) { + case QUADRATIC: + for (unsigned int j=0; j< p0.size(); ++j) { + for (unsigned int i=0; i< p1.size(); ++i) { + P::add_if_box_overlap(p1[i], p0[j]); + } + } + break; + case BBOX: + internal::bipartite_bbox_scan(p0, p1, P::get_radius_key(), + P::get_slack(), P::get_cutoff(), + internal::NBLAddPairIfNonbonded(this)); + break; + default: + IMP_failure("Can't find algorithm in BipartiteNonbondedListScoreState", + ErrorException()); + } +} + +void BipartiteNonbondedListScoreState::rebuild_nbl() { + IMP_LOG(TERSE, "Rebuilding BNBL with cutoff " + << P::get_cutoff() << " and slack " << P::get_slack() << std::endl); + process_sets(mc0_->get_particles(), mc1_->get_particles()); + P::set_nbl_is_valid(true); + IMP_LOG(TERSE, "NBL has " << P::number_of_nonbonded() + << " pairs" << std::endl); +} + +void BipartiteNonbondedListScoreState::set_algorithm(Algorithm a) { + a_=translate_algorithm(a); +} + + +void BipartiteNonbondedListScoreState::set_particles(const Particles &ps0, + const Particles &ps1) { + mc0_->clear_particles(); + mc0_->add_particles(ps0); + mc1_->clear_particles(); + mc1_->add_particles(ps1); + mcr_->clear_particles(); + mcr_->add_particles(P::particles_with_radius(ps0)); + mcr_->add_particles(P::particles_with_radius(ps1)); + P::set_nbl_is_valid(false); +} + + +void BipartiteNonbondedListScoreState::add_particles_0(const Particles &ps) { + if (P::get_nbl_is_valid()) process_sets(ps, mc1_->get_particles()); + mc0_->add_particles(ps); + mcr_->add_particles(P::particles_with_radius(ps)); + P::set_nbl_is_valid(false); +} + +void BipartiteNonbondedListScoreState::add_particles_1(const Particles &ps) { + if (P::get_nbl_is_valid()) process_sets(ps, mc0_->get_particles()); + mc1_->add_particles(ps); + mcr_->add_particles(P::particles_with_radius(ps)); + P::set_nbl_is_valid(false); +} + +void BipartiteNonbondedListScoreState::clear_particles() { + mc0_->clear_particles(); + mc1_->clear_particles(); + mcr_->clear_particles(); + P::set_nbl_is_valid(false); + P::set_nbl_is_valid(true); +} + + +void BipartiteNonbondedListScoreState::show(std::ostream &out) const +{ + out << "BipartiteNonbondedListScoreState" << std::endl; +} + +void BipartiteNonbondedListScoreState::check_nbl() const +{ + const Particles &ps0= mc0_->get_particles(); + const Particles &ps1= mc1_->get_particles(); + GetRadius gr= P::get_radius_object(); + for (unsigned int i=0; i< ps0.size(); ++i) { + XYZDecorator di= XYZDecorator::cast(ps0[i]); + for (unsigned int j=0; j< ps1.size(); ++j) { + XYZDecorator dj= XYZDecorator::cast(ps1[j]); + IMP_assert(ps0[i] != ps1[j], "Duplicate particles in list"); + if (distance(di, dj) - gr(ps0[i]) - gr(ps1[j]) + <= P::get_cutoff() && !are_bonded(ps0[i], ps1[j])) { + bool found=false; + for (NonbondedIterator nit= P::nonbonded_begin(); + nit != P::nonbonded_end(); ++nit) { + if (nit->first == ps0[i] && nit->second == ps1[j] + || nit->first == ps1[j] && nit->second == ps0[i]) { + IMP_assert(!found, "Entry is in list twice"); + found=true; + } + } + IMP_assert(found, "Nonbonded list is missing " + << ps0[i]->get_index() << " " << di + << " " << gr(ps0[i]) + << " and " << ps1[j]->get_index() << " " + << dj << gr(ps1[j]) + << " size is " << number_of_nonbonded() << std::endl); + } + } + } +} + + +unsigned int BipartiteNonbondedListScoreState::get_memory_usage() const { + return NonbondedListScoreState::get_memory_usage() + + sizeof(a_) + mc0_->get_memory_usage() + mc1_->get_memory_usage() + + mcr_->get_memory_usage(); +} + +} // namespace IMP Index: kernel/src/score_states/SConscript =================================================================== --- kernel/src/score_states/SConscript (revision 541) +++ kernel/src/score_states/SConscript (working copy) @@ -1,12 +1,10 @@ Import('env') - -files = ['NonbondedListScoreState.cpp', - 'MaxChangeScoreState.cpp', 'BondDecoratorListScoreState.cpp', - 'AllNonbondedListScoreState.cpp', - 'QuadraticBipartiteNonbondedListScoreState.cpp', - 'QuadraticAllNonbondedListScoreState.cpp', - 'QuadraticNonbondedListScoreState.cpp', 'GravityCenterScoreState.cpp', - 'ChainBondedListScoreState.cpp'] - -files = [File(x) for x in files] +files=[] +files.append(File( 'AllNonbondedListScoreState.cpp' )) +files.append(File( 'BipartiteNonbondedListScoreState.cpp' )) +files.append(File( 'BondDecoratorListScoreState.cpp' )) +files.append(File( 'CoverBondsScoreState.cpp' )) +files.append(File( 'GravityCenterScoreState.cpp' )) +files.append(File( 'MaxChangeScoreState.cpp' )) +files.append(File( 'NonbondedListScoreState.cpp' )) Return('files') Index: kernel/src/score_states/QuadraticAllNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/QuadraticAllNonbondedListScoreState.cpp (revision 541) +++ kernel/src/score_states/QuadraticAllNonbondedListScoreState.cpp (working copy) @@ -1,58 +0,0 @@ -/** - * \file QuadraticAllNonbondedListScoreState.cpp - * \brief Allow iteration through pairs of a set of s. - * - * Copyright 2007-8 Sali Lab. All rights reserved. - */ - -#include "IMP/score_states/QuadraticAllNonbondedListScoreState.h" -#include "IMP/decorators/XYZDecorator.h" - -#include - -namespace IMP -{ - - -QuadraticAllNonbondedListScoreState -::QuadraticAllNonbondedListScoreState(FloatKey radius, - const Particles &ps): P(radius) -{ - set_particles(ps); -} - -QuadraticAllNonbondedListScoreState::~QuadraticAllNonbondedListScoreState() -{ -} - -void QuadraticAllNonbondedListScoreState::set_particles(const Particles &ps) -{ - P::set_particles(ps); -} - - -void QuadraticAllNonbondedListScoreState::update() -{ - P::update(); -} - - -void QuadraticAllNonbondedListScoreState::rebuild_nbl(Float cut) -{ - IMP_LOG(TERSE, "Rebuilding QNBL with cutoff " << cut << std::endl); - const Particles &moving= P::get_particles(); - for (unsigned int j=0; j< moving.size(); ++j) { - for (unsigned int i=0; i< j; ++i) { - P::handle_nbl_pair(moving[i], moving[j], cut); - } - } - IMP_LOG(TERSE, "NBL has " << P::nbl_.size() << " pairs" << std::endl); -} - - -void QuadraticAllNonbondedListScoreState::show(std::ostream &out) const -{ - out << "QuadraticAllNonbondedListScoreState" << std::endl; -} - -} // namespace IMP Index: kernel/src/score_states/ChainBondedListScoreState.cpp =================================================================== --- kernel/src/score_states/ChainBondedListScoreState.cpp (revision 541) +++ kernel/src/score_states/ChainBondedListScoreState.cpp (working copy) @@ -1,64 +0,0 @@ -/** - * \file ChainBondedListScoreState.cpp - * \brief Allow iteration through pairs of a set of atoms. - * - * Copyright 2007-8 Sali Lab. All rights reserved. - */ - -#include "IMP/score_states/ChainBondedListScoreState.h" - -#include -#include - -namespace IMP -{ - -ChainBondedListScoreState::ChainBondedListScoreState() -{ - std::ostringstream oss; - oss << "ChainBLSS " << this; - cik_ = IntKey(oss.str().c_str()); - next_index_ = 0; -} - -void ChainBondedListScoreState::update() -{ - IMP_LOG(TERSE, "Updating ChainBondedList assumed to be static" << std::endl); -} - -void ChainBondedListScoreState::clear_chains() -{ - for (unsigned int i = 0; i < chains_.size(); ++i) { - for (unsigned int j = 0; j < chains_[i].size(); ++j) { - chains_[i][j]->set_value(cik_, -1); - } - } - chains_.clear(); - next_index_=0; -} - -void ChainBondedListScoreState::add_chain(const Particles &ps) -{ - chains_.push_back(internal::Vector(ps.begin(), ps.end())); - for (unsigned int i = 0; i < chains_.back().size(); ++i) { - Particle *p = chains_.back()[i]; - if (p->has_attribute(cik_)) { - p->set_value(cik_, next_index_); - } else { - p->add_attribute(cik_, next_index_); - } - ++next_index_; - } - ++next_index_; -} - - -bool ChainBondedListScoreState::are_bonded(Particle *a, Particle *b) const -{ - if (!a->has_attribute(cik_) || !b->has_attribute(cik_)) return false; - int ia= a->get_value(cik_); - int ib= b->get_value(cik_); - return std::abs(ia-ib) == 1; -} - -} // namespace IMP Index: kernel/src/score_states/MaxChangeScoreState.cpp =================================================================== --- kernel/src/score_states/MaxChangeScoreState.cpp (revision 541) +++ kernel/src/score_states/MaxChangeScoreState.cpp (working copy) @@ -43,7 +43,7 @@ } }, {reset();}); -void MaxChangeScoreState::update() +void MaxChangeScoreState::do_before_evaluate() { max_change_=0; // get rid of inactive particles and their stored values @@ -74,4 +74,8 @@ max_change_=0; } +void MaxChangeScoreState::show(std::ostream &out) const { + out << "MaxChangeScoreState" << std::endl; +} + } // namespace IMP Index: kernel/src/score_states/AllNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/AllNonbondedListScoreState.cpp (revision 541) +++ kernel/src/score_states/AllNonbondedListScoreState.cpp (working copy) @@ -7,101 +7,187 @@ #include "IMP/score_states/AllNonbondedListScoreState.h" #include "IMP/decorators/XYZDecorator.h" +#include "IMP/internal/utility.h" +#include "IMP/internal/bbox_nbl_helpers.h" +#include "IMP/score_states/MaxChangeScoreState.h" +#include "IMP/internal/ParticleGrid.h" -#include +#include namespace IMP { -static unsigned int min_grid_size=20; -AllNonbondedListScoreState -::AllNonbondedListScoreState(FloatKey radius, const Particles &ps): - P(radius) -{ - set_particles(ps); +// Turn the default into an actual algorithm and work around missing algorithms +static AllNonbondedListScoreState::Algorithm +translate_algorithm(AllNonbondedListScoreState::Algorithm a){ +#ifdef IMP_USE_CGAL + switch (a) { + case AllNonbondedListScoreState::DEFAULT: + return AllNonbondedListScoreState::BBOX; + default: + return a; + } +#else + switch (a) { + case AllNonbondedListScoreState::BBOX: + IMP_WARN("AllNonbondedListScoreState::BBOX requires CGAL support. " + << "GRID used instead." << std::endl); + case AllNonbondedListScoreState::DEFAULT: + return AllNonbondedListScoreState::GRID; + default: + return a; + } +#endif } -AllNonbondedListScoreState::~AllNonbondedListScoreState() +AllNonbondedListScoreState::AllNonbondedListScoreState(Float cut, + FloatKey rk, + const Particles &ps, + Algorithm a): + P(cut, rk), a_(translate_algorithm(a)) { - cleanup(bins_); + FloatKeys ks; + ks.push_back(rk); + mc_= new MaxChangeScoreState(XYZDecorator::get_xyz_keys()); + mcr_= new MaxChangeScoreState(ks); + add_particles(ps); } -float AllNonbondedListScoreState::side_from_r(float r) const +AllNonbondedListScoreState::AllNonbondedListScoreState(Float cut, + FloatKey rk, + Algorithm a): + P(cut, rk), a_(translate_algorithm(a)) { - if (r==0) return 1; - else return r*1.6; + FloatKeys ks; + ks.push_back(rk); + mc_= new MaxChangeScoreState(XYZDecorator::get_xyz_keys()); + mcr_= new MaxChangeScoreState(ks); } -Particles AllNonbondedListScoreState::get_particles() const -{ - Particles ret; - for (unsigned int i=0; i< bins_.size(); ++i) { - ret.insert(ret.end(), bins_[i].grid->get_particles().begin(), - bins_[i].grid->get_particles().end()); - } - return ret; +AllNonbondedListScoreState::~AllNonbondedListScoreState(){ } -void AllNonbondedListScoreState::set_particles(const Particles &ps) -{ - NonbondedListScoreState::clear_nbl(); - cleanup(bins_); - repartition_points(ps, bins_); +void AllNonbondedListScoreState::do_before_evaluate() { + mc_->before_evaluate(ScoreState::get_before_evaluate_iteration()); + mcr_->before_evaluate(ScoreState::get_before_evaluate_iteration()); + Float mc=mc_->get_max()+ mcr_->get_max(); + Float cost; + switch (a_){ + case QUADRATIC: + cost= 10*square(mc_->number_of_particles()); + break; + case BBOX: + cost= 1000*(mc_->number_of_particles()); + case GRID: + // completely made up + cost*=2; + break; + default: + IMP_failure("Bad algorithm", ErrorException()); + cost= 10*(mc_->number_of_particles()); + } + if (P::update(mc, cost)) { + mc_->reset(); + } + IMP_BEGIN_CHECK(EXPENSIVE); + check_nbl(); + IMP_END_CHECK; } -void AllNonbondedListScoreState::add_particles(const Particles &ps) -{ - if (bins_.empty()) set_particles(ps); - else { -#ifndef NDEBUG - for (unsigned int i=0; i< ps.size(); ++i) { - for (unsigned int j=0; j< bins_.size(); ++j) { - for (unsigned int k=0; k< bins_[j].grid->get_particles().size(); ++k) { - IMP_assert(ps[i] != bins_[j].grid->get_particles()[k], - "Can't add a particle that is already there"); - IMP_assert(ps[i]->get_index() - != bins_[j].grid->get_particles()[k]->get_index(), - "Same particle index, different particles."); - } +void AllNonbondedListScoreState::rebuild_nbl() { + IMP_LOG(TERSE, "Rebuilding AllNBL with cutoff " + << P::get_cutoff() << " and slack " << P::get_slack() << std::endl); + if (a_== QUADRATIC) { + const Particles &moving= mc_->get_particles(); + for (unsigned int j=0; j< moving.size(); ++j) { + for (unsigned int i=0; i< j; ++i) { + P::add_if_box_overlap(moving[i], moving[j]); } } -#endif + } else if (a_ == GRID) { + grid_rebuild_nbl(); + } else if (a_== BBOX) { + internal::bbox_scan(mc_->get_particles(), P::get_radius_key(), + P::get_slack(), P::get_cutoff(), + internal::NBLAddPairIfNonbonded(this)); - for (unsigned int i=0; i< ps.size(); ++i) { - float r= P::get_radius(ps[i]); - bool found=false; - for (unsigned int j=0; j< bins_.size(); ++j) { - if (bins_[j].rmax >= r) { - found=true; - IMP_LOG(VERBOSE, "Adding particle " - << ps[i]->get_index() << " to bin " << j << std::endl); - bins_[j].grid->add_particle(ps[i]); - break; + } else { + IMP_failure("Bad algorithm in AllNBL::rebuild", ErrorException()); + } + set_nbl_is_valid(true); + IMP_LOG(TERSE, "NBL has " << P::number_of_nonbonded() + << " pairs" << std::endl); +} + +void AllNonbondedListScoreState::set_particles(const Particles &ps) { + mc_->clear_particles(); + mc_->add_particles(ps); + mcr_->clear_particles(); + mcr_->add_particles(particles_with_radius(ps)); + P::set_nbl_is_valid(false); +} + + +void AllNonbondedListScoreState::add_particles(const Particles &ps) { + if (P::get_nbl_is_valid()) { + if (a_== QUADRATIC || a_ == GRID) { + const Particles &moving= mc_->get_particles(); + for (unsigned int j=0; j< moving.size(); ++j) { + for (unsigned int i=0; i< ps.size(); ++i) { + P::add_if_box_overlap(ps[i], moving[j]); } } - if (!found) { - bins_.back().rmax=r; - bins_.back().grid->add_particle(ps[i]); - } + } else if (a_== BBOX) { + internal::bipartite_bbox_scan(mc_->get_particles(), ps, + P::get_radius_key(), + P::get_slack(), P::get_cutoff(), + internal::NBLAddPairIfNonbonded(this)); } } - IMP_LOG(TERSE, "Destroying nbl in list due to additions"<< std::endl); - NonbondedListScoreState::clear_nbl(); + mcr_->add_particles(particles_with_radius(ps)); + mc_->add_particles(ps); } +void AllNonbondedListScoreState::clear_particles() { + mc_->clear_particles(); + mcr_->clear_particles(); + P::set_nbl_is_valid(false); + P::set_nbl_is_valid(true); +} -void AllNonbondedListScoreState::repartition_points(const Particles &ps, - std::vector &out) +void AllNonbondedListScoreState::show(std::ostream &out) const { - cleanup(out); - if (ps.empty()) return; + out << "AllNonbondedListScoreState" << std::endl; +} + + +void AllNonbondedListScoreState::set_algorithm(Algorithm a) { + a_=translate_algorithm(a); +} + + +// methods for grid + +static unsigned int min_grid_size=20; + +float AllNonbondedListScoreState::grid_side_from_r(float r) const +{ + if (r==0) return 1; + else return r*1.6; +} + +void AllNonbondedListScoreState +::grid_partition_points( internal::Vector &bins) +{ + if (mc_->get_particles().empty()) return; + GetRadius gr= P::get_radius_object(); float minr=std::numeric_limits::max(), maxr=0; - for (unsigned int i=0; i< ps.size(); ++i) { - float r= P::get_radius(ps[i]); + for (unsigned int i=0; i< mc_->get_particles().size(); ++i) { + float r= gr(mc_->get_particles()[i]); if ( r > maxr) maxr=r; if ( r > 0 && r < minr) minr=r; } @@ -116,13 +202,13 @@ cuts.push_back(2*maxr); std::vector ops(cuts.size()); - for (unsigned int i=0; i< ps.size(); ++i) { - float r= P::get_radius(ps[i]); + for (unsigned int i=0; i< mc_->get_particles().size(); ++i) { + float r= gr(mc_->get_particles()[i]); bool found=false; for (unsigned int j=0; ; ++j) { IMP_assert(j< cuts.size(), "Internal error in ASNBLSS"); if (cuts[j] >= r) { - ops[j].push_back(ps[i]); + ops[j].push_back(mc_->get_particles()[i]); found=true; break; } @@ -138,156 +224,122 @@ } for (unsigned int i=0; i< cuts.size(); ++i) { if (ops[i].empty()) continue; - out.push_back(Bin()); float rmax=0; for (unsigned int j=0; j< ops[i].size(); ++j) { - rmax= std::max(rmax, P::get_radius(ops[i][j])); + rmax= std::max(rmax, gr(ops[i][j])); } - out.back().rmax= rmax; - internal::ParticleGrid *pg - = new internal::ParticleGrid(side_from_r(out.back().rmax)); - out.back().grid= pg; - out.back().grid->add_particles(ops[i]); + bins.push_back(new internal::ParticleGrid(grid_side_from_r(rmax), 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 << 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 AllNonbondedListScoreState::cleanup(std::vector &bins) -{ + IMP_LOG(VERBOSE, "Created " << bins.size() << " grids" << std::endl); for (unsigned int i=0; i< bins.size(); ++i) { - delete bins[i].grid; + IMP_LOG(VERBOSE, *bins[i] << std::endl); } - bins.clear(); } -void AllNonbondedListScoreState::update() -{ - IMP_LOG(TERSE, "Updating nonbonded list"<< std::endl); - NonbondedListScoreState::update(); - bool bad=false; - for (unsigned int i=0; i< bins_.size(); ++i) { - bad = bins_[i].grid->update() || bad; - IMP_LOG(VERBOSE, bins_[i].rmax << std::endl << *bins_[i].grid << std::endl); - } - if (bad) { - IMP_LOG(TERSE, "Destroying nbl in list"<< std::endl); - NonbondedListScoreState::clear_nbl(); - } -} -void AllNonbondedListScoreState::generate_nbl(const Bin &particle_bin, - const Bin &grid_bin, - float cut) +void AllNonbondedListScoreState +::grid_generate_nbl(const internal::ParticleGrid *particle_bin, + const internal::ParticleGrid *grid_bin) { - 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); + IMP_CHECK_OBJECT(particle_bin); + IMP_CHECK_OBJECT(grid_bin); + IMP_LOG(VERBOSE, "Generate nbl for pair " << std::endl); + for (internal::ParticleGrid::ParticleVoxelIterator + it= particle_bin->particle_voxels_begin(); + it != particle_bin->particle_voxels_end(); ++it) { + Particle *p= it->first; + internal::NBLAddIfNonbonded 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())); + = grid_bin->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); + << " from " << index << std::endl); + grid_bin->apply_to_nearby(f, index, + P::get_cutoff() + 2*P::get_slack(), + false); } } -void AllNonbondedListScoreState::rebuild_nbl(Float cut) +void AllNonbondedListScoreState::grid_rebuild_nbl() { - IMP_LOG(TERSE, "Rebuilding NBL with " << bins_.size() << " bins" - << " and cutoff " << cut << 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); + IMP_LOG(TERSE, "Rebuilding NBL with Grid and cutoff " + << P::get_cutoff() << std::endl ); + internal::Vector bins; + grid_partition_points(bins); + + for (unsigned int i=0; i< bins.size(); ++i) { + for (unsigned int j=i+1; j< bins.size(); ++j) { + grid_generate_nbl(bins[i], bins[j]); } internal::ParticleGrid::Index last_index; for (internal::ParticleGrid::ParticleVoxelIterator it - = bins_[i].grid->particle_voxels_begin(); - it != bins_[i].grid->particle_voxels_end(); ++it) { + = bins[i]->particle_voxels_begin(); + it != bins[i]->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, + internal::NBLAddIfNonbonded f(this, it->first); + bins[i]->apply_to_nearby(f, it->second, + P::get_cutoff()+2*P::get_slack(), true); + if (it->second != last_index) { + internal::NBLAddPairIfNonbonded fp(this); IMP_LOG(VERBOSE, "Searching in " << it->second << std::endl); - bins_[i].grid->apply_to_cell_pairs(f, it->second); + bins[i]->apply_to_cell_pairs(fp, it->second); } last_index=it->second; } } +} - IMP_LOG(TERSE, "NBL has " << P::nbl_.size() << " pairs" << std::endl); -#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]); - IMP_assert(ps[i] != ps[j], "Duplicate particles in grid"); - if (distance(di, dj) - P::get_radius(ps[i]) - P::get_radius(ps[j]) - <= 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; - } +// debugging +void AllNonbondedListScoreState::check_nbl() const +{ + if (!get_nbl_is_valid()) return; + const Particles &ps= mc_->get_particles(); + P::GetRadius gr= P::get_radius_object(); + 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]); + IMP_assert(ps[i] != ps[j], "Duplicate particles in list"); + if (distance(di, dj) - gr(ps[i]) - gr(ps[j]) + <= P::get_cutoff() && !are_bonded(ps[i], ps[j])) { + bool found=false; + for (NonbondedIterator nit= P::nonbonded_begin(); + nit != P::nonbonded_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); } + IMP_assert(found, "Nonbonded list is missing " + << ps[i]->get_index() << " " << di + << " " << gr(ps[i]) + << " and " << ps[j]->get_index() << " " + << dj << gr(ps[j]) + << " size is " << number_of_nonbonded() << std::endl); } } } -#endif } -void AllNonbondedListScoreState::show(std::ostream &out) const -{ - out << "AllNonbondedListScoreState" << std::endl; + +unsigned int AllNonbondedListScoreState::get_memory_usage() const { + /*std::cout << "MU: " << mc_->get_memory_usage() << " " + << mcr_->get_memory_usage() << " " + << NonbondedListScoreState::get_memory_usage() << std::endl;*/ + return NonbondedListScoreState::get_memory_usage() + + sizeof(a_) + mc_->get_memory_usage() + + mcr_->get_memory_usage(); } } // namespace IMP Index: kernel/src/score_states/QuadraticNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/QuadraticNonbondedListScoreState.cpp (revision 541) +++ kernel/src/score_states/QuadraticNonbondedListScoreState.cpp (working copy) @@ -1,60 +0,0 @@ -/** - * \file QuadraticNonbondedListScoreState.cpp - * \brief Allow iteration through pairs of a set of s. - * - * Copyright 2007-8 Sali Lab. All rights reserved. - */ - -#include "IMP/score_states/QuadraticNonbondedListScoreState.h" -#include "IMP/decorators/XYZDecorator.h" - -#include - -namespace IMP -{ - -QuadraticNonbondedListScoreState -::QuadraticNonbondedListScoreState(FloatKey radius): - P(radius), slack_(.1) -{ - mc_= new MaxChangeScoreState(XYZDecorator::get_xyz_keys()); -} - -QuadraticNonbondedListScoreState::~QuadraticNonbondedListScoreState() -{ -} - - -void QuadraticNonbondedListScoreState::update() -{ - // placeholder to do tuning of slack - NonbondedListScoreState::update(); - if (mc_->get_max() > slack_) { - NonbondedListScoreState::clear_nbl(); - mc_->reset(); - } -} - -void QuadraticNonbondedListScoreState -::handle_nbl_pair( Particle *a, Particle *b, - float d) -{ - XYZDecorator da= XYZDecorator::cast(a); - XYZDecorator db= XYZDecorator::cast(b); - float ra= P::get_radius(a); - float rb= P::get_radius(b); - for (unsigned int i=0; i< 3; ++i) { - float delta=std::abs(da.get_coordinate(i) - db.get_coordinate(i)); - if (delta -ra -rb > d-slack_) { - IMP_LOG(VERBOSE, "Pair " << a->get_index() - << " and " << b->get_index() << " rejected on coordinate " - << i << std::endl); - return ; - } - } - IMP_LOG(VERBOSE, "Adding pair " << a->get_index() - << " and " << b->get_index() << std::endl); - P::add_to_nbl(a, b); -} - -} // namespace IMP Index: kernel/src/restraints/NonbondedRestraint.cpp =================================================================== --- kernel/src/restraints/NonbondedRestraint.cpp (revision 541) +++ kernel/src/restraints/NonbondedRestraint.cpp (working copy) @@ -17,9 +17,8 @@ { NonbondedRestraint::NonbondedRestraint(PairScore *ps, - NonbondedListScoreState *nbl, - Float md) : nbl_(nbl), sf_(ps), - max_dist_(md) + NonbondedListScoreState *nbl) + : nbl_(nbl), sf_(ps) { } @@ -31,12 +30,12 @@ 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_)) + << std::distance(nbl_->nonbonded_begin(), + nbl_->nonbonded_end()) << " pairs" << std::endl); for (NonbondedListScoreState::NonbondedIterator it - = nbl_->nonbonded_begin(max_dist_); - it != nbl_->nonbonded_end(max_dist_); ++it) { + = nbl_->nonbonded_begin(); + it != nbl_->nonbonded_end(); ++it) { float thisscore = sf_->evaluate(it->first, it->second, accum); if (thisscore != 0) { IMP_LOG(VERBOSE, "Pair " << it->first->get_index() @@ -57,4 +56,9 @@ out << std::endl; } +unsigned int NonbondedRestraint::get_memory_usage() const { + return sizeof(nbl_) + sizeof(sf_)+ Restraint::get_memory_usage() + + sf_->get_memory_usage(); +} + } // namespace IMP Index: kernel/src/restraints/AngleRestraint.cpp =================================================================== --- kernel/src/restraints/AngleRestraint.cpp (revision 541) +++ kernel/src/restraints/AngleRestraint.cpp (working copy) @@ -55,4 +55,9 @@ out << std::endl; } +unsigned int AngleRestraint::get_memory_usage() const { + return Restraint::get_memory_usage() + + sizeof(sf_) + sf_->get_memory_usage(); +} + } // namespace IMP Index: kernel/src/restraints/SConscript =================================================================== --- kernel/src/restraints/SConscript (revision 541) +++ kernel/src/restraints/SConscript (working copy) @@ -1,10 +1,14 @@ Import('env') - -files = ['ConnectivityRestraint.cpp', 'RestraintSet.cpp', - 'DistanceRestraint.cpp', 'AngleRestraint.cpp', 'DihedralRestraint.cpp', - 'NonbondedRestraint.cpp', 'BondDecoratorRestraint.cpp', - 'SingletonListRestraint.cpp', 'PairListRestraint.cpp', - 'TripletChainRestraint.cpp', 'PairChainRestraint.cpp'] - -files = [File(x) for x in files] +files=[] +files.append(File( 'AngleRestraint.cpp' )) +files.append(File( 'BondDecoratorRestraint.cpp' )) +files.append(File( 'ConnectivityRestraint.cpp' )) +files.append(File( 'DihedralRestraint.cpp' )) +files.append(File( 'DistanceRestraint.cpp' )) +files.append(File( 'NonbondedRestraint.cpp' )) +files.append(File( 'PairChainRestraint.cpp' )) +files.append(File( 'PairListRestraint.cpp' )) +files.append(File( 'RestraintSet.cpp' )) +files.append(File( 'SingletonListRestraint.cpp' )) +files.append(File( 'TripletChainRestraint.cpp' )) Return('files') Index: kernel/src/triplet_scores/SConscript =================================================================== --- kernel/src/triplet_scores/SConscript (revision 541) +++ kernel/src/triplet_scores/SConscript (working copy) @@ -1,6 +1,4 @@ Import('env') - -files = ['AngleTripletScore.cpp'] - -files = [File(x) for x in files] +files=[] +files.append(File( 'AngleTripletScore.cpp' )) Return('files') Index: kernel/src/optimizers/states/SConscript =================================================================== --- kernel/src/optimizers/states/SConscript (revision 541) +++ kernel/src/optimizers/states/SConscript (working copy) @@ -1,7 +1,6 @@ Import('env') - -files = ['VRMLLogOptimizerState.cpp', 'CMMLogOptimizerState.cpp', - 'VelocityScalingOptimizerState.cpp'] - -files = [File(x) for x in files] +files=[] +files.append(File( 'CMMLogOptimizerState.cpp' )) +files.append(File( 'VelocityScalingOptimizerState.cpp' )) +files.append(File( 'VRMLLogOptimizerState.cpp' )) Return('files') Index: kernel/src/optimizers/states/VRMLLogOptimizerState.cpp =================================================================== --- kernel/src/optimizers/states/VRMLLogOptimizerState.cpp (revision 541) +++ kernel/src/optimizers/states/VRMLLogOptimizerState.cpp (working copy) @@ -17,9 +17,10 @@ VRMLLogOptimizerState::VRMLLogOptimizerState(std::string filename, const Particles &pis) : - pis_(pis), filename_(filename), file_number_(0), call_number_(0), + filename_(filename), file_number_(0), call_number_(0), skip_steps_(0) { + set_particles(pis); } void VRMLLogOptimizerState::update() @@ -44,13 +45,15 @@ if (!out) { IMP_WARN("Error opening VRML log file " << buf); } else { - IMP_LOG(VERBOSE, "Writing " << pis_.size() + IMP_LOG(VERBOSE, "Writing " << number_of_particles() << " particles to file " << buf << "..." << std::flush); - write(pis_, radius_, r_, g_, b_, out); + write(out); //IMP_LOG(TERSE, "done" << std::endl); } } +IMP_LIST_IMPL(VRMLLogOptimizerState, Particle, particle, Particle*, ,); + static Float snap(Float f) { if (f < 0) return 0; @@ -58,32 +61,37 @@ return f; } -void VRMLLogOptimizerState::write(const Particles &pis, FloatKey rk, - FloatKey r, FloatKey g, FloatKey b, - std::ostream &out) +void VRMLLogOptimizerState::set_color(int c, Float r, Float g, Float b) { + colors_[c]= boost::make_tuple(snap(r),snap(g),snap(b)); +} + +void VRMLLogOptimizerState::write(std::ostream &out) const { out << "#VRML V2.0 utf8\n"; out << "Group {\n"; out << "children [\n"; - for (unsigned int i = 0; i < pis.size(); ++i) { + for (ParticleConstIterator pit=particles_begin(); + pit != particles_end(); ++pit) { + Particle *p = *pit; try { - Particle *p = pis[i]; - XYZDecorator xyz = XYZDecorator::cast(p); + XYZDecorator xyz = XYZDecorator::cast(p); float x = xyz.get_x(); float y = xyz.get_y(); float z = xyz.get_z(); Float rv = -1, gv = -1, bv = -1; - if (r != FloatKey() && b != FloatKey() && g != FloatKey() - && p->has_attribute(r) && p->has_attribute(g) - && p->has_attribute(b)) { - rv = snap(p->get_value(r)); - gv = snap(p->get_value(g)); - bv = snap(p->get_value(b)); + if (color_ != IntKey() + && p->has_attribute(color_)) { + int cv = p->get_value(color_); + if (colors_.find(cv) != colors_.end()) { + rv= colors_.find(cv)->second.get<0>(); + gv= colors_.find(cv)->second.get<1>(); + bv= colors_.find(cv)->second.get<2>(); + } } Float radius = .1; - if (rk != FloatKey() && p->has_attribute(rk)) { - radius = p->get_value(rk); + if (radius_ != FloatKey() && p->has_attribute(radius_)) { + radius = p->get_value(radius_); //oss << ".sphere " << x << " " << y << " " << z << " " << r << "\n"; } @@ -111,7 +119,7 @@ out << "}\n"; } catch (InvalidStateException &e) { - IMP_WARN("Particle " << pis[i] << " does not have " + IMP_WARN("Particle " << p << " does not have " << " cartesian coordinates"); } } Index: kernel/src/optimizers/SConscript =================================================================== --- kernel/src/optimizers/SConscript (revision 541) +++ kernel/src/optimizers/SConscript (working copy) @@ -1,10 +1,11 @@ Import('env') - -states_files = SConscript('states/SConscript') -movers_files = SConscript('movers/SConscript') - -files = [ 'SteepestDescent.cpp', 'ConjugateGradients.cpp', - 'MolecularDynamics.cpp', 'MonteCarlo.cpp', 'MoverBase.cpp']+states_files +movers_files - -files = [File(x) for x in files] +files=[] +files.append(File( 'BrownianDynamics.cpp' )) +files.append(File( 'ConjugateGradients.cpp' )) +files.append(File( 'MolecularDynamics.cpp' )) +files.append(File( 'MonteCarlo.cpp' )) +files.append(File( 'MoverBase.cpp' )) +files.append(File( 'SteepestDescent.cpp' )) +files.append(SConscript( 'movers/SConscript' )) +files.append(SConscript( 'states/SConscript' )) Return('files') Index: kernel/src/optimizers/BrownianDynamics.cpp =================================================================== --- kernel/src/optimizers/BrownianDynamics.cpp (revision 0) +++ kernel/src/optimizers/BrownianDynamics.cpp (revision 0) @@ -0,0 +1,302 @@ +/** + * \file BrownianDynamics.cpp \brief Simple brownian dynamics optimizer. + * + * Copyright 2007-8 Sali Lab. All rights reserved. + * + */ + +#include "IMP/log.h" +#include "IMP/random.h" +#include "IMP/internal/constants.h" +#include "IMP/internal/units.h" +#include "IMP/optimizers/BrownianDynamics.h" +#include "IMP/decorators/XYZDecorator.h" +#include + +#include + +namespace IMP +{ +typedef internal::MKSUnit<-3, -1, 1, 0, -1> MilliPascalSeconds; + +static MilliPascalSeconds eta(internal::Kelvin T) { + // from wikipedia + const std::pair points[] + ={ std::make_pair(internal::Kelvin(273+10.0), + MilliPascalSeconds(1.308)), + std::make_pair(internal::Kelvin(273+20.0), + MilliPascalSeconds(1.003)), + std::make_pair(internal::Kelvin(273+30.0), + MilliPascalSeconds(0.7978)), + std::make_pair(internal::Kelvin(273+40.0), + MilliPascalSeconds(0.6531)), + std::make_pair(internal::Kelvin(273+50.0), + MilliPascalSeconds(0.5471)), + std::make_pair(internal::Kelvin(273+60.0), + MilliPascalSeconds(0.4668)), + std::make_pair(internal::Kelvin(273+70.0), + MilliPascalSeconds(0.4044)), + std::make_pair(internal::Kelvin(273+80.0), + MilliPascalSeconds(0.3550)), + std::make_pair(internal::Kelvin(273+90.0), + MilliPascalSeconds(0.3150)), + std::make_pair(internal::Kelvin(273+100.0), + MilliPascalSeconds(0.2822))}; + + const unsigned int npoints= sizeof(points)/sizeof(std::pair); + if (T < points[0].first) return points[0].second; + else { + for (unsigned int i=1; i< npoints; ++i) { + if (points[i].first > T) { + internal::Scalar f= (T - points[i-1].first)/(points[i].first + - points[i-1].first); + MilliPascalSeconds ret= + (internal::Scalar(1.0)-f) *points[i-1].second + f*points[i].second; + return ret; + } + } + } + return points[npoints-1].second; +} + +BrownianDynamics::BrownianDynamics(FloatKey dkey) : + max_change_(10), max_dt_(1e7), cur_dt_(max_dt_), cur_time_(0), + T_(300.0), dkey_(dkey) +{ + IMP_check(dkey_ != FloatKey(), "BrownianDynamics needs a valid key for the " + << "diffusion coefficient", + ValueException("Bad diffusion coef key")); +} + + +BrownianDynamics::~BrownianDynamics() +{ +} + +IMP_LIST_IMPL(BrownianDynamics, Particle, particle, Particle*, + {if (0) std::cout << obj << index;},); + + +void BrownianDynamics::set_time_step(internal::FemtoSecond t) { + time_steps_.clear(); + max_dt_ = t; + cur_dt_= max_dt_; +} + +void BrownianDynamics::setup_particles() +{ + clear_particles(); + FloatKeys xyzk=XYZDecorator::get_xyz_keys(); + for (unsigned int i = 0; i < get_model()->number_of_particles(); ++i) { + Particle *p = get_model()->get_particle(i); + if (p->has_attribute(xyzk[0]) && p->get_is_optimized(xyzk[0]) + && p->has_attribute(xyzk[1]) && p->get_is_optimized(xyzk[1]) + && p->has_attribute(xyzk[2]) && p->get_is_optimized(xyzk[2]) + && p->has_attribute(dkey_) && !p->get_is_optimized(dkey_)) { + add_particle(p); + } + } +} + +// rt_dt= rt+Fdt/psi + sqrt(2kTdt/psi)omega +// omega has variance 6 phi kT + + +/* + radius + if step is xi*radius^2/(3pi E) for some E=1 + then the motion per step should be sqrt(2kT/(pi E)) + + T* is + */ + + +bool BrownianDynamics::step() { + std::vector new_pos(number_of_particles()); + bool noshrink=true; + while (!propose_step(new_pos)) { + cur_dt_= cur_dt_/internal::Scalar(2.0); + noshrink=false; + } + + for (unsigned int i=0; i< number_of_particles(); ++i) { + XYZDecorator d(get_particle(i)); + for (unsigned int j=0; j< 3; ++j) { + d.set_coordinate(j, new_pos[i][j]); + } + } + return noshrink; +} + + +//! Perform a single dynamics step. +bool BrownianDynamics::propose_step(std::vector& new_pos) +{ + // Assuming score is in kcal/mol, its derivatives in kcal/mol/angstrom, + // and mass is in g/mol, conversion factor necessary to get accelerations + // in angstrom/fs/fs from raw derivatives + // we want + //const double eta= 10;// kg/(m s) from ~1 cg/(cm s) + //const double pi=M_PI; + //const double kTnu= 1.3806504*T_.get_value(); + FloatKeys xyzk=XYZDecorator::get_xyz_keys(); + + IMP_LOG(VERBOSE, "dt is " << cur_dt_ << std::endl); + + + for (unsigned int i=0; i< number_of_particles(); ++i) { + Particle *p= get_particle(i); + XYZDecorator d(p); + //double xi= 6*pi*eta*radius; // kg/s + internal::Cm2PerSecond D(p->get_value(dkey_)); + internal::Angstrom sigma(compute_sigma_from_D(D)); +#ifndef NDEBUG + internal::Angstrom osigma(sqrt(internal::Scalar(2.0)*D*cur_dt_)); +#endif + IMP_assert(sigma + - osigma + <= internal::Scalar(.01)* sigma, + "Sigma computations don't match " << sigma + << " " + << sqrt(internal::Scalar(2.0)*D*cur_dt_)); + IMP_LOG(VERBOSE, p->get_index() << ": sigma is " + << sigma << std::endl); + boost::normal_distribution mrng(0, sigma.get_value()); + boost::variate_generator > + sampler(random_number_generator, mrng); + + //std::cout << p->get_index() << std::endl; + + internal::Angstrom delta[3]; + + for (unsigned j = 0; j < 3; ++j) { + // force originally in kcal/mol/A + // derive* 2.390e-4 gives it in J/(mol A) + // that * e10 gives it in J/(mol m) + // that / NA gives it in J/m + internal::KCalPerAMol force(-d.get_coordinate_derivative(j)); + //*4.1868e+13/NA; old conversion + internal::FemtoNewton nforce= convert_to_mks(force); + internal::Angstrom R(sampler()); + internal::Angstrom force_term=nforce*cur_dt_*D/kt(); + //std::cout << "Force term is " << force_term << " and R is " + //<< R << std::endl; + internal::Angstrom dr= force_term + R; //std::sqrt(2*kb*T_*ldt/xi) * + // get back from meters + delta[j]=dr; + } + //internal::Angstrom max_motion= internal::Scalar(4)*sigma; + /*std::cout << "delta is " << delta << " mag is " + << delta.get_magnitude() << " sigma " << sigma << std::endl;*/ + + IMP_LOG(VERBOSE, "For particle " << p->get_index() + << " delta is " << delta[0] << " " << delta[1] << " " << delta[2] + << " from a force of " + << "[" << d.get_coordinate_derivative(0) + << ", " << d.get_coordinate_derivative(1) + << ", " << d.get_coordinate_derivative(2) << "]" << std::endl); + + internal::SquaredAngstrom t= square(delta[0]); + + internal::SquaredAngstrom magnitude2=square(delta[0])+square(delta[1]) + +square(delta[2]); + + //internal::SquaredAngstrom max_motion2= square(max_motion); + if (magnitude2 > square(max_change_)) { + return false; + } + for (unsigned int j=0; j< 3; ++j) { + new_pos[i][j]= d.get_coordinate(j) + delta[j].get_value(); + } + } + return true; +} + + +//! Optimize the model. +/** \param[in] max_steps Maximum number of iterations before aborting. + \return score of the final state of the model. + */ +Float BrownianDynamics::optimize(unsigned int max_steps) +{ + setup_particles(); + IMP_LOG(SILENT, "Running brownian dynamics on " << get_particles().size() + << " particles with a step of " << cur_dt_ << std::endl); + unsigned int num_const_dt=0; + for (unsigned int i = 0; i < max_steps; ++i) { + update_states(); + get_model()->evaluate(true); + if (step()) { + ++num_const_dt; + } + + cur_time_= cur_time_+cur_dt_; + if (num_const_dt == 10) { + num_const_dt=0; + cur_dt_= min(cur_dt_*internal::Scalar(2), max_dt_); + } + { + internal::FemtoSecond lb(max_dt_/internal::Scalar(2.0)); + std::vector::size_type d=0; + while (lb > cur_dt_) { + lb = lb/internal::Scalar(2.0); + ++d; + } + time_steps_.resize(std::max(time_steps_.size(), d+1), 0); + ++time_steps_[d]; + } + } + return get_model()->evaluate(false); +} + + + +internal::FemtoJoule BrownianDynamics::kt() const { + return internal::KB*T_; +} + + +internal::Angstrom +BrownianDynamics::estimate_radius_from_mass_units(Float mass_in_kd) const { + /* Data points used: + thyroglobulin 669kD 85A + catalase 232kD 52A + aldolase 158kD 48A + */ + + internal::Angstrom r(std::pow(mass_in_kd, .3333f)*10); + return r; +} + +internal::Cm2PerSecond +BrownianDynamics::estimate_D_from_radius(internal::Angstrom r) const { + + MilliPascalSeconds e=eta(T_); + internal::MKSUnit<-13, 0, 1, 0, -1> etar( e*r); + /*std::cout << e << " " << etar << " " << kt << std::endl; + std::cout << "scalar etar " << (internal::Scalar(6*internal::PI)*etar) + << std::endl; + std::cout << "ret pre conv " << (kt/(internal::Scalar(6* internal::PI)*etar)) + << std::endl;*/ + internal::Cm2PerSecond ret( kt()/(internal::Scalar(6* internal::PI)*etar)); + //std::cout << "ret " << ret << std::endl; + return ret; +} + +internal::Angstrom +BrownianDynamics::compute_sigma_from_D(internal::Cm2PerSecond D) const { + return sqrt(internal::Scalar(2.0)*D*cur_dt_); //6*xi*kb*T_; +} + +internal::KCalPerAMol BrownianDynamics +::get_force_scale_from_D(internal::Cm2PerSecond D) const { + // force motion is f*cur_dt_*D/kT + // sigma*kT/(cur_dt_*D) + return convert_to_kcal(internal::Angstrom(compute_sigma_from_D(D))*kt() + /(max_dt_*D)); +} + + + +} // namespace IMP Index: kernel/src/optimizers/movers/SConscript =================================================================== --- kernel/src/optimizers/movers/SConscript (revision 541) +++ kernel/src/optimizers/movers/SConscript (working copy) @@ -1,6 +1,5 @@ Import('env') - -files = [ 'BallMover.cpp', 'NormalMover.cpp'] - -files = [File(x) for x in files] +files=[] +files.append(File( 'BallMover.cpp' )) +files.append(File( 'NormalMover.cpp' )) Return('files') Index: kernel/src/decorators/SConscript =================================================================== --- kernel/src/decorators/SConscript (revision 541) +++ kernel/src/decorators/SConscript (working copy) @@ -1,8 +1,10 @@ Import('env') - -files = ['AtomDecorator.cpp', 'bond_decorators.cpp', - 'HierarchyDecorator.cpp', 'NameDecorator.cpp', 'ResidueDecorator.cpp', - 'MolecularHierarchyDecorator.cpp', 'XYZDecorator.cpp'] - -files = [File(x) for x in files] +files=[] +files.append(File( 'AtomDecorator.cpp' )) +files.append(File( 'bond_decorators.cpp' )) +files.append(File( 'HierarchyDecorator.cpp' )) +files.append(File( 'MolecularHierarchyDecorator.cpp' )) +files.append(File( 'NameDecorator.cpp' )) +files.append(File( 'ResidueDecorator.cpp' )) +files.append(File( 'XYZDecorator.cpp' )) Return('files') Index: kernel/src/decorators/XYZDecorator.cpp =================================================================== --- kernel/src/decorators/XYZDecorator.cpp (revision 541) +++ kernel/src/decorators/XYZDecorator.cpp (working copy) @@ -55,6 +55,7 @@ } } + IMP_DECORATOR_INITIALIZE(XYZDecorator, DecoratorBase, { key_[0] = FloatKey("x"); @@ -74,4 +75,12 @@ return std::sqrt(d2); } +Float sphere_distance(XYZDecorator a, XYZDecorator b, FloatKey rk) +{ + double d2= d(a.get_x(), b.get_x()) + d(a.get_y(), b.get_y()) + + d(a.get_z(), b.get_z()); + double d= std::sqrt(d2); + return d - a.get_particle()->get_value(rk) - b.get_particle()->get_value(rk); +} + } // namespace IMP Index: kernel/src/pair_scores/DistancePairScore.cpp =================================================================== --- kernel/src/pair_scores/DistancePairScore.cpp (revision 541) +++ kernel/src/pair_scores/DistancePairScore.cpp (working copy) @@ -29,8 +29,13 @@ Float d2 = 0, delta[3]; Float score; - XYZDecorator d0 = XYZDecorator::cast(a); - XYZDecorator d1 = XYZDecorator::cast(b); + IMP_BEGIN_CHECK(CHEAP); + XYZDecorator::cast(a); + XYZDecorator::cast(b); + IMP_END_CHECK; + + XYZDecorator d0(a); + XYZDecorator d1(b); for (int i = 0; i < 3; ++i) { delta[i] = d0.get_coordinate(i) - d1.get_coordinate(i); d2 += square(delta[i]); Index: kernel/src/pair_scores/SConscript =================================================================== --- kernel/src/pair_scores/SConscript (revision 541) +++ kernel/src/pair_scores/SConscript (working copy) @@ -1,6 +1,5 @@ Import('env') - -files = ['DistancePairScore.cpp', 'SphereDistancePairScore.cpp'] - -files = [File(x) for x in files] +files=[] +files.append(File( 'DistancePairScore.cpp' )) +files.append(File( 'SphereDistancePairScore.cpp' )) Return('files') Index: kernel/src/pair_scores/SphereDistancePairScore.cpp =================================================================== --- kernel/src/pair_scores/SphereDistancePairScore.cpp (revision 541) +++ kernel/src/pair_scores/SphereDistancePairScore.cpp (working copy) @@ -21,6 +21,12 @@ Float SphereDistancePairScore::evaluate(Particle *a, Particle *b, DerivativeAccumulator *da) { + IMP_check(a->has_attribute(radius_), "Particle " << a->get_index() + << "missing radius in SphereDistancePairScore", + ValueException("Missing radius")); + IMP_check(b->has_attribute(radius_), "Particle " << b->get_index() + << "missing radius in SphereDistancePairScore", + ValueException("Missing radius")); Float ra = a->get_value(radius_); Float rb = b->get_value(radius_); return internal::evaluate_distance_pair_score(a,b, da, f_.get(), Index: kernel/src/internal/SConscript =================================================================== --- kernel/src/internal/SConscript (revision 541) +++ kernel/src/internal/SConscript (working copy) @@ -1,7 +1,9 @@ Import('env') - -files = ['graph_base.cpp', 'ParticleGrid.cpp', 'Object.cpp', - 'kernel_version_info.cpp', 'constants.cpp'] - -files = [File(x) for x in files] +files=[] +files.append(File( 'bbox_nbl_helpers.cpp' )) +files.append(File( 'constants.cpp' )) +files.append(File( 'graph_base.cpp' )) +files.append(File( 'kernel_version_info.cpp' )) +files.append(File( 'Object.cpp' )) +files.append(File( 'ParticleGrid.cpp' )) Return('files') Index: kernel/src/internal/bbox_nbl_helpers.cpp =================================================================== --- kernel/src/internal/bbox_nbl_helpers.cpp (revision 0) +++ kernel/src/internal/bbox_nbl_helpers.cpp (revision 0) @@ -0,0 +1,92 @@ +/** + * \file bbox_nbl_helpers.cpp + * \brief Helpers for the cgal-based NBL. + * + * Copyright 2007-8 Sali Lab. All rights reserved. + */ + +#include "IMP/internal/bbox_nbl_helpers.h" +#include "IMP/decorators/XYZDecorator.h" +#include "IMP/score_states/NonbondedListScoreState.h" + +#ifdef IMP_USE_CGAL +#include +#include +#endif + +namespace IMP +{ + +namespace internal +{ + +struct NBLBbox { + XYZDecorator d_; + typedef Float NT; + typedef void * ID; + Float r_; + NBLBbox(){} + NBLBbox(Particle *p, + Float r): d_(p), + r_(r){} + static unsigned int dimension() {return 3;} + void *id() const {return d_.get_particle();} + NT min_coord(unsigned int i) const { + return d_.get_coordinate(i)-r_; + } + NT max_coord(unsigned int i) const { + return d_.get_coordinate(i)+r_; + } + // make it so I can reused the callback provide by NBLSS + operator Particle*() {return d_.get_particle();} +}; + +#ifdef IMP_USE_CGAL +static void copy_particles_to_boxes(const Particles &ps, + FloatKey rk, Float slack, Float cutoff, + std::vector &boxes){ + boxes.resize(ps.size()); + for (unsigned int i=0; i< ps.size(); ++i) { + Particle *p= ps[i]; + + Float r= slack+cutoff/2.0; + if (rk != FloatKey() && p->has_attribute(rk)) { + r+= p->get_value(rk); + } + boxes[i]=internal::NBLBbox(p, r); + } +} +#endif + +void bipartite_bbox_scan(const Particles &ps0, const Particles &ps1, + FloatKey rk, Float slack, Float cutoff, + const NBLAddPairIfNonbonded &ap) { +#ifdef IMP_USE_CGAL + std::vector boxes0, boxes1; + copy_particles_to_boxes(ps0, rk, slack, cutoff, boxes0); + copy_particles_to_boxes(ps1, rk, slack, cutoff, boxes1); + + CGAL::box_intersection_d( boxes0.begin(), boxes0.end(), + boxes1.begin(), boxes1.end(), ap); +#else + IMP_failure( "IMP built without CGAL support.", ErrorException()); +#endif +} + +void bbox_scan(const Particles &ps, + FloatKey rk, Float slack, Float cutoff, + const NBLAddPairIfNonbonded &ap) { +#ifdef IMP_USE_CGAL + std::vector boxes; + copy_particles_to_boxes(ps, rk, slack, cutoff, boxes); + + + CGAL::box_self_intersection_d( boxes.begin(), boxes.end(), ap); +#else + IMP_failure("IMP built without CGAL support.", ErrorException()); +#endif +} + + +} +} Index: kernel/src/internal/ParticleGrid.cpp =================================================================== --- kernel/src/internal/ParticleGrid.cpp (revision 541) +++ kernel/src/internal/ParticleGrid.cpp (working copy) @@ -18,19 +18,20 @@ static const int target_cell_occupancy=10; -ParticleGrid::ParticleGrid(float tvs): target_voxel_side_(tvs), - grid_valid_(false) + ParticleGrid::ParticleGrid(float tvs, + const Particles &ps): target_voxel_side_(tvs) { IMP_assert(tvs >0, "Target voxel edge size must be positive"); - mc_= new MaxChangeScoreState(XYZDecorator::get_xyz_keys()); + build_grid(ps); } -ParticleGrid::ParticleGrid(): target_voxel_side_(0), grid_valid_(0) +ParticleGrid::ParticleGrid(): target_voxel_side_(0) { } -void ParticleGrid::build_grid() +void ParticleGrid::build_grid(const Particles &ps) { + audit_particles(ps); IMP_LOG(TERSE, "Creating nonbonded grid..." << std::flush); float mn[3]= {std::numeric_limits::max(), std::numeric_limits::max(), @@ -38,19 +39,19 @@ float mx[3]={-std::numeric_limits::max(), -std::numeric_limits::max(), -std::numeric_limits::max()}; - for (unsigned int i = 0; i < mc_->get_particles().size(); ++i) { - XYZDecorator d= XYZDecorator::cast(mc_->get_particles()[i]); + for (unsigned int i = 0; i < ps.size(); ++i) { + XYZDecorator d(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); } } - if (!mc_->get_particles().empty()) { + if (!ps.empty()) { // keep the grid size sane if things blow up float maxdim= std::max(mx[0]-mn[0], std::max(mx[1]-mn[1], mx[2]-mn[2])); float vx= std::pow(static_cast(target_cell_occupancy *(maxdim*maxdim*maxdim - /mc_->get_particles().size())), + /ps.size())), .3333f); if (vx > target_voxel_side_) { IMP_LOG(VERBOSE, "Overroade target side of " << target_voxel_side_ @@ -62,89 +63,16 @@ Vector3D(mn[0], mn[1], mn[2]), Vector3D(mx[0], mx[1], mx[2]), Particles()); - for (unsigned int i = 0; i < mc_->get_particles().size(); ++i) { - XYZDecorator d= XYZDecorator::cast(mc_->get_particles()[i]); + for (unsigned int i = 0; i < ps.size(); ++i) { + XYZDecorator d(ps[i]); Vector3D v(d.get_x(), d.get_y(), d.get_z()); - grid_.get_voxel(grid_.get_index(v)).push_back(mc_->get_particles()[i]); + grid_.get_voxel(grid_.get_index(v)).push_back(ps[i]); } - grid_valid_=true; - mc_->reset(); IMP_LOG(TERSE, "done." << std::endl); } -void ParticleGrid::add_particles(const Particles &ps) -{ - audit_particles(ps); - mc_->add_particles(ps); - for (unsigned int i=0; i< ps.size(); ++i) { - if (grid_valid_) { - add_particle_to_grid(ps[i]); - } - } -} -void ParticleGrid::add_particle(Particle *p) -{ - Particles ps(1, p); - audit_particles(ps); - mc_->add_particle(p); - - if (grid_valid_) { - add_particle_to_grid(p); - } -} - -void ParticleGrid::add_particle_to_grid(Particle *p) -{ - IMP_assert(grid_valid_, "Bad call of add particle to grid"); - XYZDecorator d= XYZDecorator::cast(p); - Vector3D v(d.get_x(), d.get_y(), d.get_z()); - Grid::VirtualIndex vi= grid_.get_virtual_index(v); - Grid::Index gi= grid_.get_index(vi); - if (gi== Grid::Index()) { - IMP_LOG(TERSE, "Adding particle off grid invalidates it " - << v << " " << vi << std::endl); - grid_valid_=false; - grid_ = Grid(); - } else { - grid_.get_voxel(gi).push_back(p); - } -} - -void ParticleGrid::clear_particles() -{ - mc_->clear_particles(); -} - -bool ParticleGrid::update() -{ - bool ret; - if (!grid_valid_ || mc_->get_max() > target_voxel_side_) { - IMP_LOG(TERSE, "Rebuilding particle grid\n"); - build_grid(); - ret= true; - } else { - IMP_LOG(TERSE, "Removing inactive particles\n"); - for (Grid::DataIterator dit= grid_.data_begin(); dit != grid_.data_end(); - ++dit) { - remove_inactive_particles(*dit); - } - ret= false; - } - unsigned int ssz=0; - for (Grid::DataIterator dit= grid_.data_begin(); dit != grid_.data_end(); - ++dit) { - ssz+= dit->size(); - } - // do this last since it has the ref counts - mc_->update(); - - IMP_assert(ssz== mc_->number_of_particles(), "Particle mismatch in PG"); - - return ret; -} - void ParticleGrid::audit_particles(const Particles &ps) const { for (unsigned int i=0; i< ps.size(); ++i) { Index: kernel/pyext/IMP/test.py =================================================================== --- kernel/pyext/IMP/test.py (revision 541) +++ kernel/pyext/IMP/test.py (working copy) @@ -6,12 +6,11 @@ class TestCase(unittest.TestCase): """Super class for IMP test cases""" - def assertInTolerance(self, num1, num2, tolerance, msg=None): + def assertInTolerance(self, num1, num2, tolerance, msg=""): """Assert that the difference between num1 and num2 is less than tolerance""" diff = abs(num1 - num2) - if msg is None: - msg = "%f != %f within %g" % (num1, num2, tolerance) + msg = msg+str("\n%f != %f within %g" % (num1, num2, tolerance)) self.assert_(diff < tolerance, msg) def create_point_particle(self, model, x, y, z): @@ -41,3 +40,32 @@ dy = p1.get_value(ykey) - p2.get_value(ykey) dz = p1.get_value(zkey) - p2.get_value(zkey) return math.sqrt(dx*dx + dy*dy + dz*dz) + + def check_unary_function_deriv(self, uf, lb, ub, step): + for i in range(0, int((ub-lb)/step)): + f= lb+ i*step + (v,d)= uf.evaluate_deriv(f) + if (i != 0): + offset= step/1024 + vmn= uf.evaluate(f-offset) + vmx= uf.evaluate(f+offset) + da= (vmx-vmn)/(2*offset) + print str(f) + ": " + str(d) + ", " + str(da) + if (abs(d-da) > abs(.1 *d)): + print f + print d + print da + self.assert_(abs(d-da) <= abs(.1 *d)+.0001, + "Deriv and approximation are not close") + def check_unary_function_min(self, uf, lb, ub, step, minf): + vmin = uf.evaluate(lb) + fmin = lb + for i in range(0, int((ub-lb)/step)): + f= lb+ i*step + (v,d)= uf.evaluate_deriv(f) + cv= uf.evaluate(f) + self.assertInTolerance(cv, v, .1*cv+.001) + if (v < vmin): + fmin= f + vmin= v + self.assert_(abs(fmin - minf) < step, "Wrong minimum") Index: kernel/pyext/SConscript =================================================================== --- kernel/pyext/SConscript (revision 541) +++ kernel/pyext/SConscript (working copy) @@ -4,7 +4,7 @@ # Get a modified build environment suitable for building Python extensions: e = get_pyext_environment(env, 'IMP', cplusplus=True) -e.Append(CPPPATH=['#/kernel/include', e['BOOST_CPPPATH']]) +e.Prepend(CPPPATH=['#/kernel/include', e['BOOST_CPPPATH']]) e.Append(LIBPATH=['../src'], LIBS='imp') # Build the Python extension from SWIG interface file: Index: kernel/pyext/IMP.i =================================================================== --- kernel/pyext/IMP.i (revision 541) +++ kernel/pyext/IMP.i (working copy) @@ -23,7 +23,7 @@ if (type(a)() == a): raise IndexError("Cannot use default Index") if (not p.has_attribute(a)): - raise IndexError("Particle does not have attribute") + raise IndexError("Particle does not have attribute "+str(a.get_string())) %} namespace IMP { @@ -136,6 +136,7 @@ %feature("director"); %include "IMP/Key.h" +%include "IMP/exception.h" %include "IMP/internal/Object.h" %include "IMP/internal/RefCountedObject.h" %include "IMP/Index.h" @@ -149,6 +150,7 @@ %include "IMP/unary_functions/ClosedCubicSpline.h" %include "IMP/unary_functions/Cosine.h" %include "IMP/unary_functions/Linear.h" +%include "IMP/unary_functions/WormLikeChain.h" %include "IMP/DerivativeAccumulator.h" %include "IMP/Restraint.h" %include "IMP/ScoreState.h" @@ -171,6 +173,7 @@ %include "IMP/Optimizer.h" %include "IMP/optimizers/SteepestDescent.h" %include "IMP/optimizers/ConjugateGradients.h" +%include "IMP/optimizers/BrownianDynamics.h" %include "IMP/optimizers/MolecularDynamics.h" %include "IMP/optimizers/Mover.h" %include "IMP/optimizers/MoverBase.h" @@ -197,15 +200,12 @@ %include "IMP/restraints/PairListRestraint.h" %include "IMP/restraints/RestraintSet.h" %include "IMP/score_states/BondedListScoreState.h" -%include "IMP/score_states/ChainBondedListScoreState.h" +%include "IMP/score_states/CoverBondsScoreState.h" %include "IMP/score_states/MaxChangeScoreState.h" %include "IMP/score_states/NonbondedListScoreState.h" %include "IMP/score_states/AllNonbondedListScoreState.h" +%include "IMP/score_states/BipartiteNonbondedListScoreState.h" %include "IMP/score_states/BondDecoratorListScoreState.h" -%include "IMP/score_states/QuadraticNonbondedListScoreState.h" -%include "IMP/score_states/QuadraticAllNonbondedListScoreState.h" -%include "IMP/score_states/QuadraticBipartiteNonbondedListScoreState.h" -%include "IMP/score_states/BipartiteNonbondedListScoreState.h" %include "IMP/score_states/GravityCenterScoreState.h" namespace IMP { @@ -217,8 +217,8 @@ %template(BondedListIndex) Index; %template(FloatKey) Key; %template(IntKey) Key; + %template(ParticleKey) Key; %template(StringKey) Key; - %template(ParticleKey) Key; %template(AtomType) Key; %template(ResidueType) Key; %template(show_named_hierarchy) show; @@ -237,4 +237,5 @@ %template(Floats) ::std::vector; %template(Strings) ::std::vector; %template(Ints) ::std::vector; -} + +} \ No newline at end of file