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 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 + +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 ts_; + std::vector 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 + +namespace IMP +{ + +ChainTripletRestraint::ChainTripletRestraint(TripletScore* ts) +{ + ts_ = std::auto_ptr(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(1.0)), + static_cast(-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(1e-12), fact_ij); + fact_kj = std::max(static_cast(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 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 - -#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(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(1.0)), - static_cast(-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(1e-12), fact_ij); - fact_kj = std::max(static_cast(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; }