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 ¢er, 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 ¢er=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 ¢er=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 ¢er,
- 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 ¢er, 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 ¢er, 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)
