1 """@namespace IMP.pmi.samplers
2 Sampling of the system.
5 from __future__
import print_function
10 class _SerialReplicaExchange(object):
11 """Dummy replica exchange class used in non-MPI builds.
12 It should act similarly to IMP.mpi.ReplicaExchange on a single processor.
16 def get_number_of_replicas(self):
18 def create_temperatures(self, tmin, tmax, nrep):
20 def get_my_index(self):
22 def set_my_parameter(self, key, val):
23 self.__params[key] = val
24 def get_my_parameter(self, key):
25 return self.__params[key]
26 def get_friend_index(self, step):
28 def get_friend_parameter(self, key, findex):
29 return self.get_my_parameter(key)
30 def do_exchange(self, myscore, fscore, findex):
32 def set_was_used(self,was_used):
33 self.was_used = was_used
37 """Sample using Monte Carlo"""
46 def __init__(self, m, objects=None, temp=1.0, filterbyname=None):
47 """Setup Monte Carlo sampling
48 @param m The IMP Model
49 @param objects What to sample. Use flat list of particles or
50 (deprecated) 'MC Sample Objects' from PMI1
51 @param temp The MC temperature
52 @param filterbyname Not used
60 self.simulated_annealing =
False
61 self.selfadaptive =
False
71 gather_objects =
False
73 objects[0].get_particles_to_sample()
80 pts = ob.get_particles_to_sample()
83 if "Rigid_Bodies" in k:
84 mvs = self.get_rigid_body_movers(
94 mvs = self.get_super_rigid_body_movers(
102 if "Floppy_Bodies" in k:
103 mvs = self.get_floppy_body_movers(pts[k][0], pts[k][1])
109 mvs = self.get_X_movers(pts[k][0], pts[k][1])
115 if not self.isd_available:
116 raise ValueError(
"isd module needed to use nuisances")
117 mvs = self.get_nuisance_movers(pts[k][0], pts[k][1])
123 if not self.isd_available:
124 raise ValueError(
"isd module needed to use weights")
125 mvs = self.get_weight_movers(pts[k][0], pts[k][1])
134 self.mc.set_scoring_function(get_restraint_set(self.m))
135 self.mc.set_return_best(
False)
136 self.mc.set_kt(self.temp)
137 self.mc.add_mover(self.smv)
139 def set_kt(self, temp):
146 def set_scoring_function(self, objectlist):
148 for ob
in objectlist:
149 rs.add_restraint(ob.get_restraint())
151 self.mc.set_scoring_function(sf)
153 def set_simulated_annealing(
159 self.simulated_annealing =
True
160 self.tempmin = min_temp
161 self.tempmax = max_temp
162 self.timemin = min_temp_time
163 self.timemax = max_temp_time
165 def set_self_adaptive(self, isselfadaptive=True):
166 self.selfadaptive = isselfadaptive
170 Return a dictionary with the mover parameters for nuisance parameters
173 for i
in range(self.get_number_of_movers()):
174 mv = self.smv.get_mover(i)
176 if "Nuisances" in name:
177 stepsize = IMP.core.NormalMover.get_from(mv).get_sigma()
178 output[name] = stepsize
181 def get_number_of_movers(self):
182 return len(self.smv.get_movers())
184 def get_particle_types():
187 def optimize(self, nstep):
189 self.mc.optimize(nstep * self.get_number_of_movers())
192 if self.simulated_annealing:
193 self.temp = self.temp_simulated_annealing()
194 self.mc.set_kt(self.temp)
197 if self.selfadaptive:
198 for i, mv
in enumerate(self.smv.get_movers()):
201 if "Nuisances" in name:
202 mvacc = mv.get_number_of_accepted()
203 mvprp = mv.get_number_of_proposed()
204 accept = float(mvacc) / float(mvprp)
205 nmv = IMP.core.NormalMover.get_from(mv)
206 stepsize = nmv.get_sigma()
208 if 0.4 > accept
or accept > 0.6:
209 nmv.set_sigma(stepsize * 2 * accept)
212 nmv.set_sigma(stepsize * 2 * accept)
215 nmv.set_sigma(stepsize * 2 * accept)
217 if "Weights" in name:
219 mvacc = mv.get_number_of_accepted()
220 mvprp = mv.get_number_of_proposed()
221 accept = float(mvacc) / float(mvprp)
222 wmv = IMP.isd.WeightMover.get_from(mv)
223 stepsize = wmv.get_radius()
225 if 0.4 > accept
or accept > 0.6:
226 wmv.set_radius(stepsize * 2 * accept)
229 wmv.set_radius(stepsize * 2 * accept)
232 wmv.set_radius(stepsize * 2 * accept)
235 def run(self, *args, **kwargs):
236 self.optimize(*args, **kwargs)
238 def get_nuisance_movers(self, nuisances, maxstep):
240 for nuisance
in nuisances:
241 print(nuisance, maxstep)
248 def get_rigid_body_movers(self, rbs, maxtrans, maxrot):
254 def get_super_rigid_body_movers(self, rbs, maxtrans, maxrot):
268 srbm.add_xyz_particle(xyz)
270 srbm.add_rigid_body_particle(rb)
274 def get_floppy_body_movers(self, fbs, maxtrans):
282 fb.set_is_optimized(fk,
True)
292 def get_X_movers(self, fbs, maxtrans):
298 raise ValueError(
"particle is part of a rigid body")
304 def get_weight_movers(self, weights, maxstep):
306 for weight
in weights:
307 if(weight.get_number_of_states() > 1):
311 def temp_simulated_annealing(self):
312 if self.nframe % (self.timemin + self.timemax) < self.timemin:
316 temp = self.tempmin + (self.tempmax - self.tempmin) * value
319 def set_label(self, label):
322 def get_frame_number(self):
325 def get_output(self):
328 for i, mv
in enumerate(self.smv.get_movers()):
329 mvname = mv.get_name()
330 mvacc = mv.get_number_of_accepted()
331 mvprp = mv.get_number_of_proposed()
333 mvacr = float(mvacc) / float(mvprp)
336 output[
"MonteCarlo_Acceptance_" +
337 mvname +
"_" + str(i)] = str(mvacr)
338 if "Nuisances" in mvname:
339 output[
"MonteCarlo_StepSize_" + mvname +
"_" +
340 str(i)] = str(IMP.core.NormalMover.get_from(mv).get_sigma())
341 if "Weights" in mvname:
342 output[
"MonteCarlo_StepSize_" + mvname +
"_" +
343 str(i)] = str(IMP.isd.WeightMover.get_from(mv).get_radius())
344 output[
"MonteCarlo_Temperature"] = str(self.mc.get_kt())
345 output[
"MonteCarlo_Nframe"] = str(self.nframe)
350 """Sample using molecular dynamics"""
352 def __init__(self,m,objects,kt,gamma=0.01,maximum_time_step=1.0):
354 @param m The IMP Model
355 @param objects What to sample. Use flat list of particles or (deprecated) 'MD Sample Objects' from PMI1
356 @param kt Temperature
357 @param gamma Viscosity parameter
358 @param maximum_time_step MD max time step
365 to_sample=obj.get_particles_to_sample()[
'Floppy_Bodies_SimplifiedModel'][0]
373 self.md.set_maximum_time_step(maximum_time_step)
374 self.md.set_scoring_function(get_restraint_set(self.m))
375 self.md.add_optimizer_state(self.ltstate)
376 self.simulated_annealing =
False
380 self.ltstate.set_temperature(temp)
381 self.md.assign_velocities(temp)
383 def set_simulated_annealing(self, min_temp, max_temp, min_temp_time,
385 self.simulated_annealing =
True
386 self.tempmin = min_temp
387 self.tempmax = max_temp
388 self.timemin = min_temp_time
389 self.timemax = max_temp_time
391 def temp_simulated_annealing(self):
392 if self.nframe % (self.timemin + self.timemax) < self.timemin:
396 temp = self.tempmin + (self.tempmax - self.tempmin) * value
399 def set_gamma(self,gamma):
400 self.ltstate.set_gamma(gamma)
402 def optimize(self,nsteps):
405 if self.simulated_annealing:
406 self.temp = self.temp_simulated_annealing()
407 self.set_kt(self.temp)
408 self.md.optimize(nsteps)
410 def get_output(self):
412 output[
"MolecularDynamics_KineticEnergy"]=str(self.md.get_kinetic_energy())
416 """Sample using conjugate gradients"""
418 def __init__(self, m, objects):
422 self.cg.set_scoring_function(get_restraint_set(self.m))
424 def set_label(self, label):
427 def get_frame_number(self):
431 def run(self, *args, **kwargs):
432 self.optimize(*args, **kwargs)
434 def optimize(self, nstep):
436 self.cg.optimize(nstep)
438 def set_scoring_function(self, objectlist):
440 for ob
in objectlist:
441 rs.add_restraint(ob.get_restraint())
443 self.cg.set_scoring_function(sf)
445 def get_output(self):
448 output[
"ConjugatedGradients_Nframe"] = str(self.nframe)
453 """Sample using replica exchange"""
462 replica_exchange_object=
None):
464 samplerobjects can be a list of MonteCarlo or MolecularDynamics
469 self.samplerobjects = samplerobjects
471 self.TEMPMIN_ = tempmin
472 self.TEMPMAX_ = tempmax
474 if replica_exchange_object
is None:
478 print(
'ReplicaExchange: MPI was found. Using Parallel Replica Exchange')
481 print(
'ReplicaExchange: Could not find MPI. Using Serial Replica Exchange')
482 self.rem = _SerialReplicaExchange()
486 print(
'got existing rex object')
487 self.rem = replica_exchange_object
490 nproc = self.rem.get_number_of_replicas()
492 if nproc % 2 != 0
and test ==
False:
493 raise Exception(
"number of replicas has to be even. set test=True to run with odd number of replicas.")
495 temp = self.rem.create_temperatures(
500 self.temperatures = temp
502 myindex = self.rem.get_my_index()
504 self.rem.set_my_parameter(
"temp", [self.temperatures[myindex]])
505 for so
in self.samplerobjects:
506 so.set_kt(self.temperatures[myindex])
512 def get_temperatures(self):
513 return self.temperatures
515 def get_my_temp(self):
516 return self.rem.get_my_parameter(
"temp")[0]
518 def get_my_index(self):
519 return self.rem.get_my_index()
521 def swap_temp(self, nframe, score=None):
523 score = self.m.evaluate(
False)
525 myindex = self.rem.get_my_index()
526 mytemp = self.rem.get_my_parameter(
"temp")[0]
528 if mytemp == self.TEMPMIN_:
531 if mytemp == self.TEMPMAX_:
535 myscore = score / mytemp
538 findex = self.rem.get_friend_index(nframe)
539 ftemp = self.rem.get_friend_parameter(
"temp", findex)[0]
541 fscore = score / ftemp
544 flag = self.rem.do_exchange(myscore, fscore, findex)
549 for so
in self.samplerobjects:
553 def get_output(self):
556 if self.nattempts != 0:
557 output[
"ReplicaExchange_SwapSuccessRatio"] = str(
558 float(self.nsuccess) / self.nattempts)
559 output[
"ReplicaExchange_MinTempFrequency"] = str(
560 float(self.nmintemp) / self.nattempts)
561 output[
"ReplicaExchange_MaxTempFrequency"] = str(
562 float(self.nmaxtemp) / self.nattempts)
564 output[
"ReplicaExchange_SwapSuccessRatio"] = str(0)
565 output[
"ReplicaExchange_MinTempFrequency"] = str(0)
566 output[
"ReplicaExchange_MaxTempFrequency"] = str(0)
567 output[
"ReplicaExchange_CurrentTemp"] = str(self.get_my_temp())
571 class PyMCMover(object):
574 def __init__(self, representation, mcchild, n_mc_steps):
579 self.rbs = representation.get_rigid_bodies()
582 self.n_mc_steps = n_mc_steps
584 def store_move(self):
587 for copy
in self.rbs:
590 crd.append(rb.get_reference_frame())
591 self.oldcoords.append(crd)
593 def propose_move(self, prob):
594 self.mc.run(self.n_mc_steps)
596 def reset_move(self):
598 for copy, crd
in zip(self.rbs, self.oldcoords):
599 for rb, ref
in zip(copy, crd):
600 rb.set_reference_frame(ref)
602 def get_number_of_steps(self):
603 return self.n_mc_steps
605 def set_number_of_steps(self, nsteps):
606 self.n_mc_steps = nsteps
616 self.restraints =
None
617 self.first_call =
True
620 def add_mover(self, mv):
623 def set_kt(self, kT):
626 def set_return_best(self, thing):
629 def set_move_probability(self, thing):
632 def get_energy(self):
634 pot = sum([r.evaluate(
False)
for r
in self.restraints])
636 pot = self.m.evaluate(
False)
639 def metropolis(self, old, new):
641 print(
": old %f new %f deltaE %f new_epot: %f" % (old, new, deltaE,
648 return exp(-deltaE / kT) > random.uniform(0, 1)
650 def optimize(self, nsteps):
653 print(
"=== new MC call")
657 self.first_call =
False
658 for i
in range(nsteps):
659 print(
"MC step %d " % i, end=
' ')
661 old = self.get_energy()
663 self.mv.propose_move(1)
665 new = self.get_energy()
666 if self.metropolis(old, new):
676 def get_number_of_forward_steps(self):
679 def set_restraints(self, restraints):
680 self.restraints = restraints
682 def set_scoring_function(self, objects):
686 rs.add_restraint(ob.get_restraint())
687 self.set_restraints([rs])
689 def get_output(self):
692 output[
"PyMC_Temperature"] = str(self.kT)
693 output[
"PyMC_Nframe"] = str(self.nframe)
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.
Modify a set of continuous variables by perturbing them within a ball.
Object used to hold a set of restraints.
Simple molecular dynamics optimizer.
def deprecated_method
Python decorator to mark a method as deprecated.
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 ...