[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[IMP-dev] More patches with a couple API changes



These huge patch includes a few API changes, so please read the next section. It includes all changes that were not checked in yet from the last patch.

The patch includes automatically generated SConscript and IMP.h files and all the tests pass.

API changes:
- The constructors for the nonbonded lists no longer take the radius key as an argument. They default to XYZRDecorator::get_radius_key and the value can be set using the set_radius_key method

- UnaryFunctions and Scores now have const evaluate methods. Things have been updated accordingly. Note that errors in python about having the wrong number of arguments for the constructor probably are a result of pure virtual methods caused by the class not matching its parent any more.

- The arguments to the constructor of Linear were swapped to make it more similar to Harmonic.

- IMP_failure takes an exception type as the third argument (rather than an exception-typed expression)

- More changes to the names of the hierarchy and molecular hierarchy decorator functions to make them more consistent. Those should be stable now

- TunnelRestraint is now TunnelSingletonScore. Not sure why I made it a restraint before.

- As before there is now a ManualBondDecoratorListScoreState which is a base class for the old one and doesn't scan for new bonds.

- PairChainListRestraint only supports one chain as was discussed previously.

- Various methods were moved out of the IMP.Test class into IMP.test so they can be used outside of the test class itself (for example in test restraints)




API additions:
- Added an XYZRDecorator which has a radius and some functionality on spheres. The various things which take radius keys now default to the XYZR radius key

- yaml save and restore is added. Only the whole model can be restored at once to avoid breaking invariants.

- There is now a PairRestraint which acts as the base for DistanceRestraint.

- various missing add_x methods were filled in.

- I added a function in log.h to control whether exception messages are dumped to standard error in addition to being embedded in the exception. The current default is yes. The python lib init code should probably set it to false for python, but I don't know how to do that.

- Pointer is now implicitly convertible to a pointer. This removes a lot of get calls and makes code look nicer. And I can't think of any reason not to do this given how we use it and the fact that it is an intrusive pointer.

- vector3Ds now have inplace operators (+=, -=, /=)


Bug fixes and cleanups:
- There are many minor cleanups to which headers are included, and comments and things like that.

- typedef pair score test was removing all references to the model an expecting the particles to still be around, which is not safe

- IMP.i was setting thisown on the wrong argument to TypedPairScore::set_pair_score

- the derivatives are all cleared at once in the attribute table

- ObjectPointer had not been removed from SVN for some reason

- Grid3D now works with bool values

- MaxChangeScoreState now stores the base values internally in order to avoid changing the set of attributes on the particles. Changing the attributes broke the yaml reader.

- nbl_helpers now doesn't complain when NDEBUG is already defined.

- Things are cleaned up for when models are reference counted but they are not at the moment due to swig's lack of real reference counting



Index: kernel/test/singleton_scores/test_tunnel.py
===================================================================
--- kernel/test/singleton_scores/test_tunnel.py	(revision 0)
+++ kernel/test/singleton_scores/test_tunnel.py	(working copy)
@@ -13,17 +13,17 @@
         m= IMP.Model()
         p= IMP.Particle()
         m.add_particle(p)
-        d= IMP.XYZDecorator.create(p)
-        rk= IMP.FloatKey("radiusk")
-        p.add_attribute(rk, 1, False)
+        d= IMP.XYZRDecorator.create(p)
+        d.set_radius(1)
         f= IMP.HarmonicLowerBound(0, 1)
-        tr= IMP.TunnelRestraint(f, rk)
-        tr.set_height(10)
-        tr.set_radius(5)
-        tr.set_center(IMP.Vector3D(10,10,10))
-        tr.add_particle(p)
-        tr.set_coordinate(2)
-        m.add_restraint(tr)
+        tss= IMP.TunnelSingletonScore(f)
+        tss.set_height(10)
+        tss.set_radius(5)
+        tss.set_center(IMP.Vector3D(10,10,10))
+        tss.set_coordinate(2)
+        sl= IMP.SingletonListRestraint(tss)
+        sl.add_particle(p)
+        m.add_restraint(sl)
         m.evaluate(True)
         d.set_coordinates(IMP.Vector3D(10,10,10))
         self.assertEqual(m.evaluate(True), 0)

Property changes on: kernel/test/singleton_scores/test_tunnel.py
___________________________________________________________________
Added: svn:mergeinfo

Index: kernel/test/unary_functions/test_linear.py
===================================================================
--- kernel/test/unary_functions/test_linear.py	(revision 626)
+++ kernel/test/unary_functions/test_linear.py	(working copy)
@@ -9,7 +9,7 @@
         """Test that linear values are correct"""
         for offset in (0.0, -1.0):
             for slope in (0.0, -5.0, 3.5):
-                func = IMP.Linear(slope, offset)
+                func = IMP.Linear(offset, slope)
                 for i in range(15):
                     val = -10.0 + 3.5 * i
                     scoreonly = func.evaluate(val)
@@ -29,7 +29,7 @@
 
     def test_show(self):
         """Check Linear::show() method"""
-        func = IMP.Linear(1.0)
+        func = IMP.Linear(0, 1.0)
         func.show()
 
 if __name__ == '__main__':
Index: kernel/test/particles/test_refcount.py
===================================================================
--- kernel/test/particles/test_refcount.py	(revision 626)
+++ kernel/test/particles/test_refcount.py	(working copy)
@@ -12,6 +12,8 @@
         IMP.set_log_level(IMP.VERBOSE)
         self.basenum= IMP.RefCountedObject.get_number_of_live_objects()
         print "The base number of objects is " + str(self.basenum)
+        # set this to true if the model is ref counted
+        self.count_model=0
 
     def _check_number(self, expected):
         print "Expected "+str(expected)\
@@ -19,7 +21,7 @@
                             - self.basenum)
         self.assertEqual(IMP.RefCountedObject.get_number_of_live_objects() - self.basenum,
                          expected,
-                         "wrong number of particles")
+                         "wrong number of objects")
 
     def __test_simple(self):
         """Check that ref counting of particles works within python"""
@@ -47,11 +49,12 @@
         m= IMP.Model()
         p= IMP.Particle()
         pi= m.add_particle(p)
-        self._check_number(1)
+        self._check_number(1+self.count_model)
         m.remove_particle(pi)
-        self._check_number(1)
         self.assert_(not p.get_is_active(), "Removed particle is still active")
         p=1
+        self._check_number(0+ self.count_model)
+        m=1
         self._check_number(0)
 
     # This test does not work since swig refcounting is broken
@@ -123,17 +126,17 @@
         print "Add particle to mc"
         mc.add_particle(p)
         # also have the score state now
-        self._check_number(2)
+        self._check_number(2+ self.count_model)
         print "Remove from model"
         m.remove_particle(pi)
-        self._check_number(2)
+        self._check_number(2+ self.count_model)
         p=1
-        self._check_number(2)
+        self._check_number(2+ self.count_model)
         print "Remove from mc"
         mc.clear_particles()
-        self._check_number(1)
+        self._check_number(1+ self.count_model)
         mc=0
-        self._check_number(0)
+        self._check_number(0+ self.count_model)
 
     def test_skip(self):
         """Check that removed particles are skipped"""
@@ -148,6 +151,7 @@
         self.assertEqual(len(ps), 0, "Should no particles particle")
 
     def test_restraints(self):
+        """Check reference counting of restraints"""
         m= IMP.Model()
         r= IMP.ConstantRestraint(1)
         s= IMP.RestraintSet()
@@ -160,7 +164,27 @@
         m=0
         self._check_number(0)
 
+    def __test_model(self):
+        """Check reference counting of models"""
+        # disabled due to the usual issues with swig
+        self._check_number(0)
+        m= IMP.Model()
+        o= IMP.ConjugateGradients()
+        self._check_number(1)
+        o.set_model(m)
+        self._check_number(1)
+        o=[]
+        self._check_number(1)
+        m=[]
+        self._check_number(0)
+        m= IMP.Model()
+        o= IMP.ConjugateGradients()
+        self._check_number(1)
+        o.set_model(m)
+        self._check_number(1)
+        m=[]
+        self._check_number(1)
+        o=[]
 
-
 if __name__ == '__main__':
     unittest.main()
Index: kernel/test/distance/test_distance.py
===================================================================
--- kernel/test/distance/test_distance.py	(revision 626)
+++ kernel/test/distance/test_distance.py	(working copy)
@@ -37,6 +37,7 @@
                    IMP.HarmonicLowerBound(mean, 0.1),
                    IMP.Harmonic(mean, 0.1)):
             r = IMP.DistanceRestraint(sf, p1, p2)
+            r.set_was_owned(True)
             self.rsrs.append(r)
 
     def _make_restraints(self):
@@ -56,6 +57,7 @@
                    IMP.HarmonicLowerBound(3.0, 0.1), IMP.Harmonic(3.0, 0.1)):
             r = IMP.DistanceRestraint(fs, self.particles[1],
                                       self.particles[0])
+            r.set_was_owned(True)
             self.rsrs.append(r)
 
         # exceed lower bound
@@ -63,6 +65,7 @@
                    IMP.HarmonicLowerBound(5.0, 0.1), IMP.Harmonic(5.0, 0.1)):
             r = IMP.DistanceRestraint(fs, self.particles[1],
                                       self.particles[2])
+            r.set_was_owned(True)
             self.rsrs.append(r)
 
         # exceed upper bound
@@ -70,6 +73,7 @@
                    IMP.HarmonicLowerBound(4.0, 0.1), IMP.Harmonic(4.0, 0.1)):
             r = IMP.DistanceRestraint(fs, self.particles[0],
                                       self.particles[2])
+            r.set_was_owned(True)
             self.rsrs.append(r)
 
     def test_show(self):
@@ -77,6 +81,7 @@
         r = IMP.DistanceRestraint( IMP.Harmonic(0.0, 0.1),
                                    self.particles[1],
                                    self.particles[0])
+        r.set_was_owned(True)
         r.show()
 
     def test_distance(self):
Index: kernel/test/restraints/test_angles.py
===================================================================
--- kernel/test/restraints/test_angles.py	(revision 626)
+++ kernel/test/restraints/test_angles.py	(working copy)
@@ -31,7 +31,7 @@
             m.add_particle(p)
             d= IMP.XYZDecorator.create(p)
             l0.append(p)
-        self.randomize_particles(l0, 20)
+        IMP.test.randomize_particles(l0, 20)
         return l0
 
     def create_angle_r(self, s, ps):
Index: kernel/test/restraints/test_list.py
===================================================================
--- kernel/test/restraints/test_list.py	(revision 626)
+++ kernel/test/restraints/test_list.py	(working copy)
@@ -11,16 +11,6 @@
     def show(self, t):
         print "One Singleton"
 
-class Linear(IMP.UnaryFunction):
-    def __init__(self):
-        IMP.UnaryFunction.__init__(self)
-    def evaluate(self, feat):
-        return feat
-    def evaluate_deriv(self, feat):
-        return feat, 1.0
-    def show(self, *args):
-        print "identity"
-
 class TestList(IMP.test.TestCase):
     def setUp(self):
         IMP.test.TestCase.setUp(self)
@@ -49,7 +39,7 @@
         d.set_z(1)
         d.set_coordinates_are_optimized(True)
         v= IMP.Vector3D(3,1,5)
-        l= Linear()
+        l= IMP.Linear(0,1)
         s= IMP.DistanceToSingletonScore(l, v)
         r= IMP.SingletonListRestraint(s, m.get_particles())
         m.add_restraint(r)
Index: kernel/test/restraints/test_pairchain.py
===================================================================
--- kernel/test/restraints/test_pairchain.py	(revision 626)
+++ kernel/test/restraints/test_pairchain.py	(working copy)
@@ -26,20 +26,14 @@
             p= IMP.Particle()
             m.add_particle(p)
             ps0.append(p)
-        ps1= IMP.Particles()
-        for i in range(0,10):
-            p= IMP.Particle()
-            m.add_particle(p)
-            ps1.append(p)
         os= IndexDiff()
         s= IMP.PairChainRestraint(os)
-        s.add_chain(ps0)
-        s.add_chain(ps1)
+        s.set_particles(ps0)
         m.add_restraint(s)
         score= m.evaluate(False)
         print str(score)
-        self.assertEqual(score, 18, "Wrong score")
-        s.clear_chains()
+        self.assertEqual(score, 9, "Wrong score")
+        s.clear_particles()
         self.assertEqual(m.evaluate(False), 0, "Should be no terms")
 
 
Index: kernel/test/restraints/test_tunnel.py
===================================================================
--- kernel/test/restraints/test_tunnel.py	(revision 626)
+++ kernel/test/restraints/test_tunnel.py	(working copy)
@@ -1,62 +0,0 @@
-import unittest
-import IMP
-import IMP.test
-import IMP.utils
-import math
-
-class TunnelTest(IMP.test.TestCase):
-    """Tests for tunnel restraints"""
-
-    def test_score(self):
-        """Test derivatives and score of tunnel restraint"""
-        IMP.set_log_level(IMP.VERBOSE)
-        m= IMP.Model()
-        p= IMP.Particle()
-        m.add_particle(p)
-        d= IMP.XYZDecorator.create(p)
-        rk= IMP.FloatKey("radiusk")
-        p.add_attribute(rk, 1, False)
-        f= IMP.HarmonicLowerBound(0, 1)
-        tr= IMP.TunnelRestraint(f, rk)
-        tr.set_height(10)
-        tr.set_radius(5)
-        tr.set_center(IMP.Vector3D(10,10,10))
-        tr.add_particle(p)
-        tr.set_coordinate(2)
-        m.add_restraint(tr)
-        m.evaluate(True)
-        d.set_coordinates(IMP.Vector3D(10,10,10))
-        self.assertEqual(m.evaluate(True), 0)
-
-        print "Test left"
-        d.set_coordinates(IMP.Vector3D(4, 10, 10))
-        self.assert_(m.evaluate(True)>0)
-        self.assertEqual(d.get_coordinate_derivative(2), 0)
-        self.assertEqual(d.get_coordinate_derivative(1), 0)
-        print d.get_coordinate_derivative(0)
-        self.assert_(d.get_coordinate_derivative(0) > 0)
-
-        print "Test above"
-        d.set_coordinates(IMP.Vector3D(10, 4, 10))
-        self.assert_(m.evaluate(True)>0)
-        self.assertEqual(d.get_coordinate_derivative(2), 0)
-        self.assertEqual(d.get_coordinate_derivative(0), 0)
-        print d.get_coordinate_derivative(1)
-        self.assert_(d.get_coordinate_derivative(1) > 0)
-
-        print "Test bottom"
-        d.set_coordinates(IMP.Vector3D(30, 30, 3))
-        self.assert_(m.evaluate(True)>0)
-        self.assertEqual(d.get_coordinate_derivative(0), 0)
-        self.assertEqual(d.get_coordinate_derivative(1), 0)
-        self.assert_(d.get_coordinate_derivative(2) < 0)
-
-        print "Test top"
-        d.set_coordinates(IMP.Vector3D(30, 30, 17))
-        self.assert_(m.evaluate(True)>0)
-        self.assertEqual(d.get_coordinate_derivative(0), 0)
-        self.assertEqual(d.get_coordinate_derivative(1), 0)
-        self.assert_(d.get_coordinate_derivative(2) > 0)
-
-if __name__ == '__main__':
-    unittest.main()
Index: kernel/test/particle_refiners/test_refine_once_ps.py
===================================================================
--- kernel/test/particle_refiners/test_refine_once_ps.py	(revision 626)
+++ kernel/test/particle_refiners/test_refine_once_ps.py	(working copy)
@@ -9,7 +9,7 @@
 
 
     def test_rops(self):
-        """Make sure that bond cover coordinates are correct"""
+        """Test refine one pair score"""
         IMP.set_log_level(IMP.VERBOSE)
         m= IMP.Model()
         pp= IMP.Particle()
Index: kernel/test/particle_refiners/test_bond_pr.py
===================================================================
--- kernel/test/particle_refiners/test_bond_pr.py	(revision 626)
+++ kernel/test/particle_refiners/test_bond_pr.py	(working copy)
@@ -10,8 +10,8 @@
     def _set_up_stuff(self, n):
         # return [model, particles, bonds]
         m= IMP.Model()
-        ps= self.create_particles_in_box(m)
-        bds= self.create_chain(ps, 10)
+        ps= IMP.test.create_particles_in_box(m)
+        bds= IMP.test.create_chain(ps, 10)
         bl= IMP.BondDecoratorListScoreState(ps)
         ss= IMP.CoverBondsScoreState(bl, rk)
         m.add_score_state(bl)
Index: kernel/test/states/test_nonbonded_list.py
===================================================================
--- kernel/test/states/test_nonbonded_list.py	(revision 626)
+++ kernel/test/states/test_nonbonded_list.py	(working copy)
@@ -11,16 +11,16 @@
 
     def test_quadratic(self):
         """Test quadratic NBL"""
-        ss='IMP.AllNonbondedListScoreState(md, self.rk, IMP.AllNonbondedListScoreState.QUADRATIC)'
+        ss='IMP.AllNonbondedListScoreState(md)'
         #eval(ss)
-        self.do_test_all(ss)
+        self.do_test_all(ss, IMP.AllNonbondedListScoreState.QUADRATIC)
 
 
     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)
+        ss='IMP.AllNonbondedListScoreState(md)'
+        self.do_test_all(ss,IMP.AllNonbondedListScoreState.GRID )
 
     def test_bbox(self):
         """Test bbox NBL"""
@@ -29,18 +29,18 @@
         except:
             print "No CGAL support"
         else:
-            ss= 'IMP.AllNonbondedListScoreState(md, self.rk, IMP.AllNonbondedListScoreState.BBOX)'
-            self.do_test_all(ss)
+            ss= 'IMP.AllNonbondedListScoreState(md)'
+            self.do_test_all(ss, IMP.AllNonbondedListScoreState.BBOX)
 
     def test_default(self):
         """Test default NBL"""
-        ss= 'IMP.AllNonbondedListScoreState(md, self.rk)'
-        self.do_test_all(ss)
+        ss= 'IMP.AllNonbondedListScoreState(md)'
+        self.do_test_all(ss,  IMP.AllNonbondedListScoreState.DEFAULT)
 
     def test_bipartite_quadratic(self):
         """Test bipartite quadratic NBL"""
-        ss='IMP.BipartiteNonbondedListScoreState(md, self.rk, IMP.BipartiteNonbondedListScoreState.QUADRATIC)'
-        self.do_test_all_bi(ss)
+        ss='IMP.BipartiteNonbondedListScoreState(md)'
+        self.do_test_all_bi(ss, IMP.BipartiteNonbondedListScoreState.QUADRATIC)
 
     def test_bipartite_bbox(self):
         """Test bipartite bbox NBL"""
@@ -49,27 +49,27 @@
         except:
             print "No CGAL support"
         else:
-            ss= 'IMP.BipartiteNonbondedListScoreState(md, self.rk, IMP.BipartiteNonbondedListScoreState.BBOX)'
-            self.do_test_all_bi(ss)
+            ss= 'IMP.BipartiteNonbondedListScoreState(md)'
+            self.do_test_all_bi(ss, IMP.BipartiteNonbondedListScoreState.BBOX)
 
 
 
-    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(self, ss, algo):
+        self.do_test_bl(ss, algo)
+        self.do_test_distfilt(ss, algo)
+        self.do_test_spheres(ss, algo)
 
-    def do_test_all_bi(self, ss):
-        self.do_test_bi(ss)
-        self.do_test_bi_update(ss)
+    def do_test_all_bi(self, ss, algo):
+        self.do_test_bi(ss, algo)
+        self.do_test_bi_update(ss, algo)
 
     def make_spheres(self, m, num, lbv, ubv, minr, maxr):
-        ps=self.create_particles_in_box(m, num, lbv, ubv)
+        ps=IMP.test.create_particles_in_box(m, num, lbv, ubv)
         for p in ps:
             p.add_attribute(self.rk, random.uniform(minr, maxr), False)
         return ps
 
-    def do_test_update(self, ss):
+    def do_test_update(self, ss, algo):
         m= IMP.Model()
         self.make_spheres(m, 20, [0,0,0], [10,10,10], .1, 1)
         ds=[]
@@ -77,6 +77,7 @@
             ds.append(IMP.XYZDecorator.cast(p))
         md=5
         s=eval(ss)
+        s.set_algorithm(algo)
         s.set_particles(m.get_particles())
         m.add_score_state(s)
         o= IMP.test.ConstPairScore(1)
@@ -88,7 +89,7 @@
             for d in ds:
                 d.randomize_in_sphere(d.get_vector(), 2.0)
 
-    def do_test_bl(self, ss):
+    def do_test_bl(self, ss, algo):
         """Test the bond decorator list"""
         m= IMP.Model()
         bds=[]
@@ -99,6 +100,7 @@
             IMP.custom_bond(bds[i-1], bds[i], 1, .1)
         md= 20
         s=eval(ss)
+        s.set_algorithm(algo)
         s.set_particles(pts)
         b= IMP.BondDecoratorListScoreState(pts)
         s.add_bonded_list(b)
@@ -114,17 +116,18 @@
         print "Score with bonds is " + str(score)
         self.assertEqual(score, 190+1900-19, "Wrong score")
 
-    def do_test_distfilt(self, ss):
+    def do_test_distfilt(self, ss, algo):
         """Test filtering based on distance in nonbonded list"""
         #IMP.set_log_level(IMP.TERSE)
         m= IMP.Model()
         ps=IMP.Particles()
-        ps= self.create_particles_in_box(m, 20, [0,0,0], [10,10,10])
-        pts= self.create_particles_in_box(m, 20, [160,160,160], [170,170,170])
+        ps= IMP.test.create_particles_in_box(m, 20, [0,0,0], [10,10,10])
+        pts= IMP.test.create_particles_in_box(m, 20, [160,160,160], [170,170,170])
         for p in pts:
             ps.append(p)
         md=15
         s=eval(ss)
+        s.set_algorithm(algo)
         s.set_particles(ps)
         m.add_score_state(s)
         o= IMP.test.ConstPairScore(1)
@@ -157,13 +160,14 @@
         self.assertEqual(score, 2*190, "Wrong score after insert")
 
 
-    def do_test_bi(self, ss):
+    def do_test_bi(self, ss, algo):
         """Test the bipartite nonbonded list and restraint which uses it"""
         m= IMP.Model()
         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)
+        s.set_algorithm(algo)
         #IMP.QuadraticBipartiteNonbondedListScoreState(IMP.FloatKey(),
         #                                              ps0, ps1)
         s.set_particles(ps0, ps1)
@@ -175,13 +179,14 @@
         print score
         self.assertEqual(score, 25, "Wrong score")
 
-    def do_test_bi_update(self, ss):
+    def do_test_bi_update(self, ss, algo):
         """Test the bipartite nonbonded list and restraint which uses it"""
         m= IMP.Model()
         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)
+        s.set_algorithm(algo)
         #IMP.QuadraticBipartiteNonbondedListScoreState(IMP.FloatKey(),
         #                                              ps0, ps1)
         s.set_particles(ps0, ps1)
@@ -199,13 +204,14 @@
                 d.randomize_in_sphere(d.get_vector(), 1)
 
 
-    def do_test_spheres(self, ss):
+    def do_test_spheres(self, ss, algo):
         """Test the nonbonded list of spheres (collision detection)"""
         m= IMP.Model()
         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_algorithm(algo)
         s.set_particles(m.get_particles())
         #IMP.AllNonbondedListScoreState(rk, m.get_particles())
         m.add_score_state(s)
Index: kernel/test/optimizers/test_sd_optimizer.py
===================================================================
--- kernel/test/optimizers/test_sd_optimizer.py	(revision 626)
+++ kernel/test/optimizers/test_sd_optimizer.py	(working copy)
@@ -70,7 +70,7 @@
 
         # Start off with all particles in close proximity (but not actually
         # colocated, as the derivative of zero distance is zero):
-        self.randomize_particles(self.particles, .01)
+        IMP.test.randomize_particles(self.particles, .01)
 
         self.steepest_descent.optimize(50)
 
Index: kernel/test/decorators/test_yaml.py
===================================================================
--- kernel/test/decorators/test_yaml.py	(revision 0)
+++ kernel/test/decorators/test_yaml.py	(revision 0)
@@ -0,0 +1,47 @@
+import unittest
+import IMP
+import IMP.test
+
+class YamlTests(IMP.test.TestCase):
+    def _create_model(self):
+        IMP.set_log_level(IMP.VERBOSE)
+        m = IMP.Model()
+        p= IMP.Particle()
+        m.add_particle(p)
+        IMP.XYZDecorator.create(p)
+        p= IMP.Particle()
+        m.add_particle(p)
+        IMP.XYZDecorator.create(p)
+        return m
+    def test_yaml(self):
+        """Check writing to yaml """
+        m= self._create_model()
+        IMP.write_yaml(m)
+    def test_read(self):
+        """Check reading from yaml"""
+        m= self._create_model()
+        st="""particle: 0
+  float-attributes:
+    x: 0
+    y: 1
+    z: 2
+  int-attributes:
+  string-attributes:
+  particle-attributes:
+particle: 1
+  float-attributes:
+    x: 3
+    y: 4
+    z: 5
+  int-attributes:
+  string-attributes:
+  particle-attributes:
+"""
+        IMP.read_yaml(st, m)
+        IMP.write_yaml(m)
+
+
+
+
+if __name__ == '__main__':
+    unittest.main()
Index: kernel/test/pair_scores/test_typed_pair_score.py
===================================================================
--- kernel/test/pair_scores/test_typed_pair_score.py	(revision 626)
+++ kernel/test/pair_scores/test_typed_pair_score.py	(working copy)
@@ -7,9 +7,8 @@
 class TypedPairScoreTests(IMP.test.TestCase):
     """Class to test TypedPairScore"""
 
-    def _make_particles(self, types):
+    def _make_particles(self, m, types):
         """Make particles with the given types"""
-        m = IMP.Model()
         ps = [IMP.Particle() for i in types]
         for p, typ in zip(ps, types):
             m.add_particle(p)
@@ -21,15 +20,23 @@
         ps = IMP.TypedPairScore(typekey)
         cps = IMP.test.ConstPairScore(5)
         ps.set_pair_score(cps, 0, 1)
-        pa, pb = self._make_particles((0, 1))
+        # particles will be destroyed when the model  goes away
+        m = IMP.Model()
+        pa, pb = self._make_particles(m, (0, 1))
         da = IMP.DerivativeAccumulator()
         # The ordering of the particles should not matter:
         self.assertEqual(ps.evaluate(pa, pb, da), 5.0)
         self.assertEqual(ps.evaluate(pb, pa, da), 5.0)
+        print "ps"
+        ps=[]
+        print "cps"
+        cps= []
+        print "done"
 
     def test_invalid_type(self):
         """Check TypedPairScore behavior with invalid particle types"""
-        pa, pb = self._make_particles((0, 1))
+        m = IMP.Model()
+        pa, pb = self._make_particles(m, (0, 1))
         da = IMP.DerivativeAccumulator()
         ps1 = IMP.TypedPairScore(typekey, True)
         self.assertEqual(ps1.evaluate(pa, pb, da), 0.0)
Index: kernel/test/pair_scores/test_transform.py
===================================================================
--- kernel/test/pair_scores/test_transform.py	(revision 626)
+++ kernel/test/pair_scores/test_transform.py	(working copy)
@@ -71,7 +71,10 @@
         d0.set_coordinates_are_optimized(True)
         d1.set_coordinates_are_optimized(True)
         tps= IMP.TransformedDistancePairScore(IMP.Harmonic(0,1))
-        tps.set_translation(0,1,0)
+        cv= IMP.Vector3D(2,1,0)
+        tv= IMP.Vector3D(0,1,0)
+        tps.set_translation(tv[0], tv[1], tv[2])
+        tps.set_center(cv[0], cv[1], cv[2])
         tps.set_rotation(1, 0, 0,
                          0, 0,-1,
                          0, 1, 0)
@@ -83,14 +86,16 @@
         cg.optimize(100)
         d0.show()
         d1.show()
-        vt= IMP.Vector3D(d1.get_vector()*IMP.Vector3D(1,0,0)+0,
-                         d1.get_vector()*IMP.Vector3D(0,0,-1)+1,
-                         d1.get_vector()*IMP.Vector3D(0,1,0)+0)
+        ncv= d1.get_vector()- cv
+        vt= IMP.Vector3D(ncv*IMP.Vector3D(1,0,0),
+                         ncv*IMP.Vector3D(0,0,-1),
+                         ncv*IMP.Vector3D(0,1,0))
+        fvt= vt+ tv+cv
         print "trans"
-        print str(vt[0]) + " " + str(vt[1])+" " + str(vt[2])
-        self.assertEqual(vt[0], d0.get_coordinate(0))
-        self.assertEqual(vt[1], d0.get_coordinate(1))
-        self.assertEqual(vt[2], d0.get_coordinate(2))
+        print str(fvt[0]) + " " + str(fvt[1])+" " + str(fvt[2])
+        self.assertEqual(fvt[0], d0.get_coordinate(0))
+        self.assertEqual(fvt[1], d0.get_coordinate(1))
+        self.assertEqual(fvt[2], d0.get_coordinate(2))
 
 if __name__ == '__main__':
     unittest.main()
Index: kernel/test/connectivity/test_connectivity.py
===================================================================
--- kernel/test/connectivity/test_connectivity.py	(revision 626)
+++ kernel/test/connectivity/test_connectivity.py	(working copy)
@@ -26,8 +26,6 @@
         IMP.set_log_level(IMP.VERBOSE)
         m = IMP.Model()
 
-        rk= IMP.FloatKey("radius")
-
         p0=IMP.Particles()
         p1=IMP.Particles()
         p2=IMP.Particles()
@@ -35,28 +33,28 @@
         for i in range(13):
             p= IMP.Particle()
             m.add_particle(p)
-            d= IMP.XYZDecorator.create(p)
+            d= IMP.XYZRDecorator.create(p)
             d.set_coordinates_are_optimized(True)
             if i % 3 == 0:
                 p0.append(p)
-                p.add_attribute(rk, 1);
+                d.set_radius(1);
             elif i % 3 == 1:
                 p1.append(p)
-                p.add_attribute(rk, 2);
+                d.set_radius(2);
             else:
                 p2.append(p)
-                p.add_attribute(rk, 3);
+                d.set_radius(3);
         p= IMP.Particle()
         m.add_particle(p)
-        d= IMP.XYZDecorator.create(p)
+        d= IMP.XYZRDecorator.create(p)
         d.set_coordinates_are_optimized(True)
         p3.append(p)
-        p.add_attribute(rk, 1);
+        d.set_radius(1);
 
         o = IMP.ConjugateGradients()
         o.set_threshold(1e-4)
         o.set_model(m)
-        self.randomize_particles(m.get_particles(), 50.0)
+        IMP.test.randomize_particles(m.get_particles(), 50.0)
 
         # add connectivity restraints
 
Index: kernel/include/IMP/singleton_scores/DistanceToSingletonScore.h
===================================================================
--- kernel/include/IMP/singleton_scores/DistanceToSingletonScore.h	(revision 626)
+++ kernel/include/IMP/singleton_scores/DistanceToSingletonScore.h	(working copy)
@@ -27,7 +27,7 @@
   DistanceToSingletonScore(UnaryFunction *f, const Vector3D& pt);
   virtual ~DistanceToSingletonScore(){}
   virtual Float evaluate(Particle *a,
-                         DerivativeAccumulator *da);
+                         DerivativeAccumulator *da) const;
   virtual void show(std::ostream &out=std::cout) const;
 };
 
Index: kernel/include/IMP/singleton_scores/SConscript
===================================================================
--- kernel/include/IMP/singleton_scores/SConscript	(revision 626)
+++ kernel/include/IMP/singleton_scores/SConscript	(working copy)
@@ -1,9 +1,10 @@
+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',
+       'TunnelSingletonScore.h',
+      ]
+includedir = os.path.join(env['includedir'], 'IMP', 'singleton_scores' )
 inst = env.Install(includedir, files)
 env.Alias('install', inst)
Index: kernel/include/IMP/singleton_scores/TunnelSingletonScore.h
===================================================================
--- kernel/include/IMP/singleton_scores/TunnelSingletonScore.h	(revision 0)
+++ kernel/include/IMP/singleton_scores/TunnelSingletonScore.h	(revision 0)
@@ -0,0 +1,79 @@
+/**
+ *  \file TunnelSingletonScore.h    \brief Tunnel restraint.
+ *
+ *  Just return a constant.
+ *
+ *  Copyright 2007-8 Sali Lab. All rights reserved.
+ *
+ */
+
+#ifndef __IMP_TUNNEL_RESTRAINT_H
+#define __IMP_TUNNEL_RESTRAINT_H
+
+#include "../IMP_config.h"
+#include "../SingletonScore.h"
+#include "../Vector3D.h"
+#include "../UnaryFunction.h"
+#include "../Pointer.h"
+#include "../internal/kernel_version_info.h"
+
+namespace IMP
+{
+
+class PairScore;
+
+//! Score particles with respect to a tunnel.
+/** Particles with x,y,z coordinates and an optional radius are
+    prevented from being in a volume destribed by a slab from
+    (center[i]-height) to (center+height) on the ith coordinate with a
+    tunnel of radius r centered at center. To set which coordinate
+    is used use the get/set _coordinate functions.
+
+    Note that the UnaryFunction should look like a lower bound with 0
+    meaning that the particle is just touching the tunnel.
+
+    \ingroup restraint
+ */
+class IMPDLLEXPORT TunnelSingletonScore : public SingletonScore
+{
+  int coordinate_;
+  Vector3D center_;
+  Float height_;
+  Float radius_;
+  Pointer<UnaryFunction> f_;
+public:
+  TunnelSingletonScore(UnaryFunction* f);
+
+  void set_center(Vector3D c){
+    center_=c;
+  }
+  void set_height(Float h){
+    IMP_check(h >= 0,
+              "Height can't be negative",
+              ValueException);
+    height_=h;
+  }
+  void set_radius(Float h){
+    IMP_check(h >= 0,
+              "Radius can't be negative",
+              ValueException);
+    radius_=h;
+  }
+  void set_coordinate(unsigned int i) {
+    IMP_check(i < 3,
+              "Invalid coordinate value",
+              ValueException);
+    coordinate_=i;
+  }
+  unsigned int get_coordinate() const {
+    return coordinate_;
+  }
+
+  virtual Float evaluate(Particle *a,
+                         DerivativeAccumulator *da) const;
+  virtual void show(std::ostream &out=std::cout) const;
+};
+
+} // namespace IMP
+
+#endif  /* __IMP_TUNNEL_RESTRAINT_H */
Index: kernel/include/IMP/base_types.h
===================================================================
--- kernel/include/IMP/base_types.h	(revision 626)
+++ kernel/include/IMP/base_types.h	(working copy)
@@ -14,7 +14,6 @@
 
 #include <string>
 #include <vector>
-#include <map>
 
 namespace IMP
 {
Index: kernel/include/IMP/SingletonScore.h
===================================================================
--- kernel/include/IMP/SingletonScore.h	(revision 626)
+++ kernel/include/IMP/SingletonScore.h	(working copy)
@@ -34,7 +34,7 @@
   virtual ~SingletonScore() {}
   //! Compute the score for the particle and the derivative if needed.
   virtual Float evaluate(Particle *a,
-                         DerivativeAccumulator *da) = 0;
+                         DerivativeAccumulator *da) const = 0;
   virtual void show(std::ostream &out=std::cout) const = 0;
 };
 
Index: kernel/include/IMP/Model.h
===================================================================
--- kernel/include/IMP/Model.h	(revision 626)
+++ kernel/include/IMP/Model.h	(working copy)
@@ -28,6 +28,7 @@
 /** All attribute data for particles is stored through indexing in the
     model_data_ structure.
     Currently no suport for constraints (e.g. rigid bodies).
+
     \ingroup kernel
  */
 class IMPDLLEXPORT Model: public Object
Index: kernel/include/IMP/UnaryFunction.h
===================================================================
--- kernel/include/IMP/UnaryFunction.h	(revision 626)
+++ kernel/include/IMP/UnaryFunction.h	(working copy)
@@ -28,7 +28,7 @@
   /** \param[in] feature Value of feature being tested.
       \return Score
    */
-  virtual Float evaluate(Float feature) = 0;
+  virtual Float evaluate(Float feature) const = 0;
 
   //! Calculate score and derivative with respect to the given feature.
   /** \param[in] feature Value of feature being tested.
@@ -36,7 +36,7 @@
                         given feaure.
       \return Score
    */
-  virtual Float evaluate_deriv(Float feature, Float& deriv) = 0;
+  virtual Float evaluate_deriv(Float feature, Float& deriv) const = 0;
 
   virtual void show(std::ostream &out=std::cout) const = 0;
 };
Index: kernel/include/IMP/score_states/ManualBondDecoratorListScoreState.h
===================================================================
--- kernel/include/IMP/score_states/ManualBondDecoratorListScoreState.h	(revision 0)
+++ kernel/include/IMP/score_states/ManualBondDecoratorListScoreState.h	(revision 0)
@@ -0,0 +1,50 @@
+/**
+ *  \file BondDecoratorListScoreState.h
+ *  \brief Allow iteration through pairs of a set of atoms.
+ *
+ *  Copyright 2007-8 Sali Lab. All rights reserved.
+ */
+
+#ifndef __IMP_MANUAL_BOND_DECORATOR_LIST_SCORE_STATE_H
+#define __IMP_MANUAL_BOND_DECORATOR_LIST_SCORE_STATE_H
+
+#include "BondedListScoreState.h"
+#include "../decorators/bond_decorators.h"
+
+#include <vector>
+
+namespace IMP
+{
+
+class ManualBondDecoratorListScoreState;
+typedef Index<ManualBondDecoratorListScoreState> BondDecoratorListIndex;
+
+//! Keep track of particles that are connected by BondDecorator bonds.
+/** This class just keeps track of the bonds it has been told about.
+    Use the non-manual version if you want auto discovery of bonds.
+    \ingroup bond
+ */
+class IMPDLLEXPORT ManualBondDecoratorListScoreState: 
+    public BondedListScoreState
+{
+ public:
+  //! Find bonds amongst the following points. 
+  /** \param [in] ps The set of particles to use.
+   */
+  ManualBondDecoratorListScoreState(const BondDecorators &ps= BondDecorators());
+  virtual ~ManualBondDecoratorListScoreState(){}
+
+  virtual void set_particles(const Particles &ps) {}
+  virtual void add_particles(const Particles &ps) {}
+
+  IMP_LIST(public, BondDecorator, bond_decorator, BondDecorator);
+
+  virtual bool are_bonded(Particle *a, Particle *b) const;
+
+protected:
+  virtual void do_before_evaluate() {}
+};
+
+} // namespace IMP
+
+#endif  /* __IMP_BOND_DECORATOR_LIST_SCORE_STATE_H */
Index: kernel/include/IMP/score_states/BondDecoratorListScoreState.h
===================================================================
--- kernel/include/IMP/score_states/BondDecoratorListScoreState.h	(revision 626)
+++ kernel/include/IMP/score_states/BondDecoratorListScoreState.h	(working copy)
@@ -8,25 +8,21 @@
 #ifndef __IMP_BOND_DECORATOR_LIST_SCORE_STATE_H
 #define __IMP_BOND_DECORATOR_LIST_SCORE_STATE_H
 
-#include "BondedListScoreState.h"
+#include "ManualBondDecoratorListScoreState.h"
 #include "../decorators/bond_decorators.h"
 
 #include <vector>
 
 namespace IMP
 {
-
-class BondDecoratorListScoreState;
-typedef Index<BondDecoratorListScoreState> BondDecoratorListIndex;
-
 //! Keep track of particles that are connected by BondDecorator bonds.
 /** We also may want to add lazy rescanning of bonds rather than doing
     it every update call and a faster lookup of bonds.
     \ingroup bond
  */
-class IMPDLLEXPORT BondDecoratorListScoreState: public BondedListScoreState
+class IMPDLLEXPORT BondDecoratorListScoreState:
+    public ManualBondDecoratorListScoreState
 {
-  std::vector<BondDecorator> bonds_;
   Particles ps_;
 public:
   //! Find bonds amongst the following points. 
@@ -36,24 +32,8 @@
   virtual ~BondDecoratorListScoreState(){}
 
   virtual void set_particles(const Particles &ps);
+  virtual void add_particles(const Particles &ps);
 
-  virtual bool are_bonded(Particle *a, Particle *b) const;
-
-  //! This iterates through the pairs of bonded particles
-  /** \note update() must be called first for this to be valid.
-   */
-  typedef std::vector<BondDecorator>::const_iterator BondIterator;
-  BondIterator bonds_begin() const {
-    return bonds_.begin();
-  }
-  BondIterator bonds_end() const {
-    return bonds_.end();
-  }
-
-  unsigned int get_number_of_bonds() const {
-    return bonds_.size();
-  }
-
 protected:
   virtual void do_before_evaluate();
 };
Index: kernel/include/IMP/score_states/NonbondedListScoreState.h
===================================================================
--- kernel/include/IMP/score_states/NonbondedListScoreState.h	(revision 626)
+++ kernel/include/IMP/score_states/NonbondedListScoreState.h	(working copy)
@@ -9,6 +9,7 @@
 #define __IMP_NONBONDED_LIST_SCORE_STATE_H
 
 #include "../ScoreState.h"
+#include "MaxChangeScoreState.h"
 #include "../internal/kernel_version_info.h"
 #include "../decorators/XYZDecorator.h"
 #include "BondedListScoreState.h"
@@ -16,7 +17,6 @@
 #include <boost/iterator/filter_iterator.hpp>
 
 #include <vector>
-#include <limits>
 
 namespace IMP
 {
@@ -51,6 +51,7 @@
    */
   unsigned int max_nbl_size_;
   ParticlePairs nbl_;
+  Pointer<MaxChangeScoreState> mcr_;
 
   struct NBLTooLargeException{};
 
@@ -184,12 +185,22 @@
   BoxesOverlap boxes_overlap_object(float cut) const {
     return BoxesOverlap(get_radius_object(), cut);
   }
+
+  Float get_max_radius_change() const;
+
+  void reset_max_radius_change();
+
+  void update_max_radius_change(unsigned int i);
+
+  void add_particles_for_max_radius_change(const Particles &ps);
+  void clear_particles_for_max_radius_change();
+
 public:
-  NonbondedListScoreState(Float cutoff, FloatKey rk);
+  NonbondedListScoreState(Float cutoff);
   ~NonbondedListScoreState();
 
   FloatKey get_radius_key() const {return rk_;}
-  void set_radius_key(FloatKey rk) {rk_=rk;} 
+  void set_radius_key(FloatKey rk);
 
   //! Set the maximum allowable size for the NBL
   /** The NBL will keep reducing the slack and trying to
Index: kernel/include/IMP/score_states/SConscript
===================================================================
--- kernel/include/IMP/score_states/SConscript	(revision 626)
+++ kernel/include/IMP/score_states/SConscript	(working copy)
@@ -1,12 +1,16 @@
+Import('env')
 import os.path
-Import('env')
-
-files = ['BondedListScoreState.h', 'NonbondedListScoreState.h',
-         'BipartiteNonbondedListScoreState.h', 'MaxChangeScoreState.h',
-         'BondDecoratorListScoreState.h', 'AllNonbondedListScoreState.h',
-         'GravityCenterScoreState.h', 'CoverBondsScoreState.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',
+       'ManualBondDecoratorListScoreState.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/CoverBondsScoreState.h
===================================================================
--- kernel/include/IMP/score_states/CoverBondsScoreState.h	(revision 626)
+++ kernel/include/IMP/score_states/CoverBondsScoreState.h	(working copy)
@@ -9,7 +9,7 @@
 #define __IMP_COVER_BONDS_SCORE_STATE_H
 
 #include "../ScoreState.h"
-#include "BondDecoratorListScoreState.h"
+#include "ManualBondDecoratorListScoreState.h"
 #include "../internal/kernel_version_info.h"
 #include "../Pointer.h"
 
@@ -30,13 +30,13 @@
  */
 class IMPDLLEXPORT CoverBondsScoreState: public ScoreState
 {
-  Pointer<BondDecoratorListScoreState> bl_;
+  Pointer<ManualBondDecoratorListScoreState> 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,
+  CoverBondsScoreState(ManualBondDecoratorListScoreState *bl,
                        FloatKey rk=FloatKey("radius"));
   ~CoverBondsScoreState();
 
Index: kernel/include/IMP/score_states/BipartiteNonbondedListScoreState.h
===================================================================
--- kernel/include/IMP/score_states/BipartiteNonbondedListScoreState.h	(revision 626)
+++ kernel/include/IMP/score_states/BipartiteNonbondedListScoreState.h	(working copy)
@@ -20,7 +20,9 @@
 /** To iterate through the list of pairs use the NonbondedListScoreState::begin,
     NonbondedListScoreState::end functions.
 
-    The radius key can be the default key.
+    By default the XYZRDecorator radius is used as the radius. The radius key
+    can be set to FloatKey to force all particles to be point particles.
+    Particles without radius attribute are treated as points.
 
     Changes in coordinates and radii are handled properly.
 
@@ -43,7 +45,7 @@
   };
 protected:
   Algorithm a_;
-  Pointer<MaxChangeScoreState> mc0_, mc1_, mcr_;
+  Pointer<MaxChangeScoreState> mc0_, mc1_;
 
   void process_sets(const Particles &p0,
                     const Particles &p1);
@@ -57,8 +59,7 @@
       \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);
+  BipartiteNonbondedListScoreState(Float cutoff);
   /** \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.
@@ -66,10 +67,8 @@
       \param[in] a Which algorithm to use. The default is the best.
    */ 
   BipartiteNonbondedListScoreState(Float cutoff,
-                                   FloatKey radius,
                                    const Particles &ps0,
-                                   const Particles &ps1,
-                                   Algorithm a=DEFAULT);
+                                   const Particles &ps1);
 
   ~BipartiteNonbondedListScoreState();
   IMP_SCORE_STATE(internal::kernel_version_info);
@@ -78,6 +77,12 @@
   void add_particles_0(const Particles &ps);
   //! Add the particles to the second set
   void add_particles_1(const Particles &ps);
+
+  //! Add the particles to the first set
+  void add_particle_0(Particle* ps);
+  //! Add the particles to the second set
+  void add_particle_1(Particle* ps);
+
   //! Remove all particles
   void clear_particles();
   //! Replace the set of particles
Index: kernel/include/IMP/score_states/AllNonbondedListScoreState.h
===================================================================
--- kernel/include/IMP/score_states/AllNonbondedListScoreState.h	(revision 626)
+++ kernel/include/IMP/score_states/AllNonbondedListScoreState.h	(working copy)
@@ -26,7 +26,9 @@
 /** To iterate through the list of pairs use the NonbondedListScoreState::begin,
     NonbondedListScoreState::end functions.
 
-    The radius key can be the default key.
+    By default the XYZRDecorator radius is used as the radius. The radius key
+    can be set to FloatKey to force all particles to be point particles.
+    Particles without radius attribute are treated as points.
 
     Changes in coordinates and radii are handled properly unless grid is used,
     then changes in radii may not be handled properly.
@@ -50,7 +52,7 @@
   };
 protected:
   Algorithm a_;
-  Pointer<MaxChangeScoreState> mc_, mcr_;
+  Pointer<MaxChangeScoreState> mc_;
 
   void check_nbl() const;
 
@@ -72,20 +74,12 @@
 
   /** \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(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);
+                             const Particles &ps= Particles());
+
   ~AllNonbondedListScoreState();
 
   IMP_SCORE_STATE(internal::kernel_version_info);
@@ -97,6 +91,11 @@
   //! Replace the set of particles
   void set_particles(const Particles &ps);
 
+  //! hack
+  void add_particle(Particle *p) {
+    add_particles(Particles(1,p));
+  }
+
   //! If there is CGAL support, a more efficient algorithm BBOX can be used
   void set_algorithm(Algorithm a);
 };
Index: kernel/include/IMP/score_states/MaxChangeScoreState.h
===================================================================
--- kernel/include/IMP/score_states/MaxChangeScoreState.h	(revision 626)
+++ kernel/include/IMP/score_states/MaxChangeScoreState.h	(working copy)
@@ -20,13 +20,14 @@
 //! Keeps track of the maximum change of a set of attributes.
 /** The score state maintains a list of particle and a list of
     float attribute keys and keeps track of the maximum amount
-    any of these have changed since the last time reset was called.
-
+    any of these have changed since the last time reset was called.    
  */
 class IMPDLLEXPORT MaxChangeScoreState: public ScoreState
 {
   FloatKeys keys_;
-  FloatKeys origkeys_;
+  typedef internal::AttributeTable<internal::FloatAttributeTableTraits>
+    Table;
+  std::vector<Table> old_values_;
   float max_change_;
 public:
   //! Track the changes with the specified keys.
Index: kernel/include/IMP/restraints/SingletonListRestraint.h
===================================================================
--- kernel/include/IMP/restraints/SingletonListRestraint.h	(revision 626)
+++ kernel/include/IMP/restraints/SingletonListRestraint.h	(working copy)
@@ -38,6 +38,7 @@
   using Restraint::add_particles;
   using Restraint::clear_particles;
   using Restraint::set_particles;
+  using Restraint::add_particle;
 
 protected:
   Pointer<SingletonScore> ss_;
Index: kernel/include/IMP/restraints/PairChainRestraint.h
===================================================================
--- kernel/include/IMP/restraints/PairChainRestraint.h	(revision 626)
+++ kernel/include/IMP/restraints/PairChainRestraint.h	(working copy)
@@ -20,31 +20,26 @@
 namespace IMP
 {
 
-//! Restrain each pair of consecutive particles in each chain.
+//! Restrain each pair of consecutive particles in chain.
 /** \ingroup restraint
  */
 class IMPDLLEXPORT PairChainRestraint : public Restraint
 {
 public:
-  //! Create the pair restraint.
+  //! Create the chain restraint.
   /** \param[in] pair_score Pair score to apply.
    */
   PairChainRestraint(PairScore* pair_score);
   virtual ~PairChainRestraint(){}
 
-  IMP_RESTRAINT(internal::kernel_version_info)
+  IMP_RESTRAINT(internal::kernel_version_info);
 
-  //! Add a chain of particles
-  /** Each two successive particles are restrained.
-   */
-  void add_chain(const Particles &ps);
+  using Restraint::add_particles;
+  using Restraint::set_particles;
+  using Restraint::clear_particles;
 
-  //! Clear all the stored chains
-  void clear_chains();
-
 protected:
   Pointer<PairScore> ts_;
-  std::vector<unsigned int> chain_splits_;
 };
 
 } // namespace IMP
Index: kernel/include/IMP/restraints/SConscript
===================================================================
--- kernel/include/IMP/restraints/SConscript	(revision 626)
+++ kernel/include/IMP/restraints/SConscript	(working copy)
@@ -1,15 +1,20 @@
+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',
-         'ConstantRestraint.h', 'TunnelRestraint.h']
-
-# Install the include files:
-includedir = os.path.join(env['includedir'], 'IMP', 'restraints')
+files=[
+       'AngleRestraint.h',
+       'BondDecoratorRestraint.h',
+       'ConnectivityRestraint.h',
+       'ConstantRestraint.h',
+       'DihedralRestraint.h',
+       'DistanceRestraint.h',
+       'NonbondedRestraint.h',
+       'PairChainRestraint.h',
+       'PairListRestraint.h',
+       'PairRestraint.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/TunnelRestraint.h
===================================================================
--- kernel/include/IMP/restraints/TunnelRestraint.h	(revision 626)
+++ kernel/include/IMP/restraints/TunnelRestraint.h	(working copy)
@@ -1,83 +0,0 @@
-/**
- *  \file TunnelRestraint.h    \brief Tunnel restraint.
- *
- *  Just return a constant.
- *
- *  Copyright 2007-8 Sali Lab. All rights reserved.
- *
- */
-
-#ifndef __IMP_TUNNEL_RESTRAINT_H
-#define __IMP_TUNNEL_RESTRAINT_H
-
-#include "../IMP_config.h"
-#include "../Restraint.h"
-#include "../Vector3D.h"
-#include "../UnaryFunction.h"
-#include "../internal/kernel_version_info.h"
-
-namespace IMP
-{
-
-class PairScore;
-
-//! Restrain particles to a tunnel.
-/** Particles with x,y,z coordinates and an optional radius are
-    prevented from being in a volume destribed by a slab from
-    (center[i]-height) to (center+height) on the ith coordinate with a
-    tunnel of radius r centered at center. To set which coordinate
-    is used use the get/set _coordinate functions.
-
-    Note that the UnaryFunction should look like a lower bound with 0
-    meaning that the particle is just touching the tunnel.
-
-    \ingroup restraint
- */
-class IMPDLLEXPORT TunnelRestraint : public Restraint
-{
-  int coordinate_;
-  Vector3D center_;
-  Float height_;
-  Float radius_;
-  Pointer<UnaryFunction> f_;
-  FloatKey rk_;
-public:
-  TunnelRestraint(UnaryFunction* f, FloatKey rk);
-
-  void set_center(Vector3D c){
-    center_=c;
-  }
-  void set_height(Float h){
-    IMP_check(h >= 0,
-              "Height can't be negative",
-              ValueException);
-    height_=h;
-  }
-  void set_radius(Float h){
-    IMP_check(h >= 0,
-              "Radius can't be negative",
-              ValueException);
-    radius_=h;
-  }
-  void set_coordinate(unsigned int i) {
-    IMP_check(i < 3,
-              "Invalid coordinate value",
-              ValueException);
-    coordinate_=i;
-  }
-  unsigned int get_coordinate() const {
-    return coordinate_;
-  }
-
-  using Restraint::add_particles;
-  using Restraint::add_particle;
-  using Restraint::set_particles;
-  using Restraint::clear_particles;
-  using Restraint::erase_particle;
-
-  IMP_RESTRAINT(internal::kernel_version_info)
-};
-
-} // namespace IMP
-
-#endif  /* __IMP_TUNNEL_RESTRAINT_H */
Index: kernel/include/IMP/restraints/NonbondedRestraint.h
===================================================================
--- kernel/include/IMP/restraints/NonbondedRestraint.h	(revision 626)
+++ kernel/include/IMP/restraints/NonbondedRestraint.h	(working copy)
@@ -16,12 +16,11 @@
 #include "../internal/kernel_version_info.h"
 #include "../Pointer.h"
 #include "../score_states/NonbondedListScoreState.h"
+#include "../PairScore.h"
 
 namespace IMP
 {
 
-class PairScore;
-
 //! Apply a PairScore to all nonbonded pairs of particles
 /**
    \ingroup restraint
Index: kernel/include/IMP/restraints/BondDecoratorRestraint.h
===================================================================
--- kernel/include/IMP/restraints/BondDecoratorRestraint.h	(revision 626)
+++ kernel/include/IMP/restraints/BondDecoratorRestraint.h	(working copy)
@@ -16,12 +16,11 @@
 #include "../internal/kernel_version_info.h"
 #include "../Pointer.h"
 #include "../UnaryFunction.h"
+#include "../score_states/ManualBondDecoratorListScoreState.h"
 
 namespace IMP
 {
 
-class BondDecoratorListScoreState;
-
 //! Restrain all pairs of non-bonded particles
 /** This restraint currently only works for bonds which have their 
     length set explicitly. Eventually we should add a table for standard
@@ -40,7 +39,8 @@
       \param[in] bl The BondDecoratorListScoreState to use to get the list
       of bonds.
    */
-  BondDecoratorRestraint(UnaryFunction *f, BondDecoratorListScoreState *bl);
+  BondDecoratorRestraint(UnaryFunction *f,
+                         ManualBondDecoratorListScoreState *bl);
   virtual ~BondDecoratorRestraint(){}
 
   IMP_RESTRAINT(internal::kernel_version_info)
@@ -48,7 +48,7 @@
   void set_function(UnaryFunction *f) {f_=f;}
 
 protected:
-  BondDecoratorListScoreState *bl_;
+  Pointer<ManualBondDecoratorListScoreState> bl_;
   Pointer<UnaryFunction> f_;
 };
 
Index: kernel/include/IMP/restraints/PairListRestraint.h
===================================================================
--- kernel/include/IMP/restraints/PairListRestraint.h	(revision 626)
+++ kernel/include/IMP/restraints/PairListRestraint.h	(working copy)
@@ -39,6 +39,10 @@
   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);
+  }
 
 protected:
   Pointer<PairScore> ss_;
Index: kernel/include/IMP/restraints/DistanceRestraint.h
===================================================================
--- kernel/include/IMP/restraints/DistanceRestraint.h	(revision 626)
+++ kernel/include/IMP/restraints/DistanceRestraint.h	(working copy)
@@ -10,7 +10,7 @@
 
 #include "../IMP_config.h"
 #include "../pair_scores/DistancePairScore.h"
-#include "../Restraint.h"
+#include "PairRestraint.h"
 #include "../internal/kernel_version_info.h"
 
 #include <iostream>
@@ -23,7 +23,7 @@
    \note If the particles are closer than a certain distance, then
    the contributions to the derivatives are set to 0.
  */
-class IMPDLLEXPORT DistanceRestraint : public Restraint
+class IMPDLLEXPORT DistanceRestraint : public PairRestraint
 {
 public:
   //! Create the distance restraint.
@@ -35,11 +35,7 @@
                     Particle* p1, Particle* p2);
   virtual ~DistanceRestraint() {}
 
-  IMP_RESTRAINT(internal::kernel_version_info)
-
-protected:
-  //! scoring function for this restraint
-  DistancePairScore dp_;
+  void show(std::ostream &out=std::cout) const;
 };
 
 } // namespace IMP
Index: kernel/include/IMP/restraints/PairRestraint.h
===================================================================
--- kernel/include/IMP/restraints/PairRestraint.h	(revision 0)
+++ kernel/include/IMP/restraints/PairRestraint.h	(revision 0)
@@ -0,0 +1,47 @@
+/**
+ *  \file PairRestraint.h   \brief Pair restraint between two particles.
+ *
+ *  Copyright 2007-8 Sali Lab. All rights reserved.
+ *
+ */
+
+#ifndef __IMP_PAIR_RESTRAINT_H
+#define __IMP_PAIR_RESTRAINT_H
+
+#include "../IMP_config.h"
+#include "../PairScore.h"
+#include "../Restraint.h"
+#include "../Pointer.h"
+#include "../internal/kernel_version_info.h"
+
+#include <iostream>
+
+namespace IMP
+{
+
+//! Pair restraint between two particles
+/**
+   Apply a pair score to a pair of particles.
+ */
+class IMPDLLEXPORT PairRestraint : public Restraint
+{
+public:
+  //! Create the pair restraint.
+  /** \param[in] score_func Scoring function for the restraint.
+      \param[in] p1 Pointer to first particle in distance restraint.
+      \param[in] p2 Pointer to second particle in distance restraint.
+   */
+  PairRestraint(PairScore* score_func,
+                Particle* p1, Particle* p2);
+  virtual ~PairRestraint() {}
+
+  IMP_RESTRAINT(internal::kernel_version_info)
+
+protected:
+  //! scoring function for this restraint
+  Pointer<PairScore> dp_;
+};
+
+} // namespace IMP
+
+#endif  /* __IMP_PAIR_RESTRAINT_H */
Index: kernel/include/IMP/Restraint.h
===================================================================
--- kernel/include/IMP/Restraint.h	(revision 626)
+++ kernel/include/IMP/Restraint.h	(working copy)
@@ -51,6 +51,11 @@
 
     \note Physical restraints should use the units of kcal/mol for restraint
     values and kcal/mol/A for derivatives.
+
+    \note Restraints will print a warning message if they are destroyed
+    without ever having been added to a model as this is an easy mistake
+    to make. To disable this warning for a particular restraint, call
+    set_was_owned(true).
  */
 class IMPDLLEXPORT Restraint : public RefCountedObject
 {
@@ -94,18 +99,31 @@
   Model *get_model() const {
     IMP_assert(model_,
                "get_model() called before set_model()");
-    return model_.get();
+    return model_;
   }
 
+  /** A warning is printed if a restraint is destroyed
+      without ever having belonged to a restraint set or a model.
+   */
+  void set_was_owned(bool tf) {
+    was_owned_=tf;
+  }
+
   IMP_LIST(protected, Particle, particle, Particle*)
 
 private:
-  Pointer<Model> model_;
+  // not ref counted as the model has a pointer to this
+  Model *model_;
 
   /* True if restraint has not been deactivated.
      If it is not active, evaluate should not be called
    */
   bool is_active_;
+
+  /* keep track of whether the restraint ever was in a model.
+     Give warnings on destruction if it was not.
+   */
+  bool was_owned_;
 };
 
 IMP_OUTPUT_OPERATOR(Restraint);
Index: kernel/include/IMP/ParticleRefiner.h
===================================================================
--- kernel/include/IMP/ParticleRefiner.h	(revision 626)
+++ kernel/include/IMP/ParticleRefiner.h	(working copy)
@@ -34,7 +34,7 @@
   //! Refine the passed particle into a set of particles.
   /** As a precondition can_refine_particle(a) should be true.
    */
-  virtual Particles get_refined(Particle *a);
+  virtual Particles get_refined(Particle *a) const;
 
   //! Cleanup after refining
   /** If da is non-NULL then the derivatives should be propagated
@@ -43,7 +43,7 @@
       can be destroyed if they are temporary.
    */
   virtual void cleanup_refined(Particle *a, Particles &b,
-                               DerivativeAccumulator *da=0) {}
+                               DerivativeAccumulator *da=0) const {}
   virtual void show(std::ostream &out=std::cout) const {
     out << "ParticleRefiner base" << std::endl;
   };
Index: kernel/include/IMP/optimizers/states/SConscript
===================================================================
--- kernel/include/IMP/optimizers/states/SConscript	(revision 626)
+++ 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/SConscript
===================================================================
--- kernel/include/IMP/optimizers/SConscript	(revision 626)
+++ 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', 'BrownianDynamics.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 626)
+++ 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/decorators/XYZDecorator.h
===================================================================
--- kernel/include/IMP/decorators/XYZDecorator.h	(revision 626)
+++ kernel/include/IMP/decorators/XYZDecorator.h	(working copy)
@@ -8,16 +8,13 @@
 #ifndef __IMP_XYZ_DECORATOR_H
 #define __IMP_XYZ_DECORATOR_H
 
-#include <vector>
-#include <deque>
-#include <limits>
-
-#include "../Particle.h"
-#include "../Model.h"
 #include "../DecoratorBase.h"
 #include "../Vector3D.h"
 #include "utility.h"
 
+#include <vector>
+#include <limits>
+
 namespace IMP
 {
 
@@ -110,10 +107,7 @@
   //! Get a vector containing the keys for x,y,z
   /** This is quite handy for initializing movers and things.
    */
-  static const FloatKeys get_xyz_keys() {
-    decorator_initialize_static_data();
-    return key_;
-  }
+  IMP_DECORATOR_GET_KEY(FloatKeys, xyz_keys, key_)
 
   //! Generate random coordinates in a sphere centered at the vector
   void randomize_in_sphere(const Vector3D &center, float radius);
Index: kernel/include/IMP/decorators/XYZRDecorator.h
===================================================================
--- kernel/include/IMP/decorators/XYZRDecorator.h	(revision 0)
+++ kernel/include/IMP/decorators/XYZRDecorator.h	(revision 0)
@@ -0,0 +1,63 @@
+/**
+ *  \file XYZRDecorator.h     \brief Simple xyz decorator.
+ *
+ *  Copyright 2007-8 Sali Lab. All rights reserved.
+ *
+ */
+
+#ifndef __IMP_XYZR_DECORATOR_H
+#define __IMP_XYZR_DECORATOR_H
+
+#include "XYZDecorator.h"
+
+#include <vector>
+#include <deque>
+#include <limits>
+
+namespace IMP
+{
+
+//! A a decorator for a particle with x,y,z coordinates.
+/** \ingroup helper
+ */
+class IMPDLLEXPORT XYZRDecorator: public XYZDecorator
+{
+  IMP_DECORATOR(XYZRDecorator, XYZDecorator,
+                {
+                  return p->has_attribute(radius_);
+                },
+                { p->add_attribute(radius_,0);
+                });
+
+protected:
+  static FloatKey radius_;
+
+public:
+  IMP_DECORATOR_GET_SET(radius, radius_, Float, Float);
+  IMP_DECORATOR_GET_KEY(FloatKey, radius_key, radius_);
+};
+
+IMP_OUTPUT_OPERATOR(XYZRDecorator);
+
+//! Compute the distance between a pair of particles
+/** \ingroup helper
+ */
+IMPDLLEXPORT Float distance(XYZRDecorator a, XYZRDecorator b);
+
+//! Set the coordinates and radius of the last to enclose the list
+/** \ingroup helper
+    \param[in] v The vector of XYZR objects to enclose
+    \param[out] b The one whose values should be set
+
+    \note this takes a Particles object rather than a vector of
+    something else since you can't easily cast vectors of
+    different things to one another. Well, you can cast vectors
+    of decorators in C++, but you can't in python without adding
+    explicit support.
+ */
+IMPDLLEXPORT void set_enclosing_sphere(const Particles &v,
+                                       XYZRDecorator b);
+
+} // namespace IMP
+
+#endif  /* __IMP_XYZR_DECORATOR_H */
Index: kernel/include/IMP/decorators/HierarchyDecorator.h
===================================================================
--- kernel/include/IMP/decorators/HierarchyDecorator.h	(revision 626)
+++ kernel/include/IMP/decorators/HierarchyDecorator.h	(working copy)
@@ -19,6 +19,7 @@
 #include <limits>
 #include <vector>
 #include <deque>
+#include <memory>
 
 namespace IMP
 {
@@ -135,14 +136,14 @@
     \ingroup hierarchy
  */
 IMPDLLEXPORT
-void breadth_first_traversal(HierarchyDecorator d,  HierarchyVisitor &v);
+void breadth_first_traversal(HierarchyDecorator d,  HierarchyVisitor *v);
 
 //! Depth first traversal of the hierarchy
 /** See breadth_first_traversal and HierarchyVisitor for more information
     \ingroup hierarchy
  */
 IMPDLLEXPORT
-void depth_first_traversal(HierarchyDecorator d,  HierarchyVisitor &v);
+void depth_first_traversal(HierarchyDecorator d,  HierarchyVisitor *v);
 
 
 //! Apply functor F to each particle, traversing the hierarchy breadth first.
@@ -313,23 +314,35 @@
 
 };
 
+/** This returns an auto_ptr so that the objects are cleaned up properly.
+    I have to pass a non-const ref (or a pointer) to the algorithms which 
+    are getting the gather object. Returning an autoptr ensures that the
+    object is alive until after the function call has completed (if gather
+    is called in the parens for the function).
+
+   \internal
+*/
+template <class F, class Out>
+std::auto_ptr<Gather<F, Out> > gather(F f, Out out) {
+  return std::auto_ptr<Gather<F, Out> >(new Gather<F, Out>(f, out));
+}
+
 } // namespace internal
 
 //! Gather all the Particle* in the hierarchy which meet some criteria
 /** \ingroup hierarchy
  */
-template <class Out, class F>
-Out hierarchy_gather(HierarchyDecorator h, F f, Out out)
+template <class F>
+Particles hierarchy_get(HierarchyDecorator h, F f)
 {
-  internal::Gather<F,Out> gather(f,out);
-  depth_first_traversal(h, gather);
-  return gather.get_out();
+  Particles ret;
+  depth_first_traversal(h, internal::gather(f,std::back_inserter(ret)).get());
+  return ret;
 }
 
 
 namespace internal
 {
-
 template <class K, class V>
 struct MatchAttribute
 {
@@ -341,20 +354,19 @@
     else return o->get_value(k_) == v_;
   }
 };
+}
 
-} // namespace internal
-
 //! Gather all the Particle* in the hierarchy which match on an attribute
 /** \ingroup hierarchy
  */
-template <class Out, class K, class V>
-Out hierarchy_gather_by_attribute(HierarchyDecorator h, K k, V v, Out out)
+template <class K, class V>
+Particles hierarchy_get_by_attribute(HierarchyDecorator h, K k,
+                                     V v)
 {
-  internal::Gather<internal::MatchAttribute<K, V>,Out>
-    gather(internal::MatchAttribute<K,V>(k,v),
-           out);
-  depth_first_traversal(h, gather);
-  return gather.get_out();
+  Particles ret;
+  depth_first_traversal(h, gather(internal::MatchAttribute<K,V>(k,v),
+                                  std::back_inserter(ret)).get());
+  return ret;
 }
 
 
@@ -380,20 +392,27 @@
   }
 };
 
+template <class K0, class V0, class K1, class V1>
+MatchAttributes<K0, V0, K1, V1>
+match_attributes(K0 k0, V0 v0,
+                 K1 k1, V1 v1) {
+  return MatchAttributes<K0, V0, K1, V1>(k0, v0, k1, v1);
+}
+
 } // namespace internal
 
 //! Gather all the Particle* in the hierarchy which match on two attributes
 /** \ingroup hierarchy
  */
-template <class Out, class K0, class V0, class K1, class V1>
-Out hierarchy_gather_by_attributes(HierarchyDecorator h, K0 k0,
-                                   V0 v0, K1 k1, V1 v1, Out out)
+template <class K0, class V0, class K1, class V1>
+Particles hierarchy_get_by_attributes(HierarchyDecorator h, K0 k0,
+                                      V0 v0, K1 k1, V1 v1)
 {
-  internal::Gather<internal::MatchAttributes<K0, V0, K1, V1>,Out>
-    gather(internal::MatchAttributes<K0,V0, K1, V1>(k0,v0, k1, v1),
-           out);
-  depth_first_traversal(h, gather);
-  return gather.get_out();
+  Particles ret;
+  depth_first_traversal(h,
+                        gather(internal::match_attributes(k0,v0, k1, v1),
+                               std::back_inserter(ret)).get());
+  return ret;
 }
 
 
Index: kernel/include/IMP/decorators/SConscript
===================================================================
--- kernel/include/IMP/decorators/SConscript	(revision 626)
+++ kernel/include/IMP/decorators/SConscript	(working copy)
@@ -1,11 +1,18 @@
+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',
+       'XYZRDecorator.h',
+       'yaml.h',
+      ]
+includedir = os.path.join(env['includedir'], 'IMP', 'decorators' )
 inst = env.Install(includedir, files)
 env.Alias('install', inst)
Index: kernel/include/IMP/decorators/yaml.h
===================================================================
--- kernel/include/IMP/decorators/yaml.h	(revision 0)
+++ kernel/include/IMP/decorators/yaml.h	(revision 0)
@@ -0,0 +1,49 @@
+#ifndef IMP_DECORATORS_YAML_H
+#define IMP_DECORATORS_YAML_H
+#include "../Particle.h"
+#include "../Model.h"
+
+namespace IMP
+{
+
+  //! Write the particle to a stream as a YAML stream
+  /** \param[in] p The particle to write
+      \param[in] out The stream to write to
+      \param[in] indent The base level of indentation
+   */
+  IMPDLLEXPORT void write_yaml(Particle *p,
+                               std::ostream &out= std::cout,
+                               std::string indent="");
+
+//! Write the model to a stream as a YAML stream
+  /** \param[in] m The model to write
+      \param[in] out The stream to write to
+      \param[in] indent The base level of indentation
+   */
+  IMPDLLEXPORT void write_yaml(Model *m,
+                               std::ostream &out= std::cout,
+                               std::string indent="");
+
+  //! Read the particle from a YAML stream
+  /** The model must already have particles matching all read particles.
+      Currently the particles must already have the same attributes
+      as are being read, but this probably should change due to lists
+      stored in attributes.
+
+      The intentended usage model is that the model (with restraints) is
+      initialized. Optimization is performed and then the model is written
+      out to a file. Some time later, when you want to reload the model, 
+      you can reuse the initialization code to set up the restraints,
+      and then read in the values for the attributes.
+
+      \note The base indent is determined from the first line. */
+  IMPDLLEXPORT void read_yaml(std::istream &in,
+                              Model *m);
+
+  //! primarily for debuggin*/
+  IMPDLLEXPORT void read_yaml(std::string contents,
+                              Model *m);
+
+}
+
+#endif
Index: kernel/include/IMP/decorators/AtomDecorator.h
===================================================================
--- kernel/include/IMP/decorators/AtomDecorator.h	(revision 626)
+++ kernel/include/IMP/decorators/AtomDecorator.h	(working copy)
@@ -14,9 +14,6 @@
 #include "utility.h"
 #include "XYZDecorator.h"
 
-#include <vector>
-#include <deque>
-
 namespace IMP
 {
 
Index: kernel/include/IMP/decorators/macros.h
===================================================================
--- kernel/include/IMP/decorators/macros.h	(revision 626)
+++ kernel/include/IMP/decorators/macros.h	(working copy)
@@ -261,5 +261,15 @@
 #define IMP_DECORATOR_ARRAY_INIT(DecoratorType, name)   \
   name##_data_.initialize();
 
+//! add a method to get a key
+/** One has to make sure to call the
+    decorator_initialize_static_data method first
+ */
+#define IMP_DECORATOR_GET_KEY(KeyType, key_name, variable_name)\
+  static KeyType get_##key_name() {                            \
+  decorator_initialize_static_data();                          \
+  return variable_name;                                        \
+  }
 
+
 #endif  /* __IMP_DECORATOR_MACROS_H */
Index: kernel/include/IMP/OptimizerState.h
===================================================================
--- kernel/include/IMP/OptimizerState.h	(revision 626)
+++ kernel/include/IMP/OptimizerState.h	(working copy)
@@ -60,11 +60,11 @@
   Optimizer *get_optimizer() const {
     IMP_assert(optimizer_,
                "Must call set_optimizer before get_optimizer on state");
-    return optimizer_.get();
+    return optimizer_;
   }
 protected:
-  //! Stored optimizer
-  Pointer<Optimizer> optimizer_;
+  //! Stored optimizer. Non-ref counted
+  Optimizer *optimizer_;
 };
 
 IMP_OUTPUT_OPERATOR(OptimizerState);
Index: kernel/include/IMP/pair_scores/SphereDistancePairScore.h
===================================================================
--- kernel/include/IMP/pair_scores/SphereDistancePairScore.h	(revision 626)
+++ kernel/include/IMP/pair_scores/SphereDistancePairScore.h	(working copy)
@@ -11,6 +11,7 @@
 #include "../PairScore.h"
 #include "../Pointer.h"
 #include "../UnaryFunction.h"
+#include "../decorators/XYZRDecorator.h"
 
 namespace IMP
 {
@@ -24,10 +25,10 @@
   FloatKey radius_;
 public:
   SphereDistancePairScore(UnaryFunction *f, 
-                          FloatKey radius=FloatKey("radius"));
+                          FloatKey radius= XYZRDecorator::get_radius_key());
   virtual ~SphereDistancePairScore(){}
   virtual Float evaluate(Particle *a, Particle *b,
-                         DerivativeAccumulator *da);
+                         DerivativeAccumulator *da) const;
   virtual void show(std::ostream &out=std::cout) const;
 };
 
Index: kernel/include/IMP/pair_scores/TypedPairScore.h
===================================================================
--- kernel/include/IMP/pair_scores/TypedPairScore.h	(revision 626)
+++ kernel/include/IMP/pair_scores/TypedPairScore.h	(working copy)
@@ -38,7 +38,7 @@
   virtual ~TypedPairScore() {}
 
   virtual Float evaluate(Particle *a, Particle *b,
-                         DerivativeAccumulator *da);
+                         DerivativeAccumulator *da) const;
 
   virtual void show(std::ostream &out=std::cout) const;
 
@@ -48,7 +48,7 @@
       be overridden in a subclass to automatically set the type of a particle,
       e.g. from other particle attributes such as an atom or residue name.
    */
-  virtual void set_particle_type(Particle *p) {}
+  virtual void set_particle_type(Particle *p) const {}
 
   //! Set the PairScore to delegate to for a given pair of particle types.
   /** \param[in] ps PairScore to use at evaluate time.
Index: kernel/include/IMP/pair_scores/SConscript
===================================================================
--- kernel/include/IMP/pair_scores/SConscript	(revision 626)
+++ kernel/include/IMP/pair_scores/SConscript	(working copy)
@@ -1,11 +1,12 @@
+Import('env')
 import os.path
-Import('env')
-
-files = ['DistancePairScore.h', 'SphereDistancePairScore.h',
-         'RefineOncePairScore.h', 'TypedPairScore.h',
-         'TransformedDistancePairScore.h']
-
-# Install the include files:
-includedir = os.path.join(env['includedir'], 'IMP', 'pair_scores')
+files=[
+       'DistancePairScore.h',
+       'RefineOncePairScore.h',
+       'SphereDistancePairScore.h',
+       'TransformedDistancePairScore.h',
+       'TypedPairScore.h',
+      ]
+includedir = os.path.join(env['includedir'], 'IMP', 'pair_scores' )
 inst = env.Install(includedir, files)
 env.Alias('install', inst)
Index: kernel/include/IMP/pair_scores/DistancePairScore.h
===================================================================
--- kernel/include/IMP/pair_scores/DistancePairScore.h	(revision 626)
+++ kernel/include/IMP/pair_scores/DistancePairScore.h	(working copy)
@@ -25,8 +25,11 @@
   DistancePairScore(UnaryFunction *f);
   virtual ~DistancePairScore(){}
   virtual Float evaluate(Particle *a, Particle *b,
-                         DerivativeAccumulator *da);
+                         DerivativeAccumulator *da) const;
   virtual void show(std::ostream &out=std::cout) const;
+  UnaryFunction *get_unary_function() const {
+    return f_;
+  }
 };
 
 } // namespace IMP
Index: kernel/include/IMP/pair_scores/RefineOncePairScore.h
===================================================================
--- kernel/include/IMP/pair_scores/RefineOncePairScore.h	(revision 626)
+++ kernel/include/IMP/pair_scores/RefineOncePairScore.h	(working copy)
@@ -34,7 +34,7 @@
   RefineOncePairScore(ParticleRefiner *r, PairScore *f);
   virtual ~RefineOncePairScore(){}
   virtual Float evaluate(Particle *a, Particle *b,
-                         DerivativeAccumulator *da);
+                         DerivativeAccumulator *da) const;
   virtual void show(std::ostream &out=std::cout) const;
 };
 
Index: kernel/include/IMP/pair_scores/TransformedDistancePairScore.h
===================================================================
--- kernel/include/IMP/pair_scores/TransformedDistancePairScore.h	(revision 626)
+++ kernel/include/IMP/pair_scores/TransformedDistancePairScore.h	(working copy)
@@ -25,6 +25,11 @@
 
     The second particle, x, is transformed as R*(x-center)+ translation+center
 
+    \note While it would be nice to have this apply a general pair
+    score to the transformed particle, this would require either creating
+    a temporary particle or rewriting and then restoring the coordinates
+    of the particle being transformed. Neither of which sound appealing.
+
     \ingroup pairscore
  */
 class IMPDLLEXPORT TransformedDistancePairScore : public PairScore
@@ -36,7 +41,7 @@
   TransformedDistancePairScore(UnaryFunction *f);
   virtual ~TransformedDistancePairScore(){}
   virtual Float evaluate(Particle *a, Particle *b,
-                         DerivativeAccumulator *da);
+                         DerivativeAccumulator *da) const;
   virtual void show(std::ostream &out=std::cout) const;
 
   void set_rotation(float r00, float r01, float r02,
Index: kernel/include/IMP/internal/SConscript
===================================================================
--- kernel/include/IMP/internal/SConscript	(revision 626)
+++ kernel/include/IMP/internal/SConscript	(working copy)
@@ -1,14 +1,23 @@
+Import('env')
 import os.path
-Import('env')
-
-files = ['AttributeTable.h',  'graph_base.h',  'Grid3D.h', 'Vector.h',
-         'ref_counting.h', 'ObjectContainer.h', 'ParticleGrid.h',
-         'kernel_version_info.h', 'constants.h', 'units.h',
-         'utility.h', 'bbox_nbl_helpers.h',
-         'ArrayOnAttributesHelper.h', 'Unit.h', 'ExponentialNumber.h',
-         'evaluate_distance_pair_score.h']
-
-# Install the include files:
-includedir = os.path.join(env['includedir'], 'IMP', 'internal')
+files=[
+       'ArrayOnAttributesHelper.h',
+       'AttributeTable.h',
+       'bbox_nbl_helpers.h',
+       'constants.h',
+       'evaluate_distance_pair_score.h',
+       'ExponentialNumber.h',
+       'graph_base.h',
+       'Grid3D.h',
+       'kernel_version_info.h',
+       'ObjectContainer.h',
+       'ParticleGrid.h',
+       'ref_counting.h',
+       'Unit.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/units.h
===================================================================
--- kernel/include/IMP/internal/units.h	(revision 626)
+++ kernel/include/IMP/internal/units.h	(working copy)
@@ -119,6 +119,12 @@
 typedef Multiply<Multiply<Centimeter, Centimeter>::type,
                       Centimeter>::type CubicCentimeter;
 typedef Divide<Gram, CubicCentimeter>::type GramPerCubicCentimeter;
+typedef Shift<Second, -9>::type Nanosecond;
+typedef Shift<Meter, -9>::type Nanometer;
+typedef Multiply<Multiply<Nanometer, Nanometer>::type,
+                      Nanometer>::type CubicNanometer;
+typedef Shift<Second, -15>::type Femtosecond;
+typedef Divide<Piconewton, Nanometer>::type PiconewtonPerNanometer;
 
 
 
Index: kernel/include/IMP/internal/AttributeTable.h
===================================================================
--- kernel/include/IMP/internal/AttributeTable.h	(revision 626)
+++ kernel/include/IMP/internal/AttributeTable.h	(working copy)
@@ -183,6 +183,11 @@
     map_[k.get_index()] = v;
   }
 
+  void set_values(Value v) {
+    for (unsigned int i=0; i< size_; ++i) {
+      map_[i]=v;
+    }
+  }
 
   void insert(Key k, Value v) {
     IMP_assert(!contains(k),
Index: kernel/include/IMP/internal/Grid3D.h
===================================================================
--- kernel/include/IMP/internal/Grid3D.h	(revision 626)
+++ kernel/include/IMP/internal/Grid3D.h	(working copy)
@@ -139,6 +139,16 @@
   }
 };
 
+
+
+
+
+
+
+
+
+
+
 //! Represent a real cell in a grid
 /** These indexes represent an actual cell in the grid with no mapping.
     They can only be constructed by the grid.
@@ -160,6 +170,16 @@
   }
 };
 
+
+
+
+
+
+
+
+
+
+
 /** \brief A voxel grid in 3D space.
     VT can be any class.
     \internal
@@ -340,12 +360,12 @@
   }
 
   //! Get the data in a particular cell
-  VoxelData& get_voxel(Index gi) {
+  typename std::vector<VT>::reference get_voxel(Index gi) {
     return data_[index(gi)];
   }
 
   //! Get the data in a particular cell
-  const VoxelData& get_voxel(Index gi) const  {
+  typename std::vector<VT>::const_reference get_voxel(Index gi) const  {
     return data_[index(gi)];
   }
 
Index: kernel/include/IMP/internal/ref_counting.h
===================================================================
--- kernel/include/IMP/internal/ref_counting.h	(revision 626)
+++ kernel/include/IMP/internal/ref_counting.h	(working copy)
@@ -27,7 +27,6 @@
   template <class O>
   static void eval(O* o) {
     BOOST_STATIC_ASSERT((!boost::is_base_of<RefCountedObject, O >::value));
-    IMP_LOG(VERBOSE, "Not refing particle " << o << std::endl);
   }
 };
 
@@ -36,8 +35,6 @@
 {
   template <class O>
   static void eval(O* o) {
-    IMP_LOG(VERBOSE, "Refing particle " << o->get_index() 
-            << o->get_ref_count() << std::endl);
     o->assert_is_valid();
     o->ref();
   }
@@ -49,7 +46,6 @@
   template <class O>
   static void eval(O* o) {
     BOOST_STATIC_ASSERT((!boost::is_base_of<RefCountedObject, O >::value));
-    IMP_LOG(VERBOSE, "Not Unrefing object " << o << std::endl);
   }
 };
 
@@ -58,8 +54,6 @@
 {
   template <class O>
   static void eval(O *o) {
-    IMP_LOG(VERBOSE, "Unrefing particle " << o->get_index()
-            << " " << o->get_ref_count() << std::endl);
     o->assert_is_valid();
     o->unref();
     if (!o->get_has_ref()) {
@@ -104,9 +98,6 @@
 template <class O>
 void disown(O* o)
 {
-  /*IMP_LOG(VERBOSE, "Disown called with " 
-          << (boost::is_base_of<RefCountedObject, O >::value)
-          << " for " << o << " " << o->get_ref_count() << std::endl);*/
   o->unref();
   if (!o->get_has_ref()) {
     delete o;
@@ -118,10 +109,6 @@
 template <class O>
 void own(O* o)
 {
-  /*IMP_LOG(VERBOSE, "Own called with "
-          << (boost::is_base_of<RefCountedObject, O >::value)
-          << " for " << o
-          << " " << o->get_ref_count() << std::endl);*/
   if (boost::is_base_of<RefCountedObject, O >::value) {
     // no checks
   } else {
@@ -133,6 +120,9 @@
 }
 
 
+
+
+
 } // namespace internal
 
 } // namespace IMP
Index: kernel/include/IMP/PairScore.h
===================================================================
--- kernel/include/IMP/PairScore.h	(revision 626)
+++ kernel/include/IMP/PairScore.h	(working copy)
@@ -32,7 +32,7 @@
   virtual ~PairScore() {}
   //! Compute the score for the pair and the derivative if needed.
   virtual Float evaluate(Particle *a, Particle *b,
-                         DerivativeAccumulator *da) = 0;
+                         DerivativeAccumulator *da) const = 0;
   virtual void show(std::ostream &out=std::cout) const = 0;
 };
 
Index: kernel/include/IMP/Key.h
===================================================================
--- kernel/include/IMP/Key.h	(revision 626)
+++ kernel/include/IMP/Key.h	(working copy)
@@ -70,7 +70,7 @@
   typedef std::vector<Name> Name##s
 
 
-/** This must occur in exactly one .o in the internal namespace. Should 
+/** This must occur in exactly one .o in the application. Should 
  be used in the IMP namespace.*/
 #define IMP_DEFINE_KEY_TYPE(Name, Tag)                  \
   namespace internal {                                  \
Index: kernel/include/IMP/Vector3D.h
===================================================================
--- kernel/include/IMP/Vector3D.h	(revision 626)
+++ kernel/include/IMP/Vector3D.h	(working copy)
@@ -120,6 +120,27 @@
                     operator[](2) + o[2]);
   }
 
+  //! Accumulate the vector
+  void operator+=(const Vector3D &o) {
+    vec_[0]+= o[0];
+    vec_[1]+= o[1];
+    vec_[2]+= o[2];
+  }
+
+  //! rescale the vector
+  void operator/=(float f) {
+    vec_[0]/= f;
+    vec_[1]/= f;
+    vec_[2]/= f;
+  }
+
+  //! rescale the vector
+  void operator*=(float f) {
+    vec_[0]*= f;
+    vec_[1]*= f;
+    vec_[2]*= f;
+  }
+
   void show(std::ostream &out=std::cout) const {
     out << "(" << operator[](0) << ", " << operator[](1) << ", "
         << operator[](2) << ")";
Index: kernel/include/IMP/exception.h
===================================================================
--- kernel/include/IMP/exception.h	(revision 626)
+++ kernel/include/IMP/exception.h	(working copy)
@@ -20,6 +20,7 @@
 
 namespace IMP
 {
+
 //! 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.
@@ -65,6 +66,7 @@
   }
 };
 
+
 //! A general exception for an error in IMP.
 /** \ingroup assert
  */
@@ -147,6 +149,10 @@
  */
 IMPDLLEXPORT void check_fail(const char *msg);
 
+
+//! Control this through functions in Log.h
+extern bool print_exceptions;
+
 } // namespace internal
 
 } // namespace IMP
@@ -199,6 +205,10 @@
     \param[in] exception Throw the object constructed by this expression.
     \ingroup assert
  */
-#define IMP_failure(message, exception) {IMP_ERROR(message); throw exception;}
+#define IMP_failure(message, exception) { \
+    std::ostringstream oss;                                             \
+    oss << message << std::endl;                                        \
+    IMP::internal::check_fail(oss.str().c_str());                       \
+    throw exception(oss.str().c_str());}
 
 #endif  /* __IMP_EXCEPTION_H */
Index: kernel/include/IMP/unary_functions/ClosedCubicSpline.h
===================================================================
--- kernel/include/IMP/unary_functions/ClosedCubicSpline.h	(revision 626)
+++ kernel/include/IMP/unary_functions/ClosedCubicSpline.h	(working copy)
@@ -36,7 +36,7 @@
       \exception ValueException Feature is out of defined range.
       \return Score
    */
-  virtual Float evaluate(Float feature);
+  virtual Float evaluate(Float feature) const;
 
   //! Calculate score and derivative with respect to the given feature.
   /** \param[in] feature Value of feature being tested.
@@ -45,7 +45,7 @@
       \exception ValueException Feature is out of defined range.
       \return Score
    */
-  virtual Float evaluate_deriv(Float feature, Float& deriv);
+  virtual Float evaluate_deriv(Float feature, Float& deriv) const;
 
   void show(std::ostream &out=std::cout) const {
     out << "Closed cubic spline of " << values_.size() << " values from "
Index: kernel/include/IMP/unary_functions/Linear.h
===================================================================
--- kernel/include/IMP/unary_functions/Linear.h	(revision 626)
+++ kernel/include/IMP/unary_functions/Linear.h	(working copy)
@@ -20,13 +20,15 @@
 class Linear : public UnaryFunction
 {
 public:
-  Linear(Float slope, Float offset=0) : slope_(slope), offset_(offset) {}
+  Linear(Float offset, Float slope) : slope_(slope), offset_(offset) {}
 
   virtual ~Linear() {}
 
-  virtual Float evaluate(Float feature) { return (feature-offset_)*slope_; }
+  virtual Float evaluate(Float feature) const {
+    return (feature-offset_)*slope_;
+  }
 
-  virtual Float evaluate_deriv(Float feature, Float& deriv) {
+  virtual Float evaluate_deriv(Float feature, Float& deriv) const {
     deriv= slope_;
     return evaluate(feature);
   }
Index: kernel/include/IMP/unary_functions/WormLikeChain.h
===================================================================
--- kernel/include/IMP/unary_functions/WormLikeChain.h	(revision 626)
+++ kernel/include/IMP/unary_functions/WormLikeChain.h	(working copy)
@@ -43,7 +43,7 @@
   /** \param[in] l Current length in Angstroms
       \return Energy in kcal/mol
    */
-  virtual Float evaluate(Float lf) {
+  virtual Float evaluate(Float lf) const {
     static const unit::Picojoule zero=eval(unit::Angstrom(0));
     unit::Angstrom l(lf);
     if (l < unit::Angstrom(0)) l=unit::Angstrom(0);
@@ -65,7 +65,7 @@
       \param[out] deriv force in kcal/angstrom mol
       \return Score
    */
-  virtual Float evaluate_deriv(Float fl, Float& deriv) {
+  virtual Float evaluate_deriv(Float fl, Float& deriv) const {
     unit::Angstrom l(fl);
     if (l < unit::Angstrom(0)) l=unit::Angstrom(0);
     unit::Piconewton doubled;
Index: kernel/include/IMP/unary_functions/Cosine.h
===================================================================
--- kernel/include/IMP/unary_functions/Cosine.h	(revision 626)
+++ kernel/include/IMP/unary_functions/Cosine.h	(working copy)
@@ -39,7 +39,7 @@
   /** \param[in] feature Value of feature being tested.
       \return Score
    */
-  virtual Float evaluate(Float feature);
+  virtual Float evaluate(Float feature) const;
 
   //! Calculate score and derivative with respect to the given feature.
   /** \param[in] feature Value of feature being tested.
@@ -47,7 +47,7 @@
                         the feature value.
       \return Score
    */
-  virtual Float evaluate_deriv(Float feature, Float& deriv);
+  virtual Float evaluate_deriv(Float feature, Float& deriv) const;
 
   void show(std::ostream &out=std::cout) const {
     out << "Cosine function with force " << force_constant_
Index: kernel/include/IMP/unary_functions/Harmonic.h
===================================================================
--- kernel/include/IMP/unary_functions/Harmonic.h	(revision 626)
+++ kernel/include/IMP/unary_functions/Harmonic.h	(working copy)
@@ -55,7 +55,7 @@
   /** \param[in] feature Value of feature being tested.
       \return Score
    */
-  virtual Float evaluate(Float feature) {
+  virtual Float evaluate(Float feature) const {
     Float d;
     return evaluate_deriv(feature, d);
   }
@@ -66,7 +66,7 @@
                         the feature value.
       \return Score
    */
-  virtual Float evaluate_deriv(Float feature, Float& deriv) {
+  virtual Float evaluate_deriv(Float feature, Float& deriv) const {
     Float e = (feature - mean_);
     deriv = k_ * e;
     return 0.5 * k_ * e * e;
Index: kernel/include/IMP/unary_functions/HarmonicLowerBound.h
===================================================================
--- kernel/include/IMP/unary_functions/HarmonicLowerBound.h	(revision 626)
+++ kernel/include/IMP/unary_functions/HarmonicLowerBound.h	(working copy)
@@ -26,7 +26,7 @@
       \param[in] feature Value of feature being tested.
       \return Score
   */
-  virtual Float evaluate(Float feature) {
+  virtual Float evaluate(Float feature) const {
     if (feature >= Harmonic::get_mean()) {
       return 0.0;
     } else {
@@ -41,7 +41,7 @@
                         the feature value.
       \return Score
    */
-  virtual Float evaluate_deriv(Float feature, Float& deriv) {
+  virtual Float evaluate_deriv(Float feature, Float& deriv) const {
     if (feature >= Harmonic::get_mean()) {
       deriv = 0.0;
       return 0.0;
Index: kernel/include/IMP/unary_functions/SConscript
===================================================================
--- kernel/include/IMP/unary_functions/SConscript	(revision 626)
+++ kernel/include/IMP/unary_functions/SConscript	(working copy)
@@ -1,11 +1,15 @@
+Import('env')
 import os.path
-Import('env')
-
-files = ['Harmonic.h', 'HarmonicLowerBound.h', 'HarmonicUpperBound.h',
-         'OpenCubicSpline.h', 'ClosedCubicSpline.h', 'Cosine.h', 'Linear.h',
-         'WormLikeChain.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/unary_functions/OpenCubicSpline.h
===================================================================
--- kernel/include/IMP/unary_functions/OpenCubicSpline.h	(revision 626)
+++ kernel/include/IMP/unary_functions/OpenCubicSpline.h	(working copy)
@@ -25,7 +25,7 @@
       \param[in] minrange Feature value at first spline point
       \param[in] spacing  Distance (in feature space) between points
    */
-  OpenCubicSpline(const std::vector<Float> &values, Float minrange,
+  OpenCubicSpline(const Floats &values, Float minrange,
                   Float spacing);
 
   virtual ~OpenCubicSpline() {}
@@ -35,7 +35,7 @@
       \exception ValueException Feature is out of defined range.
       \return Score
    */
-  virtual Float evaluate(Float feature);
+  virtual Float evaluate(Float feature) const;
 
   //! Calculate score and derivative with respect to the given feature.
   /** \param[in] feature Value of feature being tested.
@@ -44,7 +44,7 @@
       \exception ValueException Feature is out of defined range.
       \return Score
    */
-  virtual Float evaluate_deriv(Float feature, Float& deriv);
+  virtual Float evaluate_deriv(Float feature, Float& deriv) const ;
 
   void show(std::ostream &out=std::cout) const {
     out << "Open cubic spline of " << values_.size() << " values from "
Index: kernel/include/IMP/unary_functions/HarmonicUpperBound.h
===================================================================
--- kernel/include/IMP/unary_functions/HarmonicUpperBound.h	(revision 626)
+++ kernel/include/IMP/unary_functions/HarmonicUpperBound.h	(working copy)
@@ -26,7 +26,7 @@
       \param[in] feature Value of feature being tested.
       \return Score
    */
-  virtual Float evaluate(Float feature) {
+  virtual Float evaluate(Float feature) const {
     if (feature <= Harmonic::get_mean()) {
       return 0.0;
     } else {
@@ -41,7 +41,7 @@
                         the feature value.
       \return Score
    */
-  virtual Float evaluate_deriv(Float feature, Float& deriv) {
+  virtual Float evaluate_deriv(Float feature, Float& deriv) const {
     if (feature <= Harmonic::get_mean()) {
       deriv = 0.0;
       return 0.0;
Index: kernel/include/IMP/SConscript
===================================================================
--- kernel/include/IMP/SConscript	(revision 626)
+++ kernel/include/IMP/SConscript	(working copy)
@@ -1,28 +1,44 @@
+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',
-         'Object.h', 'Pointer.h', 'RefCountedObject.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',
+       'Object.h',
+       'Optimizer.h',
+       'OptimizerState.h',
+       'PairScore.h',
+       'Particle.h',
+       'ParticleRefiner.h',
+       'Pointer.h',
+       'random.h',
+       'RefCountedObject.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('particle_refiners/SConscript')
+SConscript('decorators/SConscript' )
+SConscript('internal/SConscript' )
+SConscript('optimizers/SConscript' )
+SConscript('pair_scores/SConscript' )
+SConscript('particle_refiners/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 626)
+++ kernel/include/IMP/Particle.h	(working copy)
@@ -98,7 +98,7 @@
   /** \return all particle data in the model.
    */
   Model* get_model() const {
-    return model_.get();
+    return model_;
   }
 
   //! Add a Float attribute to this particle.
@@ -369,7 +369,8 @@
   // Set pointer to model particle data.
   void set_model(Model *md, ParticleIndex pi);
 
-  Pointer<Model> model_;
+  // not ref counted as the model has a pointer to this
+  Model *model_;
 
   // true if particle is active
   bool is_active_;
Index: kernel/include/IMP/particle_refiners/BondCoverParticleRefiner.h
===================================================================
--- kernel/include/IMP/particle_refiners/BondCoverParticleRefiner.h	(revision 626)
+++ kernel/include/IMP/particle_refiners/BondCoverParticleRefiner.h	(working copy)
@@ -8,12 +8,14 @@
 #ifndef __IMP_BOND_COVER_PARTICLE_REFINER_H
 #define __IMP_BOND_COVER_PARTICLE_REFINER_H
 
+#include "../ParticleRefiner.h"
 #include "../internal/kernel_version_info.h"
-#include "../ParticleRefiner.h"
 
 namespace IMP
 {
 
+class Particle;
+
 //! Cover a bond with a constant volume set of spheres.
 /** Perhaps I want to add various custom bond types so that 
     this will only expand some custom bonds. Currently any
@@ -23,9 +25,11 @@
 {
   FloatKey rk_;
   FloatKey vk_;
+  IntKey tk_;
 public:
   BondCoverParticleRefiner(FloatKey rk,
-                           FloatKey vk);
+                           FloatKey vk,
+                           IntKey tk=IntKey());
 
   virtual ~BondCoverParticleRefiner() {}
 
Index: kernel/include/IMP/particle_refiners/SConscript
===================================================================
--- kernel/include/IMP/particle_refiners/SConscript	(revision 626)
+++ kernel/include/IMP/particle_refiners/SConscript	(working copy)
@@ -1,9 +1,9 @@
 Import('env')
 import os.path
-
-files=['BondCoverParticleRefiner.h', 'ChildrenParticleRefiner.h']
-
-# Install the include files
-includedir = os.path.join(env['includedir'], 'IMP', 'particle_refiners')
+files=[
+       'BondCoverParticleRefiner.h',
+       'ChildrenParticleRefiner.h',
+      ]
+includedir = os.path.join(env['includedir'], 'IMP', 'particle_refiners' )
 inst = env.Install(includedir, files)
 env.Alias('install', inst)
Index: kernel/include/IMP/triplet_scores/SConscript
===================================================================
--- kernel/include/IMP/triplet_scores/SConscript	(revision 626)
+++ 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/triplet_scores/AngleTripletScore.h
===================================================================
--- kernel/include/IMP/triplet_scores/AngleTripletScore.h	(revision 626)
+++ kernel/include/IMP/triplet_scores/AngleTripletScore.h	(working copy)
@@ -26,7 +26,7 @@
   AngleTripletScore(UnaryFunction *f);
   virtual ~AngleTripletScore(){}
   virtual Float evaluate(Particle *a, Particle *b, Particle *c,
-                         DerivativeAccumulator *da);
+                         DerivativeAccumulator *da) const;
   virtual void show(std::ostream &out=std::cout) const;
 };
 
Index: kernel/include/IMP/ScoreState.h
===================================================================
--- kernel/include/IMP/ScoreState.h	(revision 626)
+++ kernel/include/IMP/ScoreState.h	(working copy)
@@ -82,7 +82,7 @@
   Model *get_model() const {
     IMP_assert(model_,
                "Must call set_model before get_model on state");
-    return model_.get();
+    return model_;
   }
 
 protected:
@@ -120,8 +120,8 @@
 
   unsigned int update_iteration_;
   unsigned int after_iteration_;
-  // all of the particle data
-  Pointer<Model> model_;
+  // not ref counted or destructed
+  Model *model_;
   std::string name_;
 };
 
Index: kernel/include/IMP/TripletScore.h
===================================================================
--- kernel/include/IMP/TripletScore.h	(revision 626)
+++ kernel/include/IMP/TripletScore.h	(working copy)
@@ -33,7 +33,7 @@
   virtual ~TripletScore() {}
   //! Compute the score for the triplet and the derivative if needed.
   virtual Float evaluate(Particle *a, Particle *b, Particle *c,
-                         DerivativeAccumulator *da) = 0;
+                         DerivativeAccumulator *da) const = 0;
   virtual void show(std::ostream &out=std::cout) const = 0;
 };
 
Index: kernel/include/IMP/Pointer.h
===================================================================
--- kernel/include/IMP/Pointer.h	(revision 626)
+++ kernel/include/IMP/Pointer.h	(working copy)
@@ -103,8 +103,8 @@
     return !o_;
   }
 
-  operator unspecified_bool() const {
-    return o_ ? &This::operator! : 0;
+  operator O*() const {
+    return o_;
   }
 };
 
Index: kernel/include/IMP/log.h
===================================================================
--- kernel/include/IMP/log.h	(revision 626)
+++ kernel/include/IMP/log.h	(working copy)
@@ -183,6 +183,20 @@
 
 
 
+//! Print error messages
+/** If this is true, then failures of the IMP_check
+    macro will print their messages to std::cerr. This
+    should be true if the code is used from C++ and false
+    if used from python.
+ */
+IMPDLLEXPORT bool get_print_exception_messages();
+
+
+//! Turn on and off printing of error messages
+IMPDLLEXPORT void set_print_exception_messages(bool tf);
+
+
+
 } // namespace IMP
 
 #ifndef IMP_DISABLE_LOGGING
Index: kernel/include/IMP/macros.h
===================================================================
--- kernel/include/IMP/macros.h	(revision 626)
+++ kernel/include/IMP/macros.h	(working copy)
@@ -226,9 +226,9 @@
   virtual void show(std::ostream &out) const;                           \
   /** Destroy any created particles and propagate derivatives */        \
   virtual void cleanup_refined(Particle *a, Particles &b,               \
-                       DerivativeAccumulator *da=0);                    \
+                               DerivativeAccumulator *da=0) const;      \
   /** Return a list of particles which refines the passed particle.*/   \
-  virtual Particles get_refined(Particle *);                            \
+  virtual Particles get_refined(Particle *) const;                      \
   virtual IMP::VersionInfo get_version_info() const { return version_info; }
 
 //! Use the swap_with member function to swap two objects
Index: kernel/include/SConscript
===================================================================
--- kernel/include/SConscript	(revision 626)
+++ 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 626)
+++ kernel/include/IMP.h	(working copy)
@@ -1,94 +1,116 @@
 /**
- *  \file IMP.h   \brief IMP, an Integrative Modeling Platform.
- *
- *  Copyright 2007-8 Sali Lab. All rights reserved.
- *
- */
-
+*  ile IMP.h   rief 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/Object.h"
-#include "IMP/RefCountedObject.h"
-#include "IMP/Pointer.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/unary_functions/WormLikeChain.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/ParticleRefiner.h"
-#include "IMP/particle_refiners/BondCoverParticleRefiner.h"
-#include "IMP/particle_refiners/ChildrenParticleRefiner.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/BrownianDynamics.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/pair_scores/RefineOncePairScore.h"
-#include "IMP/pair_scores/TypedPairScore.h"
-#include "IMP/pair_scores/TransformedDistancePairScore.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/ConstantRestraint.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/restraints/TunnelRestraint.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/GravityCenterScoreState.h"
-#include "IMP/score_states/CoverBondsScoreState.h"
-
-
-/**
-   \namespace IMP The IMP namespace.
- */
-
+#include <IMP/base_types.h>
+#include <IMP/DecoratorBase.h>
+#include <IMP/DerivativeAccumulator.h>
+#include <IMP/exception.h>
+#include <IMP/IMP_config.h>
+#include <IMP/Index.h>
+#include <IMP/Key.h>
+#include <IMP/log.h>
+#include <IMP/macros.h>
+#include <IMP/Model.h>
+#include <IMP/Object.h>
+#include <IMP/Optimizer.h>
+#include <IMP/OptimizerState.h>
+#include <IMP/PairScore.h>
+#include <IMP/Particle.h>
+#include <IMP/ParticleRefiner.h>
+#include <IMP/Pointer.h>
+#include <IMP/random.h>
+#include <IMP/RefCountedObject.h>
+#include <IMP/Restraint.h>
+#include <IMP/ScoreState.h>
+#include <IMP/SingletonScore.h>
+#include <IMP/TripletScore.h>
+#include <IMP/UnaryFunction.h>
+#include <IMP/utility.h>
+#include <IMP/Vector3D.h>
+#include <IMP/VersionInfo.h>
+#include <IMP/decorators/AtomDecorator.h>
+#include <IMP/decorators/bond_decorators.h>
+#include <IMP/decorators/HierarchyDecorator.h>
+#include <IMP/decorators/macros.h>
+#include <IMP/decorators/MolecularHierarchyDecorator.h>
+#include <IMP/decorators/NameDecorator.h>
+#include <IMP/decorators/ResidueDecorator.h>
+#include <IMP/decorators/utility.h>
+#include <IMP/decorators/XYZDecorator.h>
+#include <IMP/decorators/XYZRDecorator.h>
+#include <IMP/decorators/yaml.h>
+#include <IMP/internal/ArrayOnAttributesHelper.h>
+#include <IMP/internal/AttributeTable.h>
+#include <IMP/internal/bbox_nbl_helpers.h>
+#include <IMP/internal/constants.h>
+#include <IMP/internal/evaluate_distance_pair_score.h>
+#include <IMP/internal/ExponentialNumber.h>
+#include <IMP/internal/graph_base.h>
+#include <IMP/internal/Grid3D.h>
+#include <IMP/internal/kernel_version_info.h>
+#include <IMP/internal/ObjectContainer.h>
+#include <IMP/internal/ParticleGrid.h>
+#include <IMP/internal/ref_counting.h>
+#include <IMP/internal/Unit.h>
+#include <IMP/internal/units.h>
+#include <IMP/internal/utility.h>
+#include <IMP/internal/Vector.h>
+#include <IMP/optimizers/BrownianDynamics.h>
+#include <IMP/optimizers/ConjugateGradients.h>
+#include <IMP/optimizers/MolecularDynamics.h>
+#include <IMP/optimizers/MonteCarlo.h>
+#include <IMP/optimizers/MoverBase.h>
+#include <IMP/optimizers/Mover.h>
+#include <IMP/optimizers/SteepestDescent.h>
+#include <IMP/pair_scores/DistancePairScore.h>
+#include <IMP/pair_scores/RefineOncePairScore.h>
+#include <IMP/pair_scores/SphereDistancePairScore.h>
+#include <IMP/pair_scores/TransformedDistancePairScore.h>
+#include <IMP/pair_scores/TypedPairScore.h>
+#include <IMP/particle_refiners/BondCoverParticleRefiner.h>
+#include <IMP/particle_refiners/ChildrenParticleRefiner.h>
+#include <IMP/restraints/AngleRestraint.h>
+#include <IMP/restraints/BondDecoratorRestraint.h>
+#include <IMP/restraints/ConnectivityRestraint.h>
+#include <IMP/restraints/ConstantRestraint.h>
+#include <IMP/restraints/DihedralRestraint.h>
+#include <IMP/restraints/DistanceRestraint.h>
+#include <IMP/restraints/NonbondedRestraint.h>
+#include <IMP/restraints/PairChainRestraint.h>
+#include <IMP/restraints/PairListRestraint.h>
+#include <IMP/restraints/PairRestraint.h>
+#include <IMP/restraints/RestraintSet.h>
+#include <IMP/restraints/SingletonListRestraint.h>
+#include <IMP/restraints/TripletChainRestraint.h>
+#include <IMP/score_states/AllNonbondedListScoreState.h>
+#include <IMP/score_states/BipartiteNonbondedListScoreState.h>
+#include <IMP/score_states/BondDecoratorListScoreState.h>
+#include <IMP/score_states/BondedListScoreState.h>
+#include <IMP/score_states/CoverBondsScoreState.h>
+#include <IMP/score_states/GravityCenterScoreState.h>
+#include <IMP/score_states/ManualBondDecoratorListScoreState.h>
+#include <IMP/score_states/MaxChangeScoreState.h>
+#include <IMP/score_states/NonbondedListScoreState.h>
+#include <IMP/singleton_scores/AttributeSingletonScore.h>
+#include <IMP/singleton_scores/DistanceToSingletonScore.h>
+#include <IMP/singleton_scores/TunnelSingletonScore.h>
+#include <IMP/triplet_scores/AngleTripletScore.h>
+#include <IMP/unary_functions/ClosedCubicSpline.h>
+#include <IMP/unary_functions/Cosine.h>
+#include <IMP/unary_functions/Harmonic.h>
+#include <IMP/unary_functions/HarmonicLowerBound.h>
+#include <IMP/unary_functions/HarmonicUpperBound.h>
+#include <IMP/unary_functions/Linear.h>
+#include <IMP/unary_functions/OpenCubicSpline.h>
+#include <IMP/unary_functions/WormLikeChain.h>
+#include <IMP/optimizers/movers/BallMover.h>
+#include <IMP/optimizers/movers/NormalMover.h>
+#include <IMP/optimizers/states/CMMLogOptimizerState.h>
+#include <IMP/optimizers/states/VelocityScalingOptimizerState.h>
+#include <IMP/optimizers/states/VRMLLogOptimizerState.h>
 #endif  /* __IMP_H */
Index: kernel/doc/examples/chain.py
===================================================================
--- kernel/doc/examples/chain.py	(revision 626)
+++ kernel/doc/examples/chain.py	(working copy)
@@ -9,18 +9,17 @@
 #IMP.set_log_level(IMP.VERBOSE)
 np=20
 radius =1.0
-rk= IMP.FloatKey("radius")
 m= IMP.Model()
 # The particles in the chain
 chain= IMP.Particles()
 for i in range(0,np):
     p= IMP.Particle()
     pi= m.add_particle(p)
-    d= IMP.XYZDecorator.create(p)
+    d= IMP.XYZRDecorator.create(p)
     d.randomize_in_box(IMP.Vector3D(0,0,0),
                        IMP.Vector3D(10,10,10))
     d.set_coordinates_are_optimized(True)
-    p.add_attribute(rk, radius, False)
+    d.set_radius(radius)
     chain.append(p)
 
 # create a bond between successive particles
@@ -36,15 +35,14 @@
     p.show()
 
 # Set up the nonbonded list
-nbl= IMP.AllNonbondedListScoreState(1, rk, chain)
+nbl= IMP.AllNonbondedListScoreState(1, chain)
 nbli= m.add_score_state(nbl)
 # This ScoreState uses the bonds constructed above to restrain
 bl= IMP.BondDecoratorListScoreState(chain)
 bli= nbl.add_bonded_list(bl)
 
 # Set up excluded volume
-ps= IMP.SphereDistancePairScore(IMP.HarmonicLowerBound(0,1),
-                                rk)
+ps= IMP.SphereDistancePairScore(IMP.HarmonicLowerBound(0,1))
 evr= IMP.NonbondedRestraint(ps, nbl)
 evri= m.add_restraint(evr)
 
@@ -63,7 +61,7 @@
 p= IMP.ParticlePair(chain[0], chain[-1])
 pps= IMP.ParticlePairs()
 pps.append(p)
-cr= IMP.PairListRestraint(IMP.SphereDistancePairScore(IMP.Harmonic(3,1), rk),
+cr= IMP.PairListRestraint(IMP.SphereDistancePairScore(IMP.Harmonic(3,1)),
                           pps)
 cri=m.add_restraint(cr)
 
@@ -74,7 +72,6 @@
 # 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_key(rk)
 vrml.update()
 vrml.set_skip_steps(100)
 o.add_optimizer_state(vrml)
Index: kernel/src/exception.cpp
===================================================================
--- kernel/src/exception.cpp	(revision 626)
+++ kernel/src/exception.cpp	(working copy)
@@ -10,6 +10,12 @@
 namespace IMP
 {
 
+namespace internal {
+  // The error message is already in the exception
+  bool print_exceptions=true;
+}
+
+
 static CheckLevel check_mode =
 #ifdef NDEBUG
   CHEAP;
@@ -30,12 +36,9 @@
 namespace internal
 {
 
-// The error message is already in the exception
-bool print_exceptions=false;
-
 void assert_fail(const char *msg)
 {
-  if (print_exceptions) {
+  if (internal::print_exceptions) {
     IMP_ERROR(msg);
   }
   throw ErrorException(msg);
@@ -43,11 +46,14 @@
 
 void check_fail(const char *msg)
 {
-  if (print_exceptions) {
+  if (internal::print_exceptions) {
     IMP_ERROR(msg);
   }
 }
 
 } // namespace internal
 
+
+
+
 } // namespace IMP
Index: kernel/src/singleton_scores/TunnelSingletonScore.cpp
===================================================================
--- kernel/src/singleton_scores/TunnelSingletonScore.cpp	(revision 0)
+++ kernel/src/singleton_scores/TunnelSingletonScore.cpp	(revision 0)
@@ -0,0 +1,103 @@
+/**
+ *  \file TunnelSingletonScore.cpp     \brief Tunnel restraint.
+ *
+ *  Copyright 2007-8 Sali Lab. All rights reserved.
+ *
+ */
+
+#include "IMP/singleton_scores/TunnelSingletonScore.h"
+#include "IMP/decorators/XYZRDecorator.h"
+
+namespace IMP
+{
+
+TunnelSingletonScore::TunnelSingletonScore(UnaryFunction *f) : f_(f)
+{
+  coordinate_=0;
+  center_=Vector3D(0,0,0);
+  height_=0;
+  radius_=0;
+}
+
+
+
+Float TunnelSingletonScore::evaluate(Particle *p,
+                                     DerivativeAccumulator *accum) const
+{
+  Float ret=0;
+  XYZDecorator d(p);
+  Float radius=0;
+  if ( p->has_attribute(XYZRDecorator::get_radius_key())) {
+    radius= p->get_value(XYZRDecorator::get_radius_key());
+  }
+  IMP_LOG(VERBOSE, "Tunnel restraint on particle " << d << std::endl);
+  Float hr= height_+radius;
+  if (d.get_coordinate(coordinate_) < center_[coordinate_] + hr
+      && d.get_coordinate(coordinate_) > center_[coordinate_] - hr) {
+    Float sd=0;
+    for (int i=1; i< 3; ++i) {
+      int oc= (i+coordinate_)%3;
+      sd+= square(d.get_coordinate(oc)- center_[oc]);
+    }
+    sd= std::sqrt(sd);
+    IMP_LOG(VERBOSE, "The distance is " << sd << " and radius "
+            << radius_ << std::endl);
+    if (sd > radius_-radius) {
+      Float rd= sd-radius_;
+      Float hdu= center_[coordinate_] + height_
+        - d.get_coordinate(coordinate_);
+      Float hdd=  d.get_coordinate(coordinate_)
+        + height_- center_[coordinate_];
+      Vector3D deriv(0,0,0);
+      Float score=0;
+      Float deriv_scalar=0;
+      /*! \todo Clean up these tests so I am not dependent on two expressions
+        being the same and evaluating to the same thing */
+      // look below if changed
+      Float dist= -std::min(std::min(rd, hdu), hdd) - radius;
+      if (accum) {
+        score= f_->evaluate_deriv(dist, deriv_scalar);
+      } else {
+        score= f_->evaluate(dist);
+      }
+
+      // kind if evil
+      if (dist== -rd -radius) {
+        Vector3D v= (d.get_vector() - center_).get_unit_vector();
+        for (int i=0; i< 2; ++i) {
+          int oc= (i+coordinate_+1)%3;
+          deriv[oc]= v[oc]*deriv_scalar;
+        }
+      } else {
+        // kind of evil
+        if (dist == -hdu -radius) {
+          deriv_scalar= -deriv_scalar;
+        }
+        deriv[coordinate_]= deriv_scalar;
+      }
+
+      if (accum) {
+        for (unsigned int i=0; i< 3; ++i) {
+          d.add_to_coordinate_derivative(i, deriv[i], *accum);
+        }
+      }
+      ret+= score;
+    } else {
+      IMP_LOG(VERBOSE, "Particle " << p->get_index()
+              << " is in channel" << std::endl);
+    }
+  } else {
+    IMP_LOG(VERBOSE, "Particle " << p->get_index()
+            << " is outside of slab" << std::endl);
+  }
+  return ret;
+}
+
+void TunnelSingletonScore::show(std::ostream& out) const
+{
+  out << "Tunnel score :" << *f_ 
+      << " " << center_ << " " << height_ << " " 
+      << coordinate_ << std::endl;
+}
+
+}  // namespace IMP
Index: kernel/src/singleton_scores/SConscript
===================================================================
--- kernel/src/singleton_scores/SConscript	(revision 626)
+++ kernel/src/singleton_scores/SConscript	(working copy)
@@ -1,6 +1,7 @@
 Import('env')
-
-files = ['DistanceToSingletonScore.cpp', 'AttributeSingletonScore.cpp']
-
-files = [File(x) for x in files]
+files=[
+        File( 'AttributeSingletonScore.cpp' ),
+        File( 'DistanceToSingletonScore.cpp' ),
+        File( 'TunnelSingletonScore.cpp' ),
+      ]
 Return('files')
Index: kernel/src/singleton_scores/DistanceToSingletonScore.cpp
===================================================================
--- kernel/src/singleton_scores/DistanceToSingletonScore.cpp	(revision 626)
+++ kernel/src/singleton_scores/DistanceToSingletonScore.cpp	(working copy)
@@ -27,7 +27,7 @@
                                                                        pt_(v){}
 
 Float DistanceToSingletonScore::evaluate(Particle *b,
-                                         DerivativeAccumulator *da)
+                                         DerivativeAccumulator *da) const
 {
   return internal::evaluate_distance_pair_score(XYZDecorator(b),
                                                 StaticD(pt_), da,
Index: kernel/src/Particle.cpp
===================================================================
--- kernel/src/Particle.cpp	(revision 626)
+++ kernel/src/Particle.cpp	(working copy)
@@ -15,6 +15,7 @@
 
 Particle::Particle()
 {
+  model_=NULL;
   is_active_ = true;
 }
 
@@ -43,9 +44,7 @@
 
 void Particle::zero_derivatives()
 {
-  for (FloatKeyIterator it= float_keys_begin(); it != float_keys_end(); ++it) {
-    derivatives_.set_value(*it, 0);
-  }
+  derivatives_.set_values(0);
 }
 
 
Index: kernel/src/unary_functions/Cosine.cpp
===================================================================
--- kernel/src/unary_functions/Cosine.cpp	(revision 626)
+++ kernel/src/unary_functions/Cosine.cpp	(working copy)
@@ -12,13 +12,13 @@
 namespace IMP
 {
 
-Float Cosine::evaluate(Float feature)
+Float Cosine::evaluate(Float feature) const
 {
   return std::abs(force_constant_)
          - force_constant_ * std::cos(periodicity_ * feature + phase_);
 }
 
-Float Cosine::evaluate_deriv(Float feature, Float& deriv)
+Float Cosine::evaluate_deriv(Float feature, Float& deriv) const
 {
   deriv = force_constant_ * periodicity_
           * std::sin(periodicity_ * feature + phase_);
Index: kernel/src/unary_functions/OpenCubicSpline.cpp
===================================================================
--- kernel/src/unary_functions/OpenCubicSpline.cpp	(revision 626)
+++ kernel/src/unary_functions/OpenCubicSpline.cpp	(working copy)
@@ -43,7 +43,7 @@
 }
 
 
-Float OpenCubicSpline::evaluate(Float feature)
+Float OpenCubicSpline::evaluate(Float feature) const
 {
   // check for feature in range
   if (feature < minrange_ || feature > maxrange_) {
@@ -66,7 +66,7 @@
            * (spacing_ * spacing_) / 6.;
 }
 
-Float OpenCubicSpline::evaluate_deriv(Float feature, Float& deriv)
+Float OpenCubicSpline::evaluate_deriv(Float feature, Float& deriv) const
 {
   size_t lowbin = static_cast<size_t>((feature - minrange_) / spacing_);
   // handle the case where feature ~= maxrange
Index: kernel/src/unary_functions/SConscript
===================================================================
--- kernel/src/unary_functions/SConscript	(revision 626)
+++ kernel/src/unary_functions/SConscript	(working copy)
@@ -1,6 +1,7 @@
 Import('env')
-
-files = ['OpenCubicSpline.cpp', 'ClosedCubicSpline.cpp', 'Cosine.cpp']
-
-files = [File(x) for x in files]
+files=[
+        File( 'ClosedCubicSpline.cpp' ),
+        File( 'Cosine.cpp' ),
+        File( 'OpenCubicSpline.cpp' ),
+      ]
 Return('files')
Index: kernel/src/unary_functions/ClosedCubicSpline.cpp
===================================================================
--- kernel/src/unary_functions/ClosedCubicSpline.cpp	(revision 626)
+++ kernel/src/unary_functions/ClosedCubicSpline.cpp	(working copy)
@@ -55,7 +55,7 @@
 }
 
 
-Float ClosedCubicSpline::evaluate(Float feature)
+Float ClosedCubicSpline::evaluate(Float feature) const
 {
   // check for feature in range
   if (feature < minrange_ || feature > maxrange_) {
@@ -82,7 +82,7 @@
            * (spacing_ * spacing_) / 6.;
 }
 
-Float ClosedCubicSpline::evaluate_deriv(Float feature, Float& deriv)
+Float ClosedCubicSpline::evaluate_deriv(Float feature, Float& deriv) const
 {
   size_t lowbin = static_cast<size_t>((feature - minrange_) / spacing_);
   size_t highbin = lowbin + 1;
Index: kernel/src/Restraint.cpp
===================================================================
--- kernel/src/Restraint.cpp	(revision 626)
+++ kernel/src/Restraint.cpp	(working copy)
@@ -17,6 +17,7 @@
 
 Restraint::Restraint()
 {
+  model_=NULL;
   is_active_ = true; // active by default
 }
 
@@ -24,6 +25,12 @@
 //! Destructor
 Restraint::~Restraint()
 {
+  if (!was_owned_) {
+    // can't use virtual functions in the destructor
+    std::cerr << "Restraint " << this << " is being destroyed "
+              << "without ever having been added to a model."
+              << std::endl;
+  }
 }
 
 
@@ -43,9 +50,10 @@
 {
   IMP_assert(model==NULL || get_number_of_particles()==0
              || model == get_particle(0)->get_model()
-             || (model_ && model_.get() == model),
+             || (model_ && model_ == model),
              "Model* different from Particle Model*");
   model_=model;
+  was_owned_=true;
 }
 
 void Restraint::show(std::ostream& out) const
Index: kernel/src/ScoreState.cpp
===================================================================
--- kernel/src/ScoreState.cpp	(revision 626)
+++ kernel/src/ScoreState.cpp	(working copy)
@@ -17,6 +17,7 @@
 //! Constructor
 ScoreState::ScoreState(std::string name) : name_(name)
 {
+  model_=NULL;
   update_iteration_= std::numeric_limits<unsigned int>::max();
   after_iteration_= std::numeric_limits<unsigned int>::max();
   IMP_LOG(VERBOSE, "ScoreState constructed " << name << std::endl);
Index: kernel/src/score_states/BondDecoratorListScoreState.cpp
===================================================================
--- kernel/src/score_states/BondDecoratorListScoreState.cpp	(revision 626)
+++ kernel/src/score_states/BondDecoratorListScoreState.cpp	(working copy)
@@ -21,7 +21,7 @@
 {
   IMP_LOG(TERSE, "Updating BondDecoratorList for "
           << ps_.size() << " particles" << std::endl);
-  bonds_.clear();
+  std::vector<BondDecorator> bonds;
   for (unsigned int i=0; i< ps_.size(); ++i) {
     if (!ps_[i]->get_is_active()) continue;
     BondedDecorator di(ps_[i]);
@@ -36,31 +36,26 @@
       }
       if (di < dj) {
         IMP_LOG(VERBOSE, "Found bond " << di.get_bond(j) << std::endl);
-        bonds_.push_back(di.get_bond(j));
+        bonds.push_back(di.get_bond(j));
       }
     }
   }
-  IMP_LOG(TERSE, "Found " << bonds_.size() << " bonds"<< std::endl);
+  ManualBondDecoratorListScoreState::set_bond_decorators(bonds);
+  IMP_LOG(TERSE, "Found " << bonds.size() << " bonds"<< std::endl);
 }
 
 void BondDecoratorListScoreState::set_particles(const Particles &ps)
 {
   ps_=ps;
   std::sort(ps_.begin(), ps_.end());
-  bonds_.clear();
+  clear_bond_decorators();
 }
 
-
-bool BondDecoratorListScoreState::are_bonded(Particle *a, Particle *b) const
+void BondDecoratorListScoreState::add_particles(const Particles &ps)
 {
-  try {
-    BondedDecorator da= BondedDecorator::cast(a);
-    BondedDecorator db= BondedDecorator::cast(b);
-    return get_bond(da, db) != BondDecorator();
-  } catch (...) {
-    IMP_LOG(VERBOSE, "Exception thrown in are_bonded"<< std::endl);
-  }
-  return false;
+  ps_.insert(ps_.end(), ps.begin(), ps.end());
+  std::sort(ps_.begin(), ps_.end());
+  clear_bond_decorators();
 }
 
 } // namespace IMP
Index: kernel/src/score_states/NonbondedListScoreState.cpp
===================================================================
--- kernel/src/score_states/NonbondedListScoreState.cpp	(revision 626)
+++ kernel/src/score_states/NonbondedListScoreState.cpp	(working copy)
@@ -6,7 +6,7 @@
  */
 
 #include "IMP/score_states/NonbondedListScoreState.h"
-#include "IMP/decorators/XYZDecorator.h"
+#include "IMP/decorators/XYZRDecorator.h"
 #include "IMP/internal/Grid3D.h"
 #include "IMP/score_states/MaxChangeScoreState.h"
 
@@ -32,16 +32,15 @@
 
 
 NonbondedListScoreState
-::NonbondedListScoreState(Float cut,
-                          FloatKey rk): rk_(rk),
-                                        cutoff_(cut),
-                                        nbl_is_valid_(false)
+::NonbondedListScoreState(Float cut): cutoff_(cut),
+                                      nbl_is_valid_(false)
 {
   slack_=cutoff_;
   number_of_updates_=1;
   number_of_rebuilds_=0;
   number_of_overflows_=0;
   max_nbl_size_= std::numeric_limits<unsigned int>::max();
+  set_radius_key(XYZRDecorator::get_radius_key());
 }
 
 
@@ -51,6 +50,48 @@
 {
 }
 
+void NonbondedListScoreState::set_radius_key(FloatKey rk) {
+  rk_=rk;
+  Particles old_particles;
+  if (mcr_) {
+    old_particles= mcr_->get_particles();
+  }
+  if (rk != FloatKey()) {
+    FloatKeys ks;
+    ks.push_back(rk);
+    mcr_= new MaxChangeScoreState(ks);
+  } else {
+    mcr_= new MaxChangeScoreState(FloatKeys());
+  }
+  mcr_->add_particles(old_particles);
+}
+
+
+Float NonbondedListScoreState::get_max_radius_change() const {
+  if (rk_ != FloatKey()) {
+    return mcr_->get_max();
+  } else {
+    return 0;
+  }
+}
+
+void NonbondedListScoreState::reset_max_radius_change() {
+  if (mcr_) mcr_->reset();
+}
+
+void NonbondedListScoreState::update_max_radius_change(unsigned int i) {
+  if (mcr_) mcr_->before_evaluate(i);
+}
+
+void NonbondedListScoreState
+::add_particles_for_max_radius_change(const Particles &ps) {
+  mcr_->add_particles(particles_with_radius(ps));
+}
+
+void NonbondedListScoreState::clear_particles_for_max_radius_change() {
+  mcr_->clear_particles();
+}
+
 void NonbondedListScoreState::show_statistics(std::ostream &out) const
 {
   out << "Nonbonded list averaged " 
Index: kernel/src/score_states/CoverBondsScoreState.cpp
===================================================================
--- kernel/src/score_states/CoverBondsScoreState.cpp	(revision 626)
+++ kernel/src/score_states/CoverBondsScoreState.cpp	(working copy)
@@ -14,8 +14,9 @@
 {
 
 
-CoverBondsScoreState::CoverBondsScoreState(BondDecoratorListScoreState *bl,
-                                           FloatKey rk): bl_(bl), rk_(rk)
+CoverBondsScoreState
+::CoverBondsScoreState(ManualBondDecoratorListScoreState *bl,
+                       FloatKey rk): bl_(bl), rk_(rk)
 {
 }
 
@@ -25,8 +26,9 @@
 
 void CoverBondsScoreState::do_before_evaluate()
 {
-  for (BondDecoratorListScoreState::BondIterator it= bl_->bonds_begin();
-       it != bl_->bonds_end(); ++it) {
+  for (BondDecoratorListScoreState::BondDecoratorIterator it
+         = bl_->bond_decorators_begin();
+       it != bl_->bond_decorators_end(); ++it) {
     BondDecorator bd= *it;
     BondedDecorator pa= bd.get_bonded(0);
     BondedDecorator pb= bd.get_bonded(1);
@@ -62,8 +64,9 @@
 void CoverBondsScoreState::after_evaluate(DerivativeAccumulator *dva)
 {
   if (dva) {
-    for (BondDecoratorListScoreState::BondIterator it= bl_->bonds_begin();
-         it != bl_->bonds_end(); ++it) {
+    for (BondDecoratorListScoreState::BondDecoratorIterator it
+           = bl_->bond_decorators_begin();
+         it != bl_->bond_decorators_end(); ++it) {
       XYZDecorator d(it->get_particle());
       Vector3D deriv;
       // divide derivatives equally between endpoints
@@ -84,7 +87,7 @@
 
 void CoverBondsScoreState::show(std::ostream &out) const
 {
-  out << "CoverBondsScoreState on " << bl_->get_number_of_bonds()
+  out << "CoverBondsScoreState on " << bl_->get_number_of_bond_decorators()
       << " bonds " << std::endl;
 }
 
Index: kernel/src/score_states/BipartiteNonbondedListScoreState.cpp
===================================================================
--- kernel/src/score_states/BipartiteNonbondedListScoreState.cpp	(revision 626)
+++ kernel/src/score_states/BipartiteNonbondedListScoreState.cpp	(working copy)
@@ -15,8 +15,10 @@
 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
+
+
+// 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
@@ -45,30 +47,21 @@
 
 
 BipartiteNonbondedListScoreState
-::BipartiteNonbondedListScoreState(Float cut,
-                                   FloatKey rk,
-                                   Algorithm a):
-  P(cut, rk), a_(translate_algorithm(a))
+::BipartiteNonbondedListScoreState(Float cut):
+  P(cut), a_(translate_algorithm(DEFAULT))
 {
-  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,
+::BipartiteNonbondedListScoreState(Float cut,
                                    const Particles &ps0,
-                                   const Particles &ps1,
-                                   Algorithm a):
-  P(cut, rk), a_(translate_algorithm(a))
+                                   const Particles &ps1):
+  P(cut), a_(translate_algorithm(DEFAULT))
 {
-  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);
 }
 
@@ -79,7 +72,7 @@
 {
   mc0_->before_evaluate(ScoreState::get_before_evaluate_iteration());
   mc0_->before_evaluate(ScoreState::get_before_evaluate_iteration());
-  mcr_->before_evaluate(ScoreState::get_before_evaluate_iteration());
+  update_max_radius_change(ScoreState::get_before_evaluate_iteration());
   Float mc= std::max(mc0_->get_max(), mc1_->get_max());
 
   Float cost;
@@ -88,17 +81,19 @@
     cost= 10*mc0_->get_number_of_particles()* mc1_->get_number_of_particles();
     break;
   case BBOX:
-    cost= 1000 * mcr_->get_number_of_particles();
+    cost= 1000 * (mc0_->get_number_of_particles()
+                  + mc1_->get_number_of_particles());
     break;
   default:
     IMP_assert(0, "Bad algorithm");
-    cost= 1000 * mcr_->get_number_of_particles();
+    cost= 1000 * (mc0_->get_number_of_particles()
+                  + mc1_->get_number_of_particles());
   }
 
-  if (P::update(mc+ mcr_->get_max(), cost)) {
+  if (P::update(mc+ get_max_radius_change(), cost)) {
     mc0_->reset();
     mc1_->reset();
-    mcr_->reset();
+    reset_max_radius_change();
   }
   IMP_IF_CHECK(EXPENSIVE) {
     check_nbl();
@@ -123,7 +118,7 @@
     break;
   default:
     IMP_failure("Can't find algorithm in BipartiteNonbondedListScoreState",
-                ErrorException());
+                ErrorException);
   }
 }
 
@@ -150,9 +145,9 @@
   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));
+  clear_particles_for_max_radius_change();
+  add_particles_for_max_radius_change(ps0);
+  add_particles_for_max_radius_change(ps1);
   P::set_nbl_is_valid(false);
 }
 
@@ -160,24 +155,36 @@
 void BipartiteNonbondedListScoreState::add_particles_0(const Particles &ps)
 {
   if (P::get_nbl_is_valid()) process_sets(ps, mc1_->get_particles());
+  else P::set_nbl_is_valid(false);
   mc0_->add_particles(ps);
-  mcr_->add_particles(P::particles_with_radius(ps));
-  P::set_nbl_is_valid(false);
+  add_particles_for_max_radius_change(ps);
 }
 
 void BipartiteNonbondedListScoreState::add_particles_1(const Particles &ps)
 {
   if (P::get_nbl_is_valid()) process_sets(ps, mc0_->get_particles());
+  else P::set_nbl_is_valid(false);
   mc1_->add_particles(ps);
-  mcr_->add_particles(P::particles_with_radius(ps));
-  P::set_nbl_is_valid(false);
+  add_particles_for_max_radius_change(ps);
 }
 
+
+void BipartiteNonbondedListScoreState::add_particle_0(Particle *ps)
+{
+  add_particles_0(Particles(1, ps));
+}
+
+void BipartiteNonbondedListScoreState::add_particle_1(Particle *ps)
+{
+  add_particles_1(Particles(1, ps));
+}
+
+
 void BipartiteNonbondedListScoreState::clear_particles()
 {
   mc0_->clear_particles();
   mc1_->clear_particles();
-  mcr_->clear_particles();
+  clear_particles_for_max_radius_change();
   P::set_nbl_is_valid(false);
   P::set_nbl_is_valid(true);
 }
Index: kernel/src/score_states/SConscript
===================================================================
--- kernel/src/score_states/SConscript	(revision 626)
+++ kernel/src/score_states/SConscript	(working copy)
@@ -1,10 +1,12 @@
 Import('env')
-
-files = ['NonbondedListScoreState.cpp',
-         'MaxChangeScoreState.cpp', 'BondDecoratorListScoreState.cpp',
-         'AllNonbondedListScoreState.cpp',
-         'BipartiteNonbondedListScoreState.cpp',
-         'GravityCenterScoreState.cpp', 'CoverBondsScoreState.cpp']
-
-files = [File(x) for x in files]
+files=[
+        File( 'AllNonbondedListScoreState.cpp' ),
+        File( 'BipartiteNonbondedListScoreState.cpp' ),
+        File( 'BondDecoratorListScoreState.cpp' ),
+        File( 'CoverBondsScoreState.cpp' ),
+        File( 'GravityCenterScoreState.cpp' ),
+        File( 'ManualBondDecoratorListScoreState.cpp' ),
+        File( 'MaxChangeScoreState.cpp' ),
+        File( 'NonbondedListScoreState.cpp' ),
+      ]
 Return('files')
Index: kernel/src/score_states/AllNonbondedListScoreState.cpp
===================================================================
--- kernel/src/score_states/AllNonbondedListScoreState.cpp	(revision 626)
+++ kernel/src/score_states/AllNonbondedListScoreState.cpp	(working copy)
@@ -18,7 +18,7 @@
 {
 
 
-//! Turn the default into an actual algorithm and work around missing algorithms
+// Turn the default into an actual algorithm and work around missing algorithms
 static AllNonbondedListScoreState::Algorithm
 translate_algorithm(AllNonbondedListScoreState::Algorithm a)
 {
@@ -44,29 +44,13 @@
 
 
 AllNonbondedListScoreState::AllNonbondedListScoreState(Float cut,
-                                                       FloatKey rk,
-                                                       const Particles &ps,
-                                                       Algorithm a):
-  P(cut, rk), a_(translate_algorithm(a))
+                                                       const Particles &ps):
+  P(cut), a_(translate_algorithm(DEFAULT))
 {
-  FloatKeys ks;
-  ks.push_back(rk);
   mc_= new MaxChangeScoreState(XYZDecorator::get_xyz_keys());
-  mcr_= new MaxChangeScoreState(ks);
   add_particles(ps);
 }
 
-AllNonbondedListScoreState::AllNonbondedListScoreState(Float cut,
-                                                       FloatKey rk,
-                                                       Algorithm a):
-  P(cut, rk), a_(translate_algorithm(a))
-{
-  FloatKeys ks;
-  ks.push_back(rk);
-  mc_= new MaxChangeScoreState(XYZDecorator::get_xyz_keys());
-  mcr_= new MaxChangeScoreState(ks);
-}
-
 AllNonbondedListScoreState::~AllNonbondedListScoreState()
 {
 }
@@ -74,8 +58,8 @@
 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();
+  update_max_radius_change(ScoreState::get_before_evaluate_iteration());
+  Float mc=mc_->get_max()+ get_max_radius_change();
   Float cost;
   switch (a_){
   case QUADRATIC:
@@ -89,12 +73,12 @@
     cost = 2000 * mc_->get_number_of_particles();
     break;
   default:
-    IMP_failure("Bad algorithm", ErrorException());
+    IMP_failure("Bad algorithm", ErrorException);
     cost = 10 * mc_->get_number_of_particles();
   }
   if (P::update(mc, cost)) {
     mc_->reset();
-    mcr_->reset();
+    reset_max_radius_change();
   }
   IMP_IF_CHECK(EXPENSIVE) {
     check_nbl();
@@ -121,7 +105,7 @@
                         internal::NBLAddPairIfNonbonded(this));
 
   } else {
-    IMP_failure("Bad algorithm in AllNBL::rebuild", ErrorException());
+    IMP_failure("Bad algorithm in AllNBL::rebuild", ErrorException);
   }
   set_nbl_is_valid(true);
   IMP_LOG(TERSE, "NBL has " << P::get_number_of_nonbonded()
@@ -132,8 +116,8 @@
 {
   mc_->clear_particles();
   mc_->add_particles(ps);
-  mcr_->clear_particles();
-  mcr_->add_particles(particles_with_radius(ps));
+  clear_particles_for_max_radius_change();
+  add_particles_for_max_radius_change(ps);
   P::set_nbl_is_valid(false);
 }
 
@@ -155,14 +139,14 @@
                                     internal::NBLAddPairIfNonbonded(this));
     }
   }
-  mcr_->add_particles(particles_with_radius(ps));
+  add_particles_for_max_radius_change(ps);
   mc_->add_particles(ps);
 }
 
 void AllNonbondedListScoreState::clear_particles()
 {
   mc_->clear_particles();
-  mcr_->clear_particles();
+  clear_particles_for_max_radius_change();
   P::set_nbl_is_valid(false);
   P::set_nbl_is_valid(true);
 }
Index: kernel/src/score_states/MaxChangeScoreState.cpp
===================================================================
--- kernel/src/score_states/MaxChangeScoreState.cpp	(revision 626)
+++ kernel/src/score_states/MaxChangeScoreState.cpp	(working copy)
@@ -19,12 +19,6 @@
                                          const Particles &ps): keys_(keys)
 {
   add_particles(ps);
-  origkeys_.resize(keys_.size());
-  std::ostringstream oss;
-  oss << " MCSS base " << this;
-  for (unsigned int i=0; i< keys_.size(); ++i) {
-    origkeys_[i]= FloatKey((keys_[i].get_string()+oss.str()).c_str());
-  }
 }
 
 
@@ -35,12 +29,6 @@
                             << obj,
                             ValueException);
                 };
-                for (unsigned int i=0; i< origkeys_.size(); ++i) {
-                  if (!obj->has_attribute(origkeys_[i])) {
-                    obj->add_attribute(origkeys_[i],
-                                       obj->get_value(keys_[i]), false);
-                  }
-                }
               }, {reset();});
 
 void MaxChangeScoreState::do_before_evaluate()
@@ -48,12 +36,12 @@
   max_change_=0;
   // get rid of inactive particles and their stored values
   internal::remove_inactive_particles(particle_vector_);
-  for (ParticleIterator it= particles_begin(); it != particles_end(); ++it) {
-    (*it)->assert_is_valid();
+  for (unsigned int i=0; i< get_number_of_particles(); ++i) {
+    get_particle(i)->assert_is_valid();
     for (unsigned int j=0; j < keys_.size(); ++j) {
-      Float v= (*it)->get_value(keys_[j]);
-      Float ov= (*it)->get_value(origkeys_[j]);
-      IMP_LOG(VERBOSE, "Particle " << (*it)->get_index() 
+      Float v= get_particle(i)->get_value(keys_[j]);
+      Float ov= old_values_[i].get_value(keys_[j]);
+      IMP_LOG(VERBOSE, "Particle " << get_particle(i)->get_index() 
               << " and attribute " << keys_[j]
               << " moved " << std::abs(v - ov) << std::endl);
       max_change_= std::max(max_change_,
@@ -66,9 +54,15 @@
 
 void MaxChangeScoreState::reset()
 {
-  for (ParticleIterator it= particles_begin(); it != particles_end(); ++it) {
+  Table t;
+  for (unsigned int i=0; i < keys_.size(); ++i) {
+    t.insert(keys_[i], 0);
+  }
+  old_values_.resize(get_number_of_particles(), t);
+  for (unsigned int i=0; i< get_number_of_particles(); ++i) {
     for (unsigned int j=0; j < keys_.size(); ++j) {
-      (*it)->set_value(origkeys_[j], (*it)->get_value(keys_[j]));
+      old_values_[i].set_value(keys_[j],
+                               get_particle(i)->get_value(keys_[j]));
     }
   }
   max_change_=0;
Index: kernel/src/score_states/ManualBondDecoratorListScoreState.cpp
===================================================================
--- kernel/src/score_states/ManualBondDecoratorListScoreState.cpp	(revision 0)
+++ kernel/src/score_states/ManualBondDecoratorListScoreState.cpp	(revision 0)
@@ -0,0 +1,41 @@
+#include "IMP/score_states/ManualBondDecoratorListScoreState.h"
+
+namespace IMP
+{
+  ManualBondDecoratorListScoreState
+  ::ManualBondDecoratorListScoreState(const BondDecorators &ps) {
+    set_bond_decorators(ps);
+  }
+
+  IMP_LIST_IMPL(ManualBondDecoratorListScoreState, BondDecorator,
+                bond_decorator, BondDecorator,
+                ,);
+
+  bool ManualBondDecoratorListScoreState::are_bonded(Particle *a,
+                                                     Particle *b) const
+  {
+    if (BondedDecorator::is_instance_of(a)
+        && BondedDecorator::is_instance_of(b)) {
+      BondedDecorator da(a);
+      BondedDecorator db(b);
+      return get_bond(da, db) != BondDecorator();
+    } else if (BondDecorator::is_instance_of(a)
+               && BondDecorator::is_instance_of(b)){
+      // check if they are bonds which share an endpoint
+      /**
+         \todo Decide if the check for adject bonds should go elsewhere
+       */
+      BondDecorator da(a);
+      BondDecorator db(b);
+      Particle *a0= da.get_bonded(0).get_particle();
+      Particle *a1= da.get_bonded(1).get_particle();
+      Particle *b0= db.get_bonded(0).get_particle();
+      Particle *b1= db.get_bonded(1).get_particle();
+
+      return (a0== b0 || a0 == b1 || a1== b0 || a1 == b1);
+    } else {
+      return false;
+    }
+  }
+
+}
Index: kernel/src/restraints/TunnelRestraint.cpp
===================================================================
--- kernel/src/restraints/TunnelRestraint.cpp	(revision 626)
+++ kernel/src/restraints/TunnelRestraint.cpp	(working copy)
@@ -1,103 +0,0 @@
-/**
- *  \file TunnelRestraint.cpp     \brief Tunnel restraint.
- *
- *  Copyright 2007-8 Sali Lab. All rights reserved.
- *
- */
-
-#include "IMP/restraints/TunnelRestraint.h"
-#include "IMP/decorators/XYZDecorator.h"
-
-namespace IMP
-{
-
-TunnelRestraint::TunnelRestraint(UnaryFunction *f, FloatKey r) : f_(f), rk_(r)
-{
-  coordinate_=0;
-  center_=Vector3D(0,0,0);
-  height_=0;
-  radius_=0;
-}
-
-
-
-Float TunnelRestraint::evaluate(DerivativeAccumulator *accum)
-{
-  Float ret=0;
-  for (Restraint::ParticleIterator it= Restraint::particles_begin();
-       it != Restraint::particles_end(); ++it) {
-    XYZDecorator d(*it);
-    Float radius=0;
-    if (rk_ != FloatKey() && (*it)->has_attribute(rk_)) {
-      radius= (*it)->get_value(rk_);
-    }
-    IMP_LOG(VERBOSE, "Tunnel restraint on particle " << d << std::endl);
-    Float hr= height_+radius;
-    if (d.get_coordinate(coordinate_) < center_[coordinate_] + hr
-        && d.get_coordinate(coordinate_) > center_[coordinate_] - hr) {
-      Float sd=0;
-      for (int i=1; i< 3; ++i) {
-        int oc= (i+coordinate_)%3;
-        sd+= square(d.get_coordinate(oc)- center_[oc]);
-      }
-      sd= std::sqrt(sd);
-      IMP_LOG(VERBOSE, "The distance is " << sd << " and radius "
-              << radius_ << std::endl);
-      if (sd > radius_-radius) {
-        Float rd= sd-radius_;
-        Float hdu= center_[coordinate_] + height_
-          - d.get_coordinate(coordinate_);
-        Float hdd=  d.get_coordinate(coordinate_)
-          + height_- center_[coordinate_];
-        Vector3D deriv(0,0,0);
-        Float score=0;
-        Float deriv_scalar=0;
-        /*! \todo Clean up these tests so I am not dependent on two expressions
-          being the same and evaluating to the same thing */
-        // look below if changed
-        Float dist= -std::min(std::min(rd, hdu), hdd) - radius;
-        if (accum) {
-          score= f_->evaluate_deriv(dist, deriv_scalar);
-        } else {
-          score= f_->evaluate(dist);
-        }
-
-        // kind if evil
-        if (dist== -rd -radius) {
-          Vector3D v= (d.get_vector() - center_).get_unit_vector();
-          for (int i=0; i< 2; ++i) {
-            int oc= (i+coordinate_+1)%3;
-            deriv[oc]= v[oc]*deriv_scalar;
-          }
-        } else {
-          // kind of evil
-          if (dist == -hdu -radius) deriv_scalar= -deriv_scalar;
-          deriv[coordinate_]= deriv_scalar;
-        }
-
-        if (accum) {
-          for (unsigned int i=0; i< 3; ++i) {
-            d.add_to_coordinate_derivative(i, deriv[i], *accum);
-          }
-        }
-        ret+= score;
-      } else {
-        IMP_LOG(VERBOSE, "Particle " << (*it)->get_index()
-                << " is in channel" << std::endl);
-      }
-    } else {
-      IMP_LOG(VERBOSE, "Particle " << (*it)->get_index()
-              << " is outside of slab" << std::endl);
-    }
-  }
-  return ret;
-}
-
-void TunnelRestraint::show(std::ostream& out) const
-{
-  out << "Tunnel restraint :" << f_ 
-      << " " << center_ << " " << height_ << " " 
-      << coordinate_ << std::endl;
-}
-
-}  // namespace IMP
Index: kernel/src/restraints/SConscript
===================================================================
--- kernel/src/restraints/SConscript	(revision 626)
+++ kernel/src/restraints/SConscript	(working copy)
@@ -1,11 +1,17 @@
 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',
-         'ConstantRestraint.cpp', 'TunnelRestraint.cpp']
-
-files = [File(x) for x in files]
+files=[
+        File( 'AngleRestraint.cpp' ),
+        File( 'BondDecoratorRestraint.cpp' ),
+        File( 'ConnectivityRestraint.cpp' ),
+        File( 'ConstantRestraint.cpp' ),
+        File( 'DihedralRestraint.cpp' ),
+        File( 'DistanceRestraint.cpp' ),
+        File( 'NonbondedRestraint.cpp' ),
+        File( 'PairChainRestraint.cpp' ),
+        File( 'PairListRestraint.cpp' ),
+        File( 'PairRestraint.cpp' ),
+        File( 'RestraintSet.cpp' ),
+        File( 'SingletonListRestraint.cpp' ),
+        File( 'TripletChainRestraint.cpp' ),
+      ]
 Return('files')
Index: kernel/src/restraints/ConnectivityRestraint.cpp
===================================================================
--- kernel/src/restraints/ConnectivityRestraint.cpp	(revision 626)
+++ kernel/src/restraints/ConnectivityRestraint.cpp	(working copy)
@@ -118,6 +118,8 @@
     IMP_LOG(VERBOSE, "ConnectivityRestraint edge between "
             << get_particle(i)->get_index()
             << " and " << get_particle(j)->get_index() << std::endl);
+    IMP_LOG(VERBOSE, *get_particle(i) << std::endl
+            << *get_particle(j) << std::endl);
     sum+= ps_->evaluate(get_particle(i),
                         get_particle(j),
                         accum);
Index: kernel/src/restraints/BondDecoratorRestraint.cpp
===================================================================
--- kernel/src/restraints/BondDecoratorRestraint.cpp	(revision 626)
+++ kernel/src/restraints/BondDecoratorRestraint.cpp	(working copy)
@@ -20,14 +20,15 @@
 
 BondDecoratorRestraint
 ::BondDecoratorRestraint(UnaryFunction *f,
-                         BondDecoratorListScoreState *s): bl_(s),
-                                                          f_(f){}
+                         ManualBondDecoratorListScoreState *s): bl_(s),
+                                                                f_(f){}
 
 Float BondDecoratorRestraint::evaluate(DerivativeAccumulator *accum)
 {
   Float sum=0;
-  for (BondDecoratorListScoreState::BondIterator bi= bl_->bonds_begin();
-       bi != bl_->bonds_end(); ++bi) {
+  for (BondDecoratorListScoreState::BondDecoratorIterator bi
+         = bl_->bond_decorators_begin();
+       bi != bl_->bond_decorators_end(); ++bi) {
     BondDecorator bd= *bi;
     Float l= bd.get_length();
     Float s= bd.get_stiffness();
@@ -64,7 +65,7 @@
 {
   out << "Bond decorator restraint with unary function ";
   f_->show(out);
-  out << " on " << bl_->get_number_of_bonds() << " bonds";
+  out << " on " << bl_->get_number_of_bond_decorators() << " bonds";
   out << std::endl;
 }
 
Index: kernel/src/restraints/DistanceRestraint.cpp
===================================================================
--- kernel/src/restraints/DistanceRestraint.cpp	(revision 626)
+++ kernel/src/restraints/DistanceRestraint.cpp	(working copy)
@@ -17,23 +17,10 @@
 
 DistanceRestraint::DistanceRestraint(UnaryFunction* score_func,
                                      Particle* p1, Particle* p2) :
-    dp_(score_func)
+  PairRestraint(new DistancePairScore(score_func), p1, p2)
 {
-  add_particle(p1);
-  add_particle(p2);
 }
 
-//! Calculate the score for this distance restraint.
-/** \param[in] accum If not NULL, use this object to accumulate partial first
-                     derivatives.
-    \return Current score.
- */
-Float DistanceRestraint::evaluate(DerivativeAccumulator *accum)
-{
-  return dp_.evaluate(get_particle(0), get_particle(1), accum);
-}
-
-
 //! Show the current restraint.
 /** \param[in] out Stream to send restraint description to.
  */
@@ -49,7 +36,7 @@
   out << "  particles: " << get_particle(0)->get_index();
   out << " and " << get_particle(1)->get_index();
   out << "  ";
-  dp_.show(out);
+  dp_->show(out);
   out << std::endl;
 }
 
Index: kernel/src/restraints/PairRestraint.cpp
===================================================================
--- kernel/src/restraints/PairRestraint.cpp	(revision 0)
+++ kernel/src/restraints/PairRestraint.cpp	(revision 0)
@@ -0,0 +1,51 @@
+/**
+ *  \file PairRestraint.cpp \brief Pair restraint between two particles.
+ *
+ *  Copyright 2007-8 Sali Lab. All rights reserved.
+ *
+ */
+
+#include "IMP/Particle.h"
+#include "IMP/Model.h"
+#include "IMP/log.h"
+#include "IMP/restraints/PairRestraint.h"
+#include "IMP/decorators/XYZDecorator.h"
+#include "IMP/PairScore.h"
+
+namespace IMP
+{
+
+PairRestraint::PairRestraint(PairScore* score_func,
+                             Particle* p1, Particle* p2) :
+    dp_(score_func)
+{
+  add_particle(p1);
+  add_particle(p2);
+}
+
+Float PairRestraint::evaluate(DerivativeAccumulator *accum)
+{
+  return dp_->evaluate(get_particle(0), get_particle(1), accum);
+}
+
+
+//! Show the current restraint.
+/** \param[in] out Stream to send restraint description to.
+ */
+void PairRestraint::show(std::ostream& out) const
+{
+  if (get_is_active()) {
+    out << "pair restraint (active):" << std::endl;
+  } else {
+    out << "pair restraint (inactive):" << std::endl;
+  }
+
+  get_version_info().show(out);
+  out << "  particles: " << get_particle(0)->get_index();
+  out << " and " << get_particle(1)->get_index();
+  out << "  ";
+  dp_->show(out);
+  out << std::endl;
+}
+
+}  // namespace IMP
Index: kernel/src/restraints/PairChainRestraint.cpp
===================================================================
--- kernel/src/restraints/PairChainRestraint.cpp	(revision 626)
+++ kernel/src/restraints/PairChainRestraint.cpp	(working copy)
@@ -21,51 +21,20 @@
 PairChainRestraint::PairChainRestraint(PairScore* ts)
 {
   ts_ = ts;
-  clear_chains();
 }
 
-void PairChainRestraint::add_chain(const Particles &ps)
-{
-  if (ps.size() <= 1) {
-    IMP_WARN("Adding a chain of length 1 or less to the PairChainsRestraint"
-             << " doesn't accomplish anything."<< std::endl);
-  } else {
-    Restraint::add_particles(ps);
-    chain_splits_.back()= Restraint::get_number_of_particles();
-    chain_splits_.push_back(Restraint::get_number_of_particles());
-  }
-}
 
 Float PairChainRestraint::evaluate(DerivativeAccumulator *accum)
 {
-  int cur_break=0;
-  unsigned int i=1;
   float score=0;
-  while (i < Restraint::get_number_of_particles()) {
-    /*IMP_LOG(VERBOSE, "Chain eval on " 
-            << Restraint::get_particle(i-2)->get_index()
-            << Restraint::get_particle(i-1)->get_index()
-            << Restraint::get_particle(i)->get_index()
-            << " split is " << chain_splits_[cur_break]
-            << std::endl);*/
+  for (unsigned int i=1; i< Restraint::get_number_of_particles(); ++i) {
     score += ts_->evaluate(Restraint::get_particle(i-1),
                            Restraint::get_particle(i),
                            accum);
-    if (chain_splits_[cur_break] == i) {
-      i+=2;
-      ++cur_break;
-    } else {
-      ++i;
-    }
   } 
   return score;
 }
 
-void PairChainRestraint::clear_chains() {
-  Restraint::clear_particles();
-  chain_splits_.clear();
-  chain_splits_.push_back(0);
-}
 
 
 void PairChainRestraint::show(std::ostream& out) const
@@ -77,7 +46,6 @@
   }
 
   get_version_info().show(out);
-  out << "  " << chain_splits_.size()-1 << " chains" << std::endl;
   ts_->show(out);
   out << std::endl;
 }
Index: kernel/src/ParticleRefiner.cpp
===================================================================
--- kernel/src/ParticleRefiner.cpp	(revision 626)
+++ kernel/src/ParticleRefiner.cpp	(working copy)
@@ -12,8 +12,8 @@
 
 ParticleRefiner::~ParticleRefiner(){}
 
-Particles ParticleRefiner::get_refined(Particle *p) {
-  throw ErrorException("Can't refine");
+Particles ParticleRefiner::get_refined(Particle *p) const {
+  IMP_failure("get_refined not implemented. Oops.", ErrorException);
   return Particles();
 }
 
Index: kernel/src/particle_refiners/BondCoverParticleRefiner.cpp
===================================================================
--- kernel/src/particle_refiners/BondCoverParticleRefiner.cpp	(revision 626)
+++ kernel/src/particle_refiners/BondCoverParticleRefiner.cpp	(working copy)
@@ -16,12 +16,11 @@
 {
 
 BondCoverParticleRefiner::BondCoverParticleRefiner(FloatKey rk,
-                                                   FloatKey vk): rk_(rk),
-                                                                 vk_(vk)
+                                                   FloatKey vk,
+                                                   IntKey tk): rk_(rk),
+                                                               vk_(vk),
+                                                               tk_(tk)
 {
-  if (0) {
-    BondCoverParticleRefiner t(rk, vk);
-  }
 }
 
 
@@ -35,7 +34,7 @@
 /* n (4/3) pi (d/(2n))^3 = v
    n^2= (4/3) pi (d/2)^3 / v
  */
-Particles BondCoverParticleRefiner::get_refined(Particle *p)
+Particles BondCoverParticleRefiner::get_refined(Particle *p) const
 {
   IMP_assert(get_can_refine(p), "Trying to refine the unrefinable");
 
@@ -98,6 +97,10 @@
     XYZDecorator d= XYZDecorator::create(np);
     d.set_coordinates(vb+ (1+2*i)*r* ud);
     np->add_attribute(rk_, r, false);
+    if (tk_ != IntKey() && p->has_attribute(tk_)) {
+      np->add_attribute(tk_,
+                        p->get_value(tk_));
+    }
     ret.push_back(np);
   }
   return ret;
@@ -107,7 +110,7 @@
 
 void BondCoverParticleRefiner::cleanup_refined(Particle *p,
                                                Particles &ps,
-                                               DerivativeAccumulator *da)
+                                               DerivativeAccumulator *da) const
 {
   IMP_assert(get_can_refine(p), "Cleanup called with non-refinable particle");
   BondDecorator bd(p);
Index: kernel/src/particle_refiners/SConscript
===================================================================
--- kernel/src/particle_refiners/SConscript	(revision 626)
+++ kernel/src/particle_refiners/SConscript	(working copy)
@@ -1,4 +1,6 @@
-files=['BondCoverParticleRefiner.cpp', 'ChildrenParticleRefiner.cpp']
-
-files = [File(x) for x in files]
+Import('env')
+files=[
+        File( 'BondCoverParticleRefiner.cpp' ),
+        File( 'ChildrenParticleRefiner.cpp' ),
+      ]
 Return('files')
Index: kernel/src/particle_refiners/ChildrenParticleRefiner.cpp
===================================================================
--- kernel/src/particle_refiners/ChildrenParticleRefiner.cpp	(revision 626)
+++ kernel/src/particle_refiners/ChildrenParticleRefiner.cpp	(working copy)
@@ -23,7 +23,7 @@
 
 }
 
-Particles ChildrenParticleRefiner::get_refined(Particle *p)
+Particles ChildrenParticleRefiner::get_refined(Particle *p) const
 {
   IMP_assert(get_can_refine(p), "Trying to refine the unrefinable");
   HierarchyDecorator d(p);
@@ -38,7 +38,7 @@
 
 void ChildrenParticleRefiner::cleanup_refined(Particle *,
                                               Particles &,
-                                              DerivativeAccumulator *)
+                                              DerivativeAccumulator *) const
 {
   // This space left intentionally blank
 }
Index: kernel/src/triplet_scores/AngleTripletScore.cpp
===================================================================
--- kernel/src/triplet_scores/AngleTripletScore.cpp	(revision 626)
+++ kernel/src/triplet_scores/AngleTripletScore.cpp	(working copy)
@@ -15,7 +15,7 @@
 AngleTripletScore::AngleTripletScore(UnaryFunction *f): f_(f){}
 
 Float AngleTripletScore::evaluate(Particle *a, Particle *b, Particle *c,
-                                  DerivativeAccumulator *da)
+                                  DerivativeAccumulator *da) const
 {
   IMP_CHECK_OBJECT(f_.get());
   IMP_CHECK_OBJECT(a);
Index: kernel/src/triplet_scores/SConscript
===================================================================
--- kernel/src/triplet_scores/SConscript	(revision 626)
+++ kernel/src/triplet_scores/SConscript	(working copy)
@@ -1,6 +1,5 @@
 Import('env')
-
-files = ['AngleTripletScore.cpp']
-
-files = [File(x) for x in files]
+files=[
+        File( 'AngleTripletScore.cpp' ),
+      ]
 Return('files')
Index: kernel/src/Log.cpp
===================================================================
--- kernel/src/Log.cpp	(revision 626)
+++ kernel/src/Log.cpp	(working copy)
@@ -14,4 +14,13 @@
 /* Initialize singleton pointer to NULL */
 Log* Log::logpt_ = NULL;
 
+
+bool get_print_exception_messages() {
+  return internal::print_exceptions;
+}
+
+void set_print_exception_messages(bool tf) {
+  internal::print_exceptions=tf;
+}
+
 } // namespace IMP
Index: kernel/src/OptimizerState.cpp
===================================================================
--- kernel/src/OptimizerState.cpp	(revision 626)
+++ kernel/src/OptimizerState.cpp	(working copy)
@@ -15,6 +15,7 @@
 //! Constructor
 OptimizerState::OptimizerState() 
 {
+  optimizer_=NULL;
   IMP_LOG(VERBOSE, "OptimizerState constructed " << std::endl);
 }
 
Index: kernel/src/decorators/XYZRDecorator.cpp
===================================================================
--- kernel/src/decorators/XYZRDecorator.cpp	(revision 0)
+++ kernel/src/decorators/XYZRDecorator.cpp	(revision 0)
@@ -0,0 +1,102 @@
+/**
+ *  \file XYZRDecorator.cpp   \brief Simple xyzr decorator.
+ *
+ *  Copyright 2007-8 Sali Lab. All rights reserved.
+ *
+ */
+
+#include "IMP/decorators/XYZRDecorator.h"
+
+//#undef IMP_USE_CGAL
+
+#ifdef IMP_USE_CGAL
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Min_sphere_of_spheres_d.h>
+
+#endif
+
+namespace IMP
+{
+
+// These aren't statically initialized, as that way they may be initialized
+// before the table that caches them
+FloatKey XYZRDecorator::radius_;
+
+void XYZRDecorator::show(std::ostream &out, std::string prefix) const
+{
+  out << prefix << "(" << get_x()<< ", "
+      << get_y() << ", " << get_z() << ": " << get_radius() << ")";
+
+}
+
+
+
+IMP_DECORATOR_INITIALIZE(XYZRDecorator, XYZDecorator,
+                         {
+                         radius_= FloatKey("radius");
+                         })
+
+
+Float distance(XYZRDecorator a, XYZRDecorator b)
+{
+  float d= distance(static_cast<XYZDecorator>(a),
+                    static_cast<XYZDecorator>(b));
+  return d - a.get_radius() - b.get_radius();
+}
+
+
+void set_enclosing_sphere(const Particles &v,
+                          XYZRDecorator out)
+{
+  IMP_check(!v.empty(), "Must pass some particles to have a bounding sphere",
+            ValueException);
+  for (unsigned int i=0; i< v.size(); ++i) {
+    XYZDecorator::cast(v[i]);
+  }
+  FloatKey rk= XYZRDecorator::get_radius_key();
+
+
+#ifdef IMP_USE_CGAL
+  typedef CGAL::Exact_predicates_inexact_constructions_kernel             K;
+  typedef CGAL::Min_sphere_of_spheres_d_traits_3<K, K::FT> Traits;
+  typedef CGAL::Min_sphere_of_spheres_d<Traits> Min_sphere;
+  typedef K::Point_3                        Point;
+  typedef Traits::Sphere                    Sphere;
+
+  std::vector<Sphere> spheres;
+  for (unsigned int i=0; i< v.size(); ++i) {
+    // need cast to resolve ambiguity
+    XYZDecorator d(v[i]);
+    float r=0;
+    if (v[i]->has_attribute(rk)) r= v[i]->get_value(rk);
+
+    spheres.push_back(Sphere(Point(d.get_x(),
+                                   d.get_y(),
+                                   d.get_z()),
+                             square(r)));
+  }
+  Min_sphere ms(spheres.begin(), spheres.end());
+
+  out.set_radius(ms.radius());
+  out.set_x(*ms.center_cartesian_begin());
+  out.set_y(*(ms.center_cartesian_begin()+1));
+  out.set_z(*(ms.center_cartesian_begin()+2));
+#else
+  Vector3D c(0,0,0);
+  for (unsigned int i=0; i< v.size(); ++i) {
+    XYZDecorator d(v[i]);
+    c+= d.get_vector();
+  }
+  c/=v.size();
+  out.set_coordinates(c);
+  float r=0;
+  for (unsigned int i=0; i< v.size(); ++i) {
+    float d= distance(XYZDecorator(v[i]), out);
+    if (v[i]->has_attribute(rk)) d+= v[i]->get_value(rk);
+    r= std::max(r, d);
+  }
+  out.set_radius(r);
+#endif
+}
+
+} // namespace IMP
Index: kernel/src/decorators/HierarchyDecorator.cpp
===================================================================
--- kernel/src/decorators/HierarchyDecorator.cpp	(revision 626)
+++ kernel/src/decorators/HierarchyDecorator.cpp	(working copy)
@@ -10,6 +10,7 @@
 #include "IMP/decorators/NameDecorator.h"
 
 #include <sstream>
+#include <memory>
 
 namespace IMP
 {
@@ -44,9 +45,9 @@
 
 unsigned int count_hierarchy(HierarchyDecorator h)
 {
-  HierarchyCounter hc;
-  depth_first_traversal(h,hc);
-  return hc.get_count();
+  std::auto_ptr<HierarchyCounter> hc(new HierarchyCounter());
+  depth_first_traversal(h,hc.get());
+  return hc->get_count();
 }
 
 
@@ -70,8 +71,8 @@
 void HierarchyDecorator::validate() const
 {
   //std::cerr << "Checking hierarchy" << std::endl;
-  internal::AssertHierarchy ah;
-  depth_first_traversal(*this, ah);
+  std::auto_ptr<internal::AssertHierarchy> ah(new internal::AssertHierarchy());
+  depth_first_traversal(*this, ah.get());
 }
 
 
@@ -95,7 +96,7 @@
                          })
 
 
-void breadth_first_traversal(HierarchyDecorator d, HierarchyVisitor &f)
+void breadth_first_traversal(HierarchyDecorator d, HierarchyVisitor *f)
 {
   std::deque<HierarchyDecorator> stack;
   stack.push_back(d);
@@ -103,7 +104,7 @@
   do {
     HierarchyDecorator cur= stack.front();
     stack.pop_front();
-    if (f.visit(cur.get_particle())) {
+    if (f->visit(cur.get_particle())) {
       //std::cerr << "Visiting particle " << cur.get_particle() << std::endl;
       for (int i=cur.get_number_of_children()-1; i>=0; --i) {
         stack.push_back(cur.get_child(i));
@@ -112,14 +113,14 @@
   } while (!stack.empty());
 }
 
-void depth_first_traversal(HierarchyDecorator d, HierarchyVisitor &f)
+void depth_first_traversal(HierarchyDecorator d, HierarchyVisitor *f)
 {
   std::vector<HierarchyDecorator> stack;
   stack.push_back(d);
   do {
     HierarchyDecorator cur= stack.back();
     stack.pop_back();
-    if (f.visit(cur.get_particle())) {
+    if (f->visit(cur.get_particle())) {
       for (int i=cur.get_number_of_children()-1; i>=0; --i) {
         stack.push_back(cur.get_child(i));
       }
@@ -144,16 +145,11 @@
 
 
 Particles
-hierarchy_get_leaves(HierarchyDecorator mhd)
-{
-  Particles out;
-  hierarchy_gather(mhd, MHDMatchingLeaves(),
-                   std::back_inserter(out));
-  return out;
+hierarchy_get_leaves(HierarchyDecorator mhd) {
+  return hierarchy_get(mhd, MHDMatchingLeaves());
 }
 
-BondDecorators hierarchy_get_internal_bonds(HierarchyDecorator mhd)
-{
+BondDecorators hierarchy_get_internal_bonds(HierarchyDecorator mhd) {
   Particles ps= hierarchy_get_all_descendants(mhd);
   std::set<Particle*> sps(ps.begin(), ps.end());
   BondDecorators ret;
@@ -186,12 +182,8 @@
 } // namespace
 
 Particles
-hierarchy_get_all_descendants(HierarchyDecorator mhd)
-{
-  Particles out;
-  hierarchy_gather(mhd, MHDMatchingAll(),
-                   std::back_inserter(out));
-  return out;
+hierarchy_get_all_descendants(HierarchyDecorator mhd) {
+  return  hierarchy_get(mhd, MHDMatchingAll());
 }
 
 } // namespace IMP
Index: kernel/src/decorators/yaml.cpp
===================================================================
--- kernel/src/decorators/yaml.cpp	(revision 0)
+++ kernel/src/decorators/yaml.cpp	(revision 0)
@@ -0,0 +1,272 @@
+#include <IMP/decorators/yaml.h>
+
+namespace IMP
+{
+
+  namespace {
+    struct DefaultWrite {
+      template <class T>
+      std::ostream & operator()(const T t,
+                                std::ostream &out) const {
+        return out << t;
+      }
+    };
+    struct ParticleWrite {
+      std::ostream & operator()(Particle *p,
+                                std::ostream &out) const {
+        return out << p->get_index().get_index();
+      }
+    };
+    struct FloatWrite {
+      std::ostream & operator()(Float t,
+                                std::ostream &out) const {
+        // set format better
+        return out << t;
+      }
+    };
+
+    template <class It, class Write>
+    void write_attributes(std::string indent,
+                          Particle *p,
+                          It b, It e, Write w,
+                          std::ostream &out) {
+      for (It c= b; c != e; ++c) {
+        /** \todo should escape things properly */
+        out << indent << c->get_string() << ": ";
+        w(p->get_value(*c), out) << "\n";
+      }
+    }
+
+    // skips empty lines
+    struct LineStream{
+      std::istream &in;
+      std::vector<std::string> line;
+      bool has_indent(std::string str, unsigned int indent) {
+        for (unsigned int i=0; i< indent; ++i) {
+          if (i== str.size()) return false;
+          if (str[i] != ' ') return false;
+        }
+        return true;
+      }
+      bool not_white(char buf[]) const {
+        for (int i=0; buf[i] != '\0'; ++i) {
+          if (buf[i] != ' ') return true;
+        }
+        return false;
+      }
+
+      LineStream(std::istream &init): in(init){}
+      operator bool() const {return !line.empty() || static_cast<bool>(in);}
+      std::string get_line(unsigned int min_indent) {
+        while (line.empty()) {
+          char buf[3000];
+          in.getline(buf, 3000);
+          if (!in) return std::string();
+          if (buf[0] == '#') continue;
+          if (not_white(buf)) {
+            line.push_back(buf);
+          }
+        }
+        if (has_indent(line.back(), min_indent)) {
+          std::string ret(line.back(), min_indent);
+          line.pop_back();
+          return ret;
+        } else {
+          IMP_LOG(VERBOSE, "Line \"" << line.back() << "\" lacks "
+                  << min_indent << " spaces" << std::endl);
+          return std::string();
+        }
+      }
+      void push_line(std::string s) {
+        if (s.empty()) return;
+        for (unsigned int i=0; i< s.size(); ++i) {
+          if (s[i] != ' ') {
+            line.push_back(s);
+            return;
+          }
+        }
+      }
+    };
+
+    template <class K, class V>
+    struct DefaultRead {
+      void operator()(Particle *p, std::string key, std::string value) const {
+        IMP_LOG(VERBOSE,
+                "Reading values from pair " << key << " "
+                << value << std::endl);
+        K k(key.c_str());
+        std::istringstream iss(value.c_str());
+        V v; 
+        iss >> v;
+        IMP_check(iss, "Error reading value. Got " << v , ValueException);
+        p->set_value(k, v);
+      }
+    };
+
+    struct ParticleRead {
+      void operator()(Particle *p, std::string key, std::string value) const {
+        IMP_LOG(VERBOSE,
+                "Reading values from pair " << key << " "
+                << value << std::endl);
+        ParticleKey k(key.c_str());
+        std::istringstream iss(value.c_str());
+        int i;
+        iss >> i;
+        IMP_check(iss, "Error reading value" , ValueException);
+        Particle *op= p->get_model()->get_particle(ParticleIndex(i));
+        p->set_value(k, op);
+      }
+    };
+
+    int get_next_indent(LineStream &in) {
+      std::string buf= in.get_line(0);
+      if (buf.empty()) return 0;
+      unsigned int i=0;
+      for (; i < buf.size() && buf[i] == ' '; ++i) {
+      }
+      in.push_line(buf);
+      return i;
+    }
+
+
+    template <class Read>
+    void read_attributes(Particle *p, LineStream &in,
+                         int indent,
+                         Read read) {
+      IMP_LOG(VERBOSE, "Reading attributes " << indent << std::endl);
+      int nindent= get_next_indent(in);
+      if (nindent <= indent) return;
+      indent=nindent;
+      IMP_LOG(VERBOSE, "Required indent is " << indent<< std::endl);
+      do {
+        std::string buf = in.get_line(indent);
+        if (buf.empty()) {
+          IMP_LOG(VERBOSE, "Done reading attributes" << std::endl);
+          return;
+        }
+        IMP_check(buf[0] != ' ', "Extra white space on line " 
+                  << buf, InvalidStateException);
+        std::istringstream iss(buf.c_str());
+        char key[2000];
+        iss.get(key, 2000, ':');
+        IMP_check(iss, "no : found in line " << buf,
+                  ValueException);
+        char colon;
+        iss >> colon;
+        IMP_check(colon == ':', "No colon found" << buf,
+                  ValueException);
+
+        char value[2000];
+        iss.getline(value, 2000);
+        IMP_check(iss, "Error reading line " << buf,
+                  ValueException);
+        read(p, key, value);
+
+      } while (true);
+    }
+
+    void read_yaml(Model *m,
+                   LineStream &in,
+                   unsigned int indent) {
+      std::string buf=in.get_line(indent);
+      if (buf.empty()) return;
+      //IMP_LOG(VERBOSE, "Got line " << buf << std::endl);
+      //IMP_check(in, "Error reading particle line from yaml", ValueException);
+      int id;
+      int nread=sscanf(buf.c_str(), "particle: %d", &id);
+      IMP_check(nread==1, "Couldn't read id", InvalidStateException);
+      Particle *p= m->get_particle(id);
+      IMP_LOG(VERBOSE, "Reading particle " << id << std::endl);
+      unsigned int nindent= get_next_indent(in);
+      if (nindent <= indent) return;
+      indent=nindent;
+      while (in) {
+        std::string buf=in.get_line(indent);
+        if (buf.empty()) break;
+        IMP_check(buf[0] != ' ', "Indent error" << buf, InvalidStateException);
+
+        IMP_LOG(VERBOSE, "Looking for attributes in line " << buf << std::endl);
+        std::istringstream iss(buf);
+        std::string type;
+        iss >> type;
+        if (type.compare("float-attributes:")==0) {
+          read_attributes(p, in, indent, DefaultRead<FloatKey, Float>()); 
+        } else if (type.compare("int-attributes:")==0) {
+          read_attributes(p, in, indent, DefaultRead<IntKey, Int>()); 
+        } else if (type.compare("string-attributes:")==0) {
+          read_attributes(p, in, indent, DefaultRead<StringKey, String>()); 
+        } else if (type.compare("particle-attributes:")==0) {
+          read_attributes(p, in, indent, ParticleRead());  
+        } else {
+          break;
+        }
+      }
+      IMP_LOG(VERBOSE, "Done reading particle " << id << std::endl);
+    }
+  }
+
+  static std::string indent_level="  ";
+
+  void write_yaml(Particle *p,
+                  std::ostream &out,
+                  std::string indent) {
+    out << indent << "particle: " << p->get_index().get_index() << "\n";
+    out << indent << indent_level << "float-attributes:\n";
+    write_attributes(indent+indent_level+"  ",
+                     p,
+                     p->float_keys_begin(),
+                     p->float_keys_end(),
+                     FloatWrite(),
+                     out);
+    out << indent << indent_level << "int-attributes:\n";
+    write_attributes(indent+indent_level+"  ",
+                     p,
+                     p->int_keys_begin(),
+                     p->int_keys_end(),
+                     DefaultWrite(),
+                     out);
+    out << indent << indent_level << "string-attributes:\n";
+    write_attributes(indent+indent_level+"  ",
+                     p,
+                     p->string_keys_begin(),
+                     p->string_keys_end(),
+                     DefaultWrite(),
+                     out);
+    out << indent << indent_level << "particle-attributes:\n";
+    write_attributes(indent+indent_level+"  ",
+                     p,
+                     p->particle_keys_begin(),
+                     p->particle_keys_end(),
+                     ParticleWrite(),
+                     out);
+  }
+
+  void write_yaml(Model *m,
+                  std::ostream &out,
+                  std::string indent) {
+    for (Model::ParticleIterator pit= m->particles_begin(); 
+         pit != m->particles_end(); ++pit) {
+      write_yaml(*pit, out, indent);
+    }
+  }
+
+  void read_yaml(std::istream &in,
+                 Model *m) {
+    LineStream r(in);
+    unsigned int nread=0;
+    do {
+      read_yaml(m, r, get_next_indent(r));
+      ++nread;
+    } while (r);
+    IMP_check(nread== m->get_number_of_particles(),
+              "Read wrong number of particles. Model is corrupted. Bye.",
+              ErrorException);
+  }
+
+  void read_yaml(std::string contents,
+                 Model *m) {
+    std::istringstream iss(contents.c_str());
+    IMP_LOG(VERBOSE, "Got string:\n" << contents << std::endl);
+    read_yaml(iss, m);
+  }
+}
Index: kernel/src/decorators/AtomDecorator.cpp
===================================================================
--- kernel/src/decorators/AtomDecorator.cpp	(revision 626)
+++ kernel/src/decorators/AtomDecorator.cpp	(working copy)
@@ -5,13 +5,13 @@
  *
  */
 
+#include "IMP/decorators/AtomDecorator.h"
+#include "IMP/log.h"
+
 #include <sstream>
 #include <vector>
 #include <limits>
 
-#include "IMP/decorators/AtomDecorator.h"
-#include "IMP/log.h"
-
 namespace IMP
 {
 
Index: kernel/src/decorators/ResidueDecorator.cpp
===================================================================
--- kernel/src/decorators/ResidueDecorator.cpp	(revision 626)
+++ kernel/src/decorators/ResidueDecorator.cpp	(working copy)
@@ -12,6 +12,7 @@
 
 #include <sstream>
 #include <vector>
+#include <iterator>
 
 namespace IMP
 {
Index: kernel/src/decorators/MolecularHierarchyDecorator.cpp
===================================================================
--- kernel/src/decorators/MolecularHierarchyDecorator.cpp	(revision 626)
+++ kernel/src/decorators/MolecularHierarchyDecorator.cpp	(working copy)
@@ -54,34 +54,12 @@
                          })
 
 
-namespace
-{
-
-struct MHDMatchingType
-{
-  MHDMatchingType(MolecularHierarchyDecorator::Type t): t_(t){}
-
-  bool operator()(Particle *p) const {
-    MolecularHierarchyDecorator mhd= MolecularHierarchyDecorator::cast(p);
-    if (mhd== MolecularHierarchyDecorator()) {
-      return false;
-    } else {
-      return mhd.get_type()==t_;
-    }
-  }
-
-  MolecularHierarchyDecorator::Type t_;
-};
-
-} // namespace
-
 Particles molecular_hierarchy_get_by_type(MolecularHierarchyDecorator mhd, 
                                           MolecularHierarchyDecorator::Type t)
 {
-  Particles out;
-  hierarchy_gather(mhd, MHDMatchingType(t),
-                   std::back_inserter(out));
-  return out;
+  return hierarchy_get_by_attribute(mhd,
+                                    MolecularHierarchyDecorator::get_type_key(),
+                                    t);
 }
 
 
@@ -156,5 +134,4 @@
   return fd;
 }
 
-
 } // namespace IMP
Index: kernel/src/decorators/SConscript
===================================================================
--- kernel/src/decorators/SConscript	(revision 626)
+++ kernel/src/decorators/SConscript	(working copy)
@@ -1,8 +1,13 @@
 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=[
+        File( 'AtomDecorator.cpp' ),
+        File( 'bond_decorators.cpp' ),
+        File( 'HierarchyDecorator.cpp' ),
+        File( 'MolecularHierarchyDecorator.cpp' ),
+        File( 'NameDecorator.cpp' ),
+        File( 'ResidueDecorator.cpp' ),
+        File( 'XYZDecorator.cpp' ),
+        File( 'XYZRDecorator.cpp' ),
+        File( 'yaml.cpp' ),
+      ]
 Return('files')
Index: kernel/src/decorators/XYZDecorator.cpp
===================================================================
--- kernel/src/decorators/XYZDecorator.cpp	(revision 626)
+++ kernel/src/decorators/XYZDecorator.cpp	(working copy)
@@ -62,16 +62,9 @@
                          key_[2] = FloatKey("z");
                          })
 
-namespace {
-  template <class T>
-  T d(T a, T b){T d=a-b; return d*d;}
-}
-
 Float distance(XYZDecorator a, XYZDecorator b)
 {
-  double d2= d(a.get_x(), b.get_x()) + d(a.get_y(), b.get_y())
-    + d(a.get_z(), b.get_z());
-  return std::sqrt(d2);
+  return (a.get_vector()-b.get_vector()).get_magnitude();
 }
 
 } // namespace IMP
Index: kernel/src/optimizers/states/SConscript
===================================================================
--- kernel/src/optimizers/states/SConscript	(revision 626)
+++ kernel/src/optimizers/states/SConscript	(working copy)
@@ -1,7 +1,7 @@
 Import('env')
-
-files = ['VRMLLogOptimizerState.cpp', 'CMMLogOptimizerState.cpp',
-         'VelocityScalingOptimizerState.cpp']
-
-files = [File(x) for x in files]
+files=[
+        File( 'CMMLogOptimizerState.cpp' ),
+        File( 'VelocityScalingOptimizerState.cpp' ),
+        File( 'VRMLLogOptimizerState.cpp' ),
+      ]
 Return('files')
Index: kernel/src/optimizers/states/VRMLLogOptimizerState.cpp
===================================================================
--- kernel/src/optimizers/states/VRMLLogOptimizerState.cpp	(revision 626)
+++ kernel/src/optimizers/states/VRMLLogOptimizerState.cpp	(working copy)
@@ -10,7 +10,7 @@
 #include <sstream>
 
 #include "IMP/optimizers/states/VRMLLogOptimizerState.h"
-#include "IMP/decorators/XYZDecorator.h"
+#include "IMP/decorators/XYZRDecorator.h"
 
 namespace IMP
 {
@@ -18,7 +18,7 @@
 VRMLLogOptimizerState::VRMLLogOptimizerState(std::string filename,
                                              const Particles &pis) :
     filename_(filename), file_number_(0), call_number_(0),
-    skip_steps_(0)
+    skip_steps_(0), radius_(XYZRDecorator::get_radius_key())
 {
   set_particles(pis);
 }
Index: kernel/src/optimizers/MonteCarlo.cpp
===================================================================
--- kernel/src/optimizers/MonteCarlo.cpp	(revision 626)
+++ kernel/src/optimizers/MonteCarlo.cpp	(working copy)
@@ -39,7 +39,7 @@
 Float MonteCarlo::optimize(unsigned int max_steps)
 {
   IMP_CHECK_OBJECT(this);
-  if (cg_.get() != NULL) {
+  if (cg_) {
     IMP_CHECK_OBJECT(cg_.get());
     IMP_check(cg_->get_model() == get_model(),
                "The model used by the local optimizer does not match "\
@@ -58,7 +58,7 @@
       (*it)->propose_move(probability_);
     }
     Float next_energy;
-    if (cg_.get() != NULL && num_local_steps_!= 0) {
+    if (cg_ && num_local_steps_!= 0) {
       IMP_LOG(VERBOSE,
               "MC Performing local optimization "<< std::flush);
       IMP_CHECK_OBJECT(cg_.get());
Index: kernel/src/optimizers/SConscript
===================================================================
--- kernel/src/optimizers/SConscript	(revision 626)
+++ kernel/src/optimizers/SConscript	(working copy)
@@ -1,11 +1,12 @@
 Import('env')
-
-states_files = SConscript('states/SConscript')
-movers_files = SConscript('movers/SConscript')
-
-files = ['SteepestDescent.cpp', 'ConjugateGradients.cpp',
-         'MolecularDynamics.cpp', 'MonteCarlo.cpp', 'MoverBase.cpp',
-         'BrownianDynamics.cpp'] + states_files + movers_files
-
-files = [File(x) for x in files]
+files=[
+        File( 'BrownianDynamics.cpp' ),
+        File( 'ConjugateGradients.cpp' ),
+        File( 'MolecularDynamics.cpp' ),
+        File( 'MonteCarlo.cpp' ),
+        File( 'MoverBase.cpp' ),
+        File( 'SteepestDescent.cpp' ),
+        SConscript( 'movers/SConscript' ),
+        SConscript( 'states/SConscript' ),
+      ]
 Return('files')
Index: kernel/src/optimizers/movers/SConscript
===================================================================
--- kernel/src/optimizers/movers/SConscript	(revision 626)
+++ kernel/src/optimizers/movers/SConscript	(working copy)
@@ -1,6 +1,6 @@
 Import('env')
-
-files = [ 'BallMover.cpp', 'NormalMover.cpp']
-
-files = [File(x) for x in files]
+files=[
+        File( 'BallMover.cpp' ),
+        File( 'NormalMover.cpp' ),
+      ]
 Return('files')
Index: kernel/src/pair_scores/DistancePairScore.cpp
===================================================================
--- kernel/src/pair_scores/DistancePairScore.cpp	(revision 626)
+++ kernel/src/pair_scores/DistancePairScore.cpp	(working copy)
@@ -22,7 +22,7 @@
 };
 
 Float DistancePairScore::evaluate(Particle *a, Particle *b,
-                                  DerivativeAccumulator *da)
+                                  DerivativeAccumulator *da) const
 {
   return internal::evaluate_distance_pair_score(XYZDecorator(a),
                                                 XYZDecorator(b),
Index: kernel/src/pair_scores/RefineOncePairScore.cpp
===================================================================
--- kernel/src/pair_scores/RefineOncePairScore.cpp	(revision 626)
+++ kernel/src/pair_scores/RefineOncePairScore.cpp	(working copy)
@@ -20,7 +20,7 @@
                                          PairScore *f): r_(r), f_(f) {}
 
 Float RefineOncePairScore::evaluate(Particle *a, Particle *b,
-                                    DerivativeAccumulator *da)
+                                    DerivativeAccumulator *da) const
 {
   Particle* p[2]={a,b};
   Particles ps[2];
Index: kernel/src/pair_scores/SConscript
===================================================================
--- kernel/src/pair_scores/SConscript	(revision 626)
+++ kernel/src/pair_scores/SConscript	(working copy)
@@ -1,8 +1,9 @@
 Import('env')
-
-files = ['DistancePairScore.cpp', 'SphereDistancePairScore.cpp',
-         'RefineOncePairScore.cpp', 'TypedPairScore.cpp',
-         'TransformedDistancePairScore.cpp']
-
-files = [File(x) for x in files]
+files=[
+        File( 'DistancePairScore.cpp' ),
+        File( 'RefineOncePairScore.cpp' ),
+        File( 'SphereDistancePairScore.cpp' ),
+        File( 'TransformedDistancePairScore.cpp' ),
+        File( 'TypedPairScore.cpp' ),
+      ]
 Return('files')
Index: kernel/src/pair_scores/TransformedDistancePairScore.cpp
===================================================================
--- kernel/src/pair_scores/TransformedDistancePairScore.cpp	(revision 626)
+++ kernel/src/pair_scores/TransformedDistancePairScore.cpp	(working copy)
@@ -83,12 +83,13 @@
 };
 
 Float TransformedDistancePairScore::evaluate(Particle *a, Particle *b,
-                                             DerivativeAccumulator *da)
+                                             DerivativeAccumulator *da) const
 {
   TransformParticle tb(r_,ri_, tc_,c_,b);
   IMP_LOG(VERBOSE, "Transformed particle is " 
           << tb.get_coordinate(0) << " " << tb.get_coordinate(1)
-          << " " << tb.get_coordinate(2) << std::endl);
+          << " " << tb.get_coordinate(2) 
+          << " from " << XYZDecorator(b) << std::endl);
   return internal::evaluate_distance_pair_score(XYZDecorator(a),
                                                 tb,
                                                 da, f_.get(), 
@@ -136,6 +137,14 @@
 {
   out << "TransformedDistancePairScore using ";
   f_->show(out);
+  out << "With rotation:\n";
+  out << r_[0] << std::endl;
+  out << r_[1] << std::endl;
+  out << r_[2] << std::endl;
+  out << "and inverse rotation:\n";
+  out << ri_[0] << std::endl;
+  out << ri_[1] << std::endl;
+  out << ri_[2] << std::endl;
 }
 
 } // namespace IMP
Index: kernel/src/pair_scores/SphereDistancePairScore.cpp
===================================================================
--- kernel/src/pair_scores/SphereDistancePairScore.cpp	(revision 626)
+++ kernel/src/pair_scores/SphereDistancePairScore.cpp	(working copy)
@@ -29,7 +29,7 @@
 };
 
 Float SphereDistancePairScore::evaluate(Particle *a, Particle *b,
-                                        DerivativeAccumulator *da)
+                                        DerivativeAccumulator *da) const
 {
   IMP_check(a->has_attribute(radius_), "Particle " << a->get_index() 
             << "missing radius in SphereDistancePairScore",
Index: kernel/src/pair_scores/TypedPairScore.cpp
===================================================================
--- kernel/src/pair_scores/TypedPairScore.cpp	(revision 626)
+++ kernel/src/pair_scores/TypedPairScore.cpp	(working copy)
@@ -11,7 +11,7 @@
 {
 
 Float TypedPairScore::evaluate(Particle *a, Particle *b,
-                               DerivativeAccumulator *da)
+                               DerivativeAccumulator *da) const
 {
   if (!a->has_attribute(typekey_)) {
     set_particle_type(a);
@@ -22,7 +22,7 @@
   Int atype = a->get_value(typekey_);
   Int btype = b->get_value(typekey_);
 
-  ScoreMap::iterator psit =
+  ScoreMap::const_iterator psit =
       score_map_.find(std::pair<Int,Int>(std::min(atype, btype),
                                          std::max(atype, btype)));
   if (psit == score_map_.end()) {
Index: kernel/src/internal/SConscript
===================================================================
--- kernel/src/internal/SConscript	(revision 626)
+++ kernel/src/internal/SConscript	(working copy)
@@ -1,7 +1,9 @@
 Import('env')
-
-files = ['graph_base.cpp', 'ParticleGrid.cpp',
-         'kernel_version_info.cpp', 'constants.cpp', 'bbox_nbl_helpers.cpp']
-
-files = [File(x) for x in files]
+files=[
+        File( 'bbox_nbl_helpers.cpp' ),
+        File( 'constants.cpp' ),
+        File( 'graph_base.cpp' ),
+        File( 'kernel_version_info.cpp' ),
+        File( 'ParticleGrid.cpp' ),
+      ]
 Return('files')
Index: kernel/src/internal/bbox_nbl_helpers.cpp
===================================================================
--- kernel/src/internal/bbox_nbl_helpers.cpp	(revision 626)
+++ kernel/src/internal/bbox_nbl_helpers.cpp	(working copy)
@@ -12,7 +12,9 @@
 /* compile the CGAL code with NDEBUG since it doesn't have the 
    same level of control over errors as IMP
 */
+#ifndef NDEBUG
 #define NDEBUG
+#endif
 #ifdef IMP_USE_CGAL
 #include <CGAL/box_intersection_d.h>
 #include <vector>
@@ -76,7 +78,7 @@
   CGAL::box_intersection_d( boxes0.begin(), boxes0.end(), 
                             boxes1.begin(), boxes1.end(), ap);
 #else
-  IMP_failure( "IMP built without CGAL support.", ErrorException());
+  IMP_failure( "IMP built without CGAL support.", ErrorException);
 #endif
 }
 
@@ -91,7 +93,7 @@
 
   CGAL::box_self_intersection_d( boxes.begin(), boxes.end(), ap);
 #else
-  IMP_failure("IMP built without CGAL support.", ErrorException());
+  IMP_failure("IMP built without CGAL support.", ErrorException);
 #endif
 }
 
Index: kernel/pyext/IMP/test.py
===================================================================
--- kernel/pyext/IMP/test.py	(revision 626)
+++ kernel/pyext/IMP/test.py	(working copy)
@@ -3,6 +3,40 @@
 import IMP
 
 
+def create_particles_in_box( model, num=10,
+                            lb= [0,0,0],
+                            ub= [10,10,10]):
+    """Create a bunch of particles in a box"""
+    lbv=IMP.Vector3D(lb[0],lb[1],lb[2])
+    ubv=IMP.Vector3D(ub[0],ub[1],ub[2])
+    ps= IMP.Particles()
+    for i in range(0,num):
+        p= IMP.Particle()
+        model.add_particle(p)
+        d= IMP.XYZDecorator.create(p)
+        d.randomize_in_box(lbv, ubv)
+        ps.append(p)
+        d.set_coordinates_are_optimized(True)
+    return ps
+
+def randomize_particles(particles, deviation):
+    """Randomize the xyz coordinates of a list of particles"""
+    for p in particles:
+        d= IMP.XYZDecorator.cast(p)
+        d.set_x(random.uniform(-deviation, deviation))
+        d.set_y(random.uniform(-deviation, deviation))
+        d.set_z(random.uniform(-deviation, deviation))
+
+def create_chain(ps, length=1, stiffness=1):
+    bds= []
+    IMP.BondedDecorator.create(ps[0])
+    for i in range(1,len(ps)):
+        ba= IMP.BondedDecorator.cast(ps[i-1])
+        bb= IMP.BondedDecorator.create(ps[i])
+        bds.append(IMP.custom_bond(ba, bb, length, stiffness))
+    return bds
+
+
 class TestCase(unittest.TestCase):
     """Super class for IMP test cases"""
 
@@ -10,6 +44,8 @@
         self.__check_level = IMP.get_check_level()
         # Turn on expensive runtime checks while running the test suite:
         IMP.set_check_level(IMP.EXPENSIVE)
+        # should go in library init code
+        IMP.set_print_exception_messages(False)
 
     def tearDown(self):
         # Restore original check level
@@ -33,14 +69,6 @@
         p.add_attribute(IMP.FloatKey("z"), z, True)
         return p
 
-    def randomize_particles(self, particles, deviation):
-        """Randomize the xyz coordinates of a list of particles"""
-        for p in particles:
-            d= IMP.XYZDecorator.cast(p)
-            d.set_x(random.uniform(-deviation, deviation))
-            d.set_y(random.uniform(-deviation, deviation))
-            d.set_z(random.uniform(-deviation, deviation))
-
     def particle_distance(self, p1, p2):
         """Return distance between two given particles"""
         xkey = IMP.FloatKey("x")
@@ -73,30 +101,7 @@
                 fmin, vmin = f, v
         self.assertInTolerance(fmin, expected_fmin, step)
 
-    def create_particles_in_box(self, model, num=10,
-                                lb= [0,0,0],
-                                ub= [10,10,10]):
-        """Create a bunch of particles in a box"""
-        lbv=IMP.Vector3D(lb[0],lb[1],lb[2])
-        ubv=IMP.Vector3D(ub[0],ub[1],ub[2])
-        ps= IMP.Particles()
-        for i in range(0,num):
-            p= IMP.Particle()
-            model.add_particle(p)
-            d= IMP.XYZDecorator.create(p)
-            d.randomize_in_box(lbv, ubv)
-            ps.append(p)
-            d.set_coordinates_are_optimized(True)
-        return ps
 
-    def create_chain(self, ps, length=1, stiffness=1):
-        bds= []
-        IMP.BondedDecorator.create(ps[0])
-        for i in range(1,len(ps)):
-            ba= IMP.BondedDecorator.cast(ps[i-1])
-            bb= IMP.BondedDecorator.create(ps[i])
-            bds.append(IMP.custom_bond(ba, bb, length, stiffness))
-        return bds
 
 
 class ConstPairScore(IMP.PairScore):
Index: kernel/pyext/IMP/modeller_intf.py
===================================================================
--- kernel/pyext/IMP/modeller_intf.py	(revision 626)
+++ kernel/pyext/IMP/modeller_intf.py	(working copy)
@@ -241,7 +241,7 @@
 
 def _LinearGenerator(parameters, modalities):
     (scale,) = parameters
-    return IMP.Linear(scale)
+    return IMP.Linear(0,scale)
 
 def _SplineGenerator(parameters, modalities):
     (open, low, high, delta, lowderiv, highderiv) = parameters[:6]
Index: kernel/pyext/IMP.i
===================================================================
--- kernel/pyext/IMP.i	(revision 626)
+++ kernel/pyext/IMP.i	(working copy)
@@ -36,13 +36,16 @@
   %pythonprepend Optimizer::add_optimizer_state %{
         args[1].thisown=0
   %}
+  %pythonprepend Optimizer::set_model %{
+        args[1].thisown=0
+  %}
   %pythonprepend RestraintSet::add_restraint %{
         args[1].thisown=0
   %}
   %pythonprepend NonbondedListScoreState::add_bonded_list %{
         args[1].thisown=0
   %}
-  %pythonprepend TunnelRestraint::TunnelRestraint %{
+  %pythonprepend TunnelSingletonScore::TunnelSingletonScore %{
         args[0].thisown=0
   %}
   %pythonprepend DistanceRestraint::DistanceRestraint %{
@@ -75,6 +78,9 @@
   %pythonprepend PairChainRestraint::PairChainRestraint %{
         args[0].thisown=0
   %}
+  %pythonprepend PairRestraint::PairRestraint %{
+        args[0].thisown=0
+  %}
   %pythonprepend ConnectivityRestraint::ConnectivityRestraint %{
         args[0].thisown=0
   %}
@@ -113,7 +119,7 @@
         args[1].thisown=0
   %}
   %pythonprepend TypedPairScore::set_pair_score %{
-        args[0].thisown=0
+        args[1].thisown=0
   %}
   %pythonprepend Particle::get_value %{
         check_particle(args[0], args[1])
@@ -148,7 +154,10 @@
 %feature("ref")   Particle "$this->ref();"
 %feature("unref") Particle "$this->unref(); if (! $this->get_has_ref()) delete $this;"
 
+/*%feature("ref")   Model "$this->ref();"
+%feature("unref") Model "$this->unref(); if (! $this->get_has_ref()) delete $this;"*/
 
+
 /* Don't wrap internal functions */
 %ignore IMP::internal::evaluate_distance_pair_score;
 %ignore IMP::internal::check_particles_active;
@@ -171,14 +180,6 @@
 %include "IMP/base_types.h"
 %include "IMP/VersionInfo.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/unary_functions/WormLikeChain.h"
 %include "IMP/DerivativeAccumulator.h"
 %include "IMP/Restraint.h"
 %include "IMP/ScoreState.h"
@@ -190,6 +191,7 @@
 %include "IMP/SingletonScore.h"
 %include "IMP/TripletScore.h"
 %include "IMP/Particle.h"
+%include "IMP/ParticleRefiner.h"
 %include "Vector3D.i"
 %include "IMP/DecoratorBase.h"
 %include "IMP/decorators/bond_decorators.h"
@@ -198,8 +200,9 @@
 %include "IMP/decorators/NameDecorator.h"
 %include "IMP/decorators/ResidueDecorator.h"
 %include "IMP/decorators/XYZDecorator.h"
+%include "IMP/decorators/XYZRDecorator.h"
 %include "IMP/decorators/AtomDecorator.h"
-%include "IMP/ParticleRefiner.h"
+%include "IMP/decorators/yaml.h"
 %include "IMP/particle_refiners/BondCoverParticleRefiner.h"
 %include "IMP/particle_refiners/ChildrenParticleRefiner.h"
 %include "IMP/Optimizer.h"
@@ -224,16 +227,18 @@
 %include "IMP/particle_refiners/ChildrenParticleRefiner.h"
 %include "IMP/singleton_scores/DistanceToSingletonScore.h"
 %include "IMP/singleton_scores/AttributeSingletonScore.h"
-%include "IMP/triplet_scores/AngleTripletScore.h"
+%include "IMP/singleton_scores/TunnelSingletonScore.h"
 %include "IMP/score_states/BondedListScoreState.h"
 %include "IMP/score_states/MaxChangeScoreState.h"
 %include "IMP/score_states/NonbondedListScoreState.h"
 %include "IMP/score_states/AllNonbondedListScoreState.h"
+%include "IMP/score_states/ManualBondDecoratorListScoreState.h"
 %include "IMP/score_states/BondDecoratorListScoreState.h"
 %include "IMP/score_states/BipartiteNonbondedListScoreState.h"
 %include "IMP/score_states/GravityCenterScoreState.h"
 %include "IMP/score_states/CoverBondsScoreState.h"
 %include "IMP/restraints/AngleRestraint.h"
+%include "IMP/restraints/PairRestraint.h"
 %include "IMP/restraints/BondDecoratorRestraint.h"
 %include "IMP/restraints/ConnectivityRestraint.h"
 %include "IMP/restraints/ConstantRestraint.h"
@@ -245,7 +250,15 @@
 %include "IMP/restraints/RestraintSet.h"
 %include "IMP/restraints/SingletonListRestraint.h"
 %include "IMP/restraints/TripletChainRestraint.h"
-%include "IMP/restraints/TunnelRestraint.h"
+%include "IMP/triplet_scores/AngleTripletScore.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/unary_functions/WormLikeChain.h"
 
 namespace IMP {
   %template(ParticleIndex) Index<ParticleTag>;