IMP logo
IMP Reference Guide  2.6.1
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 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.
13  """
14  def __init__(self):
15  self.__params = {}
16  def get_number_of_replicas(self):
17  return 1
18  def create_temperatures(self, tmin, tmax, nrep):
19  return [tmin]
20  def get_my_index(self):
21  return 0
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):
27  return 0
28  def get_friend_parameter(self, key, findex):
29  return self.get_my_parameter(key)
30  def do_exchange(self, myscore, fscore, findex):
31  return False
32  def set_was_used(self,was_used):
33  self.was_used = was_used
34 
35 
36 class MonteCarlo(object):
37  """Sample using Monte Carlo"""
38 
39  # check that isd is installed
40  try:
41  import IMP.isd
42  isd_available = True
43  except ImportError:
44  isd_available = False
45 
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
53  """
54  self.losp = [
55  "Rigid_Bodies",
56  "Floppy_Bodies",
57  "Nuisances",
58  "X_coord",
59  "Weights"]
60  self.simulated_annealing = False
61  self.selfadaptive = False
62  # that is -1 because mc has not yet run
63  self.nframe = -1
64  self.temp = temp
65  self.mvs = []
66  self.mvslabels = []
67  self.label = "None"
68  self.m = m
69 
70  # check if using PMI1 or just passed a list of movers
71  gather_objects = False
72  try:
73  objects[0].get_particles_to_sample()
74  gather_objects = True
75  except:
76  self.mvs = objects
77 
78  if gather_objects:
79  for ob in objects:
80  pts = ob.get_particles_to_sample()
81  for k in pts.keys():
82 
83  if "Rigid_Bodies" in k:
84  mvs = self.get_rigid_body_movers(
85  pts[k][0],
86  pts[k][1],
87  pts[k][2])
88  for mv in mvs:
89  mv.set_name(k)
90  self.mvs += mvs
91 
92  if "SR_Bodies" in k:
93  print(len(pts[k]))
94  mvs = self.get_super_rigid_body_movers(
95  pts[k][0],
96  pts[k][1],
97  pts[k][2])
98  for mv in mvs:
99  mv.set_name(k)
100  self.mvs += mvs
101 
102  if "Floppy_Bodies" in k:
103  mvs = self.get_floppy_body_movers(pts[k][0], pts[k][1])
104  for mv in mvs:
105  mv.set_name(k)
106  self.mvs += mvs
107 
108  if "X_coord" in k:
109  mvs = self.get_X_movers(pts[k][0], pts[k][1])
110  for mv in mvs:
111  mv.set_name(k)
112  self.mvs += mvs
113 
114  if "Nuisances" in k:
115  if not self.isd_available:
116  raise ValueError("isd module needed to use nuisances")
117  mvs = self.get_nuisance_movers(pts[k][0], pts[k][1])
118  for mv in mvs:
119  mv.set_name(k)
120  self.mvs += mvs
121 
122  if "Weights" in k:
123  if not self.isd_available:
124  raise ValueError("isd module needed to use weights")
125  mvs = self.get_weight_movers(pts[k][0], pts[k][1])
126  for mv in mvs:
127  mv.set_name(k)
128  self.mvs += mvs
129 
130  # SerialMover
131  self.smv = IMP.core.SerialMover(self.mvs)
132 
133  self.mc = IMP.core.MonteCarlo(self.m)
134  self.mc.set_scoring_function(get_restraint_set(self.m))
135  self.mc.set_return_best(False)
136  self.mc.set_kt(self.temp)
137  self.mc.add_mover(self.smv)
138 
139  def set_kt(self, temp):
140  self.temp = temp
141  self.mc.set_kt(temp)
142 
143  def get_mc(self):
144  return self.mc
145 
146  def set_scoring_function(self, objectlist):
147  rs = IMP.RestraintSet(self.m, 1.0, 'sfo')
148  for ob in objectlist:
149  rs.add_restraint(ob.get_restraint())
151  self.mc.set_scoring_function(sf)
152 
153  def set_simulated_annealing(
154  self,
155  min_temp,
156  max_temp,
157  min_temp_time,
158  max_temp_time):
159  self.simulated_annealing = True
160  self.tempmin = min_temp
161  self.tempmax = max_temp
162  self.timemin = min_temp_time
163  self.timemax = max_temp_time
164 
165  def set_self_adaptive(self, isselfadaptive=True):
166  self.selfadaptive = isselfadaptive
167 
169  '''
170  Return a dictionary with the mover parameters for nuisance parameters
171  '''
172  output = {}
173  for i in range(self.get_number_of_movers()):
174  mv = self.smv.get_mover(i)
175  name = mv.get_name()
176  if "Nuisances" in name:
177  stepsize = IMP.core.NormalMover.get_from(mv).get_sigma()
178  output[name] = stepsize
179  return output
180 
181  def get_number_of_movers(self):
182  return len(self.smv.get_movers())
183 
184  def get_particle_types():
185  return self.losp
186 
187  def optimize(self, nstep):
188  self.nframe += 1
189  self.mc.optimize(nstep * self.get_number_of_movers())
190 
191  # apply simulated annealing protocol
192  if self.simulated_annealing:
193  self.temp = self.temp_simulated_annealing()
194  self.mc.set_kt(self.temp)
195 
196  # apply self adaptive protocol
197  if self.selfadaptive:
198  for i, mv in enumerate(self.smv.get_movers()):
199  name = mv.get_name()
200 
201  if "Nuisances" in name:
202  mvacc = mv.get_number_of_accepted()
203  mvprp = mv.get_number_of_proposed()
204  accept = float(mvacc) / float(mvprp)
205  nmv = IMP.core.NormalMover.get_from(mv)
206  stepsize = nmv.get_sigma()
207 
208  if 0.4 > accept or accept > 0.6:
209  nmv.set_sigma(stepsize * 2 * accept)
210  if accept < 0.05:
211  accept = 0.05
212  nmv.set_sigma(stepsize * 2 * accept)
213  if accept > 1.0:
214  accept = 1.0
215  nmv.set_sigma(stepsize * 2 * accept)
216 
217  if "Weights" in name:
218 
219  mvacc = mv.get_number_of_accepted()
220  mvprp = mv.get_number_of_proposed()
221  accept = float(mvacc) / float(mvprp)
222  wmv = IMP.isd.WeightMover.get_from(mv)
223  stepsize = wmv.get_radius()
224 
225  if 0.4 > accept or accept > 0.6:
226  wmv.set_radius(stepsize * 2 * accept)
227  if accept < 0.05:
228  accept = 0.05
229  wmv.set_radius(stepsize * 2 * accept)
230  if accept > 1.0:
231  accept = 1.0
232  wmv.set_radius(stepsize * 2 * accept)
233 
234  @IMP.deprecated_method("2.5", "Use optimize() instead.")
235  def run(self, *args, **kwargs):
236  self.optimize(*args, **kwargs)
237 
238  def get_nuisance_movers(self, nuisances, maxstep):
239  mvs = []
240  for nuisance in nuisances:
241  print(nuisance, maxstep)
242  mvs.append(
243  IMP.core.NormalMover([nuisance],
244  IMP.FloatKeys([IMP.FloatKey("nuisance")]),
245  maxstep))
246  return mvs
247 
248  def get_rigid_body_movers(self, rbs, maxtrans, maxrot):
249  mvs = []
250  for rb in rbs:
251  mvs.append(IMP.core.RigidBodyMover(rb, maxtrans, maxrot))
252  return mvs
253 
254  def get_super_rigid_body_movers(self, rbs, maxtrans, maxrot):
255  mvs = []
256  for rb in rbs:
257  if len(rb) == 2:
258  # normal Super Rigid Body
259  srbm = IMP.pmi.TransformMover(self.m, maxtrans, maxrot)
260  if len(rb) == 3:
261  # super rigid body with 2D rotation, rb[2] is the axis
262  srbm = IMP.pmi.TransformMover(
263  self.m,
264  IMP.algebra.Vector3D(rb[2]),
265  maxtrans,
266  maxrot)
267  for xyz in rb[0]:
268  srbm.add_xyz_particle(xyz)
269  for rb in rb[1]:
270  srbm.add_rigid_body_particle(rb)
271  mvs.append(srbm)
272  return mvs
273 
274  def get_floppy_body_movers(self, fbs, maxtrans):
275  mvs = []
276  for fb in fbs:
277  # check is that is a rigid body member:
279  # if so force the particles to move anyway
280  floatkeys = [IMP.FloatKey(4), IMP.FloatKey(5), IMP.FloatKey(6)]
281  for fk in floatkeys:
282  fb.set_is_optimized(fk, True)
283  mvs.append(
284  IMP.core.BallMover([fb],
285  IMP.FloatKeys(floatkeys),
286  maxtrans))
287  else:
288  # otherwise use the normal ball mover
289  mvs.append(IMP.core.BallMover([fb], maxtrans))
290  return mvs
291 
292  def get_X_movers(self, fbs, maxtrans):
293  mvs = []
294  Xfloatkey = IMP.core.XYZ.get_xyz_keys()[0]
295  for fb in fbs:
296  # check is that is a rigid body member:
298  raise ValueError("particle is part of a rigid body")
299  else:
300  # otherwise use the normal ball mover
301  mvs.append(IMP.core.NormalMover([fb], [Xfloatkey], maxtrans))
302  return mvs
303 
304  def get_weight_movers(self, weights, maxstep):
305  mvs = []
306  for weight in weights:
307  if(weight.get_number_of_states() > 1):
308  mvs.append(IMP.isd.WeightMover(weight, maxstep))
309  return mvs
310 
311  def temp_simulated_annealing(self):
312  if self.nframe % (self.timemin + self.timemax) < self.timemin:
313  value = 0.0
314  else:
315  value = 1.0
316  temp = self.tempmin + (self.tempmax - self.tempmin) * value
317  return temp
318 
319  def set_label(self, label):
320  self.label = label
321 
322  def get_frame_number(self):
323  return self.nframe
324 
325  def get_output(self):
326  output = {}
327  acceptances = []
328  for i, mv in enumerate(self.smv.get_movers()):
329  mvname = mv.get_name()
330  mvacc = mv.get_number_of_accepted()
331  mvprp = mv.get_number_of_proposed()
332  try:
333  mvacr = float(mvacc) / float(mvprp)
334  except:
335  mvacr = 0.0
336  output["MonteCarlo_Acceptance_" +
337  mvname + "_" + str(i)] = str(mvacr)
338  if "Nuisances" in mvname:
339  output["MonteCarlo_StepSize_" + mvname + "_" +
340  str(i)] = str(IMP.core.NormalMover.get_from(mv).get_sigma())
341  if "Weights" in mvname:
342  output["MonteCarlo_StepSize_" + mvname + "_" +
343  str(i)] = str(IMP.isd.WeightMover.get_from(mv).get_radius())
344  output["MonteCarlo_Temperature"] = str(self.mc.get_kt())
345  output["MonteCarlo_Nframe"] = str(self.nframe)
346  return output
347 
348 
349 class MolecularDynamics(object):
350  """Sample using molecular dynamics"""
351 
352  def __init__(self,m,objects,kt,gamma=0.01,maximum_time_step=1.0):
353  """Setup MD
354  @param m The IMP Model
355  @param objects What to sample. Use flat list of particles or (deprecated) 'MD Sample Objects' from PMI1
356  @param kt Temperature
357  @param gamma Viscosity parameter
358  @param maximum_time_step MD max time step
359  """
360  self.m=m
361 
362  # check if using PMI1 objects dictionary, or just list of particles
363  try:
364  for obj in objects:
365  to_sample=obj.get_particles_to_sample()['Floppy_Bodies_SimplifiedModel'][0]
366  except:
367  to_sample = objects
368 
369  self.ltstate=IMP.atom.LangevinThermostatOptimizerState(self.m,to_sample,
370  kt/0.0019872041,
371  gamma)
372  self.md = IMP.atom.MolecularDynamics(self.m)
373  self.md.set_maximum_time_step(maximum_time_step)
374  self.md.set_scoring_function(get_restraint_set(self.m))
375  self.md.add_optimizer_state(self.ltstate)
376  self.simulated_annealing = False
377  self.nframe = -1
378  def set_kt(self,kt):
379  temp=kt/0.0019872041
380  self.ltstate.set_temperature(temp)
381  self.md.assign_velocities(temp)
382 
383  def set_simulated_annealing(self, min_temp, max_temp, min_temp_time,
384  max_temp_time):
385  self.simulated_annealing = True
386  self.tempmin = min_temp
387  self.tempmax = max_temp
388  self.timemin = min_temp_time
389  self.timemax = max_temp_time
390 
391  def temp_simulated_annealing(self):
392  if self.nframe % (self.timemin + self.timemax) < self.timemin:
393  value = 0.0
394  else:
395  value = 1.0
396  temp = self.tempmin + (self.tempmax - self.tempmin) * value
397  return temp
398 
399  def set_gamma(self,gamma):
400  self.ltstate.set_gamma(gamma)
401 
402  def optimize(self,nsteps):
403  # apply simulated annealing protocol
404  self.nframe+=1
405  if self.simulated_annealing:
406  self.temp = self.temp_simulated_annealing()
407  self.set_kt(self.temp)
408  self.md.optimize(nsteps)
409 
410  def get_output(self):
411  output={}
412  output["MolecularDynamics_KineticEnergy"]=str(self.md.get_kinetic_energy())
413  return output
414 
415 class ConjugateGradients(object):
416  """Sample using conjugate gradients"""
417 
418  def __init__(self, m, objects):
419  self.m = m
420  self.nframe = -1
421  self.cg = IMP.core.ConjugateGradients(self.m)
422  self.cg.set_scoring_function(get_restraint_set(self.m))
423 
424  def set_label(self, label):
425  self.label = label
426 
427  def get_frame_number(self):
428  return self.nframe
429 
430  @IMP.deprecated_method("2.5", "Use optimize() instead.")
431  def run(self, *args, **kwargs):
432  self.optimize(*args, **kwargs)
433 
434  def optimize(self, nstep):
435  self.nframe += 1
436  self.cg.optimize(nstep)
437 
438  def set_scoring_function(self, objectlist):
439  rs = IMP.RestraintSet(self.m, 1.0, 'sfo')
440  for ob in objectlist:
441  rs.add_restraint(ob.get_restraint())
443  self.cg.set_scoring_function(sf)
444 
445  def get_output(self):
446  output = {}
447  acceptances = []
448  output["ConjugatedGradients_Nframe"] = str(self.nframe)
449  return output
450 
451 
452 class ReplicaExchange(object):
453  """Sample using replica exchange"""
454 
455  def __init__(
456  self,
457  model,
458  tempmin,
459  tempmax,
460  samplerobjects,
461  test=True,
462  replica_exchange_object=None):
463  '''
464  samplerobjects can be a list of MonteCarlo or MolecularDynamics
465  '''
466 
467 
468  self.m = model
469  self.samplerobjects = samplerobjects
470  # min and max temperature
471  self.TEMPMIN_ = tempmin
472  self.TEMPMAX_ = tempmax
473 
474  if replica_exchange_object is None:
475  # initialize Replica Exchange class
476  try:
477  import IMP.mpi
478  print('ReplicaExchange: MPI was found. Using Parallel Replica Exchange')
479  self.rem = IMP.mpi.ReplicaExchange()
480  except ImportError:
481  print('ReplicaExchange: Could not find MPI. Using Serial Replica Exchange')
482  self.rem = _SerialReplicaExchange()
483 
484  else:
485  # get the replica exchange class instance from elsewhere
486  print('got existing rex object')
487  self.rem = replica_exchange_object
488 
489  # get number of replicas
490  nproc = self.rem.get_number_of_replicas()
491 
492  if nproc % 2 != 0 and test == False:
493  raise Exception("number of replicas has to be even. set test=True to run with odd number of replicas.")
494  # create array of temperatures, in geometric progression
495  temp = self.rem.create_temperatures(
496  self.TEMPMIN_,
497  self.TEMPMAX_,
498  nproc)
499  # get replica index
500  self.temperatures = temp
501 
502  myindex = self.rem.get_my_index()
503  # set initial value of the parameter (temperature) to exchange
504  self.rem.set_my_parameter("temp", [self.temperatures[myindex]])
505  for so in self.samplerobjects:
506  so.set_kt(self.temperatures[myindex])
507  self.nattempts = 0
508  self.nmintemp = 0
509  self.nmaxtemp = 0
510  self.nsuccess = 0
511 
512  def get_temperatures(self):
513  return self.temperatures
514 
515  def get_my_temp(self):
516  return self.rem.get_my_parameter("temp")[0]
517 
518  def get_my_index(self):
519  return self.rem.get_my_index()
520 
521  def swap_temp(self, nframe, score=None):
522  if score is None:
523  score = self.m.evaluate(False)
524  # get my replica index and temperature
525  myindex = self.rem.get_my_index()
526  mytemp = self.rem.get_my_parameter("temp")[0]
527 
528  if mytemp == self.TEMPMIN_:
529  self.nmintemp += 1
530 
531  if mytemp == self.TEMPMAX_:
532  self.nmaxtemp += 1
533 
534  # score divided by kbt
535  myscore = score / mytemp
536 
537  # get my friend index and temperature
538  findex = self.rem.get_friend_index(nframe)
539  ftemp = self.rem.get_friend_parameter("temp", findex)[0]
540  # score divided by kbt
541  fscore = score / ftemp
542 
543  # try exchange
544  flag = self.rem.do_exchange(myscore, fscore, findex)
545 
546  self.nattempts += 1
547  # if accepted, change temperature
548  if (flag):
549  for so in self.samplerobjects:
550  so.set_kt(ftemp)
551  self.nsuccess += 1
552 
553  def get_output(self):
554  output = {}
555  acceptances = []
556  if self.nattempts != 0:
557  output["ReplicaExchange_SwapSuccessRatio"] = str(
558  float(self.nsuccess) / self.nattempts)
559  output["ReplicaExchange_MinTempFrequency"] = str(
560  float(self.nmintemp) / self.nattempts)
561  output["ReplicaExchange_MaxTempFrequency"] = str(
562  float(self.nmaxtemp) / self.nattempts)
563  else:
564  output["ReplicaExchange_SwapSuccessRatio"] = str(0)
565  output["ReplicaExchange_MinTempFrequency"] = str(0)
566  output["ReplicaExchange_MaxTempFrequency"] = str(0)
567  output["ReplicaExchange_CurrentTemp"] = str(self.get_my_temp())
568  return output
569 
570 
571 class PyMCMover(object):
572  # only works if the sampled particles are rigid bodies
573 
574  def __init__(self, representation, mcchild, n_mc_steps):
575 
576  # mcchild must be pmi.samplers.MonteCarlo
577  # representation must be pmi.representation
578 
579  self.rbs = representation.get_rigid_bodies()
580 
581  self.mc = mcchild
582  self.n_mc_steps = n_mc_steps
583 
584  def store_move(self):
585  # get all xyz coordinates of all rigid bodies of all copies
586  self.oldcoords = []
587  for copy in self.rbs:
588  crd = []
589  for rb in copy:
590  crd.append(rb.get_reference_frame())
591  self.oldcoords.append(crd)
592 
593  def propose_move(self, prob):
594  self.mc.run(self.n_mc_steps)
595 
596  def reset_move(self):
597  # reset each copy to crd
598  for copy, crd in zip(self.rbs, self.oldcoords):
599  for rb, ref in zip(copy, crd):
600  rb.set_reference_frame(ref)
601 
602  def get_number_of_steps(self):
603  return self.n_mc_steps
604 
605  def set_number_of_steps(self, nsteps):
606  self.n_mc_steps = nsteps
607 
608 
609 class PyMC(object):
610 
611  def __init__(self, model):
612  from math import exp
613  import random
614 
615  self.m = model
616  self.restraints = None
617  self.first_call = True
618  self.nframe = -1
619 
620  def add_mover(self, mv):
621  self.mv = mv
622 
623  def set_kt(self, kT):
624  self.kT = kT
625 
626  def set_return_best(self, thing):
627  pass
628 
629  def set_move_probability(self, thing):
630  pass
631 
632  def get_energy(self):
633  if self.restraints:
634  pot = sum([r.evaluate(False) for r in self.restraints])
635  else:
636  pot = self.m.evaluate(False)
637  return pot
638 
639  def metropolis(self, old, new):
640  deltaE = new - old
641  print(": old %f new %f deltaE %f new_epot: %f" % (old, new, deltaE,
642  self.m.evaluate(
643  False)), end=' ')
644  kT = self.kT
645  if deltaE < 0:
646  return True
647  else:
648  return exp(-deltaE / kT) > random.uniform(0, 1)
649 
650  def optimize(self, nsteps):
651  self.naccept = 0
652  self.nframe += 1
653  print("=== new MC call")
654  # store initial coordinates
655  if self.first_call:
656  self.mv.store_move()
657  self.first_call = False
658  for i in range(nsteps):
659  print("MC step %d " % i, end=' ')
660  # get total energy
661  old = self.get_energy()
662  # make a MD move
663  self.mv.propose_move(1)
664  # get new total energy
665  new = self.get_energy()
666  if self.metropolis(old, new):
667  # move was accepted: store new conformation
668  self.mv.store_move()
669  self.naccept += 1
670  print("accepted ")
671  else:
672  # move rejected: restore old conformation
673  self.mv.reset_move()
674  print(" ")
675 
676  def get_number_of_forward_steps(self):
677  return self.naccept
678 
679  def set_restraints(self, restraints):
680  self.restraints = restraints
681 
682  def set_scoring_function(self, objects):
683  # objects should be pmi.restraints
684  rs = IMP.RestraintSet(self.m, 1.0, 'sfo')
685  for ob in objects:
686  rs.add_restraint(ob.get_restraint())
687  self.set_restraints([rs])
688 
689  def get_output(self):
690  output = {}
691  acceptances = []
692  output["PyMC_Temperature"] = str(self.kT)
693  output["PyMC_Nframe"] = str(self.nframe)
694  return output
695 
696 
697 # 3
def __init__
samplerobjects can be a list of MonteCarlo or MolecularDynamics
Definition: samplers.py:455
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:349
def get_nuisance_movers_parameters
Return a dictionary with the mover parameters for nuisance parameters.
Definition: samplers.py:168
Miscellaneous utilities.
Definition: tools.py:1
Modify the transformation of a rigid body.
Simple conjugate gradients optimizer.
Sample using conjugate gradients.
Definition: samplers.py:415
Modify a set of continuous variables by perturbing them within a ball.
Object used to hold a set of restraints.
Definition: RestraintSet.h:36
Simple molecular dynamics optimizer.
def deprecated_method
Python decorator to mark a method as deprecated.
Definition: __init__.py:9537
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:46
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:36
The general base class for IMP exceptions.
Definition: exception.h:49
VectorD< 3 > Vector3D
Definition: VectorD.h:395
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:490
Sample using replica exchange.
Definition: samplers.py:452
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...