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)
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(self):
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)
229 if isinstance(mv, MP.isd.WeightMover):
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)
253 def get_nuisance_movers(self, nuisances, maxstep):
255 for nuisance
in nuisances:
256 print(nuisance, maxstep)
263 def get_rigid_body_movers(self, rbs, maxtrans, maxrot):
270 def get_super_rigid_body_movers(self, rbs, maxtrans, maxrot):
277 if isinstance(rb[2], tuple)
and len(rb[2]) == 3 \
278 and isinstance(rb[2][0], float) \
279 and isinstance(rb[2][1], float) \
280 and isinstance(rb[2][2], float):
297 print(
"Setting up a super rigid body with wrong parameters")
301 srbm.add_xyz_particle(xyz)
303 srbm.add_rigid_body_particle(rb)
307 def get_floppy_body_movers(self, fbs, maxtrans):
315 fb.set_is_optimized(fk,
True)
325 def get_X_movers(self, fbs, maxtrans):
331 raise ValueError(
"particle is part of a rigid body")
337 def get_weight_movers(self, weights, maxstep):
339 for weight
in weights:
340 if(weight.get_number_of_states() > 1):
344 def get_surface_movers(self, surfaces, maxtrans, maxrot, refprob):
346 for surface
in surfaces:
351 def temp_simulated_annealing(self):
352 if self.nframe % (self.timemin + self.timemax) < self.timemin:
356 temp = self.tempmin + (self.tempmax - self.tempmin) * value
359 def set_label(self, label):
362 def get_frame_number(self):
365 def get_output(self):
368 for i, mv
in enumerate(self.smv.get_movers()):
369 mvname = mv.get_name()
370 mvacc = mv.get_number_of_accepted()
371 mvprp = mv.get_number_of_proposed()
373 mvacr = float(mvacc) / float(mvprp)
376 output[
"MonteCarlo_Acceptance_" +
377 mvname +
"_" + str(i)] = str(mvacr)
378 if "Nuisances" in mvname:
379 output[
"MonteCarlo_StepSize_" + mvname +
"_" +
380 str(i)] = str(IMP.core.NormalMover.get_from(mv).get_sigma())
381 if "Weights" in mvname:
382 output[
"MonteCarlo_StepSize_" + mvname +
"_" +
383 str(i)] = str(IMP.isd.WeightMover.get_from(mv).get_radius())
384 output[
"MonteCarlo_Temperature"] = str(self.mc.get_kt())
385 output[
"MonteCarlo_Nframe"] = str(self.nframe)
390 """Sample using molecular dynamics"""
392 def __init__(self,model,objects,kt,gamma=0.01,maximum_time_step=1.0,sf=None):
394 @param model The IMP Model
395 @param objects What to sample. Use flat list of particles or (deprecated) 'MD Sample Objects' from PMI1
396 @param kt Temperature
397 @param gamma Viscosity parameter
398 @param maximum_time_step MD max time step
405 to_sample=obj.get_particles_to_sample()[
'Floppy_Bodies_SimplifiedModel'][0]
413 self.md.set_maximum_time_step(maximum_time_step)
415 self.md.set_scoring_function(sf)
417 self.md.set_scoring_function(get_restraint_set(self.model))
418 self.md.add_optimizer_state(self.ltstate)
419 self.simulated_annealing =
False
424 self.ltstate.set_temperature(temp)
425 self.md.assign_velocities(temp)
427 def set_simulated_annealing(self, min_temp, max_temp, min_temp_time,
429 self.simulated_annealing =
True
430 self.tempmin = min_temp
431 self.tempmax = max_temp
432 self.timemin = min_temp_time
433 self.timemax = max_temp_time
435 def temp_simulated_annealing(self):
436 if self.nframe % (self.timemin + self.timemax) < self.timemin:
440 temp = self.tempmin + (self.tempmax - self.tempmin) * value
443 def set_gamma(self,gamma):
444 self.ltstate.set_gamma(gamma)
446 def optimize(self,nsteps):
449 if self.simulated_annealing:
450 self.temp = self.temp_simulated_annealing()
451 self.set_kt(self.temp)
452 self.md.optimize(nsteps)
454 def get_output(self):
456 output[
"MolecularDynamics_KineticEnergy"]=str(self.md.get_kinetic_energy())
460 """Sample using conjugate gradients"""
462 def __init__(self, model, objects):
466 self.cg.set_scoring_function(get_restraint_set(self.model))
468 def set_label(self, label):
471 def get_frame_number(self):
474 def optimize(self, nstep):
476 self.cg.optimize(nstep)
478 def set_scoring_function(self, objectlist):
480 for ob
in objectlist:
481 rs.add_restraint(ob.get_restraint())
483 self.cg.set_scoring_function(sf)
485 def get_output(self):
488 output[
"ConjugatedGradients_Nframe"] = str(self.nframe)
497 """Sample using replica exchange"""
506 replica_exchange_object=
None):
508 samplerobjects can be a list of MonteCarlo or MolecularDynamics
513 self.samplerobjects = samplerobjects
515 self.TEMPMIN_ = tempmin
516 self.TEMPMAX_ = tempmax
518 if replica_exchange_object
is None:
522 print(
'ReplicaExchange: MPI was found. Using Parallel Replica Exchange')
525 print(
'ReplicaExchange: Could not find MPI. Using Serial Replica Exchange')
526 self.rem = _SerialReplicaExchange()
530 print(
'got existing rex object')
531 self.rem = replica_exchange_object
534 nproc = self.rem.get_number_of_replicas()
536 if nproc % 2 != 0
and test ==
False:
537 raise Exception(
"number of replicas has to be even. set test=True to run with odd number of replicas.")
539 temp = self.rem.create_temperatures(
544 self.temperatures = temp
546 myindex = self.rem.get_my_index()
548 self.rem.set_my_parameter(
"temp", [self.temperatures[myindex]])
549 for so
in self.samplerobjects:
550 so.set_kt(self.temperatures[myindex])
556 def get_temperatures(self):
557 return self.temperatures
559 def get_my_temp(self):
560 return self.rem.get_my_parameter(
"temp")[0]
562 def get_my_index(self):
563 return self.rem.get_my_index()
565 def swap_temp(self, nframe, score=None):
567 score = self.model.evaluate(
False)
569 myindex = self.rem.get_my_index()
570 mytemp = self.rem.get_my_parameter(
"temp")[0]
572 if mytemp == self.TEMPMIN_:
575 if mytemp == self.TEMPMAX_:
579 myscore = score / mytemp
582 findex = self.rem.get_friend_index(nframe)
583 ftemp = self.rem.get_friend_parameter(
"temp", findex)[0]
585 fscore = score / ftemp
588 flag = self.rem.do_exchange(myscore, fscore, findex)
593 for so
in self.samplerobjects:
597 def get_output(self):
600 if self.nattempts != 0:
601 output[
"ReplicaExchange_SwapSuccessRatio"] = str(
602 float(self.nsuccess) / self.nattempts)
603 output[
"ReplicaExchange_MinTempFrequency"] = str(
604 float(self.nmintemp) / self.nattempts)
605 output[
"ReplicaExchange_MaxTempFrequency"] = str(
606 float(self.nmaxtemp) / self.nattempts)
608 output[
"ReplicaExchange_SwapSuccessRatio"] = str(0)
609 output[
"ReplicaExchange_MinTempFrequency"] = str(0)
610 output[
"ReplicaExchange_MaxTempFrequency"] = str(0)
611 output[
"ReplicaExchange_CurrentTemp"] = str(self.get_my_temp())
616 class MPI_values(object):
617 def __init__(self,replica_exchange_object=None):
618 """Query values (ie score, and others)
619 from a set of parallel jobs"""
621 if replica_exchange_object
is None:
625 print(
'MPI_values: MPI was found. Using Parallel Replica Exchange')
628 print(
'MPI_values: Could not find MPI. Using Serial Replica Exchange')
629 self.rem = _SerialReplicaExchange()
633 print(
'got existing rex object')
634 self.rem = replica_exchange_object
636 def set_value(self,name,value):
637 self.rem.set_my_parameter(name,[value])
639 def get_values(self,name):
641 for i
in range(self.rem.get_number_of_replicas()):
642 v=self.rem.get_friend_parameter(name, i)[0]
646 def get_percentile(self,name):
647 value=self.rem.get_my_parameter(name)[0]
648 values=sorted(self.get_values(name))
649 ind=values.index(value)
650 percentile=float(ind)/len(values)
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.
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 ...