Index: kernel/include/IMP/score_states/QuadraticNonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/QuadraticNonbondedListScoreState.h (revision 0) +++ kernel/include/IMP/score_states/QuadraticNonbondedListScoreState.h (revision 0) @@ -0,0 +1,53 @@ +/** + * \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 +#include + +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 IMP::NonbondedListScoreState +{ + typedef NonbondedListScoreState P; + internal::ObjectPointer 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(); +}; + +//! This class maintains a list of non-bonded pairs of spheres + +} // 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 489) +++ kernel/include/IMP/score_states/BondDecoratorListScoreState.h (working copy) @@ -31,7 +31,7 @@ //! Find bonds amongst the following points. /** \param [in] ps The set of particles to use. */ - BondDecoratorListScoreState(const Particles &ps); + BondDecoratorListScoreState(const Particles &ps= Particles()); virtual ~BondDecoratorListScoreState(){} virtual void set_particles(const Particles &ps); Index: kernel/include/IMP/score_states/NonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/NonbondedListScoreState.h (revision 489) +++ kernel/include/IMP/score_states/NonbondedListScoreState.h (working copy) @@ -22,18 +22,19 @@ typedef std::vector BondedListScoreStates; //! A base class for classes that maintain a list of non-bonded pairs. -/** \ingroup restraint +/** */ class IMPDLLEXPORT NonbondedListScoreState: public ScoreState { private: + FloatKey rk_; + protected: // made protected for debugging code, do not use otherwise typedef std::vector > NBL; NBL nbl_; float last_cutoff_; - unsigned int size_nbl() const {return nbl_.size();} //! rebuild the nonbonded list @@ -46,14 +47,25 @@ nbl_.clear(); } + bool are_bonded(Particle *a, Particle *b) const { + for (BondedListConstIterator bli= bonded_lists_begin(); + bli != bonded_lists_end(); ++bli) { + if ((*bli)->are_bonded(a, b)) { + return true; + } + } + return false; + } + struct AddToNBL; friend struct AddToNBL; + // 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){} + p_(p){} void operator()(Particle *p) { operator()(p_, p); } @@ -62,16 +74,6 @@ } }; - bool are_bonded(Particle *a, Particle *b) const { - for (BondedListConstIterator bli= bonded_lists_begin(); - bli != bonded_lists_end(); ++bli) { - if ((*bli)->are_bonded(a, b)) { - return true; - } - } - return false; - } - //! tell the bonded lists what particles to pay attention to void propagate_particles(const Particles &ps); @@ -82,14 +84,29 @@ IMP_LOG(VERBOSE, "Found pair " << a->get_index() << " " << b->get_index() << std::endl); nbl_.push_back(std::make_pair(a, b)); + } else { + IMP_LOG(VERBOSE, "Pair " << a->get_index() + << " and " << b->get_index() << " rejected on bond" + <has_attribute(rk_)) { + return a->get_value(rk_); + } else { + return 0; + } + } + public: /** */ - NonbondedListScoreState(); + NonbondedListScoreState(FloatKey rk=FloatKey()); + FloatKey get_radius_key() const {return rk_;} + void set_radius_key(FloatKey rk) {rk_=rk;} + IMP_CONTAINER(BondedListScoreState, bonded_list, BondedListIndex); // kind of evil hack to make the names better Index: kernel/include/IMP/score_states/AllSphereNonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/AllSphereNonbondedListScoreState.h (revision 489) +++ kernel/include/IMP/score_states/AllSphereNonbondedListScoreState.h (working copy) @@ -1,76 +0,0 @@ -/** - * \file AllSphereNonbondedListScoreState.h - * \brief Allow iteration through pairs of a set of spheres. - * - * Copyright 2007-8 Sali Lab. All rights reserved. - */ - -#ifndef __IMP_ALL_SPHERE_NONBONDED_LIST_SCORE_STATE_H -#define __IMP_ALL_SPHERE_NONBONDED_LIST_SCORE_STATE_H - -#include "NonbondedListScoreState.h" -#include "../internal/ParticleGrid.h" -#include "../internal/kernel_version_info.h" -#include "MaxChangeScoreState.h" - -#include -#include - -namespace IMP -{ - -//! This class maintains a list of non-bonded pairs of spheres -/** \note The radius is currently assumed not to change. This could - be fixed later. - \note Points whose coordinates are not optimized are assumed to - stay that way (and vice - versa, although that direction doesn't matter so much). - \note 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). - - \ingroup restraint - */ -class IMPDLLEXPORT AllSphereNonbondedListScoreState: - public NonbondedListScoreState -{ -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 bins_; - std::vector fixed_bins_; - FloatKey rk_; - - //! \internal - void rebuild_nbl(Float cut); - - void repartition_points(const Particles &ps, std::vector &out); - - float side_from_r(float r) const {return r*1.6;} - - void generate_nbl(const Bin &particle_bin, const Bin &grid_bin, float cut); - - void cleanup(std::vector &bins); - -public: - /** - \param[in] ps A list of particles to use. - \param[in] radius The key to use to get the radius - */ - AllSphereNonbondedListScoreState(const Particles &ps, FloatKey radius); - ~AllSphereNonbondedListScoreState(); - IMP_SCORE_STATE(internal::kernel_version_info) - - void set_particles(const Particles &ps); -}; - -} // namespace IMP - -#endif /* __IMP_ALL_SPHERE_NONBONDED_LIST_SCORE_STATE_H */ Index: kernel/include/IMP/score_states/SConscript =================================================================== --- kernel/include/IMP/score_states/SConscript (revision 489) +++ kernel/include/IMP/score_states/SConscript (working copy) @@ -1,13 +1,15 @@ +Import('env') import os.path -Import('env') - -files = ['BondedListScoreState.h', 'NonbondedListScoreState.h', - 'BipartiteNonbondedListScoreState.h', 'MaxChangeScoreState.h', - 'BondDecoratorListScoreState.h', 'AllNonbondedListScoreState.h', - 'AllSphereNonbondedListScoreState.h'] - - -# Install the include files: -includedir = os.path.join(env['includedir'], 'IMP', 'score_states') +files=[] +files.append( 'AllNonbondedListScoreState.h' ) +files.append( 'BipartiteNonbondedListScoreState.h' ) +files.append( 'BondDecoratorListScoreState.h' ) +files.append( 'BondedListScoreState.h' ) +files.append( 'MaxChangeScoreState.h' ) +files.append( 'NonbondedListScoreState.h' ) +files.append( 'QuadraticAllNonbondedListScoreState.h' ) +files.append( 'QuadraticBipartiteNonbondedListScoreState.h' ) +files.append( 'QuadraticNonbondedListScoreState.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 0) +++ kernel/include/IMP/score_states/QuadraticBipartiteNonbondedListScoreState.h (revision 0) @@ -0,0 +1,51 @@ +/** + * \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 +#include + +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(const Particles &ps0, + const Particles &ps1); + QuadraticBipartiteNonbondedListScoreState(FloatKey rk=FloatKey()); + + 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 489) +++ kernel/include/IMP/score_states/BipartiteNonbondedListScoreState.h (working copy) @@ -8,39 +8,27 @@ #ifndef __IMP_BIPARTITE_NONBONDED_LIST_SCORE_STATE_H #define __IMP_BIPARTITE_NONBONDED_LIST_SCORE_STATE_H -#include "NonbondedListScoreState.h" -#include "../internal/ParticleGrid.h" -#include "../internal/kernel_version_info.h" +#include "QuadraticBipartiteNonbondedListScoreState.h" -#include -#include - 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. -class BondedListScoreState; - -//! 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. - \ingroup restraint + \ingroup restraint */ class IMPDLLEXPORT BipartiteNonbondedListScoreState: - public NonbondedListScoreState -{ - typedef NonbondedListScoreState P; - internal::ObjectPointer mc_; - internal::ParticleGrid grid_; - - virtual void rebuild_nbl(float cut); -public: + public QuadraticBipartiteNonbondedListScoreState { + typedef QuadraticBipartiteNonbondedListScoreState P; + public: + BipartiteNonbondedListScoreState(FloatKey rk, + const Particles &ps0, + const Particles &ps1): P(rk, ps0, ps1){} BipartiteNonbondedListScoreState(const Particles &ps0, - const Particles &ps1, - float target_side); - - IMP_SCORE_STATE(internal::kernel_version_info) - - void set_particles(const Particles &ps0, const Particles &ps1); + const Particles &ps1): P(ps0, ps1){} + BipartiteNonbondedListScoreState(FloatKey rk=FloatKey()): P(rk){} }; } // namespace IMP Index: kernel/include/IMP/score_states/QuadraticAllNonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/QuadraticAllNonbondedListScoreState.h (revision 0) +++ kernel/include/IMP/score_states/QuadraticAllNonbondedListScoreState.h (revision 0) @@ -0,0 +1,54 @@ +/** + * \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_SPHERE_NONBONDED_LIST_SCORE_STATE_H +#define __IMP_QUADRATIC_ALL_SPHERE_NONBONDED_LIST_SCORE_STATE_H + +#include "NonbondedListScoreState.h" +#include "../internal/kernel_version_info.h" +#include "QuadraticNonbondedListScoreState.h" + +#include +#include + +namespace IMP +{ + +//! This class maintains a list of non-bonded pairs of spheres +/** \note Points whose coordinates are not optimized are assumed to + stay that way and are pairs involving only fixed points are skipped + + \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; + Particles fixed_; + + //! \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= FloatKey(), + const Particles &ps= Particles()); + QuadraticAllNonbondedListScoreState(const Particles &ps); + ~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/MaxChangeScoreState.h =================================================================== --- kernel/include/IMP/score_states/MaxChangeScoreState.h (revision 489) +++ kernel/include/IMP/score_states/MaxChangeScoreState.h (working copy) @@ -18,20 +18,31 @@ { //! Keeps track of the maximum change of a set of attributes. +/** The score state maintains a list of particle and a list of + float attribute keys and keeps track of the maximum amount + any of these have changed since the last time reset was called. + + */ class IMPDLLEXPORT MaxChangeScoreState: public ScoreState { FloatKeys keys_; - std::vector > orig_; + std::vector orig_; Particles ps_; float max_change_; public: - MaxChangeScoreState(const FloatKeys &keys):keys_(keys){} + //! Track the changes with the specified keys. + MaxChangeScoreState(const FloatKeys &keys, + const Particles &ps= Particles()):keys_(keys){ + add_particles(ps); + } virtual ~MaxChangeScoreState(){} virtual void update(); + //! Measure differences from the current value. void reset(); + //! Return the maximum amount any attribute has changed. float get_max() const { return max_change_; } Index: kernel/include/IMP/score_states/AllNonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/AllNonbondedListScoreState.h (revision 489) +++ kernel/include/IMP/score_states/AllNonbondedListScoreState.h (working copy) @@ -1,6 +1,6 @@ /** * \file AllNonbondedListScoreState.h - * \brief Allow iteration through pairs of a set of particles. + * \brief Allow iteration through pairs of a set of s. * * Copyright 2007-8 Sali Lab. All rights reserved. */ @@ -9,7 +9,9 @@ #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 #include @@ -17,24 +19,59 @@ namespace IMP { -//! This class maintains a list of non-bonded pairs. -/** The nonbonded list returns all nonboned pairs in the set of points. +//! 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. + + \note The radius is currently assumed not to change. This could + be fixed later. + + \note Points whose coordinates are not optimized are assumed to + stay that way (and vice versa, although that direction doesn't + matter so much). + \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). + \ingroup restraint */ -class IMPDLLEXPORT AllNonbondedListScoreState: public NonbondedListScoreState +class IMPDLLEXPORT AllNonbondedListScoreState: + public NonbondedListScoreState { + typedef NonbondedListScoreState P; protected: - internal::ParticleGrid grid_; + //! \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 bins_; + std::vector fixed_bins_; - virtual void rebuild_nbl(float cut); + //! \internal + void rebuild_nbl(Float cut); + void repartition_points(const Particles &ps, std::vector &out); + + float side_from_r(float r) const; + + void generate_nbl(const Bin &particle_bin, const Bin &grid_bin, float cut); + + void cleanup(std::vector &bins); + public: /** \param[in] ps A list of particles to use. - \param[in] tvs A suggested size for the voxel side. - */ - AllNonbondedListScoreState(const Particles &ps, float tvs=1); - + \param[in] radius The key to use to get the radius + */ + AllNonbondedListScoreState(FloatKey radius= FloatKey(), + const Particles &ps=Particles()); + AllNonbondedListScoreState(const Particles &ps); + ~AllNonbondedListScoreState(); IMP_SCORE_STATE(internal::kernel_version_info) void set_particles(const Particles &ps); Index: kernel/test/states/test_nonbonded_list.py =================================================================== --- kernel/test/states/test_nonbonded_list.py (revision 489) +++ kernel/test/states/test_nonbonded_list.py (working copy) @@ -33,14 +33,17 @@ p= IMP.Particle() m.add_particle(p) d=IMP.XYZDecorator.create(p) + d.set_coordinates_are_optimized(True) d.randomize_in_box(IMP.Vector3D(0,0,0), IMP.Vector3D(10,10,10)); - s= IMP.AllNonbondedListScoreState(m.get_particles()) + s= IMP.AllNonbondedListScoreState(IMP.FloatKey(), + m.get_particles()) m.add_score_state(s) o= OnePair() r= IMP.NonbondedRestraint(o, s, 15) m.add_restraint(r) score= m.evaluate(False) + print score self.assertEqual(score, 4950, "Wrong score") def test_bl(self): @@ -51,6 +54,7 @@ p= IMP.Particle() m.add_particle(p) d=IMP.XYZDecorator.create(p) + d.set_coordinates_are_optimized(True) d.set_x(0) d.set_y(10) d.set_z(10) @@ -60,7 +64,8 @@ pts.append(p) for i in range(1,100): IMP.custom_bond(bds[i-1], bds[i], 1, .1) - s= IMP.AllNonbondedListScoreState(pts, 1) + s= IMP.AllNonbondedListScoreState(IMP.FloatKey(), + pts) b= IMP.BondDecoratorListScoreState(pts) s.add_bonded_list(b) m.add_score_state(s) @@ -84,19 +89,21 @@ m.add_particle(p) d=IMP.XYZDecorator.create(p) ps.append(p) + d.set_coordinates_are_optimized(True) if (i < 25): d.randomize_in_box(IMP.Vector3D(0,0,0), IMP.Vector3D(10,10,10)); else: d.randomize_in_box(IMP.Vector3D(60,60,60), IMP.Vector3D(70,70,70)); - s= IMP.AllNonbondedListScoreState(ps) + s= IMP.AllNonbondedListScoreState(IMP.FloatKey(), ps) m.add_score_state(s) o= OnePair() r= IMP.NonbondedRestraint(o, s, 15) m.add_restraint(r) score= m.evaluate(False) self.assertEqual(score, 1225, "Wrong score") + def test_bi(self): """Test the bipartite nonbonded list and restraint which uses it""" m= IMP.Model() @@ -106,13 +113,15 @@ p= IMP.Particle() m.add_particle(p) d=IMP.XYZDecorator.create(p) + d.set_coordinates_are_optimized(True) d.randomize_in_box(IMP.Vector3D(0,0,0), IMP.Vector3D(10,10,10)); if (i < 5): ps0.append(p) else: ps1.append(p) - s= IMP.BipartiteNonbondedListScoreState(ps0, ps1, 1000) + s= IMP.QuadraticBipartiteNonbondedListScoreState(IMP.FloatKey(), + ps0, ps1) m.add_score_state(s) o= OnePair() r= IMP.NonbondedRestraint(o, s, 15) @@ -120,7 +129,7 @@ score= m.evaluate(False) self.assertEqual(score, 25, "Wrong score") def test_spheres2(self): - """Test the nonbonded list of spheres (num pairs)""" + """Test the (quadratic and optimized) nonbonded lists of spheres (num pairs)""" m= IMP.Model() rk= IMP.FloatKey("radius") for i in range(0,100): @@ -131,12 +140,37 @@ IMP.Vector3D(10,10,10)); p.add_attribute(rk, i, False) d.set_coordinates_are_optimized(True) - s= IMP.AllSphereNonbondedListScoreState(m.get_particles(), rk) + s= IMP.AllNonbondedListScoreState(rk, + m.get_particles()) + s2= IMP.QuadraticAllNonbondedListScoreState(rk, m.get_particles()) m.add_score_state(s) + m.add_score_state(s2) o= OnePair() + o2= OnePair() r= IMP.NonbondedRestraint(o, s, 15) + r2= IMP.NonbondedRestraint(o2, s2, 15) m.add_restraint(r) + m.add_restraint(r2) score= m.evaluate(False) + self.assertEqual(score, 2*4950, "Wrong score") + def test_spheres3(self): + """Test the quadratic nonbonded list of spheres (num pairs)""" + m= IMP.Model() + rk= IMP.FloatKey("radius") + for i in range(0,100): + p= IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.randomize_in_box(IMP.Vector3D(0,0,0), + IMP.Vector3D(10,10,10)); + p.add_attribute(rk, i, False) + d.set_coordinates_are_optimized(True) + s= IMP.QuadraticAllNonbondedListScoreState(rk, m.get_particles()) + m.add_score_state(s) + o= OnePair() + r= IMP.NonbondedRestraint(o, s, 15) + m.add_restraint(r) + score= m.evaluate(False) self.assertEqual(score, 4950, "Wrong score") def test_frido_spheres(self): """Test the nonbonded list with frido's spheres""" @@ -262,7 +296,7 @@ print ".sphere "+str(d.get_x())+ " " + str(d.get_y())\ + " " + str(d.get_z()) + " " + str(p.get_value(rk)) - s= IMP.AllSphereNonbondedListScoreState(m.get_particles(), rk) + s= IMP.AllNonbondedListScoreState(rk, m.get_particles()) m.add_score_state(s) sd= IMP.SphereDistancePairScore(IMP.HarmonicLowerBound(0,1), rk) @@ -283,7 +317,7 @@ IMP.Vector3D(20,20,20)); p.add_attribute(rk, random.uniform(0,100), False) d.set_coordinates_are_optimized(True) - s= IMP.AllSphereNonbondedListScoreState(m.get_particles(), rk) + s= IMP.AllNonbondedListScoreState(rk, m.get_particles()) m.add_score_state(s) sd= IMP.SphereDistancePairScore(IMP.HarmonicLowerBound(0,1), rk) Index: kernel/src/score_states/BondDecoratorListScoreState.cpp =================================================================== --- kernel/src/score_states/BondDecoratorListScoreState.cpp (revision 489) +++ kernel/src/score_states/BondDecoratorListScoreState.cpp (working copy) @@ -19,7 +19,7 @@ void BondDecoratorListScoreState::update() { - IMP_LOG(VERBOSE, "Updating BondDecoratorList for " + IMP_LOG(TERSE, "Updating BondDecoratorList for " << ps_.size() << " particles" << std::endl); bonds_.clear(); for (unsigned int i=0; i< ps_.size(); ++i) { @@ -35,11 +35,12 @@ continue; } if (di < dj) { + IMP_LOG(VERBOSE, "Found bond " << di.get_bond(j) << std::endl); bonds_.push_back(di.get_bond(j)); } } } - IMP_LOG(VERBOSE, "Found " << bonds_.size() << " bonds"<< std::endl); + IMP_LOG(TERSE, "Found " << bonds_.size() << " bonds"<< std::endl); } void BondDecoratorListScoreState::set_particles(const Particles &ps) Index: kernel/src/score_states/NonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/NonbondedListScoreState.cpp (revision 489) +++ kernel/src/score_states/NonbondedListScoreState.cpp (working copy) @@ -13,7 +13,7 @@ namespace IMP { -NonbondedListScoreState::NonbondedListScoreState() +NonbondedListScoreState::NonbondedListScoreState(FloatKey rk): rk_(rk) { last_cutoff_=-1; } Index: kernel/src/score_states/AllSphereNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/AllSphereNonbondedListScoreState.cpp (revision 489) +++ kernel/src/score_states/AllSphereNonbondedListScoreState.cpp (working copy) @@ -1,242 +0,0 @@ -/** - * \file AllSphereNonbondedListScoreState.cpp - * \brief Allow iteration through pairs of a set of spheres. - * - * Copyright 2007-8 Sali Lab. All rights reserved. - */ - -#include "IMP/score_states/AllSphereNonbondedListScoreState.h" -#include "IMP/decorators/XYZDecorator.h" - -#include - -namespace IMP -{ - -static unsigned int min_grid_size=20; - -AllSphereNonbondedListScoreState -::AllSphereNonbondedListScoreState(const Particles &ps, FloatKey radius): - rk_(radius) -{ - set_particles(ps); -} - -AllSphereNonbondedListScoreState::~AllSphereNonbondedListScoreState() -{ - cleanup(bins_); - cleanup(fixed_bins_); -} - -void AllSphereNonbondedListScoreState::set_particles(const Particles &ps) -{ - NonbondedListScoreState::clear_nbl(); - // split into moving and fixed and call repartition twice - Particles moving, fixed; - for (unsigned int i=0; i< ps.size(); ++i) { - XYZDecorator d= XYZDecorator::cast(ps[i]); - if (d.get_coordinates_are_optimized()) { - moving.push_back(ps[i]); - } else { - fixed.push_back(ps[i]); - } - } - cleanup(bins_); - cleanup(fixed_bins_); - repartition_points(moving, bins_); - repartition_points(fixed, fixed_bins_); - for (unsigned int i=0; i< fixed_bins_.size(); ++i) { - fixed_bins_[i].grid->update(); - } -} - -void AllSphereNonbondedListScoreState::repartition_points(const Particles &ps, - std::vector &out) -{ - cleanup(out); - if (ps.empty()) return; - float minr=std::numeric_limits::max(), maxr=0; - for (unsigned int i=0; i< ps.size(); ++i) { - float r= ps[i]->get_value(rk_); - if ( r > maxr) maxr=r; - if ( r > 0 && r < minr) minr=r; - } - float curr=minr*2; - Floats cuts; - do { - cuts.push_back(curr); - curr *= 2; - } while (curr < maxr); - cuts.push_back(2*maxr); - - std::vector ops(cuts.size()); - for (unsigned int i=0; i< ps.size(); ++i) { - float r= ps[i]->get_value(rk_); - 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]); - break; - } - } - } - // consolidate - for (unsigned int i=1; i< ops.size(); ++i) { - if (ops[i-1].size() + ops[i].size() < min_grid_size) { - ops[i].insert(ops[i].end(), ops[i-1].begin(), ops[i-1].end()); - ops[i-1].clear(); - } - } - 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, ops[i][j]->get_value(rk_)); - } - 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]); - } - 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); - } - -#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 AllSphereNonbondedListScoreState::cleanup(std::vector &bins) -{ - for (unsigned int i=0; i< bins.size(); ++i) { - delete bins[i].grid; - } - bins.clear(); -} - -void AllSphereNonbondedListScoreState::update() -{ - NonbondedListScoreState::update(); - bool bad=false; - for (unsigned int i=0; i< bins_.size(); ++i) { - if (bins_[i].grid->update()) bad=true; - IMP_LOG(VERBOSE, bins_[i].rmax << std::endl << *bins_[i].grid << std::endl); - } - if (bad) { - IMP_LOG(VERBOSE, "Destroying nbl in Sphere list"<< std::endl); - NonbondedListScoreState::clear_nbl(); - } -} - -void AllSphereNonbondedListScoreState::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); - 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())); - 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); - } - -} - -void AllSphereNonbondedListScoreState::rebuild_nbl(Float cut) -{ - IMP_LOG(VERBOSE, "Rebuilding NBL with " << bins_.size() << " dynamic and " - << fixed_bins_.size() << " fixed" << std::endl); - 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); - } - - for (unsigned int j=0; j< fixed_bins_.size(); ++j) { - generate_nbl(bins_[i], fixed_bins_[j], cut); - } - - // same code as AllNonbonded. Would be nice to consolidate - internal::ParticleGrid::Index last_index; - for (internal::ParticleGrid::ParticleVoxelIterator it - = bins_[i].grid->particle_voxels_begin(); - it != bins_[i].grid->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, - true); - if (it->second != last_index) { - IMP_LOG(VERBOSE, "Searching in " << it->second - << std::endl); - bins_[i].grid->apply_to_cell_pairs(f, it->second); - } - last_index=it->second; - } - } - -#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]); - if (distance(di, dj) - ps[i]->get_value(rk_) - ps[j]->get_value(rk_) - <= 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; - } - } - IMP_assert(found, "Nonbonded list is missing " - << ps[i]->get_index() << " " << di - << " and " << ps[j]->get_index() << " " << dj << std::endl); - } - } - } -#endif -} - - -void AllSphereNonbondedListScoreState::show(std::ostream &out) const -{ - out << "AllSphereNonbondedListScoreState" << std::endl; -} - -} // namespace IMP Index: kernel/src/score_states/QuadraticBipartiteNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/QuadraticBipartiteNonbondedListScoreState.cpp (revision 0) +++ kernel/src/score_states/QuadraticBipartiteNonbondedListScoreState.cpp (revision 0) @@ -0,0 +1,74 @@ +/** + * \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(const Particles &ps0, + const Particles &ps1): P(FloatKey()) +{ + 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 489) +++ kernel/src/score_states/BipartiteNonbondedListScoreState.cpp (working copy) @@ -1,94 +0,0 @@ -/** - * \file BipartiteNonbondedListScoreState.cpp - * \brief Allow iteration through pairs of a set of atoms. - * - * Copyright 2007-8 Sali Lab. All rights reserved. - */ - -#include "IMP/score_states/BipartiteNonbondedListScoreState.h" -#include "IMP/decorators/XYZDecorator.h" -#include "IMP/score_states/MaxChangeScoreState.h" - -namespace IMP -{ - - /* - Iterator needs: - iterator into ps_, end iterator into ps_ - cell for ps_[i] - -> actually has value current iterator around cell - -put in logic for current - number of bins to search - - - nonbipartite: - iterator through all cells - index in current cell - -> actual value in iterator around cell - num bins - */ - -BipartiteNonbondedListScoreState -::BipartiteNonbondedListScoreState(const Particles &ps0, - const Particles &ps1, float target): - grid_(target) -{ - FloatKeys fks; - fks.push_back(FloatKey("x")); - fks.push_back(FloatKey("y")); - fks.push_back(FloatKey("z")); - mc_= new MaxChangeScoreState(fks); - set_particles(ps0, ps1); -} - -void BipartiteNonbondedListScoreState::rebuild_nbl(float cut) -{ - for (unsigned int i=0; i< mc_->get_particles().size(); ++i) { - Particle *p= mc_->get_particles()[i]; - NonbondedListScoreState::AddToNBL f(this, p); - XYZDecorator d= XYZDecorator::cast(p); - internal::ParticleGrid::VirtualIndex index - = grid_.get_virtual_index(Vector3D(d.get_x(), d.get_y(), d.get_z())); - IMP_LOG(VERBOSE, "Searching for " << p->get_index() - << " from cell " << index << std::endl); - grid_.apply_to_nearby(f, index, cut, false); - } - IMP_LOG(VERBOSE, "Found " << P::size_nbl() - << " bipartite nonbonded pairs" - << " for " << mc_->get_particles().size() - << " nongridded particles" << std::endl); -} - -void BipartiteNonbondedListScoreState::set_particles(const Particles &ps0, - const Particles &ps1) -{ - IMP_LOG(VERBOSE, "Setting bipartite particles " << ps0.size() - << " and " << ps1.size() << std::endl); - mc_->clear_particles(); - mc_->add_particles(ps1); - grid_.clear_particles(); - grid_.add_particles(ps0); - Particles all(ps0); - all.insert(all.end(), ps1.begin(), ps1.end()); - NonbondedListScoreState::propagate_particles(all); - /*IMP_assert(ps0.size() == P::get_particles().size(), - "Where did they go?");*/ - IMP_assert(ps1.size() == mc_->get_particles().size(), - "Where did we go?"); -} - -void BipartiteNonbondedListScoreState::update() -{ - NonbondedListScoreState::update(); - if (mc_->get_max() > grid_.get_voxel_size() - || grid_.update()) { - NonbondedListScoreState::clear_nbl(); - } -} - -void BipartiteNonbondedListScoreState::show(std::ostream &out) const -{ - out << "BipartiteNonbondedList" << std::endl; -} - -} // namespace IMP Index: kernel/src/score_states/SConscript =================================================================== --- kernel/src/score_states/SConscript (revision 489) +++ kernel/src/score_states/SConscript (working copy) @@ -1,9 +1,10 @@ Import('env') - -files = ['NonbondedListScoreState.cpp', 'BipartiteNonbondedListScoreState.cpp', - 'MaxChangeScoreState.cpp', 'BondDecoratorListScoreState.cpp', - 'AllNonbondedListScoreState.cpp', - 'AllSphereNonbondedListScoreState.cpp'] - -files = [File(x) for x in files] +files=[] +files.append(File( 'AllNonbondedListScoreState.cpp' )) +files.append(File( 'BondDecoratorListScoreState.cpp' )) +files.append(File( 'MaxChangeScoreState.cpp' )) +files.append(File( 'NonbondedListScoreState.cpp' )) +files.append(File( 'QuadraticAllNonbondedListScoreState.cpp' )) +files.append(File( 'QuadraticBipartiteNonbondedListScoreState.cpp' )) +files.append(File( 'QuadraticNonbondedListScoreState.cpp' )) Return('files') Index: kernel/src/score_states/QuadraticAllNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/QuadraticAllNonbondedListScoreState.cpp (revision 0) +++ kernel/src/score_states/QuadraticAllNonbondedListScoreState.cpp (revision 0) @@ -0,0 +1,77 @@ +/** + * \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 + +namespace IMP +{ + + +QuadraticAllNonbondedListScoreState +::QuadraticAllNonbondedListScoreState(FloatKey radius, + const Particles &ps): P(radius) +{ + set_particles(ps); +} + +QuadraticAllNonbondedListScoreState +::QuadraticAllNonbondedListScoreState(const Particles &ps): P(FloatKey()) +{ + set_particles(ps); +} + +QuadraticAllNonbondedListScoreState::~QuadraticAllNonbondedListScoreState() +{ +} + +void QuadraticAllNonbondedListScoreState::set_particles(const Particles &ps) +{ + Particles moving; + fixed_.clear(); + for (unsigned int i=0; i< ps.size(); ++i) { + XYZDecorator d= XYZDecorator::cast(ps[i]); + if (d.get_coordinates_are_optimized()) { + moving.push_back(ps[i]); + } else { + fixed_.push_back(ps[i]); + } + } + P::set_particles(moving); +} + + +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< fixed_.size(); ++i) { + P::handle_nbl_pair(fixed_[i], moving[j], cut); + } + 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/MaxChangeScoreState.cpp =================================================================== --- kernel/src/score_states/MaxChangeScoreState.cpp (revision 489) +++ kernel/src/score_states/MaxChangeScoreState.cpp (working copy) @@ -19,11 +19,16 @@ max_change_=0; for (unsigned int i=0; i < orig_.size(); ++i) { for (unsigned int j=0; j < keys_.size(); ++j) { + IMP_LOG(VERBOSE, "Particle " << ps_[i]->get_index() + << " and attribute " << keys_[j] + << " moved " << std::abs(ps_[i]->get_value(keys_[j]) + - orig_[i][j]) << std::endl); max_change_= std::max(max_change_, std::abs(ps_[i]->get_value(keys_[j]) - orig_[i][j])); } } + IMP_LOG(TERSE, "MaxChange update got " << max_change_ << std::endl); } Index: kernel/src/score_states/AllNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/AllNonbondedListScoreState.cpp (revision 489) +++ kernel/src/score_states/AllNonbondedListScoreState.cpp (working copy) @@ -1,54 +1,243 @@ /** * \file AllNonbondedListScoreState.cpp - * \brief Allow iteration through pairs of a set of atoms. + * \brief Allow iteration through pairs of a set of s. * * Copyright 2007-8 Sali Lab. All rights reserved. */ #include "IMP/score_states/AllNonbondedListScoreState.h" #include "IMP/decorators/XYZDecorator.h" -#include "IMP/internal/Grid3D.h" -#include "IMP/score_states/MaxChangeScoreState.h" +#include + namespace IMP { -AllNonbondedListScoreState::AllNonbondedListScoreState(const Particles &ps, - float tvs): - grid_(tvs) +static unsigned int min_grid_size=20; + +AllNonbondedListScoreState +::AllNonbondedListScoreState(FloatKey radius, const Particles &ps): + P(radius) { set_particles(ps); } -void AllNonbondedListScoreState::rebuild_nbl(float cut) +AllNonbondedListScoreState +::AllNonbondedListScoreState(const Particles &ps): + P(FloatKey()) { - IMP_LOG(VERBOSE, "Scanning for nonbonded pairs..." << std::flush); - internal::ParticleGrid::Index last_index; - for (internal::ParticleGrid::ParticleVoxelIterator it - = grid_.particle_voxels_begin(); - it != grid_.particle_voxels_end(); ++it) { - IMP_LOG(VERBOSE, "Searching with particle " << it->first->get_index() - << std::endl); - NonbondedListScoreState::AddToNBL f(this, it->first); - grid_.apply_to_nearby(f, it->second, cut, true); - if (it->second != last_index) { - grid_.apply_to_cell_pairs(f, it->second); + set_particles(ps); +} + + +AllNonbondedListScoreState::~AllNonbondedListScoreState() +{ + cleanup(bins_); + cleanup(fixed_bins_); +} + +float AllNonbondedListScoreState::side_from_r(float r) const { + if (r==0) return 1; + else return r*1.6; +} + + + +void AllNonbondedListScoreState::set_particles(const Particles &ps) +{ + NonbondedListScoreState::clear_nbl(); + // split into moving and fixed and call repartition twice + Particles moving, fixed; + for (unsigned int i=0; i< ps.size(); ++i) { + XYZDecorator d= XYZDecorator::cast(ps[i]); + if (d.get_coordinates_are_optimized()) { + moving.push_back(ps[i]); + } else { + fixed.push_back(ps[i]); } - last_index=it->second; } + cleanup(bins_); + cleanup(fixed_bins_); + repartition_points(moving, bins_); + repartition_points(fixed, fixed_bins_); + for (unsigned int i=0; i< fixed_bins_.size(); ++i) { + fixed_bins_[i].grid->update(); + } +} - IMP_LOG(VERBOSE, "done" << std::endl); - IMP_LOG(VERBOSE, "Found " << NonbondedListScoreState::size_nbl() - << " nonbonded pairs" - << std::endl); +void AllNonbondedListScoreState::repartition_points(const Particles &ps, + std::vector &out) +{ + cleanup(out); + if (ps.empty()) return; + float minr=std::numeric_limits::max(), maxr=0; + for (unsigned int i=0; i< ps.size(); ++i) { + float r= P::get_radius(ps[i]); + if ( r > maxr) maxr=r; + if ( r > 0 && r < minr) minr=r; + } + minr= std::min(maxr, minr); + float curr=minr*2; + Floats cuts; + cuts.push_back(0); + do { + cuts.push_back(curr); + curr *= 2; + } while (curr < maxr); + cuts.push_back(2*maxr); + std::vector ops(cuts.size()); + 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) { + IMP_assert(j< cuts.size(), "Internal error in ASNBLSS"); + if (cuts[j] >= r) { + ops[j].push_back(ps[i]); + found=true; + break; + } + } + IMP_assert(found, "Didn't put particle anywhere"); + } + // consolidate + for (unsigned int i=1; i< ops.size(); ++i) { + if (ops[i-1].size() + ops[i].size() < min_grid_size) { + ops[i].insert(ops[i].end(), ops[i-1].begin(), ops[i-1].end()); + ops[i-1].clear(); + } + } + 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]); + } + 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); + } + #ifndef NDEBUG - Particles ps= grid_.get_particles(); + 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 &bins) +{ + for (unsigned int i=0; i< bins.size(); ++i) { + delete bins[i].grid; + } + bins.clear(); +} + +void AllNonbondedListScoreState::update() +{ + NonbondedListScoreState::update(); + bool bad=false; + for (unsigned int i=0; i< bins_.size(); ++i) { + if (bins_[i].grid->update()) bad=true; + IMP_LOG(VERBOSE, bins_[i].rmax << std::endl << *bins_[i].grid << std::endl); + } + if (bad) { + IMP_LOG(VERBOSE, "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); + 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())); + 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); + } + +} + +void AllNonbondedListScoreState::rebuild_nbl(Float cut) +{ + IMP_LOG(TERSE, "Rebuilding NBL with " << bins_.size() << " dynamic and " + << fixed_bins_.size() << " fixed " + << " and cutoff " << cut << std::endl); + 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); + } + + for (unsigned int j=0; j< fixed_bins_.size(); ++j) { + generate_nbl(bins_[i], fixed_bins_[j], cut); + } + + // same code as AllNonbonded. Would be nice to consolidate + internal::ParticleGrid::Index last_index; + for (internal::ParticleGrid::ParticleVoxelIterator it + = bins_[i].grid->particle_voxels_begin(); + it != bins_[i].grid->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, + true); + if (it->second != last_index) { + IMP_LOG(VERBOSE, "Searching in " << it->second + << std::endl); + bins_[i].grid->apply_to_cell_pairs(f, it->second); + } + last_index=it->second; + } + } + + IMP_LOG(TERSE, "NBL has " << P::nbl_.size() << " 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]); - if (distance(di, dj) <= cut && !are_bonded(ps[i], ps[j]) ) { + 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) { @@ -67,25 +256,10 @@ #endif } -void AllNonbondedListScoreState::set_particles(const Particles &ps) -{ - grid_.clear_particles(); - grid_.add_particles(ps); - NonbondedListScoreState::propagate_particles(ps); -} -void AllNonbondedListScoreState::update() -{ - NonbondedListScoreState::update(); - if (grid_.update()) { - NonbondedListScoreState::clear_nbl(); - } - IMP_LOG(VERBOSE, grid_); -} - void AllNonbondedListScoreState::show(std::ostream &out) const { - out << "NonbondedList" << std::endl; + out << "AllNonbondedListScoreState" << std::endl; } } // namespace IMP Index: kernel/src/score_states/QuadraticNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/QuadraticNonbondedListScoreState.cpp (revision 0) +++ kernel/src/score_states/QuadraticNonbondedListScoreState.cpp (revision 0) @@ -0,0 +1,55 @@ +/** + * \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 + +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