May 2008 - IMP-dev - salilab.org

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

Turning on more refcounting
by Daniel Russel 30 May '08

30 May '08

Here is a patch which turns on refcounting for most of the objects (Optimizers and Model being the two exceptions, although, now that I think about, model should have it (but I would have to change ObjectPointer in order to do so, since lots of things have pointers back to model which can't keep the model alive). It shouldn't change anything for python users other than you can add restraints to multiple restraint sets. We should add a check for dependency cycles in restraint sets, but that takes some work. It also adds a trivial restraint ConstantRestraint. Mostly there for testing (since otherwise the ref count check gets buried under the work of creating all the bits needed for a restraint).

1 0

valgrind suppressions for 64 bit machines
by Daniel Russel 30 May '08

30 May '08

I get a bunch in modeller too which I can add. Index: tools/valgrind-python.supp =================================================================== --- tools/valgrind-python.supp (revision 584) +++ tools/valgrind-python.supp (working copy) @@ -13,6 +13,18 @@ } { + ADDRESS_IN_RANGE/Invalid read of size 8 + Memcheck:Addr8 + fun:PyObject_Free +} + +{ + ADDRESS_IN_RANGE/Invalid read of size 8 + Memcheck:Value8 + fun:PyObject_Free +} + +{ ADDRESS_IN_RANGE/Conditional jump or move depends on uninitialised value Memcheck:Cond fun:PyObject_Free @@ -31,6 +43,18 @@ } { + ADDRESS_IN_RANGE/Invalid read of size 8 + Memcheck:Addr8 + fun:PyObject_Realloc +} + +{ + ADDRESS_IN_RANGE/Invalid read of size 8 + Memcheck:Value8 + fun:PyObject_Realloc +} + +{ ADDRESS_IN_RANGE/Conditional jump or move depends on uninitialised value Memcheck:Cond fun:PyObject_Realloc

3 6

Improving the test framework
by Daniel Russel 29 May '08

29 May '08

Is it possible to have unittest run each test in a separate process? That we we don't have to worry about bad memory accesses in one test crashing python and messing up later tests (and hence can leave failed tests in until we can actually fix them). We are kind of abusing unittest (by using it to test C++ code) so it may not be easy to do.

3 3

Actual patches
by Daniel Russel 28 May '08

28 May '08

First, the patch fixes a bug in remove_always and removes the workaround in Particle. The second fixes a bug in reference counting and removes the workaround in the test code. In addition, for some reason all the reference counting tests were disabled so that fact that they all were failing didn't come up :-) Anyway, now things work and are tested. The containers now should deal with exceptions better (for example adding a non-reference counted Object to two containers and catching the resulting exception doesn't result in memory corruption any more). Not that the rest of the IMP code is all exception safe, but these changes should handle most errors. Index: kernel/include/IMP/internal/AttributeTable.h =================================================================== --- kernel/include/IMP/internal/AttributeTable.h (revision 584) +++ kernel/include/IMP/internal/AttributeTable.h (working copy) @@ -199,7 +199,9 @@ } void remove_always(Key k) { - map_[k.get_index()]= Traits::get_invalid(); + if (k.get_index() < size_) { + map_[k.get_index()]= Traits::get_invalid(); + } } Index: kernel/include/IMP/Particle.h =================================================================== --- kernel/include/IMP/Particle.h (revision 584) +++ kernel/include/IMP/Particle.h (working copy) @@ -265,6 +265,9 @@ \return true it the particle is active. */ bool get_is_active() const { + IMP_IF_CHECK(EXPENSIVE) { + assert_is_valid(); + } return is_active_; } @@ -543,9 +546,7 @@ { floats_.remove(name); derivatives_.remove(name); - if (optimizeds_.contains(name)) { - optimizeds_.remove_always(name); - } + optimizeds_.remove_always(name); } Index: kernel/include/IMP/internal/ObjectContainer.h =================================================================== --- kernel/include/IMP/internal/ObjectContainer.h (revision 584) +++ kernel/include/IMP/internal/ObjectContainer.h (working copy) @@ -9,6 +9,7 @@ #define __IMP_OBJECT_CONTAINER_H #include "Object.h" +#include "RefCountedObject.h" #include <boost/iterator/filter_iterator.hpp> @@ -31,8 +32,9 @@ couldn't get the casts to work out right. */ template <class O, class I> -class ObjectContainer: public std::vector<O*> +class ObjectContainer { + std::vector<O*> data_; std::vector<int> free_; struct OK { bool operator()(const O*a) const { @@ -44,108 +46,111 @@ unsigned int get_index(II i) const {return i.get_index();} unsigned int get_index(unsigned int i) const {return i;} - void ref(O*o) { - if (o) o->ref(); - } - - void unref(O* o) { - if (o) { - o->unref(); - if (!o->get_has_ref()) { - delete o; - } - } - } // hide it - void erase(){} + typedef std::vector<O*> Vector; public: - typedef std::vector<O*> Vector; - using Vector::size; - using Vector::empty; + ObjectContainer(){} ~ObjectContainer() { clear(); } + bool empty() const { + return data_.empty() || data_.size() == free_.size(); + } + unsigned int size() const { + return data_.size() - free_.size(); + } + void clear() { - for (typename Vector::iterator it= Vector::begin(); - it != Vector::end(); ++it) { - unref(*it); + for (typename Vector::iterator it= data_.begin(); + it != data_.end(); ++it) { + O* t= *it; + *it=NULL; + if (t) disown(t); } free_.clear(); - Vector::clear(); + data_.clear(); } typedef boost::filter_iterator<OK, typename Vector::iterator> iterator; - iterator begin() {return iterator(Vector::begin(), Vector::end());} - iterator end() {return iterator(Vector::end(), Vector::end());} + iterator begin() {return iterator(data_.begin(), data_.end());} + iterator end() {return iterator(data_.end(), data_.end());} typedef boost::filter_iterator<OK, typename Vector::const_iterator> const_iterator; const_iterator begin() const { - return const_iterator(Vector::begin(), Vector::end()); + return const_iterator(data_.begin(), data_.end()); } const_iterator end() const { - return const_iterator(Vector::end(), Vector::end()); + return const_iterator(data_.end(), data_.end()); } void remove(I i) { unsigned int id= get_index(i); - IMP_assert(Vector::operator[](id) != NULL, "Nothing there to remove"); - unref(Vector::operator[](id)); - Vector::operator[](id)=NULL; + IMP_assert(id < data_.size(), + "Trying to remove invalid element in container"); + IMP_assert(data_[id] != NULL, "Nothing there to remove"); + O* t= data_[id]; free_.push_back(id); + data_[id]=NULL; + disown(t); } O* operator[](I i) const { - IMP_check(get_index(i) < Vector::size(), + IMP_check(get_index(i) < data_.size(), "Index " << i << " out of range", IndexException("Out of range")); - IMP_assert(Vector::operator[](get_index(i)) != NULL, + IMP_assert(data_.operator[](get_index(i)) != NULL, "Attempting to access invalid slot in container"); - return Vector::operator[](get_index(i)); + return data_[get_index(i)]; } + I push_back(O* d) { IMP_CHECK_OBJECT(d); - ref(d); + own(d); IMP_IF_CHECK(EXPENSIVE) { - for (typename Vector::const_iterator it= Vector::begin(); - it != Vector::end(); ++it) { + for (typename Vector::const_iterator it= data_.begin(); + it != data_.end(); ++it) { IMP_assert(*it != d, "IMP Containers can only have one copy of " << " each object"); } } if (free_.empty()) { - Vector::push_back(d); - unsigned int idx= Vector::size()-1; + data_.push_back(d); + unsigned int idx= data_.size()-1; return I(idx); } else { unsigned int i= free_.back(); free_.pop_back(); - Vector::operator[](i)= d; + data_[i]= d; return I(i); } } + template <class It> void insert(iterator c, It b, It e) { + IMP_assert(c== end(), "Insert position is ignored in ObjectContainer"); IMP_IF_CHECK(EXPENSIVE) { for (It cc= b; cc != e; ++cc) { IMP_CHECK_OBJECT(*cc); - for (typename Vector::const_iterator it= Vector::begin(); - it != Vector::end(); ++it) { + for (typename Vector::const_iterator it= data_.begin(); + it != data_.end(); ++it) { IMP_assert(*it != *cc, "IMP Containers can only have one copy of " << " each object"); } } } for (It cc= b; cc != e; ++cc) { - ref(*cc); + own(*cc); } while (!free_.empty()) { - push_back(*b); + int i= free_.back(); + free_.pop_back(); + data_[i]= *b; ++b; } - Vector::insert(c.base(), b, e); + data_.insert(data_.end(), b, e); } }; Index: kernel/include/IMP/internal/Object.h =================================================================== --- kernel/include/IMP/internal/Object.h (revision 584) +++ kernel/include/IMP/internal/Object.h (working copy) @@ -9,7 +9,6 @@ #ifndef __IMP_OBJECT_H #define __IMP_OBJECT_H -#include "../log.h" #include "../exception.h" namespace IMP @@ -42,24 +41,27 @@ void assert_is_valid() const; bool get_has_ref() const {return count_ != 0;} + void ref() { assert_is_valid(); - IMP_assert(count_== 0, - "Non-reference counted objects can only belong to " \ - "one container."); ++count_; } + void unref() { assert_is_valid(); - IMP_assert(count_ ==1, "Too many unrefs on object"); + IMP_assert(count_ !=0, "Too many unrefs on object"); --count_; } - protected: - int count_; -private: + + unsigned int get_ref_count() const { + return count_; + } + + private: Object(const Object &o){} const internal::Object& operator=(const Object &o) {return *this;} + int count_; double check_value_; }; Index: kernel/include/IMP/internal/ObjectPointer.h =================================================================== --- kernel/include/IMP/internal/ObjectPointer.h (revision 584) +++ kernel/include/IMP/internal/ObjectPointer.h (working copy) @@ -11,6 +11,7 @@ #include "../log.h" #include "Object.h" +#include "RefCountedObject.h" #include "../macros.h" #include "../exception.h" @@ -35,6 +36,16 @@ typedef ObjectPointer<O, RC> This; O* o_; + void set_pointer(O* p) { + if (RC) { + if (o_) disown(o_); + if (p) own(p); + o_=p; + } else { + o_=p; + } + } + // Enforce that ref counted objects are ref counted BOOST_STATIC_ASSERT((RC || !boost::is_base_of<RefCountedObject, O>::value)); @@ -48,34 +59,21 @@ } typedef bool (This::*unspecified_bool)() const; - void ref() { - if (RC && o_) o_->ref(); - } - - void unref() { - if (RC && o_) { - o_->unref(); - if (!o_->get_has_ref()) delete o_; - } - } - public: - ObjectPointer(const ObjectPointer &o): o_(o.o_) { - ref(); + ObjectPointer(const ObjectPointer &o): o_(NULL) { + set_pointer(o.o_); } ObjectPointer& operator=(const ObjectPointer &o){ - unref(); - o_=o.o_; - ref(); + set_pointer(o.o_); return *this; } ObjectPointer(): o_(NULL) {} - explicit ObjectPointer(O* o): o_(o) { - IMP_assert(o != NULL, "Can't initialize with NULL pointer"); - ref(); + explicit ObjectPointer(O* o): o_(NULL) { + IMP_assert(o, "Can't initialize with NULL pointer"); + set_pointer(o); } ~ObjectPointer(){ - unref(); + set_pointer(NULL); } const O& operator*() const { audit(); @@ -98,9 +96,7 @@ return o_; } void operator=(O* o) { - unref(); - o_=o; - ref(); + set_pointer(o); } IMP_COMPARISONS_1(o_); Index: kernel/include/IMP/internal/RefCountedObject.h =================================================================== --- kernel/include/IMP/internal/RefCountedObject.h (revision 584) +++ kernel/include/IMP/internal/RefCountedObject.h (working copy) @@ -22,12 +22,8 @@ //! Common base class for ref counted objects. /** Currently the only ref counted objects are particles. + This class acts as a tag rather than providing any functionality. - \note Due to weirdness in SWIG, the external objects - are responsible for deleting the ref counted object when - the ref count goes to zero. This will change once we have - a real solution for Python. - \internal */ class IMPDLLEXPORT RefCountedObject: public Object @@ -43,21 +39,12 @@ public: virtual ~RefCountedObject() { - IMP_assert(count_==0, "Deleting object which still has references"); + IMP_assert(!get_has_ref(), "Deleting object which still has references"); --live_objects_; } - void ref() { - assert_is_valid(); - ++P::count_; - } - void unref() { - assert_is_valid(); - IMP_assert(count_ != 0, "Too many unrefs on object"); - --P::count_; - } - static unsigned int get_number_of_live_objects() { + // for debugging purposes only return live_objects_; } }; @@ -68,9 +55,9 @@ struct Ref { template <class O> - static void eval(O) { - BOOST_STATIC_ASSERT((!boost::is_pointer<O>::value - || !boost::is_base_of<RefCountedObject, O >::value)); + static void eval(O* o) { + BOOST_STATIC_ASSERT((!boost::is_base_of<RefCountedObject, O >::value)); + IMP_LOG(VERBOSE, "Not refing particle " << o << std::endl); } }; @@ -79,7 +66,10 @@ { template <class O> static void eval(O* o) { - if (o) o->ref(); + IMP_LOG(VERBOSE, "Refing particle " << o->get_index() + << o->get_ref_count() << std::endl); + o->assert_is_valid(); + o->ref(); } }; @@ -87,9 +77,9 @@ struct UnRef { template <class O> - static void eval(O) { - BOOST_STATIC_ASSERT((!boost::is_pointer<O>::value - || !boost::is_base_of<RefCountedObject, O >::value)); + static void eval(O* o) { + BOOST_STATIC_ASSERT((!boost::is_base_of<RefCountedObject, O >::value)); + IMP_LOG(VERBOSE, "Not Unrefing object " << o << std::endl); } }; @@ -98,22 +88,24 @@ { template <class O> static void eval(O *o) { - if (o) { - o->unref(); - if (!o->get_has_ref()) { - delete o; - } + IMP_LOG(VERBOSE, "Unrefing particle " << o->get_index() + << " " << o->get_ref_count() << std::endl); + o->assert_is_valid(); + o->unref(); + if (!o->get_has_ref()) { + delete o; } } -}; + }; //! Can be called on any object and will only unref it if appropriate template <class O> void unref(O o) { - UnRef<(boost::is_pointer<O>::value - && boost::is_base_of<RefCountedObject, O >::value)>::eval(o); + BOOST_STATIC_ASSERT(!boost::is_pointer<O>::value); + /*IMP_LOG(VERBOSE, "NonUnRef called with nonpointer for " + << o << std::endl);*/ } @@ -121,11 +113,69 @@ template <class O> void ref(O o) { - Ref<(boost::is_pointer<O>::value - && boost::is_base_of<RefCountedObject, O >::value)>::eval(o); + BOOST_STATIC_ASSERT(!boost::is_pointer<O>::value); + /*IMP_LOG(VERBOSE, "NonRef count called with nonpointer for " + << o << std::endl);*/ } +//! Can be called on any object and will only unref it if appropriate +template <class O> +void unref(O* o) +{ + /*IMP_LOG(VERBOSE, "Unref count called with " + << (boost::is_base_of<RefCountedObject, O >::value) + << " for " << o << std::endl);*/ + UnRef<(boost::is_base_of<RefCountedObject, O >::value)>::eval(o); +} + +//! Can be called on any object and will only ref it if appropriate +template <class O> +void ref(O* o) +{ + /*IMP_LOG(VERBOSE, "Ref called with " + << (boost::is_base_of<RefCountedObject, O >::value) + << " for " << o << std::endl);*/ + Ref<(boost::is_base_of<RefCountedObject, O >::value)>::eval(o); +} + + +//! Can be called on any object and will only unref it if appropriate +template <class O> +void disown(O* o) +{ + /*IMP_LOG(VERBOSE, "Disown called with " + << (boost::is_base_of<RefCountedObject, O >::value) + << " for " << o << " " << o->get_ref_count() << std::endl);*/ + o->unref(); + if (!o->get_has_ref()) { + delete o; + } +} + + +//! Can be called on any object and will only ref it if appropriate +template <class O> +void own(O* o) +{ + /*IMP_LOG(VERBOSE, "Own called with " + << (boost::is_base_of<RefCountedObject, O >::value) + << " for " << o + << " " << o->get_ref_count() << std::endl);*/ + if (boost::is_base_of<RefCountedObject, O >::value) { + // no checks + } else { + IMP_check(!o->get_has_ref(), "Trying to own already owned but " + << "non-reference-counted object.", + ValueException("Already owned object")); + } + o->ref(); +} + + + + + } // namespace internal } // namespace IMP Index: kernel/test/states/test_nonbonded_list.py =================================================================== --- kernel/test/states/test_nonbonded_list.py (revision 584) +++ kernel/test/states/test_nonbonded_list.py (working copy) @@ -165,7 +165,7 @@ score= m.evaluate(False) self.assertEqual(score, 2*190, "Wrong score") - ps[3].set_is_active(False) + m.remove_particle(ps[3].get_index()) self.assert_(not ps[3].get_is_active(), "Particle not inactive") ps=None score= m.evaluate(False) Index: kernel/test/particles/test_refcount.py =================================================================== --- kernel/test/particles/test_refcount.py (revision 584) +++ kernel/test/particles/test_refcount.py (working copy) @@ -7,15 +7,23 @@ class RefCountTests(IMP.test.TestCase): """Test refcounting of particles""" + def setUp(self): + IMP.test.TestCase.setUp(self) + IMP.set_log_level(IMP.VERBOSE) + self.basenum= IMP.RefCountedObject.get_number_of_live_objects() + print "The base number of objects is " + str(self.basenum) + def _check_number(self, expected): print "Expected "+str(expected)\ - + " got " + str(IMP.RefCountedObject.get_number_of_live_objects()) - self.assertEqual(IMP.RefCountedObject.get_number_of_live_objects(), + + " got " + str(IMP.RefCountedObject.get_number_of_live_objects() + - self.basenum) + self.assertEqual(IMP.RefCountedObject.get_number_of_live_objects() - self.basenum, expected, "wrong number of particles") - def _test_simple(self): + def __test_simple(self): """Check that ref counting of particles works within python""" + # swig is broken so this needs to be skipped self._check_number(0) p= IMP.Particle() @@ -33,7 +41,7 @@ m=1 self._check_number(0) - def _test_removal(self): + def test_removal(self): """Check that ref counting works with removing particles""" self._check_number(0) m= IMP.Model() @@ -47,7 +55,7 @@ self._check_number(0) # This test does not work since swig refcounting is broken - def _test_round_trip(self): + def __test_round_trip(self): """test that tracking survives the round trip""" print "test that the round trip object doesn't delete it" @@ -101,24 +109,30 @@ self._check_number(0) - def _test_shared(self): - """Check that ref counting works shared particles""" + def test_shared(self): + """Check that ref counting works with shared particles""" print "max change" self._check_number(0) m= IMP.Model() p= IMP.Particle() + print "Add particle" pi= m.add_particle(p) + d= IMP.XYZDecorator.create(p) + d=0 mc= IMP.MaxChangeScoreState(IMP.XYZDecorator.get_xyz_keys()) + print "Add particle to mc" mc.add_particle(p) self._check_number(1) + print "Remove from model" m.remove_particle(pi) self._check_number(1) p=1 self._check_number(1) + print "Remove from mc" mc.clear_particles() self._check_number(0) - def _test_skip(self): + def test_skip(self): """Check that removed particles are skipped""" print "skipped" m= IMP.Model()

1 0

Particle removal and deactivation
by Daniel Russel 28 May '08

28 May '08

I can't find the imp-commit email for some reason. Anyway, there was a recent patch (581) changing a model::remove_particle to deactivating the particle. This change is wrong as the particles are currently reference counted in C++ and removing the particle from the model should do the right thing assuming there are no python references hanging around. If it does not, it is a bug in kernel code which should be fixed. In fact, the "Particle::set_is_active" method should probably go away.

3 7

Patch 583
by Daniel Russel 28 May '08

28 May '08

Patch 583 doesn't really make sense. Either it should become "if (table.contains()) table.remove" and remove always should be removed or remove_always should be fixed to do bounds checking. I think prefer the first.

2 1

Removing particle attributes
by Daniel Russel 25 May '08

25 May '08

The attached patch adds the ability to remove particle attributes. It also adds more thorough checking of NaNs in the float attributes. In addition, the rearrangement creates a natural iterator through the optimized float keys which the Optimizer now uses. Internally, it changes how attributes are store, replacing the valid bit with a sentinel value. This works nicely for floats and derivs where we can can use NaN and for particles where we can use NULL. It also reduces the memory usage for those by a bit since the bools couldn't be packed well before. For ints, it is a little less nice as it currently uses numeric_limits<int>::max as a sentinel. For strings it is a particular weird string. For both of those we can easily improve things if it becomes a problem (either due to people wanting to use int max as an attribute value or the memory allocation associated with the strings becoming significant). Index: kernel/include/IMP/internal/AttributeTable.h =================================================================== --- kernel/include/IMP/internal/AttributeTable.h (revision 576) +++ kernel/include/IMP/internal/AttributeTable.h (working copy) @@ -11,12 +11,14 @@ #include "../base_types.h" #include "../utility.h" #include "../log.h" +#include "ObjectPointer.h" #include <boost/iterator/filter_iterator.hpp> #include <boost/iterator/counting_iterator.hpp> #include <boost/scoped_array.hpp> #include <vector> +#include <limits> namespace IMP { @@ -25,25 +27,101 @@ namespace internal { +struct FloatAttributeTableTraits { + 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(); + } else if (std::numeric_limits<Float>::has_infinity) { + return std::numeric_limits<Float>::infinity(); + } else { + return std::numeric_limits<Float>::max(); + } + } + static bool get_is_valid(Float f) { + if (std::numeric_limits<Float>::has_quiet_NaN) { + return f==f; + } else { + return f!= get_invalid(); + } + } +}; + +struct IntAttributeTableTraits { + typedef Int Value; + typedef KeyBase<Int> Key; + static Int get_invalid() { + return std::numeric_limits<Int>::max(); + } + static bool get_is_valid(Int f) { + return f!= get_invalid(); + } +}; + +struct BoolAttributeTableTraits { + typedef bool Value; + typedef KeyBase<Float> Key; + static bool get_invalid() { + return false; + } + static bool get_is_valid(Int f) { + return f; + } +}; + +struct StringAttributeTableTraits { + typedef String Value; + typedef KeyBase<String> Key; + static Value get_invalid() { + return "This is an invalid string in IMP"; + } + static bool get_is_valid(String f) { + return f!= get_invalid(); + } +}; + +struct ParticleAttributeTableTraits { + typedef internal::ObjectPointer<Particle,true> Value; + typedef KeyBase<Particle*> Key; + static Value get_invalid() { + return Value(); + } + static bool get_is_valid(Value f) { + return f!= Value(); + } +}; + /** \internal If memory becomes a problem then either use a char table to look up actual entries or use perfect hashing. http://burtleburtle.net/bob/hash/perfect.html + + Actuall Cuckoo hashing is probably a better bet as that gets + high occupancy without large tables for the hash function. + + \note This version of the table uses certain values of the data + to singal that the entry is invalid. Setting an entry to these + values is a checked error. The values are specified by the + Traits::invalid entry. */ -template <class T, class VT> +template <class Traits> class AttributeTable { - typedef AttributeTable<T, VT> This; - struct Bin - { - bool first; - VT second; - Bin(): first(false){} - }; - typedef boost::scoped_array<Bin> Map; + typedef AttributeTable<Traits> This; + + typedef boost::scoped_array<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"); @@ -60,15 +138,25 @@ map_.reset(); } else { size_=std::max(nlen, 6U); - map_.reset(new Bin[size_]); + 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(), + "Attribute \"" << k.get_string() + << "\" not found in table."); + IMP_check(Traits::get_is_valid(map_[k.get_index()]), + "Attribute \"" << k.get_string() + << "\" not found in table.", + IndexException("Invalid attribute")); + } + public: - typedef VT Value; - typedef KeyBase<T> Key; + 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; @@ -87,28 +175,43 @@ return *this; } const Value get_value(Key k) const { - IMP_assert(contains(k), - "Attribute \"" << k.get_string() - << "\" not found in table."); - return map_[k.get_index()].second; + check_contains(k); + return map_[k.get_index()]; } - Value& get_value(Key k) { - IMP_assert(contains(k), - "Attribute \"" << k.get_string() - << "\" not found in table."); - return map_[k.get_index()].second; + void set_value(Key k, Value v) { + check_contains(k); + IMP_check(Traits::get_is_valid(v), + "Cannot set value of attribute to " << v, + ValueException("Invalid value for attribute")); + map_[k.get_index()] = v; } - void insert(Key k, Value v); + void insert(Key k, Value v) { + IMP_assert(!contains(k), + "Trying to add attribute \"" << k.get_string() + << "\" twice"); + insert_always(k, v); + } + void insert_always(Key k, Value v); + void remove(Key k) { + check_contains(k); + map_[k.get_index()]= Traits::get_invalid(); + } + + void remove_always(Key k) { + map_[k.get_index()]= Traits::get_invalid(); + } + + bool contains(Key k) const { IMP_assert(k != Key(), "Can't search for default key"); return k.get_index() < size_ - && map_[k.get_index()].first; + && Traits::get_is_valid(map_[k.get_index()]); } @@ -143,57 +246,52 @@ KeyIterator(Key(size_))); } - - unsigned int get_heap_memory_usage() const { - return size_*sizeof(Bin); - } - }; -IMP_OUTPUT_OPERATOR_2(AttributeTable) +IMP_OUTPUT_OPERATOR_1(AttributeTable) -template <class T, class VT> -inline void AttributeTable<T, VT>::insert(Key k, Value v) +template <class Traits> +inline void AttributeTable<Traits>::insert_always(Key k, Value v) { /*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 " + << v << " into attribute " << k, + ValueException("Invalid attribute value")); if (size_ <= k.get_index()) { - boost::scoped_array<Bin> old; + boost::scoped_array<Value> old; swap(old, map_); realloc(size_, old, k.get_index()+1); } - IMP_assert(!map_[k.get_index()].first, - "Trying to add attribute \"" << k.get_string() - << "\" twice"); - map_[k.get_index()].second= v; - map_[k.get_index()].first= true; - IMP_assert(contains(k), "Something is broken"); + map_[k.get_index()]= v; } - template <class T, class VT> - inline void AttributeTable<T, VT>::show(std::ostream &out, - const char *prefix) const + +template <class Traits> +inline void AttributeTable<Traits>::show(std::ostream &out, + const char *prefix) const { for (unsigned int i=0; i< size_; ++i) { - if (map_[i].first) { + if (Traits::get_is_valid(map_[i])) { out << prefix; out << Key(i).get_string() << ": "; - out << map_[i].second; + out << map_[i]; out << std::endl; } } } -template <class T, class VT> -inline std::vector<typename AttributeTable<T, VT>::Key> - AttributeTable<T, VT>::get_keys() const +template <class Traits> +inline std::vector<typename Traits::Key> + AttributeTable<Traits>::get_keys() const { std::vector<Key> ret(attribute_keys_begin(), attribute_keys_end()); return ret; Index: kernel/include/IMP/Particle.h =================================================================== --- kernel/include/IMP/Particle.h (revision 576) +++ kernel/include/IMP/Particle.h (working copy) @@ -26,25 +26,6 @@ namespace internal { -struct IMPDLLEXPORT FloatData -{ - bool is_optimized; - Float value; - Float derivative; - FloatData(Float v, bool opt): is_optimized(opt), value(v), derivative(0){} - FloatData():is_optimized(false), value(std::numeric_limits<Float>::max()), - derivative(std::numeric_limits<Float>::max()){} -}; - -inline std::ostream &operator<<(std::ostream &out, const FloatData &d) -{ - out << d.value; - if (d.is_optimized) { - out << " (" << d.derivative << ")"; - } - return out; -} - template <class It> void check_particles_active(It b, It e, std::string msg) { @@ -61,17 +42,29 @@ //! Class to handle individual model particles. /** This class contains particle methods and indexes to particle attributes. - Particles cannot be deleted once they are added to a model, but they can - be deactivated (with their set_is_active method) after which they play no - role in the scoring (it is illegal to try to evaluate a restraint on an - inactive particle). To merely prevent a particle from moving during + To merely prevent a particle from moving during optimization, mark all of its attributes as being non-optimizable (set_is_optimized method). + \ingroup kernel */ class IMPDLLEXPORT Particle : public internal::RefCountedObject { friend class Model; + + typedef internal::AttributeTable<internal::FloatAttributeTableTraits> + FloatTable; + typedef internal::AttributeTable<internal::FloatAttributeTableTraits> + DerivativeTable; + typedef internal::AttributeTable<internal::BoolAttributeTableTraits> + OptimizedTable; + typedef internal::AttributeTable<internal::IntAttributeTableTraits> + IntTable; + typedef internal::AttributeTable<internal::StringAttributeTableTraits> + StringTable; + typedef internal::AttributeTable<internal::ParticleAttributeTableTraits> + ParticleTable; + public: @@ -100,6 +93,12 @@ void add_attribute(FloatKey name, const Float value, const bool is_optimized = false); + //! Remove a Float attribute from this particle. + /** \param[in] name Name of the attribute being added. + */ + void remove_attribute(FloatKey name); + + //! Does particle have a Float attribute with the given name. /** \param[in] name Name of the attribute being checked. \return true if Float attribute exists in this particle. @@ -147,6 +146,11 @@ */ void add_attribute(IntKey name, const Int value); + //! Remove a Int attribute from this particle. + /** \param[in] name Name of the attribute being added. + */ + void remove_attribute(IntKey name); + //! Does particle have an Int attribute with the given name. /** \param[in] name Name of the attribute being checked. \return true if Int attribute exists in this particle. @@ -173,6 +177,12 @@ */ void add_attribute(StringKey name, const String value); + //! Remove a String attribute from this particle. + /** \param[in] name Name of the attribute being added. + */ + void remove_attribute(StringKey name); + + //! Does particle have a String attribute with the given name. /** \param[in] name Name of the attribute being checked. \return true if Int attribute exists in this particle. @@ -200,6 +210,11 @@ */ void add_attribute(ParticleKey name, Particle* value); + //! Remove a Particle attribute from this particle. + /** \param[in] name Name of the attribute being added. + */ + void remove_attribute(ParticleKey name); + //! Does particle have a Particle attribute with the given name. /** \param[in] name Name of the attribute being checked. \return true if Particle attribute exists in this particle. @@ -256,72 +271,78 @@ which seems inelegant. */ std::vector<FloatKey> get_float_attributes() const { - return float_indexes_.get_keys(); + return floats_.get_keys(); } //! See get_float_attributes std::vector<IntKey> get_int_attributes() const { - return int_indexes_.get_keys(); + return ints_.get_keys(); } //! See get_float_attributes std::vector<StringKey> get_string_attributes() const { - return string_indexes_.get_keys(); + return strings_.get_keys(); } //! See get_particle_attributes std::vector<ParticleKey> get_particle_attributes() const { - return particle_indexes_.get_keys(); + return particles_.get_keys(); } //! An iterator through the keys of the float attributes of this particle - typedef internal::AttributeTable<Float, - internal::FloatData>::AttributeKeyIterator + typedef FloatTable::AttributeKeyIterator FloatKeyIterator; //! Iterate through the keys of float attributes of the particle FloatKeyIterator float_keys_begin() const { - return float_indexes_.attribute_keys_begin(); + return floats_.attribute_keys_begin(); } FloatKeyIterator float_keys_end() const { - return float_indexes_.attribute_keys_end(); + return floats_.attribute_keys_end(); } + //! An iterator through the keys of the derivative attributes of this particle + typedef OptimizedTable::AttributeKeyIterator + OptimizedKeyIterator; + //! Iterate through the keys of float attributes of the particle + OptimizedKeyIterator optimized_keys_begin() const { + return optimizeds_.attribute_keys_begin(); + } + OptimizedKeyIterator optimized_keys_end() const { + return optimizeds_.attribute_keys_end(); + } + + //! An iterator through the keys of the int attributes of this particle - typedef internal::AttributeTable<Int, Int> - ::AttributeKeyIterator IntKeyIterator; + typedef IntTable::AttributeKeyIterator IntKeyIterator; //! Iterate through the keys of int attributes of the particle IntKeyIterator int_keys_begin() const { - return int_indexes_.attribute_keys_begin(); + return ints_.attribute_keys_begin(); } IntKeyIterator int_keys_end() const { - return int_indexes_.attribute_keys_end(); + return ints_.attribute_keys_end(); } //! An iterator through the keys of the string attributes of this particle - typedef internal::AttributeTable<String, String> - ::AttributeKeyIterator StringKeyIterator; + typedef StringTable::AttributeKeyIterator StringKeyIterator; //! Iterate through the keys of string attributes of the particle StringKeyIterator string_keys_begin() const { - return string_indexes_.attribute_keys_begin(); + return strings_.attribute_keys_begin(); } StringKeyIterator string_keys_end() const { - return string_indexes_.attribute_keys_end(); + return strings_.attribute_keys_end(); } //! An iterator through the keys of the Particle attributes of this particle - typedef internal::AttributeTable<Particle*, - internal::ObjectPointer<Particle, true> > - ::AttributeKeyIterator ParticleKeyIterator; + typedef ParticleTable::AttributeKeyIterator ParticleKeyIterator; //! Iterate through the keys of Particle attributes of the particle ParticleKeyIterator particle_keys_begin() const { - return particle_indexes_.attribute_keys_begin(); + return particles_.attribute_keys_begin(); } ParticleKeyIterator particle_keys_end() const { - return particle_indexes_.attribute_keys_end(); + return particles_.attribute_keys_end(); } - protected: void zero_derivatives(); @@ -335,14 +356,18 @@ bool is_active_; // float attributes associated with the particle - internal::AttributeTable<Float, internal::FloatData> float_indexes_; + FloatTable floats_; + // float attributes associated with the particle + FloatTable derivatives_; + // Whether a given float is optimized or not + OptimizedTable optimizeds_; + // int attributes associated with the particle - internal::AttributeTable<Int, Int> int_indexes_; + IntTable ints_; // string attributes associated with the particle - internal::AttributeTable<String, String> string_indexes_; + StringTable strings_; // particle attributes associated with the particle - internal::AttributeTable<Particle*, internal::ObjectPointer<Particle,true> > - particle_indexes_; + ParticleTable particles_; ParticleIndex pi_; }; @@ -354,7 +379,7 @@ inline bool Particle::has_attribute(FloatKey name) const { - return float_indexes_.contains(name); + return floats_.contains(name); } @@ -363,36 +388,45 @@ { IMP_check(get_is_active(), "Do not touch inactive particles", InactiveParticleException()); - return float_indexes_.get_value(name).value; + return floats_.get_value(name); } inline Float Particle::get_derivative(FloatKey name) const { IMP_check(get_is_active(), "Do not touch inactive particles", InactiveParticleException()); - return float_indexes_.get_value(name).derivative; + return derivatives_.get_value(name); } inline void Particle::set_value(FloatKey name, Float value) { IMP_check(get_is_active(), "Do not touch inactive particles", InactiveParticleException()); - IMP_assert(value==value, "Can't set value to NaN"); - float_indexes_.get_value(name).value= value; + floats_.set_value(name, value); } inline bool Particle::get_is_optimized(FloatKey name) const { IMP_check(get_is_active(), "Do not touch inactive particles", InactiveParticleException()); - return float_indexes_.get_value(name).is_optimized; + IMP_check(floats_.contains(name), "get_is_optimized called " + << "with invalid attribute" << name, + IndexException("Invalid float attribute")); + return optimizeds_.contains(name); } inline void Particle::set_is_optimized(FloatKey name, bool tf) { IMP_check(get_is_active(), "Do not touch inactive particles", InactiveParticleException()); - float_indexes_.get_value(name).is_optimized=tf; + IMP_check(floats_.contains(name), "set_is_optimized called " + << "with invalid attribute" << name, + IndexException("Invalid float attribute")); + if (tf) { + optimizeds_.insert_always(name, true); + } else { + optimizeds_.remove_always(name); + } } inline void Particle::add_to_derivative(FloatKey name, Float value, @@ -401,14 +435,14 @@ IMP_check(get_is_active(), "Do not touch inactive particles", InactiveParticleException()); IMP_assert(value==value, "Can't add NaN to derivative in particle " << *this); - float_indexes_.get_value(name).derivative+= da(value); + derivatives_.set_value(name, derivatives_.get_value(name)+ da(value)); } inline bool Particle::has_attribute(IntKey name) const { IMP_check(get_is_active(), "Do not touch inactive particles", InactiveParticleException()); - return int_indexes_.contains(name); + return ints_.contains(name); } @@ -417,7 +451,7 @@ { IMP_check(get_is_active(), "Do not touch inactive particles", InactiveParticleException()); - return int_indexes_.get_value(name); + return ints_.get_value(name); } @@ -425,14 +459,14 @@ { IMP_check(get_is_active(), "Do not touch inactive particles", InactiveParticleException()); - int_indexes_.get_value(name)= value; + ints_.set_value(name, value); } inline bool Particle::has_attribute(StringKey name) const { IMP_check(get_is_active(), "Do not touch inactive particles", InactiveParticleException()); - return string_indexes_.contains(name); + return strings_.contains(name); } @@ -441,14 +475,14 @@ { IMP_check(get_is_active(), "Do not touch inactive particles", InactiveParticleException()); - return string_indexes_.get_value(name); + return strings_.get_value(name); } inline void Particle::set_value(StringKey name, String value) { IMP_check(get_is_active(), "Do not touch inactive particles", InactiveParticleException()); - string_indexes_.get_value(name)= value; + strings_.set_value(name, value); } @@ -456,7 +490,7 @@ { IMP_check(get_is_active(), "Do not touch inactive particles", InactiveParticleException()); - return particle_indexes_.contains(name); + return particles_.contains(name); } @@ -465,7 +499,7 @@ { IMP_check(get_is_active(), "Do not touch inactive particles", InactiveParticleException()); - return particle_indexes_.get_value(name).get(); + return particles_.get_value(name).get(); } @@ -473,7 +507,7 @@ { IMP_check(get_is_active(), "Do not touch inactive particles", InactiveParticleException()); - particle_indexes_.get_value(name)= value; + particles_.set_value(name, internal::ObjectPointer<Particle, true>(value)); } @@ -482,34 +516,60 @@ { IMP_assert(model_ , "Particle must be added to Model before attributes are added"); - float_indexes_.insert(name, internal::FloatData(value, is_optimized)); + floats_.insert(name, value); + derivatives_.insert(name, 0); + if (is_optimized) { + optimizeds_.insert(name, true); + } } +void inline Particle::remove_attribute(FloatKey name) +{ + floats_.remove(name); + derivatives_.remove(name); + optimizeds_.remove_always(name); +} + void inline Particle::add_attribute(IntKey name, const Int value) { IMP_assert(model_, "Particle must be added to Model before attributes are added"); - int_indexes_.insert(name, value); + ints_.insert(name, value); } +void inline Particle::remove_attribute(IntKey name) +{ + ints_.remove(name); +} + void inline Particle::add_attribute(StringKey name, const String value) { IMP_assert(model_, "Particle must be added to Model before attributes are added"); - string_indexes_.insert(name, value); + strings_.insert(name, value); } +void inline Particle::remove_attribute(StringKey name) +{ + strings_.remove(name); +} + void inline Particle::add_attribute(ParticleKey name, Particle* value) { IMP_assert(model_, "Particle must be added to Model before attributes are added"); - particle_indexes_.insert(name, + particles_.insert(name, internal::ObjectPointer<Particle, true>(value)); } +void inline Particle::remove_attribute(ParticleKey name) +{ + particles_.remove(name); +} + } // namespace IMP #endif /* __IMP_PARTICLE_H */ Index: kernel/include/IMP/Optimizer.h =================================================================== --- kernel/include/IMP/Optimizer.h (revision 576) +++ kernel/include/IMP/Optimizer.h (working copy) @@ -70,7 +70,6 @@ void update_states(); //! An index to an optimized particle attribute - class FloatIndexInderator; struct FloatIndex { /** @@ -79,7 +78,7 @@ friend class Optimizer; friend class FloatIndexIterator; Model::ParticleIterator p_; - Particle::FloatKeyIterator fk_; + Particle::OptimizedKeyIterator fk_; FloatIndex(Model::ParticleIterator p): p_(p){} public: FloatIndex() {} @@ -96,20 +95,19 @@ mutable FloatIndex i_; void search_valid() const { - while (i_.fk_ == (*i_.p_)->float_keys_end() - || !(*i_.p_)->get_is_optimized(*i_.fk_)) { - if (i_.fk_ == (*i_.p_)->float_keys_end()) { + while (i_.fk_ == (*i_.p_)->optimized_keys_end()) { + if (i_.fk_ == (*i_.p_)->optimized_keys_end()) { ++i_.p_; if (i_.p_== pe_) return; else { - i_.fk_= (*i_.p_)->float_keys_begin(); + i_.fk_= (*i_.p_)->optimized_keys_begin(); } } else { ++i_.fk_; } } IMP_assert(i_.p_ != pe_, "Should have just returned"); - IMP_assert(i_.fk_ != (*i_.p_)->float_keys_end(), + IMP_assert(i_.fk_ != (*i_.p_)->optimized_keys_end(), "Broken iterator end"); IMP_assert((*i_.p_)->get_is_optimized(*i_.fk_), "Why did the loop end?"); @@ -122,7 +120,7 @@ FloatIndexIterator(Model::ParticleIterator pc, Model::ParticleIterator pe): pe_(pe), i_(pc) { if (pc != pe) { - i_.fk_= (*pc)->float_keys_begin(); + i_.fk_= (*pc)->optimized_keys_begin(); search_valid(); } } @@ -162,7 +160,7 @@ FloatIndexIterator float_indexes_end() { return FloatIndexIterator(model_->particles_end(), - model_->particles_end()); + model_->particles_end()); } //! Set the value of an optimized attribute Index: kernel/test/particles/test_particles.py =================================================================== --- kernel/test/particles/test_particles.py (revision 576) +++ kernel/test/particles/test_particles.py (working copy) @@ -68,6 +68,28 @@ p.set_is_active(False) self.assertEqual(p.get_is_active(), False) + def _test_add_remove(self, p, ak, v): + p.add_attribute(ak, v) + self.assert_(p.has_attribute(ak)) + p.remove_attribute(ak) + self.assert_(not p.has_attribute(ak)) + + def test_remove_attributes(self): + """Test that attributes can be removed""" + p=self.particles[0] + fk= IMP.FloatKey("to_remove") + p.add_attribute(fk, 0, False) + self.assert_(p.has_attribute(fk)) + self.assert_(not p.get_is_optimized(fk)) + p.set_is_optimized(fk, True) + self.assert_(p.get_is_optimized(fk)) + p.set_is_optimized(fk, False) + self.assert_(not p.get_is_optimized(fk)) + self._test_add_remove(p, IMP.FloatKey("something"), 1.0) + self._test_add_remove(p, IMP.StringKey("something"), "Hello") + self._test_add_remove(p, IMP.IntKey("something"), 1) + self._test_add_remove(p, IMP.ParticleKey("something"), p) + def test_derivatives(self): """Test get/set of derivatives""" p = self.particles[0] Index: kernel/src/Model.cpp =================================================================== --- kernel/src/Model.cpp (revision 576) +++ kernel/src/Model.cpp (working copy) @@ -75,6 +75,7 @@ "Begin evaluate restraints " << (calc_derivs?"with derivatives":"without derivatives") << std::endl); + for (RestraintIterator it = restraints_begin(); it != restraints_end(); ++it) { IMP_CHECK_OBJECT(*it); @@ -87,34 +88,23 @@ IMP_LOG(TERSE, "Restraint score is " << tscore << std::endl); score+= tscore; } + IMP_LOG(TERSE, "End evaluate restraints." << std::endl); IMP_LOG(TERSE, "Begin after_evaluate of ScoreStates " << std::endl); + for (ScoreStateIterator it = score_states_begin(); it != score_states_end(); ++it) { IMP_CHECK_OBJECT(*it); (*it)->after_evaluate(iteration_, accpt); IMP_LOG(VERBOSE, "." << std::flush); } + IMP_LOG(TERSE, "End after_evaluate of ScoreStates." << std::endl); IMP_LOG(TERSE, "End Model::evaluate. Final score: " << score << std::endl); - for (ParticleIterator pit= particles_begin(); - pit != particles_end(); ++pit) { - for (Particle::FloatKeyIterator fkit = (*pit)->float_keys_begin(); - fkit != (*pit)->float_keys_end(); ++fkit) { - Float v= (*pit)->get_value(*fkit); - IMP_check(v==v, "NaN found in particle " << (*pit)->get_index() - << " attribute " << *fkit, - ValueException("NaN found")); - Float d= (*pit)->get_derivative(*fkit); - IMP_check(d==d, "NaN found in particle derivative " << (*pit)->get_index() - << " attribute " << *fkit, - ValueException("NaN found")); - } - } ++iteration_; return score; } Index: kernel/src/Particle.cpp =================================================================== --- kernel/src/Particle.cpp (revision 576) +++ kernel/src/Particle.cpp (working copy) @@ -44,7 +44,7 @@ void Particle::zero_derivatives() { for (FloatKeyIterator it= float_keys_begin(); it != float_keys_end(); ++it) { - float_indexes_.get_value(*it).derivative=0; + derivatives_.set_value(*it, 0); } } @@ -63,13 +63,17 @@ if (get_model() != NULL) { out << inset << inset << "float attributes:" << std::endl; - float_indexes_.show(out, " "); + floats_.show(out, " "); out << inset << inset << "int attributes:" << std::endl; - int_indexes_.show(out, " "); + ints_.show(out, " "); out << inset << inset << "string attributes:" << std::endl; - string_indexes_.show(out, " "); + strings_.show(out, " "); + + out << inset << inset << "particle attributes:" << std::endl; + particles_.show(out, " "); + } } Index: kernel/src/internal/graph_base.cpp =================================================================== --- kernel/src/internal/graph_base.cpp (revision 576) +++ kernel/src/internal/graph_base.cpp (working copy) @@ -74,7 +74,7 @@ p[i]->set_value(graph_get_edge_key(j+shift, d), v); } } - p[i]->set_value(graph_get_edge_key(nc-1, d), NULL); + p[i]->remove_attribute(graph_get_edge_key(nc-1, d)); IMP_assert(shift==-1, "no edge found"); IMP_assert(nc > 0, "Too few edges"); p[i]->set_value(d.num_edges_key_, nc-1);

2 3

Changing white space in internal/Vector.h
by Daniel Russel 23 May '08

23 May '08

Added some white space to make it easier to read. The email to commit system is definitely overkill for this :-) Index: kernel/include/IMP/internal/Vector.h =================================================================== --- kernel/include/IMP/internal/Vector.h (revision 573) +++ kernel/include/IMP/internal/Vector.h (working copy) @@ -33,30 +33,37 @@ public: typedef std::vector<D> P; Vector(){} + template <class It> Vector(It b, It e) { insert(P::end(), b, e); } + Vector(const P &p) { insert(P::end(), p.begin(), p.end()); } + ~Vector(){clear();} + const D& operator[](unsigned int i) const { IMP_check(i < P::size(), "Index " << i << " out of range", IndexException("")); return P::operator[](i); } + D& operator[](unsigned int i) { IMP_check(i < P::size(), "Index " << i << " out of range", IndexException("")); return P::operator[](i); } + void erase(typename P::iterator it) { unref(*it); P::erase(it); } + void erase(typename P::iterator b, typename P::iterator e) { for (typename P::iterator c= b; c != e; ++c) { @@ -64,21 +71,25 @@ } P::erase(b,e); } + unsigned int push_back(D d) { ref(d); P::push_back(d); return P::size()-1; } + void pop_back() { unref(P::back()); P::pop_back(); } + void clear() { for (typename P::iterator it= P::begin(); it != P::end(); ++it) { unref(*it); } P::clear(); } + template <class It> void insert(typename P::iterator it, It b, It e) { for (It c= b; c != e; ++c) {

1 0

Move exception code
by Daniel Russel 23 May '08

23 May '08

This patch moves the code called on check and exception failures to exception.cpp rather than Log.cpp. Index: kernel/include/IMP/exception.h =================================================================== --- kernel/include/IMP/exception.h (revision 573) +++ kernel/include/IMP/exception.h (working copy) @@ -134,13 +134,13 @@ { //! This is just here so you can catch errors more easily in the debugger -/** Break on Log.cpp:19 to catch assertion failures. +/** Break on exception.cpp:31 to catch assertion failures. \ingroup assert */ IMPDLLEXPORT void assert_fail(); //! Here so you can catch check failures more easily in the debugger -/** Break on Log.cpp:22 to catch check failures. +/** Break on exception.cpp:35 to catch check failures. \ingroup assert */ IMPDLLEXPORT void check_fail(); @@ -178,12 +178,13 @@ \param[in] exception Throw the object constructed by this expression. \ingroup assert */ -#define IMP_check(expr, message, exception) \ - do { \ +#define IMP_check(expr, message, exception) \ + do { \ if (IMP::get_check_level() >= IMP::CHEAP && !(expr)) { \ - IMP_ERROR(message); \ - throw exception; \ - } \ + IMP_ERROR(message); \ + IMP::internal::check_fail(); \ + throw exception; \ + } \ } while (false) //! A runtime failure for IMP. Index: kernel/src/exception.cpp =================================================================== --- kernel/src/exception.cpp (revision 573) +++ kernel/src/exception.cpp (working copy) @@ -27,4 +27,13 @@ return check_mode; } +namespace internal { + void assert_fail() { + throw ErrorException(); + } + + void check_fail() { + } +} + } // namespace IMP Index: kernel/src/Log.cpp =================================================================== --- kernel/src/Log.cpp (revision 573) +++ kernel/src/Log.cpp (working copy) @@ -14,13 +14,4 @@ /* Initialize singleton pointer to NULL */ Log* Log::logpt_ = NULL; - -namespace internal { - void assert_fail() { - throw ErrorException(); - } - - void check_fail() { - } -} } // namespace IMP

1 0