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
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.m))
147 self.mc.set_return_best(
False)
148 self.mc.set_kt(self.temp)
149 self.mc.add_mover(self.smv)
151 def set_kt(self, temp):
158 def set_scoring_function(self, objectlist):
160 for ob
in objectlist:
161 rs.add_restraint(ob.get_restraint())
163 self.mc.set_scoring_function(sf)
165 def set_simulated_annealing(
171 self.simulated_annealing =
True
172 self.tempmin = min_temp
173 self.tempmax = max_temp
174 self.timemin = min_temp_time
175 self.timemax = max_temp_time
177 def set_self_adaptive(self, isselfadaptive=True):
178 self.selfadaptive = isselfadaptive
182 Return a dictionary with the mover parameters for nuisance parameters
185 for i
in range(self.get_number_of_movers()):
186 mv = self.smv.get_mover(i)
188 if "Nuisances" in name:
189 stepsize = IMP.core.NormalMover.get_from(mv).get_sigma()
190 output[name] = stepsize
193 def get_number_of_movers(self):
194 return len(self.smv.get_movers())
196 def get_particle_types():
199 def optimize(self, nstep):
201 self.mc.optimize(nstep * self.get_number_of_movers())
204 if self.simulated_annealing:
205 self.temp = self.temp_simulated_annealing()
206 self.mc.set_kt(self.temp)
209 if self.selfadaptive:
210 for i, mv
in enumerate(self.mvs):
212 mvacc = mv.get_number_of_accepted()
213 mvprp = mv.get_number_of_proposed()
214 if mv
not in self.movers_data:
215 accept = float(mvacc) / float(mvprp)
216 self.movers_data[mv]=(mvacc,mvprp)
218 oldmvacc,oldmvprp=self.movers_data[mv]
219 accept = float(mvacc-oldmvacc) / float(mvprp-oldmvprp)
220 self.movers_data[mv]=(mvacc,mvprp)
221 if accept < 0.05: accept = 0.05
222 if accept > 1.0: accept = 1.0
225 stepsize = mv.get_sigma()
226 if 0.4 > accept
or accept > 0.6:
227 mv.set_sigma(stepsize * 2 * accept)
230 stepsize = mv.get_radius()
231 if 0.4 > accept
or accept > 0.6:
232 mv.set_radius(stepsize * 2 * accept)
235 mr=mv.get_maximum_rotation()
236 mt=mv.get_maximum_translation()
237 if 0.4 > accept
or accept > 0.6:
238 mv.set_maximum_rotation(mr * 2 * accept)
239 mv.set_maximum_translation(mt * 2 * accept)
242 mr=mv.get_maximum_rotation()
243 mt=mv.get_maximum_translation()
244 if 0.4 > accept
or accept > 0.6:
245 mv.set_maximum_rotation(mr * 2 * accept)
246 mv.set_maximum_translation(mt * 2 * accept)
250 if 0.4 > accept
or accept > 0.6:
251 mv.set_radius(mr * 2 * accept)
254 def run(self, *args, **kwargs):
255 self.optimize(*args, **kwargs)
257 def get_nuisance_movers(self, nuisances, maxstep):
259 for nuisance
in nuisances:
260 print(nuisance, maxstep)
267 def get_rigid_body_movers(self, rbs, maxtrans, maxrot):
274 def get_super_rigid_body_movers(self, rbs, maxtrans, maxrot):
281 if type(rb[2]) == tuple
and type(rb[2][0]) == float \
282 and type(rb[2][1]) == float
and type(rb[2][2]) == float \
300 print(
"Setting up a super rigid body with wrong parameters")
304 srbm.add_xyz_particle(xyz)
306 srbm.add_rigid_body_particle(rb)
310 def get_floppy_body_movers(self, fbs, maxtrans):
318 fb.set_is_optimized(fk,
True)
328 def get_X_movers(self, fbs, maxtrans):
334 raise ValueError(
"particle is part of a rigid body")
340 def get_weight_movers(self, weights, maxstep):
342 for weight
in weights:
343 if(weight.get_number_of_states() > 1):
347 def get_surface_movers(self, surfaces, maxtrans, maxrot, refprob):
349 for surface
in surfaces:
354 def temp_simulated_annealing(self):
355 if self.nframe % (self.timemin + self.timemax) < self.timemin:
359 temp = self.tempmin + (self.tempmax - self.tempmin) * value
362 def set_label(self, label):
365 def get_frame_number(self):
368 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 +
"_" +
383 str(i)] = str(IMP.core.NormalMover.get_from(mv).get_sigma())
384 if "Weights" in mvname:
385 output[
"MonteCarlo_StepSize_" + mvname +
"_" +
386 str(i)] = 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,m,objects,kt,gamma=0.01,maximum_time_step=1.0,sf=None):
397 @param m The IMP Model
398 @param objects What to sample. Use flat list of particles or (deprecated) 'MD Sample Objects' from PMI1
399 @param kt Temperature
400 @param gamma Viscosity parameter
401 @param maximum_time_step MD max time step
408 to_sample=obj.get_particles_to_sample()[
'Floppy_Bodies_SimplifiedModel'][0]
416 self.md.set_maximum_time_step(maximum_time_step)
418 self.md.set_scoring_function(sf)
420 self.md.set_scoring_function(get_restraint_set(self.m))
421 self.md.add_optimizer_state(self.ltstate)
422 self.simulated_annealing =
False
426 self.ltstate.set_temperature(temp)
427 self.md.assign_velocities(temp)
429 def set_simulated_annealing(self, min_temp, max_temp, min_temp_time,
431 self.simulated_annealing =
True
432 self.tempmin = min_temp
433 self.tempmax = max_temp
434 self.timemin = min_temp_time
435 self.timemax = max_temp_time
437 def temp_simulated_annealing(self):
438 if self.nframe % (self.timemin + self.timemax) < self.timemin:
442 temp = self.tempmin + (self.tempmax - self.tempmin) * value
445 def set_gamma(self,gamma):
446 self.ltstate.set_gamma(gamma)
448 def optimize(self,nsteps):
451 if self.simulated_annealing:
452 self.temp = self.temp_simulated_annealing()
453 self.set_kt(self.temp)
454 self.md.optimize(nsteps)
456 def get_output(self):
458 output[
"MolecularDynamics_KineticEnergy"]=str(self.md.get_kinetic_energy())
462 """Sample using conjugate gradients"""
464 def __init__(self, m, objects):
468 self.cg.set_scoring_function(get_restraint_set(self.m))
470 def set_label(self, label):
473 def get_frame_number(self):
477 def run(self, *args, **kwargs):
478 self.optimize(*args, **kwargs)
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):
494 output[
"ConjugatedGradients_Nframe"] = str(self.nframe)
503 """Sample using replica exchange"""
512 replica_exchange_object=
None):
514 samplerobjects can be a list of MonteCarlo or MolecularDynamics
519 self.samplerobjects = samplerobjects
521 self.TEMPMIN_ = tempmin
522 self.TEMPMAX_ = tempmax
524 if replica_exchange_object
is None:
528 print(
'ReplicaExchange: MPI was found. Using Parallel Replica Exchange')
531 print(
'ReplicaExchange: Could not find MPI. Using Serial Replica Exchange')
532 self.rem = _SerialReplicaExchange()
536 print(
'got existing rex object')
537 self.rem = replica_exchange_object
540 nproc = self.rem.get_number_of_replicas()
542 if nproc % 2 != 0
and test ==
False:
543 raise Exception(
"number of replicas has to be even. set test=True to run with odd number of replicas.")
545 temp = self.rem.create_temperatures(
550 self.temperatures = temp
552 myindex = self.rem.get_my_index()
554 self.rem.set_my_parameter(
"temp", [self.temperatures[myindex]])
555 for so
in self.samplerobjects:
556 so.set_kt(self.temperatures[myindex])
562 def get_temperatures(self):
563 return self.temperatures
565 def get_my_temp(self):
566 return self.rem.get_my_parameter(
"temp")[0]
568 def get_my_index(self):
569 return self.rem.get_my_index()
571 def swap_temp(self, nframe, score=None):
573 score = self.m.evaluate(
False)
575 myindex = self.rem.get_my_index()
576 mytemp = self.rem.get_my_parameter(
"temp")[0]
578 if mytemp == self.TEMPMIN_:
581 if mytemp == self.TEMPMAX_:
585 myscore = score / mytemp
588 findex = self.rem.get_friend_index(nframe)
589 ftemp = self.rem.get_friend_parameter(
"temp", findex)[0]
591 fscore = score / ftemp
594 flag = self.rem.do_exchange(myscore, fscore, findex)
599 for so
in self.samplerobjects:
603 def get_output(self):
606 if self.nattempts != 0:
607 output[
"ReplicaExchange_SwapSuccessRatio"] = str(
608 float(self.nsuccess) / self.nattempts)
609 output[
"ReplicaExchange_MinTempFrequency"] = str(
610 float(self.nmintemp) / self.nattempts)
611 output[
"ReplicaExchange_MaxTempFrequency"] = str(
612 float(self.nmaxtemp) / self.nattempts)
614 output[
"ReplicaExchange_SwapSuccessRatio"] = str(0)
615 output[
"ReplicaExchange_MinTempFrequency"] = str(0)
616 output[
"ReplicaExchange_MaxTempFrequency"] = str(0)
617 output[
"ReplicaExchange_CurrentTemp"] = str(self.get_my_temp())
622 class MPI_values(object):
623 def __init__(self,replica_exchange_object=None):
624 """Query values (ie score, and others)
625 from a set of parallel jobs"""
627 if replica_exchange_object
is None:
631 print(
'MPI_values: MPI was found. Using Parallel Replica Exchange')
634 print(
'MPI_values: Could not find MPI. Using Serial Replica Exchange')
635 self.rem = _SerialReplicaExchange()
639 print(
'got existing rex object')
640 self.rem = replica_exchange_object
642 def set_value(self,name,value):
643 self.rem.set_my_parameter(name,[value])
645 def get_values(self,name):
647 for i
in range(self.rem.get_number_of_replicas()):
648 v=self.rem.get_friend_parameter(name, i)[0]
652 def get_percentile(self,name):
653 value=self.rem.get_my_parameter(name)[0]
654 values=sorted(self.get_values(name))
655 ind=values.index(value)
656 percentile=float(ind)/len(values)
661 class PyMCMover(object):
664 def __init__(self, representation, mcchild, n_mc_steps):
669 self.rbs = representation.get_rigid_bodies()
672 self.n_mc_steps = n_mc_steps
674 def store_move(self):
677 for copy
in self.rbs:
680 crd.append(rb.get_reference_frame())
681 self.oldcoords.append(crd)
683 def propose_move(self, prob):
684 self.mc.run(self.n_mc_steps)
686 def reset_move(self):
688 for copy, crd
in zip(self.rbs, self.oldcoords):
689 for rb, ref
in zip(copy, crd):
690 rb.set_reference_frame(ref)
692 def get_number_of_steps(self):
693 return self.n_mc_steps
695 def set_number_of_steps(self, nsteps):
696 self.n_mc_steps = nsteps
706 self.restraints =
None
707 self.first_call =
True
710 def add_mover(self, mv):
713 def set_kt(self, kT):
716 def set_return_best(self, thing):
719 def set_move_probability(self, thing):
722 def get_energy(self):
724 pot = sum([r.evaluate(
False)
for r
in self.restraints])
726 pot = self.m.evaluate(
False)
729 def metropolis(self, old, new):
731 print(
": old %f new %f deltaE %f new_epot: %f" % (old, new, deltaE,
738 return exp(-deltaE / kT) > random.uniform(0, 1)
740 def optimize(self, nsteps):
743 print(
"=== new MC call")
747 self.first_call =
False
748 for i
in range(nsteps):
749 print(
"MC step %d " % i, end=
' ')
751 old = self.get_energy()
753 self.mv.propose_move(1)
755 new = self.get_energy()
756 if self.metropolis(old, new):
766 def get_number_of_forward_steps(self):
769 def set_restraints(self, restraints):
770 self.restraints = restraints
772 def set_scoring_function(self, objects):
776 rs.add_restraint(ob.get_restraint())
777 self.set_restraints([rs])
779 def get_output(self):
782 output[
"PyMC_Temperature"] = str(self.kT)
783 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 ...