IMP logo
IMP Reference Guide  2.10.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, 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("2.10", "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.pmi.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  def get_nuisance_movers(self, nuisances, maxstep):
259  mvs = []
260  for nuisance in nuisances:
261  print(nuisance, maxstep)
262  mvs.append(
263  IMP.core.NormalMover([nuisance],
264  IMP.FloatKeys([IMP.FloatKey("nuisance")]),
265  maxstep))
266  return mvs
267 
268  def get_rigid_body_movers(self, rbs, maxtrans, maxrot):
269  mvs = []
270  for rb in rbs:
271  mvs.append(IMP.core.RigidBodyMover(rb.get_model(), rb,
272  maxtrans, maxrot))
273  return mvs
274 
275  def get_super_rigid_body_movers(self, rbs, maxtrans, maxrot):
276  mvs = []
277  for rb in rbs:
278  if len(rb) == 2:
279  # normal Super Rigid Body
280  srbm = IMP.pmi.TransformMover(self.model, maxtrans, maxrot)
281  elif len(rb) == 3:
282  if type(rb[2]) == tuple and type(rb[2][0]) == float \
283  and type(rb[2][1]) == float and type(rb[2][2]) == float \
284  and len(rb[2])== 3:
285  # super rigid body with 2D rotation, rb[2] is the axis
286  srbm = IMP.pmi.TransformMover(
287  self.model,
288  IMP.algebra.Vector3D(rb[2]),
289  maxtrans,
290  maxrot)
291  #elif type(rb[2]) == tuple and type(rb[2][0]) == IMP.Particle \
292  # and type(rb[2][1]) == IMP.Particle and len(rb[2])== 2:
293  # # super rigid body with bond rotation
294 
295  # srbm = IMP.pmi.TransformMover(
296  # self.model,
297  # rb[2][0],rb[2][1],
298  # 0, #no translation
299  # maxrot)
300  else:
301  print("Setting up a super rigid body with wrong parameters")
302  raise
303 
304  for xyz in rb[0]:
305  srbm.add_xyz_particle(xyz)
306  for rb in rb[1]:
307  srbm.add_rigid_body_particle(rb)
308  mvs.append(srbm)
309  return mvs
310 
311  def get_floppy_body_movers(self, fbs, maxtrans):
312  mvs = []
313  for fb in fbs:
314  # check is that is a rigid body member:
316  # if so force the particles to move anyway
317  floatkeys = [IMP.FloatKey(4), IMP.FloatKey(5), IMP.FloatKey(6)]
318  for fk in floatkeys:
319  fb.set_is_optimized(fk, True)
320  mvs.append(
321  IMP.core.BallMover(fb.get_model(), fb,
322  IMP.FloatKeys(floatkeys),
323  maxtrans))
324  else:
325  # otherwise use the normal ball mover
326  mvs.append(IMP.core.BallMover(fb.get_model(), fb, maxtrans))
327  return mvs
328 
329  def get_X_movers(self, fbs, maxtrans):
330  mvs = []
331  Xfloatkey = IMP.core.XYZ.get_xyz_keys()[0]
332  for fb in fbs:
333  # check is that is a rigid body member:
335  raise ValueError("particle is part of a rigid body")
336  else:
337  # otherwise use the normal ball mover
338  mvs.append(IMP.core.NormalMover([fb], [Xfloatkey], maxtrans))
339  return mvs
340 
341  def get_weight_movers(self, weights, maxstep):
342  mvs = []
343  for weight in weights:
344  if(weight.get_number_of_states() > 1):
345  mvs.append(IMP.isd.WeightMover(weight, maxstep))
346  return mvs
347 
348  def get_surface_movers(self, surfaces, maxtrans, maxrot, refprob):
349  mvs = []
350  for surface in surfaces:
351  mvs.append(IMP.core.SurfaceMover(surface, maxtrans, maxrot,
352  refprob))
353  return mvs
354 
355  def temp_simulated_annealing(self):
356  if self.nframe % (self.timemin + self.timemax) < self.timemin:
357  value = 0.0
358  else:
359  value = 1.0
360  temp = self.tempmin + (self.tempmax - self.tempmin) * value
361  return temp
362 
363  def set_label(self, label):
364  self.label = label
365 
366  def get_frame_number(self):
367  return self.nframe
368 
369  def get_output(self):
370  output = {}
371  acceptances = []
372  for i, mv in enumerate(self.smv.get_movers()):
373  mvname = mv.get_name()
374  mvacc = mv.get_number_of_accepted()
375  mvprp = mv.get_number_of_proposed()
376  try:
377  mvacr = float(mvacc) / float(mvprp)
378  except:
379  mvacr = 0.0
380  output["MonteCarlo_Acceptance_" +
381  mvname + "_" + str(i)] = str(mvacr)
382  if "Nuisances" in mvname:
383  output["MonteCarlo_StepSize_" + mvname + "_" +
384  str(i)] = str(IMP.core.NormalMover.get_from(mv).get_sigma())
385  if "Weights" in mvname:
386  output["MonteCarlo_StepSize_" + mvname + "_" +
387  str(i)] = str(IMP.isd.WeightMover.get_from(mv).get_radius())
388  output["MonteCarlo_Temperature"] = str(self.mc.get_kt())
389  output["MonteCarlo_Nframe"] = str(self.nframe)
390  return output
391 
392 
393 class MolecularDynamics(object):
394  """Sample using molecular dynamics"""
395 
396  def __init__(self,model,objects,kt,gamma=0.01,maximum_time_step=1.0,sf=None):
397  """Setup MD
398  @param model The IMP Model
399  @param objects What to sample. Use flat list of particles or (deprecated) 'MD Sample Objects' from PMI1
400  @param kt Temperature
401  @param gamma Viscosity parameter
402  @param maximum_time_step MD max time step
403  """
404  self.model=model
405 
406  # check if using PMI1 objects dictionary, or just list of particles
407  try:
408  for obj in objects:
409  to_sample=obj.get_particles_to_sample()['Floppy_Bodies_SimplifiedModel'][0]
410  except:
411  to_sample = objects
412 
413  self.ltstate=IMP.atom.LangevinThermostatOptimizerState(self.model,to_sample,
414  kt/0.0019872041,
415  gamma)
416  self.md = IMP.atom.MolecularDynamics(self.model)
417  self.md.set_maximum_time_step(maximum_time_step)
418  if sf:
419  self.md.set_scoring_function(sf)
420  else:
421  self.md.set_scoring_function(get_restraint_set(self.model))
422  self.md.add_optimizer_state(self.ltstate)
423  self.simulated_annealing = False
424  self.nframe = -1
425 
426  @property
427  @IMP.deprecated_method("2.10", "Model should be accessed with `.model`.")
428  def m(self):
429  return self.model
430 
431  def set_kt(self,kt):
432  temp=kt/0.0019872041
433  self.ltstate.set_temperature(temp)
434  self.md.assign_velocities(temp)
435 
436  def set_simulated_annealing(self, min_temp, max_temp, min_temp_time,
437  max_temp_time):
438  self.simulated_annealing = True
439  self.tempmin = min_temp
440  self.tempmax = max_temp
441  self.timemin = min_temp_time
442  self.timemax = max_temp_time
443 
444  def temp_simulated_annealing(self):
445  if self.nframe % (self.timemin + self.timemax) < self.timemin:
446  value = 0.0
447  else:
448  value = 1.0
449  temp = self.tempmin + (self.tempmax - self.tempmin) * value
450  return temp
451 
452  def set_gamma(self,gamma):
453  self.ltstate.set_gamma(gamma)
454 
455  def optimize(self,nsteps):
456  # apply simulated annealing protocol
457  self.nframe+=1
458  if self.simulated_annealing:
459  self.temp = self.temp_simulated_annealing()
460  self.set_kt(self.temp)
461  self.md.optimize(nsteps)
462 
463  def get_output(self):
464  output={}
465  output["MolecularDynamics_KineticEnergy"]=str(self.md.get_kinetic_energy())
466  return output
467 
468 class ConjugateGradients(object):
469  """Sample using conjugate gradients"""
470 
471  def __init__(self, model, objects):
472  self.model = model
473  self.nframe = -1
474  self.cg = IMP.core.ConjugateGradients(self.model)
475  self.cg.set_scoring_function(get_restraint_set(self.model))
476 
477  @property
478  @IMP.deprecated_method("2.10", "Model should be accessed with `.model`.")
479  def m(self):
480  return 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  acceptances = []
502  output["ConjugatedGradients_Nframe"] = str(self.nframe)
503  return output
504 
505 
506 
507 
508 
509 
510 class ReplicaExchange(object):
511  """Sample using replica exchange"""
512 
513  def __init__(
514  self,
515  model,
516  tempmin,
517  tempmax,
518  samplerobjects,
519  test=True,
520  replica_exchange_object=None):
521  '''
522  samplerobjects can be a list of MonteCarlo or MolecularDynamics
523  '''
524 
525 
526  self.model = model
527  self.samplerobjects = samplerobjects
528  # min and max temperature
529  self.TEMPMIN_ = tempmin
530  self.TEMPMAX_ = tempmax
531 
532  if replica_exchange_object is None:
533  # initialize Replica Exchange class
534  try:
535  import IMP.mpi
536  print('ReplicaExchange: MPI was found. Using Parallel Replica Exchange')
537  self.rem = IMP.mpi.ReplicaExchange()
538  except ImportError:
539  print('ReplicaExchange: Could not find MPI. Using Serial Replica Exchange')
540  self.rem = _SerialReplicaExchange()
541 
542  else:
543  # get the replica exchange class instance from elsewhere
544  print('got existing rex object')
545  self.rem = replica_exchange_object
546 
547  # get number of replicas
548  nproc = self.rem.get_number_of_replicas()
549 
550  if nproc % 2 != 0 and test == False:
551  raise Exception("number of replicas has to be even. set test=True to run with odd number of replicas.")
552  # create array of temperatures, in geometric progression
553  temp = self.rem.create_temperatures(
554  self.TEMPMIN_,
555  self.TEMPMAX_,
556  nproc)
557  # get replica index
558  self.temperatures = temp
559 
560  myindex = self.rem.get_my_index()
561  # set initial value of the parameter (temperature) to exchange
562  self.rem.set_my_parameter("temp", [self.temperatures[myindex]])
563  for so in self.samplerobjects:
564  so.set_kt(self.temperatures[myindex])
565  self.nattempts = 0
566  self.nmintemp = 0
567  self.nmaxtemp = 0
568  self.nsuccess = 0
569 
570  @property
571  @IMP.deprecated_method("2.10", "Model should be accessed with `.model`.")
572  def m(self):
573  return self.model
574 
575  def get_temperatures(self):
576  return self.temperatures
577 
578  def get_my_temp(self):
579  return self.rem.get_my_parameter("temp")[0]
580 
581  def get_my_index(self):
582  return self.rem.get_my_index()
583 
584  def swap_temp(self, nframe, score=None):
585  if score is None:
586  score = self.model.evaluate(False)
587  # get my replica index and temperature
588  myindex = self.rem.get_my_index()
589  mytemp = self.rem.get_my_parameter("temp")[0]
590 
591  if mytemp == self.TEMPMIN_:
592  self.nmintemp += 1
593 
594  if mytemp == self.TEMPMAX_:
595  self.nmaxtemp += 1
596 
597  # score divided by kbt
598  myscore = score / mytemp
599 
600  # get my friend index and temperature
601  findex = self.rem.get_friend_index(nframe)
602  ftemp = self.rem.get_friend_parameter("temp", findex)[0]
603  # score divided by kbt
604  fscore = score / ftemp
605 
606  # try exchange
607  flag = self.rem.do_exchange(myscore, fscore, findex)
608 
609  self.nattempts += 1
610  # if accepted, change temperature
611  if (flag):
612  for so in self.samplerobjects:
613  so.set_kt(ftemp)
614  self.nsuccess += 1
615 
616  def get_output(self):
617  output = {}
618  acceptances = []
619  if self.nattempts != 0:
620  output["ReplicaExchange_SwapSuccessRatio"] = str(
621  float(self.nsuccess) / self.nattempts)
622  output["ReplicaExchange_MinTempFrequency"] = str(
623  float(self.nmintemp) / self.nattempts)
624  output["ReplicaExchange_MaxTempFrequency"] = str(
625  float(self.nmaxtemp) / self.nattempts)
626  else:
627  output["ReplicaExchange_SwapSuccessRatio"] = str(0)
628  output["ReplicaExchange_MinTempFrequency"] = str(0)
629  output["ReplicaExchange_MaxTempFrequency"] = str(0)
630  output["ReplicaExchange_CurrentTemp"] = str(self.get_my_temp())
631  return output
632 
633 
634 
635 class MPI_values(object):
636  def __init__(self,replica_exchange_object=None):
637  """Query values (ie score, and others)
638  from a set of parallel jobs"""
639 
640  if replica_exchange_object is None:
641  # initialize Replica Exchange class
642  try:
643  import IMP.mpi
644  print('MPI_values: MPI was found. Using Parallel Replica Exchange')
645  self.rem = IMP.mpi.ReplicaExchange()
646  except ImportError:
647  print('MPI_values: Could not find MPI. Using Serial Replica Exchange')
648  self.rem = _SerialReplicaExchange()
649 
650  else:
651  # get the replica exchange class instance from elsewhere
652  print('got existing rex object')
653  self.rem = replica_exchange_object
654 
655  def set_value(self,name,value):
656  self.rem.set_my_parameter(name,[value])
657 
658  def get_values(self,name):
659  values=[]
660  for i in range(self.rem.get_number_of_replicas()):
661  v=self.rem.get_friend_parameter(name, i)[0]
662  values.append(v)
663  return values
664 
665  def get_percentile(self,name):
666  value=self.rem.get_my_parameter(name)[0]
667  values=sorted(self.get_values(name))
668  ind=values.index(value)
669  percentile=float(ind)/len(values)
670  return percentile
671 
672 
673 
674 class PyMCMover(object):
675  # only works if the sampled particles are rigid bodies
676 
677  def __init__(self, representation, mcchild, n_mc_steps):
678 
679  # mcchild must be pmi.samplers.MonteCarlo
680  # representation must be pmi.representation
681 
682  self.rbs = representation.get_rigid_bodies()
683 
684  self.mc = mcchild
685  self.n_mc_steps = n_mc_steps
686 
687  def store_move(self):
688  # get all xyz coordinates of all rigid bodies of all copies
689  self.oldcoords = []
690  for copy in self.rbs:
691  crd = []
692  for rb in copy:
693  crd.append(rb.get_reference_frame())
694  self.oldcoords.append(crd)
695 
696  def propose_move(self, prob):
697  self.mc.run(self.n_mc_steps)
698 
699  def reset_move(self):
700  # reset each copy to crd
701  for copy, crd in zip(self.rbs, self.oldcoords):
702  for rb, ref in zip(copy, crd):
703  rb.set_reference_frame(ref)
704 
705  def get_number_of_steps(self):
706  return self.n_mc_steps
707 
708  def set_number_of_steps(self, nsteps):
709  self.n_mc_steps = nsteps
710 
711 
712 class PyMC(object):
713 
714  def __init__(self, model):
715  from math import exp
716  import random
717 
718  self.model = model
719  self.restraints = None
720  self.first_call = True
721  self.nframe = -1
722 
723  @property
724  @IMP.deprecated_method("2.10", "Model should be accessed with `.model`.")
725  def m(self):
726  return self.model
727 
728  def add_mover(self, mv):
729  self.mv = mv
730 
731  def set_kt(self, kT):
732  self.kT = kT
733 
734  def set_return_best(self, thing):
735  pass
736 
737  def set_move_probability(self, thing):
738  pass
739 
740  def get_energy(self):
741  if self.restraints:
742  pot = sum([r.evaluate(False) for r in self.restraints])
743  else:
744  pot = self.model.evaluate(False)
745  return pot
746 
747  def metropolis(self, old, new):
748  deltaE = new - old
749  print(": old %f new %f deltaE %f new_epot: %f" % (old, new, deltaE,
750  self.model.evaluate(
751  False)), end=' ')
752  kT = self.kT
753  if deltaE < 0:
754  return True
755  else:
756  return exp(-deltaE / kT) > random.uniform(0, 1)
757 
758  def optimize(self, nsteps):
759  self.naccept = 0
760  self.nframe += 1
761  print("=== new MC call")
762  # store initial coordinates
763  if self.first_call:
764  self.mv.store_move()
765  self.first_call = False
766  for i in range(nsteps):
767  print("MC step %d " % i, end=' ')
768  # get total energy
769  old = self.get_energy()
770  # make a MD move
771  self.mv.propose_move(1)
772  # get new total energy
773  new = self.get_energy()
774  if self.metropolis(old, new):
775  # move was accepted: store new conformation
776  self.mv.store_move()
777  self.naccept += 1
778  print("accepted ")
779  else:
780  # move rejected: restore old conformation
781  self.mv.reset_move()
782  print(" ")
783 
784  def get_number_of_forward_steps(self):
785  return self.naccept
786 
787  def set_restraints(self, restraints):
788  self.restraints = restraints
789 
790  def set_scoring_function(self, objects):
791  # objects should be pmi.restraints
792  rs = IMP.RestraintSet(self.model, 1.0, 'sfo')
793  for ob in objects:
794  rs.add_restraint(ob.get_restraint())
795  self.set_restraints([rs])
796 
797  def get_output(self):
798  output = {}
799  acceptances = []
800  output["PyMC_Temperature"] = str(self.kT)
801  output["PyMC_Nframe"] = str(self.nframe)
802  return output
803 
804 
805 # 3
def __init__
samplerobjects can be a list of MonteCarlo or MolecularDynamics
Definition: samplers.py:513
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:393
def get_nuisance_movers_parameters
Return a dictionary with the mover parameters for nuisance parameters.
Definition: samplers.py:185
Miscellaneous utilities.
Definition: tools.py:1
Modify the transformation of a rigid body.
Simple conjugate gradients optimizer.
Sample using conjugate gradients.
Definition: samplers.py:468
Move continuous particle variables by perturbing them within a ball.
Modify a surface orientation.
Definition: SurfaceMover.h:34
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:9703
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:711
Sample using replica exchange.
Definition: samplers.py:510
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...