1 """@namespace IMP.pmi1.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)
259 def run(self, *args, **kwargs):
260 self.optimize(*args, **kwargs)
262 def get_nuisance_movers(self, nuisances, maxstep):
264 for nuisance
in nuisances:
265 print(nuisance, maxstep)
272 def get_rigid_body_movers(self, rbs, maxtrans, maxrot):
279 def get_super_rigid_body_movers(self, rbs, maxtrans, maxrot):
286 if type(rb[2]) == tuple
and type(rb[2][0]) == float \
287 and type(rb[2][1]) == float
and type(rb[2][2]) == float \
305 print(
"Setting up a super rigid body with wrong parameters")
309 srbm.add_xyz_particle(xyz)
311 srbm.add_rigid_body_particle(rb)
315 def get_floppy_body_movers(self, fbs, maxtrans):
323 fb.set_is_optimized(fk,
True)
333 def get_X_movers(self, fbs, maxtrans):
339 raise ValueError(
"particle is part of a rigid body")
345 def get_weight_movers(self, weights, maxstep):
347 for weight
in weights:
348 if weight.get_number_of_weights() > 1:
352 def get_surface_movers(self, surfaces, maxtrans, maxrot, refprob):
354 for surface
in surfaces:
359 def temp_simulated_annealing(self):
360 if self.nframe % (self.timemin + self.timemax) < self.timemin:
364 temp = self.tempmin + (self.tempmax - self.tempmin) * value
367 def set_label(self, label):
370 def get_frame_number(self):
373 def get_output(self):
376 for i, mv
in enumerate(self.smv.get_movers()):
377 mvname = mv.get_name()
378 mvacc = mv.get_number_of_accepted()
379 mvprp = mv.get_number_of_proposed()
381 mvacr = float(mvacc) / float(mvprp)
384 output[
"MonteCarlo_Acceptance_" +
385 mvname +
"_" + str(i)] = str(mvacr)
386 if "Nuisances" in mvname:
387 output[
"MonteCarlo_StepSize_" + mvname +
"_" +
388 str(i)] = str(IMP.core.NormalMover.get_from(mv).get_sigma())
389 if "Weights" in mvname:
390 output[
"MonteCarlo_StepSize_" + mvname +
"_" +
391 str(i)] = str(IMP.isd.WeightMover.get_from(mv).get_radius())
392 output[
"MonteCarlo_Temperature"] = str(self.mc.get_kt())
393 output[
"MonteCarlo_Nframe"] = str(self.nframe)
398 """Sample using molecular dynamics"""
400 def __init__(self,model,objects,kt,gamma=0.01,maximum_time_step=1.0,sf=None):
402 @param model The IMP Model
403 @param objects What to sample. Use flat list of particles or (deprecated) 'MD Sample Objects' from PMI1
404 @param kt Temperature
405 @param gamma Viscosity parameter
406 @param maximum_time_step MD max time step
413 to_sample=obj.get_particles_to_sample()[
'Floppy_Bodies_SimplifiedModel'][0]
421 self.md.set_maximum_time_step(maximum_time_step)
423 self.md.set_scoring_function(sf)
425 self.md.set_scoring_function(get_restraint_set(self.model))
426 self.md.add_optimizer_state(self.ltstate)
427 self.simulated_annealing =
False
437 self.ltstate.set_temperature(temp)
438 self.md.assign_velocities(temp)
440 def set_simulated_annealing(self, min_temp, max_temp, min_temp_time,
442 self.simulated_annealing =
True
443 self.tempmin = min_temp
444 self.tempmax = max_temp
445 self.timemin = min_temp_time
446 self.timemax = max_temp_time
448 def temp_simulated_annealing(self):
449 if self.nframe % (self.timemin + self.timemax) < self.timemin:
453 temp = self.tempmin + (self.tempmax - self.tempmin) * value
456 def set_gamma(self,gamma):
457 self.ltstate.set_gamma(gamma)
459 def optimize(self,nsteps):
462 if self.simulated_annealing:
463 self.temp = self.temp_simulated_annealing()
464 self.set_kt(self.temp)
465 self.md.optimize(nsteps)
467 def get_output(self):
469 output[
"MolecularDynamics_KineticEnergy"]=str(self.md.get_kinetic_energy())
473 """Sample using conjugate gradients"""
475 def __init__(self, model, objects):
479 self.cg.set_scoring_function(get_restraint_set(self.model))
486 def set_label(self, label):
489 def get_frame_number(self):
493 def run(self, *args, **kwargs):
494 self.optimize(*args, **kwargs)
496 def optimize(self, nstep):
498 self.cg.optimize(nstep)
500 def set_scoring_function(self, objectlist):
502 for ob
in objectlist:
503 rs.add_restraint(ob.get_restraint())
505 self.cg.set_scoring_function(sf)
507 def get_output(self):
510 output[
"ConjugatedGradients_Nframe"] = str(self.nframe)
519 """Sample using replica exchange"""
528 replica_exchange_object=
None):
530 samplerobjects can be a list of MonteCarlo or MolecularDynamics
535 self.samplerobjects = samplerobjects
537 self.TEMPMIN_ = tempmin
538 self.TEMPMAX_ = tempmax
540 if replica_exchange_object
is None:
544 print(
'ReplicaExchange: MPI was found. Using Parallel Replica Exchange')
547 print(
'ReplicaExchange: Could not find MPI. Using Serial Replica Exchange')
548 self.rem = _SerialReplicaExchange()
552 print(
'got existing rex object')
553 self.rem = replica_exchange_object
556 nproc = self.rem.get_number_of_replicas()
558 if nproc % 2 != 0
and test ==
False:
559 raise Exception(
"number of replicas has to be even. set test=True to run with odd number of replicas.")
561 temp = self.rem.create_temperatures(
566 self.temperatures = temp
568 myindex = self.rem.get_my_index()
570 self.rem.set_my_parameter(
"temp", [self.temperatures[myindex]])
571 for so
in self.samplerobjects:
572 so.set_kt(self.temperatures[myindex])
583 def get_temperatures(self):
584 return self.temperatures
586 def get_my_temp(self):
587 return self.rem.get_my_parameter(
"temp")[0]
589 def get_my_index(self):
590 return self.rem.get_my_index()
592 def swap_temp(self, nframe, score=None):
594 score = self.model.evaluate(
False)
596 myindex = self.rem.get_my_index()
597 mytemp = self.rem.get_my_parameter(
"temp")[0]
599 if mytemp == self.TEMPMIN_:
602 if mytemp == self.TEMPMAX_:
606 myscore = score / mytemp
609 findex = self.rem.get_friend_index(nframe)
610 ftemp = self.rem.get_friend_parameter(
"temp", findex)[0]
612 fscore = score / ftemp
615 flag = self.rem.do_exchange(myscore, fscore, findex)
620 for so
in self.samplerobjects:
624 def get_output(self):
627 if self.nattempts != 0:
628 output[
"ReplicaExchange_SwapSuccessRatio"] = str(
629 float(self.nsuccess) / self.nattempts)
630 output[
"ReplicaExchange_MinTempFrequency"] = str(
631 float(self.nmintemp) / self.nattempts)
632 output[
"ReplicaExchange_MaxTempFrequency"] = str(
633 float(self.nmaxtemp) / self.nattempts)
635 output[
"ReplicaExchange_SwapSuccessRatio"] = str(0)
636 output[
"ReplicaExchange_MinTempFrequency"] = str(0)
637 output[
"ReplicaExchange_MaxTempFrequency"] = str(0)
638 output[
"ReplicaExchange_CurrentTemp"] = str(self.get_my_temp())
643 class MPI_values(object):
644 def __init__(self,replica_exchange_object=None):
645 """Query values (ie score, and others)
646 from a set of parallel jobs"""
648 if replica_exchange_object
is None:
652 print(
'MPI_values: MPI was found. Using Parallel Replica Exchange')
655 print(
'MPI_values: Could not find MPI. Using Serial Replica Exchange')
656 self.rem = _SerialReplicaExchange()
660 print(
'got existing rex object')
661 self.rem = replica_exchange_object
663 def set_value(self,name,value):
664 self.rem.set_my_parameter(name,[value])
666 def get_values(self,name):
668 for i
in range(self.rem.get_number_of_replicas()):
669 v=self.rem.get_friend_parameter(name, i)[0]
673 def get_percentile(self,name):
674 value=self.rem.get_my_parameter(name)[0]
675 values=sorted(self.get_values(name))
676 ind=values.index(value)
677 percentile=float(ind)/len(values)
682 class PyMCMover(object):
685 def __init__(self, representation, mcchild, n_mc_steps):
690 self.rbs = representation.get_rigid_bodies()
693 self.n_mc_steps = n_mc_steps
695 def store_move(self):
698 for copy
in self.rbs:
701 crd.append(rb.get_reference_frame())
702 self.oldcoords.append(crd)
704 def propose_move(self, prob):
705 self.mc.run(self.n_mc_steps)
707 def reset_move(self):
709 for copy, crd
in zip(self.rbs, self.oldcoords):
710 for rb, ref
in zip(copy, crd):
711 rb.set_reference_frame(ref)
713 def get_number_of_steps(self):
714 return self.n_mc_steps
716 def set_number_of_steps(self, nsteps):
717 self.n_mc_steps = nsteps
727 self.restraints =
None
728 self.first_call =
True
736 def add_mover(self, mv):
739 def set_kt(self, kT):
742 def set_return_best(self, thing):
745 def set_move_probability(self, thing):
748 def get_energy(self):
750 pot = sum([r.evaluate(
False)
for r
in self.restraints])
752 pot = self.model.evaluate(
False)
755 def metropolis(self, old, new):
757 print(
": old %f new %f deltaE %f new_epot: %f" % (old, new, deltaE,
764 return exp(-deltaE / kT) > random.uniform(0, 1)
766 def optimize(self, nsteps):
769 print(
"=== new MC call")
773 self.first_call =
False
774 for i
in range(nsteps):
775 print(
"MC step %d " % i, end=
' ')
777 old = self.get_energy()
779 self.mv.propose_move(1)
781 new = self.get_energy()
782 if self.metropolis(old, new):
792 def get_number_of_forward_steps(self):
795 def set_restraints(self, restraints):
796 self.restraints = restraints
798 def set_scoring_function(self, objects):
802 rs.add_restraint(ob.get_restraint())
803 self.set_restraints([rs])
805 def get_output(self):
808 output[
"PyMC_Temperature"] = str(self.kT)
809 output[
"PyMC_Nframe"] = str(self.nframe)
A class to implement Hamiltonian Replica Exchange.
Maintains temperature during molecular dynamics.
Sample using replica exchange.
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.
Create a scoring function on a list of restraints.
Move continuous particle variables by perturbing them within a ball.
Modify a surface orientation.
Sample using conjugate gradients.
Sample using Monte Carlo.
Object used to hold a set of restraints.
Simple molecular dynamics simulator.
def deprecated_method
Python decorator to mark a method as deprecated.
def __init__
Setup Monte Carlo sampling.
Code that uses the MPI parallel library.
A mover that perturbs a Weight particle.
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.
The general base class for IMP exceptions.
def __init__
samplerobjects can be a list of MonteCarlo or MolecularDynamics
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)
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...