March 2008 - IMP-dev - salilab.org

GFit
by Daniel Russel 24 Mar '08

24 Mar '08

There was a poster at the TCNP meeting from the developers of gfit. It is a very general package for building models to fit experimental data via optimization. You basically write a matlab program for your model and describe a set of experiments (as values for some of the input and output parameters for your model) and gfit tries to optimize the unknowns in your model to fid the experimental data. http://gfit.sourceforge.net It is from the same group as wrote Virtual Cell. Something we should be aware of.

2 1

A few small patches
by Daniel Russel 17 Mar '08

17 Mar '08

First, ObjectPointer checks to in operator= to make sure another ObjectPointer doesn't own the object already. BondDecoratorRestraint was catching too many exceptions which hid some real errors. Now it only catches exceptions caused by the internal structure of the bonds being screwed up. exception.h was missing some includes And Vector.h was missing an include of exception.h

1 0

Remove Vector3D::get_component()
by Ben Webb 17 Mar '08

17 Mar '08

The attached patch removes Vector3D::get_component(). We don't need it any more, since we have operator[] now which does the same thing. Ben -- ben(a)salilab.org http://salilab.org/~ben/ "It is a capital mistake to theorize before one has data." - Sir Arthur Conan Doyle

1 0

Logging
by Daniel Russel 15 Mar '08

15 Mar '08

Logging is currently a bit of a mess in IMP as it basically output whatever I was last interested in looking at. One way of systematizing it would be VERBOSE: each restraint/state should output enough information to determine the flow through the code. For example, a restraint should print out each term of its energy (so the energy for each pair of particles for example). TERSE: each restraint/state should print a constant number of lines per evaluate call

2 1

Removing deprecated restraints
by Daniel Russel 12 Mar '08

12 Mar '08

Since Marc asked me about using one of the restraints that was scheduled for destruction, I figured I should just go ahead and remove them all. The patch removes PairConnectivety, Torus, and Proximity restraints and their test code and the xml reading code which we already had to warn him away from.

2 1

Minor patches
by Daniel Russel 07 Mar '08

07 Mar '08

First, I changed vector so that it is get_magnitude rather than just magnitude to be consistent with the rest of IMP. I also added operator* for the dot product and operator- (and difference for python) since that is pretty standard and added a method to XYZDecorator to get a vector for the current position. Second, I added a check for NaN derivative values in Particle. Third, I added a comment about the internal checks in nonbonded list test. Index: kernel/src/restraints/DihedralRestraint.cpp =================================================================== --- kernel/src/restraints/DihedralRestraint.cpp (revision 457) +++ kernel/src/restraints/DihedralRestraint.cpp (working copy) @@ -58,7 +58,7 @@ Vector3D v1 = rij.vector_product(rkj); Vector3D v2 = rkj.vector_product(rkl); Float scalar_product = v1.scalar_product(v2); - Float mag_product = v1.magnitude() * v2.magnitude(); + Float mag_product = v1.get_magnitude() * v2.get_magnitude(); // avoid division by zero Float cosangle = std::abs(mag_product) > 1e-12 @@ -86,9 +86,9 @@ // J. Mol. Biol. 234, 751-762 (1993) Vector3D vijkj = rij.vector_product(rkj); Vector3D vkjkl = rkj.vector_product(rkl); - Float sijkj2 = vijkj.squared_magnitude(); - Float skjkl2 = vkjkl.squared_magnitude(); - Float skj = rkj.magnitude(); + Float sijkj2 = vijkj.get_squared_magnitude(); + Float skjkl2 = vkjkl.get_squared_magnitude(); + Float skj = rkj.get_magnitude(); Float rijkj = rij.scalar_product(rkj); Float rkjkl = rkj.scalar_product(rkl); Index: kernel/include/IMP/decorators/XYZDecorator.h =================================================================== --- kernel/include/IMP/decorators/XYZDecorator.h (revision 457) +++ kernel/include/IMP/decorators/XYZDecorator.h (working copy) @@ -82,6 +82,13 @@ b.get_coordinate(2) - get_coordinate(2)); } + //! Convert it to a vector. + /** Somewhat suspect based on wanting a Point/Vector differentiation + but we don't have points*/ + Vector3D get_vector() const { + return Vector3D(get_x(), get_y(), get_z()); + } + //! Get a vector containing the keys for x,y,z /** This is quite handy for initializing movers and things. */ Index: kernel/include/IMP/Vector3D.h =================================================================== --- kernel/include/IMP/Vector3D.h (revision 457) +++ kernel/include/IMP/Vector3D.h (working copy) @@ -52,6 +52,11 @@ return vec_[0] * vec2.vec_[0] + vec_[1] * vec2.vec_[1] + vec_[2] * vec2.vec_[2]; } + + //! scalar product + Float operator*(const Vector3D &o) const { + return scalar_product(o); + } //! \return the vector product of two vectors. /** \param[in] vec2 The other vector to use in the product. @@ -63,23 +68,33 @@ } //! \return The square of the magnitude of this vector. - Float squared_magnitude() const { + Float get_squared_magnitude() const { return vec_[0] * vec_[0] + vec_[1] * vec_[1] + vec_[2] * vec_[2]; } //! \return The magnitude of this vector. - Float magnitude() const { - return std::sqrt(squared_magnitude()); + Float get_magnitude() const { + return std::sqrt(get_squared_magnitude()); } //! \return This vector normalized to unit length. Vector3D get_unit_vector() const { - Float mag = magnitude(); + Float mag = get_magnitude(); // avoid division by zero mag = std::max(mag, static_cast<Float>(1e-12)); return Vector3D(vec_[0] / mag, vec_[1] / mag, vec_[2] / mag); } + //! for C++ + Vector3D operator-(const Vector3D &o) const { + return Vector3D(operator[](0)-o[0], + operator[](1)-o[1], + operator[](2)-o[2]); + } + //! for python + Vector3D difference(const Vector3D &o) const { + return operator-(o); + } private: Float vec_[3]; }; Index: kernel/include/IMP/Particle.h =================================================================== --- kernel/include/IMP/Particle.h (revision 457) +++ kernel/include/IMP/Particle.h (working copy) @@ -336,6 +336,7 @@ { IMP_check(get_is_active(), "Do not touch inactive particles", InactiveParticleException()); + IMP_assert(value==value, "Can't add NaN to derivative"); float_indexes_.get_value(name).derivative+= da(value); } Index: kernel/test/states/test_nonbonded_list.py =================================================================== --- kernel/test/states/test_nonbonded_list.py (revision 457) +++ kernel/test/states/test_nonbonded_list.py (working copy) @@ -119,6 +119,8 @@ self.assertEqual(score, 25, "Wrong score") def test_spheres2(self): """Test the nonbonded list of spheres (num pairs)""" + # This uses the internal checks of the nonbonded list to + # verify correctness m= IMP.Model() rk= IMP.FloatKey("radius") for i in range(0,100):

2 2

Patch to AllSpheresNonbonded
by Daniel Russel 05 Mar '08

05 Mar '08

Here are some changes to all spheres nonbonded list to fix a bug that Frido found as well as make things faster. The nonbonded lists now have checks built into them so asserts should fail if anything goes wrong. Note that these checks are very expensive for large numbers of particles (quadratic) and so any real runs with IMP should use builds with NDEBUG defined. We may want to define macros like IMP_CHECK_EXPENSIVE to provide finer grain control over debugging checks. Index: kernel/include/IMP/score_states/NonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/NonbondedListScoreState.h (revision 456) +++ kernel/include/IMP/score_states/NonbondedListScoreState.h (working copy) @@ -27,11 +27,13 @@ class IMPDLLEXPORT NonbondedListScoreState: public ScoreState { private: +protected: + // made protected for debuggin code, do not use otherwise typedef std::vector<std::pair<Particle*, Particle*> > NBL; NBL nbl_; float last_cutoff_; -protected: + unsigned int size_nbl() const {return nbl_.size();} //! rebuild the nonbonded list @@ -60,20 +62,23 @@ } }; + 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); void add_to_nbl(Particle *a, Particle *b) { if (!a->get_is_active() || !b->get_is_active()) return; - bool found=false; - for (BondedListIterator bli= bonded_lists_begin(); - bli != bonded_lists_end(); ++bli) { - if ((*bli)->are_bonded(a, b)) { - found = true; - break; - } - } - if (!found) { + + if (!are_bonded(a,b)) { IMP_LOG(VERBOSE, "Found pair " << a->get_index() << " " << b->get_index() << std::endl); nbl_.push_back(std::make_pair(a, b)); @@ -90,6 +95,7 @@ // kind of evil hack to make the names better // perhaps the macro should be made more flexible typedef BondedListScoreStateIterator BondedListIterator; + typedef BondedListScoreStateConstIterator BondedListConstIterator; IMP_SCORE_STATE(internal::kernel_version_info) Index: kernel/test/states/test_nonbonded_list.py =================================================================== --- kernel/test/states/test_nonbonded_list.py (revision 456) +++ kernel/test/states/test_nonbonded_list.py (working copy) @@ -27,7 +27,7 @@ def test_it(self): """Test the nonbonded list and restraint which uses it""" m= IMP.Model() - for i in range(0,10): + for i in range(0,100): p= IMP.Particle() m.add_particle(p) d=IMP.XYZDecorator.create(p) @@ -39,13 +39,13 @@ r= IMP.NonbondedRestraint(o, s, 15) m.add_restraint(r) score= m.evaluate(False) - self.assertEqual(score, 45, "Wrong score") + self.assertEqual(score, 4950, "Wrong score") def test_bl(self): """Test the bonded list""" m= IMP.Model() bds=[] - for i in range(0,10): + for i in range(0,100): p= IMP.Particle() m.add_particle(p) d=IMP.XYZDecorator.create(p) @@ -56,7 +56,7 @@ pts=IMP.Particles() for p in m.get_particles(): pts.append(p) - for i in range(1,10): + for i in range(1,100): IMP.custom_bond(bds[i-1], bds[i], 1, .1) s= IMP.AllNonbondedListScoreState(pts, 1) b= IMP.BondDecoratorListScoreState(pts) @@ -70,7 +70,8 @@ m.add_restraint(r) m.add_restraint(br) score= m.evaluate( False ) - self.assertEqual(score, 900+45-9, "Wrong score") + print "Score with bonds is " + str(score) + self.assertEqual(score, 9900+4950-99, "Wrong score") def test_distfilt(self): """Test filtering based on distance in nonbonded list""" @@ -120,7 +121,7 @@ """Test the nonbonded list of spheres (num pairs)""" m= IMP.Model() rk= IMP.FloatKey("radius") - for i in range(0,10): + for i in range(0,100): p= IMP.Particle() m.add_particle(p) d=IMP.XYZDecorator.create(p) @@ -134,18 +135,149 @@ r= IMP.NonbondedRestraint(o, s, 15) m.add_restraint(r) score= m.evaluate(False) - self.assertEqual(score, 45, "Wrong score") + self.assertEqual(score, 4950, "Wrong score") + def test_frido_spheres(self): + """Test the nonbonded list with frido's spheres""" + m= IMP.Model() + rk= IMP.FloatKey("radius") + print "Frido begin" + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(87.0621490479) + d.set_y(140.299957275) + d.set_z(76.7119979858) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 20.1244087219, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(80.805770874) + d.set_y(99.9667434692) + d.set_z(66.9167098999) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 19.8160457611, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(112.992515564) + d.set_y(119.602111816) + d.set_z(53.6557235718) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 20.3428726196, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(118.784294128) + d.set_y(86.842376709) + d.set_z(65.2264938354) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 20.1794681549, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(132.930664062) + d.set_y(85.9250793457) + d.set_z(62.0034713745) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 19.730260849, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(97.6803894043) + d.set_y(77.0734939575) + d.set_z(45.9188423157) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 20.0133743286, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(74.9787445068) + d.set_y(103.48789978) + d.set_z(65.8464813232) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 20.1519756317, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(74.5077972412) + d.set_y(114.997116089) + d.set_z(103.185188293) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 18.6368904114, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(110.806472778) + d.set_y(100.937240601) + d.set_z(91.8201828003) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 19.2294254303, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(125.188072205) + d.set_y(135.360610962) + d.set_z(84.0617752075) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 18.9221935272, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(122.894897461) + d.set_y(114.53062439) + d.set_z(92.5285186768) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 18.9533672333, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(124.656105042) + d.set_y(88.4534988403) + d.set_z(99.5341949463) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 18.8280506134, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(109.428604126) + d.set_y(60.1015129089) + d.set_z(79.1919250488) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 19.1991424561, False) + p=IMP.Particle() + m.add_particle(p) + d=IMP.XYZDecorator.create(p) + d.set_x(79.7605895996) + d.set_y(78.1874771118) + d.set_z(95.7365951538) + d.set_coordinates_are_optimized(True) + p.add_attribute(rk, 19.3197059631, False) + for p in m.get_particles(): + d= IMP.XYZDecorator.cast(p) + 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) + m.add_score_state(s) + sd= IMP.SphereDistancePairScore(IMP.HarmonicLowerBound(0,1), + rk) + r= IMP.NonbondedRestraint(sd, s, 1) + m.add_restraint(r) + score= m.evaluate(False) + print "Frido score is" + print score def test_spheres(self): """Test the nonbonded list of spheres (collision detection)""" m= IMP.Model() rk= IMP.FloatKey("radius") - for i in range(0,10): + 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, random.uniform(0,1000), False) + 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) m.add_score_state(s) @@ -156,6 +288,7 @@ score= m.evaluate(False) opt= IMP.ConjugateGradients() opt.set_model(m) + IMP.set_log_level(IMP.TERSE) score =opt.optimize(10000) print score for p in m.get_particles(): Index: kernel/include/IMP/internal/ParticleGrid.h =================================================================== --- kernel/include/IMP/internal/ParticleGrid.h (revision 456) +++ kernel/include/IMP/internal/ParticleGrid.h (working copy) @@ -42,6 +42,8 @@ const Particles& get_particles() const {return mc_->get_particles();} bool update(); + void show(std::ostream &out) const; + typedef Grid::VirtualIndex VirtualIndex; typedef Grid::Index Index; Grid::VirtualIndex get_virtual_index(Vector3D pt) const { @@ -148,6 +150,8 @@ if (curp_== grid_->get_voxel(*cvoxel_).end()) { ++cvoxel_; find_voxel(); + } else { + temp_= std::make_pair(*curp_, *cvoxel_); } } @@ -178,6 +182,8 @@ } }; +IMP_OUTPUT_OPERATOR(ParticleGrid); + } // namespace internal } // namespace IMP Index: kernel/src/score_states/AllSphereNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/AllSphereNonbondedListScoreState.cpp (revision 456) +++ kernel/src/score_states/AllSphereNonbondedListScoreState.cpp (working copy) @@ -13,6 +13,8 @@ namespace IMP { +static unsigned int min_grid_size=20; + AllSphereNonbondedListScoreState ::AllSphereNonbondedListScoreState(const Particles &ps, FloatKey radius): rk_(radius) @@ -59,29 +61,40 @@ if ( r > maxr) maxr=r; if ( r > 0 && r < minr) minr=r; } - float curr=minr; + float curr=minr*2; Floats cuts; do { cuts.push_back(curr); curr *= 2; } while (curr < maxr); - cuts.push_back(maxr); + cuts.push_back(2*maxr); std::vector<Particles> 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) { - if (cuts[j] >= r) { + 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()); - out.back().rmax= cuts[i]; + 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; @@ -90,7 +103,7 @@ 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->get_particles().size() << std::endl); + << ": " << *out[i].grid << std::endl); } #ifndef NDEBUG @@ -121,6 +134,7 @@ 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); @@ -186,6 +200,37 @@ 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 } Index: kernel/src/score_states/AllNonbondedListScoreState.cpp =================================================================== --- kernel/src/score_states/AllNonbondedListScoreState.cpp (revision 456) +++ kernel/src/score_states/AllNonbondedListScoreState.cpp (working copy) @@ -41,6 +41,30 @@ IMP_LOG(VERBOSE, "Found " << NonbondedListScoreState::size_nbl() << " nonbonded pairs" << std::endl); + +#ifndef NDEBUG + Particles ps= grid_.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]); + if (distance(di, dj) <= 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 AllNonbondedListScoreState::set_particles(const Particles &ps) @@ -56,6 +80,7 @@ if (grid_.update()) { NonbondedListScoreState::clear_nbl(); } + IMP_LOG(VERBOSE, grid_); } void AllNonbondedListScoreState::show(std::ostream &out) const Index: kernel/src/internal/ParticleGrid.cpp =================================================================== --- kernel/src/internal/ParticleGrid.cpp (revision 456) +++ kernel/src/internal/ParticleGrid.cpp (working copy) @@ -108,6 +108,19 @@ } } + +void ParticleGrid::show(std::ostream &out) const { + for (Grid::IndexIterator it= grid_.all_indexes_begin(); + it != grid_.all_indexes_end(); ++it) { + out << *it << ": "; + //Grid::Index + for (unsigned int i=0; i< grid_.get_voxel(*it).size(); ++i) { + out << grid_.get_voxel(*it)[i]->get_index() << " "; + } + out << std::endl; + } +} + } // namespace internal } // namespace IMP

2 1

Remove uniform_01 from Monte Carlo code
by Ben Webb 05 Mar '08

05 Mar '08

The attached patch removes the code from the Monte Carlo NormalMover::generate_move method that uses boost::uniform_01, and replaces it with code that uses boost::variate_generator instead. This is because uniform_01 takes the random number generator by value rather than reference; thus, it creates a temporary copy of it and doesn't update the generator state (so that every time you call the method, you get the same result, unless you call some other method somewhere which _does_ update the RNG). At least, this is what was happening in the MD code, which also used uniform_01, until I fixed it in r454. variate_generator in the attached patch takes a reference to the underlying RNG, so *does* update the state. variate_generator seems to be the "official" Boost way to tie a distribution with a random number generator these days anyway. Comments? Ben -- ben(a)salilab.org http://salilab.org/~ben/ "It is a capital mistake to theorize before one has data." - Sir Arthur Conan Doyle Index: kernel/src/optimizers/movers/NormalMover.cpp =================================================================== --- kernel/src/optimizers/movers/NormalMover.cpp (revision 452) +++ kernel/src/optimizers/movers/NormalMover.cpp (working copy) @@ -26,13 +26,15 @@ void NormalMover::generate_move(float scale) { - std::vector<Float> center(number_of_float_keys()); - boost::uniform_01<RandomNumberGenerator> u01(random_number_generator); boost::normal_distribution<double> mrng(0, stddev_); + boost::variate_generator<RandomNumberGenerator&, + boost::normal_distribution<double> > + sampler(random_number_generator, mrng); + for (unsigned int i = 0; i < number_of_particles(); ++i) { for (unsigned int j = 0; j < number_of_float_keys(); ++j) { float c = get_float(i, j); - float r = mrng(u01); + float r = sampler(); // Check for NaN (x!=x when x==NaN) (can only use std::isnan with C99) IMP_assert(r == r, "Bad random"); IMP_assert(c == c, "Bad stored");

1 0

Change constructor argument orders
by Daniel Russel 04 Mar '08

04 Mar '08

As I threatened, I changed the argument orders for the nice restraints and for the scoring functions so that the score function or unary function is always first. Notes about these conventions have been put in various base classes. XYZDecorator.h should be committed separately as the only change to there is to make the get_xyz_keys function static rather than const.

1 0

Proposed rearrangement
by Daniel Russel 04 Mar '08

04 Mar '08

I find it currently somewhat confusing what the order of arguments is to restraints and score function constructors. I propose establishing the convention that the unary function or the score function is always first. Such a convention has several advantages: - it makes it easier to get the thisown line in IMP.i right - the remainder of the arguments are less consistent, but everything takes and owns a score function or a unary function - the score function is needed for the restraints, other things are less important - we can then make the other arguments option. For example, I generally create restraints before I know the particles they will act on (since those change), so I create them with an empty list of particles. It would be a minor convenience to be able to omit this I'll make the changes if no one objects. It should take only a few minutes.

2 1