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 on a single processor.
17 def get_number_of_replicas(self):
19 def create_temperatures(self, tmin, tmax, nrep):
21 def get_my_index(self):
23 def set_my_parameter(self, key, val):
24 self.__params[key] = val
25 def get_my_parameter(self, key):
26 return self.__params[key]
27 def get_friend_index(self, step):
29 def get_friend_parameter(self, key, findex):
30 return self.get_my_parameter(key)
31 def do_exchange(self, myscore, fscore, findex):
36 """Sample using Monte Carlo"""
45 def __init__(self, m, objects, temp, filterbyname=None,dof=None):
46 """Setup Monte Carlo sampling
47 @param m The IMP Model
48 @param objects Sample Objects. Alternatively pass the degrees of freedom object
49 @param temp The initial MC temperature
50 @param filterbyname Not implemented
51 @param dof DegreesOfFreedom object (replaces "objects")
63 self.simulated_annealing =
False
64 self.selfadaptive =
False
74 self.mvs = dof.get_all_movers()
79 ob.get_particles_to_sample()
81 print(
"MonteCarlo: object ", ob,
" doesn't have get_particles_to_sample() method")
83 pts = ob.get_particles_to_sample()
86 if "Rigid_Bodies" in k:
87 mvs = self.get_rigid_body_movers(
97 mvs = self.get_super_rigid_body_movers(
105 if "Floppy_Bodies" in k:
106 mvs = self.get_floppy_body_movers(pts[k][0], pts[k][1])
112 mvs = self.get_X_movers(pts[k][0], pts[k][1])
118 if not self.isd_available:
119 raise ValueError(
"isd module needed to use nuisances")
120 mvs = self.get_nuisance_movers(pts[k][0], pts[k][1])
126 if not self.isd_available:
127 raise ValueError(
"isd module needed to use weights")
128 mvs = self.get_weight_movers(pts[k][0], pts[k][1])
137 self.mc.set_scoring_function(get_restraint_set(self.m))
138 self.mc.set_return_best(
False)
139 self.mc.set_kt(self.temp)
140 self.mc.add_mover(self.smv)
142 def set_kt(self, temp):
149 def set_scoring_function(self, objectlist):
151 for ob
in objectlist:
152 rs.add_restraint(ob.get_restraint())
154 self.mc.set_scoring_function(sf)
156 def set_simulated_annealing(
162 self.simulated_annealing =
True
163 self.tempmin = min_temp
164 self.tempmax = max_temp
165 self.timemin = min_temp_time
166 self.timemax = max_temp_time
168 def set_self_adaptive(self, isselfadaptive=True):
169 self.selfadaptive = isselfadaptive
173 Return a dictionary with the mover parameters for nuisance parameters
176 for i
in range(self.get_number_of_movers()):
177 mv = self.smv.get_mover(i)
179 if "Nuisances" in name:
180 stepsize = IMP.core.NormalMover.get_from(mv).get_sigma()
181 output[name] = stepsize
184 def get_number_of_movers(self):
185 return len(self.smv.get_movers())
187 def get_particle_types():
190 def optimize(self, nstep):
192 self.mc.optimize(nstep * self.get_number_of_movers())
195 if self.simulated_annealing:
196 self.temp = self.temp_simulated_annealing()
197 self.mc.set_kt(self.temp)
200 if self.selfadaptive:
201 for i, mv
in enumerate(self.smv.get_movers()):
204 if "Nuisances" in name:
205 mvacc = mv.get_number_of_accepted()
206 mvprp = mv.get_number_of_proposed()
207 accept = float(mvacc) / float(mvprp)
208 nmv = IMP.core.NormalMover.get_from(mv)
209 stepsize = nmv.get_sigma()
211 if 0.4 > accept
or accept > 0.6:
212 nmv.set_sigma(stepsize * 2 * accept)
215 nmv.set_sigma(stepsize * 2 * accept)
218 nmv.set_sigma(stepsize * 2 * accept)
220 if "Weights" in name:
222 mvacc = mv.get_number_of_accepted()
223 mvprp = mv.get_number_of_proposed()
224 accept = float(mvacc) / float(mvprp)
225 wmv = IMP.isd.WeightMover.get_from(mv)
226 stepsize = wmv.get_radius()
228 if 0.4 > accept
or accept > 0.6:
229 wmv.set_radius(stepsize * 2 * accept)
232 wmv.set_radius(stepsize * 2 * accept)
235 wmv.set_radius(stepsize * 2 * accept)
238 def run(self, *args, **kwargs):
239 self.optimize(*args, **kwargs)
241 def get_nuisance_movers(self, nuisances, maxstep):
243 for nuisance
in nuisances:
244 print(nuisance, maxstep)
251 def get_rigid_body_movers(self, rbs, maxtrans, maxrot):
257 def get_super_rigid_body_movers(self, rbs, maxtrans, maxrot):
271 srbm.add_xyz_particle(xyz)
273 srbm.add_rigid_body_particle(rb)
277 def get_floppy_body_movers(self, fbs, maxtrans):
285 fb.set_is_optimized(fk,
True)
295 def get_X_movers(self, fbs, maxtrans):
301 raise ValueError(
"particle is part of a rigid body")
307 def get_weight_movers(self, weights, maxstep):
309 for weight
in weights:
310 if(weight.get_number_of_states() > 1):
314 def temp_simulated_annealing(self):
315 if self.nframe % (self.timemin + self.timemax) < self.timemin:
319 temp = self.tempmin + (self.tempmax - self.tempmin) * value
322 def set_label(self, label):
325 def get_frame_number(self):
328 def get_output(self):
331 for i, mv
in enumerate(self.smv.get_movers()):
332 mvname = mv.get_name()
333 mvacc = mv.get_number_of_accepted()
334 mvprp = mv.get_number_of_proposed()
336 mvacr = float(mvacc) / float(mvprp)
339 output[
"MonteCarlo_Acceptance_" +
340 mvname +
"_" + str(i)] = str(mvacr)
341 if "Nuisances" in mvname:
342 output[
"MonteCarlo_StepSize_" + mvname +
"_" +
343 str(i)] = str(IMP.core.NormalMover.get_from(mv).get_sigma())
344 if "Weights" in mvname:
345 output[
"MonteCarlo_StepSize_" + mvname +
"_" +
346 str(i)] = str(IMP.isd.WeightMover.get_from(mv).get_radius())
347 output[
"MonteCarlo_Temperature"] = str(self.mc.get_kt())
348 output[
"MonteCarlo_Nframe"] = str(self.nframe)
353 """Sample using molecular dynamics"""
355 def __init__(self,m,objects,kt,gamma=0.01,maximum_time_step=1.0):
359 to_sample+=obj.get_particles_to_sample()[
'Floppy_Bodies_SimplifiedModel'][0]
364 self.md.set_maximum_time_step(maximum_time_step)
365 self.md.add_optimizer_state(self.ltstate)
366 self.simulated_annealing =
False
370 self.ltstate.set_temperature(temp)
371 self.md.assign_velocities(temp)
373 def set_simulated_annealing(self, min_temp, max_temp, min_temp_time,
375 self.simulated_annealing =
True
376 self.tempmin = min_temp
377 self.tempmax = max_temp
378 self.timemin = min_temp_time
379 self.timemax = max_temp_time
381 def temp_simulated_annealing(self):
382 if self.nframe % (self.timemin + self.timemax) < self.timemin:
386 temp = self.tempmin + (self.tempmax - self.tempmin) * value
389 def set_gamma(self,gamma):
390 self.ltstate.set_gamma(gamma)
392 def optimize(self,nsteps):
395 if self.simulated_annealing:
396 self.temp = self.temp_simulated_annealing()
397 self.set_kt(self.temp)
398 self.md.optimize(nsteps)
400 def get_output(self):
402 output[
"MolecularDynamics_KineticEnergy"]=str(self.md.get_kinetic_energy())
406 """Sample using conjugate gradients"""
408 def __init__(self, m, objects):
412 self.cg.set_scoring_function(get_restraint_set(self.m))
414 def set_label(self, label):
417 def get_frame_number(self):
421 def run(self, *args, **kwargs):
422 self.optimize(*args, **kwargs)
424 def optimize(self, nstep):
426 self.cg.optimize(nstep)
428 def set_scoring_function(self, objectlist):
430 for ob
in objectlist:
431 rs.add_restraint(ob.get_restraint())
433 self.cg.set_scoring_function(sf)
435 def get_output(self):
438 output[
"ConjugatedGradients_Nframe"] = str(self.nframe)
443 """Sample using replica exchange"""
452 replica_exchange_object=
None):
454 samplerobjects can be a list of MonteCarlo or MolecularDynamics
459 self.samplerobjects = samplerobjects
461 self.TEMPMIN_ = tempmin
462 self.TEMPMAX_ = tempmax
464 if replica_exchange_object
is None:
468 print(
'ReplicaExchange: MPI was found. Using Parallel Replica Exchange')
471 print(
'ReplicaExchange: Could not find MPI. Using Serial Replica Exchange')
472 self.rem = _SerialReplicaExchange()
476 print(
'got existing rex object')
477 self.rem = replica_exchange_object
480 nproc = self.rem.get_number_of_replicas()
482 if nproc % 2 != 0
and test ==
False:
483 raise Exception(
"number of replicas has to be even. set test=True to run with odd number of replicas.")
485 temp = self.rem.create_temperatures(
490 self.temperatures = temp
492 myindex = self.rem.get_my_index()
494 self.rem.set_my_parameter(
"temp", [self.temperatures[myindex]])
495 for so
in self.samplerobjects:
496 so.set_kt(self.temperatures[myindex])
502 def get_temperatures(self):
503 return self.temperatures
505 def get_my_temp(self):
506 return self.rem.get_my_parameter(
"temp")[0]
508 def get_my_index(self):
509 return self.rem.get_my_index()
511 def swap_temp(self, nframe, score=None):
513 score = self.m.evaluate(
False)
515 myindex = self.rem.get_my_index()
516 mytemp = self.rem.get_my_parameter(
"temp")[0]
518 if mytemp == self.TEMPMIN_:
521 if mytemp == self.TEMPMAX_:
525 myscore = score / mytemp
528 findex = self.rem.get_friend_index(nframe)
529 ftemp = self.rem.get_friend_parameter(
"temp", findex)[0]
531 fscore = score / ftemp
534 flag = self.rem.do_exchange(myscore, fscore, findex)
539 for so
in self.samplerobjects:
543 def get_output(self):
546 if self.nattempts != 0:
547 output[
"ReplicaExchange_SwapSuccessRatio"] = str(
548 float(self.nsuccess) / self.nattempts)
549 output[
"ReplicaExchange_MinTempFrequency"] = str(
550 float(self.nmintemp) / self.nattempts)
551 output[
"ReplicaExchange_MaxTempFrequency"] = str(
552 float(self.nmaxtemp) / self.nattempts)
554 output[
"ReplicaExchange_SwapSuccessRatio"] = str(0)
555 output[
"ReplicaExchange_MinTempFrequency"] = str(0)
556 output[
"ReplicaExchange_MaxTempFrequency"] = str(0)
557 output[
"ReplicaExchange_CurrentTemp"] = str(self.get_my_temp())
561 class PyMCMover(object):
564 def __init__(self, representation, mcchild, n_mc_steps):
569 self.rbs = representation.get_rigid_bodies()
572 self.n_mc_steps = n_mc_steps
574 def store_move(self):
577 for copy
in self.rbs:
580 crd.append(rb.get_reference_frame())
581 self.oldcoords.append(crd)
583 def propose_move(self, prob):
584 self.mc.run(self.n_mc_steps)
586 def reset_move(self):
588 for copy, crd
in zip(self.rbs, self.oldcoords):
589 for rb, ref
in zip(copy, crd):
590 rb.set_reference_frame(ref)
592 def get_number_of_steps(self):
593 return self.n_mc_steps
595 def set_number_of_steps(self, nsteps):
596 self.n_mc_steps = nsteps
606 self.restraints =
None
607 self.first_call =
True
610 def add_mover(self, mv):
613 def set_kt(self, kT):
616 def set_return_best(self, thing):
619 def set_move_probability(self, thing):
622 def get_energy(self):
624 pot = sum([r.evaluate(
False)
for r
in self.restraints])
626 pot = self.m.evaluate(
False)
629 def metropolis(self, old, new):
631 print(
": old %f new %f deltaE %f new_epot: %f" % (old, new, deltaE,
638 return exp(-deltaE / kT) > random.uniform(0, 1)
640 def optimize(self, nsteps):
643 print(
"=== new MC call")
647 self.first_call =
False
648 for i
in range(nsteps):
649 print(
"MC step %d " % i, end=
' ')
651 old = self.get_energy()
653 self.mv.propose_move(1)
655 new = self.get_energy()
656 if self.metropolis(old, new):
666 def get_number_of_forward_steps(self):
669 def set_restraints(self, restraints):
670 self.restraints = restraints
672 def set_scoring_function(self, objects):
676 rs.add_restraint(ob.get_restraint())
677 self.set_restraints([rs])
679 def get_output(self):
682 output[
"PyMC_Temperature"] = str(self.kT)
683 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.
Setup constraints and create movers for an IMP Hierarchy.
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 ...