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, model, objects=None, temp=1.0, filterbyname=None):
47 """Setup Monte Carlo sampling
48 @param model 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
61 self.simulated_annealing =
False
62 self.selfadaptive =
False
73 gather_objects =
False
75 objects[0].get_particles_to_sample()
82 pts = ob.get_particles_to_sample()
85 if "Rigid_Bodies" in k:
86 mvs = self.get_rigid_body_movers(
96 mvs = self.get_super_rigid_body_movers(
104 if "Floppy_Bodies" in k:
105 mvs = self.get_floppy_body_movers(pts[k][0], pts[k][1])
111 mvs = self.get_X_movers(pts[k][0], pts[k][1])
117 if not self.isd_available:
118 raise ValueError(
"isd module needed to use nuisances")
119 mvs = self.get_nuisance_movers(pts[k][0], pts[k][1])
125 if not self.isd_available:
126 raise ValueError(
"isd module needed to use weights")
127 mvs = self.get_weight_movers(pts[k][0], pts[k][1])
133 mvs = self.get_surface_movers(
146 self.mc.set_scoring_function(get_restraint_set(self.model))
147 self.mc.set_return_best(
False)
148 self.mc.set_kt(self.temp)
149 self.mc.add_mover(self.smv)
156 def set_kt(self, temp):
163 def set_scoring_function(self, objectlist):
165 for ob
in objectlist:
166 rs.add_restraint(ob.get_restraint())
168 self.mc.set_scoring_function(sf)
170 def set_simulated_annealing(
176 self.simulated_annealing =
True
177 self.tempmin = min_temp
178 self.tempmax = max_temp
179 self.timemin = min_temp_time
180 self.timemax = max_temp_time
182 def set_self_adaptive(self, isselfadaptive=True):
183 self.selfadaptive = isselfadaptive
187 Return a dictionary with the mover parameters for nuisance parameters
190 for i
in range(self.get_number_of_movers()):
191 mv = self.smv.get_mover(i)
193 if "Nuisances" in name:
194 stepsize = IMP.core.NormalMover.get_from(mv).get_sigma()
195 output[name] = stepsize
198 def get_number_of_movers(self):
199 return len(self.smv.get_movers())
201 def get_particle_types():
204 def optimize(self, nstep):
206 self.mc.optimize(nstep * self.get_number_of_movers())
209 if self.simulated_annealing:
210 self.temp = self.temp_simulated_annealing()
211 self.mc.set_kt(self.temp)
214 if self.selfadaptive:
215 for i, mv
in enumerate(self.mvs):
217 mvacc = mv.get_number_of_accepted()
218 mvprp = mv.get_number_of_proposed()
219 if mv
not in self.movers_data:
220 accept = float(mvacc) / float(mvprp)
221 self.movers_data[mv]=(mvacc,mvprp)
223 oldmvacc,oldmvprp=self.movers_data[mv]
224 accept = float(mvacc-oldmvacc) / float(mvprp-oldmvprp)
225 self.movers_data[mv]=(mvacc,mvprp)
226 if accept < 0.05: accept = 0.05
227 if accept > 1.0: accept = 1.0
230 stepsize = mv.get_sigma()
231 if 0.4 > accept
or accept > 0.6:
232 mv.set_sigma(stepsize * 2 * accept)
235 stepsize = mv.get_radius()
236 if 0.4 > accept
or accept > 0.6:
237 mv.set_radius(stepsize * 2 * accept)
240 mr=mv.get_maximum_rotation()
241 mt=mv.get_maximum_translation()
242 if 0.4 > accept
or accept > 0.6:
243 mv.set_maximum_rotation(mr * 2 * accept)
244 mv.set_maximum_translation(mt * 2 * accept)
247 mr=mv.get_maximum_rotation()
248 mt=mv.get_maximum_translation()
249 if 0.4 > accept
or accept > 0.6:
250 mv.set_maximum_rotation(mr * 2 * accept)
251 mv.set_maximum_translation(mt * 2 * accept)
255 if 0.4 > accept
or accept > 0.6:
256 mv.set_radius(mr * 2 * accept)
258 def get_nuisance_movers(self, nuisances, maxstep):
260 for nuisance
in nuisances:
261 print(nuisance, maxstep)
268 def get_rigid_body_movers(self, rbs, maxtrans, maxrot):
275 def get_super_rigid_body_movers(self, rbs, maxtrans, maxrot):
282 if type(rb[2]) == tuple
and type(rb[2][0]) == float \
283 and type(rb[2][1]) == float
and type(rb[2][2]) == float \
301 print(
"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):
372 for i, mv
in enumerate(self.smv.get_movers()):
373 mvname = mv.get_name()
374 mvacc = mv.get_number_of_accepted()
375 mvprp = mv.get_number_of_proposed()
377 mvacr = float(mvacc) / float(mvprp)
380 output[
"MonteCarlo_Acceptance_" +
381 mvname +
"_" + str(i)] = str(mvacr)
382 if "Nuisances" in mvname:
383 output[
"MonteCarlo_StepSize_" + mvname +
"_" +
384 str(i)] = str(IMP.core.NormalMover.get_from(mv).get_sigma())
385 if "Weights" in mvname:
386 output[
"MonteCarlo_StepSize_" + mvname +
"_" +
387 str(i)] = str(IMP.isd.WeightMover.get_from(mv).get_radius())
388 output[
"MonteCarlo_Temperature"] = str(self.mc.get_kt())
389 output[
"MonteCarlo_Nframe"] = str(self.nframe)
394 """Sample using molecular dynamics"""
396 def __init__(self,model,objects,kt,gamma=0.01,maximum_time_step=1.0,sf=None):
398 @param model The IMP Model
399 @param objects What to sample. Use flat list of particles or (deprecated) 'MD Sample Objects' from PMI1
400 @param kt Temperature
401 @param gamma Viscosity parameter
402 @param maximum_time_step MD max time step
409 to_sample=obj.get_particles_to_sample()[
'Floppy_Bodies_SimplifiedModel'][0]
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
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"]=str(self.md.get_kinetic_energy())
469 """Sample using conjugate gradients"""
471 def __init__(self, model, objects):
475 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):
502 output[
"ConjugatedGradients_Nframe"] = str(self.nframe)
511 """Sample using replica exchange"""
520 replica_exchange_object=
None):
522 samplerobjects can be a list of MonteCarlo or MolecularDynamics
527 self.samplerobjects = samplerobjects
529 self.TEMPMIN_ = tempmin
530 self.TEMPMAX_ = tempmax
532 if replica_exchange_object
is None:
536 print(
'ReplicaExchange: MPI was found. Using Parallel Replica Exchange')
539 print(
'ReplicaExchange: Could not find MPI. Using Serial Replica Exchange')
540 self.rem = _SerialReplicaExchange()
544 print(
'got existing rex object')
545 self.rem = replica_exchange_object
548 nproc = self.rem.get_number_of_replicas()
550 if nproc % 2 != 0
and test ==
False:
551 raise Exception(
"number of replicas has to be even. set test=True to run with odd number of replicas.")
553 temp = self.rem.create_temperatures(
558 self.temperatures = temp
560 myindex = self.rem.get_my_index()
562 self.rem.set_my_parameter(
"temp", [self.temperatures[myindex]])
563 for so
in self.samplerobjects:
564 so.set_kt(self.temperatures[myindex])
575 def get_temperatures(self):
576 return self.temperatures
578 def get_my_temp(self):
579 return self.rem.get_my_parameter(
"temp")[0]
581 def get_my_index(self):
582 return self.rem.get_my_index()
584 def swap_temp(self, nframe, score=None):
586 score = self.model.evaluate(
False)
588 myindex = self.rem.get_my_index()
589 mytemp = self.rem.get_my_parameter(
"temp")[0]
591 if mytemp == self.TEMPMIN_:
594 if mytemp == self.TEMPMAX_:
598 myscore = score / mytemp
601 findex = self.rem.get_friend_index(nframe)
602 ftemp = self.rem.get_friend_parameter(
"temp", findex)[0]
604 fscore = score / ftemp
607 flag = self.rem.do_exchange(myscore, fscore, findex)
612 for so
in self.samplerobjects:
616 def get_output(self):
619 if self.nattempts != 0:
620 output[
"ReplicaExchange_SwapSuccessRatio"] = str(
621 float(self.nsuccess) / self.nattempts)
622 output[
"ReplicaExchange_MinTempFrequency"] = str(
623 float(self.nmintemp) / self.nattempts)
624 output[
"ReplicaExchange_MaxTempFrequency"] = str(
625 float(self.nmaxtemp) / self.nattempts)
627 output[
"ReplicaExchange_SwapSuccessRatio"] = str(0)
628 output[
"ReplicaExchange_MinTempFrequency"] = str(0)
629 output[
"ReplicaExchange_MaxTempFrequency"] = str(0)
630 output[
"ReplicaExchange_CurrentTemp"] = str(self.get_my_temp())
635 class MPI_values(object):
636 def __init__(self,replica_exchange_object=None):
637 """Query values (ie score, and others)
638 from a set of parallel jobs"""
640 if replica_exchange_object
is None:
644 print(
'MPI_values: MPI was found. Using Parallel Replica Exchange')
647 print(
'MPI_values: Could not find MPI. Using Serial Replica Exchange')
648 self.rem = _SerialReplicaExchange()
652 print(
'got existing rex object')
653 self.rem = replica_exchange_object
655 def set_value(self,name,value):
656 self.rem.set_my_parameter(name,[value])
658 def get_values(self,name):
660 for i
in range(self.rem.get_number_of_replicas()):
661 v=self.rem.get_friend_parameter(name, i)[0]
665 def get_percentile(self,name):
666 value=self.rem.get_my_parameter(name)[0]
667 values=sorted(self.get_values(name))
668 ind=values.index(value)
669 percentile=float(ind)/len(values)
674 class PyMCMover(object):
677 def __init__(self, representation, mcchild, n_mc_steps):
682 self.rbs = representation.get_rigid_bodies()
685 self.n_mc_steps = n_mc_steps
687 def store_move(self):
690 for copy
in self.rbs:
693 crd.append(rb.get_reference_frame())
694 self.oldcoords.append(crd)
696 def propose_move(self, prob):
697 self.mc.run(self.n_mc_steps)
699 def reset_move(self):
701 for copy, crd
in zip(self.rbs, self.oldcoords):
702 for rb, ref
in zip(copy, crd):
703 rb.set_reference_frame(ref)
705 def get_number_of_steps(self):
706 return self.n_mc_steps
708 def set_number_of_steps(self, nsteps):
709 self.n_mc_steps = nsteps
719 self.restraints =
None
720 self.first_call =
True
728 def add_mover(self, mv):
731 def set_kt(self, kT):
734 def set_return_best(self, thing):
737 def set_move_probability(self, thing):
740 def get_energy(self):
742 pot = sum([r.evaluate(
False)
for r
in self.restraints])
744 pot = self.model.evaluate(
False)
747 def metropolis(self, old, new):
749 print(
": old %f new %f deltaE %f new_epot: %f" % (old, new, deltaE,
756 return exp(-deltaE / kT) > random.uniform(0, 1)
758 def optimize(self, nsteps):
761 print(
"=== new MC call")
765 self.first_call =
False
766 for i
in range(nsteps):
767 print(
"MC step %d " % i, end=
' ')
769 old = self.get_energy()
771 self.mv.propose_move(1)
773 new = self.get_energy()
774 if self.metropolis(old, new):
784 def get_number_of_forward_steps(self):
787 def set_restraints(self, restraints):
788 self.restraints = restraints
790 def set_scoring_function(self, objects):
794 rs.add_restraint(ob.get_restraint())
795 self.set_restraints([rs])
797 def get_output(self):
800 output[
"PyMC_Temperature"] = str(self.kT)
801 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.
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.
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 ...