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::Vectorinternal::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) {