1 """@namespace IMP.pmi.samplers
2 Sampling of the system.
5 from __future__
import print_function
11 class _SerialReplicaExchange(object):
12 """Dummy replica exchange class used in non-MPI builds.
13 It should act similarly to IMP.mpi.ReplicaExchange
14 on a single processor.
19 def get_number_of_replicas(self):
22 def create_temperatures(self, tmin, tmax, nrep):
25 def get_my_index(self):
28 def set_my_parameter(self, key, val):
29 self.__params[key] = val
31 def get_my_parameter(self, key):
32 return self.__params[key]
34 def get_friend_index(self, step):
37 def get_friend_parameter(self, key, findex):
38 return self.get_my_parameter(key)
40 def do_exchange(self, myscore, fscore, findex):
43 def set_was_used(self, was_used):
44 self.was_used = was_used
48 """Sample using Monte Carlo"""
57 def __init__(self, model, objects=None, temp=1.0, filterbyname=None,
59 """Setup Monte Carlo sampling
60 @param model The IMP Model
61 @param objects What to sample (a list of Movers)
62 @param temp The MC temperature
63 @param filterbyname Not used
64 @param score_moved If True, attempt to speed up sampling by
65 caching scoring function terms on particles that didn't move
74 self.simulated_annealing =
False
75 self.selfadaptive =
False
86 gather_objects =
False
88 objects[0].get_particles_to_sample()
90 except AttributeError:
95 "Pass Movers to MonteCarlo(), not objects that implement "
96 "get_particles_to_sample()")
98 pts = ob.get_particles_to_sample()
101 if "Rigid_Bodies" in k:
102 mvs = self.get_rigid_body_movers(
112 mvs = self.get_super_rigid_body_movers(
120 if "Floppy_Bodies" in k:
121 mvs = self.get_floppy_body_movers(pts[k][0], pts[k][1])
127 mvs = self.get_X_movers(pts[k][0], pts[k][1])
133 if not self.isd_available:
135 "isd module needed to use nuisances")
136 mvs = self.get_nuisance_movers(pts[k][0], pts[k][1])
142 if not self.isd_available:
144 "isd module needed to use weights")
145 mvs = self.get_weight_movers(pts[k][0], pts[k][1])
151 mvs = self.get_surface_movers(
164 self.mc.set_scoring_function(get_restraint_set(self.model))
165 self.mc.set_return_best(
False)
166 self.mc.set_score_moved(score_moved)
167 self.mc.set_kt(self.temp)
168 self.mc.add_mover(self.smv)
170 def set_kt(self, temp):
177 def set_scoring_function(self, objectlist):
179 for ob
in objectlist:
180 rs.add_restraint(ob.get_restraint())
182 self.mc.set_scoring_function(sf)
184 def set_simulated_annealing(
190 self.simulated_annealing =
True
191 self.tempmin = min_temp
192 self.tempmax = max_temp
193 self.timemin = min_temp_time
194 self.timemax = max_temp_time
196 def set_self_adaptive(self, isselfadaptive=True):
197 self.selfadaptive = isselfadaptive
201 Return a dictionary with the mover parameters for nuisance parameters
204 for i
in range(self.get_number_of_movers()):
205 mv = self.smv.get_mover(i)
207 if "Nuisances" in name:
208 stepsize = IMP.core.NormalMover.get_from(mv).get_sigma()
209 output[name] = stepsize
212 def get_number_of_movers(self):
213 return len(self.smv.get_movers())
215 def get_particle_types(self):
218 def optimize(self, nstep):
220 self.mc.optimize(nstep * self.get_number_of_movers())
223 if self.simulated_annealing:
224 self.temp = self.temp_simulated_annealing()
225 self.mc.set_kt(self.temp)
228 if self.selfadaptive:
229 for i, mv
in enumerate(self.mvs):
231 mvacc = mv.get_number_of_accepted()
232 mvprp = mv.get_number_of_proposed()
233 if mv
not in self.movers_data:
234 accept = float(mvacc) / float(mvprp)
235 self.movers_data[mv] = (mvacc, mvprp)
237 oldmvacc, oldmvprp = self.movers_data[mv]
238 accept = float(mvacc-oldmvacc) / float(mvprp-oldmvprp)
239 self.movers_data[mv] = (mvacc, mvprp)
246 stepsize = mv.get_sigma()
247 if 0.4 > accept
or accept > 0.6:
248 mv.set_sigma(stepsize * 2 * accept)
251 stepsize = mv.get_radius()
252 if 0.4 > accept
or accept > 0.6:
253 mv.set_radius(stepsize * 2 * accept)
256 mr = mv.get_maximum_rotation()
257 mt = mv.get_maximum_translation()
258 if 0.4 > accept
or accept > 0.6:
259 mv.set_maximum_rotation(mr * 2 * accept)
260 mv.set_maximum_translation(mt * 2 * accept)
263 mr = mv.get_maximum_rotation()
264 mt = mv.get_maximum_translation()
265 if 0.4 > accept
or accept > 0.6:
266 mv.set_maximum_rotation(mr * 2 * accept)
267 mv.set_maximum_translation(mt * 2 * accept)
271 if 0.4 > accept
or accept > 0.6:
272 mv.set_radius(mr * 2 * accept)
274 def get_nuisance_movers(self, nuisances, maxstep):
276 for nuisance
in nuisances:
277 print(nuisance, maxstep)
284 def get_rigid_body_movers(self, rbs, maxtrans, maxrot):
291 def get_super_rigid_body_movers(self, rbs, maxtrans, maxrot):
298 if isinstance(rb[2], tuple)
and len(rb[2]) == 3 \
299 and isinstance(rb[2][0], float) \
300 and isinstance(rb[2][1], float) \
301 and isinstance(rb[2][2], float):
308 "Setting up a super rigid body with wrong parameters")
312 srbm.add_xyz_particle(xyz)
314 srbm.add_rigid_body_particle(rb)
318 def get_floppy_body_movers(self, fbs, maxtrans):
327 fb.set_is_optimized(fk,
True)
337 def get_X_movers(self, fbs, maxtrans):
343 raise ValueError(
"particle is part of a rigid body")
349 def get_weight_movers(self, weights, maxstep):
351 for weight
in weights:
352 if weight.get_number_of_weights() > 1:
356 def get_surface_movers(self, surfaces, maxtrans, maxrot, refprob):
358 for surface
in surfaces:
363 def temp_simulated_annealing(self):
364 if self.nframe % (self.timemin + self.timemax) < self.timemin:
368 temp = self.tempmin + (self.tempmax - self.tempmin) * value
371 def set_label(self, label):
374 def get_frame_number(self):
377 def get_output(self):
379 for i, mv
in enumerate(self.smv.get_movers()):
380 mvname = mv.get_name()
381 mvacc = mv.get_number_of_accepted()
382 mvprp = mv.get_number_of_proposed()
384 mvacr = float(mvacc) / float(mvprp)
387 output[
"MonteCarlo_Acceptance_" +
388 mvname +
"_" + str(i)] = str(mvacr)
389 if "Nuisances" in mvname:
390 output[
"MonteCarlo_StepSize_" + mvname +
"_" + str(i)] = \
391 str(IMP.core.NormalMover.get_from(mv).get_sigma())
392 if "Weights" in mvname:
393 output[
"MonteCarlo_StepSize_" + mvname +
"_" + str(i)] = \
394 str(IMP.isd.WeightMover.get_from(mv).get_radius())
395 output[
"MonteCarlo_Temperature"] = str(self.mc.get_kt())
396 output[
"MonteCarlo_Nframe"] = str(self.nframe)
401 """Sample using molecular dynamics"""
403 def __init__(self, model, objects, kt, gamma=0.01, maximum_time_step=1.0,
406 @param model The IMP Model
407 @param objects What to sample. Use flat list of particles
408 @param kt Temperature
409 @param gamma Viscosity parameter
410 @param maximum_time_step MD max time step
417 psamp = obj.get_particles_to_sample()
418 to_sample = psamp[
'Floppy_Bodies_SimplifiedModel'][0]
423 self.model, to_sample, kt/0.0019872041, gamma)
425 self.md.set_maximum_time_step(maximum_time_step)
427 self.md.set_scoring_function(sf)
429 self.md.set_scoring_function(get_restraint_set(self.model))
430 self.md.add_optimizer_state(self.ltstate)
431 self.simulated_annealing =
False
434 def set_kt(self, kt):
435 temp = kt/0.0019872041
436 self.ltstate.set_temperature(temp)
437 self.md.assign_velocities(temp)
439 def set_simulated_annealing(self, min_temp, max_temp, min_temp_time,
441 self.simulated_annealing =
True
442 self.tempmin = min_temp
443 self.tempmax = max_temp
444 self.timemin = min_temp_time
445 self.timemax = max_temp_time
447 def temp_simulated_annealing(self):
448 if self.nframe % (self.timemin + self.timemax) < self.timemin:
452 temp = self.tempmin + (self.tempmax - self.tempmin) * value
455 def set_gamma(self, gamma):
456 self.ltstate.set_gamma(gamma)
458 def optimize(self, nsteps):
461 if self.simulated_annealing:
462 self.temp = self.temp_simulated_annealing()
463 self.set_kt(self.temp)
464 self.md.optimize(nsteps)
466 def get_output(self):
468 output[
"MolecularDynamics_KineticEnergy"] = \
469 str(self.md.get_kinetic_energy())
474 """Sample using conjugate gradients"""
476 def __init__(self, model, objects):
480 self.cg.set_scoring_function(get_restraint_set(self.model))
482 def set_label(self, label):
485 def get_frame_number(self):
488 def optimize(self, nstep):
490 self.cg.optimize(nstep)
492 def set_scoring_function(self, objectlist):
494 for ob
in objectlist:
495 rs.add_restraint(ob.get_restraint())
497 self.cg.set_scoring_function(sf)
499 def get_output(self):
501 output[
"ConjugatedGradients_Nframe"] = str(self.nframe)
506 """Sample using replica exchange"""
508 def __init__(self, model, tempmin, tempmax, samplerobjects, test=True,
509 replica_exchange_object=
None):
511 samplerobjects can be a list of MonteCarlo or MolecularDynamics
515 self.samplerobjects = samplerobjects
517 self.TEMPMIN_ = tempmin
518 self.TEMPMAX_ = tempmax
520 if replica_exchange_object
is None:
524 print(
'ReplicaExchange: MPI was found. '
525 'Using Parallel Replica Exchange')
528 print(
'ReplicaExchange: Could not find MPI. '
529 'Using Serial Replica Exchange')
530 self.rem = _SerialReplicaExchange()
534 print(
'got existing rex object')
535 self.rem = replica_exchange_object
538 nproc = self.rem.get_number_of_replicas()
540 if nproc % 2 != 0
and not test:
542 "number of replicas has to be even. "
543 "set test=True to run with odd number of replicas.")
545 temp = self.rem.create_temperatures(
550 self.temperatures = temp
552 myindex = self.rem.get_my_index()
554 self.rem.set_my_parameter(
"temp", [self.temperatures[myindex]])
555 for so
in self.samplerobjects:
556 so.set_kt(self.temperatures[myindex])
562 def get_temperatures(self):
563 return self.temperatures
565 def get_my_temp(self):
566 return self.rem.get_my_parameter(
"temp")[0]
568 def get_my_index(self):
569 return self.rem.get_my_index()
571 def swap_temp(self, nframe, score=None):
573 score = self.model.evaluate(
False)
575 _ = self.rem.get_my_index()
576 mytemp = self.rem.get_my_parameter(
"temp")[0]
578 if mytemp == self.TEMPMIN_:
581 if mytemp == self.TEMPMAX_:
585 myscore = score / mytemp
588 findex = self.rem.get_friend_index(nframe)
589 ftemp = self.rem.get_friend_parameter(
"temp", findex)[0]
591 fscore = score / ftemp
594 flag = self.rem.do_exchange(myscore, fscore, findex)
599 for so
in self.samplerobjects:
603 def get_output(self):
605 if self.nattempts != 0:
606 output[
"ReplicaExchange_SwapSuccessRatio"] = str(
607 float(self.nsuccess) / self.nattempts)
608 output[
"ReplicaExchange_MinTempFrequency"] = str(
609 float(self.nmintemp) / self.nattempts)
610 output[
"ReplicaExchange_MaxTempFrequency"] = str(
611 float(self.nmaxtemp) / self.nattempts)
613 output[
"ReplicaExchange_SwapSuccessRatio"] = str(0)
614 output[
"ReplicaExchange_MinTempFrequency"] = str(0)
615 output[
"ReplicaExchange_MaxTempFrequency"] = str(0)
616 output[
"ReplicaExchange_CurrentTemp"] = str(self.get_my_temp())
620 class MPI_values(object):
621 def __init__(self, replica_exchange_object=None):
622 """Query values (ie score, and others)
623 from a set of parallel jobs"""
625 if replica_exchange_object
is None:
629 print(
'MPI_values: MPI was found. '
630 'Using Parallel Replica Exchange')
633 print(
'MPI_values: Could not find MPI. '
634 'Using Serial Replica Exchange')
635 self.rem = _SerialReplicaExchange()
639 print(
'got existing rex object')
640 self.rem = replica_exchange_object
642 def set_value(self, name, value):
643 self.rem.set_my_parameter(name, [value])
645 def get_values(self, name):
647 for i
in range(self.rem.get_number_of_replicas()):
648 v = self.rem.get_friend_parameter(name, i)[0]
652 def get_percentile(self, name):
653 value = self.rem.get_my_parameter(name)[0]
654 values = sorted(self.get_values(name))
655 ind = values.index(value)
656 percentile = float(ind)/len(values)
def __init__
samplerobjects can be a list of MonteCarlo or MolecularDynamics
A class to implement Hamiltonian Replica Exchange.
Maintains temperature during molecular dynamics.
Sample using molecular dynamics.
def get_nuisance_movers_parameters
Return a dictionary with the mover parameters for nuisance parameters.
void handle_use_deprecated(std::string message)
Break in this method in gdb to find deprecated uses at runtime.
Modify the transformation of a rigid body.
Simple conjugate gradients optimizer.
Sample using conjugate gradients.
Create a scoring function on a list of restraints.
Move continuous particle variables by perturbing them within a ball.
Modify a surface orientation.
static FloatKeys get_internal_coordinate_keys()
Object used to hold a set of restraints.
Simple molecular dynamics optimizer.
Code that uses the MPI parallel library.
A mover that perturbs a Weight particle.
def __init__
Setup Monte Carlo sampling.
Modify a set of continuous variables using a normal distribution.
Basic functionality that is expected to be used by a wide variety of IMP users.
Sample using Monte Carlo.
The general base class for IMP exceptions.
static const FloatKeys & get_xyz_keys()
Get a vector containing the keys for x,y,z.
Applies a list of movers one at a time.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Sample using replica exchange.
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...