Index: kernel/include/IMP/DecoratorBase.h
===================================================================
--- kernel/include/IMP/DecoratorBase.h (revision 358)
+++ kernel/include/IMP/DecoratorBase.h (working copy)
@@ -9,6 +9,7 @@
#define __IMP_DECORATOR_BASE_H
#include "Object.h"
+#include "Particle.h"
namespace IMP
{
Index: kernel/include/IMP/TripletScore.h
===================================================================
--- kernel/include/IMP/TripletScore.h (revision 0)
+++ kernel/include/IMP/TripletScore.h (revision 0)
@@ -0,0 +1,39 @@
+/**
+ * \file TripletScore.h \brief A Score on a triplet of particles.
+ *
+ * Copyright 2007-8 Sali Lab. All rights reserved.
+ */
+
+#ifndef __IMP_TRIPLET_SCORE_H
+#define __IMP_TRIPLET_SCORE_H
+
+#include "IMP_config.h"
+#include "base_types.h"
+#include "Object.h"
+#include "Particle.h"
+#include "DerivativeAccumulator.h"
+
+namespace IMP
+{
+/**
+ \ingroup restraint
+ \addtogroup tripletscore Score functions on three particles
+ Score functions to by applied to a triplet of particles. These can be
+ used to make more flexible restraints.
+ */
+
+//! Abstract score function for a triplet of particles.
+class IMPDLLEXPORT TripletScore : public Object
+{
+public:
+ TripletScore() {}
+ virtual ~TripletScore() {}
+ //! Compute the score for the triplet and the derivative if needed.
+ virtual Float evaluate(Particle *a, Particle *b, Particle *c,
+ DerivativeAccumulator *da) = 0;
+ virtual void show(std::ostream &out=std::cout) const = 0;
+};
+
+} // namespace IMP
+
+#endif /* __IMP_TRIPLET_SCORE_H */
Index: kernel/include/IMP/triplet_scores/SConscript
===================================================================
--- kernel/include/IMP/triplet_scores/SConscript (revision 0)
+++ kernel/include/IMP/triplet_scores/SConscript (revision 0)
@@ -0,0 +1,7 @@
+Import('env')
+import os.path
+files=[]
+files.append( 'AngleTripletScore.h' )
+includedir = os.path.join(env['includedir'], 'IMP', 'triplet_scores' )
+inst = env.Install(includedir, files)
+env.Alias('install', inst)
Index: kernel/include/IMP/triplet_scores/AngleTripletScore.h
===================================================================
--- kernel/include/IMP/triplet_scores/AngleTripletScore.h (revision 0)
+++ kernel/include/IMP/triplet_scores/AngleTripletScore.h (revision 0)
@@ -0,0 +1,34 @@
+/**
+ * \file AngleTripletScore.h
+ * \brief A Score on the angle between three of particles.
+ *
+ * Copyright 2007-8 Sali Lab. All rights reserved.
+ */
+
+#ifndef __IMP_ANGLE_TRIPLET_SCORE_H
+#define __IMP_ANGLE_TRIPLET_SCORE_H
+
+#include "../TripletScore.h"
+#include "../UnaryFunction.h"
+
+namespace IMP
+{
+
+
+//! Apply a function to the angle between three particles.
+/** \ingroup triplet
+ */
+class IMPDLLEXPORT AngleTripletScore : public TripletScore
+{
+ std::auto_ptr<UnaryFunction> f_;
+public:
+ AngleTripletScore(UnaryFunction *f);
+ virtual ~AngleTripletScore(){}
+ virtual Float evaluate(Particle *a, Particle *b, Particle *c,
+ DerivativeAccumulator *da);
+ virtual void show(std::ostream &out=std::cout) const;
+};
+
+} // namespace IMP
+
+#endif /* __IMP_ANGLE_TRIPLET_SCORE_H */
Index: kernel/include/IMP/restraints/ChainTripletRestraint.h
===================================================================
--- kernel/include/IMP/restraints/ChainTripletRestraint.h (revision 0)
+++ kernel/include/IMP/restraints/ChainTripletRestraint.h (revision 0)
@@ -0,0 +1,49 @@
+/**
+ * \file ChainTripletRestraint.h
+ * \brief Restraint triplets of particles in chains.
+ *
+ * Copyright 2007-8 Sali Lab. All rights reserved.
+ *
+ */
+
+#ifndef __IMP_ANGLES_RESTRAINT_H
+#define __IMP_ANGLES_RESTRAINT_H
+
+#include "../IMP_config.h"
+#include "../Restraint.h"
+#include <vector>
+
+namespace IMP
+{
+
+class TripletScore;
+
+//! Restrain each three consecutive particles in each chain.
+/** \ingroup restraint
+ */
+class IMPDLLEXPORT ChainTripletRestraint : public Restraint
+{
+public:
+ //! Create the triplet restraint.
+ /** \param[in] trip_score Triplet score to apply.
+ */
+ ChainTripletRestraint(TripletScore* trip_score);
+ virtual ~ChainTripletRestraint(){}
+
+ IMP_RESTRAINT("0.5", "Daniel Russel");
+
+ //! Add a chain of particles
+ /** Each three successive particles are restrained*/
+ void add_chain(const Particles &ps);
+
+ //! Clear all the stored chains
+ void clear_chains();
+
+protected:
+ std::auto_ptr<TripletScore> ts_;
+ std::vector<unsigned int> chain_splits_;
+};
+
+} // namespace IMP
+
+#endif /* __IMP_ANGLES_RESTRAINT_H */
Index: kernel/src/restraints/ChainTripletRestraint.cpp
===================================================================
--- kernel/src/restraints/ChainTripletRestraint.cpp (revision 0)
+++ kernel/src/restraints/ChainTripletRestraint.cpp (revision 0)
@@ -0,0 +1,87 @@
+/**
+ * \file ChainTripletRestraint.cpp \brief Triplet restraint.
+ *
+ * Copyright 2007-8 Sali Lab. All rights reserved.
+ *
+ */
+
+#include "IMP/restraints/ChainTripletRestraint.h"
+#include "IMP/decorators/XYZDecorator.h"
+#include "IMP/Particle.h"
+#include "IMP/Model.h"
+#include "IMP/log.h"
+#include "IMP/Vector3D.h"
+#include "IMP/TripletScore.h"
+
+#include <cmath>
+
+namespace IMP
+{
+
+ChainTripletRestraint::ChainTripletRestraint(TripletScore* ts)
+{
+ ts_ = std::auto_ptr<TripletScore>(ts);
+ clear_chains();
+}
+
+void ChainTripletRestraint::add_chain(const Particles &ps)
+{
+ if (ps.size() <3) {
+ IMP_WARN("Adding a chain of length 2 or less to the AnglesRestraint"
+ << " doesn't accomplish anything."<< std::endl);
+ } else {
+ Restraint::add_particles(ps);
+ chain_splits_.back()= Restraint::number_of_particles();
+ chain_splits_.push_back(Restraint::number_of_particles());
+ }
+}
+
+Float ChainTripletRestraint::evaluate(DerivativeAccumulator *accum)
+{
+ int cur_break=0;
+ unsigned int i=2;
+ float score=0;
+ while (i < Restraint::number_of_particles()) {
+ /*IMP_LOG(VERBOSE, "Chain eval on "
+ << Restraint::get_particle(i-2)->get_index()
+ << Restraint::get_particle(i-1)->get_index()
+ << Restraint::get_particle(i)->get_index()
+ << " split is " << chain_splits_[cur_break]
+ << std::endl);*/
+ score += ts_->evaluate(Restraint::get_particle(i-2),
+ Restraint::get_particle(i-1),
+ Restraint::get_particle(i),
+ accum);
+ if (chain_splits_[cur_break] == i) {
+ i+=3;
+ ++cur_break;
+ } else {
+ ++i;
+ }
+ }
+ return score;
+}
+
+void ChainTripletRestraint::clear_chains() {
+ Restraint::clear_particles();
+ chain_splits_.clear();
+ chain_splits_.push_back(0);
+}
+
+
+void ChainTripletRestraint::show(std::ostream& out) const
+{
+ if (get_is_active()) {
+ out << "Chain triplet restraint (active):" << std::endl;
+ } else {
+ out << "Chain triplet restraint (inactive):" << std::endl;
+ }
+
+ out << " version: " << version() << " ";
+ out << " last_modified_by: " << last_modified_by() << std::endl;
+ out << " " << chain_splits_.size()-1 << " chains" << std::endl;
+ ts_->show(out);
+ out << std::endl;
+}
+
+} // namespace IMP
Index: kernel/src/triplet_scores/AngleTripletScore.cpp
===================================================================
--- kernel/src/triplet_scores/AngleTripletScore.cpp (revision 0)
+++ kernel/src/triplet_scores/AngleTripletScore.cpp (revision 0)
@@ -0,0 +1,83 @@
+/**
+ * \file AngleTripletScore.cpp
+ * \brief A Score on the angle between a triplet of particles.
+ *
+ * Copyright 2007-8 Sali Lab. All rights reserved.
+ */
+
+#include "IMP/triplet_scores/AngleTripletScore.h"
+#include "IMP/decorators/XYZDecorator.h"
+#include "IMP/UnaryFunction.h"
+
+namespace IMP
+{
+
+AngleTripletScore::AngleTripletScore(UnaryFunction *f): f_(f){}
+
+Float AngleTripletScore::evaluate(Particle *a, Particle *b, Particle *c,
+ DerivativeAccumulator *da)
+{
+ IMP_CHECK_OBJECT(f_.get());
+ IMP_CHECK_OBJECT(a);
+ IMP_CHECK_OBJECT(b);
+ IMP_CHECK_OBJECT(c);
+ XYZDecorator d0 = XYZDecorator::cast(a);
+ XYZDecorator d1 = XYZDecorator::cast(b);
+ XYZDecorator d2 = XYZDecorator::cast(c);
+
+ Vector3D rij = d1.get_vector_to(d0);
+ Vector3D rkj = d1.get_vector_to(d2);
+
+ Float scalar_product = rij.scalar_product(rkj);
+ Float mag_rij = rij.magnitude();
+ Float mag_rkj = rkj.magnitude();
+ Float mag_product = mag_rij * mag_rkj;
+
+ // avoid division by zero
+ Float cosangle = std::abs(mag_product) > 1e-12 ? scalar_product / mag_product
+ : 0.0;
+
+ // avoid range error for acos
+ cosangle = std::max(std::min(cosangle, static_cast<Float>(1.0)),
+ static_cast<Float>(-1.0));
+
+ Float angle = std::acos(cosangle);
+ Float score;
+
+ if (da) {
+ Float deriv;
+ score = (*f_)(angle, deriv);
+
+ Vector3D unit_rij = rij.get_unit_vector();
+ Vector3D unit_rkj = rkj.get_unit_vector();
+
+ Float sinangle = std::abs(std::sin(angle));
+
+ Float fact_ij = sinangle * mag_rij;
+ Float fact_kj = sinangle * mag_rkj;
+ // avoid division by zero
+ fact_ij = std::max(static_cast<float>(1e-12), fact_ij);
+ fact_kj = std::max(static_cast<float>(1e-12), fact_kj);
+
+ for (int i = 0; i < 3; ++i) {
+ Float derv0 = deriv * (rij.get_component(i) * cosangle
+ - rkj.get_component(i)) / fact_ij;
+ Float derv2 = deriv * (rkj.get_component(i) * cosangle
+ - rij.get_component(i)) / fact_kj;
+ d0.add_to_coordinate_derivative(i, derv0, *da);
+ d1.add_to_coordinate_derivative(i, -derv0 - derv2, *da);
+ d2.add_to_coordinate_derivative(i, derv2, *da);
+ }
+ } else {
+ score = (*f_)(angle);
+ }
+ return score;
+}
+
+void AngleTripletScore::show(std::ostream &out) const
+{
+ out << "AngleTripletScore using ";
+ f_->show(out);
+}
+
+} // namespace IMP
Index: kernel/test/restraints/test_angles.py
===================================================================
--- kernel/test/restraints/test_angles.py (revision 0)
+++ kernel/test/restraints/test_angles.py (revision 0)
@@ -0,0 +1,77 @@
+import unittest
+import IMP
+import IMP.test
+import IMP.utils
+import math
+
+class One(IMP.UnaryFunction):
+ def __init__(self):
+ IMP.UnaryFunction.__init__(self)
+ def __call__(self, *args):
+ return 1
+ def show(self, *args):
+ print "identity"
+
+
+class AngleRestraintTests(IMP.test.TestCase):
+ """Tests for angle restraints"""
+ def setUp(self):
+ # for some reason objects get destroyed even though
+ # python doesn't own them
+ self.stupid_hack=[]
+
+ def create_particles(self, m, n):
+ l0=IMP.Particles()
+ for i in range(0,n):
+ p= IMP.Particle()
+ m.add_particle(p)
+ d= IMP.XYZDecorator.create(p)
+ l0.append(p)
+ self.randomize_particles(l0, 20)
+ return l0
+
+ def create_angle_r(self, s, ps):
+ print s.thisown
+ for i in range(2, len(ps)):
+ print str(i)
+ l= One()
+ self.stupid_hack.append(l)
+ r= IMP.AngleRestraint(ps[i-2],ps[i-1], ps[i], l)
+ s.add_restraint(r)
+ #print r.evaluate(None)
+ #print r.thisown
+ #print l.thisown
+ #print s.evaluate(None)
+ print s.evaluate(None)
+
+ def test_score(self):
+ """Check score of chain triples restraints"""
+ IMP.set_log_level(IMP.VERBOSE)
+ m=IMP.Model()
+ l0= self.create_particles(m, 3)
+ l1= self.create_particles(m, 10)
+ l2= self.create_particles(m, 2)
+ l= One()
+ t= IMP.AngleTripletScore(l)
+ r= IMP.ChainTripletRestraint(t)
+ r.add_chain(l0)
+ r.add_chain(l1)
+ r.add_chain(l2)
+ m.add_restraint(r)
+ s= IMP.RestraintSet("angle restraints")
+ m.add_restraint(s)
+ print "creating angle restraints"
+ self.create_angle_r(s, l0)
+ self.create_angle_r(s, l1)
+ self.create_angle_r(s, l2)
+ ss= s.evaluate(None)
+ rs= r.evaluate(None)
+ diff = ss-rs
+ print ss
+ print rs
+ if (diff <0): diff=-diff
+ self.assert_(diff < .001, "The restraints are not equal")
+
+
+if __name__ == '__main__':
+ unittest.main()
Index: kernel/include/IMP/restraints/AngleRestraint.h
===================================================================
--- kernel/include/IMP/restraints/AngleRestraint.h (revision 358)
+++ kernel/include/IMP/restraints/AngleRestraint.h (working copy)
@@ -9,12 +9,13 @@
#define __IMP_ANGLE_RESTRAINT_H
#include "../IMP_config.h"
-#include "../UnaryFunction.h"
#include "../Restraint.h"
namespace IMP
{
+class AngleTripletScore;
+class UnaryFunction;
//! Angle restraint between three particles
class IMPDLLEXPORT AngleRestraint : public Restraint
@@ -28,13 +29,12 @@
*/
AngleRestraint(Particle* p1, Particle* p2, Particle* p3,
UnaryFunction* score_func);
- virtual ~AngleRestraint();
+ virtual ~AngleRestraint(){}
- IMP_RESTRAINT("0.1", "Ben Webb")
+ IMP_RESTRAINT("0.2", "Daniel Russel")
protected:
- //! scoring function for this restraint
- UnaryFunction* score_func_;
+ std::auto_ptr<AngleTripletScore> sf_;
};
} // namespace IMP
Index: kernel/src/restraints/AngleRestraint.cpp
===================================================================
--- kernel/src/restraints/AngleRestraint.cpp (revision 358)
+++ kernel/src/restraints/AngleRestraint.cpp (working copy)
@@ -5,14 +5,8 @@
*
*/
-#include <cmath>
-
-#include "IMP/Particle.h"
-#include "IMP/Model.h"
-#include "IMP/log.h"
-#include "IMP/Vector3D.h"
#include "IMP/restraints/AngleRestraint.h"
-#include "IMP/decorators/XYZDecorator.h"
+#include "IMP/triplet_scores/AngleTripletScore.h"
namespace IMP
{
@@ -24,17 +18,9 @@
add_particle(p2);
add_particle(p3);
- score_func_ = score_func;
+ sf_= std::auto_ptr<AngleTripletScore>(new AngleTripletScore(score_func));
}
-
-//! Destructor
-AngleRestraint::~AngleRestraint()
-{
- delete score_func_;
-}
-
-
//! Calculate the score for this angle restraint.
/** \param[in] accum If not NULL, use this object to accumulate partial first
derivatives.
@@ -42,58 +28,10 @@
*/
Float AngleRestraint::evaluate(DerivativeAccumulator *accum)
{
- IMP_CHECK_OBJECT(score_func_);
- XYZDecorator d0 = XYZDecorator::cast(get_particle(0));
- XYZDecorator d1 = XYZDecorator::cast(get_particle(1));
- XYZDecorator d2 = XYZDecorator::cast(get_particle(2));
-
- Vector3D rij = d1.get_vector_to(d0);
- Vector3D rkj = d1.get_vector_to(d2);
-
- Float scalar_product = rij.scalar_product(rkj);
- Float mag_rij = rij.magnitude();
- Float mag_rkj = rkj.magnitude();
- Float mag_product = mag_rij * mag_rkj;
-
- // avoid division by zero
- Float cosangle = std::abs(mag_product) > 1e-12 ? scalar_product / mag_product
- : 0.0;
-
- // avoid range error for acos
- cosangle = std::max(std::min(cosangle, static_cast<Float>(1.0)),
- static_cast<Float>(-1.0));
-
- Float angle = std::acos(cosangle);
- Float score;
-
- if (accum) {
- Float deriv;
- score = (*score_func_)(angle, deriv);
-
- Vector3D unit_rij = rij.get_unit_vector();
- Vector3D unit_rkj = rkj.get_unit_vector();
-
- Float sinangle = std::abs(std::sin(angle));
-
- Float fact_ij = sinangle * mag_rij;
- Float fact_kj = sinangle * mag_rkj;
- // avoid division by zero
- fact_ij = std::max(static_cast<float>(1e-12), fact_ij);
- fact_kj = std::max(static_cast<float>(1e-12), fact_kj);
-
- for (int i = 0; i < 3; ++i) {
- Float derv0 = deriv * (rij.get_component(i) * cosangle
- - rkj.get_component(i)) / fact_ij;
- Float derv2 = deriv * (rkj.get_component(i) * cosangle
- - rij.get_component(i)) / fact_kj;
- d0.add_to_coordinate_derivative(i, derv0, *accum);
- d1.add_to_coordinate_derivative(i, -derv0 - derv2, *accum);
- d2.add_to_coordinate_derivative(i, derv2, *accum);
- }
- } else {
- score = (*score_func_)(angle);
- }
- return score;
+ return sf_->evaluate(get_particle(0),
+ get_particle(1),
+ get_particle(2),
+ accum);
}
@@ -114,7 +52,7 @@
out << ", " << get_particle(1)->get_index();
out << " and " << get_particle(2)->get_index();
out << " ";
- score_func_->show(out);
+ sf_->show(out);
out << std::endl;
}