August 2008 - IMP-dev

warning on unused restraints
by Daniel Russel 26 Aug '08

26 Aug '08

Since I got bitten by forgetting to add a restraint to a model I added a warning. Restraint will print out a message if it is destroyed without ever having been added to a Model. This warning can be disabled for a restraint by calling Restraint::set_was_owned(true). Index: kernel/include/IMP/Restraint.h =================================================================== --- kernel/include/IMP/Restraint.h (revision 672) +++ kernel/include/IMP/Restraint.h (working copy) @@ -43,7 +43,7 @@ constructor and should provide methods so that the set of particles can be modified after construction. - A restraint can be added to the model multiple times or to multiple + A restraint can be added to the model multiple times or to multiple restraint sets in the same model. \note When logging is VERBOSE, restraints should print enough information @@ -53,6 +53,13 @@ \note Physical restraints should use the units of kcal/mol for restraint values and kcal/mol/A for derivatives. + + \note Restraints will print a warning message if they are destroyed + without ever having been added to a model as this is an easy mistake + to make. To disable this warning for a particular restraint, call + set_was_owned(true). + + \ingroup kernel */ class IMPDLLEXPORT Restraint : public RefCountedObject { @@ -96,9 +103,16 @@ Model *get_model() const { IMP_assert(model_, "get_model() called before set_model()"); - return model_.get(); + return model_; } + /** A warning is printed if a restraint is destroyed + without ever having belonged to a restraint set or a model. + */ + void set_was_owned(bool tf) { + was_owned_=tf; + } + //! Return a list of sets of particles that are restrained by this restraint /** This function returns a list of sets of particles that are interacting within this restraint. Particles can appear in more @@ -114,12 +128,20 @@ IMP_LIST(protected, Particle, particle, Particle*) private: - Pointer<Model> model_; + /* This pointer should never be ref counted as Model has a + pointer to this object. Not that Model is refcounted yet. + */ + Model *model_; /* True if restraint has not been deactivated. If it is not active, evaluate should not be called */ bool is_active_; + + /* keep track of whether the restraint ever was in a model. + Give warnings on destruction if it was not. + */ + bool was_owned_; }; IMP_OUTPUT_OPERATOR(Restraint); Index: kernel/src/Restraint.cpp =================================================================== --- kernel/src/Restraint.cpp (revision 672) +++ kernel/src/Restraint.cpp (working copy) @@ -17,6 +17,7 @@ Restraint::Restraint() { + model_=NULL; is_active_ = true; // active by default } @@ -24,6 +25,12 @@ //! Destructor Restraint::~Restraint() { + if (!was_owned_) { + // can't use virtual functions in the destructor + std::cerr << "Restraint " << this << " is being destroyed " + << "without ever having been added to a model." + << std::endl; + } } @@ -43,9 +50,10 @@ { IMP_assert(model==NULL || get_number_of_particles()==0 || model == get_particle(0)->get_model() - || (model_ && model_.get() == model), + || (model_ && model_ == model), "Model* different from Particle Model*"); model_=model; + was_owned_=true; } void Restraint::show(std::ostream& out) const

1 0

Re: [IMP-dev] [IMP-commits] r672 - in trunk/kernel: pyext test/optimizers
by Daniel Russel 26 Aug '08

26 Aug '08

Why did you skip the MonteCarlo and NonbondedList parts? Notification of IMP commits wrote: > Author: ben(a)SALILAB.ORG > Date: 2008-08-26 16:59:29 -0700 (Tue, 26 Aug 2008) > New Revision: 672 > > Modified: > trunk/kernel/pyext/IMP.i > trunk/kernel/pyext/IMP_macros.i > trunk/kernel/test/optimizers/test_md_optimizer.py > Log: > Patch from Daniel, with minor changes and testcases by me: > - Add macros to handle the thisown stuff in IMP.i (and add a few missed ones). > - Enable Optimizer::get_optimizer_states() in Python code. > > > Modified: trunk/kernel/pyext/IMP.i > =================================================================== > --- trunk/kernel/pyext/IMP.i 2008-08-26 23:40:30 UTC (rev 671) > +++ trunk/kernel/pyext/IMP.i 2008-08-26 23:59:29 UTC (rev 672) > @@ -27,94 +27,30 @@ > %} > > namespace IMP { > - %pythonprepend Model::add_restraint %{ > - args[1].thisown=0 > - %} > - %pythonprepend Model::add_score_state %{ > - args[1].thisown=0 > - %} > - %pythonprepend Optimizer::add_optimizer_state %{ > - args[1].thisown=0 > - %} > - %pythonprepend RestraintSet::add_restraint %{ > - args[1].thisown=0 > - %} > - %pythonprepend NonbondedListScoreState::add_bonded_list %{ > - args[1].thisown=0 > - %} > - %pythonprepend DistanceRestraint::DistanceRestraint %{ > - args[0].thisown=0 > - %} > - %pythonprepend AngleRestraint::AngleRestraint %{ > - args[0].thisown=0 > - %} > - %pythonprepend DihedralRestraint::DihedralRestraint %{ > - args[0].thisown=0 > - %} > - %pythonprepend TorusRestraint::TorusRestraint %{ > - args[3].thisown=0 > - %} > - %pythonprepend NonbondedRestraint::NonbondedRestraint %{ > - args[0].thisown=0 > - %} > - %pythonprepend BondDecoratorRestraint::BondDecoratorRestraint %{ > - args[0].thisown=0 > - %} > - %pythonprepend SingletonListRestraint::SingletonListRestraint %{ > - args[0].thisown=0 > - %} > - %pythonprepend PairListRestraint::PairListRestraint %{ > - args[0].thisown=0 > - %} > - %pythonprepend TripletChainRestraint::TripletChainRestraint %{ > - args[0].thisown=0 > - %} > - %pythonprepend PairChainRestraint::PairChainRestraint %{ > - args[0].thisown=0 > - %} > - %pythonprepend ConnectivityRestraint::ConnectivityRestraint %{ > - args[0].thisown=0 > - %} > - %pythonprepend DistancePairScore::DistancePairScore %{ > - args[0].thisown=0 > - %} > - %pythonprepend TransformedDistancePairScore::TransformedDistancePairScore %{ > - args[0].thisown=0 > - %} > - %pythonprepend BondCoverPairScore::BondCoverPairScore %{ > - args[0].thisown=0 > - %} > - %pythonprepend SphereDistancePairScore::SphereDistancePairScore %{ > - args[0].thisown=0 > - %} > - %pythonprepend RefineOncePairScore::RefineOncePairScore %{ > - args[0].thisown=0 > - args[1].thisown=0 > - %} > - %pythonprepend DistanceToSingletonScore::DistanceToSingletonScore %{ > - args[0].thisown=0 > - %} > - %pythonprepend AttributeSingletonScore::AttributeSingletonScore %{ > - args[0].thisown=0 > - %} > - %pythonprepend TunnelSingletonScore::TunnelSingletonScore %{ > - args[0].thisown=0 > - %} > - %pythonprepend AngleTripletScore::AngleTripletScore %{ > - args[0].thisown=0 > - %} > - %pythonprepend MonteCarlo::add_mover %{ > - args[1].thisown=0 > - %} > - %pythonprepend MonteCarlo::set_local_optimizer %{ > - args[1].thisown=0 > - %} > - %pythonprepend VRMLLogOptimizerState::add_particle_refiner %{ > - args[1].thisown=0 > - %} > - %pythonprepend TypedPairScore::set_pair_score %{ > - args[1].thisown=0 > - %} > + // need to special case particle so can't add this to macro > + IMP_OWN_FIRST_CONSTRUCTOR(DistanceRestraint) > + IMP_OWN_FIRST_CONSTRUCTOR(AngleRestraint) > + IMP_OWN_FIRST_CONSTRUCTOR(DihedralRestraint) > + IMP_OWN_FIRST_CONSTRUCTOR(TorusRestraint) > + IMP_OWN_FIRST_CONSTRUCTOR(NonbondedRestraint) > + IMP_OWN_FIRST_CONSTRUCTOR(BondDecoratorRestraint) > + IMP_OWN_FIRST_CONSTRUCTOR(SingletonListRestraint) > + IMP_OWN_FIRST_CONSTRUCTOR(PairListRestraint) > + IMP_OWN_FIRST_CONSTRUCTOR(TripletChainRestraint) > + IMP_OWN_FIRST_CONSTRUCTOR(PairChainRestraint) > + IMP_OWN_FIRST_CONSTRUCTOR(ConnectivityRestraint) > + IMP_OWN_FIRST_CONSTRUCTOR(DistancePairScore) > + IMP_OWN_FIRST_CONSTRUCTOR(TransformedDistancePairScore) > + IMP_OWN_FIRST_CONSTRUCTOR(BondCoverPairScore) > + IMP_OWN_FIRST_CONSTRUCTOR(SphereDistancePairScore) > + IMP_OWN_FIRST_SECOND_CONSTRUCTOR(RefineOncePairScore) > + IMP_OWN_FIRST_CONSTRUCTOR(DistanceToSingletonScore) > + IMP_OWN_FIRST_CONSTRUCTOR(AttributeSingletonScore) > + IMP_OWN_FIRST_CONSTRUCTOR(TunnelSingletonScore) > + IMP_OWN_FIRST_CONSTRUCTOR(AngleTripletScore) > + IMP_SET_OBJECT(MonteCarlo, set_local_optimizer) > + IMP_SET_OBJECT(TypedPairScore, set_pair_score) > + > %pythonprepend Particle::get_value %{ > check_particle(args[0], args[1]) > %} > @@ -144,10 +80,17 @@ > > %} > > - IMP_CONTAINER_SWIG(Model, Particle, particle); > + // special case since particles are ref counted > + %extend Model { > + Particles get_particles() const { > + IMP::Particles ret(self->particles_begin(), self->particles_end()); > + 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); > } > > %feature("ref") Particle "$this->ref();" > > Modified: trunk/kernel/pyext/IMP_macros.i > =================================================================== > --- trunk/kernel/pyext/IMP_macros.i 2008-08-26 23:40:30 UTC (rev 671) > +++ trunk/kernel/pyext/IMP_macros.i 2008-08-26 23:59:29 UTC (rev 672) > @@ -1,6 +1,38 @@ > #define IMPDLLEXPORT > #define IMPDLLLOCAL > > +%define IMP_OWN_FIRST_CONSTRUCTOR(Ucname) > +%pythonprepend Ucname::Ucname %{ > + args[0].thisown=0 > +%} > +%enddef > + > +%define IMP_OWN_FIRST_SECOND_CONSTRUCTOR(Ucname) > +%pythonprepend Ucname::Ucname %{ > + args[0].thisown=0 > + args[1].thisown=0 > +%} > +%enddef > + > +%define IMP_SET_OBJECT(Ucname, methodname) > +%pythonprepend Ucname::methodname %{ > + args[1].thisown=0 > +%} > +%enddef > + > +%define IMP_ADD_OBJECT(Ucname, methodname) > +%pythonprepend Ucname::methodname %{ > + args[1].thisown=0 > +%} > +%enddef > + > +%define IMP_ADD_OBJECTS(Ucname, methodname) > +%pythonprepend Ucname::methodname %{ > + for p in args[1]: > + p.thisown=0 > +%} > +%enddef > + > /* Add additional IMP_CONTAINER methods for scripting languages */ > %define IMP_CONTAINER_SWIG(type, Ucname, lcname) > %extend type { > @@ -9,6 +41,8 @@ > return ret; > } > } > +IMP_ADD_OBJECT(type, add_##lcname) > +IMP_ADD_OBJECTS(type, add_##lcname##s) > %enddef > > %include "IMP/macros.h" > > Modified: trunk/kernel/test/optimizers/test_md_optimizer.py > =================================================================== > --- trunk/kernel/test/optimizers/test_md_optimizer.py 2008-08-26 23:40:30 UTC (rev 671) > +++ trunk/kernel/test/optimizers/test_md_optimizer.py 2008-08-26 23:59:29 UTC (rev 672) > @@ -159,6 +159,18 @@ > ekinetic2 = self.md.get_kinetic_energy() > self.assertNotEqual(ekinetic, ekinetic2) > > + def test_get_optimizer_states(self): > + """Test get_optimizer_states() method""" > + wrtraj = WriteTrajState([]) > + scaler = IMP.VelocityScalingOptimizerState( > + IMP.Particles(self.particles), 298.0, 10) > + self.md.add_optimizer_state(wrtraj) > + self.md.add_optimizer_state(scaler) > + m = self.md.get_optimizer_states() > + self.assertEqual(len(m), 2) > + for a in m: > + self.assert_(isinstance(a, IMP.OptimizerState)) > + > def test_rescaling(self): > """Test thermostatting by velocity rescaling""" > for i in range(100): > > _______________________________________________ > IMP-commits mailing list > IMP-commits(a)salilab.org > https://salilab.org/mailman/listinfo/imp-commits >

2 2

attribute table bug fix
by Daniel Russel 26 Aug '08

26 Aug '08

Here is a simplification of the attribute tables which also fixes an assertion failure when there are few attributes around. It adds a new file IMP_keys.i which defines some base classes for the Keys. This file should probably be somehow included in the other module's .is. I stuck it in IMPEM.i, but Ben needs to make sure such things are correct. The patch includes the previous patches since they can't easily be separated at this point. Index: em/pyext/IMPEM.i =================================================================== --- em/pyext/IMPEM.i (revision 669) +++ em/pyext/IMPEM.i (working copy) @@ -7,6 +7,7 @@ %} %include "kernel/pyext/IMP_macros.i" +%include "kernel/pyext/IMP_keys.i" /* Ignore shared object import/export stuff */ #define EMDLLEXPORT Index: kernel/pyext/IMP.i =================================================================== --- kernel/pyext/IMP.i (revision 669) +++ kernel/pyext/IMP.i (working copy) @@ -11,11 +11,21 @@ %include "IMP_macros.i" %include "IMP_exceptions.i" +%include "IMP_keys.i" /* Return derivatives from unary functions */ %include "typemaps.i" -%apply double &OUTPUT { IMP::Float& deriv }; +namespace IMP { + %typemap(out) std::pair<Float,Float> { + PyObject *tup= PyTuple_New(2); + PyTuple_SetItem(tup, 0, PyFloat_FromDouble($1.first)); + PyTuple_SetItem(tup, 1, PyFloat_FromDouble($1.second)); + $result= tup; + } +} + + %pythoncode %{ def check_particle(p, a): if (not p.get_is_active()): @@ -27,94 +37,30 @@ %} namespace IMP { - %pythonprepend Model::add_restraint %{ - args[1].thisown=0 - %} - %pythonprepend Model::add_score_state %{ - args[1].thisown=0 - %} - %pythonprepend Optimizer::add_optimizer_state %{ - args[1].thisown=0 - %} - %pythonprepend RestraintSet::add_restraint %{ - args[1].thisown=0 - %} - %pythonprepend NonbondedListScoreState::add_bonded_list %{ - args[1].thisown=0 - %} - %pythonprepend DistanceRestraint::DistanceRestraint %{ - args[0].thisown=0 - %} - %pythonprepend AngleRestraint::AngleRestraint %{ - args[0].thisown=0 - %} - %pythonprepend DihedralRestraint::DihedralRestraint %{ - args[0].thisown=0 - %} - %pythonprepend TorusRestraint::TorusRestraint %{ - args[3].thisown=0 - %} - %pythonprepend NonbondedRestraint::NonbondedRestraint %{ - args[0].thisown=0 - %} - %pythonprepend BondDecoratorRestraint::BondDecoratorRestraint %{ - args[0].thisown=0 - %} - %pythonprepend SingletonListRestraint::SingletonListRestraint %{ - args[0].thisown=0 - %} - %pythonprepend PairListRestraint::PairListRestraint %{ - args[0].thisown=0 - %} - %pythonprepend TripletChainRestraint::TripletChainRestraint %{ - args[0].thisown=0 - %} - %pythonprepend PairChainRestraint::PairChainRestraint %{ - args[0].thisown=0 - %} - %pythonprepend ConnectivityRestraint::ConnectivityRestraint %{ - args[0].thisown=0 - %} - %pythonprepend DistancePairScore::DistancePairScore %{ - args[0].thisown=0 - %} - %pythonprepend TransformedDistancePairScore::TransformedDistancePairScore %{ - args[0].thisown=0 - %} - %pythonprepend BondCoverPairScore::BondCoverPairScore %{ - args[0].thisown=0 - %} - %pythonprepend SphereDistancePairScore::SphereDistancePairScore %{ - args[0].thisown=0 - %} - %pythonprepend RefineOncePairScore::RefineOncePairScore %{ - args[0].thisown=0 - args[1].thisown=0 - %} - %pythonprepend DistanceToSingletonScore::DistanceToSingletonScore %{ - args[0].thisown=0 - %} - %pythonprepend AttributeSingletonScore::AttributeSingletonScore %{ - args[0].thisown=0 - %} - %pythonprepend TunnelSingletonScore::TunnelSingletonScore %{ - args[0].thisown=0 - %} - %pythonprepend AngleTripletScore::AngleTripletScore %{ - args[0].thisown=0 - %} - %pythonprepend MonteCarlo::add_mover %{ - args[1].thisown=0 - %} - %pythonprepend MonteCarlo::set_local_optimizer %{ - args[1].thisown=0 - %} - %pythonprepend VRMLLogOptimizerState::add_particle_refiner %{ - args[1].thisown=0 - %} - %pythonprepend TypedPairScore::set_pair_score %{ - args[1].thisown=0 - %} + // need to special case particle so can't add this to macro + IMP_OWN_FIRST_CONSTRUCTOR(DistanceRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(AngleRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(DihedralRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(TorusRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(NonbondedRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(BondDecoratorRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(SingletonListRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(PairListRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(TripletChainRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(PairChainRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(ConnectivityRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(DistancePairScore) + IMP_OWN_FIRST_CONSTRUCTOR(TransformedDistancePairScore) + IMP_OWN_FIRST_CONSTRUCTOR(BondCoverPairScore) + IMP_OWN_FIRST_CONSTRUCTOR(SphereDistancePairScore) + IMP_OWN_FIRST_SECOND_CONSTRUCTOR(RefineOncePairScore) + IMP_OWN_FIRST_CONSTRUCTOR(DistanceToSingletonScore) + IMP_OWN_FIRST_CONSTRUCTOR(AttributeSingletonScore) + IMP_OWN_FIRST_CONSTRUCTOR(TunnelSingletonScore) + IMP_OWN_FIRST_CONSTRUCTOR(AngleTripletScore) + IMP_SET_OBJECT(MonteCarlo, set_local_optimizer) + IMP_SET_OBJECT(TypedPairScore, set_pair_score) + %pythonprepend Particle::get_value %{ check_particle(args[0], args[1]) %} @@ -144,10 +90,19 @@ %} - IMP_CONTAINER_SWIG(Model, Particle, particle); + // special case since particles are ref counted + %extend Model { + Particles get_particles() const { + IMP::Particles ret(self->particles_begin(), self->particles_end()); + return ret; + } + } 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();" @@ -168,8 +123,10 @@ %feature("director") IMP::TripletScore; %feature("director") IMP::Optimizer; %feature("director") IMP::ParticleRefiner; +%feature("director") IMP::Mover; -%include "IMP/Key.h" + +%include "IMP/base_types.h" %include "IMP/Object.h" %include "IMP/RefCountedObject.h" %include "IMP/Index.h" @@ -259,12 +216,6 @@ %template(OptimizerStateIndex) Index<OptimizerStateTag>; %template(MoverIndex) Index<Mover>; %template(BondedListIndex) Index<BondedListScoreState>; - %template(FloatKey) KeyBase<Float>; - %template(IntKey) KeyBase<Int>; - %template(StringKey) KeyBase<String>; - %template(ParticleKey) KeyBase<Particle*>; - %template(AtomType) KeyBase<AtomTypeTag>; - %template(ResidueType) KeyBase<ResidueTypeTag>; %template(show_named_hierarchy) show<NameDecorator>; %template(show_molecular_hierarchy) show<MolecularHierarchyDecorator>; %template(Particles) ::std::vector<Particle*>; Index: kernel/pyext/IMP_keys.i =================================================================== --- kernel/pyext/IMP_keys.i (revision 0) +++ kernel/pyext/IMP_keys.i (revision 0) @@ -0,0 +1,11 @@ +%include "IMP/Key.h" + +namespace IMP +{ + %template(FloatKeyBase) ::IMP::KeyBase<0>; + %template(IntKeyBase) ::IMP::KeyBase<1>; + %template(StringKeyBase) ::IMP::KeyBase<2>; + %template(ParticleKeyBase) ::IMP::KeyBase<3>; + %template(AtomTypeBase) ::IMP::KeyBase<IMP_ATOM_TYPE_INDEX>; + %template(ResidueTypeBase) ::IMP::KeyBase<IMP_RESIDUE_TYPE_INDEX>; +} Index: kernel/src/Key.cpp =================================================================== --- kernel/src/Key.cpp (revision 669) +++ kernel/src/Key.cpp (working copy) @@ -31,6 +31,9 @@ } } + +std::map<unsigned int, KeyData> key_data; + } // namespace internal } // namespace IMP Index: kernel/src/decorators/ResidueDecorator.cpp =================================================================== --- kernel/src/decorators/ResidueDecorator.cpp (revision 669) +++ kernel/src/decorators/ResidueDecorator.cpp (working copy) @@ -16,8 +16,6 @@ namespace IMP { -IMP_DEFINE_KEY_TYPE(ResidueType, ResidueTypeTag); - IntKey ResidueDecorator::type_key_; IntKey ResidueDecorator::index_key_; Index: kernel/src/decorators/AtomDecorator.cpp =================================================================== --- kernel/src/decorators/AtomDecorator.cpp (revision 669) +++ kernel/src/decorators/AtomDecorator.cpp (working copy) @@ -15,8 +15,6 @@ namespace IMP { -IMP_DEFINE_KEY_TYPE(AtomType, AtomTypeTag); - #define TYPE_INIT(STR) AT_##STR= AtomType(#STR); #define TYPE_INIT2(NAME, STR) AT_##NAME = AtomType(#STR); Index: kernel/src/Particle.cpp =================================================================== --- kernel/src/Particle.cpp (revision 669) +++ kernel/src/Particle.cpp (working copy) @@ -43,17 +43,15 @@ void Particle::zero_derivatives() { - for (FloatKeyIterator it= float_keys_begin(); it != float_keys_end(); ++it) { - derivatives_.set_value(*it, 0); - } + derivatives_.set_values(0); } void Particle::show(std::ostream& out) const { - const char* inset = " "; + const std::string inset(" "); out << std::endl; - out << "--" << get_index() << "--" << std::endl; + out << "Particle: " << get_index() << std::endl; if (is_active_) { out << inset << inset << "active"; } else { @@ -62,18 +60,24 @@ out << std::endl; if (get_model() != NULL) { - out << inset << inset << "float attributes:" << std::endl; - floats_.show(out, " "); + out << inset << "float attributes:" << std::endl; + floats_.show(out, inset+inset); - out << inset << inset << "int attributes:" << std::endl; - ints_.show(out, " "); + out << inset << "float derivatives:" << std::endl; + derivatives_.show(out, inset+inset); - out << inset << inset << "string attributes:" << std::endl; - strings_.show(out, " "); + out << inset << "optimizeds:" << std::endl; + optimizeds_.show(out, inset+inset); - out << inset << inset << "particle attributes:" << std::endl; - particles_.show(out, " "); + out << inset << "int attributes:" << std::endl; + ints_.show(out, inset+inset); + out << inset << "string attributes:" << std::endl; + strings_.show(out, inset+inset); + + out << inset << "particle attributes:" << std::endl; + particles_.show(out, inset+inset); + } } Index: kernel/include/IMP/utility.h =================================================================== --- kernel/include/IMP/utility.h (revision 669) +++ 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,15 @@ return t*t; } +//! put in check for G++ +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/include/IMP/Particle.h =================================================================== --- kernel/include/IMP/Particle.h (revision 669) +++ kernel/include/IMP/Particle.h (working copy) @@ -60,7 +60,7 @@ struct ParticleAttributeTableTraits { typedef Pointer<Particle> Value; - typedef KeyBase<Particle*> Key; + typedef ParticleKey Key; static Value get_invalid() { return Value(); } @@ -69,13 +69,13 @@ } }; - typedef internal::AttributeTable<internal::FloatAttributeTableTraits> + typedef internal::AttributeTable<internal::FloatAttributeTableTraits> FloatTable; - typedef internal::AttributeTable<internal::FloatAttributeTableTraits> + typedef internal::AttributeTable<internal::FloatAttributeTableTraits> DerivativeTable; - typedef internal::AttributeTable<internal::BoolAttributeTableTraits> + typedef internal::AttributeTable<internal::BoolAttributeTableTraits> OptimizedTable; - typedef internal::AttributeTable<internal::IntAttributeTableTraits> + typedef internal::AttributeTable<internal::IntAttributeTableTraits> IntTable; typedef internal::AttributeTable<internal::StringAttributeTableTraits> StringTable; Index: kernel/src/base_types.cpp =================================================================== --- kernel/src/base_types.cpp (revision 669) +++ kernel/src/base_types.cpp (working copy) @@ -10,9 +10,4 @@ namespace IMP { -IMP_DEFINE_KEY_TYPE(FloatKey, Float); -IMP_DEFINE_KEY_TYPE(IntKey, Int); -IMP_DEFINE_KEY_TYPE(StringKey, String); -IMP_DEFINE_KEY_TYPE(ParticleKey, Particle*) - } // namespace IMP Index: kernel/include/IMP/Key.h =================================================================== --- kernel/include/IMP/Key.h (revision 669) +++ kernel/include/IMP/Key.h (working copy) @@ -21,13 +21,11 @@ The keys in IMP maintain a cached mapping between strings and indexes. This mapping is global--that is all IMP Models and Particles in the - same program use the same mapping. + same program use the same mapping for each type of key. The type of + the key is determined by an integer which should be unique for + each type. If the integer is not unique, everything works, just + more memory is wasted and types are interconvertible - The mapping uses a static table which is defined with the - IMP_DEFINE_KEY_TYPE macro. As a result, the macro must be used in exactly - one .o. If it appears more than once or that .o is linked in multiple times, - then bad things can happen. - Since the order of initialization of static data is undefined between .os, it is important that no other static data not in the same .o uses the key data. Specifically, typedes defined by @@ -37,54 +35,43 @@ indices and would make the set of indices less dense. */ -// Swig chokes on the specialization -#ifndef SWIG -//! \internal -#define IMP_KEY_DATA_HOLDER(Name, Tag) \ - namespace internal { \ - template <> \ - struct IMPDLLEXPORT KeyDataHolder<Tag> { \ - static KeyData data; \ - }; \ - } -#else -#define IMP_KEY_DATA_HOLDER(Name, Tag) -#endif - /** - Define a new key type. There must be an accompanying IMP_DEFINE_KEY_TYPE - located in some .o file. This macro should be used in the IMP namespace. + Define a new key type. It defines two public types Name, which is an instantiation of KeyBase and Names which is a vector of Name. \param[in] Name The name for the new type. - \param[in] Tag A class which is unique for this type. For attributes - use the type of the attributes. For other Keys, declare an empty - class with a unique name and use it. + \param[in] Tag A (hopefully) unique integer to define this key type. + + \note We define a new class rather than use a typedef since SWIG has a + bug dealing with names that start with ::. A fix has been commited to SVN + for swig. + + \note The name in the typedef would have to start with ::IMP so it + could be used out of the IMP namespace. */ #define IMP_DECLARE_KEY_TYPE(Name, Tag) \ - IMP_KEY_DATA_HOLDER(Name, Tag); \ - typedef KeyBase<Tag> Name; \ - typedef std::vector<Name> Name##s + struct Name: public ::IMP::KeyBase<Tag> { \ + typedef ::IMP::KeyBase<Tag> P; \ + typedef Name This; \ + Name(){}; \ + Name(unsigned int i): P(i){} \ + Name(const char *nm): P(nm){} \ +}; \ +typedef std::vector<Name> Name##s -/** This must occur in exactly one .o in the internal namespace. Should - be used in the IMP namespace.*/ -#define IMP_DEFINE_KEY_TYPE(Name, Tag) \ - namespace internal { \ - KeyData KeyDataHolder<Tag>::data; \ - } - - namespace IMP { namespace internal { -/** The data concerning a particular type of key. + + +/** The data concerning keys. \internal */ struct IMPDLLEXPORT KeyData @@ -109,16 +96,10 @@ double heuristic_; Map map_; RMap rmap_; - }; -/** A dummy class. Actual keys types create specializations +IMPDLLEXPORT extern std::map<unsigned int, KeyData> key_data; - \internal -*/ -template <class T> -struct KeyDataHolder {}; - } // namespace internal @@ -132,22 +113,45 @@ attribute_table_index. Yes, this is an evil hack, but I couldn't get linking to work with static members of the template class. */ -template <class T> +template <int ID> class KeyBase { int str_; - typedef T Type; + static const internal::KeyData::Map& get_map() + { + return IMP::internal::key_data[ID].get_map(); + } + static const internal::KeyData::RMap& get_rmap() { + return IMP::internal::key_data[ID].get_rmap(); + } - bool is_default() const; - static const internal::KeyData::Map& get_map(); - static const internal::KeyData::RMap& get_rmap(); - + static unsigned int find_index(std::string sc) { + if (get_map().find(sc) == get_map().end()) { + return IMP::internal::key_data[ID].add_key(sc); + } else { + return get_map().find(sc)->second; + } + } +protected: + bool is_default() const; public: - static const std::string &get_string(int i); + static const std::string get_string(int i) + { + if (static_cast<unsigned int>(i) + < get_rmap().size()) { + return get_rmap()[i]; + } else { + IMP_WARN("Corrupted Key Table asking for key " << i + << " with a table of size " << get_rmap().size()); + std::ostringstream oss; + oss<< "Invalid index " << i; + return oss.str(); + } + } - typedef KeyBase<T> This; + typedef KeyBase<ID> This; //! make a default (uninitalized) key KeyBase(): str_(-1) {} @@ -155,16 +159,12 @@ //! Generate a key from the given string KeyBase(const char *c) { - std::string sc(c); - if (get_map().find(sc) == get_map().end()) { - str_= internal::KeyDataHolder<T>::data.add_key(sc); - } else { - str_= get_map().find(sc)->second; - } + str_= find_index(c); } explicit KeyBase(unsigned int i): str_(i) { - //IMP_assert(get_rmap().size() > i, "There is no such attribute " << i); + IMP_assert(str_ >= 0, "Invalid initializer " << i); + // cannot check here as we need a past end iterator } //! Turn a key into a pretty string @@ -223,46 +223,28 @@ #endif }; -IMP_OUTPUT_OPERATOR_1(KeyBase) - -template <class T> -inline const internal::KeyData::Map& KeyBase<T>::get_map() -{ - return internal::KeyDataHolder<T>::data.get_map(); +template <int ID> +std::ostream &operator<<(std::ostream &out, KeyBase<ID> k) { + k.show(out); + return out; } -template <class T> -inline const internal::KeyData::RMap& KeyBase<T>::get_rmap() +template <int ID> +inline bool KeyBase<ID>::is_default() const { - return internal::KeyDataHolder<T>::data.get_rmap(); -} - -template <class T> -inline bool KeyBase<T>::is_default() const -{ return str_==-1; } -template <class T> -inline const std::string &KeyBase<T>::get_string(int i) -{ - IMP_assert(static_cast<unsigned int>(i) - < get_rmap().size(), - "Corrupted " << " Key " << i - << " vs " << get_rmap().size()); - return get_rmap()[i]; -} - -template <class T> -inline void KeyBase<T>::show_all(std::ostream &out) +template <int ID> +inline void KeyBase<ID>::show_all(std::ostream &out) { - internal::KeyDataHolder<T>::data.show(out); + internal::key_data[ID].show(out); } -template <class T> -std::vector<std::string> KeyBase<T>::get_all_strings() +template <int ID> +std::vector<std::string> KeyBase<ID>::get_all_strings() { std::vector<std::string> str; for (internal::KeyData::Map::const_iterator it= get_map().begin(); Index: kernel/include/IMP/internal/AttributeTable.h =================================================================== --- kernel/include/IMP/internal/AttributeTable.h (revision 669) +++ kernel/include/IMP/internal/AttributeTable.h (working copy) @@ -27,10 +27,14 @@ namespace internal { -struct FloatAttributeTableTraits +template <class T, class K > +struct DefaultTraits { + typedef T Value; + typedef K Key; +}; + + struct FloatAttributeTableTraits: public DefaultTraits<Float, FloatKey> { - typedef Float Value; - typedef KeyBase<Float> Key; static Float get_invalid() { if (std::numeric_limits<Float>::has_quiet_NaN) { return std::numeric_limits<Float>::quiet_NaN(); @@ -42,17 +46,15 @@ } 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(); + return f != get_invalid(); } } }; -struct IntAttributeTableTraits + struct IntAttributeTableTraits: public DefaultTraits<Int, IntKey> { - typedef Int Value; - typedef KeyBase<Int> Key; static Int get_invalid() { return std::numeric_limits<Int>::max(); } @@ -61,10 +63,8 @@ } }; -struct BoolAttributeTableTraits +struct BoolAttributeTableTraits: public DefaultTraits<bool, FloatKey > { - typedef bool Value; - typedef KeyBase<Float> Key; static bool get_invalid() { return false; } @@ -73,10 +73,8 @@ } }; -struct StringAttributeTableTraits + struct StringAttributeTableTraits: public DefaultTraits<String, StringKey> { - typedef String Value; - typedef KeyBase<String> Key; static Value get_invalid() { return "This is an invalid string in IMP"; } @@ -87,9 +85,9 @@ // The traits for the particle class are declared in the Particle.h -/** \internal +/** \internal If memory becomes a problem then either use a char table to look up - actual entries or use perfect hashing. + actual entries or use perfect hashing. http://burtleburtle.net/bob/hash/perfect.html Actuall Cuckoo hashing is probably a better bet as that gets @@ -105,42 +103,11 @@ { typedef AttributeTable<Traits> This; - typedef boost::scoped_array<typename Traits::Value> Map; + typedef std::vector<typename Traits::Value> Map; Map map_; - unsigned int size_; - typename Traits::Value* new_array(unsigned int asz) const { - typename Traits::Value* ret= new typename Traits::Value[asz]; - for (unsigned int i=0; i< asz; ++i) { - ret[i]= Traits::get_invalid(); - } - return ret; - } - - void copy_from(unsigned int len, const Map &o) { - IMP_assert(map_, "Internal error in attribute table"); - IMP_assert(size_ >= len, "Table too small"); - //std::cout << "Copy from " << o << " to " << map_ << std::endl; - for (unsigned int i=0; i< len; ++i) { - map_[i]= o[i]; - } - } - - void realloc(unsigned int olen, const Map &o, unsigned int nlen) { - //std::cout << "Realloc from " << size_ << " " << map_ << " "; - if (nlen==0) { - size_=0; - map_.reset(); - } else { - size_=std::max(nlen, 6U); - map_.reset(new_array(size_)); - copy_from(olen, o); - } - //std::cout << " to " << size_ << " " << map_ << std::endl; - } - void check_contains(typename Traits::Key k) const { - IMP_assert(size_ > k.get_index(), + IMP_assert(map_.size() > k.get_index(), "Attribute \"" << k.get_string() << "\" not found in table."); IMP_check(Traits::get_is_valid(map_[k.get_index()]), @@ -149,26 +116,18 @@ IndexException); } + void clear() { + map_.clear(); + } + public: typedef typename Traits::Value Value; typedef typename Traits::Key Key; - AttributeTable(): size_(0){} - AttributeTable(const This &o): size_(0) { - //std::cout << "Copy constructor called from " << o.map_ << std::endl; - realloc(o.size_, o.map_, o.size_); - } + AttributeTable(){} ~AttributeTable() { - //std::cout << "Deleting " << map_ << std::endl; + //std::cout << "Deleting " << map_ << std::endl; } - This &operator=(const This &o) { - //std::cout << "Operator= called from " << o.map_ << std::endl; - if (&o == this) { - //std::cout << "Self assignment" << std::endl; - return *this; - } - realloc(o.size_, o.map_, o.size_); - return *this; - } + const Value get_value(Key k) const { check_contains(k); return map_[k.get_index()]; @@ -183,6 +142,11 @@ map_[k.get_index()] = v; } + void set_values(Value v) { + for (unsigned int i=0; i< map_.size(); ++i) { + map_[i]=v; + } + } void insert(Key k, Value v) { IMP_assert(!contains(k), @@ -195,24 +159,27 @@ void remove(Key k) { check_contains(k); - map_[k.get_index()]= Traits::get_invalid(); + remove_always(k); } void remove_always(Key k) { - if (k.get_index() < size_) { + IMP_assert(k != Key(), "Can't remove invalid key"); + if (k.get_index() < map_.size()) { map_[k.get_index()]= Traits::get_invalid(); + while (!map_.empty() + && map_.back()== Traits::get_invalid()) map_.pop_back(); } } bool contains(Key k) const { IMP_assert(k != Key(), "Can't search for default key"); - return k.get_index() < size_ + return k.get_index() < map_.size() && Traits::get_is_valid(map_[k.get_index()]); } - void show(std::ostream &out, const char *prefix="") const; + void show(std::ostream &out, const std::string prefix="") const; std::vector<Key> get_keys() const; @@ -228,19 +195,26 @@ } }; - typedef boost::counting_iterator<Key, boost::random_access_traversal_tag, - std::size_t> KeyIterator; + typedef boost::counting_iterator<Key, boost::forward_traversal_tag, + unsigned int> KeyIterator; typedef boost::filter_iterator<IsAttribute, KeyIterator> AttributeKeyIterator; AttributeKeyIterator attribute_keys_begin() const { + KeyIterator b(0U); + KeyIterator e(map_.size()); + IMP_assert(std::distance(b,e) + == map_.size(), "Something is broken with the iterators"); + IMP_assert(std::distance(AttributeKeyIterator(IsAttribute(this), b,e), + AttributeKeyIterator(IsAttribute(this), e,e)) + <= map_.size(), "Broken in filter"); return AttributeKeyIterator(IsAttribute(this), KeyIterator(Key(0U)), - KeyIterator(Key(size_))); + KeyIterator(Key(map_.size()))); } AttributeKeyIterator attribute_keys_end() const { return AttributeKeyIterator(IsAttribute(this), - KeyIterator(Key(size_)), - KeyIterator(Key(size_))); + KeyIterator(Key(map_.size())), + KeyIterator(Key(map_.size()))); } }; @@ -253,19 +227,20 @@ template <class Traits> inline void AttributeTable<Traits>::insert_always(Key k, Value v) { - /*std::cout << "Insert " << k << " in v of size " + /*std::cout << "Insert " << k << " in v of size " << size_ << " " << map_ << " " << k.get_index() << std::endl;*/ IMP_assert(k != Key(), "Can't insert default key"); IMP_check(Traits::get_is_valid(v), - "Trying to insert invalid value for attribute " + "Trying to insert invalid value for attribute " << v << " into attribute " << k, ValueException); - if (size_ <= k.get_index()) { - boost::scoped_array<Value> old; - swap(old, map_); - realloc(size_, old, k.get_index()+1); - } + typename Map::size_type val + =static_cast<typename Map::size_type>(k.get_index()); + IMP_assert(val <1000, "Bad key index"); + map_.resize(std::max(map_.size(), + val+1), + Traits::get_invalid()); map_[k.get_index()]= v; } @@ -273,9 +248,9 @@ template <class Traits> inline void AttributeTable<Traits>::show(std::ostream &out, - const char *prefix) const + const std::string prefix) const { - for (unsigned int i=0; i< size_; ++i) { + for (unsigned int i=0; i< map_.size(); ++i) { if (Traits::get_is_valid(map_[i])) { out << prefix; out << Key(i).get_string() << ": "; @@ -287,12 +262,12 @@ template <class Traits> -inline std::vector<typename Traits::Key> +inline std::vector<typename Traits::Key> AttributeTable<Traits>::get_keys() const { std::vector<Key> ret(attribute_keys_begin(), attribute_keys_end()); return ret; -} +} } // namespace internal Index: kernel/include/IMP/decorators/macros.h =================================================================== --- kernel/include/IMP/decorators/macros.h (revision 669) +++ kernel/include/IMP/decorators/macros.h (working copy) @@ -262,4 +262,7 @@ name##_data_.initialize(); +#define IMP_ATOM_TYPE_INDEX 8974343 +#define IMP_RESIDUE_TYPE_INDEX 90784334 + #endif /* __IMP_DECORATOR_MACROS_H */ Index: kernel/include/IMP/decorators/ResidueDecorator.h =================================================================== --- kernel/include/IMP/decorators/ResidueDecorator.h (revision 669) +++ kernel/include/IMP/decorators/ResidueDecorator.h (working copy) @@ -17,8 +17,7 @@ namespace IMP { -struct ResidueTypeTag{}; -IMP_DECLARE_KEY_TYPE(ResidueType, ResidueTypeTag); +IMP_DECLARE_KEY_TYPE(ResidueType, IMP_RESIDUE_TYPE_INDEX); //! A decorator for a residue. Index: kernel/include/IMP/decorators/AtomDecorator.h =================================================================== --- kernel/include/IMP/decorators/AtomDecorator.h (revision 669) +++ kernel/include/IMP/decorators/AtomDecorator.h (working copy) @@ -13,6 +13,7 @@ #include "../Model.h" #include "utility.h" #include "XYZDecorator.h" +#include "macros.h" #include <vector> #include <deque> @@ -20,8 +21,7 @@ namespace IMP { -struct AtomTypeTag{}; -IMP_DECLARE_KEY_TYPE(AtomType, AtomTypeTag); +IMP_DECLARE_KEY_TYPE(AtomType, IMP_ATOM_TYPE_INDEX); //! A decorator for a particle representing an atom. Index: kernel/include/IMP/base_types.h =================================================================== --- kernel/include/IMP/base_types.h (revision 669) +++ kernel/include/IMP/base_types.h (working copy) @@ -14,7 +14,6 @@ #include <string> #include <vector> -#include <map> namespace IMP { @@ -63,7 +62,7 @@ class Particle; //! A class which is used for representing collections of particles /** - We need this to have a uniform return type for python. + We need this to have a uniform return type for python. \todo It would be nice to use internal::Vector instead, but that is not as pretty for Python. */ @@ -77,13 +76,13 @@ class Particle; //! The type used to identify float attributes in the Particles -IMP_DECLARE_KEY_TYPE(FloatKey, Float); +IMP_DECLARE_KEY_TYPE(FloatKey, 0); //! The type used to identify int attributes in the Particles -IMP_DECLARE_KEY_TYPE(IntKey, Int); +IMP_DECLARE_KEY_TYPE(IntKey, 1); //! The type used to identify string attributes in the Particles -IMP_DECLARE_KEY_TYPE(StringKey, String); +IMP_DECLARE_KEY_TYPE(StringKey, 2); //! The type used to identify a particle attribute in the Particles -IMP_DECLARE_KEY_TYPE(ParticleKey, Particle*); +IMP_DECLARE_KEY_TYPE(ParticleKey, 3); } // namespace IMP

1 0

Quick cleanup of IMP.i
by Daniel Russel 26 Aug '08

26 Aug '08

Here is a simple patch to IMP.i which adds macros to handle the thisowns and adds a few missed ones. Index: kernel/include/IMP/optimizers/MonteCarlo.h =================================================================== --- kernel/include/IMP/optimizers/MonteCarlo.h (revision 669) +++ 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/Mover.h =================================================================== --- kernel/include/IMP/optimizers/Mover.h (revision 669) +++ 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/pyext/IMP_macros.i =================================================================== --- kernel/pyext/IMP_macros.i (revision 669) +++ kernel/pyext/IMP_macros.i (working copy) @@ -1,6 +1,38 @@ #define IMPDLLEXPORT #define IMPDLLLOCAL +%define IMP_OWN_FIRST_CONSTRUCTOR(Ucname) +%pythonprepend Ucname::Ucname %{ + args[0].thisown=0 +%} +%enddef + +%define IMP_OWN_FIRST_SECOND_CONSTRUCTOR(Ucname) +%pythonprepend Ucname::Ucname %{ + args[0].thisown=0 + args[1].thisown=0 +%} +%enddef + +%define IMP_SET_OBJECT(Ucname, methodname) +%pythonprepend Ucname::methodname %{ + args[1].thisown=0 +%} +%enddef + +%define IMP_ADD_OBJECT(Ucname, methodname) +%pythonprepend Ucname::methodname %{ + args[1].thisown=0 +%} +%enddef + +%define IMP_ADD_OBJECTS(Ucname, methodname) +%pythonprepend Ucname::methodname %{ + for p in args[1]: + p.thisown=0 +%} +%enddef + /* Add additional IMP_CONTAINER methods for scripting languages */ %define IMP_CONTAINER_SWIG(type, Ucname, lcname) %extend type { @@ -9,6 +41,8 @@ return ret; } } +IMP_ADD_OBJECT(type, add_##lcname) +IMP_ADD_OBJECTS(type, add_##lcname##s) %enddef %include "IMP/macros.h" Index: kernel/pyext/IMP.i =================================================================== --- kernel/pyext/IMP.i (revision 669) +++ kernel/pyext/IMP.i (working copy) @@ -27,94 +27,30 @@ %} namespace IMP { - %pythonprepend Model::add_restraint %{ - args[1].thisown=0 - %} - %pythonprepend Model::add_score_state %{ - args[1].thisown=0 - %} - %pythonprepend Optimizer::add_optimizer_state %{ - args[1].thisown=0 - %} - %pythonprepend RestraintSet::add_restraint %{ - args[1].thisown=0 - %} - %pythonprepend NonbondedListScoreState::add_bonded_list %{ - args[1].thisown=0 - %} - %pythonprepend DistanceRestraint::DistanceRestraint %{ - args[0].thisown=0 - %} - %pythonprepend AngleRestraint::AngleRestraint %{ - args[0].thisown=0 - %} - %pythonprepend DihedralRestraint::DihedralRestraint %{ - args[0].thisown=0 - %} - %pythonprepend TorusRestraint::TorusRestraint %{ - args[3].thisown=0 - %} - %pythonprepend NonbondedRestraint::NonbondedRestraint %{ - args[0].thisown=0 - %} - %pythonprepend BondDecoratorRestraint::BondDecoratorRestraint %{ - args[0].thisown=0 - %} - %pythonprepend SingletonListRestraint::SingletonListRestraint %{ - args[0].thisown=0 - %} - %pythonprepend PairListRestraint::PairListRestraint %{ - args[0].thisown=0 - %} - %pythonprepend TripletChainRestraint::TripletChainRestraint %{ - args[0].thisown=0 - %} - %pythonprepend PairChainRestraint::PairChainRestraint %{ - args[0].thisown=0 - %} - %pythonprepend ConnectivityRestraint::ConnectivityRestraint %{ - args[0].thisown=0 - %} - %pythonprepend DistancePairScore::DistancePairScore %{ - args[0].thisown=0 - %} - %pythonprepend TransformedDistancePairScore::TransformedDistancePairScore %{ - args[0].thisown=0 - %} - %pythonprepend BondCoverPairScore::BondCoverPairScore %{ - args[0].thisown=0 - %} - %pythonprepend SphereDistancePairScore::SphereDistancePairScore %{ - args[0].thisown=0 - %} - %pythonprepend RefineOncePairScore::RefineOncePairScore %{ - args[0].thisown=0 - args[1].thisown=0 - %} - %pythonprepend DistanceToSingletonScore::DistanceToSingletonScore %{ - args[0].thisown=0 - %} - %pythonprepend AttributeSingletonScore::AttributeSingletonScore %{ - args[0].thisown=0 - %} - %pythonprepend TunnelSingletonScore::TunnelSingletonScore %{ - args[0].thisown=0 - %} - %pythonprepend AngleTripletScore::AngleTripletScore %{ - args[0].thisown=0 - %} - %pythonprepend MonteCarlo::add_mover %{ - args[1].thisown=0 - %} - %pythonprepend MonteCarlo::set_local_optimizer %{ - args[1].thisown=0 - %} - %pythonprepend VRMLLogOptimizerState::add_particle_refiner %{ - args[1].thisown=0 - %} - %pythonprepend TypedPairScore::set_pair_score %{ - args[1].thisown=0 - %} + // need to special case particle so can't add this to macro + IMP_OWN_FIRST_CONSTRUCTOR(DistanceRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(AngleRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(DihedralRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(TorusRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(NonbondedRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(BondDecoratorRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(SingletonListRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(PairListRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(TripletChainRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(PairChainRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(ConnectivityRestraint) + IMP_OWN_FIRST_CONSTRUCTOR(DistancePairScore) + IMP_OWN_FIRST_CONSTRUCTOR(TransformedDistancePairScore) + IMP_OWN_FIRST_CONSTRUCTOR(BondCoverPairScore) + IMP_OWN_FIRST_CONSTRUCTOR(SphereDistancePairScore) + IMP_OWN_FIRST_SECOND_CONSTRUCTOR(RefineOncePairScore) + IMP_OWN_FIRST_CONSTRUCTOR(DistanceToSingletonScore) + IMP_OWN_FIRST_CONSTRUCTOR(AttributeSingletonScore) + IMP_OWN_FIRST_CONSTRUCTOR(TunnelSingletonScore) + IMP_OWN_FIRST_CONSTRUCTOR(AngleTripletScore) + IMP_SET_OBJECT(MonteCarlo, set_local_optimizer) + IMP_SET_OBJECT(TypedPairScore, set_pair_score) + %pythonprepend Particle::get_value %{ check_particle(args[0], args[1]) %} @@ -144,10 +80,19 @@ %} - IMP_CONTAINER_SWIG(Model, Particle, particle); + // special case since particles are ref counted + %extend Model { + Particles get_particles() const { + IMP::Particles ret(self->particles_begin(), self->particles_end()); + return ret; + } + } 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();" @@ -168,6 +113,7 @@ %feature("director") IMP::TripletScore; %feature("director") IMP::Optimizer; %feature("director") IMP::ParticleRefiner; +%feature("director") IMP::Mover; %include "IMP/Key.h" %include "IMP/Object.h" Index: SConstruct =================================================================== --- SConstruct (revision 669) +++ 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)

1 0

One final attempt to stir discussion before I give up on IMP svn
by Daniel Russel 21 Aug '08

21 Aug '08

The current IMP rcs setup is completely idiotic. As a result, in order to be able to actually track the changes I make to the IMP kernel I have to maintain my own fork and then move patches back to the main svn. I am tired of doing this and am somewhat inclined just to maintain my own (public) fork and forget about contributing back to the main IMP svn. As there are occasionally a few global changes, reconciling things will get increasingly hard over time. Restricting who can commit to svn offers no benefits over open submissions. - code review as it now stands appears to simply consist of checking formatting of the source files. Such can changes can be more easily done as a modification of already checked in code. That would have the added advantage of not having svn report everything as a conflict (svn gets confused by my changes being committed from another working copy). - No one is counting on IMP svn being constantly in a working state. People wait to check things out until they know they have some time to clean things up. As a result, there the invariant that all versions of the repository work doesn't buy us anything over ensuring that it works almost all of the time. The latter would make it easier to submit patches since you don't have to ensure all files are submitted at once and you can more easily check that your patch is OK by checking IMP out and building it. Unless the automatic test scripts for IMP are a mess it should be easy enough to have them create a "last stable" branch in svn upon successful completion of the tests. - svn is a revision control system. It allows us to roll back changes we don't like after they are committed. That is the point. And the current system has many disadvantages: - submitting patches takes a reasonable amount of effort. As a result minor changes to documentation and things which result in IMP overall being a nicer experience never get submitted. In addition, it makes people wait longer to submit things making the eventually submission that much more complicated and reviewing of changes that much more involved. - since it can take days (sometimes weeks if Ben is away) to get even the simplest patches committed, it takes a lot of work to make focused patches. Our time is better spent elsewhere. - IMP is supposed to be a collaborative effort. It is hard enough to get people to share their code without added hurdles. As it is, I seem to be the only one inclined to go through the effort. - we can't actually use svn as a revision control system. If everything in in sync most of the time, then once can easily try speculative changes and then revert them if they don't work. Since I can't do this, I have to maintain my own svn repository with a copy of IMP. - we don't have any easy record of who last changed each file without reading the whole log making it harder to know who to ask about changes and breakages. Creating yet another project in IMP svn doesn't solve these problems and results in another complication for people who want to use IMP since they always have to think about which library to get it from/which namespace to use.

4 19

returning references (regarding get_foo.begin(), .end())
by Daniel Russel 20 Aug '08

20 Aug '08

As a general guideline in IMP get_ methods never returns a reference and so doing things like get_foo.begin() is never a safe thing. Any such method which happens to return a reference is just an optimization for efficiency and might change later. The only way I know to nicely catch the problem that Keren encountered is to turn on iterator debugging in gcc, but that changes the C++ ABI and so is a bit awkward if you link against any other C++ libraries.

1 0

get_interacting_particles patch
by Daniel Russel 19 Aug '08

19 Aug '08

Here is a new patch which includes the get_interacting_particles code. It is not implemented for NonbondedRestraint as that did not appear to be currently needed and required a few changes to the nonbonded lists and so was not entirely trivial. The other changes from the last version of the patch are: - TripletChainRestraint now only has one chain (as with PairChainRestraint) - There is a new PairScore RefineHighPairScore which refines high scoring pairs for use in hierarchical collision detection.

2 1

NBL cleanup and other small patches
by Daniel Russel 18 Aug '08

18 Aug '08

Frido was complaining that the nbl was too slow, so I thought I would check in my local changes. The big thing is that the fastest implementation I have uses the CGAL library. It is part of fedora extras so getting it is not that big a deal (and I think there are installers for macs). The nbl code is controlled by the macro IMP_USE_CGAL so everything works fine without it. I have a patched scons init file which adds support (but does not do a very good job since I don't know how to write link tests and such). First some small patches: The first adds a constructor to decorators which assumes that the particles have the needed fields (this is checked if debug code is used). We should add checks to the python interface eventually. The second patches the init file as above. It also adds an argument extraroot which adds to the include search path and the lib path (and to the runtime link path since it is needed to make the boost tests pass currently). There are probably lots of ways to make this better as I don't really know scons. Also, as a general point, the SConscript files which build things should all be changed to add the local include and lib paths to the front so that local headers mask globally installed headers. The third prints a nicer error message if you forget to set the model before CG::optimize (this is independent of the other patches) The fourth uses the decorator constructors in the bond decorators. The fifth adds set_particles to the the list restraint I added before (this is independent of anything else). And the fifth patches the nonbonded lists. There is now just AllNonbondedListScoreState and BipartiteNonbondedListScore state with set_algorithm methods (and arguments to the constructor) to determine which algorithm is used. The default is either grid or CGAL if the latter is available. The other interface change is that the nbl now takes the distance cutoff in its constructor (rather than the cutoff being passed to the NBRestraint as before). This simplified things a great deal. The patch also brings in a BondedList for chains which should be faster than using the bond decorators. Some timing results for the various nbls follow. The results are normalized by the number of particles involved. Nothing: (just run the python code to set up the particles) For 50 particles got 3.24961185455 For 100 particles got 3.17380809784 For 1000 particles got 3.14385008812 For 10000 particles got 3.13716006279 For 100000 particles got 3.16522884369 Bbox- (CGAL) For 50 particles got 3.32897806168 For 100 particles got 3.29577898979 For 1000 particles got 3.42112493515 For 10000 particles got 4.14875411987 For 100000 particles got 5.81269311905 Grid For 50 particles got 3.63301706314 For 100 particles got 3.84657287598 For 1000 particles got 6.80274581909 For 10000 particles got 11.4696931839 Quadratic For 50 particles got 3.39585018158 For 100 particles got 3.51785087585 For 1000 particles got 6.79182696342 For 10000 particles got 44.5892150402 Index: kernel/include/IMP/decorators/macros.h =================================================================== --- kernel/include/IMP/decorators/macros.h (revision 524) +++ kernel/include/IMP/decorators/macros.h (working copy) @@ -44,15 +44,17 @@ add_required; \ } \ friend class DecoratorBase; \ - Name(::IMP::Particle* p): Parent(p) { \ - IMP_assert(has_required_attributes(p), \ - "This is not a particle of type " \ - << #Name << *p); \ - } \ public: \ typedef Name This; \ - /** The default constructor. This is used as a null value */ \ + /** \short The default constructor. This is used as a null value */ \ Name(): Parent(){} \ + /** \short Construct from a Particle which has all needed attributes */\ + Name(::IMP::Particle *p): Parent(p){ \ + if (!decorator_keys_initialized_) decorator_initialize_static_data();\ + IMP_assert(has_required_attributes(p), \ + "Particle missing required attributes for decorator" \ + << #Name << *p << std::endl); \ + } \ /** Add the necessary attributes to p and return a decorator. */ \ static Name create(::IMP::Particle *p) { \ return IMP::DecoratorBase::create<Name>(p); \ Index: tools/__init__.py =================================================================== --- tools/__init__.py (revision 524) +++ tools/__init__.py (working copy) @@ -53,7 +53,10 @@ def _add_release_flags(env): """Add compiler flags for release builds, if requested""" if env.get('release', False): - env.Append(CPPDEFINES='NDEBUG') + env.Append(CPPDEFINES=['NDEBUG']) + env.Append(CCFLAGS=["-O3"]) + else: + env.Append(CCFLAGS=["-g"]) def CheckGNUHash(context): """Disable GNU_HASH-style linking (if found) for backwards compatibility""" @@ -186,16 +189,30 @@ pass env.Prepend(SCANNERS = _SWIGScanner) if env['CC'] == 'gcc': - env.Append(CCFLAGS="-Wall -g -O3") + env.Append(CCFLAGS=["-Wall"]) _add_release_flags(env) sys = platform.system() if sys == 'SunOS': # Find locally-installed libraries in /usr/local (e.g. for SWIG) env['ENV']['LD_LIBRARY_PATH'] = '/usr/local/lib' - # Make Modeller exetype variable available: - if os.environ.has_key('EXECUTABLE_TYPESVN'): - env['ENV']['EXECUTABLE_TYPESVN'] = os.environ['EXECUTABLE_TYPESVN'] + + if (ARGUMENTS.get('extraroot', '') != ''): + env.Append(CPPPATH=['$extraroot/include']) + env.Append(LIBPATH=['$extraroot/lib']) + env.Append(LINKFLAGS=['-Wl,-rpath,$extraroot/lib']) + env.Append(LD_LIBRARY_PATH=['$extraroot/lib']) + env['EXTRA_LDPATH']=ARGUMENTS.get('extraroot', '')+'/lib' + env['EXTRA_PYTHONPATH']=ARGUMENTS.get('extraroot', '')+'/lib/python2.5/site-packages' + + + if (ARGUMENTS.get('usecgal')): + # should do tests here to make sure it is usable, but I don't know how + env.Append(CPPDEFINES=['IMP_USE_CGAL']) + env.Append(LIBS= ['CGAL']) + if env['CC'] == 'gcc': + env.Append(CCFLAGS=["-frounding-math"]) + # Set empty variables in case the Modeller check fails: for mod in ('MODPY', 'EXETYPE'): env['MODELLER_' + mod] = '' @@ -309,6 +326,8 @@ # Remove options that don't work with C++ code: if cplusplus: opt = opt.replace("-Wstrict-prototypes", "") + if cplusplus: + opt = opt.replace("-O2", "") e.Replace(CC=cc, CXX=cxx, LDMODULESUFFIX=so) e.Replace(CPPFLAGS=basecflags.split() + opt.split()) @@ -372,3 +391,10 @@ opts.Add(BoolOption('release', 'Disable most runtime checks (e.g. for releases)', False)) + opts.Add(PathOption('extraroot', + 'Extra place to search for includes and libs', + '', PathOption.PathAccept)) + # I would like to have a check, but it doesn't seem to turn it off it no value is passed + #PathOption.PathIsDir)) + opts.Add(BoolOption('usecgal', + 'Use the CGAL library for some algorithms', False)) Index: kernel/src/optimizers/ConjugateGradients.cpp =================================================================== --- kernel/src/optimizers/ConjugateGradients.cpp (revision 524) +++ kernel/src/optimizers/ConjugateGradients.cpp (working copy) @@ -214,6 +214,9 @@ Float ConjugateGradients::optimize(unsigned int max_steps) { + IMP_check(get_model(), + "Must set the model on the optimizer before optimizing", + ValueException("Must set the model before optimizing")); std::vector<Float> x, dx; int i; //ModelData* model_data = get_model()->get_model_data(); Index: kernel/src/decorators/bond_decorators.cpp =================================================================== --- kernel/src/decorators/bond_decorators.cpp (revision 524) +++ kernel/src/decorators/bond_decorators.cpp (working copy) @@ -77,7 +77,7 @@ { Particle *p= internal::graph_connect(a.get_particle(), b.get_particle(), internal::bond_graph_data_); - BondDecorator bd= BondDecorator::cast(p); + BondDecorator bd(p); bd.set_type(t); return bd; } Index: kernel/include/IMP/decorators/bond_decorators.h =================================================================== --- kernel/include/IMP/decorators/bond_decorators.h (revision 524) +++ kernel/include/IMP/decorators/bond_decorators.h (working copy) @@ -102,7 +102,7 @@ BondDecorator get_bond(unsigned int i) const { Particle *p= graph_get_edge(get_particle(), i, internal::bond_graph_data_); - return BondDecorator::cast(p); + return BondDecorator(p); } //! Get a BondedDecorator of the ith bonded particle @@ -117,7 +117,7 @@ */ BondedDecorator get_bonded(unsigned int i) const { Particle *p= graph_get_edge(get_particle(), i, internal::bond_graph_data_); - BondDecorator bd= BondDecorator::cast(p); + BondDecorator bd(p); if (bd.get_bonded(0) == *this) return bd.get_bonded(1); else return bd.get_bonded(0); } @@ -130,7 +130,7 @@ { Particle *p= graph_get_node(get_particle(), i, internal::bond_graph_data_); - return BondedDecorator::cast(p); + return BondedDecorator(p); } //! Connect the two wrapped particles by a bond. Index: kernel/include/IMP/restraints/SingletonListRestraint.h =================================================================== --- kernel/include/IMP/restraints/SingletonListRestraint.h (revision 524) +++ kernel/include/IMP/restraints/SingletonListRestraint.h (working copy) @@ -38,6 +38,7 @@ using Restraint::add_particles; using Restraint::clear_particles; + using Restraint::set_particles; protected: internal::ObjectPointer<SingletonScore, true> ss_; Index: kernel/include/IMP/restraints/NonbondedRestraint.h =================================================================== --- kernel/include/IMP/restraints/NonbondedRestraint.h (revision 524) +++ kernel/include/IMP/restraints/NonbondedRestraint.h (working copy) @@ -34,10 +34,8 @@ object is deleted upon destruction. \param[in] nbl The non-bonded list to use to get the list of particles. - \param[in] max_dist Pairs beyond this distance may be dropped. */ - NonbondedRestraint(PairScore *ps, NonbondedListScoreState *nbl, - Float max_dist= std::numeric_limits<Float>::max()); + NonbondedRestraint(PairScore *ps, NonbondedListScoreState *nbl); virtual ~NonbondedRestraint(){} IMP_RESTRAINT(internal::kernel_version_info) @@ -45,7 +43,6 @@ protected: NonbondedListScoreState *nbl_; internal::ObjectPointer<PairScore, true> sf_; - Float max_dist_; }; } // namespace IMP Index: kernel/src/restraints/NonbondedRestraint.cpp =================================================================== --- kernel/src/restraints/NonbondedRestraint.cpp (revision 524) +++ kernel/src/restraints/NonbondedRestraint.cpp (working copy) @@ -17,9 +17,8 @@ { NonbondedRestraint::NonbondedRestraint(PairScore *ps, - NonbondedListScoreState *nbl, - Float md) : nbl_(nbl), sf_(ps), - max_dist_(md) + NonbondedListScoreState *nbl) + : nbl_(nbl), sf_(ps) { } @@ -31,12 +30,12 @@ IMP_CHECK_OBJECT(nbl_); Float score=0; IMP_LOG(VERBOSE, "Nonbonded restraint on " - << std::distance(nbl_->nonbonded_begin(max_dist_), - nbl_->nonbonded_end(max_dist_)) + << std::distance(nbl_->nonbonded_begin(), + nbl_->nonbonded_end()) << " pairs" << std::endl); for (NonbondedListScoreState::NonbondedIterator it - = nbl_->nonbonded_begin(max_dist_); - it != nbl_->nonbonded_end(max_dist_); ++it) { + = nbl_->nonbonded_begin(); + it != nbl_->nonbonded_end(); ++it) { float thisscore = sf_->evaluate(it->first, it->second, accum); if (thisscore != 0) { IMP_LOG(VERBOSE, "Pair " << it->first->get_index() Index: kernel/include/IMP/score_states/QuadraticNonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/QuadraticNonbondedListScoreState.h (revision 524) +++ kernel/include/IMP/score_states/QuadraticNonbondedListScoreState.h (working copy) @@ -1,51 +0,0 @@ -/** - * \file QuadraticNonbondedListScoreState.h - * \brief Allow iteration through pairs of a set of spheres. - * - * Copyright 2007-8 Sali Lab. All rights reserved. - */ - -#ifndef __IMP_QUADRATIC_NONBONDED_LIST_SCORE_STATE_H -#define __IMP_QUADRATIC_NONBONDED_LIST_SCORE_STATE_H - -#include "NonbondedListScoreState.h" -#include "../internal/kernel_version_info.h" -#include "MaxChangeScoreState.h" - -#include <vector> -#include <limits> - -namespace IMP -{ - -//! A base class for nonbonded lists which test all pairs of particles -/** This should not be used by end users. But since it needs to be in the - inheritance hierarchy, it should be in the IMP namespace. - */ -class IMPDLLEXPORT QuadraticNonbondedListScoreState: - public NonbondedListScoreState -{ - typedef NonbondedListScoreState P; - internal::ObjectPointer<MaxChangeScoreState, true> mc_; - float slack_; - protected: - void handle_nbl_pair( Particle *a, Particle *b, float d); - const Particles &get_particles() const { - return mc_->get_particles(); - } - void set_particles(const Particles &ps) { - P::clear_nbl(); - mc_->clear_particles(); - mc_->add_particles(ps); - } - - QuadraticNonbondedListScoreState(FloatKey radius); - ~QuadraticNonbondedListScoreState(); - - public: - void update(); -}; - -} // namespace IMP - -#endif /* __IMP_QUADRATIC_NONBONDED_LIST_SCORE_STATE_H */ Index: kernel/include/IMP/score_states/BondDecoratorListScoreState.h =================================================================== --- kernel/include/IMP/score_states/BondDecoratorListScoreState.h (revision 524) +++ kernel/include/IMP/score_states/BondDecoratorListScoreState.h (working copy) @@ -8,10 +8,11 @@ #ifndef __IMP_BOND_DECORATOR_LIST_SCORE_STATE_H #define __IMP_BOND_DECORATOR_LIST_SCORE_STATE_H -#include <vector> #include "BondedListScoreState.h" #include "../decorators/bond_decorators.h" +#include <vector> + namespace IMP { @@ -25,6 +26,7 @@ */ class IMPDLLEXPORT BondDecoratorListScoreState: public BondedListScoreState { + // This is only used to provide the iterators. Probably should be lazy. std::vector<BondDecorator> bonds_; Particles ps_; public: Index: kernel/include/IMP/score_states/NonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/NonbondedListScoreState.h (revision 524) +++ kernel/include/IMP/score_states/NonbondedListScoreState.h (working copy) @@ -11,6 +11,7 @@ #include "../ScoreState.h" #include "../internal/ParticleGrid.h" #include "../internal/kernel_version_info.h" +#include "../decorators/XYZDecorator.h" #include "BondedListScoreState.h" #include <vector> @@ -19,34 +20,52 @@ namespace IMP { +namespace internal +{ +class NBLAddPairIfNonbonded; +class NBLAddIfNonbonded; +} + typedef std::vector<BondedListScoreState*> BondedListScoreStates; -//! A base class for classes that maintain a list of non-bonded pairs. -/** +//! An abstract base class for classes that maintain a list of non-bonded pairs. +/** \note If no value for the radius key is specified, all radii are + considered to be zero. */ class IMPDLLEXPORT NonbondedListScoreState: public ScoreState { private: FloatKey rk_; - -protected: + // the distance cutoff to count a pair as adjacent + float cutoff_; + /* How much to add to the size of particles to allow them to move + without rebuilding*/ + Float slack_; + /* An estimate of what the slack should be next time they are recomputed*/ + Float next_slack_; + int num_steps_; + bool nbl_is_valid_; // made protected for debugging code, do not use otherwise typedef std::vector<std::pair<Particle*, Particle*> > NBL; NBL nbl_; - float last_cutoff_; +protected: + + Float get_slack() const {return slack_;} + Float get_cutoff() const {return cutoff_;} + unsigned int size_nbl() const {return nbl_.size();} //! rebuild the nonbonded list - /** \internal + /** This is used by classes which inherit from this to rebuild the NBL + + \internal */ - virtual void rebuild_nbl(float cut)=0; + virtual void rebuild_nbl()=0; - void clear_nbl() { - last_cutoff_=-1; - nbl_.clear(); - } + + // Return true if two particles are bonded bool are_bonded(Particle *a, Particle *b) const { IMP_assert(a->get_is_active() && b->get_is_active(), "Inactive particles should have been removed" @@ -60,27 +79,13 @@ return false; } - struct AddToNBL; - friend struct AddToNBL; + friend struct internal::NBLAddPairIfNonbonded; + friend struct internal::NBLAddIfNonbonded; - // these should be two separate classes, but they are not - struct AddToNBL{ - NonbondedListScoreState *state_; - Particle *p_; - AddToNBL(NonbondedListScoreState *s, Particle *p): state_(s), - p_(p){} - void operator()(Particle *p) { - operator()(p_, p); - } - void operator()(Particle *a, Particle *b) { - state_->add_to_nbl(a,b); - } - }; - //! tell the bonded lists what particles to pay attention to - void propagate_particles(const Particles &ps); - void add_to_nbl(Particle *a, Particle *b) { + // Add the pair to the nbl if the particles are not bonded + void add_if_nonbonded(Particle *a, Particle *b) { IMP_assert(a->get_is_active() && b->get_is_active(), "Inactive particles should have been stripped"); @@ -95,6 +100,28 @@ } } + // Check if the bounding boxes of the spheres overlap + void add_if_box_overlap(Particle *a, Particle *b) { + XYZDecorator da(a); + XYZDecorator db(b); + float ra= get_radius(a); + float rb= get_radius(b); + for (unsigned int i=0; i< 3; ++i) { + float delta=std::abs(da.get_coordinate(i) - db.get_coordinate(i)); + if (delta -ra -rb > cutoff_+slack_) { + IMP_LOG(VERBOSE, "Pair " << a->get_index() + << " and " << b->get_index() << " rejected on coordinate " + << i << " " << da.get_coordinate(i) << " and " + << db.get_coordinate(i) << " " << ra << " " + << rb << " " << slack_ << " " << slack_ << std::endl); + return ; + } + } + IMP_LOG(VERBOSE, "Adding pair " << a->get_index() + << " and " << b->get_index() << std::endl); + add_if_nonbonded(a, b); + } + Float get_radius(Particle *a) const { if (rk_ != FloatKey() && a->has_attribute(rk_)) { return a->get_value(rk_); @@ -102,11 +129,26 @@ return 0; } } + // return true if the nbl was invalidated by a move of mc + bool update(Float mc); + // Return a list of all the particles in in with a radius field + Particles particles_with_radius(const Particles &in) const; + + + // Return true if the nbl is valid + bool get_nbl_is_valid() const {return nbl_is_valid_;} + + void set_nbl_is_valid(bool tf) { + nbl_is_valid_= tf; + if (!nbl_is_valid_) { + nbl_.clear(); + } + } public: /** */ - NonbondedListScoreState(FloatKey rk); + NonbondedListScoreState(Float cutoff, FloatKey rk); FloatKey get_radius_key() const {return rk_;} void set_radius_key(FloatKey rk) {rk_=rk;} @@ -118,48 +160,56 @@ typedef BondedListScoreStateIterator BondedListIterator; typedef BondedListScoreStateConstIterator BondedListConstIterator; - IMP_SCORE_STATE(internal::kernel_version_info) + virtual void show(std::ostream &out=std::cout) const; + virtual IMP::VersionInfo get_version_info() const { + return internal::kernel_version_info; + } //! An iterator through nonbonded particles - /** The value type is an std::pair<Particle*, Particle*> + /** The value type is an ParticlePair. */ typedef NBL::const_iterator NonbondedIterator; //! This iterates through the pairs of non-bonded particles - /** \param[in] cutoff The state may ignore pairs which are futher - apart than the cutoff. - \note that this is highly unsafe and iteration can only be - done once at a time. I will fix that eventually. - - \note The distance cutoff is the l2 norm between the 3D coordinates - of the Particles. It ignore any size that may be associated with - the particles. - */ - NonbondedIterator nonbonded_begin(Float cutoff - =std::numeric_limits<Float>::max()) { - IMP_assert(last_cutoff_== cutoff || last_cutoff_==-1, - "Bad things are happening with the iterators in " - << "NonbondedListScoreState"); - if (last_cutoff_ < cutoff) { - IMP_LOG(VERBOSE, "Rebuilding NBL cutoff " << cutoff << std::endl); - clear_nbl(); - rebuild_nbl(cutoff); - last_cutoff_=cutoff; - } + NonbondedIterator nonbonded_begin() const { + IMP_check(get_nbl_is_valid(), "Must call update first", + ValueException("Must call update first")); return nbl_.begin(); } - NonbondedIterator nonbonded_end(Float cutoff - =std::numeric_limits<Float>::max()) { - if (last_cutoff_ < cutoff) { - IMP_LOG(VERBOSE, "Rebuilding NBL cutoff " << cutoff << std::endl); - clear_nbl(); - rebuild_nbl(cutoff); - last_cutoff_=cutoff; - } + NonbondedIterator nonbonded_end() const { + IMP_check(get_nbl_is_valid(), "Must call update first", + ValueException("Must call update first")); return nbl_.end(); } + + + unsigned int number_of_nonbonded() const { + return nbl_.size(); + } }; + +namespace internal +{ + struct NBLAddPairIfNonbonded{ + NonbondedListScoreState *state_; + NBLAddPairIfNonbonded(NonbondedListScoreState *s): state_(s){} + void operator()(Particle *a, Particle *b) { + state_->add_if_nonbonded(a,b); + } + }; + + struct NBLAddIfNonbonded{ + NonbondedListScoreState *state_; + Particle *p_; + NBLAddIfNonbonded(NonbondedListScoreState *s, Particle *p): state_(s), + p_(p){} + void operator()(Particle *p) { + state_->add_if_nonbonded(p_, p); + } + }; +} + } // namespace IMP #endif /* __IMP_NONBONDED_LIST_SCORE_STATE_H */ Index: kernel/include/IMP/score_states/SConscript =================================================================== --- kernel/include/IMP/score_states/SConscript (revision 524) +++ kernel/include/IMP/score_states/SConscript (working copy) @@ -1,14 +1,15 @@ +Import('env') import os.path -Import('env') - -files = ['BondedListScoreState.h', 'NonbondedListScoreState.h', - 'BipartiteNonbondedListScoreState.h', 'MaxChangeScoreState.h', - 'BondDecoratorListScoreState.h', 'AllNonbondedListScoreState.h', - 'QuadraticBipartiteNonbondedListScoreState.h', - 'QuadraticAllNonbondedListScoreState.h', - 'QuadraticNonbondedListScoreState.h', 'GravityCenterScoreState.h'] - -# Install the include files: -includedir = os.path.join(env['includedir'], 'IMP', 'score_states') +files=[ + 'AllNonbondedListScoreState.h', + 'BipartiteNonbondedListScoreState.h', + 'BondDecoratorListScoreState.h', + 'BondedListScoreState.h', + 'ChainBondedListScoreState.h', + 'GravityCenterScoreState.h', + 'MaxChangeScoreState.h', + 'NonbondedListScoreState.h', + ] +includedir = os.path.join(env['includedir'], 'IMP', 'score_states' ) inst = env.Install(includedir, files) env.Alias('install', inst) Index: kernel/include/IMP/score_states/QuadraticBipartiteNonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/QuadraticBipartiteNonbondedListScoreState.h (revision 524) +++ kernel/include/IMP/score_states/QuadraticBipartiteNonbondedListScoreState.h (working copy) @@ -1,49 +0,0 @@ -/** - * \file QuadraticBipartiteNonbondedListScoreState.h - * \brief Allow iteration through pairs of a set of atoms. - * - * Copyright 2007-8 Sali Lab. All rights reserved. - */ - -#ifndef __IMP_QUADRATIC_BIPARTITE_NONBONDED_LIST_SCORE_STATE_H -#define __IMP_QUADRATIC_BIPARTITE_NONBONDED_LIST_SCORE_STATE_H - -#include "QuadraticNonbondedListScoreState.h" -#include "../internal/ParticleGrid.h" -#include "../internal/kernel_version_info.h" - -#include <vector> -#include <limits> - -namespace IMP -{ - -//! This class maintains a list of non-bonded pairs between two sets. -/** The class works roughly like the NonbondedListScoreState except - only pairs where one particle is taken from each set are returned. - - \note QuadraticBipartiteNonbondedListScoreState is basically an - implementation detail for performance analysis and should not - be used by end users. - */ -class IMPDLLEXPORT QuadraticBipartiteNonbondedListScoreState: - public QuadraticNonbondedListScoreState -{ - typedef QuadraticNonbondedListScoreState P; - unsigned int na_; - - virtual void rebuild_nbl(float cut); -public: - QuadraticBipartiteNonbondedListScoreState(FloatKey rk, - const Particles &ps0, - const Particles &ps1); - QuadraticBipartiteNonbondedListScoreState(FloatKey rk); - - IMP_SCORE_STATE(internal::kernel_version_info) - - void set_particles(const Particles &ps0, const Particles &ps1); -}; - -} // namespace IMP - -#endif /* __IMP_QUADRATIC_BIPARTITE_NONBONDED_LIST_SCORE_STATE_H */ Index: kernel/include/IMP/score_states/BipartiteNonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/BipartiteNonbondedListScoreState.h (revision 524) +++ kernel/include/IMP/score_states/BipartiteNonbondedListScoreState.h (working copy) @@ -1,32 +1,92 @@ /** * \file BipartiteNonbondedListScoreState.h - * \brief Allow iteration through pairs of a set of atoms. + * \brief Bipartiteow iteration through pairs of a set of s. * - * Copyright 2007-8 Sali Lab. All rights reserved. + * Copyright 2007-8 Sali Lab. Bipartite rights reserved. */ #ifndef __IMP_BIPARTITE_NONBONDED_LIST_SCORE_STATE_H #define __IMP_BIPARTITE_NONBONDED_LIST_SCORE_STATE_H -#include "QuadraticBipartiteNonbondedListScoreState.h" +#include "NonbondedListScoreState.h" +#include "../internal/kernel_version_info.h" +#include "MaxChangeScoreState.h" + namespace IMP { -//! Maintain a nonbonded list between two disjoint sets. -/** - \note If no value for the radius key is specified, all radii are - considered to be zero. - \ingroup restraint +//! This class maintains a list of non-bonded pairs of spheres between two sets +/** To iterate through the list of pairs use the NonbondedListScoreState::begin, + NonbondedListScoreState::end functions. + + The radius key can be the default key. + + Changes in coordinates and radii are handled properly. + + \ingroup restraint */ -class BipartiteNonbondedListScoreState: - public QuadraticBipartiteNonbondedListScoreState { - typedef QuadraticBipartiteNonbondedListScoreState P; +class IMPDLLEXPORT BipartiteNonbondedListScoreState: + public NonbondedListScoreState +{ + typedef NonbondedListScoreState P; public: - BipartiteNonbondedListScoreState(FloatKey rk, + //! What algorithm to use to perform the computations + enum Algorithm { + //! Check all pairs of particles to see if they are close enough + QUADRATIC, +#ifdef IMP_USE_CGAL + //! Sweep space looking for intersecting bounding boxes. + BBOX, + DEFAULT= BBOX +#else + DEFAULT= QUADRATIC +#endif + }; +protected: + Algorithm a_; + internal::ObjectPointer<MaxChangeScoreState, true> mc0_, mc1_, mcr_; + + void process_sets(const Particles &p0, + const Particles &p1); + + //! \internal + void rebuild_nbl(); +public: + /**\param[in] cutoff The distance cutoff to use. + \param[in] radius The key to use to get the radius + \param[in] a Which algorithm to use. The default is the best. + */ + BipartiteNonbondedListScoreState(Float cutoff, FloatKey radius, + Algorithm a=DEFAULT); + /**\param[in] cutoff The distance cutoff to use. + \param[in] radius The key to use to get the radius + \param[in] ps0 The first set. + \param[in] ps1 The second set. + \param[in] a Which algorithm to use. The default is the best. + */ + BipartiteNonbondedListScoreState(Float cutoff, + FloatKey radius, const Particles &ps0, - const Particles &ps1): P(rk, ps0, ps1){} - BipartiteNonbondedListScoreState(FloatKey rk): P(rk){} + const Particles &ps1, + Algorithm a=DEFAULT); + + ~BipartiteNonbondedListScoreState(); + IMP_SCORE_STATE(internal::kernel_version_info); + + //! Add the particles to the first set + void add_particles_0(const Particles &ps); + //! Add the particles to the second set + void add_particles_1(const Particles &ps); + //! Remove all particles + void clear_particles(); + //! Replace the set of particles + void set_particles(const Particles &ps0, const Particles &ps1); + + /** If there is CGAL support, a more efficient algorithm BBOX can be used*/ + void set_algorithm(Algorithm a) { + a_=a; + } }; } // namespace IMP Index: kernel/include/IMP/score_states/BondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/BondedListScoreState.h (revision 524) +++ kernel/include/IMP/score_states/BondedListScoreState.h (working copy) @@ -25,9 +25,6 @@ BondedListScoreState(){} virtual ~BondedListScoreState(){} - //! Set the set of particles used - virtual void set_particles(const Particles &ps)=0; - //! Return true if the two particles are bonded virtual bool are_bonded(Particle *a, Particle *b) const =0; }; Index: kernel/include/IMP/score_states/QuadraticAllNonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/QuadraticAllNonbondedListScoreState.h (revision 524) +++ kernel/include/IMP/score_states/QuadraticAllNonbondedListScoreState.h (working copy) @@ -1,49 +0,0 @@ -/** - * \file QuadraticAllNonbondedListScoreState.h - * \brief Allow iteration through pairs of a set of spheres. - * - * Copyright 2007-8 Sali Lab. All rights reserved. - */ - -#ifndef __IMP_QUADRATIC_ALL_NONBONDED_LIST_SCORE_STATE_H -#define __IMP_QUADRATIC_ALL_NONBONDED_LIST_SCORE_STATE_H - -#include "NonbondedListScoreState.h" -#include "../internal/kernel_version_info.h" -#include "QuadraticNonbondedListScoreState.h" - -#include <vector> -#include <limits> - -namespace IMP -{ - -//! This class maintains a list of non-bonded pairs of spheres -/** \note QuadraticBipartiteNonbondedListScoreState is basically an - implementation detail for performance analysis and should not - be used by end users. - */ -class IMPDLLEXPORT QuadraticAllNonbondedListScoreState: - public QuadraticNonbondedListScoreState -{ - typedef QuadraticNonbondedListScoreState P; - - //! \internal - void rebuild_nbl(Float cut); - -public: - /** - \param[in] ps A list of particles to use. - \param[in] radius The key to use to get the radius - */ - QuadraticAllNonbondedListScoreState(FloatKey radius, - const Particles &ps= Particles()); - ~QuadraticAllNonbondedListScoreState(); - IMP_SCORE_STATE(internal::kernel_version_info) - - void set_particles(const Particles &ps); -}; - -} // namespace IMP - -#endif /* __IMP_QUADRATIC_ALL_NONBONDED_LIST_SCORE_STATE_H */ Index: kernel/include/IMP/score_states/ChainBondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/ChainBondedListScoreState.h (revision 0) +++ kernel/include/IMP/score_states/ChainBondedListScoreState.h (revision 0) @@ -0,0 +1,44 @@ +/** + * \file ChainBondedListScoreState.h + * \brief Allow iteration through pairs of a set of atoms. + * + * Copyright 2007-8 Sali Lab. All rights reserved. + */ + +#ifndef __IMP_CHAIN_BONDED_LIST_SCORE_STATE_H +#define __IMP_CHAIN_BONDED_LIST_SCORE_STATE_H + +#include "BondedListScoreState.h" +#include "../decorators/bond_decorators.h" +#include "../internal/Vector.h" + +namespace IMP +{ + +//! Exclude consecutive particles in a chain. +/** \ingroup bond + */ +class IMPDLLEXPORT ChainBondedListScoreState: public BondedListScoreState +{ + std::vector<internal::Vector<Particle*> > chains_; + IntKey cik_; + unsigned int next_index_; +public: + //! Find bonds amongst the following points. + /** \param [in] ps The set of particles to use. + */ + ChainBondedListScoreState(); + virtual ~ChainBondedListScoreState(){} + + void add_chain(const Particles &ps); + + void clear_chains(); + + virtual bool are_bonded(Particle *a, Particle *b) const; + + virtual void update(); +}; + +} // namespace IMP + +#endif /* __IMP_BOND_DECORATOR_LIST_SCORE_STATE_H */ Index: kernel/include/IMP/score_states/AllNonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/AllNonbondedListScoreState.h (revision 524) +++ kernel/include/IMP/score_states/AllNonbondedListScoreState.h (working copy) @@ -9,76 +9,93 @@ #define __IMP_ALL_NONBONDED_LIST_SCORE_STATE_H #include "NonbondedListScoreState.h" -#include "../internal/ParticleGrid.h" #include "../internal/kernel_version_info.h" -#include "MaxChangeScoreState.h" +#include "../internal/Vector.h" -#include <vector> -#include <limits> namespace IMP { -//! This class maintains a list of non-bonded pairs of particles -/** \note If no value for the radius key is specified, all radii are - considered to be zero. +class MaxChangeScoreState; +namespace internal +{ +class ParticleGrid; +} - \note The radius is currently assumed not to change. This could - be fixed later. +//! This class maintains a list of non-bonded pairs of spheres +/** To iterate through the list of pairs use the NonbondedListScoreState::begin, + NonbondedListScoreState::end functions. - \todo The structure is slightly dumb about rebuilding and will - rebuild the whole list of any of the grids become invalidated. - This could be improved as each piece is computed separately (so - they could be cached). + The radius key can be the default key. + Changes in coordinates and radii are handled properly unless grid is used, + then changes in radii may not be handled properly. + \ingroup restraint */ class IMPDLLEXPORT AllNonbondedListScoreState: public NonbondedListScoreState { typedef NonbondedListScoreState P; + public: + //! What algorithm to use to perform the computations + enum Algorithm { + //! Check all pairs of particles to see if they are close enough + QUADRATIC, + GRID, + //! Sweep space looking for intersecting bounding boxes. + BBOX, + DEFAULT + }; protected: - //! \internal - struct Bin - { - internal::ParticleGrid *grid; - Float rmax; - Bin(): grid(NULL), rmax(-1){} - Bin(const Bin &o): grid(o.grid), rmax(o.rmax){} - }; - std::vector<Bin> bins_; + Algorithm a_; + internal::ObjectPointer<MaxChangeScoreState, true> mc_, mcr_; + internal::Vector<internal::ParticleGrid*> bins_; + void check_nbl() const; + //! \internal - void rebuild_nbl(Float cut); + void rebuild_nbl(); - void repartition_points(const Particles &ps, std::vector<Bin> &out); - float side_from_r(float r) const; + // methods for grid + void grid_rebuild_nbl(); - void generate_nbl(const Bin &particle_bin, const Bin &grid_bin, float cut); + void grid_partition_points(); - void cleanup(std::vector<Bin> &bins); + void grid_generate_nbl(const internal::ParticleGrid *particle_bin, + const internal::ParticleGrid *grid_bin); + float grid_side_from_r(float r) const; public: - /** - \param[in] ps A list of particles to use. + /**\param[in] cutoff The distance cutoff to use. \param[in] radius The key to use to get the radius + \param[in] a Which algorithm to use. The default is the best. */ - AllNonbondedListScoreState(FloatKey radius, - const Particles &ps=Particles()); + AllNonbondedListScoreState(Float cutoff, + FloatKey radius, + Algorithm a= DEFAULT); + /**\param[in] cutoff The distance cutoff to use. + \param[in] radius The key to use to get the radius + \param[in] ps A list of particles to use. + \param[in] a Which algorithm to use. The default is the best. + */ + AllNonbondedListScoreState(Float cutoff, + FloatKey radius, + const Particles &ps, + Algorithm a= DEFAULT); ~AllNonbondedListScoreState(); - IMP_SCORE_STATE(internal::kernel_version_info) + IMP_SCORE_STATE(internal::kernel_version_info); + //! Add the particles and add them to the NBL + void add_particles(const Particles &ps); + //! Remove all particles + void clear_particles(); + //! Replace the set of particles void set_particles(const Particles &ps); - //! Add a few particles to the nonbonded list - /** Note that this invalidates the nonbonded list. - \todo We could just add newly created pairs to the nonbonded list. - */ - void add_particles(const Particles &ps); - - //! Return a list of all the particles used - Particles get_particles() const; + /** If there is CGAL support, a more efficient algorithm BBOX can be used*/ + void set_algorithm(Algorithm a); }; } // namespace IMP Index: kernel/src/score_states/BondDecoratorListScoreState.cpp =================================================================== --- kernel/src/score_states/BondDecoratorListScoreState.cpp (revision 524) +++ kernel/src/score_states/BondDecoratorListScoreState.cpp (working copy) @@ -24,7 +24,7 @@ bonds_.clear(); for (unsigned int i=0; i< ps_.size(); ++i) { if (!ps_[i]->get_is_active()) continue; - BondedDecorator di= BondedDecorator::cast(ps_[i]); + BondedDecorator di(ps_[i]); ParticleIndex pi= ps_[i]->get_index(); for (unsigned int j=0; j< di.get_number_of_bonds(); ++j) { BondedDecorator dj= di.get_bonded(j); Index: kernel/src/score_states/NonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/NonbondedListScoreState.cpp (revision 524) +++ kernel/src/score_states/NonbondedListScoreState.cpp (working copy) @@ -15,22 +15,27 @@ namespace IMP { -NonbondedListScoreState::NonbondedListScoreState(FloatKey rk): rk_(rk) +NonbondedListScoreState +::NonbondedListScoreState(Float cut, + FloatKey rk): rk_(rk), + cutoff_(cut), + nbl_is_valid_(false) { - last_cutoff_=-1; + slack_=10; + num_steps_=0; } - - - - -void NonbondedListScoreState::propagate_particles(const Particles&ps) -{ - clear_nbl(); - for (BondedListScoreStateIterator bli= bonded_lists_begin(); - bli != bonded_lists_end(); ++bli) { - (*bli)->set_particles(ps); +Particles NonbondedListScoreState +::particles_with_radius(const Particles &in) const { + Particles ret; + if (rk_== FloatKey()) return ret; + ret.reserve(in.size()); + for (unsigned int i=0; i< in.size(); ++i) { + if (in[i]->has_attribute(rk_)) { + ret.push_back(in[i]); + } } + return ret; } namespace internal @@ -45,7 +50,7 @@ } // namespace internal -void NonbondedListScoreState::update() +bool NonbondedListScoreState::update(Float mc) { IMP_LOG(VERBOSE, "Updating non-bonded list" << std::endl); for (BondedListScoreStateIterator bli= bonded_lists_begin(); @@ -53,9 +58,21 @@ (*bli)->update(); } - // if the list is not deleted, we need to scan for inactive particles - nbl_.erase(std::remove_if(nbl_.begin(), nbl_.end(), internal::HasInactive()), - nbl_.end()); + if (mc > slack_ || nbl_.size() ==0) { + num_steps_=0; + set_nbl_is_valid(false); + slack_= next_slack_; + rebuild_nbl(); + return true; + } else { + // if the list is not deleted, we need to scan for inactive particles + nbl_.erase(std::remove_if(nbl_.begin(), nbl_.end(), + internal::HasInactive()), + nbl_.end()); + ++num_steps_; + next_slack_= mc/num_steps_ * 100.0; + return false; + } } void NonbondedListScoreState::show(std::ostream &out) const Index: kernel/src/score_states/QuadraticBipartiteNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/QuadraticBipartiteNonbondedListScoreState.cpp (revision 524) +++ kernel/src/score_states/QuadraticBipartiteNonbondedListScoreState.cpp (working copy) @@ -1,66 +0,0 @@ -/** - * \file QuadraticBipartiteNonbondedListScoreState.cpp - * \brief Allow iteration through pairs of a set of atoms. - * - * Copyright 2007-8 Sali Lab. All rights reserved. - */ - -#include "IMP/score_states/QuadraticBipartiteNonbondedListScoreState.h" -#include "IMP/decorators/XYZDecorator.h" -#include "IMP/score_states/MaxChangeScoreState.h" - -namespace IMP -{ - -QuadraticBipartiteNonbondedListScoreState -::QuadraticBipartiteNonbondedListScoreState(FloatKey rk, - const Particles &ps0, - const Particles &ps1): P(rk) -{ - set_particles(ps0, ps1); -} - - -QuadraticBipartiteNonbondedListScoreState -::QuadraticBipartiteNonbondedListScoreState(FloatKey rk): P(rk) -{ -} - -void QuadraticBipartiteNonbondedListScoreState::rebuild_nbl(float cut) -{ - IMP_LOG(TERSE, "Rebuilding QBNBL with lists of size " << na_ - << " and " << P::get_particles().size() -na_ - << " and cutoff " << cut << std::endl); - for (unsigned int i=0; i< na_; ++i) { - for (unsigned int j=na_; j < P::get_particles().size(); ++j) { - P::handle_nbl_pair(P::get_particles()[i], - P::get_particles()[j], - cut); - } - } - IMP_LOG(TERSE, "NBL has " << P::nbl_.size() << " pairs" << std::endl); -} - -void QuadraticBipartiteNonbondedListScoreState -::set_particles(const Particles &ps0, - const Particles &ps1) -{ - IMP_LOG(TERSE, "Setting QBNBLSS particles " << ps0.size() - << " and " << ps1.size() << std::endl); - Particles all(ps0); - all.insert(all.end(), ps1.begin(), ps1.end()); - P::set_particles(all); - na_= ps0.size(); -} - -void QuadraticBipartiteNonbondedListScoreState::update() -{ - P::update(); -} - -void QuadraticBipartiteNonbondedListScoreState::show(std::ostream &out) const -{ - out << "QuadraticBipartiteNonbondedList" << std::endl; -} - -} // namespace IMP Index: kernel/src/score_states/BipartiteNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/BipartiteNonbondedListScoreState.cpp (revision 0) +++ kernel/src/score_states/BipartiteNonbondedListScoreState.cpp (revision 0) @@ -0,0 +1,128 @@ +/** + * \file BipartiteNonbondedListScoreState.cpp + * \brief Bipartiteow iteration through pairs of a set of s. + * + * Copyright 2007-8 Sali Lab. Bipartite rights reserved. + */ + +#include "IMP/score_states/BipartiteNonbondedListScoreState.h" +#include "IMP/decorators/XYZDecorator.h" +#include "IMP/internal/utility.h" +#include "IMP/internal/bbox_nbl_helpers.h" + +#include <vector> + +namespace IMP +{ + +BipartiteNonbondedListScoreState +::BipartiteNonbondedListScoreState(Float cut, + FloatKey rk, + Algorithm a): + P(cut, rk), a_(a) +{ + FloatKeys ks; + ks.push_back(rk); + mc0_= new MaxChangeScoreState(XYZDecorator::get_xyz_keys()); + mc1_= new MaxChangeScoreState(XYZDecorator::get_xyz_keys()); + mcr_= new MaxChangeScoreState(ks); +} + +BipartiteNonbondedListScoreState +::BipartiteNonbondedListScoreState(Float cut, FloatKey rk, + const Particles &ps0, + const Particles &ps1, + Algorithm a): + P(cut, rk), a_(a) +{ + FloatKeys ks; + ks.push_back(rk); + mc0_= new MaxChangeScoreState(XYZDecorator::get_xyz_keys()); + mc1_= new MaxChangeScoreState(XYZDecorator::get_xyz_keys()); + mcr_= new MaxChangeScoreState(ks); + set_particles(ps0, ps1); +} + + +BipartiteNonbondedListScoreState::~BipartiteNonbondedListScoreState(){} + +void BipartiteNonbondedListScoreState::update() { + mc0_->update(); + mc0_->update(); + mcr_->update(); + Float mc= std::max(mc0_->get_max(), mc1_->get_max()); + if (P::update(mc+ mcr_->get_max())) { + mc0_->reset(); + mc1_->reset(); + mcr_->reset(); + } +} + +void BipartiteNonbondedListScoreState::process_sets(const Particles &p0, + const Particles &p1) { + if (a_== QUADRATIC) { + for (unsigned int j=0; j< p0.size(); ++j) { + for (unsigned int i=0; i< p1.size(); ++i) { + P::add_if_box_overlap(p1[i], p0[j]); + } + } + } +#ifdef IMP_USE_CGAL + else if (a_== BBOX) { + internal::bipartite_bbox_scan(p0, p1, P::get_radius_key(), + P::get_slack(), P::get_cutoff(), + internal::NBLAddPairIfNonbonded(this)); + + } +#endif + } + +void BipartiteNonbondedListScoreState::rebuild_nbl() { + IMP_LOG(TERSE, "Rebuilding QNBL with cutoff " + << P::get_cutoff() << std::endl); + process_sets(mc0_->get_particles(), mc1_->get_particles()); + P::set_nbl_is_valid(true); + IMP_LOG(TERSE, "NBL has " << P::number_of_nonbonded() + << " pairs" << std::endl); +} + +void BipartiteNonbondedListScoreState::set_particles(const Particles &ps0, + const Particles &ps1) { + mc0_->clear_particles(); + mc0_->add_particles(ps0); + mc1_->clear_particles(); + mc1_->add_particles(ps1); + mcr_->clear_particles(); + mcr_->add_particles(P::particles_with_radius(ps0)); + mcr_->add_particles(P::particles_with_radius(ps1)); + P::set_nbl_is_valid(false); +} + + +void BipartiteNonbondedListScoreState::add_particles_0(const Particles &ps) { + if (P::get_nbl_is_valid()) process_sets(ps, mc1_->get_particles()); + mc0_->add_particles(ps); + mcr_->add_particles(P::particles_with_radius(ps)); +} + +void BipartiteNonbondedListScoreState::add_particles_1(const Particles &ps) { + if (P::get_nbl_is_valid()) process_sets(ps, mc0_->get_particles()); + mc1_->add_particles(ps); + mcr_->add_particles(P::particles_with_radius(ps)); +} + +void BipartiteNonbondedListScoreState::clear_particles() { + mc0_->clear_particles(); + mc1_->clear_particles(); + mcr_->clear_particles(); + P::set_nbl_is_valid(false); + P::set_nbl_is_valid(true); +} + + +void BipartiteNonbondedListScoreState::show(std::ostream &out) const +{ + out << "BipartiteNonbondedListScoreState" << std::endl; +} + +} // namespace IMP Index: kernel/src/score_states/SConscript =================================================================== --- kernel/src/score_states/SConscript (revision 524) +++ kernel/src/score_states/SConscript (working copy) @@ -1,11 +1,10 @@ Import('env') - -files = ['NonbondedListScoreState.cpp', - 'MaxChangeScoreState.cpp', 'BondDecoratorListScoreState.cpp', - 'AllNonbondedListScoreState.cpp', - 'QuadraticBipartiteNonbondedListScoreState.cpp', - 'QuadraticAllNonbondedListScoreState.cpp', - 'QuadraticNonbondedListScoreState.cpp', 'GravityCenterScoreState.cpp'] - -files = [File(x) for x in files] +files=[] +files.append(File( 'AllNonbondedListScoreState.cpp' )) +files.append(File( 'BipartiteNonbondedListScoreState.cpp' )) +files.append(File( 'BondDecoratorListScoreState.cpp' )) +files.append(File( 'ChainBondedListScoreState.cpp' )) +files.append(File( 'GravityCenterScoreState.cpp' )) +files.append(File( 'MaxChangeScoreState.cpp' )) +files.append(File( 'NonbondedListScoreState.cpp' )) Return('files') Index: kernel/src/score_states/QuadraticAllNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/QuadraticAllNonbondedListScoreState.cpp (revision 524) +++ kernel/src/score_states/QuadraticAllNonbondedListScoreState.cpp (working copy) @@ -1,58 +0,0 @@ -/** - * \file QuadraticAllNonbondedListScoreState.cpp - * \brief Allow iteration through pairs of a set of s. - * - * Copyright 2007-8 Sali Lab. All rights reserved. - */ - -#include "IMP/score_states/QuadraticAllNonbondedListScoreState.h" -#include "IMP/decorators/XYZDecorator.h" - -#include <algorithm> - -namespace IMP -{ - - -QuadraticAllNonbondedListScoreState -::QuadraticAllNonbondedListScoreState(FloatKey radius, - const Particles &ps): P(radius) -{ - set_particles(ps); -} - -QuadraticAllNonbondedListScoreState::~QuadraticAllNonbondedListScoreState() -{ -} - -void QuadraticAllNonbondedListScoreState::set_particles(const Particles &ps) -{ - P::set_particles(ps); -} - - -void QuadraticAllNonbondedListScoreState::update() -{ - P::update(); -} - - -void QuadraticAllNonbondedListScoreState::rebuild_nbl(Float cut) -{ - IMP_LOG(TERSE, "Rebuilding QNBL with cutoff " << cut << std::endl); - const Particles &moving= P::get_particles(); - for (unsigned int j=0; j< moving.size(); ++j) { - for (unsigned int i=0; i< j; ++i) { - P::handle_nbl_pair(moving[i], moving[j], cut); - } - } - IMP_LOG(TERSE, "NBL has " << P::nbl_.size() << " pairs" << std::endl); -} - - -void QuadraticAllNonbondedListScoreState::show(std::ostream &out) const -{ - out << "QuadraticAllNonbondedListScoreState" << std::endl; -} - -} // namespace IMP Index: kernel/src/score_states/ChainBondedListScoreState.cpp =================================================================== --- kernel/src/score_states/ChainBondedListScoreState.cpp (revision 0) +++ kernel/src/score_states/ChainBondedListScoreState.cpp (revision 0) @@ -0,0 +1,64 @@ +/** + * \file ChainBondedListScoreState.cpp + * \brief Allow iteration through pairs of a set of atoms. + * + * Copyright 2007-8 Sali Lab. All rights reserved. + */ + +#include "IMP/score_states/ChainBondedListScoreState.h" + +#include <algorithm> +#include <sstream> + +namespace IMP +{ + +ChainBondedListScoreState::ChainBondedListScoreState() +{ + std::ostringstream oss; + oss << "ChainBLSS " << this; + cik_= IntKey(oss.str().c_str()); + next_index_=0; +} + +void ChainBondedListScoreState::update() +{ + IMP_LOG(TERSE, "Updating ChainBondedList assumed to be static" << std::endl); +} + +void ChainBondedListScoreState::clear_chains() +{ + for (unsigned int i=0; i< chains_.size(); ++i) { + for (unsigned int j=0; j< chains_[i].size(); ++j) { + chains_[i][j]->set_value(cik_, -1); + } + } + chains_.clear(); + next_index_=0; +} + +void ChainBondedListScoreState::add_chain(const Particles &ps) +{ + chains_.push_back(internal::Vector<Particle*>(ps.begin(), ps.end())); + for (unsigned int i=0; i< chains_.back().size(); ++i) { + Particle *p= chains_.back()[i]; + if (p->has_attribute(cik_)) { + p->set_value(cik_, next_index_); + } else { + p->add_attribute(cik_, next_index_); + } + ++next_index_; + } + ++next_index_; +} + + +bool ChainBondedListScoreState::are_bonded(Particle *a, Particle *b) const +{ + if (!a->has_attribute(cik_) || !b->has_attribute(cik_)) return false; + int ia= a->get_value(cik_); + int ib= b->get_value(cik_); + return std::abs(ia-ib) ==1; +} + +} // namespace IMP Index: kernel/src/score_states/MaxChangeScoreState.cpp =================================================================== --- kernel/src/score_states/MaxChangeScoreState.cpp (revision 524) +++ kernel/src/score_states/MaxChangeScoreState.cpp (working copy) @@ -31,12 +31,15 @@ IMP_LIST_IMPL(MaxChangeScoreState, Particle, particle, Particle*, {for (unsigned int i=0; i< keys_.size(); ++i) { IMP_check(obj->has_attribute(keys_[i]), - "Particle missing needed attribute", + "Particle missing needed attribute " << keys_[i] + << obj, ValueException("Particle missing attribute")); }; for (unsigned int i=0; i< origkeys_.size(); ++i) { - obj->add_attribute(origkeys_[i], - obj->get_value(keys_[i]), false); + if (!obj->has_attribute(origkeys_[i])) { + obj->add_attribute(origkeys_[i], + obj->get_value(keys_[i]), false); + } } }, {reset();}); Index: kernel/src/score_states/AllNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/AllNonbondedListScoreState.cpp (revision 524) +++ kernel/src/score_states/AllNonbondedListScoreState.cpp (working copy) @@ -7,101 +7,170 @@ #include "IMP/score_states/AllNonbondedListScoreState.h" #include "IMP/decorators/XYZDecorator.h" +#include "IMP/internal/utility.h" +#include "IMP/internal/bbox_nbl_helpers.h" -#include <algorithm> +#include <vector> namespace IMP { -static unsigned int min_grid_size=20; +// Turn the default into an actuall algorithm and work around missing algorithms +static AllNonbondedListScoreState::Algorithm +translate_algorithm(AllNonbondedListScoreState::Algorithm a){ +#ifdef IMP_USE_CGAL + switch (a) { + case AllNonbondedListScoreState::DEFAULT: + return AllNonbondedListScoreState::BBOX; + default: + return a; + } +#else + switch (a) { + case AllNonbondedListScoreState::BBOX: + IMP_WARN("AllNonbondedListScoreState::BBOX requires CGAL support. " + << "GRID used instead." << std::endl); + case AllNonbondedListScoreState::DEFAULT: + return AllNonbondedListScoreState::GRID; + default: + return a; + } +#endif +} -AllNonbondedListScoreState -::AllNonbondedListScoreState(FloatKey radius, const Particles &ps): - P(radius) +AllNonbondedListScoreState::AllNonbondedListScoreState(Float cut, + FloatKey rk, + const Particles &ps, + Algorithm a): + P(cut, rk), a_(translate_algorithm(a)) { - set_particles(ps); + FloatKeys ks; + ks.push_back(rk); + mc_= new MaxChangeScoreState(XYZDecorator::get_xyz_keys()); + mcr_= new MaxChangeScoreState(ks); + add_particles(ps); } - -AllNonbondedListScoreState::~AllNonbondedListScoreState() +AllNonbondedListScoreState::AllNonbondedListScoreState(Float cut, + FloatKey rk, + Algorithm a): + P(cut, rk), a_(translate_algorithm(a)) { - cleanup(bins_); + FloatKeys ks; + ks.push_back(rk); + mc_= new MaxChangeScoreState(XYZDecorator::get_xyz_keys()); + mcr_= new MaxChangeScoreState(ks); } -float AllNonbondedListScoreState::side_from_r(float r) const -{ - if (r==0) return 1; - else return r*1.6; +AllNonbondedListScoreState::~AllNonbondedListScoreState(){ } -Particles AllNonbondedListScoreState::get_particles() const -{ - Particles ret; - for (unsigned int i=0; i< bins_.size(); ++i) { - ret.insert(ret.end(), bins_[i].grid->get_particles().begin(), - bins_[i].grid->get_particles().end()); +void AllNonbondedListScoreState::update() { + mc_->update(); + mcr_->update(); + Float mc=mc_->get_max()+ mcr_->get_max(); + if (P::update(mc)) { + mc_->reset(); } - return ret; } -void AllNonbondedListScoreState::set_particles(const Particles &ps) -{ - NonbondedListScoreState::clear_nbl(); - cleanup(bins_); - repartition_points(ps, bins_); -} -void AllNonbondedListScoreState::add_particles(const Particles &ps) -{ - if (bins_.empty()) set_particles(ps); - else { -#ifndef NDEBUG - for (unsigned int i=0; i< ps.size(); ++i) { - for (unsigned int j=0; j< bins_.size(); ++j) { - for (unsigned int k=0; k< bins_[j].grid->get_particles().size(); ++k) { - IMP_assert(ps[i] != bins_[j].grid->get_particles()[k], - "Can't add a particle that is already there"); - - IMP_assert(ps[i]->get_index() - != bins_[j].grid->get_particles()[k]->get_index(), - "Same particle index, different particles."); - } +void AllNonbondedListScoreState::rebuild_nbl() { + IMP_LOG(TERSE, "Rebuilding QNBL with cutoff " + << P::get_cutoff() << std::endl); + if (a_== QUADRATIC) { + const Particles &moving= mc_->get_particles(); + for (unsigned int j=0; j< moving.size(); ++j) { + for (unsigned int i=0; i< j; ++i) { + P::add_if_box_overlap(moving[i], moving[j]); } } + } else if (a_ == GRID) { + grid_rebuild_nbl(); + } else if (a_== BBOX) { + internal::bbox_scan(mc_->get_particles(), P::get_radius_key(), + P::get_slack(), P::get_cutoff(), + internal::NBLAddPairIfNonbonded(this)); + + } else { + IMP_failure("Bad algorithm in AllNBL::rebuild", ErrorException()); + } + set_nbl_is_valid(true); +#ifndef NDEBUG + check_nbl(); #endif + IMP_LOG(TERSE, "NBL has " << P::number_of_nonbonded() + << " pairs" << std::endl); +} - for (unsigned int i=0; i< ps.size(); ++i) { - float r= P::get_radius(ps[i]); - bool found=false; - for (unsigned int j=0; j< bins_.size(); ++j) { - if (bins_[j].rmax >= r) { - found=true; - IMP_LOG(VERBOSE, "Adding particle " - << ps[i]->get_index() << " to bin " << j << std::endl); - bins_[j].grid->add_particle(ps[i]); - break; +void AllNonbondedListScoreState::set_particles(const Particles &ps) { + mc_->clear_particles(); + mc_->add_particles(ps); + mcr_->clear_particles(); + mcr_->add_particles(particles_with_radius(ps)); + P::set_nbl_is_valid(false); +} + + +void AllNonbondedListScoreState::add_particles(const Particles &ps) { + if (P::get_nbl_is_valid()) { + if (a_== QUADRATIC || a_ == GRID) { + const Particles &moving= mc_->get_particles(); + for (unsigned int j=0; j< moving.size(); ++j) { + for (unsigned int i=0; i< ps.size(); ++i) { + P::add_if_box_overlap(ps[i], moving[j]); } } - if (!found) { - bins_.back().rmax=r; - bins_.back().grid->add_particle(ps[i]); - } + } else if (a_== BBOX) { + internal::bipartite_bbox_scan(mc_->get_particles(), ps, + P::get_radius_key(), + P::get_slack(), P::get_cutoff(), + internal::NBLAddPairIfNonbonded(this)); } } - IMP_LOG(TERSE, "Destroying nbl in list due to additions"<< std::endl); - NonbondedListScoreState::clear_nbl(); + mcr_->add_particles(particles_with_radius(ps)); + mc_->add_particles(ps); } +void AllNonbondedListScoreState::clear_particles() { + mc_->clear_particles(); + mcr_->clear_particles(); + P::set_nbl_is_valid(false); + P::set_nbl_is_valid(true); +} -void AllNonbondedListScoreState::repartition_points(const Particles &ps, - std::vector<Bin> &out) +void AllNonbondedListScoreState::show(std::ostream &out) const { - cleanup(out); - if (ps.empty()) return; + out << "AllNonbondedListScoreState" << std::endl; +} + + +void AllNonbondedListScoreState::set_algorithm(Algorithm a) { + a_=translate_algorithm(a); + if (a_!= GRID) { + bins_.clear(); + } +} + + +// methods for grid + +static unsigned int min_grid_size=20; + +float AllNonbondedListScoreState::grid_side_from_r(float r) const +{ + if (r==0) return 1; + else return r*1.6; +} + +void AllNonbondedListScoreState::grid_partition_points() +{ + bins_.clear(); + if (mc_->get_particles().empty()) return; float minr=std::numeric_limits<float>::max(), maxr=0; - for (unsigned int i=0; i< ps.size(); ++i) { - float r= P::get_radius(ps[i]); + for (unsigned int i=0; i< mc_->get_particles().size(); ++i) { + float r= P::get_radius(mc_->get_particles()[i]); if ( r > maxr) maxr=r; if ( r > 0 && r < minr) minr=r; } @@ -116,13 +185,13 @@ cuts.push_back(2*maxr); std::vector<Particles> ops(cuts.size()); - for (unsigned int i=0; i< ps.size(); ++i) { - float r= P::get_radius(ps[i]); + for (unsigned int i=0; i< mc_->get_particles().size(); ++i) { + float r= P::get_radius(mc_->get_particles()[i]); bool found=false; for (unsigned int j=0; ; ++j) { IMP_assert(j< cuts.size(), "Internal error in ASNBLSS"); if (cuts[j] >= r) { - ops[j].push_back(ps[i]); + ops[j].push_back(mc_->get_particles()[i]); found=true; break; } @@ -138,156 +207,110 @@ } for (unsigned int i=0; i< cuts.size(); ++i) { if (ops[i].empty()) continue; - out.push_back(Bin()); float rmax=0; for (unsigned int j=0; j< ops[i].size(); ++j) { rmax= std::max(rmax, P::get_radius(ops[i][j])); } - out.back().rmax= rmax; - internal::ParticleGrid *pg - = new internal::ParticleGrid(side_from_r(out.back().rmax)); - out.back().grid= pg; - out.back().grid->add_particles(ops[i]); + bins_.push_back(new internal::ParticleGrid(grid_side_from_r(rmax), ops[i])); } - IMP_LOG(VERBOSE, "Created " << out.size() << " grids" << std::endl); - for (unsigned int i=0; i< out.size(); ++i) { - IMP_LOG(VERBOSE, out[i].rmax - << ": " << *out[i].grid << std::endl); + IMP_LOG(VERBOSE, "Created " << bins_.size() << " grids" << std::endl); + for (unsigned int i=0; i< bins_.size(); ++i) { + IMP_LOG(VERBOSE, *bins_[i] << std::endl); } - -#ifndef NDEBUG - Particles all; - for (unsigned int i=0; i< out.size(); ++i) { - all.insert(all.end(), out[i].grid->get_particles().begin(), - out[i].grid->get_particles().end()); - } - std::sort(all.begin(), all.end()); - Particles::iterator it = std::unique(all.begin(), all.end()); - IMP_assert(it == all.end(), "Duplicate particles " << all.size()); - IMP_assert(all.size() == ps.size(), "Wrong number of particles at end " - << all.size() << " " << ps.size()); -#endif } -void AllNonbondedListScoreState::cleanup(std::vector<Bin> &bins) -{ - for (unsigned int i=0; i< bins.size(); ++i) { - delete bins[i].grid; - } - bins.clear(); -} -void AllNonbondedListScoreState::update() +void AllNonbondedListScoreState +::grid_generate_nbl(const internal::ParticleGrid *particle_bin, + const internal::ParticleGrid *grid_bin) { - IMP_LOG(TERSE, "Updating nonbonded list"<< std::endl); - NonbondedListScoreState::update(); - bool bad=false; - for (unsigned int i=0; i< bins_.size(); ++i) { - bad = bins_[i].grid->update() || bad; - IMP_LOG(VERBOSE, bins_[i].rmax << std::endl << *bins_[i].grid << std::endl); - } - if (bad) { - IMP_LOG(TERSE, "Destroying nbl in list"<< std::endl); - NonbondedListScoreState::clear_nbl(); - } -} - -void AllNonbondedListScoreState::generate_nbl(const Bin &particle_bin, - const Bin &grid_bin, - float cut) -{ - IMP_CHECK_OBJECT(particle_bin.grid); - IMP_CHECK_OBJECT(grid_bin.grid); - IMP_LOG(VERBOSE, "Generate nbl for " << particle_bin.rmax - << " and " << grid_bin.rmax << std::endl); - for (unsigned int k=0; k< particle_bin.grid->get_particles().size(); ++k) { - Particle *p= particle_bin.grid->get_particles()[k]; - NonbondedListScoreState::AddToNBL f(this, p); + IMP_CHECK_OBJECT(particle_bin); + IMP_CHECK_OBJECT(grid_bin); + IMP_LOG(VERBOSE, "Generate nbl for pair " << std::endl); + for (internal::ParticleGrid::ParticleVoxelIterator + it= particle_bin->particle_voxels_begin(); + it != particle_bin->particle_voxels_end(); ++it) { + Particle *p= it->first; + internal::NBLAddIfNonbonded f(this, p); XYZDecorator d= XYZDecorator::cast(p); internal::ParticleGrid::VirtualIndex index - = grid_bin.grid->get_virtual_index(Vector3D(d.get_x(), - d.get_y(), - d.get_z())); + = grid_bin->get_virtual_index(Vector3D(d.get_x(), + d.get_y(), + d.get_z())); IMP_LOG(VERBOSE, "Searching for " << p->get_index() - << " from " << index - << " of bin " << grid_bin.rmax << std::endl); - grid_bin.grid->apply_to_nearby(f, index, - cut+2*particle_bin.rmax + grid_bin.rmax, - false); + << " from " << index << std::endl); + grid_bin->apply_to_nearby(f, index, + P::get_cutoff() + 2*P::get_slack(), + false); } } -void AllNonbondedListScoreState::rebuild_nbl(Float cut) +void AllNonbondedListScoreState::grid_rebuild_nbl() { - IMP_LOG(TERSE, "Rebuilding NBL with " << bins_.size() << " bins" - << " and cutoff " << cut << std::endl); + IMP_LOG(TERSE, "Rebuilding NBL with Grid and cutoff " + << P::get_cutoff() << std::endl ); + grid_partition_points(); + for (unsigned int i=0; i< bins_.size(); ++i) { for (unsigned int j=i+1; j< bins_.size(); ++j) { - generate_nbl(bins_[i], bins_[j], cut); + grid_generate_nbl(bins_[i], bins_[j]); } internal::ParticleGrid::Index last_index; for (internal::ParticleGrid::ParticleVoxelIterator it - = bins_[i].grid->particle_voxels_begin(); - it != bins_[i].grid->particle_voxels_end(); ++it) { + = bins_[i]->particle_voxels_begin(); + it != bins_[i]->particle_voxels_end(); ++it) { IMP_LOG(VERBOSE, "Searching with particle " << it->first->get_index() << std::endl); - NonbondedListScoreState::AddToNBL f(this, it->first); - bins_[i].grid->apply_to_nearby(f, it->second, - cut+3*bins_[i].rmax, + internal::NBLAddIfNonbonded f(this, it->first); + bins_[i]->apply_to_nearby(f, it->second, + P::get_cutoff()+2*P::get_slack(), true); + if (it->second != last_index) { + internal::NBLAddPairIfNonbonded fp(this); IMP_LOG(VERBOSE, "Searching in " << it->second << std::endl); - bins_[i].grid->apply_to_cell_pairs(f, it->second); + bins_[i]->apply_to_cell_pairs(fp, it->second); } last_index=it->second; } } - IMP_LOG(TERSE, "NBL has " << P::nbl_.size() << " pairs" << std::endl); + bins_.clear(); + IMP_LOG(TERSE, "NBL has " << P::number_of_nonbonded() + << " pairs" << std::endl); +} -#ifndef NDEBUG - { - Particles ps; - for (unsigned int i=0; i< bins_.size(); ++i) { - if (bins_[i].grid) { - ps.insert(ps.end(), bins_[i].grid->get_particles().begin(), - bins_[i].grid->get_particles().end()); - } - } - for (unsigned int i=0; i< ps.size(); ++i) { - XYZDecorator di= XYZDecorator::cast(ps[i]); - for (unsigned int j=0; j< i; ++j) { - XYZDecorator dj= XYZDecorator::cast(ps[j]); - IMP_assert(ps[i] != ps[j], "Duplicate particles in grid"); - if (distance(di, dj) - P::get_radius(ps[i]) - P::get_radius(ps[j]) - <= cut && !are_bonded(ps[i], ps[j])) { - bool found=false; - for (NonbondedIterator nit= nbl_.begin(); - nit != nbl_.end(); ++nit) { - if (nit->first == ps[i] && nit->second == ps[j] - || nit->first == ps[j] && nit->second == ps[i]) { - IMP_assert(!found, "Entry is in list twice"); - found=true; - } + +// debugging +void AllNonbondedListScoreState::check_nbl() const +{ + const Particles &ps= mc_->get_particles(); + for (unsigned int i=0; i< ps.size(); ++i) { + XYZDecorator di= XYZDecorator::cast(ps[i]); + for (unsigned int j=0; j< i; ++j) { + XYZDecorator dj= XYZDecorator::cast(ps[j]); + IMP_assert(ps[i] != ps[j], "Duplicate particles in grid"); + if (distance(di, dj) - P::get_radius(ps[i]) - P::get_radius(ps[j]) + <= P::get_cutoff() && !are_bonded(ps[i], ps[j])) { + bool found=false; + for (NonbondedIterator nit= P::nonbonded_begin(); + nit != P::nonbonded_end(); ++nit) { + if (nit->first == ps[i] && nit->second == ps[j] + || nit->first == ps[j] && nit->second == ps[i]) { + IMP_assert(!found, "Entry is in list twice"); + found=true; } - IMP_assert(found, "Nonbonded list is missing " - << ps[i]->get_index() << " " << di - << " and " << ps[j]->get_index() << " " - << dj << std::endl); } + IMP_assert(found, "Nonbonded list is missing " + << ps[i]->get_index() << " " << di + << " and " << ps[j]->get_index() << " " + << dj << std::endl); } } } -#endif } - -void AllNonbondedListScoreState::show(std::ostream &out) const -{ - out << "AllNonbondedListScoreState" << std::endl; -} - } // namespace IMP Index: kernel/src/score_states/QuadraticNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/QuadraticNonbondedListScoreState.cpp (revision 524) +++ kernel/src/score_states/QuadraticNonbondedListScoreState.cpp (working copy) @@ -1,60 +0,0 @@ -/** - * \file QuadraticNonbondedListScoreState.cpp - * \brief Allow iteration through pairs of a set of s. - * - * Copyright 2007-8 Sali Lab. All rights reserved. - */ - -#include "IMP/score_states/QuadraticNonbondedListScoreState.h" -#include "IMP/decorators/XYZDecorator.h" - -#include <algorithm> - -namespace IMP -{ - -QuadraticNonbondedListScoreState -::QuadraticNonbondedListScoreState(FloatKey radius): - P(radius), slack_(.1) -{ - mc_= new MaxChangeScoreState(XYZDecorator::get_xyz_keys()); -} - -QuadraticNonbondedListScoreState::~QuadraticNonbondedListScoreState() -{ -} - - -void QuadraticNonbondedListScoreState::update() -{ - // placeholder to do tuning of slack - NonbondedListScoreState::update(); - if (mc_->get_max() > slack_) { - NonbondedListScoreState::clear_nbl(); - mc_->reset(); - } -} - -void QuadraticNonbondedListScoreState -::handle_nbl_pair( Particle *a, Particle *b, - float d) -{ - XYZDecorator da= XYZDecorator::cast(a); - XYZDecorator db= XYZDecorator::cast(b); - float ra= P::get_radius(a); - float rb= P::get_radius(b); - for (unsigned int i=0; i< 3; ++i) { - float delta=std::abs(da.get_coordinate(i) - db.get_coordinate(i)); - if (delta -ra -rb > d-slack_) { - IMP_LOG(VERBOSE, "Pair " << a->get_index() - << " and " << b->get_index() << " rejected on coordinate " - << i << std::endl); - return ; - } - } - IMP_LOG(VERBOSE, "Adding pair " << a->get_index() - << " and " << b->get_index() << std::endl); - P::add_to_nbl(a, b); -} - -} // namespace IMP Index: kernel/include/IMP/internal/ParticleGrid.h =================================================================== --- kernel/include/IMP/internal/ParticleGrid.h (revision 524) +++ kernel/include/IMP/internal/ParticleGrid.h (working copy) @@ -25,26 +25,18 @@ // don't need ref counting since mc_ has the same set of points typedef internal::Grid3D<Particles> Grid; Grid grid_; - internal::ObjectPointer<MaxChangeScoreState, true> mc_; Float target_voxel_side_; - bool grid_valid_; - void build_grid(); + void build_grid(const Particles &ps); void audit_particles(const Particles &ps) const; void add_particle_to_grid(Particle *p); public: ParticleGrid(); //! suggested grid edge size. - ParticleGrid(Float sz); + ParticleGrid(Float sz, const Particles &ps); Float get_voxel_size() const {return target_voxel_side_;} - void add_particles(const Particles &ps); - void add_particle(Particle *p); - void clear_particles(); - const Particles& get_particles() const {return mc_->get_particles();} - bool update(); - void show(std::ostream &out) const; typedef Grid::VirtualIndex VirtualIndex; Index: kernel/src/internal/ParticleGrid.cpp =================================================================== --- kernel/src/internal/ParticleGrid.cpp (revision 524) +++ kernel/src/internal/ParticleGrid.cpp (working copy) @@ -18,19 +18,20 @@ static const int target_cell_occupancy=10; -ParticleGrid::ParticleGrid(float tvs): target_voxel_side_(tvs), - grid_valid_(false) + ParticleGrid::ParticleGrid(float tvs, + const Particles &ps): target_voxel_side_(tvs) { IMP_assert(tvs >0, "Target voxel edge size must be positive"); - mc_= new MaxChangeScoreState(XYZDecorator::get_xyz_keys()); + build_grid(ps); } -ParticleGrid::ParticleGrid(): target_voxel_side_(0), grid_valid_(0) +ParticleGrid::ParticleGrid(): target_voxel_side_(0) { } -void ParticleGrid::build_grid() +void ParticleGrid::build_grid(const Particles &ps) { + audit_particles(ps); IMP_LOG(TERSE, "Creating nonbonded grid..." << std::flush); float mn[3]= {std::numeric_limits<float>::max(), std::numeric_limits<float>::max(), @@ -38,19 +39,19 @@ float mx[3]={-std::numeric_limits<float>::max(), -std::numeric_limits<float>::max(), -std::numeric_limits<float>::max()}; - for (unsigned int i = 0; i < mc_->get_particles().size(); ++i) { - XYZDecorator d= XYZDecorator::cast(mc_->get_particles()[i]); + for (unsigned int i = 0; i < ps.size(); ++i) { + XYZDecorator d(ps[i]); for (unsigned int j=0; j<3; ++j) { if (d.get_coordinate(j)< mn[j]) mn[j]= d.get_coordinate(j); if (d.get_coordinate(j)> mx[j]) mx[j]= d.get_coordinate(j); } } - if (!mc_->get_particles().empty()) { + if (!ps.empty()) { // keep the grid size sane if things blow up float maxdim= std::max(mx[0]-mn[0], std::max(mx[1]-mn[1], mx[2]-mn[2])); float vx= std::pow(static_cast<float>(target_cell_occupancy *(maxdim*maxdim*maxdim - /mc_->get_particles().size())), + /ps.size())), .3333f); if (vx > target_voxel_side_) { IMP_LOG(VERBOSE, "Overroade target side of " << target_voxel_side_ @@ -62,89 +63,16 @@ Vector3D(mn[0], mn[1], mn[2]), Vector3D(mx[0], mx[1], mx[2]), Particles()); - for (unsigned int i = 0; i < mc_->get_particles().size(); ++i) { - XYZDecorator d= XYZDecorator::cast(mc_->get_particles()[i]); + for (unsigned int i = 0; i < ps.size(); ++i) { + XYZDecorator d(ps[i]); Vector3D v(d.get_x(), d.get_y(), d.get_z()); - grid_.get_voxel(grid_.get_index(v)).push_back(mc_->get_particles()[i]); + grid_.get_voxel(grid_.get_index(v)).push_back(ps[i]); } - grid_valid_=true; - mc_->reset(); IMP_LOG(TERSE, "done." << std::endl); } -void ParticleGrid::add_particles(const Particles &ps) -{ - audit_particles(ps); - mc_->add_particles(ps); - for (unsigned int i=0; i< ps.size(); ++i) { - if (grid_valid_) { - add_particle_to_grid(ps[i]); - } - } -} -void ParticleGrid::add_particle(Particle *p) -{ - Particles ps(1, p); - audit_particles(ps); - mc_->add_particle(p); - - if (grid_valid_) { - add_particle_to_grid(p); - } -} - -void ParticleGrid::add_particle_to_grid(Particle *p) -{ - IMP_assert(grid_valid_, "Bad call of add particle to grid"); - XYZDecorator d= XYZDecorator::cast(p); - Vector3D v(d.get_x(), d.get_y(), d.get_z()); - Grid::VirtualIndex vi= grid_.get_virtual_index(v); - Grid::Index gi= grid_.get_index(vi); - if (gi== Grid::Index()) { - IMP_LOG(TERSE, "Adding particle off grid invalidates it " - << v << " " << vi << std::endl); - grid_valid_=false; - grid_ = Grid(); - } else { - grid_.get_voxel(gi).push_back(p); - } -} - -void ParticleGrid::clear_particles() -{ - mc_->clear_particles(); -} - -bool ParticleGrid::update() -{ - bool ret; - if (!grid_valid_ || mc_->get_max() > target_voxel_side_) { - IMP_LOG(TERSE, "Rebuilding particle grid\n"); - build_grid(); - ret= true; - } else { - IMP_LOG(TERSE, "Removing inactive particles\n"); - for (Grid::DataIterator dit= grid_.data_begin(); dit != grid_.data_end(); - ++dit) { - remove_inactive_particles(*dit); - } - ret= false; - } - unsigned int ssz=0; - for (Grid::DataIterator dit= grid_.data_begin(); dit != grid_.data_end(); - ++dit) { - ssz+= dit->size(); - } - // do this last since it has the ref counts - mc_->update(); - - IMP_assert(ssz== mc_->number_of_particles(), "Particle mismatch in PG"); - - return ret; -} - void ParticleGrid::audit_particles(const Particles &ps) const { for (unsigned int i=0; i< ps.size(); ++i) {

5 57

Grandson of mongo patch
by Daniel Russel 14 Aug '08

14 Aug '08

I have fixed a few bugs so I thought I should submit things. This email lists the changes in addition to the un checked in ones from the last patch. API Changes: - UnaryFunctions::evaluate_derivative now returns a pair of Floats alleviating the need to use the swig hack and making the C++ and swig syntax identical. Python code doesn't need to change but C++ code does. I suggest using boost::tie(value, deriv) (from boost/tuple/ tuple.hpp) to capture the return values for minimal code change. This one glues lots of the disparate changes together. - the TransformPairScore has changed a bit to now functions take Vector3Ds instead of three floats - The randomize functions have been moved out of XYZDecorator into Vector3D.h. This is a more generally useful location for them (which did not exist when I first wrote them). Calls to XYZDecorator::randomize_in_box should be replaced with calls to XYZDecorator::set_coordinates(random_vector_in_box(lb, ub)) - per Ben's suggestion, the hierarchy_ and molecular_hierarchy_ prefixes were removed from some function names Bug fixes: - fixed the get_foos documentation - fixed the hierarchy_get_internal_bonds function, which was not finding all bonds before - circular ref counted pointers with Optimizers/optimizer states. Likewise for Model/Particle, Model/ScoreState (for when Model is refcounted) - BallMover was broken in the IMP repository. I shouldn't return things by passing references to functions as one can easily loose the return value. Additions: - There is a restraint which uses the Lowest n of a set of pair scores. - an is_nan function which uses the built in one when compiling with gcc and the lack of self equality check otherwise. Various places use it now. - a ParticleRefinerGuard which makes sure the refine and cleanup functions are called in pairs by calling cleanup in its destructor. - The BondCoverPartirticleRefiner takes an IntKey whose value is propagated to the particles for use with the typed pair score Not submitted things that might be useful: - A mover to do kmeans clustering. This perhaps should be split into a score state for storing the clustering and a mover to move things to the center. - a particle refiner that uses a lookup table - a pair score which refines pairs of particles which are close under some metric (can be used for hierarchical collision detection - a pair score on the distance between two cylinders - an optimization harness which runs various copies of your optimization (in parallel if you have multiple cores at your disposal) and keeps track of the best results. The rest is various cleanups and new test cases and some extra checks to avoid misuse. I also think that the BondCoverParticleRefiner should be removed as it is not as interesting as hoped.

1 0

bug iterating over restraints
by Keren Lasker 13 Aug '08

13 Aug '08

hi, For some reason iterating over restraints like: for(Restraints::iterator it = opt_mdl->get_restraints().begin(); it != opt_mdl->get_restraints().end(); it++) { } fails as *it == NULL, even though it contains data. The problem is solved if iterating like: for(unsigned int i =0;i<opt_mdl->get_restraints().size();i++){} Keren.

3 5