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):
58 """Setup Monte Carlo sampling
59 @param model The IMP Model
60 @param objects What to sample. Use flat list of particles
61 @param temp The MC temperature
62 @param filterbyname Not used
71 self.simulated_annealing =
False
72 self.selfadaptive =
False
83 gather_objects =
False
85 objects[0].get_particles_to_sample()
92 pts = ob.get_particles_to_sample()
95 if "Rigid_Bodies" in k:
96 mvs = self.get_rigid_body_movers(
106 mvs = self.get_super_rigid_body_movers(
114 if "Floppy_Bodies" in k:
115 mvs = self.get_floppy_body_movers(pts[k][0], pts[k][1])
121 mvs = self.get_X_movers(pts[k][0], pts[k][1])
127 if not self.isd_available:
129 "isd module needed to use nuisances")
130 mvs = self.get_nuisance_movers(pts[k][0], pts[k][1])
136 if not self.isd_available:
138 "isd module needed to use weights")
139 mvs = self.get_weight_movers(pts[k][0], pts[k][1])
145 mvs = self.get_surface_movers(
158 self.mc.set_scoring_function(get_restraint_set(self.model))
159 self.mc.set_return_best(
False)
160 self.mc.set_kt(self.temp)
161 self.mc.add_mover(self.smv)
163 def set_kt(self, temp):
170 def set_scoring_function(self, objectlist):
172 for ob
in objectlist:
173 rs.add_restraint(ob.get_restraint())
175 self.mc.set_scoring_function(sf)
177 def set_simulated_annealing(
183 self.simulated_annealing =
True
184 self.tempmin = min_temp
185 self.tempmax = max_temp
186 self.timemin = min_temp_time
187 self.timemax = max_temp_time
189 def set_self_adaptive(self, isselfadaptive=True):
190 self.selfadaptive = isselfadaptive
194 Return a dictionary with the mover parameters for nuisance parameters
197 for i
in range(self.get_number_of_movers()):
198 mv = self.smv.get_mover(i)
200 if "Nuisances" in name:
201 stepsize = IMP.core.NormalMover.get_from(mv).get_sigma()
202 output[name] = stepsize
205 def get_number_of_movers(self):
206 return len(self.smv.get_movers())
208 def get_particle_types(self):
211 def optimize(self, nstep):
213 self.mc.optimize(nstep * self.get_number_of_movers())
216 if self.simulated_annealing:
217 self.temp = self.temp_simulated_annealing()
218 self.mc.set_kt(self.temp)
221 if self.selfadaptive:
222 for i, mv
in enumerate(self.mvs):
224 mvacc = mv.get_number_of_accepted()
225 mvprp = mv.get_number_of_proposed()
226 if mv
not in self.movers_data:
227 accept = float(mvacc) / float(mvprp)
228 self.movers_data[mv] = (mvacc, mvprp)
230 oldmvacc, oldmvprp = self.movers_data[mv]
231 accept = float(mvacc-oldmvacc) / float(mvprp-oldmvprp)
232 self.movers_data[mv] = (mvacc, mvprp)
239 stepsize = mv.get_sigma()
240 if 0.4 > accept
or accept > 0.6:
241 mv.set_sigma(stepsize * 2 * accept)
244 stepsize = mv.get_radius()
245 if 0.4 > accept
or accept > 0.6:
246 mv.set_radius(stepsize * 2 * accept)
249 mr = mv.get_maximum_rotation()
250 mt = mv.get_maximum_translation()
251 if 0.4 > accept
or accept > 0.6:
252 mv.set_maximum_rotation(mr * 2 * accept)
253 mv.set_maximum_translation(mt * 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)
264 if 0.4 > accept
or accept > 0.6:
265 mv.set_radius(mr * 2 * accept)
267 def get_nuisance_movers(self, nuisances, maxstep):
269 for nuisance
in nuisances:
270 print(nuisance, maxstep)
277 def get_rigid_body_movers(self, rbs, maxtrans, maxrot):
284 def get_super_rigid_body_movers(self, rbs, maxtrans, maxrot):
291 if isinstance(rb[2], tuple)
and len(rb[2]) == 3 \
292 and isinstance(rb[2][0], float) \
293 and isinstance(rb[2][1], float) \
294 and isinstance(rb[2][2], float):
301 "Setting up a super rigid body with wrong parameters")
305 srbm.add_xyz_particle(xyz)
307 srbm.add_rigid_body_particle(rb)
311 def get_floppy_body_movers(self, fbs, maxtrans):
319 fb.set_is_optimized(fk,
True)
329 def get_X_movers(self, fbs, maxtrans):
335 raise ValueError(
"particle is part of a rigid body")
341 def get_weight_movers(self, weights, maxstep):
343 for weight
in weights:
344 if(weight.get_number_of_states() > 1):
348 def get_surface_movers(self, surfaces, maxtrans, maxrot, refprob):
350 for surface
in surfaces:
355 def temp_simulated_annealing(self):
356 if self.nframe % (self.timemin + self.timemax) < self.timemin:
360 temp = self.tempmin + (self.tempmax - self.tempmin) * value
363 def set_label(self, label):
366 def get_frame_number(self):
369 def get_output(self):
371 for i, mv
in enumerate(self.smv.get_movers()):
372 mvname = mv.get_name()
373 mvacc = mv.get_number_of_accepted()
374 mvprp = mv.get_number_of_proposed()
376 mvacr = float(mvacc) / float(mvprp)
379 output[
"MonteCarlo_Acceptance_" +
380 mvname +
"_" + str(i)] = str(mvacr)
381 if "Nuisances" in mvname:
382 output[
"MonteCarlo_StepSize_" + mvname +
"_" + str(i)] = \
383 str(IMP.core.NormalMover.get_from(mv).get_sigma())
384 if "Weights" in mvname:
385 output[
"MonteCarlo_StepSize_" + mvname +
"_" + str(i)] = \
386 str(IMP.isd.WeightMover.get_from(mv).get_radius())
387 output[
"MonteCarlo_Temperature"] = str(self.mc.get_kt())
388 output[
"MonteCarlo_Nframe"] = str(self.nframe)
393 """Sample using molecular dynamics"""
395 def __init__(self, model, objects, kt, gamma=0.01, maximum_time_step=1.0,
398 @param model The IMP Model
399 @param objects What to sample. Use flat list of particles
400 @param kt Temperature
401 @param gamma Viscosity parameter
402 @param maximum_time_step MD max time step
409 psamp = obj.get_particles_to_sample()
410 to_sample = psamp[
'Floppy_Bodies_SimplifiedModel'][0]
415 self.model, to_sample, kt/0.0019872041, gamma)
417 self.md.set_maximum_time_step(maximum_time_step)
419 self.md.set_scoring_function(sf)
421 self.md.set_scoring_function(get_restraint_set(self.model))
422 self.md.add_optimizer_state(self.ltstate)
423 self.simulated_annealing =
False
426 def set_kt(self, kt):
427 temp = kt/0.0019872041
428 self.ltstate.set_temperature(temp)
429 self.md.assign_velocities(temp)
431 def set_simulated_annealing(self, min_temp, max_temp, min_temp_time,
433 self.simulated_annealing =
True
434 self.tempmin = min_temp
435 self.tempmax = max_temp
436 self.timemin = min_temp_time
437 self.timemax = max_temp_time
439 def temp_simulated_annealing(self):
440 if self.nframe % (self.timemin + self.timemax) < self.timemin:
444 temp = self.tempmin + (self.tempmax - self.tempmin) * value
447 def set_gamma(self, gamma):
448 self.ltstate.set_gamma(gamma)
450 def optimize(self, nsteps):
453 if self.simulated_annealing:
454 self.temp = self.temp_simulated_annealing()
455 self.set_kt(self.temp)
456 self.md.optimize(nsteps)
458 def get_output(self):
460 output[
"MolecularDynamics_KineticEnergy"] = \
461 str(self.md.get_kinetic_energy())
466 """Sample using conjugate gradients"""
468 def __init__(self, model, objects):
472 self.cg.set_scoring_function(get_restraint_set(self.model))
474 def set_label(self, label):
477 def get_frame_number(self):
480 def optimize(self, nstep):
482 self.cg.optimize(nstep)
484 def set_scoring_function(self, objectlist):
486 for ob
in objectlist:
487 rs.add_restraint(ob.get_restraint())
489 self.cg.set_scoring_function(sf)
491 def get_output(self):
493 output[
"ConjugatedGradients_Nframe"] = str(self.nframe)
498 """Sample using replica exchange"""
500 def __init__(self, model, tempmin, tempmax, samplerobjects, test=True,
501 replica_exchange_object=
None):
503 samplerobjects can be a list of MonteCarlo or MolecularDynamics
507 self.samplerobjects = samplerobjects
509 self.TEMPMIN_ = tempmin
510 self.TEMPMAX_ = tempmax
512 if replica_exchange_object
is None:
516 print(
'ReplicaExchange: MPI was found. '
517 'Using Parallel Replica Exchange')
520 print(
'ReplicaExchange: Could not find MPI. '
521 'Using Serial Replica Exchange')
522 self.rem = _SerialReplicaExchange()
526 print(
'got existing rex object')
527 self.rem = replica_exchange_object
530 nproc = self.rem.get_number_of_replicas()
532 if nproc % 2 != 0
and not test:
534 "number of replicas has to be even. "
535 "set test=True to run with odd number of replicas.")
537 temp = self.rem.create_temperatures(
542 self.temperatures = temp
544 myindex = self.rem.get_my_index()
546 self.rem.set_my_parameter(
"temp", [self.temperatures[myindex]])
547 for so
in self.samplerobjects:
548 so.set_kt(self.temperatures[myindex])
554 def get_temperatures(self):
555 return self.temperatures
557 def get_my_temp(self):
558 return self.rem.get_my_parameter(
"temp")[0]
560 def get_my_index(self):
561 return self.rem.get_my_index()
563 def swap_temp(self, nframe, score=None):
565 score = self.model.evaluate(
False)
567 _ = self.rem.get_my_index()
568 mytemp = self.rem.get_my_parameter(
"temp")[0]
570 if mytemp == self.TEMPMIN_:
573 if mytemp == self.TEMPMAX_:
577 myscore = score / mytemp
580 findex = self.rem.get_friend_index(nframe)
581 ftemp = self.rem.get_friend_parameter(
"temp", findex)[0]
583 fscore = score / ftemp
586 flag = self.rem.do_exchange(myscore, fscore, findex)
591 for so
in self.samplerobjects:
595 def get_output(self):
597 if self.nattempts != 0:
598 output[
"ReplicaExchange_SwapSuccessRatio"] = str(
599 float(self.nsuccess) / self.nattempts)
600 output[
"ReplicaExchange_MinTempFrequency"] = str(
601 float(self.nmintemp) / self.nattempts)
602 output[
"ReplicaExchange_MaxTempFrequency"] = str(
603 float(self.nmaxtemp) / self.nattempts)
605 output[
"ReplicaExchange_SwapSuccessRatio"] = str(0)
606 output[
"ReplicaExchange_MinTempFrequency"] = str(0)
607 output[
"ReplicaExchange_MaxTempFrequency"] = str(0)
608 output[
"ReplicaExchange_CurrentTemp"] = str(self.get_my_temp())
612 class MPI_values(object):
613 def __init__(self, replica_exchange_object=None):
614 """Query values (ie score, and others)
615 from a set of parallel jobs"""
617 if replica_exchange_object
is None:
621 print(
'MPI_values: MPI was found. '
622 'Using Parallel Replica Exchange')
625 print(
'MPI_values: Could not find MPI. '
626 'Using Serial Replica Exchange')
627 self.rem = _SerialReplicaExchange()
631 print(
'got existing rex object')
632 self.rem = replica_exchange_object
634 def set_value(self, name, value):
635 self.rem.set_my_parameter(name, [value])
637 def get_values(self, name):
639 for i
in range(self.rem.get_number_of_replicas()):
640 v = self.rem.get_friend_parameter(name, i)[0]
644 def get_percentile(self, name):
645 value = self.rem.get_my_parameter(name)[0]
646 values = sorted(self.get_values(name))
647 ind = values.index(value)
648 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.
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 ...