IMP logo
IMP Reference Guide  2.18.0
The Integrative Modeling Platform
samplers.py
1 """@namespace IMP.pmi.samplers
2  Sampling of the system.
3 """
4 
5 from __future__ import print_function
6 import IMP
7 import IMP.core
8 from IMP.pmi.tools import get_restraint_set
9 
10 
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.
15  """
16  def __init__(self):
17  self.__params = {}
18 
19  def get_number_of_replicas(self):
20  return 1
21 
22  def create_temperatures(self, tmin, tmax, nrep):
23  return [tmin]
24 
25  def get_my_index(self):
26  return 0
27 
28  def set_my_parameter(self, key, val):
29  self.__params[key] = val
30 
31  def get_my_parameter(self, key):
32  return self.__params[key]
33 
34  def get_friend_index(self, step):
35  return 0
36 
37  def get_friend_parameter(self, key, findex):
38  return self.get_my_parameter(key)
39 
40  def do_exchange(self, myscore, fscore, findex):
41  return False
42 
43  def set_was_used(self, was_used):
44  self.was_used = was_used
45 
46 
47 class MonteCarlo(object):
48  """Sample using Monte Carlo"""
49 
50  # check that isd is installed
51  try:
52  import IMP.isd
53  isd_available = True
54  except ImportError:
55  isd_available = False
56 
57  def __init__(self, model, objects=None, temp=1.0, filterbyname=None,
58  score_moved=False):
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
66  """
67  self.losp = [
68  "Rigid_Bodies",
69  "Floppy_Bodies",
70  "Nuisances",
71  "X_coord",
72  "Weights"
73  "Surfaces"]
74  self.simulated_annealing = False
75  self.selfadaptive = False
76  # that is -1 because mc has not yet run
77  self.nframe = -1
78  self.temp = temp
79  self.mvs = []
80  self.mvslabels = []
81  self.label = "None"
82  self.model = model
83  self.movers_data = {}
84 
85  # check if using PMI1 or just passed a list of movers
86  gather_objects = False
87  try:
88  objects[0].get_particles_to_sample()
89  gather_objects = True
90  except AttributeError:
91  self.mvs = objects
92 
93  if gather_objects:
95  "Pass Movers to MonteCarlo(), not objects that implement "
96  "get_particles_to_sample()")
97  for ob in objects:
98  pts = ob.get_particles_to_sample()
99  for k in pts.keys():
100 
101  if "Rigid_Bodies" in k:
102  mvs = self.get_rigid_body_movers(
103  pts[k][0],
104  pts[k][1],
105  pts[k][2])
106  for mv in mvs:
107  mv.set_name(k)
108  self.mvs += mvs
109 
110  if "SR_Bodies" in k:
111  print(len(pts[k]))
112  mvs = self.get_super_rigid_body_movers(
113  pts[k][0],
114  pts[k][1],
115  pts[k][2])
116  for mv in mvs:
117  mv.set_name(k)
118  self.mvs += mvs
119 
120  if "Floppy_Bodies" in k:
121  mvs = self.get_floppy_body_movers(pts[k][0], pts[k][1])
122  for mv in mvs:
123  mv.set_name(k)
124  self.mvs += mvs
125 
126  if "X_coord" in k:
127  mvs = self.get_X_movers(pts[k][0], pts[k][1])
128  for mv in mvs:
129  mv.set_name(k)
130  self.mvs += mvs
131 
132  if "Nuisances" in k:
133  if not self.isd_available:
134  raise ValueError(
135  "isd module needed to use nuisances")
136  mvs = self.get_nuisance_movers(pts[k][0], pts[k][1])
137  for mv in mvs:
138  mv.set_name(k)
139  self.mvs += mvs
140 
141  if "Weights" in k:
142  if not self.isd_available:
143  raise ValueError(
144  "isd module needed to use weights")
145  mvs = self.get_weight_movers(pts[k][0], pts[k][1])
146  for mv in mvs:
147  mv.set_name(k)
148  self.mvs += mvs
149 
150  if "Surfaces" in k:
151  mvs = self.get_surface_movers(
152  pts[k][0],
153  pts[k][1],
154  pts[k][2],
155  pts[k][3])
156  for mv in mvs:
157  mv.set_name(k)
158  self.mvs += mvs
159 
160  # SerialMover
161  self.smv = IMP.core.SerialMover(self.mvs)
162 
163  self.mc = IMP.core.MonteCarlo(self.model)
164  self.mc.set_scoring_function(get_restraint_set(self.model))
165  self.mc.set_return_best(False)
166  self.mc.set_score_moved(score_moved)
167  self.mc.set_kt(self.temp)
168  self.mc.add_mover(self.smv)
169 
170  def set_kt(self, temp):
171  self.temp = temp
172  self.mc.set_kt(temp)
173 
174  def get_mc(self):
175  return self.mc
176 
177  def set_scoring_function(self, objectlist):
178  rs = IMP.RestraintSet(self.model, 1.0, 'sfo')
179  for ob in objectlist:
180  rs.add_restraint(ob.get_restraint())
182  self.mc.set_scoring_function(sf)
183 
184  def set_simulated_annealing(
185  self,
186  min_temp,
187  max_temp,
188  min_temp_time,
189  max_temp_time):
190  self.simulated_annealing = True
191  self.tempmin = min_temp
192  self.tempmax = max_temp
193  self.timemin = min_temp_time
194  self.timemax = max_temp_time
195 
196  def set_self_adaptive(self, isselfadaptive=True):
197  self.selfadaptive = isselfadaptive
198 
200  '''
201  Return a dictionary with the mover parameters for nuisance parameters
202  '''
203  output = {}
204  for i in range(self.get_number_of_movers()):
205  mv = self.smv.get_mover(i)
206  name = mv.get_name()
207  if "Nuisances" in name:
208  stepsize = IMP.core.NormalMover.get_from(mv).get_sigma()
209  output[name] = stepsize
210  return output
211 
212  def get_number_of_movers(self):
213  return len(self.smv.get_movers())
214 
215  def get_particle_types(self):
216  return self.losp
217 
218  def optimize(self, nstep):
219  self.nframe += 1
220  self.mc.optimize(nstep * self.get_number_of_movers())
221 
222  # apply simulated annealing protocol
223  if self.simulated_annealing:
224  self.temp = self.temp_simulated_annealing()
225  self.mc.set_kt(self.temp)
226 
227  # apply self adaptive protocol
228  if self.selfadaptive:
229  for i, mv in enumerate(self.mvs):
230 
231  mvacc = mv.get_number_of_accepted()
232  mvprp = mv.get_number_of_proposed()
233  if mv not in self.movers_data:
234  accept = float(mvacc) / float(mvprp)
235  self.movers_data[mv] = (mvacc, mvprp)
236  else:
237  oldmvacc, oldmvprp = self.movers_data[mv]
238  accept = float(mvacc-oldmvacc) / float(mvprp-oldmvprp)
239  self.movers_data[mv] = (mvacc, mvprp)
240  if accept < 0.05:
241  accept = 0.05
242  if accept > 1.0:
243  accept = 1.0
244 
245  if isinstance(mv, IMP.core.NormalMover):
246  stepsize = mv.get_sigma()
247  if 0.4 > accept or accept > 0.6:
248  mv.set_sigma(stepsize * 2 * accept)
249 
250  if isinstance(mv, IMP.isd.WeightMover):
251  stepsize = mv.get_radius()
252  if 0.4 > accept or accept > 0.6:
253  mv.set_radius(stepsize * 2 * accept)
254 
255  if isinstance(mv, IMP.core.RigidBodyMover):
256  mr = mv.get_maximum_rotation()
257  mt = mv.get_maximum_translation()
258  if 0.4 > accept or accept > 0.6:
259  mv.set_maximum_rotation(mr * 2 * accept)
260  mv.set_maximum_translation(mt * 2 * accept)
261 
262  if isinstance(mv, IMP.pmi.TransformMover):
263  mr = mv.get_maximum_rotation()
264  mt = mv.get_maximum_translation()
265  if 0.4 > accept or accept > 0.6:
266  mv.set_maximum_rotation(mr * 2 * accept)
267  mv.set_maximum_translation(mt * 2 * accept)
268 
269  if isinstance(mv, IMP.core.BallMover):
270  mr = mv.get_radius()
271  if 0.4 > accept or accept > 0.6:
272  mv.set_radius(mr * 2 * accept)
273 
274  def get_nuisance_movers(self, nuisances, maxstep):
275  mvs = []
276  for nuisance in nuisances:
277  print(nuisance, maxstep)
278  mvs.append(
279  IMP.core.NormalMover([nuisance],
280  IMP.FloatKeys([IMP.FloatKey("nuisance")]),
281  maxstep))
282  return mvs
283 
284  def get_rigid_body_movers(self, rbs, maxtrans, maxrot):
285  mvs = []
286  for rb in rbs:
287  mvs.append(IMP.core.RigidBodyMover(rb.get_model(), rb,
288  maxtrans, maxrot))
289  return mvs
290 
291  def get_super_rigid_body_movers(self, rbs, maxtrans, maxrot):
292  mvs = []
293  for rb in rbs:
294  if len(rb) == 2:
295  # normal Super Rigid Body
296  srbm = IMP.pmi.TransformMover(self.model, maxtrans, maxrot)
297  elif len(rb) == 3:
298  if isinstance(rb[2], tuple) and len(rb[2]) == 3 \
299  and isinstance(rb[2][0], float) \
300  and isinstance(rb[2][1], float) \
301  and isinstance(rb[2][2], float):
302  # super rigid body with 2D rotation, rb[2] is the axis
303  srbm = IMP.pmi.TransformMover(
304  self.model, IMP.algebra.Vector3D(rb[2]), maxtrans,
305  maxrot)
306  else:
307  print(
308  "Setting up a super rigid body with wrong parameters")
309  raise
310 
311  for xyz in rb[0]:
312  srbm.add_xyz_particle(xyz)
313  for rb in rb[1]:
314  srbm.add_rigid_body_particle(rb)
315  mvs.append(srbm)
316  return mvs
317 
318  def get_floppy_body_movers(self, fbs, maxtrans):
319  mvs = []
320  for fb in fbs:
321  # check is that is a rigid body member:
323  # if so force the particles to move anyway
324  floatkeys = \
326  for fk in floatkeys:
327  fb.set_is_optimized(fk, True)
328  mvs.append(
329  IMP.core.BallMover(fb.get_model(), fb,
330  IMP.FloatKeys(floatkeys),
331  maxtrans))
332  else:
333  # otherwise use the normal ball mover
334  mvs.append(IMP.core.BallMover(fb.get_model(), fb, maxtrans))
335  return mvs
336 
337  def get_X_movers(self, fbs, maxtrans):
338  mvs = []
339  Xfloatkey = IMP.core.XYZ.get_xyz_keys()[0]
340  for fb in fbs:
341  # check is that is a rigid body member:
343  raise ValueError("particle is part of a rigid body")
344  else:
345  # otherwise use the normal ball mover
346  mvs.append(IMP.core.NormalMover([fb], [Xfloatkey], maxtrans))
347  return mvs
348 
349  def get_weight_movers(self, weights, maxstep):
350  mvs = []
351  for weight in weights:
352  if weight.get_number_of_weights() > 1:
353  mvs.append(IMP.isd.WeightMover(weight, maxstep))
354  return mvs
355 
356  def get_surface_movers(self, surfaces, maxtrans, maxrot, refprob):
357  mvs = []
358  for surface in surfaces:
359  mvs.append(IMP.core.SurfaceMover(surface, maxtrans, maxrot,
360  refprob))
361  return mvs
362 
363  def temp_simulated_annealing(self):
364  if self.nframe % (self.timemin + self.timemax) < self.timemin:
365  value = 0.0
366  else:
367  value = 1.0
368  temp = self.tempmin + (self.tempmax - self.tempmin) * value
369  return temp
370 
371  def set_label(self, label):
372  self.label = label
373 
374  def get_frame_number(self):
375  return self.nframe
376 
377  def get_output(self):
378  output = {}
379  for i, mv in enumerate(self.smv.get_movers()):
380  mvname = mv.get_name()
381  mvacc = mv.get_number_of_accepted()
382  mvprp = mv.get_number_of_proposed()
383  try:
384  mvacr = float(mvacc) / float(mvprp)
385  except: # noqa: E722
386  mvacr = 0.0
387  output["MonteCarlo_Acceptance_" +
388  mvname + "_" + str(i)] = str(mvacr)
389  if "Nuisances" in mvname:
390  output["MonteCarlo_StepSize_" + mvname + "_" + str(i)] = \
391  str(IMP.core.NormalMover.get_from(mv).get_sigma())
392  if "Weights" in mvname:
393  output["MonteCarlo_StepSize_" + mvname + "_" + str(i)] = \
394  str(IMP.isd.WeightMover.get_from(mv).get_radius())
395  output["MonteCarlo_Temperature"] = str(self.mc.get_kt())
396  output["MonteCarlo_Nframe"] = str(self.nframe)
397  return output
398 
399 
400 class MolecularDynamics(object):
401  """Sample using molecular dynamics"""
402 
403  def __init__(self, model, objects, kt, gamma=0.01, maximum_time_step=1.0,
404  sf=None):
405  """Setup MD
406  @param model The IMP Model
407  @param objects What to sample. Use flat list of particles
408  @param kt Temperature
409  @param gamma Viscosity parameter
410  @param maximum_time_step MD max time step
411  """
412  self.model = model
413 
414  # check if using PMI1 objects dictionary, or just list of particles
415  try:
416  for obj in objects:
417  psamp = obj.get_particles_to_sample()
418  to_sample = psamp['Floppy_Bodies_SimplifiedModel'][0]
419  except: # noqa: E722
420  to_sample = objects
421 
423  self.model, to_sample, kt/0.0019872041, gamma)
424  self.md = IMP.atom.MolecularDynamics(self.model)
425  self.md.set_maximum_time_step(maximum_time_step)
426  if sf:
427  self.md.set_scoring_function(sf)
428  else:
429  self.md.set_scoring_function(get_restraint_set(self.model))
430  self.md.add_optimizer_state(self.ltstate)
431  self.simulated_annealing = False
432  self.nframe = -1
433 
434  def set_kt(self, kt):
435  temp = kt/0.0019872041
436  self.ltstate.set_temperature(temp)
437  self.md.assign_velocities(temp)
438 
439  def set_simulated_annealing(self, min_temp, max_temp, min_temp_time,
440  max_temp_time):
441  self.simulated_annealing = True
442  self.tempmin = min_temp
443  self.tempmax = max_temp
444  self.timemin = min_temp_time
445  self.timemax = max_temp_time
446 
447  def temp_simulated_annealing(self):
448  if self.nframe % (self.timemin + self.timemax) < self.timemin:
449  value = 0.0
450  else:
451  value = 1.0
452  temp = self.tempmin + (self.tempmax - self.tempmin) * value
453  return temp
454 
455  def set_gamma(self, gamma):
456  self.ltstate.set_gamma(gamma)
457 
458  def optimize(self, nsteps):
459  # apply simulated annealing protocol
460  self.nframe += 1
461  if self.simulated_annealing:
462  self.temp = self.temp_simulated_annealing()
463  self.set_kt(self.temp)
464  self.md.optimize(nsteps)
465 
466  def get_output(self):
467  output = {}
468  output["MolecularDynamics_KineticEnergy"] = \
469  str(self.md.get_kinetic_energy())
470  return output
471 
472 
473 class ConjugateGradients(object):
474  """Sample using conjugate gradients"""
475 
476  def __init__(self, model, objects):
477  self.model = model
478  self.nframe = -1
479  self.cg = IMP.core.ConjugateGradients(self.model)
480  self.cg.set_scoring_function(get_restraint_set(self.model))
481 
482  def set_label(self, label):
483  self.label = label
484 
485  def get_frame_number(self):
486  return self.nframe
487 
488  def optimize(self, nstep):
489  self.nframe += 1
490  self.cg.optimize(nstep)
491 
492  def set_scoring_function(self, objectlist):
493  rs = IMP.RestraintSet(self.model, 1.0, 'sfo')
494  for ob in objectlist:
495  rs.add_restraint(ob.get_restraint())
497  self.cg.set_scoring_function(sf)
498 
499  def get_output(self):
500  output = {}
501  output["ConjugatedGradients_Nframe"] = str(self.nframe)
502  return output
503 
504 
505 class ReplicaExchange(object):
506  """Sample using replica exchange"""
507 
508  def __init__(self, model, tempmin, tempmax, samplerobjects, test=True,
509  replica_exchange_object=None):
510  '''
511  samplerobjects can be a list of MonteCarlo or MolecularDynamics
512  '''
513 
514  self.model = model
515  self.samplerobjects = samplerobjects
516  # min and max temperature
517  self.TEMPMIN_ = tempmin
518  self.TEMPMAX_ = tempmax
519 
520  if replica_exchange_object is None:
521  # initialize Replica Exchange class
522  try:
523  import IMP.mpi
524  print('ReplicaExchange: MPI was found. '
525  'Using Parallel Replica Exchange')
526  self.rem = IMP.mpi.ReplicaExchange()
527  except ImportError:
528  print('ReplicaExchange: Could not find MPI. '
529  'Using Serial Replica Exchange')
530  self.rem = _SerialReplicaExchange()
531 
532  else:
533  # get the replica exchange class instance from elsewhere
534  print('got existing rex object')
535  self.rem = replica_exchange_object
536 
537  # get number of replicas
538  nproc = self.rem.get_number_of_replicas()
539 
540  if nproc % 2 != 0 and not test:
541  raise Exception(
542  "number of replicas has to be even. "
543  "set test=True to run with odd number of replicas.")
544  # create array of temperatures, in geometric progression
545  temp = self.rem.create_temperatures(
546  self.TEMPMIN_,
547  self.TEMPMAX_,
548  nproc)
549  # get replica index
550  self.temperatures = temp
551 
552  myindex = self.rem.get_my_index()
553  # set initial value of the parameter (temperature) to exchange
554  self.rem.set_my_parameter("temp", [self.temperatures[myindex]])
555  for so in self.samplerobjects:
556  so.set_kt(self.temperatures[myindex])
557  self.nattempts = 0
558  self.nmintemp = 0
559  self.nmaxtemp = 0
560  self.nsuccess = 0
561 
562  def get_temperatures(self):
563  return self.temperatures
564 
565  def get_my_temp(self):
566  return self.rem.get_my_parameter("temp")[0]
567 
568  def get_my_index(self):
569  return self.rem.get_my_index()
570 
571  def swap_temp(self, nframe, score=None):
572  if score is None:
573  score = self.model.evaluate(False)
574  # get my replica index and temperature
575  _ = self.rem.get_my_index()
576  mytemp = self.rem.get_my_parameter("temp")[0]
577 
578  if mytemp == self.TEMPMIN_:
579  self.nmintemp += 1
580 
581  if mytemp == self.TEMPMAX_:
582  self.nmaxtemp += 1
583 
584  # score divided by kbt
585  myscore = score / mytemp
586 
587  # get my friend index and temperature
588  findex = self.rem.get_friend_index(nframe)
589  ftemp = self.rem.get_friend_parameter("temp", findex)[0]
590  # score divided by kbt
591  fscore = score / ftemp
592 
593  # try exchange
594  flag = self.rem.do_exchange(myscore, fscore, findex)
595 
596  self.nattempts += 1
597  # if accepted, change temperature
598  if (flag):
599  for so in self.samplerobjects:
600  so.set_kt(ftemp)
601  self.nsuccess += 1
602 
603  def get_output(self):
604  output = {}
605  if self.nattempts != 0:
606  output["ReplicaExchange_SwapSuccessRatio"] = str(
607  float(self.nsuccess) / self.nattempts)
608  output["ReplicaExchange_MinTempFrequency"] = str(
609  float(self.nmintemp) / self.nattempts)
610  output["ReplicaExchange_MaxTempFrequency"] = str(
611  float(self.nmaxtemp) / self.nattempts)
612  else:
613  output["ReplicaExchange_SwapSuccessRatio"] = str(0)
614  output["ReplicaExchange_MinTempFrequency"] = str(0)
615  output["ReplicaExchange_MaxTempFrequency"] = str(0)
616  output["ReplicaExchange_CurrentTemp"] = str(self.get_my_temp())
617  return output
618 
619 
620 class MPI_values(object):
621  def __init__(self, replica_exchange_object=None):
622  """Query values (ie score, and others)
623  from a set of parallel jobs"""
624 
625  if replica_exchange_object is None:
626  # initialize Replica Exchange class
627  try:
628  import IMP.mpi
629  print('MPI_values: MPI was found. '
630  'Using Parallel Replica Exchange')
631  self.rem = IMP.mpi.ReplicaExchange()
632  except ImportError:
633  print('MPI_values: Could not find MPI. '
634  'Using Serial Replica Exchange')
635  self.rem = _SerialReplicaExchange()
636 
637  else:
638  # get the replica exchange class instance from elsewhere
639  print('got existing rex object')
640  self.rem = replica_exchange_object
641 
642  def set_value(self, name, value):
643  self.rem.set_my_parameter(name, [value])
644 
645  def get_values(self, name):
646  values = []
647  for i in range(self.rem.get_number_of_replicas()):
648  v = self.rem.get_friend_parameter(name, i)[0]
649  values.append(v)
650  return values
651 
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)
657  return percentile
def __init__
samplerobjects can be a list of MonteCarlo or MolecularDynamics
Definition: samplers.py:508
A Monte Carlo optimizer.
Definition: MonteCarlo.h:45
A class to implement Hamiltonian Replica Exchange.
Maintains temperature during molecular dynamics.
Sample using molecular dynamics.
Definition: samplers.py:400
def get_nuisance_movers_parameters
Return a dictionary with the mover parameters for nuisance parameters.
Definition: samplers.py:199
Miscellaneous utilities.
Definition: tools.py:1
void handle_use_deprecated(std::string message)
Break in this method in gdb to find deprecated uses at runtime.
Modify the transformation of a rigid body.
Simple conjugate gradients optimizer.
Sample using conjugate gradients.
Definition: samplers.py:473
Create a scoring function on a list of restraints.
Move continuous particle variables by perturbing them within a ball.
Modify a surface orientation.
Definition: SurfaceMover.h:34
static FloatKeys get_internal_coordinate_keys()
Definition: rigid_bodies.h:709
Object used to hold a set of restraints.
Definition: RestraintSet.h:37
Simple molecular dynamics optimizer.
Code that uses the MPI parallel library.
A mover that perturbs a Weight particle.
Definition: WeightMover.h:20
Modify the transformation of a rigid body.
def __init__
Setup Monte Carlo sampling.
Definition: samplers.py:57
Modify a set of continuous variables using a normal distribution.
Definition: NormalMover.h:20
Basic functionality that is expected to be used by a wide variety of IMP users.
Sample using Monte Carlo.
Definition: samplers.py:47
The general base class for IMP exceptions.
Definition: exception.h:48
VectorD< 3 > Vector3D
Definition: VectorD.h:421
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.
Definition: SerialMover.h:23
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:752
Sample using replica exchange.
Definition: samplers.py:505
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...