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