April 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

IMP todo list
by Daniel Russel 30 Apr '08

30 Apr '08

I added a link to a todo page on the IMP wiki to keep track of the things which should get done eventually.

1 0

PairChainRestraint
by Daniel Russel 15 Apr '08

15 Apr '08

Added an analog of TripletChainRestraint which acts on pairs in the chain.

1 0

decorators
by Daniel Russel 13 Apr '08

13 Apr '08

Another source of excess checks in the IMP code is the decorators. There are three reasons to create them - you have a new particle and need to add fields to it - you have a particle which may or may not have the needed fields - you have a particle which has the needed fields (you know already) We currently support the first two and have no way of expressing the third. In addition, the current syntax is kind of annoying. I propose adding constructors to the decorators from Particle* to express the third. The constructors would have an assertion that they particles have the correct attributes (but this assertion would go away in non debugged code). Thoughts?

2 2

Nightly builds failed
by Ben Webb 12 Apr '08

12 Apr '08

See https://salilab.org/internal/imp/tests.html Or more specifically, https://salilab.org/internal/imp/logs/imp/bin.i386-w32.log (since I've already addressed the other failures). Looks like MSVC doesn't like some part of Daniel's latest patch. I'm not super-familiar with the Boost traits stuff (which looks like the culprit). Anybody have any idea what may be wrong? Ben -- ben(a)salilab.org http://salilab.org/~ben/ "It is a capital mistake to theorize before one has data." - Sir Arthur Conan Doyle

2 2

pyext sconscript
by Daniel Russel 10 Apr '08

10 Apr '08

make the libs change use a list so that it plays nicely with other changes to the lib path. Before if the LIBS var was non-empty, the last entry and ../src became one argument.

1 0

IMP patches
by Daniel Russel 09 Apr '08

09 Apr '08

First some minor cleanups - adding operator<< to Vector3D and fixing some header inclusion issues Second adding particle removal as previous discussed: - There is a new base class RefCountedObject which is used for Particle and is ref counted. Vector and ObjectContainer ref count ref counted objects that are contained in them - particles now have a Particle attribute type - the usual caveats about changes to SConscript files being ignore and IMP.i containing some extra changes apply (extra %includes) - if a particle is removed from the model it is made inactive also. We should probably have the particle constructor make the particle inactive until it is added to a model. Bonds use the Particle attribute: - no visible changes Hierarchy uses the particle attribute - no visible changes The nonbonded list is more dynamic: - it drops inactive particles silently - you can add particles and they are handled correctly, if not terribly efficiently I can't guarantee that the patches all work independently

1 0

deactivating particles
by Daniel Russel 08 Apr '08

08 Apr '08

To make my previous email more concrete (and revise things slightly), deactivating particles would involve: - telling the Model to delete the particle. That removes the pointer in the model and frees the index up to be recycled - the ParticleIterators in the model skip delete particles (skip null pointers) - the particle is marked as invalid and any accesses to it throw exceptions - all pointers to the particle are reference counted and the particle object is deleted when all of them go away A slightly more complicated alternative which would delay recycling ParticleIndexes would be - marking them as inactive - The Model ParticleInterators skip inactive (and NULL) particles. Python people will just have to deal unless someone knows how to get custom iterators working in python - all pointers held by restraints are reference counted. When all of the pointers held by everything go away, the Particle entry is set to NULL in the model (and can be reused) and the particle is deleted. - the Model maintains a list of free particle indexes and reuses them once they are gone from the system These proposals leave open the issue of making sure that particle indexes stored in the attribute table always point to the right generation of particle (i.e. make sure that the index doesn't get recycled before they get cleaned up). One solution to this would be to add a new attribute type, that of a particle pointer so these have full reference counted semantics. This would not be hard.

3 2

Inactive particles
by Daniel Russel 08 Apr '08

08 Apr '08

What is means for a particle to be inactive is not well defined at this point. My proposal would be that we define them to be permanently dead-- particles are made inactive by removing them from the Model. When restraints encounter inactive particles, they should either: - if they act on bulk sets of particles, just remove the inactive particle from their list and go on. This would apply to the nonbonded list and list-based restraints - if they act on a constant number of particles, the InactiveParticleException should propagate. The thing holding the restraint can use the exception to destroy the restraint (maybe) or, probably better, just to kill the optimization. The iterators provided by the model would skip inactive particles (and reference counting would ensure that the memory is reclaimed when everything is cleaned up), thereby allowing the optimizers to work when there are inactive particles. Using this, you could, for example, maintain a dynamic set of particles by adding all new particles to the nonbonded list and know that all inactive particles will get cleaned up properly. In order to make sure the reference counting works, Particles would get a new attribute type which is a pointer to another particle. Bonds and hierarchies would use this new attribute type. Note that swig reference counting is pretty much useless, so there is no trivial way to ensure that all references held in python keep the respective bit of memory alive. Namely, if you get a particle back from the model or a restraint, this python object will not keep the corresponding C++ object alive. There are solutions to this problem, but I don't want to implement them now as they are a bit complicated. For now, the best solution might be to have the python wrapper print warnings when the dangerous functions are used.

1 0

Moving particle checks to python
by Daniel Russel 04 Apr '08

04 Apr '08

The first patch moves the checks that particles have a requested attribute out of the C++ runtime pathway. They are kept in the C++ debug builds and are added to the function calls from python The second adds an offset to Linear. It is not really necessary, but makes printouts of the energy must easier to understand. And shouldn't add significant overhead.

1 0