Index: kernel/test/cg_optimizer/test_cg.py =================================================================== --- kernel/test/cg_optimizer/test_cg.py (revision 352) +++ kernel/test/cg_optimizer/test_cg.py (working copy) @@ -1,64 +0,0 @@ -import unittest -import IMP, IMP.test - -class WoodsFunc(IMP.Restraint): - """Woods function for four input values, defined as an IMP restraint""" - def __init__(self, model, particles): - IMP.Restraint.__init__(self) - self.particles= particles - self.index= IMP.FloatKey("x") - def show(self, junk): - print "Woods function" - def last_modified_by(self): - return "Daniel Russel" - - - def evaluate(self, accum): - (x1, x2, x3, x4) = [p.get_value(self.index) for p in self.particles] - a = x2 - x1 * x1 - b = x4 - x3 * x3 - e = 100.0 * a * a + (1.0 - x1) ** 2+ 90.0 * b * b + (1.0 - x3) ** 2 \ - + 10.1 * ((x2 - 1.0) ** 2 + (x4 - 1.0) ** 2) \ - + 19.8 * (x2 - 1.0) * (x4 - 1.0) - if accum: - dx = [-2.0 * (200.0 * x1 * a + 1.0 - x1), - 2.0 * (100.0 * a + 10.1 * (x2 - 1.0) + 9.9 * (x4 - 1.0)), - -2.0 * (180.0 * x3 * b + 1.0 - x3), - 2.0 * (90.0 * b + 10.1 * (x4 - 1.0) + 9.9 * (x2 - 1.0))] - for (p,d) in zip(self.particles, dx): - p.add_to_derivative(self.index, d, accum); - #for (i, d) in zip(self.indices, dx): - # accum.add_to_deriv(i, d) - return e - -class CGOptimizerTests(IMP.test.TestCase): - def test_cg_woods_func(self): - """Check that we can optimize the Woods function with CG""" - self._test_starting_conditions((-3.0, -1.0, -3.0, -1.0)) - self._test_starting_conditions((2.0, 3.0, 8.0, -5.0)) - - def _test_starting_conditions(self, starting_values): - """Test the optimizer with given starting conditions""" - model = IMP.Model() - particles = [] - - for value in starting_values: - p = IMP.Particle() - model.add_particle(p) - particles.append(p) - p.add_attribute(IMP.FloatKey("x"), value, True) - rsr = WoodsFunc(model, particles) - model.add_restraint(rsr) - model.set_up_trajectory('', False, False) - opt = IMP.ConjugateGradients() - opt.set_model(model) - opt.set_threshold(1e-5) - e = opt.optimize(100) - model_data = model.get_model_data() - for p in particles: - val = p.get_value(IMP.FloatKey("x")) - self.assertAlmostEqual(val, 1.0, places=1) - self.assertAlmostEqual(e, 0.0, places=2) - -if __name__ == '__main__': - unittest.main() Index: kernel/test/md_optimizer/test_md_optimizer.py =================================================================== --- kernel/test/md_optimizer/test_md_optimizer.py (revision 352) +++ kernel/test/md_optimizer/test_md_optimizer.py (working copy) @@ -1,118 +0,0 @@ -import unittest -import IMP.utils -import IMP.test, IMP - -xkey = IMP.FloatKey('x') -ykey = IMP.FloatKey('y') -zkey = IMP.FloatKey('z') -vxkey = IMP.FloatKey('vx') - -# Conversion from derivatives (in kcal/mol/A) to acceleration (A/fs/fs) -kcal2mod = 4.1868e-7 -# Mass of Carbon-12 (kg/mol) -cmass = 0.012011 - -class XTransRestraint(IMP.Restraint): - """Attempt to move the whole system along the x axis""" - def __init__(self, strength): - IMP.Restraint.__init__(self) - self.strength = strength - - def evaluate(self, accum): - e = 0. - for p in self.get_model().get_particles(): - e += p.get_value(xkey) * self.strength - if accum: - for p in self.get_model().get_particles(): - p.add_to_derivative(xkey, self.strength, accum) - p.add_to_derivative(ykey, 0.0, accum) - p.add_to_derivative(zkey, 0.0, accum) - return e - def version(self): - return "0.5" - def last_modified_by(self): - return "Daniel Russel" - - -class WriteTrajState(IMP.OptimizerState): - """Write system coordinates (trajectory) into a Python list""" - def __init__(self, traj): - IMP.OptimizerState.__init__(self) - self.traj = traj - def update(self): - model = self.get_optimizer().get_model() - self.traj.append([(p.get_value(xkey), p.get_value(ykey), - p.get_value(zkey), p.get_value(vxkey)) \ - for p in model.get_particles()]) - - -class MolecularDynamicsTests(IMP.test.TestCase): - """Test molecular dynamics optimizer""" - - def setUp(self): - """Set up particles and optimizer""" - - self.model = IMP.Model() - self.particles = [] - self.particles.append(IMP.utils.XYZParticle(self.model, - -43.0, 65.0, 93.0)) - self.md = IMP.MolecularDynamics() - self.md.set_model(self.model) - self.md.set_temperature(300) - - def _check_trajectory(self, coor, traj, timestep, vxfunc): - """Check generated trajectory against that predicted using vxfunc""" - vx = 0. - msg = "Predicted coordinate %.5f doesn't match generated %.5f, " + \ - "for step %d, coordinate %d[%d]" - velmsg = "Predicted velocity %.5f doesn't match generated %.5f, " + \ - "for step %d, particle %d" - for (num, step) in enumerate(traj): - for n in range(len(coor)): - coor[n][0] += vx * timestep - self.assert_(abs(vx - step[n][3]) <= 1e-3, - velmsg % (vx, step[n][3], num, n)) - for d in range(3): - self.assert_(abs(coor[n][d] - step[n][d]) <= 1e-3, - msg % (coor[n][d], step[n][d], num, n, d)) - vx = vxfunc(vx) - - def _optimize_model(self, timestep): - """Run a short MD optimization on the model.""" - traj = [] - start = [[p.get_value(xkey), p.get_value(ykey), p.get_value(zkey)] \ - for p in self.model.get_particles()] - state = WriteTrajState(traj) - self.md.add_optimizer_state(state) - self.md.set_time_step(timestep) - self.md.optimize(50) - return start, traj - - def test_nonrigid_translation(self): - """Check that non-rigid MD translation is Newtonian""" - timestep = 4.0 - strength = 50.0 - r = XTransRestraint(strength) - self.model.add_restraint(r) - (start, traj) = self._optimize_model(timestep) - delttm = -timestep * kcal2mod / cmass - self._check_trajectory(start, traj, timestep, - lambda a: a + strength * delttm) - - def test_non_xyz(self): - """Should be unable to do MD on optimizable non-xyz attributes""" - p = IMP.Particle() - self.model.add_particle(p) - p.add_attribute(IMP.FloatKey("attr"), 0.0, True) - self.assertRaises(ValueError, self.md.optimize, 50) - - def test_make_velocities(self): - """Test that MD generates particle velocities""" - self.md.optimize(0) - keys = [IMP.FloatKey(x) for x in ("vx", "vy", "vz")] - for p in self.particles: - for key in keys: - self.assert_(p.has_attribute(key)) - -if __name__ == '__main__': - unittest.main() Index: kernel/test/optimizers/test_mc.py =================================================================== --- kernel/test/optimizers/test_mc.py (revision 0) +++ kernel/test/optimizers/test_mc.py (revision 0) @@ -0,0 +1,101 @@ +import unittest +import IMP, IMP.test + +class WoodsFunc(IMP.Restraint): + """Woods function for four input values, defined as an IMP restraint""" + def __init__(self): + IMP.Restraint.__init__(self) + self.index= IMP.FloatKey("x") + def evaluate(self, accum): + #print "Evaluating in python\n" + (x1, x2, x3, x4) = [p.get_value(self.index) for p in self.get_model().get_particles()] + #print "Evaluating at " +str(x1) + " " + str(x2) + " " \ + # + str(x3) + " " + str(x4) + a = x2 - x1 * x1 + b = x4 - x3 * x3 + e = 100.0 * a * a + (1.0 - x1) ** 2+ 90.0 * b * b + (1.0 - x3) ** 2 \ + + 10.1 * ((x2 - 1.0) ** 2 + (x4 - 1.0) ** 2) \ + + 19.8 * (x2 - 1.0) * (x4 - 1.0) + if accum: + dx = [-2.0 * (200.0 * x1 * a + 1.0 - x1), + 2.0 * (100.0 * a + 10.1 * (x2 - 1.0) + 9.9 * (x4 - 1.0)), + -2.0 * (180.0 * x3 * b + 1.0 - x3), + 2.0 * (90.0 * b + 10.1 * (x4 - 1.0) + 9.9 * (x2 - 1.0))] + for (p,d) in zip(self.get_model().get_particles(), dx): + p.add_to_derivative(self.index, d, accum); + return e + def version(self): + return "0.5"; + def last_modified_by(self): + return "Daniel Russel" + def show(self, fakestream): + print "WoodsFunc" + +class MCOptimizerTest(IMP.test.TestCase): + def setUp(self): + IMP.set_log_level(IMP.TERSE) + self.xkey= IMP.FloatKey("x") + self.rsrs=[] + + def test_c1(self): + """test montecarlo 1""" + (model, opt)= self._setup_opt() + lopt= IMP.ConjugateGradients() + opt.set_local_optimizer(lopt) + self._test_starting_conditions(model, opt, (-3.0, -1.0, -3.0, -1.0), 5) + def test_c2(self): + """test montecarlo 2""" + (model, opt)= self._setup_opt() + lopt= IMP.ConjugateGradients() + opt.set_local_optimizer(lopt) + self._test_starting_conditions(model, opt, (2.0, 3.0, 8.0, -5.0), 5) + def test_c3(self): + """test montecarlo 3""" + (model, opt)= self._setup_opt() + lopt= IMP.ConjugateGradients() + opt.set_local_optimizer(lopt) + self._test_starting_conditions(model, opt, (2.0, 3.0, 8.0, -5.0), 5) + + def _setup_opt(self): + model = IMP.Model() + opt = IMP.MonteCarlo() + opt.set_model(model) + for value in (-3.0, -1.0, -3.0, -1.0): + p = IMP.Particle() + model.add_particle(p) + p.add_attribute(self.xkey, value, True) + fk=IMP.FloatKeys([self.xkey]) + mod= IMP.BallMover(model.get_particles(), + fk, + .25) + opt.add_mover(mod) + + rsr = WoodsFunc() + self.rsrs.append(rsr) + model.add_restraint(rsr) + #print rsr.evaluate(None) + #print model.evaluate(False) + #model.show() + #print "Restaint TO is "+str(rsr.thisown) + #opt.show() + return (model,opt) + + def _test_starting_conditions(self, model, opt, starting_values, nrounds): + """Test the optimizer with given starting conditions""" + + print "Start energy is " + str(model.evaluate(False)) + for i in range(0,nrounds): + e=opt.optimize(100) + print "Energy after step is " + str(e) + for p in model.get_particles(): + val = p.get_value(self.xkey) + #print "Particle " + str(p.get_index().get_index()) +\ + # " is at " + str(val) + + for p in model.get_particles(): + val = p.get_value(self.xkey) + self.assertAlmostEqual(val, 1.0, places=1) + self.assertAlmostEqual(e, 0.0, places=2) + +if __name__ == '__main__': + unittest.main() Index: kernel/test/sd_optimizer/test_sd_optimizer.py =================================================================== --- kernel/test/sd_optimizer/test_sd_optimizer.py (revision 352) +++ kernel/test/sd_optimizer/test_sd_optimizer.py (working copy) @@ -1,88 +0,0 @@ -import unittest -import IMP.utils -import IMP.test, IMP -import random - -class SteepestDescentTests(IMP.test.TestCase): - """Test steepest descent optimizer""" - - def setUp(self): - """set up distance restraints and create optimizer """ - - self.imp_model = IMP.Model() - self.particles = [] - self.restraint_sets = [] - self.rsrs = [] - - # create particles 0 - 1 - self.particles.append(IMP.utils.XYZParticle(self.imp_model, - -43.0, 65.0, 93.0)) - self.particles.append(IMP.utils.XYZParticle(self.imp_model, - 20.0, 74.0, -80.0)) - self.particles.append(IMP.utils.XYZParticle(self.imp_model, - 4.0, -39.0, 26.0)) - radkey= IMP.FloatKey("radius") - - p1 = self.particles[0] - p1.add_attribute(radkey, 1.0, False) - p1 = self.particles[1] - p1.add_attribute(radkey, 2.0, False) - p1 = self.particles[2] - p1.add_attribute(radkey, 3.0, False) - - # separate 3 particles by their radii - for pairs in ((0,1), (1,2), (0,2)): - p1 = self.particles[pairs[0]] - p2 = self.particles[pairs[1]] - mean = p1.get_value(radkey) + p2.get_value(radkey) - sf = IMP.Harmonic(mean, 0.1) - rsr = IMP.DistanceRestraint(p1, p2, sf) - self.rsrs.append(rsr) - - # add restraints - rs = IMP.RestraintSet("distance_rsrs") - self.imp_model.add_restraint(rs) - self.restraint_sets.append(rs) - for i in range(len(self.rsrs)): - rs.add_restraint(self.rsrs[i]) - - self.steepest_descent = IMP.SteepestDescent() - self.steepest_descent.set_model(self.imp_model) - - - - def test_sd_optimizer1(self): - """Test that optimizer brings particles together""" - - self.steepest_descent.optimize(50) - - for i in range(0, 2): - for j in range(i+1, 3): - dist = self.particle_distance(self.particles, i, j) \ - - self.particles[i].get_value(IMP.FloatKey("radius")) \ - - self.particles[j].get_value(IMP.FloatKey("radius")) - self.assertAlmostEqual(0.0, dist, places=2) - - - def test_sd_optimizer2(self): - """Test that optimizer spreads particles apart""" - - # Start off with all particles in close proximity (but not actually - # colocated, as the derivative of zero distance is zero): - for p in self.particles: - p.set_x(random.uniform(-0.01, 0.01)) - p.set_y(random.uniform(-0.01, 0.01)) - p.set_z(random.uniform(-0.01, 0.01)) - - self.steepest_descent.optimize(50) - - for i in range(0, 2): - for j in range(i+1, 3): - dist = self.particle_distance(self.particles, i, j) \ - - self.particles[i].get_value(IMP.FloatKey("radius")) \ - - self.particles[j].get_value(IMP.FloatKey("radius")) - self.assertAlmostEqual(0.0, dist, places=2) - - -if __name__ == '__main__': - unittest.main() Index: kernel/include/IMP/optimizers/SConscript =================================================================== --- kernel/include/IMP/optimizers/SConscript (revision 352) +++ kernel/include/IMP/optimizers/SConscript (working copy) @@ -1,7 +1,8 @@ import os.path Import('env') -files = ['ConjugateGradients.h', 'SteepestDescent.h', 'MolecularDynamics.h'] +files = ['ConjugateGradients.h', 'SteepestDescent.h', 'MolecularDynamics.h', + 'MoverBase.h'] # Install the include files: includedir = os.path.join(env['includedir'], 'IMP', 'optimizers') Index: kernel/include/IMP/optimizers/MonteCarlo.h =================================================================== --- kernel/include/IMP/optimizers/MonteCarlo.h (revision 0) +++ kernel/include/IMP/optimizers/MonteCarlo.h (revision 0) @@ -0,0 +1,90 @@ +/** + * \file MonteCarlo.h \brief Simple montecarlo optimizer. + * + * Copyright 2007 Sali Lab. All rights reserved. + * + */ + +#ifndef __IMP_MONTECARLO_H +#define __IMP_MONTECARLO_H + +#include "../Optimizer.h" +#include "Mover.h" + +namespace IMP +{ + +//! Simple monte carlo optimizer. +/** The optimizer uses a set of Mover objects to propose steps. Currently + each Mover is called at each MonteCarlo iteration. This may change in + the future. The movers propose some modification, which is then + accepted or rejected based on the Metropolis criteria. Optionally, a + number of local optimization steps are taken before the step is + accepted or rejected. + */ +class IMPDLLEXPORT MonteCarlo: public Optimizer +{ +public: + MonteCarlo(); + ~MonteCarlo(); + + IMP_OPTIMIZER("Daniel Russel", "0.1"); + + IMP_CONTAINER(Mover, mover, MoverIndex); + public: + //! Return the local optimizer used or NULL + Optimizer *get_local_optimizer() const { + return cg_; + } + //! Set the local optimizer + /** The number of local steps must be nonzero for the + local optimizer to be used. + */ + void set_local_optimizer(Optimizer* cg) { + cg_=cg; + cg_->set_model(get_model()); + } + //! Set the temperature for the metropolis criteria + void set_temperature(Float t) { + IMP_assert(t>0, "Temperature must be positive"); + temp_=t; + } + Float get_temperature() const { + return temp_; + } + //! Stop if the optimization falls below this energy + void set_energy_threshold(Float t) { + stop_energy_=t; + } + Float get_energy_threshold() const { + return stop_energy_; + } + //! Take this many steps of the local optimizer for each MC step + int get_local_steps() const { + return num_local_steps_; + } + //! Take this many steps of the local optimizer for each MC step + void set_local_steps(unsigned int n) { + num_local_steps_=n; + } + //! Return how many times the optimizer has succeeded in taking a step + unsigned int get_number_of_forward_steps() const { + return stat_forward_steps_taken_; + } + + void show(std::ostream &out= std::cout) const; + private: + Float temp_; + Float prior_energy_; + Float stop_energy_; + Optimizer *cg_; + unsigned int num_local_steps_; + unsigned int stat_forward_steps_taken_; + unsigned int stat_num_failures_; +}; + +IMP_OUTPUT_OPERATOR(MonteCarlo); + +} // namespace IMP + +#endif /* __IMP_MONTECARLO_H */ Index: kernel/include/IMP/optimizers/MoverBase.h =================================================================== --- kernel/include/IMP/optimizers/MoverBase.h (revision 0) +++ kernel/include/IMP/optimizers/MoverBase.h (revision 0) @@ -0,0 +1,140 @@ +/** + * \file MoverBase.h \brief A class to help implement movers. + * + * Copyright 2007 Sali Lab. All rights reserved. + * + */ + +#ifndef __IMP_OPT_MOVERHELPER_H +#define __IMP_OPT_MOVERHELPER_H + +#include "Mover.h" +#include "../Particle.h" + +#include + +namespace IMP +{ + + //! A class to help implement movers + /** + This class helps in implementing Movers by allowing changes to be easily + rolled back. It maintains a list of particles and a list of attributes. + All changes to the product of those two lists will be rolled back + when reject_move() is called. + + See NormalMover for a simple example using this class. + */ + class IMPDLLEXPORT MoverBase: public Mover { + std::vector floats_; + std::vector ints_; + FloatKeys fkeys_; + IntKeys ikeys_; + Particles ps_; + public: + virtual void accept_move(){} + virtual void reject_move(); + + /** This sets everything up and then calls the generate_move method. + */ + virtual void propose_move(float f); + + protected: + //! implement this method to propose a move + /** + See NormalMover for a simple example. + */ + virtual void generate_move(float f)=0; + + + //! Add to the set of particles being controlled + void add_particles(const Particles &ps) { + ps_.insert(ps_.end(), ps.begin(), ps.end()); + } + //! Clear the set of particles being controlled + void clear_particles() { + ps_.clear(); + } + + //! Get the number of particles + unsigned int number_of_particles() const { + return ps_.size(); + } + Particle *get_particle(unsigned int i) const { + IMP_check(i< ps_.size(), "Bad index in MoverBase::get_particle", + IndexException()); + return ps_[i]; + } + + //! Add an attribute + unsigned int add_key(FloatKey k){ + fkeys_.push_back(k); + return fkeys_.size()-1; + } + //! Add an attribute + unsigned int add_key(IntKey k){ + ikeys_.push_back(k); + return ikeys_.size()-1; + } + + unsigned int number_of_float_keys() const {return fkeys_.size();} + unsigned int number_of_int_keys() const {return ikeys_.size();} + + + + //! Get the value of a controlled attribute + /** + \param [in] i The index of the particle. + \param [in] j The index of the attribute. + */ + Float get_float(unsigned int i, unsigned int j) const { + IMP_assert(i < ps_.size(), + "Out of range controlled attribute in guard"); + IMP_assert(j < fkeys_.size(), + "Out of range controlled attribute in guard"); + IMP_assert(ps_.size() == floats_.size(), + "Only call get_float from within generate_proposal"); + return get_particle(i)->get_value(fkeys_[j]); + } + + //! Get an int attribute value + /** + \param [in] i The index of the particle. + \param [in] j The index of the attribute. + */ + Int get_int(unsigned int i, unsigned int j) const { + IMP_assert(i < ps_.size(), + "Out of range controlled attribute in guard"); + IMP_assert(j < ikeys_.size(), + "Out of range controlled attribute in guard"); + IMP_assert(ps_.size() == ints_.size(), + "Only call get_float from within generate_proposal"); + + return get_particle(i)->get_value(ikeys_[j]); + } + //! Propose a value + /** + \param[in] i The index of the particle. + \param[in] j The index of the key + \param[in] t The value to propose + */ + void propose_value(unsigned int i, unsigned int j, Float t) { + get_particle(i)->set_value(fkeys_[j], t); + } + //! Propose a value + /** + \param[in] i The index of the particle. + \param[in] j The index of the key + \param[in] t The value to propose + */ + void propose_value(unsigned int i, unsigned int j, Int t) { + get_particle(i)->set_value(ikeys_[j], t); + } + + MoverBase(){} + ~MoverBase(){} + }; + +} + +#endif Index: kernel/include/IMP/optimizers/movers/SConscript =================================================================== --- kernel/include/IMP/optimizers/movers/SConscript (revision 0) +++ kernel/include/IMP/optimizers/movers/SConscript (revision 0) @@ -0,0 +1,9 @@ +Import('env') +import os.path +files=[] +files.append( 'BallMover.h' ) +files.append( 'DiscreteBallMover.h' ) +files.append( 'NormalMover.h' ) +includedir = os.path.join(env['includedir'], 'IMP', 'optimizers', 'movers' ) +inst = env.Install(includedir, files) +env.Alias('install', inst) Index: kernel/include/IMP/optimizers/movers/NormalMover.h =================================================================== --- kernel/include/IMP/optimizers/movers/NormalMover.h (revision 0) +++ kernel/include/IMP/optimizers/movers/NormalMover.h (revision 0) @@ -0,0 +1,41 @@ +/** + * \file NormalMover.h + * \brief A modifier which perturbs a point with a gaussian. + * + * Copyright 2007 Sali Lab. All rights reserved. + * + */ + +#ifndef __IMP_NORMALMODIFIER_H +#define __IMP_NORMALMODIFIER_H + +#include + +#include "../MoverBase.h" + +namespace IMP +{ + +//! Modify a set of continuous variables. +/** + The variables are perturbed within a ball of the + give radius. + */ +class IMPDLLEXPORT NormalMover :public MoverBase +{ +public: + NormalMover(const Particles &pis, const FloatKeys &vars, + Float stdev); + void set_particles(const Particles &ps) { + MoverBase::clear_particles(); + MoverBase::add_particles(ps); + } + protected: + virtual void generate_move(float f); + private: + Float stddev_; +}; + +} // namespace IMP + +#endif Index: kernel/include/IMP/optimizers/movers/BallMover.h =================================================================== --- kernel/include/IMP/optimizers/movers/BallMover.h (revision 0) +++ kernel/include/IMP/optimizers/movers/BallMover.h (revision 0) @@ -0,0 +1,38 @@ +/** + * \file BallMover.h \brief A modifier which perturbs a discrete variable. + * + * Copyright 2007 Sali Lab. All rights reserved. + * + */ + +#ifndef __IMP_CONTINUOUSMODIFIER_H +#define __IMP_CONTINUOUSMODIFIER_H + +#include + +#include "../MonteCarlo.h" +#include "../MoverBase.h" + +namespace IMP +{ + +//! Modify a set of continuous variables. +/** + The variables are perturbed within a ball of the + give radius. + */ +class IMPDLLEXPORT BallMover :public MoverBase +{ +public: + BallMover(const Particles &pis, const FloatKeys &vars, + Float max); + protected: + void generate_move(float a); + + private: + Float max_step_; +}; + +} // namespace IMP + +#endif Index: kernel/include/IMP/optimizers/Mover.h =================================================================== --- kernel/include/IMP/optimizers/Mover.h (revision 0) +++ kernel/include/IMP/optimizers/Mover.h (revision 0) @@ -0,0 +1,76 @@ +/** + * \file Mover.h \brief The base class for movers for MC optimization. + * + * Copyright 2007 Sali Lab. All rights reserved. + * + */ + +#ifndef __IMP_OPT_MOVER_H +#define __IMP_OPT_MOVER_H + +#include "../IMP_config.h" +#include "../base_types.h" +#include "../Object.h" + +#include + +namespace IMP +{ + class MonteCarlo; + class Mover; + typedef Index MoverIndex; + + //! A class to make a monte carlo move. + /** You probably want to use MoverBase if you are implementing a Mover. . + */ + class IMPDLLEXPORT Mover: public Object + { + friend class MonteCarlo; + void set_optimizer(MonteCarlo *c, MoverIndex i) { + opt_=c; + index_=i; + } + + MonteCarlo *opt_; + MoverIndex index_; + public: + Mover(); + virtual ~Mover(){}; + + //! propose a modification + /** + \param[in] size A number between 0 and 1 used to scale the proposed + moves. This number can be either used to scale a continuous move + or affect the probability of a discrete move. + */ + virtual void propose_move(float size)=0; + //! set whether the proposed modification is accepted + /** + \note Accepting should not change the Particles at all. + */ + virtual void accept_move()=0; + + //! Roll back any changes made to the Particles + virtual void reject_move()=0; + + //! Get a pointer to the optimizer which has this mover. + MonteCarlo *get_optimizer() const { + IMP_CHECK_OBJECT(this); + IMP_assert(opt_ != NULL, "Call set_optimizer first"); + return opt_; + } + MoverIndex get_index() const { + IMP_assert(index_!= MoverIndex(), "Call set_optimizer first"); + return index_; + } + virtual void show(std::ostream&out= std::cout) const { + out << "Mover doesn't implement show " << std::endl; + } + }; + + + IMP_OUTPUT_OPERATOR(Mover); + +} + +#endif Index: kernel/include/IMP.h =================================================================== --- kernel/include/IMP.h (revision 352) +++ kernel/include/IMP.h (working copy) @@ -36,6 +36,11 @@ #include "IMP/optimizers/SteepestDescent.h" #include "IMP/optimizers/ConjugateGradients.h" #include "IMP/optimizers/MolecularDynamics.h" +#include "IMP/optimizers/MonteCarlo.h" +#include "IMP/optimizers/Mover.h" +#include "IMP/optimizers/MoverBase.h" +#include "IMP/optimizers/movers/BallMover.h" +#include "IMP/optimizers/movers/NormalMover.h" #include "IMP/optimizers/states/VRMLLogOptimizerState.h" #include "IMP/pair_scores/DistancePairScore.h" #include "IMP/pair_scores/SphereDistancePairScore.h" Index: kernel/src/SConscript =================================================================== --- kernel/src/SConscript (revision 352) +++ kernel/src/SConscript (working copy) @@ -17,7 +17,7 @@ files = ['base_types.cpp', 'Model.cpp', 'ModelData.cpp', 'Particle.cpp', 'RigidBody.cpp', 'ScoreState.cpp', 'OptimizerState.cpp', 'Log.cpp', 'Restraint.cpp', 'Optimizer.cpp', - 'Object.cpp', 'BasicScoreFuncParams.cpp' + 'Object.cpp', 'BasicScoreFuncParams.cpp', 'random.cpp' ] + decorators_files + restraints_files + optimizers_files \ + unary_functions_files + pair_scores_files + score_states_files Index: kernel/src/random.cpp =================================================================== --- kernel/src/random.cpp (revision 0) +++ kernel/src/random.cpp (revision 0) @@ -0,0 +1,4 @@ +#include +namespace IMP { + ::boost::rand48 random_number_generator; +} Index: kernel/src/optimizers/MonteCarlo.cpp =================================================================== --- kernel/src/optimizers/MonteCarlo.cpp (revision 0) +++ kernel/src/optimizers/MonteCarlo.cpp (revision 0) @@ -0,0 +1,105 @@ +#include "IMP/optimizers/MonteCarlo.h" +#include "IMP/random.h" +#include "IMP/Model.h" + +#include +#include +#include + + +namespace IMP +{ + + + Mover::Mover(): opt_(NULL){}; + + IMP_CONTAINER_IMPL(MonteCarlo, Mover, mover, MoverIndex, + obj->set_optimizer(this, index)); + + MonteCarlo::MonteCarlo(): temp_(1), + prior_energy_(std::numeric_limits::max()), + stop_energy_(-std::numeric_limits::max()), + cg_(NULL), + num_local_steps_(50), + stat_forward_steps_taken_(0), + stat_num_failures_(0){} + + MonteCarlo::~MonteCarlo() { + IMP_CONTAINER_DELETE(Mover, mover); + if (cg_ != NULL) delete cg_; + } + + Float MonteCarlo::optimize(unsigned int max_steps) { + IMP_CHECK_OBJECT(this); + for (OptimizerStateIterator it= optimizer_states_begin(); + it != optimizer_states_end(); ++it) { + (*it)->update(); + } + prior_energy_ =get_model()->evaluate(0); + IMP_LOG(VERBOSE, "MC Initial energy is " << prior_energy_ << std::endl); + ::boost::uniform_real<> rand(0,1); + for (unsigned int i=0; i< max_steps; ++i) { + //make it a parameter + for (MoverIterator it = movers_begin(); it != movers_end(); ++it) { + IMP_LOG(VERBOSE, "MC Trying move " << **it << std::endl); + IMP_CHECK_OBJECT(*it); + (*it)->propose_move(.5); + } + Float next_energy; + if (cg_ != NULL && num_local_steps_!= 0) { + IMP_LOG(VERBOSE, + "MC Performing local optimization "<< std::flush); + IMP_CHECK_OBJECT(cg_); + next_energy =cg_->optimize(num_local_steps_); + IMP_LOG(VERBOSE, next_energy << " done "<< std::endl); + } else { + next_energy = get_model()->evaluate(0); + } + + bool accept= (next_energy < prior_energy_); + if (!accept) { + Float diff= next_energy- prior_energy_; + Float e= std::exp(-diff/temp_); + Float r= rand(random_number_generator); + IMP_LOG(VERBOSE, diff << " " << temp_ << " " << e << " " << r + << std::endl); + if (e > r) { + accept=true; + } + } + IMP_LOG(VERBOSE, "MC Prior energy is " << prior_energy_ + << " and next is " << next_energy << " "); + if (accept) { + IMP_LOG(VERBOSE, " accept" << std::endl); + for (MoverIterator it = movers_begin(); it != movers_end(); ++it) { + (*it)->accept_move(); + } + } else { + IMP_LOG(VERBOSE, " reject" << std::endl); + for (MoverIterator it = movers_begin(); it != movers_end(); ++it) { + (*it)->reject_move(); + } + } + + if (accept) { + ++stat_forward_steps_taken_; + prior_energy_= next_energy; + for (OptimizerStateIterator it= optimizer_states_begin(); + it != optimizer_states_end(); ++it) { + (*it)->update(); + } + } else { + ++stat_num_failures_; + } + if (prior_energy_ < stop_energy_) break; + } + IMP_LOG(VERBOSE, "MC Final energy is " << prior_energy_ << std::endl); + return prior_energy_; + } + + + void MonteCarlo::show(std::ostream &out) const { + out << "MonteCarlo +" << stat_forward_steps_taken_ + << " -" << stat_num_failures_ << std::endl; + } +} Index: kernel/src/optimizers/MoverBase.cpp =================================================================== --- kernel/src/optimizers/MoverBase.cpp (revision 0) +++ kernel/src/optimizers/MoverBase.cpp (revision 0) @@ -0,0 +1,36 @@ +#include "IMP/optimizers/MoverBase.h" + +namespace IMP +{ + + void MoverBase::propose_move(float f) { + if (!fkeys_.empty()) { + floats_.resize(number_of_particles(), Floats(fkeys_.size(), 0)); + } + if (!ikeys_.empty()) { + ints_.resize(number_of_particles(), Ints(ikeys_.size(), 0)); + } + for (unsigned int i=0; i< number_of_particles(); ++i) { + Particle *p = get_particle(i); + for (unsigned int j=0; j< fkeys_.size(); ++j) { + floats_[i][j]= p->get_value(fkeys_[j]); + } + for (unsigned int j=0; j< ikeys_.size(); ++j) { + ints_[i][j]= p->get_value(ikeys_[j]); + } + } + generate_move(f); + } + + void MoverBase::reject_move() { + for (unsigned int i=0; i< number_of_particles(); ++i) { + Particle *p = get_particle(i); + for (unsigned int j=0; j< fkeys_.size(); ++j) { + p->set_value(fkeys_[j], floats_[i][j]); + } + for (unsigned int j=0; j< ikeys_.size(); ++j) { + p->set_value(fkeys_[j], floats_[i][j]); + } + } + } +} Index: kernel/src/optimizers/SConscript =================================================================== --- kernel/src/optimizers/SConscript (revision 352) +++ kernel/src/optimizers/SConscript (working copy) @@ -1,9 +1,10 @@ Import('env') states_files = SConscript('states/SConscript') +movers_files = SConscript('movers/SConscript') files = [ 'SteepestDescent.cpp', 'ConjugateGradients.cpp', - 'MolecularDynamics.cpp']+states_files + 'MolecularDynamics.cpp', 'MonteCarlo.cpp', 'MoverBase.cpp']+states_files +movers_files files = [File(x) for x in files] Return('files') Index: kernel/src/optimizers/movers/NormalMover.cpp =================================================================== --- kernel/src/optimizers/movers/NormalMover.cpp (revision 0) +++ kernel/src/optimizers/movers/NormalMover.cpp (revision 0) @@ -0,0 +1,33 @@ +#include +#include +#include + +namespace IMP { + + NormalMover::NormalMover(const Particles &pis, + const FloatKeys &vars, + Float max) { + IMP_assert(max != 0, "Must have some width"); + add_particles(pis); + for (unsigned int i=0; i< vars.size(); ++i) { + add_key(vars[i]); + } + stddev_=max; + } + + void NormalMover::generate_move(float scale) { + std::vector center(number_of_float_keys()); + boost::uniform_01 u01(random_number_generator); + boost::normal_distribution mrng(0, stddev_); + 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); + IMP_assert(!std::isnan(r), "Bad random"); + IMP_assert(!std::isnan(c), "Bad stored"); + propose_value(i, j, c + r); + } + } + } + +} Index: kernel/src/optimizers/movers/BallMover.cpp =================================================================== --- kernel/src/optimizers/movers/BallMover.cpp (revision 0) +++ kernel/src/optimizers/movers/BallMover.cpp (revision 0) @@ -0,0 +1,60 @@ +#include +#include +#include + +namespace IMP { + // These functions probably should be exposed at some point + + static void random_point_in_sphere(unsigned int D, + Float radius, + Float v[]) { + IMP_assert(radius > 0, "No volume there"); + ::boost::uniform_real<> rand(-radius, radius); + Float norm; + do { + norm=0; + for (unsigned int i=0; i< D; ++i) { + v[i]= rand(random_number_generator); + norm+= v[i]*v[i]; + } + } while (norm > radius*radius); + } + + static std::vector + random_point_in_sphere(const std::vector ¢er, + Float radius) { + IMP_assert(radius > 0, "No volume there"); + Float v[center.size()]; + random_point_in_sphere(center.size(), radius, v); + std::vector r(center.size()); + for (unsigned int i=0; i< center.size(); ++i) { + r[i]=center[i]+v[i]; + } + return r; + } + + + + BallMover::BallMover(const Particles &pis, + const FloatKeys &vars, + Float max) { + add_particles(pis); + for (unsigned int i=0; i< vars.size(); ++i) { + add_key(vars[i]); + } + max_step_=max; + } + + void BallMover::generate_move( float scale ) { + std::vector center(number_of_float_keys()); + for (unsigned int i=0; i< number_of_particles(); ++i) { + for (unsigned int j=0; j< number_of_float_keys(); ++j) { + center[j]= get_float(i,j); + } + std::vector npos = random_point_in_sphere(center, scale*max_step_); + for (unsigned int j=0; j< number_of_float_keys(); ++j) { + propose_value(i, j, npos[j]); + } + } + } +} Index: kernel/src/optimizers/movers/SConscript =================================================================== --- kernel/src/optimizers/movers/SConscript (revision 0) +++ kernel/src/optimizers/movers/SConscript (revision 0) @@ -0,0 +1,6 @@ +Import('env') + +files = [ 'BallMover.cpp', 'NormalMover.cpp'] + +files = [File(x) for x in files] +Return('files') Index: kernel/pyext/IMP.i =================================================================== --- kernel/pyext/IMP.i (revision 352) +++ kernel/pyext/IMP.i (working copy) @@ -58,6 +58,12 @@ %pythonprepend SphericalDistancePairScore::SphericalDistancePairScore %{ args[1].thisown=0 %} + %pythonprepend MonteCarlo::add_mover %{ + args[1].thisown=0 + %} + %pythonprepend MonteCarlo::set_local_optimizer %{ + args[1].thisown=0 + %} } /* Don't wrap internal functions */ @@ -101,6 +107,11 @@ %include "IMP/optimizers/SteepestDescent.h" %include "IMP/optimizers/ConjugateGradients.h" %include "IMP/optimizers/MolecularDynamics.h" +%include "IMP/optimizers/Mover.h" +%include "IMP/optimizers/MoverBase.h" +%include "IMP/optimizers/MonteCarlo.h" +%include "IMP/optimizers/movers/BallMover.h" +%include "IMP/optimizers/movers/NormalMover.h" %include "IMP/optimizers/states/VRMLLogOptimizerState.h" %include "IMP/pair_scores/DistancePairScore.h" %include "IMP/pair_scores/SphereDistancePairScore.h" @@ -129,6 +140,7 @@ %template(RestraintIndex) Index; %template(ScoreStateIndex) Index; %template(OptimizerStateIndex) Index; + %template(MoverIndex) Index; %template(BondeListIndex) Index; %template(FloatKey) Key; %template(IntKey) Key; Index: tools/__init__.py =================================================================== --- tools/__init__.py (revision 352) +++ tools/__init__.py (working copy) @@ -186,7 +186,7 @@ pass env.Prepend(SCANNERS = _SWIGScanner) if env['CC'] == 'gcc': - env.Append(CCFLAGS="-Wall -Werror -g -O3") + env.Append(CCFLAGS="-Wall -g") _add_release_flags(env) sys = platform.system()