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)