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