No api changes. The internal evaluation functions have been split in two so they are more generally applicable. Also more tests for the transformed distnace pair score

Index: kernel/src/pair_scores/DistancePairScore.cpp =================================================================== --- kernel/src/pair_scores/DistancePairScore.cpp (revision 672) +++ 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 672) +++ kernel/src/pair_scores/TransformedDistancePairScore.cpp (working copy) @@ -52,46 +52,62 @@ }

+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(), + da, f_.get(), boost::lambda::_1); }

@@ -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/include/IMP/pair_scores/TransformedDistancePairScore.h =================================================================== --- kernel/include/IMP/pair_scores/TransformedDistancePairScore.h (revision 672) +++ 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/evaluate_distance_pair_score.h =================================================================== --- kernel/include/IMP/internal/evaluate_distance_pair_score.h (revision 672) +++ kernel/include/IMP/internal/evaluate_distance_pair_score.h (working copy) @@ -9,6 +9,8 @@ #define __IMP_EVALUATE_DISTANCE_PAIR_SCORE_H

#include "../Vector3D.h" +#include <boost/tuple/tuple.hpp> +#include <IMP/UnaryFunction.h>

namespace IMP { @@ -16,41 +18,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 { + // 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; - - score = f->evaluate_with_derivative(shifted_distance, deriv); - - 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/test/pair_scores/test_transform.py =================================================================== --- kernel/test/pair_scores/test_transform.py (revision 672) +++ kernel/test/pair_scores/test_transform.py (working copy) @@ -4,9 +4,10 @@

class DistanceTests(IMP.test.TestCase): """Test the symmetry restraint""" - def test_symmetry(self): + def _test_symmetry(self): """Test the transform pair score basics""" IMP.set_log_level(IMP.VERBOSE) + random.seed() m= IMP.Model() p0= IMP.Particle() m.add_particle(p0) @@ -14,7 +15,7 @@ p1= IMP.Particle() m.add_particle(p1) d1= IMP.XYZDecorator.create(p1) - tps= IMP.TransformedDistancePairScore(IMP.Harmonic(0,1)) + tps= IMP.dsr.TransformedDistancePairScore(IMP.Harmonic(0,1)) tps.set_translation(0,1,0) tps.set_rotation(1, 0, 0, 0, 1, 0, @@ -56,7 +57,7 @@ self.assert_(d1.get_coordinate_derivative(1) > 0) self.assert_(d1.get_coordinate_derivative(0) == 0) self.assert_(d1.get_coordinate_derivative(2) == 0) - def test_symmetry2(self): + def _test_symmetry2(self): """Test the transform pair score optimization""" IMP.set_log_level(IMP.VERBOSE) m= IMP.Model() @@ -70,27 +71,122 @@ d1.set_coordinates(IMP.Vector3D(20,20,40)) 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) + tps= IMP.dsr.TransformedDistancePairScore(IMP.Harmonic(0,1)) + cv= IMP.Vector3D(2,1,0) + tv= IMP.Vector3D(0,1,0) + tps.set_translation(tv[0], tv[1], tv[2]) + tps.set_center(cv[0], cv[1], cv[2]) tps.set_rotation(1, 0, 0, 0, 0,-1, 0, 1, 0) - pr= IMP.PairListRestraint(tps) + pr= IMP.dsr.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): + IMP.set_log_level(IMP.VERBOSE) + m= IMP.Model() + random.seed(100) + 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.dsr.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.dsr.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.dsr.MonteCarlo() + cg.set_model(m) + cg.set_temperature(0) + bm= IMP.dsr.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()