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
14 on a single processor.
19 def get_number_of_replicas(self):
22 def create_temperatures(self, tmin, tmax, nrep):
25 def get_my_index(self):
28 def set_my_parameter(self, key, val):
29 self.__params[key] = val
31 def get_my_parameter(self, key):
32 return self.__params[key]
34 def get_friend_index(self, step):
37 def get_friend_parameter(self, key, findex):
38 return self.get_my_parameter(key)
40 def do_exchange(self, myscore, fscore, findex):
43 def set_was_used(self, was_used):
44 self.was_used = was_used
48 """Sample using Monte Carlo"""
57 def __init__(self, model, objects=None, temp=1.0, filterbyname=None,
59 """Setup Monte Carlo sampling
60 @param model The IMP Model
61 @param objects What to sample (a list of Movers)
62 @param temp The MC temperature
63 @param filterbyname Not used
64 @param score_moved If True, attempt to speed up sampling by
65 caching scoring function terms on particles that didn't move
74 self.simulated_annealing =
False
75 self.selfadaptive =
False
91 self.mc.set_scoring_function(get_restraint_set(self.model))
92 self.mc.set_return_best(
False)
93 self.mc.set_score_moved(score_moved)
94 self.mc.set_kt(self.temp)
95 self.mc.add_mover(self.smv)
97 def set_kt(self, temp):
104 def set_scoring_function(self, objectlist):
106 for ob
in objectlist:
107 rs.add_restraint(ob.get_restraint())
109 self.mc.set_scoring_function(sf)
111 def set_simulated_annealing(
117 self.simulated_annealing =
True
118 self.tempmin = min_temp
119 self.tempmax = max_temp
120 self.timemin = min_temp_time
121 self.timemax = max_temp_time
123 def set_self_adaptive(self, isselfadaptive=True):
124 self.selfadaptive = isselfadaptive
128 Return a dictionary with the mover parameters for nuisance parameters
131 for i
in range(self.get_number_of_movers()):
132 mv = self.smv.get_mover(i)
134 if "Nuisances" in name:
135 stepsize = IMP.core.NormalMover.get_from(mv).get_sigma()
136 output[name] = stepsize
139 def get_number_of_movers(self):
140 return len(self.smv.get_movers())
142 def get_particle_types(self):
145 def optimize(self, nstep):
147 self.mc.optimize(nstep * self.get_number_of_movers())
150 if self.simulated_annealing:
151 self.temp = self.temp_simulated_annealing()
152 self.mc.set_kt(self.temp)
155 if self.selfadaptive:
156 for i, mv
in enumerate(self.mvs):
158 mvacc = mv.get_number_of_accepted()
159 mvprp = mv.get_number_of_proposed()
160 if mv
not in self.movers_data:
161 accept = float(mvacc) / float(mvprp)
162 self.movers_data[mv] = (mvacc, mvprp)
164 oldmvacc, oldmvprp = self.movers_data[mv]
165 accept = float(mvacc-oldmvacc) / float(mvprp-oldmvprp)
166 self.movers_data[mv] = (mvacc, mvprp)
173 stepsize = mv.get_sigma()
174 if 0.4 > accept
or accept > 0.6:
175 mv.set_sigma(stepsize * 2 * accept)
178 stepsize = mv.get_radius()
179 if 0.4 > accept
or accept > 0.6:
180 mv.set_radius(stepsize * 2 * accept)
183 mr = mv.get_maximum_rotation()
184 mt = mv.get_maximum_translation()
185 if 0.4 > accept
or accept > 0.6:
186 mv.set_maximum_rotation(mr * 2 * accept)
187 mv.set_maximum_translation(mt * 2 * accept)
190 mr = mv.get_maximum_rotation()
191 mt = mv.get_maximum_translation()
192 if 0.4 > accept
or accept > 0.6:
193 mv.set_maximum_rotation(mr * 2 * accept)
194 mv.set_maximum_translation(mt * 2 * accept)
198 if 0.4 > accept
or accept > 0.6:
199 mv.set_radius(mr * 2 * accept)
201 def get_nuisance_movers(self, nuisances, maxstep):
203 for nuisance
in nuisances:
204 print(nuisance, maxstep)
211 def get_rigid_body_movers(self, rbs, maxtrans, maxrot):
218 def get_super_rigid_body_movers(self, rbs, maxtrans, maxrot):
225 if isinstance(rb[2], tuple)
and len(rb[2]) == 3 \
226 and isinstance(rb[2][0], float) \
227 and isinstance(rb[2][1], float) \
228 and isinstance(rb[2][2], float):
235 "Setting up a super rigid body with wrong parameters")
239 srbm.add_xyz_particle(xyz)
241 srbm.add_rigid_body_particle(rb)
245 def get_floppy_body_movers(self, fbs, maxtrans):
254 fb.set_is_optimized(fk,
True)
264 def get_X_movers(self, fbs, maxtrans):
270 raise ValueError(
"particle is part of a rigid body")
276 def get_weight_movers(self, weights, maxstep):
278 for weight
in weights:
279 if weight.get_number_of_weights() > 1:
283 def get_surface_movers(self, surfaces, maxtrans, maxrot, refprob):
285 for surface
in surfaces:
290 def temp_simulated_annealing(self):
291 if self.nframe % (self.timemin + self.timemax) < self.timemin:
295 temp = self.tempmin + (self.tempmax - self.tempmin) * value
298 def set_label(self, label):
301 def get_frame_number(self):
304 def get_output(self):
306 for i, mv
in enumerate(self.smv.get_movers()):
307 mvname = mv.get_name()
308 mvacc = mv.get_number_of_accepted()
309 mvprp = mv.get_number_of_proposed()
311 mvacr = float(mvacc) / float(mvprp)
314 output[
"MonteCarlo_Acceptance_" +
315 mvname +
"_" + str(i)] = str(mvacr)
316 if "Nuisances" in mvname:
317 output[
"MonteCarlo_StepSize_" + mvname +
"_" + str(i)] = \
318 str(IMP.core.NormalMover.get_from(mv).get_sigma())
319 if "Weights" in mvname:
320 output[
"MonteCarlo_StepSize_" + mvname +
"_" + str(i)] = \
321 str(IMP.isd.WeightMover.get_from(mv).get_radius())
322 output[
"MonteCarlo_Temperature"] = str(self.mc.get_kt())
323 output[
"MonteCarlo_Nframe"] = str(self.nframe)
328 """Sample using molecular dynamics"""
330 def __init__(self, model, objects, kt, gamma=0.01, maximum_time_step=1.0,
333 @param model The IMP Model
334 @param objects What to sample. Use flat list of particles
335 @param kt Temperature
336 @param gamma Viscosity parameter
337 @param maximum_time_step MD max time step
344 psamp = obj.get_particles_to_sample()
345 to_sample = psamp[
'Floppy_Bodies_SimplifiedModel'][0]
350 self.model, to_sample, kt/0.0019872041, gamma)
352 self.md.set_maximum_time_step(maximum_time_step)
354 self.md.set_scoring_function(sf)
356 self.md.set_scoring_function(get_restraint_set(self.model))
357 self.md.add_optimizer_state(self.ltstate)
358 self.simulated_annealing =
False
361 def set_kt(self, kt):
362 temp = kt/0.0019872041
363 self.ltstate.set_temperature(temp)
364 self.md.assign_velocities(temp)
366 def set_simulated_annealing(self, min_temp, max_temp, min_temp_time,
368 self.simulated_annealing =
True
369 self.tempmin = min_temp
370 self.tempmax = max_temp
371 self.timemin = min_temp_time
372 self.timemax = max_temp_time
374 def temp_simulated_annealing(self):
375 if self.nframe % (self.timemin + self.timemax) < self.timemin:
379 temp = self.tempmin + (self.tempmax - self.tempmin) * value
382 def set_gamma(self, gamma):
383 self.ltstate.set_gamma(gamma)
385 def optimize(self, nsteps):
388 if self.simulated_annealing:
389 self.temp = self.temp_simulated_annealing()
390 self.set_kt(self.temp)
391 self.md.optimize(nsteps)
393 def get_output(self):
395 output[
"MolecularDynamics_KineticEnergy"] = \
396 str(self.md.get_kinetic_energy())
401 """Sample using conjugate gradients"""
403 def __init__(self, model, objects):
407 self.cg.set_scoring_function(get_restraint_set(self.model))
409 def set_label(self, label):
412 def get_frame_number(self):
415 def optimize(self, nstep):
417 self.cg.optimize(nstep)
419 def set_scoring_function(self, objectlist):
421 for ob
in objectlist:
422 rs.add_restraint(ob.get_restraint())
424 self.cg.set_scoring_function(sf)
426 def get_output(self):
428 output[
"ConjugatedGradients_Nframe"] = str(self.nframe)
433 """Sample using replica exchange"""
435 def __init__(self, model, tempmin, tempmax, samplerobjects, test=True,
436 replica_exchange_object=
None):
438 samplerobjects can be a list of MonteCarlo or MolecularDynamics
442 self.samplerobjects = samplerobjects
444 self.TEMPMIN_ = tempmin
445 self.TEMPMAX_ = tempmax
447 if replica_exchange_object
is None:
451 print(
'ReplicaExchange: MPI was found. '
452 'Using Parallel Replica Exchange')
455 print(
'ReplicaExchange: Could not find MPI. '
456 'Using Serial Replica Exchange')
457 self.rem = _SerialReplicaExchange()
461 print(
'got existing rex object')
462 self.rem = replica_exchange_object
465 nproc = self.rem.get_number_of_replicas()
467 if nproc % 2 != 0
and not test:
469 "number of replicas has to be even. "
470 "set test=True to run with odd number of replicas.")
472 temp = self.rem.create_temperatures(
477 self.temperatures = temp
479 myindex = self.rem.get_my_index()
481 self.rem.set_my_parameter(
"temp", [self.temperatures[myindex]])
482 for so
in self.samplerobjects:
483 so.set_kt(self.temperatures[myindex])
489 def get_temperatures(self):
490 return self.temperatures
492 def get_my_temp(self):
493 return self.rem.get_my_parameter(
"temp")[0]
495 def get_my_index(self):
496 return self.rem.get_my_index()
498 def swap_temp(self, nframe, score=None):
500 score = self.model.evaluate(
False)
502 _ = self.rem.get_my_index()
503 mytemp = self.rem.get_my_parameter(
"temp")[0]
505 if mytemp == self.TEMPMIN_:
508 if mytemp == self.TEMPMAX_:
512 myscore = score / mytemp
515 findex = self.rem.get_friend_index(nframe)
516 ftemp = self.rem.get_friend_parameter(
"temp", findex)[0]
518 fscore = score / ftemp
521 flag = self.rem.do_exchange(myscore, fscore, findex)
526 for so
in self.samplerobjects:
530 def get_output(self):
532 if self.nattempts != 0:
533 output[
"ReplicaExchange_SwapSuccessRatio"] = str(
534 float(self.nsuccess) / self.nattempts)
535 output[
"ReplicaExchange_MinTempFrequency"] = str(
536 float(self.nmintemp) / self.nattempts)
537 output[
"ReplicaExchange_MaxTempFrequency"] = str(
538 float(self.nmaxtemp) / self.nattempts)
540 output[
"ReplicaExchange_SwapSuccessRatio"] = str(0)
541 output[
"ReplicaExchange_MinTempFrequency"] = str(0)
542 output[
"ReplicaExchange_MaxTempFrequency"] = str(0)
543 output[
"ReplicaExchange_CurrentTemp"] = str(self.get_my_temp())
547 class MPI_values(object):
548 def __init__(self, replica_exchange_object=None):
549 """Query values (ie score, and others)
550 from a set of parallel jobs"""
552 if replica_exchange_object
is None:
556 print(
'MPI_values: MPI was found. '
557 'Using Parallel Replica Exchange')
560 print(
'MPI_values: Could not find MPI. '
561 'Using Serial Replica Exchange')
562 self.rem = _SerialReplicaExchange()
566 print(
'got existing rex object')
567 self.rem = replica_exchange_object
569 def set_value(self, name, value):
570 self.rem.set_my_parameter(name, [value])
572 def get_values(self, name):
574 for i
in range(self.rem.get_number_of_replicas()):
575 v = self.rem.get_friend_parameter(name, i)[0]
579 def get_percentile(self, name):
580 value = self.rem.get_my_parameter(name)[0]
581 values = sorted(self.get_values(name))
582 ind = values.index(value)
583 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.
Create a scoring function on a list of restraints.
Move continuous particle variables by perturbing them within a ball.
Modify a surface orientation.
static FloatKeys get_internal_coordinate_keys()
Object used to hold a set of restraints.
Simple molecular dynamics simulator.
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 ...