IMP logo
IMP Reference Guide  2.14.0
The Integrative Modeling Platform
samplers.py
1 """@namespace IMP.pmi.samplers
2  Sampling of the system.
3 """
4 
5 from __future__ import print_function
6 import IMP
7 import IMP.core
8 from IMP.pmi.tools import get_restraint_set
9 
10 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  def set_kt(self, temp):
152  self.temp = temp
153  self.mc.set_kt(temp)
154 
155  def get_mc(self):
156  return self.mc
157 
158  def set_scoring_function(self, objectlist):
159  rs = IMP.RestraintSet(self.model, 1.0, 'sfo')
160  for ob in objectlist:
161  rs.add_restraint(ob.get_restraint())
163  self.mc.set_scoring_function(sf)
164 
165  def set_simulated_annealing(
166  self,
167  min_temp,
168  max_temp,
169  min_temp_time,
170  max_temp_time):
171  self.simulated_annealing = True
172  self.tempmin = min_temp
173  self.tempmax = max_temp
174  self.timemin = min_temp_time
175  self.timemax = max_temp_time
176 
177  def set_self_adaptive(self, isselfadaptive=True):
178  self.selfadaptive = isselfadaptive
179 
181  '''
182  Return a dictionary with the mover parameters for nuisance parameters
183  '''
184  output = {}
185  for i in range(self.get_number_of_movers()):
186  mv = self.smv.get_mover(i)
187  name = mv.get_name()
188  if "Nuisances" in name:
189  stepsize = IMP.core.NormalMover.get_from(mv).get_sigma()
190  output[name] = stepsize
191  return output
192 
193  def get_number_of_movers(self):
194  return len(self.smv.get_movers())
195 
196  def get_particle_types(self):
197  return self.losp
198 
199  def optimize(self, nstep):
200  self.nframe += 1
201  self.mc.optimize(nstep * self.get_number_of_movers())
202 
203  # apply simulated annealing protocol
204  if self.simulated_annealing:
205  self.temp = self.temp_simulated_annealing()
206  self.mc.set_kt(self.temp)
207 
208  # apply self adaptive protocol
209  if self.selfadaptive:
210  for i, mv in enumerate(self.mvs):
211 
212  mvacc = mv.get_number_of_accepted()
213  mvprp = mv.get_number_of_proposed()
214  if mv not in self.movers_data:
215  accept = float(mvacc) / float(mvprp)
216  self.movers_data[mv]=(mvacc,mvprp)
217  else:
218  oldmvacc,oldmvprp=self.movers_data[mv]
219  accept = float(mvacc-oldmvacc) / float(mvprp-oldmvprp)
220  self.movers_data[mv]=(mvacc,mvprp)
221  if accept < 0.05: accept = 0.05
222  if accept > 1.0: accept = 1.0
223 
224  if isinstance(mv, IMP.core.NormalMover):
225  stepsize = mv.get_sigma()
226  if 0.4 > accept or accept > 0.6:
227  mv.set_sigma(stepsize * 2 * accept)
228 
229  if isinstance(mv, MP.isd.WeightMover):
230  stepsize = mv.get_radius()
231  if 0.4 > accept or accept > 0.6:
232  mv.set_radius(stepsize * 2 * accept)
233 
234  if isinstance(mv, IMP.core.RigidBodyMover):
235  mr=mv.get_maximum_rotation()
236  mt=mv.get_maximum_translation()
237  if 0.4 > accept or accept > 0.6:
238  mv.set_maximum_rotation(mr * 2 * accept)
239  mv.set_maximum_translation(mt * 2 * accept)
240 
241  if isinstance(mv, IMP.pmi.TransformMover):
242  mr=mv.get_maximum_rotation()
243  mt=mv.get_maximum_translation()
244  if 0.4 > accept or accept > 0.6:
245  mv.set_maximum_rotation(mr * 2 * accept)
246  mv.set_maximum_translation(mt * 2 * accept)
247 
248  if isinstance(mv, IMP.core.BallMover):
249  mr=mv.get_radius()
250  if 0.4 > accept or accept > 0.6:
251  mv.set_radius(mr * 2 * accept)
252 
253  def get_nuisance_movers(self, nuisances, maxstep):
254  mvs = []
255  for nuisance in nuisances:
256  print(nuisance, maxstep)
257  mvs.append(
258  IMP.core.NormalMover([nuisance],
259  IMP.FloatKeys([IMP.FloatKey("nuisance")]),
260  maxstep))
261  return mvs
262 
263  def get_rigid_body_movers(self, rbs, maxtrans, maxrot):
264  mvs = []
265  for rb in rbs:
266  mvs.append(IMP.core.RigidBodyMover(rb.get_model(), rb,
267  maxtrans, maxrot))
268  return mvs
269 
270  def get_super_rigid_body_movers(self, rbs, maxtrans, maxrot):
271  mvs = []
272  for rb in rbs:
273  if len(rb) == 2:
274  # normal Super Rigid Body
275  srbm = IMP.pmi.TransformMover(self.model, maxtrans, maxrot)
276  elif len(rb) == 3:
277  if isinstance(rb[2], tuple) and len(rb[2]) == 3 \
278  and isinstance(rb[2][0], float) \
279  and isinstance(rb[2][1], float) \
280  and isinstance(rb[2][2], float):
281  # super rigid body with 2D rotation, rb[2] is the axis
282  srbm = IMP.pmi.TransformMover(
283  self.model,
284  IMP.algebra.Vector3D(rb[2]),
285  maxtrans,
286  maxrot)
287  #elif type(rb[2]) == tuple and type(rb[2][0]) == IMP.Particle \
288  # and type(rb[2][1]) == IMP.Particle and len(rb[2])== 2:
289  # # super rigid body with bond rotation
290 
291  # srbm = IMP.pmi.TransformMover(
292  # self.model,
293  # rb[2][0],rb[2][1],
294  # 0, #no translation
295  # maxrot)
296  else:
297  print("Setting up a super rigid body with wrong parameters")
298  raise
299 
300  for xyz in rb[0]:
301  srbm.add_xyz_particle(xyz)
302  for rb in rb[1]:
303  srbm.add_rigid_body_particle(rb)
304  mvs.append(srbm)
305  return mvs
306 
307  def get_floppy_body_movers(self, fbs, maxtrans):
308  mvs = []
309  for fb in fbs:
310  # check is that is a rigid body member:
312  # if so force the particles to move anyway
313  floatkeys = [IMP.FloatKey(4), IMP.FloatKey(5), IMP.FloatKey(6)]
314  for fk in floatkeys:
315  fb.set_is_optimized(fk, True)
316  mvs.append(
317  IMP.core.BallMover(fb.get_model(), fb,
318  IMP.FloatKeys(floatkeys),
319  maxtrans))
320  else:
321  # otherwise use the normal ball mover
322  mvs.append(IMP.core.BallMover(fb.get_model(), fb, maxtrans))
323  return mvs
324 
325  def get_X_movers(self, fbs, maxtrans):
326  mvs = []
327  Xfloatkey = IMP.core.XYZ.get_xyz_keys()[0]
328  for fb in fbs:
329  # check is that is a rigid body member:
331  raise ValueError("particle is part of a rigid body")
332  else:
333  # otherwise use the normal ball mover
334  mvs.append(IMP.core.NormalMover([fb], [Xfloatkey], maxtrans))
335  return mvs
336 
337  def get_weight_movers(self, weights, maxstep):
338  mvs = []
339  for weight in weights:
340  if(weight.get_number_of_states() > 1):
341  mvs.append(IMP.isd.WeightMover(weight, maxstep))
342  return mvs
343 
344  def get_surface_movers(self, surfaces, maxtrans, maxrot, refprob):
345  mvs = []
346  for surface in surfaces:
347  mvs.append(IMP.core.SurfaceMover(surface, maxtrans, maxrot,
348  refprob))
349  return mvs
350 
351  def temp_simulated_annealing(self):
352  if self.nframe % (self.timemin + self.timemax) < self.timemin:
353  value = 0.0
354  else:
355  value = 1.0
356  temp = self.tempmin + (self.tempmax - self.tempmin) * value
357  return temp
358 
359  def set_label(self, label):
360  self.label = label
361 
362  def get_frame_number(self):
363  return self.nframe
364 
365  def get_output(self):
366  output = {}
367  acceptances = []
368  for i, mv in enumerate(self.smv.get_movers()):
369  mvname = mv.get_name()
370  mvacc = mv.get_number_of_accepted()
371  mvprp = mv.get_number_of_proposed()
372  try:
373  mvacr = float(mvacc) / float(mvprp)
374  except:
375  mvacr = 0.0
376  output["MonteCarlo_Acceptance_" +
377  mvname + "_" + str(i)] = str(mvacr)
378  if "Nuisances" in mvname:
379  output["MonteCarlo_StepSize_" + mvname + "_" +
380  str(i)] = str(IMP.core.NormalMover.get_from(mv).get_sigma())
381  if "Weights" in mvname:
382  output["MonteCarlo_StepSize_" + mvname + "_" +
383  str(i)] = str(IMP.isd.WeightMover.get_from(mv).get_radius())
384  output["MonteCarlo_Temperature"] = str(self.mc.get_kt())
385  output["MonteCarlo_Nframe"] = str(self.nframe)
386  return output
387 
388 
389 class MolecularDynamics(object):
390  """Sample using molecular dynamics"""
391 
392  def __init__(self,model,objects,kt,gamma=0.01,maximum_time_step=1.0,sf=None):
393  """Setup MD
394  @param model The IMP Model
395  @param objects What to sample. Use flat list of particles or (deprecated) 'MD Sample Objects' from PMI1
396  @param kt Temperature
397  @param gamma Viscosity parameter
398  @param maximum_time_step MD max time step
399  """
400  self.model=model
401 
402  # check if using PMI1 objects dictionary, or just list of particles
403  try:
404  for obj in objects:
405  to_sample=obj.get_particles_to_sample()['Floppy_Bodies_SimplifiedModel'][0]
406  except:
407  to_sample = objects
408 
409  self.ltstate=IMP.atom.LangevinThermostatOptimizerState(self.model,to_sample,
410  kt/0.0019872041,
411  gamma)
412  self.md = IMP.atom.MolecularDynamics(self.model)
413  self.md.set_maximum_time_step(maximum_time_step)
414  if sf:
415  self.md.set_scoring_function(sf)
416  else:
417  self.md.set_scoring_function(get_restraint_set(self.model))
418  self.md.add_optimizer_state(self.ltstate)
419  self.simulated_annealing = False
420  self.nframe = -1
421 
422  def set_kt(self,kt):
423  temp=kt/0.0019872041
424  self.ltstate.set_temperature(temp)
425  self.md.assign_velocities(temp)
426 
427  def set_simulated_annealing(self, min_temp, max_temp, min_temp_time,
428  max_temp_time):
429  self.simulated_annealing = True
430  self.tempmin = min_temp
431  self.tempmax = max_temp
432  self.timemin = min_temp_time
433  self.timemax = max_temp_time
434 
435  def temp_simulated_annealing(self):
436  if self.nframe % (self.timemin + self.timemax) < self.timemin:
437  value = 0.0
438  else:
439  value = 1.0
440  temp = self.tempmin + (self.tempmax - self.tempmin) * value
441  return temp
442 
443  def set_gamma(self,gamma):
444  self.ltstate.set_gamma(gamma)
445 
446  def optimize(self,nsteps):
447  # apply simulated annealing protocol
448  self.nframe+=1
449  if self.simulated_annealing:
450  self.temp = self.temp_simulated_annealing()
451  self.set_kt(self.temp)
452  self.md.optimize(nsteps)
453 
454  def get_output(self):
455  output={}
456  output["MolecularDynamics_KineticEnergy"]=str(self.md.get_kinetic_energy())
457  return output
458 
459 class ConjugateGradients(object):
460  """Sample using conjugate gradients"""
461 
462  def __init__(self, model, objects):
463  self.model = model
464  self.nframe = -1
465  self.cg = IMP.core.ConjugateGradients(self.model)
466  self.cg.set_scoring_function(get_restraint_set(self.model))
467 
468  def set_label(self, label):
469  self.label = label
470 
471  def get_frame_number(self):
472  return self.nframe
473 
474  def optimize(self, nstep):
475  self.nframe += 1
476  self.cg.optimize(nstep)
477 
478  def set_scoring_function(self, objectlist):
479  rs = IMP.RestraintSet(self.model, 1.0, 'sfo')
480  for ob in objectlist:
481  rs.add_restraint(ob.get_restraint())
483  self.cg.set_scoring_function(sf)
484 
485  def get_output(self):
486  output = {}
487  acceptances = []
488  output["ConjugatedGradients_Nframe"] = str(self.nframe)
489  return output
490 
491 
492 
493 
494 
495 
496 class ReplicaExchange(object):
497  """Sample using replica exchange"""
498 
499  def __init__(
500  self,
501  model,
502  tempmin,
503  tempmax,
504  samplerobjects,
505  test=True,
506  replica_exchange_object=None):
507  '''
508  samplerobjects can be a list of MonteCarlo or MolecularDynamics
509  '''
510 
511 
512  self.model = model
513  self.samplerobjects = samplerobjects
514  # min and max temperature
515  self.TEMPMIN_ = tempmin
516  self.TEMPMAX_ = tempmax
517 
518  if replica_exchange_object is None:
519  # initialize Replica Exchange class
520  try:
521  import IMP.mpi
522  print('ReplicaExchange: MPI was found. Using Parallel Replica Exchange')
523  self.rem = IMP.mpi.ReplicaExchange()
524  except ImportError:
525  print('ReplicaExchange: Could not find MPI. Using Serial Replica Exchange')
526  self.rem = _SerialReplicaExchange()
527 
528  else:
529  # get the replica exchange class instance from elsewhere
530  print('got existing rex object')
531  self.rem = replica_exchange_object
532 
533  # get number of replicas
534  nproc = self.rem.get_number_of_replicas()
535 
536  if nproc % 2 != 0 and test == False:
537  raise Exception("number of replicas has to be even. set test=True to run with odd number of replicas.")
538  # create array of temperatures, in geometric progression
539  temp = self.rem.create_temperatures(
540  self.TEMPMIN_,
541  self.TEMPMAX_,
542  nproc)
543  # get replica index
544  self.temperatures = temp
545 
546  myindex = self.rem.get_my_index()
547  # set initial value of the parameter (temperature) to exchange
548  self.rem.set_my_parameter("temp", [self.temperatures[myindex]])
549  for so in self.samplerobjects:
550  so.set_kt(self.temperatures[myindex])
551  self.nattempts = 0
552  self.nmintemp = 0
553  self.nmaxtemp = 0
554  self.nsuccess = 0
555 
556  def get_temperatures(self):
557  return self.temperatures
558 
559  def get_my_temp(self):
560  return self.rem.get_my_parameter("temp")[0]
561 
562  def get_my_index(self):
563  return self.rem.get_my_index()
564 
565  def swap_temp(self, nframe, score=None):
566  if score is None:
567  score = self.model.evaluate(False)
568  # get my replica index and temperature
569  myindex = self.rem.get_my_index()
570  mytemp = self.rem.get_my_parameter("temp")[0]
571 
572  if mytemp == self.TEMPMIN_:
573  self.nmintemp += 1
574 
575  if mytemp == self.TEMPMAX_:
576  self.nmaxtemp += 1
577 
578  # score divided by kbt
579  myscore = score / mytemp
580 
581  # get my friend index and temperature
582  findex = self.rem.get_friend_index(nframe)
583  ftemp = self.rem.get_friend_parameter("temp", findex)[0]
584  # score divided by kbt
585  fscore = score / ftemp
586 
587  # try exchange
588  flag = self.rem.do_exchange(myscore, fscore, findex)
589 
590  self.nattempts += 1
591  # if accepted, change temperature
592  if (flag):
593  for so in self.samplerobjects:
594  so.set_kt(ftemp)
595  self.nsuccess += 1
596 
597  def get_output(self):
598  output = {}
599  acceptances = []
600  if self.nattempts != 0:
601  output["ReplicaExchange_SwapSuccessRatio"] = str(
602  float(self.nsuccess) / self.nattempts)
603  output["ReplicaExchange_MinTempFrequency"] = str(
604  float(self.nmintemp) / self.nattempts)
605  output["ReplicaExchange_MaxTempFrequency"] = str(
606  float(self.nmaxtemp) / self.nattempts)
607  else:
608  output["ReplicaExchange_SwapSuccessRatio"] = str(0)
609  output["ReplicaExchange_MinTempFrequency"] = str(0)
610  output["ReplicaExchange_MaxTempFrequency"] = str(0)
611  output["ReplicaExchange_CurrentTemp"] = str(self.get_my_temp())
612  return output
613 
614 
615 
616 class MPI_values(object):
617  def __init__(self,replica_exchange_object=None):
618  """Query values (ie score, and others)
619  from a set of parallel jobs"""
620 
621  if replica_exchange_object is None:
622  # initialize Replica Exchange class
623  try:
624  import IMP.mpi
625  print('MPI_values: MPI was found. Using Parallel Replica Exchange')
626  self.rem = IMP.mpi.ReplicaExchange()
627  except ImportError:
628  print('MPI_values: Could not find MPI. Using Serial Replica Exchange')
629  self.rem = _SerialReplicaExchange()
630 
631  else:
632  # get the replica exchange class instance from elsewhere
633  print('got existing rex object')
634  self.rem = replica_exchange_object
635 
636  def set_value(self,name,value):
637  self.rem.set_my_parameter(name,[value])
638 
639  def get_values(self,name):
640  values=[]
641  for i in range(self.rem.get_number_of_replicas()):
642  v=self.rem.get_friend_parameter(name, i)[0]
643  values.append(v)
644  return values
645 
646  def get_percentile(self,name):
647  value=self.rem.get_my_parameter(name)[0]
648  values=sorted(self.get_values(name))
649  ind=values.index(value)
650  percentile=float(ind)/len(values)
651  return percentile
def __init__
samplerobjects can be a list of MonteCarlo or MolecularDynamics
Definition: samplers.py:499
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:389
def get_nuisance_movers_parameters
Return a dictionary with the mover parameters for nuisance parameters.
Definition: samplers.py:180
Miscellaneous utilities.
Definition: tools.py:1
Modify the transformation of a rigid body.
Simple conjugate gradients optimizer.
Sample using conjugate gradients.
Definition: samplers.py:459
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.
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:421
static const FloatKeys & get_xyz_keys()
Get a vector containing the keys for x,y,z.
Applies a list of movers one at a time.
Definition: SerialMover.h:23
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:752
Sample using replica exchange.
Definition: samplers.py:496
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...