[IMP-dev] previously broken patches

27 Aug 2008

In order to avoid going back and forth a few times, here is a patch with 
all the changed files
The contained patches are:
kernel/include/IMP/optimizers/movers/BallMover.h
kernel/src/IMP/optimizers/movers/BallMover.cpp
Somewhere between my copy of IMP and the general copy of IMP BallMover 
became completely broken due to a missing &, showing why returning by 
reference arguments is problematic :-)
----
kernel/pyext/IMP.i
For uniformity we should make all containers look the same. And Mover 
has virtual methods so they
should get directors.
-----
kernel/include/IMP/optimizers/MonteCarlo.h
kernel/include/IMP/optimizers/Mover.h
Movers should probably be declared in Mover.h.
-----
kernel/include/IMP/utility.h
kernel/include/IMP/internal/AttributeTable.h
isnan function to make it clearer what we are testing for when we test 
things
----
kernel/test/states/test_nonbonded_list.py
kernel/test/states/test_cover_bonds.py
kernel/pyext/IMP/test.py
kernel/src/Vector3D.cpp
kernel/src/decorators/XYZDecorator.cpp
kernel/src/SConscript
kernel/doc/examples/chain.py
kernel/include/IMP/decorators/macros.h
kernel/include/IMP/decorators/XYZDecorator.h
Move the random functions from XYZDecorator to Vector3D.
----
kernel/src/optimizers/MonteCarlo.cpp
don't try to get a null pointer
-----
kernel/test/pair_scores/test_transform.py
kernel/src/pair_scores/TransformedDistancePairScore.cpp
kernel/src/pair_scores/DistancePairScore.cpp
kernel/include/IMP/pair_scores/TransformedDistancePairScore.h
kernel/include/IMP/internal/evaluate_distance_pair_score.h
Clean up the evaluate distance pair score functions.
Index: kernel/test/states/test_nonbonded_list.py
===================================================================

--- kernel/test/states/test_nonbonded_list.py	(revision 681)
+++ kernel/test/states/test_nonbonded_list.py	(working copy)
@@ -86,7 +86,8 @@
         for i in range(0,10):
             score= m.evaluate(False)
             for d in ds:
-                d.randomize_in_sphere(d.get_vector(), 2.0)
+                d.set_coordinates(IMP.random_vector_in_sphere(d.get_vector(),
+                                                              2.0))
def do_test_bl(self, ss):
         """Test the bond decorator list"""
@@ -148,8 +149,8 @@
         print "Index is " +str(p.get_index().get_index())
         d=IMP.XYZDecorator.create(p)
         d.set_coordinates_are_optimized(True)
-        d.randomize_in_box(IMP.Vector3D(0,0,0),
-                           IMP.Vector3D(10,10,10));
+        d.set_coordinates(IMP.random_vector_in_box(IMP.Vector3D(0,0,0),
+                                                   IMP.Vector3D(10,10,10)))
         nps= IMP.Particles([p])
         s.add_particles(nps)
         score= m.evaluate(False)
@@ -193,10 +194,12 @@
             score= m.evaluate(False)
             for p in ps0:
                 d= IMP.XYZDecorator.cast(p)
-                d.randomize_in_sphere(d.get_vector(), 1)
+                d.set_coordinates(IMP.random_vector_in_sphere(d.get_vector(),
+                                                              1))
             for p in ps1:
                 d= IMP.XYZDecorator.cast(p)
-                d.randomize_in_sphere(d.get_vector(), 1)
+                d.set_coordinates(IMP.random_vector_in_sphere(d.get_vector(),
+                                                              1))
def do_test_spheres(self, ss):
Index: kernel/test/states/test_cover_bonds.py
===================================================================
--- kernel/test/states/test_cover_bonds.py	(revision 681)
+++ kernel/test/states/test_cover_bonds.py	(working copy)
@@ -15,8 +15,8 @@
             p= IMP.Particle()
             m.add_particle(p)
             d= IMP.XYZDecorator.create(p)
-            d.randomize_in_box(IMP.Vector3D(0,0,0),
-                               IMP.Vector3D(10,10,10))
+            d.set_coordinates(IMP.random_vector_in_box(IMP.Vector3D(0,0,0),
+                                                       IMP.Vector3D(10,10,10)))
             ps.append(p)
         bds= []
         bb= IMP.BondedDecorator.create(ps[0])
Index: kernel/test/pair_scores/test_transform.py
===================================================================
--- kernel/test/pair_scores/test_transform.py	(revision 681)
+++ kernel/test/pair_scores/test_transform.py	(working copy)
@@ -1,6 +1,8 @@
 import unittest
 import IMP.utils
 import IMP.test, IMP
+import math
+import random
class DistanceTests(IMP.test.TestCase):
     """Test the symmetry restraint"""
@@ -15,7 +17,7 @@
         m.add_particle(p1)
         d1= IMP.XYZDecorator.create(p1)
         tps= IMP.TransformedDistancePairScore(IMP.Harmonic(0,1))
-        tps.set_translation(0,1,0)
+        tps.set_translation(IMP.Vector3D(0,1,0))
         tps.set_rotation(1, 0, 0,
                          0, 1, 0,
                          0, 0, 1)
@@ -23,10 +25,10 @@
         d1.set_coordinates(IMP.Vector3D(2,2,4))
         self.assertEqual(tps.evaluate(p0, p1, None), 0)
         self.assert_(tps.evaluate(p1, p0, None) != 0)
-        tps.set_center(10,13,4)
+        tps.set_center(IMP.Vector3D(10,13,4))
         self.assertEqual(tps.evaluate(p0, p1, None), 0)
         self.assert_(tps.evaluate(p1, p0, None) != 0)
-        tps.set_center(0,0,0)
+        tps.set_center(IMP.Vector3D(0,0,0))
         print "test rotation"
         tps.set_rotation(0, 0,-1,
                          0, 1, 0,
@@ -34,7 +36,7 @@
         d1.set_coordinates(IMP.Vector3D(4, 2, -2))
         self.assertEqual(tps.evaluate(p0, p1, None), 0)
         self.assert_(tps.evaluate(p1, p0, None) != 0)
-        tps.set_translation(0,0,0)
+        tps.set_translation(IMP.Vector3D(0,0,0))
         tps.set_rotation(0,-1, 0,
                          1, 0, 0,
                          0, 0, 1)
@@ -71,26 +73,121 @@
         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)
+        tps.set_center(cv)
         tps.set_rotation(1, 0, 0,
                          0, 0,-1,
                          0, 1, 0)
         pr= IMP.PairListRestraint(tps)
         pr.add_particle_pair(IMP.ParticlePair(p0, p1))
         m.add_restraint(pr)
-        cg= IMP.ConjugateGradients()
+        cg= IMP.SteepestDescent()
         cg.set_model(m)
+        IMP.set_log_level(IMP.SILENT)
         cg.optimize(100)
+        IMP.set_log_level(IMP.VERBOSE)
         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)
-        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))
+        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 "itrans: "
+        tps.get_transformed(d1.get_vector()).show()
+        print "trans: "
+        fvt.show()
+        self.assertInTolerance(fvt[0], d0.get_coordinate(0), .1)
+        self.assertInTolerance(fvt[1], d0.get_coordinate(1), .1)
+        self.assertInTolerance(fvt[2], d0.get_coordinate(2), .1)
+    def _setup_system(self):
+        random.seed()
+        IMP.set_log_level(IMP.VERBOSE)
+        m= IMP.Model()
+        p0= IMP.Particle()
+        m.add_particle(p0)
+        d0= IMP.XYZDecorator.create(p0)
+        p1= IMP.Particle()
+        m.add_particle(p1)
+        d1= IMP.XYZDecorator.create(p1)
+        d0.set_coordinates(IMP.random_vector_in_box(IMP.Vector3D(-10, -10, -10),
+                                                    IMP.Vector3D(10,10,10)))
+        d1.set_coordinates(IMP.random_vector_in_box(IMP.Vector3D(-10, -10, -10),
+                                                    IMP.Vector3D(10,10,10)))
+        d0.set_coordinates_are_optimized(True)
+        d1.set_coordinates_are_optimized(True)
+        tps= IMP.TransformedDistancePairScore(IMP.Harmonic(0,1))
+        cv= IMP.random_vector_in_box(IMP.Vector3D(-10, -10, -10),
+                                     IMP.Vector3D(10,10,10))
+        tv= IMP.random_vector_in_box(IMP.Vector3D(-10, -10, -10),
+                                     IMP.Vector3D(10,10,10))
+        tps.set_translation(tv)
+        tps.set_center(cv)
+        tps.set_rotation(IMP.random_vector_on_sphere(), 2.0*math.pi*random.random())
+        pr= IMP.PairListRestraint(tps)
+        pr.add_particle_pair(IMP.ParticlePair(p0, p1))
+        m.add_restraint(pr)
+        return (m, tps, d0, d1)
+
+
+    def test_transoformish(self):
+        """Test optimizing with monte carlo to test the non-deriv part"""
+        IMP.set_log_level(IMP.VERBOSE)
+        (m, tps, d0, d1)= self._setup_system()
+
+        cg= IMP.MonteCarlo()
+        cg.set_model(m)
+        cg.set_temperature(0.0001)
+        bm= IMP.BallMover(IMP.XYZDecorator.get_xyz_keys(),
+                          1.0, IMP.Particles(1, d0.get_particle()))
+        cg.add_mover(bm)
+        IMP.set_log_level(IMP.SILENT)
+        cg.optimize(2000)
+        IMP.set_log_level(IMP.VERBOSE)
+        d0.get_vector().show(); print
+        d1.get_vector().show(); print
+        self.assertInTolerance((d0.get_vector()
+                                -tps.get_transformed(d1.get_vector())).get_magnitude(), 0,
+                               1)
+
+    def test_transoformish2(self):
+        """Test optimizing with steepest descent for testing derivatives part"""
+        IMP.set_log_level(IMP.VERBOSE)
+        (m, tps, d0, d1)= self._setup_system()
+        tps.show()
+        cg= IMP.SteepestDescent()
+        cg.set_model(m)
+        IMP.set_log_level(IMP.SILENT)
+        cg.optimize(2000)
+        IMP.set_log_level(IMP.VERBOSE)
+        print "p0 is "
+        d0.get_vector().show(); print
+        print "p1 is "
+        d1.get_vector().show(); print
+        print "target is"
+        tps.get_transformed(d1.get_vector()).show(); print
+        self.assertInTolerance((d0.get_vector()
+                                -tps.get_transformed(d1.get_vector())).get_magnitude(), 0,
+                               1)
+
+
+    def test_set_rot(self):
+        """Testing setting the rotation in tps"""
+        (m, tps, d0, d1)= self._setup_system()
+        v= IMP.random_vector_on_sphere()
+        a= random.random()*2.0* math.pi
+        tps.set_rotation(v, a)
+        tv= tps.get_transformed(d0.get_vector())
+        tps.set_rotation(-v, -a)
+        tv2= tps.get_transformed(d0.get_vector())
+        self.assertInTolerance((tv-tv2).get_magnitude(), 0, .1);
+
+
+
+
+
 if __name__ == '__main__':
     unittest.main()
Index: kernel/include/IMP/optimizers/MonteCarlo.h
===================================================================
--- kernel/include/IMP/optimizers/MonteCarlo.h	(revision 681)
+++ kernel/include/IMP/optimizers/MonteCarlo.h	(working copy)
@@ -15,8 +15,6 @@
 namespace IMP
 {
-typedef std::vector<Mover*> Movers;
-
 //! A Monte Carlo optimizer.
 /** The optimizer uses a set of Mover objects to propose steps. Currently
     each Mover is called at each Monte Carlo iteration. This may change in
Index: kernel/include/IMP/optimizers/movers/BallMover.h
===================================================================
--- kernel/include/IMP/optimizers/movers/BallMover.h	(revision 681)
+++ kernel/include/IMP/optimizers/movers/BallMover.h	(working copy)
@@ -42,6 +42,8 @@
   Float get_radius() const {
     return radius_;
   }
+
+  void show(std::ostream &out= std::cout) const;
 protected:
   /** \internal */
   void generate_move(float a);
Index: kernel/include/IMP/optimizers/Mover.h
===================================================================
--- kernel/include/IMP/optimizers/Mover.h	(revision 681)
+++ kernel/include/IMP/optimizers/Mover.h	(working copy)
@@ -76,6 +76,8 @@
IMP_OUTPUT_OPERATOR(Mover);
+typedef std::vector<Mover*> Movers;
+
 } // namespace IMP
#endif  /* __IMP_MOVER_H */
Index: kernel/include/IMP/decorators/XYZDecorator.h
===================================================================
--- kernel/include/IMP/decorators/XYZDecorator.h	(revision 681)
+++ 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
 {
@@ -104,23 +101,23 @@
   /** Somewhat suspect based on wanting a Point/Vector differentiation
       but we don't have points */
   Vector3D get_vector() const {
-    return Vector3D(get_x(), get_y(), get_z());
+       return Vector3D(get_x(), get_y(), get_z());
   }
+  //! Get the vector of derivatives.
+  /** Somewhat suspect based on wanting a Point/Vector differentiation
+      but we don't have points */
+  Vector3D get_derivative_vector() const {
+    return Vector3D(get_coordinate_derivative(0),
+                    get_coordinate_derivative(1),
+                    get_coordinate_derivative(2));
+  }
+
   //! 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);
-
-  //! Generate random coordinates in a box defined by the vectors
-  void randomize_in_box(const Vector3D &lower_corner,
-                        const Vector3D &upper_corner);
 protected:
   static FloatKey get_coordinate_key(unsigned int i) {
     IMP_check(i <3, "Out of range coordinate",
Index: kernel/include/IMP/decorators/macros.h
===================================================================
--- kernel/include/IMP/decorators/macros.h	(revision 681)
+++ kernel/include/IMP/decorators/macros.h	(working copy)
@@ -44,7 +44,7 @@
     }                                                                   \
     add_required;                                                       \
   }                                                                     \
-  friend class DecoratorBase;                                           \
+  friend class ::IMP::DecoratorBase;                                    \
 public:                                                                 \
  typedef Name This;                                                     \
  /** \short The default constructor. This is used as a null value */    \
@@ -58,14 +58,14 @@
  }                                                                      \
  /** Add the necessary attributes to p and return a decorator. */       \
  static Name create(::IMP::Particle *p) {                               \
-   return IMP::DecoratorBase::create<Name>(p);                          \
+   return ::IMP::DecoratorBase::create<Name>(p);                        \
  }                                                                      \
  /** Check that p has the necessary attributes and return a decorator.  \
      \throws InvalidStateException if some required attributes are      \
      missing                                                            \
  */                                                                     \
  static Name cast(::IMP::Particle *p) {                                 \
-   return IMP::DecoratorBase::cast<Name>(p);                            \
+   return ::IMP::DecoratorBase::cast<Name>(p);                          \
  }                                                                      \
  static bool is_instance_of(::IMP::Particle *p) {                       \
    decorator_initialize_static_data();                                  \
@@ -262,6 +262,17 @@
   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;                                        \
+  }
+
+
 #define IMP_ATOM_TYPE_INDEX 8974343
 #define IMP_RESIDUE_TYPE_INDEX 90784334
Index: kernel/include/IMP/pair_scores/TransformedDistancePairScore.h
===================================================================
--- kernel/include/IMP/pair_scores/TransformedDistancePairScore.h	(revision 681)
+++ kernel/include/IMP/pair_scores/TransformedDistancePairScore.h	(working copy)
@@ -25,13 +25,28 @@
The second particle, x, is transformed as R*(x-center)+ translation+center
+    \note This score is asymetric, that is which point is passed first or
+    second makes a difference. If you don't know which direction the symetry
+    goes make sure you try both directions.
+
+    \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
 {
+  class TransformParticle;
+  friend class TransformParticle;
   Pointer<UnaryFunction> f_;
   Vector3D tc_, c_;
   Vector3D r_[3], ri_[3];
+
+  Float get_transformed(const Vector3D &v, unsigned int i) const ;
+  // Transform a derivative back
+  Vector3D get_back_rotated(const Vector3D &v) const ;
 public:
   TransformedDistancePairScore(UnaryFunction *f);
   virtual ~TransformedDistancePairScore(){}
@@ -39,11 +54,23 @@
                          DerivativeAccumulator *da) const;
   virtual void show(std::ostream &out=std::cout) const;
+  //! Set the rotation from a rotation matrix
+  /** The matrix is multiplied by a vector from the right.
+   */
   void set_rotation(float r00, float r01, float r02,
                     float r10, float r11, float r12,
                     float r20, float r21, float r22);
-  void set_translation(float t0, float t1, float t2);
-  void set_center(float t0, float t1, float t2);
+  //! Set the rotation from an axis and an angle (radians)
+  /** The axis should be close to a unit vector.
+   */
+  void set_rotation(const Vector3D &axis, Float angle);
+
+  void set_translation(const Vector3D &t);
+  //! Set the rotation center
+  void set_center(const Vector3D &t);
+
+  //! Apply the current transformation to the given vector
+  Vector3D get_transformed(const Vector3D &v) const ;
 };
} // namespace IMP
Index: kernel/include/IMP/internal/AttributeTable.h
===================================================================
--- kernel/include/IMP/internal/AttributeTable.h	(revision 681)
+++ kernel/include/IMP/internal/AttributeTable.h	(working copy)
@@ -47,7 +47,7 @@
   }
   static bool get_is_valid(Float f) {
     if (std::numeric_limits<Float>::has_quiet_NaN) {
-      return f==f;
+      return !is_nan(f);
     } else {
       return f != get_invalid();
     }
Index: kernel/include/IMP/internal/evaluate_distance_pair_score.h
===================================================================
--- kernel/include/IMP/internal/evaluate_distance_pair_score.h	(revision 681)
+++ kernel/include/IMP/internal/evaluate_distance_pair_score.h	(working copy)
@@ -17,41 +17,51 @@
 namespace internal
 {
+template <class SD>
+Float compute_distance_pair_score(const Vector3D &delta,
+                                  const UnaryFunction *f,
+                                  Vector3D *d,
+                                  SD sd) {
+  static const Float MIN_DISTANCE = .00001;
+  Float distance= delta.get_magnitude();
+  Float shifted_distance = sd(distance);
+
+  // if needed, calculate the partial derivatives of the scores with respect
+  // to the particle attributes
+  // avoid division by zero if the distance is too small
+  Float score, deriv;
+  if (d && distance >= MIN_DISTANCE) {
+    boost::tie(score, deriv) = f->evaluate_with_derivative(shifted_distance);
+
+    *d= delta/distance *deriv;
+  } else {
+    if (d) *d= Vector3D(0,0,0);
+    // calculate the score based on the distance feature
+    score = f->evaluate(shifted_distance);
+  }
+  return score;
+}
+
+
 template <class W0, class W1, class SD>
 Float evaluate_distance_pair_score(W0 d0, W1 d1,
                                    DerivativeAccumulator *da,
-                                   UnaryFunction *f, SD sd)
+                                   const UnaryFunction *f, SD sd)
 {
-  static const Float MIN_DISTANCE = .00001;
   IMP_CHECK_OBJECT(f);
-  Float d2 = 0;
   Vector3D delta;
-  Float score;
for (int i = 0; i < 3; ++i) {
     delta[i] = d0.get_coordinate(i) - d1.get_coordinate(i);
-    d2 += square(delta[i]);
   }
-  Float distance = std::sqrt(d2);
+  Vector3D d;
+  Float score= compute_distance_pair_score(delta, f, (da? &d : NULL), sd);
-  Float shifted_distance = sd(distance); //scale*(distance - offset);
-
-  // if needed, calculate the partial derivatives of the scores with respect
-  // to the particle attributes
-  // avoid division by zero if the distance is too small
-  if (da && distance >= MIN_DISTANCE) {
-    Float deriv;
-
-    boost::tie(score, deriv) = f->evaluate_with_derivative(shifted_distance);
-
-    Vector3D d= delta/distance *deriv;
+  if (da) {
     d0.add_to_coordinates_derivative(d, *da);
     d1.add_to_coordinates_derivative(-d, *da);
-  } else {
-    // calculate the score based on the distance feature
-    score = f->evaluate(shifted_distance);
   }
return score;
Index: kernel/include/IMP/Vector3D.h
===================================================================
--- kernel/include/IMP/Vector3D.h	(revision 681)
+++ kernel/include/IMP/Vector3D.h	(working copy)
@@ -11,6 +11,7 @@
 #include "IMP_config.h"
 #include "base_types.h"
 #include "macros.h"
+#include "exception.h"
#include <cmath>
@@ -22,7 +23,11 @@
  */
 class IMPDLLEXPORT Vector3D
 {
+  bool is_default() const {return false;}
 public:
+  // public for swig
+  typedef Vector3D This;
+
   //! Initialize the vector from separate x,y,z values.
   Vector3D(Float x, Float y, Float z) {
     vec_[0] = x;
@@ -33,6 +38,8 @@
   //! Default constructor
   Vector3D() {}
+  IMP_COMPARISONS_3(vec_[0], vec_[1], vec_[2]);
+
   //! \return A single component of this vector (0-2).
   Float operator[](unsigned int i) const {
     IMP_assert(i < 3, "Invalid component of vector requested");
@@ -149,14 +156,6 @@
         << operator[](2) << ")";
   }
-  bool operator<(const Vector3D &o) const {
-    for (unsigned int i=0; i< 3; ++i) {
-      if (operator[](i) < o[i]) return true;
-      else if (operator[](i) > o[i]) return false;
-    }
-    return false;
-  }
-
 private:
   Float vec_[3];
 };
@@ -171,6 +170,53 @@
                   o[2]*s);
 }
+//! Generate a random vector in a box with uniform density
+/**
+   \ingroup uncommitted
+ */
+IMPDLLEXPORT Vector3D
+random_vector_in_box(const Vector3D &lb=Vector3D(0,0,0),
+                     const Vector3D &ub=Vector3D(10,10,10));
+
+//! Generate a random vector in a sphere with uniform density
+/**
+   \ingroup uncommitted
+ */
+IMPDLLEXPORT Vector3D
+random_vector_in_sphere(const Vector3D &center=Vector3D(0,0,0),
+                        Float radius=1);
+
+//! Generate a random vector on a sphere with uniform density
+/**
+   \ingroup uncommitted
+ */
+IMPDLLEXPORT Vector3D
+random_vector_on_sphere(const Vector3D &center=Vector3D(0,0,0),
+                        Float radius=1);
+
+
+struct SpacesIO {
+  const Vector3D &v_;
+  SpacesIO(const Vector3D &v): v_(v){}
+};
+
+
+inline std::ostream &operator<<(std::ostream &out, const SpacesIO &s) {
+  out << s.v_[0] << " " << s.v_[1] << " " << s.v_[2];
+  return out;
+}
+
+struct CommasIO {
+  const Vector3D &v_;
+  CommasIO(const Vector3D &v): v_(v){}
+};
+
+
+inline std::ostream &operator<<(std::ostream &out, const CommasIO &s) {
+  out << s.v_[0] << ", " << s.v_[1] << ", " << s.v_[2];
+  return out;
+}
+
 } // namespace IMP
#endif  /* __IMP_VECTOR_3D_H */
Index: kernel/include/IMP/utility.h
===================================================================
--- kernel/include/IMP/utility.h	(revision 681)
+++ kernel/include/IMP/utility.h	(working copy)
@@ -8,8 +8,16 @@
 #ifndef __IMP_UTILITY_H
 #define __IMP_UTILITY_H
+#include "base_types.h"
 #include "macros.h"
+#ifdef __GNUC__
+/** \todo Use an sconscript test for the isnan function directly
+ */
+#include <cmath>
+#endif
+
+
 namespace IMP
 {
@@ -20,6 +28,18 @@
   return t*t;
 }
+//! put in check for G++
+/** This should not be called isnan as C99 says isnan is a macro. It's
+    all a bit messed up, but it is nice to have something grepable.
+ */
+inline bool is_nan(Float a) {
+#ifdef __GNUC__
+  return std::isnan(a);
+#else
+  return a != a;
+#endif
+}
+
 } // namespace IMP
#endif  /* __IMP_UTILITY_H */
Index: kernel/doc/examples/chain.py
===================================================================
--- kernel/doc/examples/chain.py	(revision 681)
+++ kernel/doc/examples/chain.py	(working copy)
@@ -17,8 +17,8 @@
     p= IMP.Particle()
     pi= m.add_particle(p)
     d= IMP.XYZDecorator.create(p)
-    d.randomize_in_box(IMP.Vector3D(0,0,0),
-                       IMP.Vector3D(10,10,10))
+    d.set_coordinates(IMP.random_vector_in_box(IMP.Vector3D(0,0,0),
+                                               IMP.Vector3D(10,10,10)))
     d.set_coordinates_are_optimized(True)
     p.add_attribute(rk, radius, False)
     chain.append(p)
Index: kernel/src/SConscript
===================================================================
--- kernel/src/SConscript	(revision 681)
+++ kernel/src/SConscript	(working copy)
@@ -22,7 +22,8 @@
 files = ['base_types.cpp', 'Model.cpp',
          'Particle.cpp', 'ScoreState.cpp', 'Object.cpp',
          'OptimizerState.cpp', 'Log.cpp', 'Restraint.cpp', 'Optimizer.cpp',
-         'random.cpp', 'Key.cpp', 'exception.cpp', 'ParticleRefiner.cpp'
+         'random.cpp', 'Key.cpp', 'exception.cpp', 'ParticleRefiner.cpp',
+         'Vector3D.cpp'
         ] + decorators_files + restraints_files + optimizers_files \
           + unary_functions_files + pair_scores_files + singleton_scores_files \
           + triplet_scores_files + score_states_files + internal_files \
Index: kernel/src/optimizers/MonteCarlo.cpp
===================================================================
--- kernel/src/optimizers/MonteCarlo.cpp	(revision 681)
+++ 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/movers/BallMover.cpp
===================================================================
--- kernel/src/optimizers/movers/BallMover.cpp	(revision 681)
+++ kernel/src/optimizers/movers/BallMover.cpp	(working copy)
@@ -16,7 +16,7 @@
static void random_point_in_sphere(unsigned int D,
                                    Float radius,
-                                   std::vector<Float> v)
+                                   std::vector<Float>& v)
 {
   IMP_assert(radius > 0, "No volume there");
   ::boost::uniform_real<> rand(-radius, radius);
@@ -66,7 +66,21 @@
     for (unsigned int j = 0; j < get_number_of_float_keys(); ++j) {
       propose_value(i, j, npos[j]);
     }
+    IMP_LOG(VERBOSE, "Proposing move to ");
+    for (unsigned int j=0; j< npos.size(); ++j) {
+      IMP_LOG(VERBOSE, npos[j] << " ");
+    }
+    IMP_LOG(VERBOSE, " for particle " << get_particle(i)->get_index()
+            << std::endl);
   }
 }
+
+
+void BallMover::show(std::ostream &out) const
+{
+  out << "Ball Mover with radius " << radius_ << std::endl;
+
+}
+
 } // namespace IMP
Index: kernel/src/decorators/XYZDecorator.cpp
===================================================================
--- kernel/src/decorators/XYZDecorator.cpp	(revision 681)
+++ kernel/src/decorators/XYZDecorator.cpp	(working copy)
@@ -6,9 +6,7 @@
  */
#include "IMP/decorators/XYZDecorator.h"
-#include "IMP/random.h"
-#include <boost/random/uniform_real.hpp>
#include <cmath>
@@ -26,35 +24,6 @@
}
-void XYZDecorator::randomize_in_sphere(const Vector3D &center,
-                                       float radius)
-{
-  IMP_check(radius > 0, "Radius in randomize must be postive",
-            ValueException);
-  Vector3D min(center[0]-radius, center[1]-radius, center[2]-radius);
-  Vector3D max(center[0]+radius, center[1]+radius, center[2]+radius);
-  float norm;
-  do {
-    randomize_in_box(min, max);
-    norm=0;
-    for (int i=0; i< 3; ++i) {
-      norm+= square(center[i]-get_coordinate(i));
-    }
-    norm = std::sqrt(norm);
-  } while (norm > radius);
-}
-
-void XYZDecorator::randomize_in_box(const Vector3D &min,
-                                    const Vector3D &max)
-{
-  for (unsigned int i=0; i< 3; ++i) {
-    IMP_check(min[i] < max[i], "Box for randomize must be non-empty",
-              ValueException);
-    ::boost::uniform_real<> rand(min[i], max[i]);
-    set_coordinate(i, rand(random_number_generator));
-  }
-}
-
 IMP_DECORATOR_INITIALIZE(XYZDecorator, DecoratorBase,
                          {
                          key_[0] = FloatKey("x");
@@ -62,16 +31,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/pair_scores/DistancePairScore.cpp
===================================================================
--- kernel/src/pair_scores/DistancePairScore.cpp	(revision 681)
+++ kernel/src/pair_scores/DistancePairScore.cpp	(working copy)
@@ -16,10 +16,6 @@
DistancePairScore::DistancePairScore(UnaryFunction *f): f_(f){}
-struct Identity
-{
-  Float operator()(Float t) const {return t;}
-};
Float DistancePairScore::evaluate(Particle *a, Particle *b,
                                   DerivativeAccumulator *da) const
Index: kernel/src/pair_scores/TransformedDistancePairScore.cpp
===================================================================
--- kernel/src/pair_scores/TransformedDistancePairScore.cpp	(revision 681)
+++ kernel/src/pair_scores/TransformedDistancePairScore.cpp	(working copy)
@@ -52,43 +52,59 @@
 }
+Vector3D TransformedDistancePairScore
+::get_transformed(const Vector3D &v) const {
+  return Vector3D(get_transformed(v,0),
+                  get_transformed(v,1),
+                  get_transformed(v,2));
+}
+
+
+Float TransformedDistancePairScore
+::get_transformed(const Vector3D &v, unsigned int i) const {
+  IMP_assert(i<3, "Bad coordinate");
+  return (v-c_)*r_[i] + tc_[i];
+}
+
+
+Vector3D TransformedDistancePairScore
+::get_back_rotated(const Vector3D &v) const {
+  return Vector3D(v*ri_[0],
+                  v*ri_[1],
+                  v*ri_[2]);
+}
+
+
 /*
   Compute R(x-c)+tc and the reverse
   dt= R(d-c)+tc-> Rt(dt-tc)+c
  */
-struct TransformParticle
+struct TransformedDistancePairScore::TransformParticle
 {
-  const Vector3D *r_, *ri_;
-  const Vector3D &tc_;
-  const Vector3D &c_;
+  const TransformedDistancePairScore *trd_;
   XYZDecorator d_;
-  TransformParticle(const Vector3D *r,
-                    const Vector3D *ri,
-                    const Vector3D &tc,
-                    const Vector3D &c,
-                    Particle *p): r_(r), ri_(ri),
-                                  tc_(tc), c_(c), d_(p){}
+  TransformParticle(const TransformedDistancePairScore *trd,
+                    Particle *p): trd_(trd), d_(p){}
Float get_coordinate(unsigned int i) const {
-    return (d_.get_vector()-c_)*r_[i] + tc_[i];
+    return trd_->get_transformed(d_.get_vector(), i);
   }
void add_to_coordinates_derivative(const Vector3D& f,
                                      DerivativeAccumulator &da) {
-    Vector3D r(f*ri_[0],
-               f*ri_[1],
-               f*ri_[2]);
-    d_.add_to_coordinates_derivative(r, da);
+
+    d_.add_to_coordinates_derivative(trd_->get_back_rotated(f), da);
   }
 };
Float TransformedDistancePairScore::evaluate(Particle *a, Particle *b,
                                              DerivativeAccumulator *da) const
 {
-  TransformParticle tb(r_,ri_, tc_,c_,b);
-  IMP_LOG(VERBOSE, "Transformed particle is "
+  TransformParticle tb(this,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(),
@@ -118,24 +134,55 @@
           << ri_[2] << std::endl);
 }
-void TransformedDistancePairScore::set_translation(float t0, float t1, float t2)
+void TransformedDistancePairScore::set_translation(const Vector3D &t)
 {
-  tc_= Vector3D(t0, t1, t2) + c_;
+  tc_= t + c_;
 }
-void TransformedDistancePairScore::set_center(float t0, float t1, float t2)
+void TransformedDistancePairScore::set_center(const Vector3D &t)
 {
   tc_= tc_-c_;
-  c_= Vector3D(t0, t1, t2);
+  c_= t;
   tc_= tc_+c_;
 }
+void TransformedDistancePairScore::set_rotation(const Vector3D &v, Float theta)
+{
+  IMP_check(v.get_magnitude() <1.1 && v.get_magnitude() > .9,
+            "Vector must be normalized", ValueException);
+  IMP_LOG(TERSE, "Computing rotation from " << v << " and "
+          << theta << std::endl);
+  Float c= std::cos(theta);
+  Float s= std::sin(theta);
+  Float C= 1-c;
+  Vector3D vs = v*s;
+  Vector3D vC = v*C;
+
+  Float xyC = v[0]*vC[1];
+  Float yzC = v[1]*vC[2];
+  Float zxC = v[2]*vC[0];
+  set_rotation(v[0]*vC[0] + c,   xyC - vs[2],   zxC + vs[1],
+               xyC + vs[2],      v[1]*vC[1]+c,  yzC - vs[0],
+               zxC - vs[1],      yzC + vs[0],   v[2]*vC[2] + c);
+}
+
+
 void TransformedDistancePairScore::show(std::ostream &out) const
 {
   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;
+  out << "and center: " << c_ << std::endl;
+  out << "and translation: " << tc_-c_ << std::endl;
 }
} // namespace IMP
Index: kernel/src/Vector3D.cpp
===================================================================
--- kernel/src/Vector3D.cpp	(revision 0)
+++ kernel/src/Vector3D.cpp	(revision 0)
@@ -0,0 +1,63 @@
+/**
+ *  \file Vector3D.cpp   \brief Classes to handle individual a vector.
+ *
+ *  Copyright 2007-8 Sali Lab. All rights reserved.
+ *
+ */
+
+#include "IMP/Vector3D.h"
+#include "IMP/random.h"
+#include "IMP/internal/constants.h"
+#include "IMP/utility.h"
+
+#include <boost/random/uniform_real.hpp>
+
+namespace IMP
+{
+
+Vector3D random_vector_in_box(const Vector3D &min, const Vector3D &max)
+{
+  Vector3D ret;
+  for (unsigned int i=0; i< 3; ++i) {
+    IMP_check(min[i] < max[i], "Box for randomize must be non-empty",
+              ValueException);
+    ::boost::uniform_real<> rand(min[i], max[i]);
+    ret[i]=rand(random_number_generator);
+  }
+  return ret;
+}
+
+
+Vector3D random_vector_in_sphere(const Vector3D &center, Float radius) {
+  IMP_check(radius > 0, "Radius in randomize must be postive",
+            ValueException);
+  Vector3D min(center[0]-radius, center[1]-radius, center[2]-radius);
+  Vector3D max(center[0]+radius, center[1]+radius, center[2]+radius);
+  float norm;
+  Vector3D ret;
+  do {
+    ret=random_vector_in_box(min, max);
+    norm= (center- ret).get_magnitude();
+  } while (norm > radius);
+  return ret;
+}
+
+Vector3D random_vector_on_sphere(const Vector3D &center, Float radius) {
+  IMP_check(radius > 0, "Radius in randomize must be postive",
+            ValueException);
+  ::boost::uniform_real<> rand(-1,1);
+  Vector3D up;
+  up[2]= rand(random_number_generator);
+  ::boost::uniform_real<> trand(0, 2*internal::PI);
+  Float theta= trand(random_number_generator);
+  // radius of circle
+  Float r= std::sqrt(1-square(up[2]));
+  up[0]= std::sin(theta)*r;
+  up[1]= std::cos(theta)*r;
+  IMP_assert(std::abs(up.get_magnitude() -1) < .1,
+             "Error generating unit vector on sphere");
+  IMP_LOG(VERBOSE, "Random vector on sphere is " << up << std::endl);
+  return center+ up*radius;
+}
+
+}  // namespace IMP
Index: kernel/pyext/IMP/test.py
===================================================================
--- kernel/pyext/IMP/test.py	(revision 681)
+++ kernel/pyext/IMP/test.py	(working copy)
@@ -84,7 +84,7 @@
             p= IMP.Particle()
             model.add_particle(p)
             d= IMP.XYZDecorator.create(p)
-            d.randomize_in_box(lbv, ubv)
+            d.set_coordinates(IMP.random_vector_in_box(lbv, ubv))
             ps.append(p)
             d.set_coordinates_are_optimized(True)
         return ps
Index: kernel/pyext/IMP.i
===================================================================
--- kernel/pyext/IMP.i	(revision 681)
+++ kernel/pyext/IMP.i	(working copy)
@@ -96,14 +96,12 @@
       return ret;
     }
   }
-  IMP_CONTAINER_SWIG(Model, ScoreState, score_state)
-  IMP_CONTAINER_SWIG(Model, Restraint, restraint)
-  IMP_CONTAINER_SWIG(RestraintSet, Restraint, restraint)
-  IMP_CONTAINER_SWIG(Optimizer, OptimizerState, optimizer_state)
-  IMP_ADD_OBJECT(MonteCarlo, add_mover)
-  IMP_ADD_OBJECTS(MonteCarlo, add_movers)
-  IMP_ADD_OBJECT(NonbondedListScoreState, add_bonded_list)
-  IMP_ADD_OBJECTS(NonbondedListScoreState, add_bonded_lists)
+  IMP_CONTAINER_SWIG(Model, ScoreState, score_state);
+  IMP_CONTAINER_SWIG(Model, Restraint, restraint);
+  IMP_CONTAINER_SWIG(RestraintSet, Restraint, restraint);
+  IMP_CONTAINER_SWIG(MonteCarlo, Mover, mover);
+  IMP_CONTAINER_SWIG(Optimizer, OptimizerState, optimizer_state);
+  IMP_CONTAINER_SWIG(NonbondedListScoreState, BondedListScoreState, bonded_list);
 }
%feature("ref")   Particle "$this->ref();"
@@ -124,6 +122,7 @@
 %feature("director") IMP::TripletScore;
 %feature("director") IMP::Optimizer;
 %feature("director") IMP::ParticleRefiner;
+%feature("director") IMP::Mover;
%include "IMP/base_types.h"
 %include "IMP/Object.h"
Index: SConstruct
===================================================================
--- SConstruct	(revision 681)
+++ SConstruct	(working copy)
@@ -52,7 +52,7 @@
 #SConscript('domino/SConscript')
# bin script first requires kernel libraries to be built:
-env.Depends(bin, [src, pyext])
+env.Depends(bin, [pyext, src])
# Build the binaries by default:
 env.Default(bin)

    

[IMP-dev] previously broken patches

Daniel Russel