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. Use flat list of particles
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()
95 pts = ob.get_particles_to_sample()
98 if "Rigid_Bodies" in k:
99 mvs = self.get_rigid_body_movers(
109 mvs = self.get_super_rigid_body_movers(
117 if "Floppy_Bodies" in k:
118 mvs = self.get_floppy_body_movers(pts[k][0], pts[k][1])
124 mvs = self.get_X_movers(pts[k][0], pts[k][1])
130 if not self.isd_available:
132 "isd module needed to use nuisances")
133 mvs = self.get_nuisance_movers(pts[k][0], pts[k][1])
139 if not self.isd_available:
141 "isd module needed to use weights")
142 mvs = self.get_weight_movers(pts[k][0], pts[k][1])
148 mvs = self.get_surface_movers(
161 self.mc.set_scoring_function(get_restraint_set(self.model))
162 self.mc.set_return_best(
False)
163 self.mc.set_score_moved(score_moved)
164 self.mc.set_kt(self.temp)
165 self.mc.add_mover(self.smv)
167 def set_kt(self, temp):
174 def set_scoring_function(self, objectlist):
176 for ob
in objectlist:
177 rs.add_restraint(ob.get_restraint())
179 self.mc.set_scoring_function(sf)
181 def set_simulated_annealing(
187 self.simulated_annealing =
True
188 self.tempmin = min_temp
189 self.tempmax = max_temp
190 self.timemin = min_temp_time
191 self.timemax = max_temp_time
193 def set_self_adaptive(self, isselfadaptive=True):
194 self.selfadaptive = isselfadaptive
198 Return a dictionary with the mover parameters for nuisance parameters
201 for i
in range(self.get_number_of_movers()):
202 mv = self.smv.get_mover(i)
204 if "Nuisances" in name:
205 stepsize = IMP.core.NormalMover.get_from(mv).get_sigma()
206 output[name] = stepsize
209 def get_number_of_movers(self):
210 return len(self.smv.get_movers())
212 def get_particle_types(self):
215 def optimize(self, nstep):
217 self.mc.optimize(nstep * self.get_number_of_movers())
220 if self.simulated_annealing:
221 self.temp = self.temp_simulated_annealing()
222 self.mc.set_kt(self.temp)
225 if self.selfadaptive:
226 for i, mv
in enumerate(self.mvs):
228 mvacc = mv.get_number_of_accepted()
229 mvprp = mv.get_number_of_proposed()
230 if mv
not in self.movers_data:
231 accept = float(mvacc) / float(mvprp)
232 self.movers_data[mv] = (mvacc, mvprp)
234 oldmvacc, oldmvprp = self.movers_data[mv]
235 accept = float(mvacc-oldmvacc) / float(mvprp-oldmvprp)
236 self.movers_data[mv] = (mvacc, mvprp)
243 stepsize = mv.get_sigma()
244 if 0.4 > accept
or accept > 0.6:
245 mv.set_sigma(stepsize * 2 * accept)
248 stepsize = mv.get_radius()
249 if 0.4 > accept
or accept > 0.6:
250 mv.set_radius(stepsize * 2 * accept)
253 mr = mv.get_maximum_rotation()
254 mt = mv.get_maximum_translation()
255 if 0.4 > accept
or accept > 0.6:
256 mv.set_maximum_rotation(mr * 2 * accept)
257 mv.set_maximum_translation(mt * 2 * accept)
260 mr = mv.get_maximum_rotation()
261 mt = mv.get_maximum_translation()
262 if 0.4 > accept
or accept > 0.6:
263 mv.set_maximum_rotation(mr * 2 * accept)
264 mv.set_maximum_translation(mt * 2 * accept)
268 if 0.4 > accept
or accept > 0.6:
269 mv.set_radius(mr * 2 * accept)
271 def get_nuisance_movers(self, nuisances, maxstep):
273 for nuisance
in nuisances:
274 print(nuisance, maxstep)
281 def get_rigid_body_movers(self, rbs, maxtrans, maxrot):
288 def get_super_rigid_body_movers(self, rbs, maxtrans, maxrot):
295 if isinstance(rb[2], tuple)
and len(rb[2]) == 3 \
296 and isinstance(rb[2][0], float) \
297 and isinstance(rb[2][1], float) \
298 and isinstance(rb[2][2], float):
305 "Setting up a super rigid body with wrong parameters")
309 srbm.add_xyz_particle(xyz)
311 srbm.add_rigid_body_particle(rb)
315 def get_floppy_body_movers(self, fbs, maxtrans):
324 fb.set_is_optimized(fk,
True)
334 def get_X_movers(self, fbs, maxtrans):
340 raise ValueError(
"particle is part of a rigid body")
346 def get_weight_movers(self, weights, maxstep):
348 for weight
in weights:
349 if weight.get_number_of_weights() > 1:
353 def get_surface_movers(self, surfaces, maxtrans, maxrot, refprob):
355 for surface
in surfaces:
360 def temp_simulated_annealing(self):
361 if self.nframe % (self.timemin + self.timemax) < self.timemin:
365 temp = self.tempmin + (self.tempmax - self.tempmin) * value
368 def set_label(self, label):
371 def get_frame_number(self):
374 def get_output(self):
376 for i, mv
in enumerate(self.smv.get_movers()):
377 mvname = mv.get_name()
378 mvacc = mv.get_number_of_accepted()
379 mvprp = mv.get_number_of_proposed()
381 mvacr = float(mvacc) / float(mvprp)
384 output[
"MonteCarlo_Acceptance_" +
385 mvname +
"_" + str(i)] = str(mvacr)
386 if "Nuisances" in mvname:
387 output[
"MonteCarlo_StepSize_" + mvname +
"_" + str(i)] = \
388 str(IMP.core.NormalMover.get_from(mv).get_sigma())
389 if "Weights" in mvname:
390 output[
"MonteCarlo_StepSize_" + mvname +
"_" + str(i)] = \
391 str(IMP.isd.WeightMover.get_from(mv).get_radius())
392 output[
"MonteCarlo_Temperature"] = str(self.mc.get_kt())
393 output[
"MonteCarlo_Nframe"] = str(self.nframe)
398 """Sample using molecular dynamics"""
400 def __init__(self, model, objects, kt, gamma=0.01, maximum_time_step=1.0,
403 @param model The IMP Model
404 @param objects What to sample. Use flat list of particles
405 @param kt Temperature
406 @param gamma Viscosity parameter
407 @param maximum_time_step MD max time step
414 psamp = obj.get_particles_to_sample()
415 to_sample = psamp[
'Floppy_Bodies_SimplifiedModel'][0]
420 self.model, to_sample, kt/0.0019872041, gamma)
422 self.md.set_maximum_time_step(maximum_time_step)
424 self.md.set_scoring_function(sf)
426 self.md.set_scoring_function(get_restraint_set(self.model))
427 self.md.add_optimizer_state(self.ltstate)
428 self.simulated_annealing =
False
431 def set_kt(self, kt):
432 temp = kt/0.0019872041
433 self.ltstate.set_temperature(temp)
434 self.md.assign_velocities(temp)
436 def set_simulated_annealing(self, min_temp, max_temp, min_temp_time,
438 self.simulated_annealing =
True
439 self.tempmin = min_temp
440 self.tempmax = max_temp
441 self.timemin = min_temp_time
442 self.timemax = max_temp_time
444 def temp_simulated_annealing(self):
445 if self.nframe % (self.timemin + self.timemax) < self.timemin:
449 temp = self.tempmin + (self.tempmax - self.tempmin) * value
452 def set_gamma(self, gamma):
453 self.ltstate.set_gamma(gamma)
455 def optimize(self, nsteps):
458 if self.simulated_annealing:
459 self.temp = self.temp_simulated_annealing()
460 self.set_kt(self.temp)
461 self.md.optimize(nsteps)
463 def get_output(self):
465 output[
"MolecularDynamics_KineticEnergy"] = \
466 str(self.md.get_kinetic_energy())
471 """Sample using conjugate gradients"""
473 def __init__(self, model, objects):
477 self.cg.set_scoring_function(get_restraint_set(self.model))
479 def set_label(self, label):
482 def get_frame_number(self):
485 def optimize(self, nstep):
487 self.cg.optimize(nstep)
489 def set_scoring_function(self, objectlist):
491 for ob
in objectlist:
492 rs.add_restraint(ob.get_restraint())
494 self.cg.set_scoring_function(sf)
496 def get_output(self):
498 output[
"ConjugatedGradients_Nframe"] = str(self.nframe)
503 """Sample using replica exchange"""
505 def __init__(self, model, tempmin, tempmax, samplerobjects, test=True,
506 replica_exchange_object=
None):
508 samplerobjects can be a list of MonteCarlo or MolecularDynamics
512 self.samplerobjects = samplerobjects
514 self.TEMPMIN_ = tempmin
515 self.TEMPMAX_ = tempmax
517 if replica_exchange_object
is None:
521 print(
'ReplicaExchange: MPI was found. '
522 'Using Parallel Replica Exchange')
525 print(
'ReplicaExchange: Could not find MPI. '
526 'Using Serial Replica Exchange')
527 self.rem = _SerialReplicaExchange()
531 print(
'got existing rex object')
532 self.rem = replica_exchange_object
535 nproc = self.rem.get_number_of_replicas()
537 if nproc % 2 != 0
and not test:
539 "number of replicas has to be even. "
540 "set test=True to run with odd number of replicas.")
542 temp = self.rem.create_temperatures(
547 self.temperatures = temp
549 myindex = self.rem.get_my_index()
551 self.rem.set_my_parameter(
"temp", [self.temperatures[myindex]])
552 for so
in self.samplerobjects:
553 so.set_kt(self.temperatures[myindex])
559 def get_temperatures(self):
560 return self.temperatures
562 def get_my_temp(self):
563 return self.rem.get_my_parameter(
"temp")[0]
565 def get_my_index(self):
566 return self.rem.get_my_index()
568 def swap_temp(self, nframe, score=None):
570 score = self.model.evaluate(
False)
572 _ = self.rem.get_my_index()
573 mytemp = self.rem.get_my_parameter(
"temp")[0]
575 if mytemp == self.TEMPMIN_:
578 if mytemp == self.TEMPMAX_:
582 myscore = score / mytemp
585 findex = self.rem.get_friend_index(nframe)
586 ftemp = self.rem.get_friend_parameter(
"temp", findex)[0]
588 fscore = score / ftemp
591 flag = self.rem.do_exchange(myscore, fscore, findex)
596 for so
in self.samplerobjects:
600 def get_output(self):
602 if self.nattempts != 0:
603 output[
"ReplicaExchange_SwapSuccessRatio"] = str(
604 float(self.nsuccess) / self.nattempts)
605 output[
"ReplicaExchange_MinTempFrequency"] = str(
606 float(self.nmintemp) / self.nattempts)
607 output[
"ReplicaExchange_MaxTempFrequency"] = str(
608 float(self.nmaxtemp) / self.nattempts)
610 output[
"ReplicaExchange_SwapSuccessRatio"] = str(0)
611 output[
"ReplicaExchange_MinTempFrequency"] = str(0)
612 output[
"ReplicaExchange_MaxTempFrequency"] = str(0)
613 output[
"ReplicaExchange_CurrentTemp"] = str(self.get_my_temp())
617 class MPI_values(object):
618 def __init__(self, replica_exchange_object=None):
619 """Query values (ie score, and others)
620 from a set of parallel jobs"""
622 if replica_exchange_object
is None:
626 print(
'MPI_values: MPI was found. '
627 'Using Parallel Replica Exchange')
630 print(
'MPI_values: Could not find MPI. '
631 'Using Serial Replica Exchange')
632 self.rem = _SerialReplicaExchange()
636 print(
'got existing rex object')
637 self.rem = replica_exchange_object
639 def set_value(self, name, value):
640 self.rem.set_my_parameter(name, [value])
642 def get_values(self, name):
644 for i
in range(self.rem.get_number_of_replicas()):
645 v = self.rem.get_friend_parameter(name, i)[0]
649 def get_percentile(self, name):
650 value = self.rem.get_my_parameter(name)[0]
651 values = sorted(self.get_values(name))
652 ind = values.index(value)
653 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.
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 ...