IMP logo
IMP Reference Guide  2.10.0
The Integrative Modeling Platform
tools.py
1 #!/usr/bin/env python
2 
3 """@namespace IMP.pmi.tools
4  Miscellaneous utilities.
5 """
6 
7 from __future__ import print_function
8 import IMP
9 import IMP.algebra
10 import IMP.isd
11 import IMP.pmi
12 import IMP.pmi.topology
13 import collections
14 import itertools
15 from math import log,pi,sqrt,exp
16 import sys,os
17 import random
18 import ast
19 import time
20 import RMF
21 import IMP.rmf
22 from collections import defaultdict
23 try:
24  from collections import OrderedDict
25 except ImportError:
26  from IMP.pmi._compat_collections import OrderedDict
27 
28 def _get_system_for_hier(hier):
29  """Given a top-level hierarchy, return the System that created it,
30  or None"""
31  if hasattr(hier, '_pmi2_system'):
32  return hier._pmi2_system()
33 
34 def _all_protocol_outputs(representations, root_hier):
35  """Iterate over all (ProtocolOutput, State) pairs for the given
36  representations (PMI1) or root hier (PMI2)"""
37  if root_hier:
38  system = _get_system_for_hier(root_hier)
39  if system:
40  for state in system.states:
41  for p in state._protocol_output:
42  yield p
43  else:
44  for p in representations[0]._protocol_output:
45  yield p
46 
47 def _add_pmi_provenance(p):
48  """Tag the given particle as being created by the current version of PMI."""
50  IMP.core.add_software_provenance(p, name="IMP PMI module",
52  location="https://integrativemodeling.org")
54 
55 def _get_restraint_set_keys():
56  if not hasattr(_get_restraint_set_keys, 'pmi_rs_key'):
57  _get_restraint_set_keys.pmi_rs_key = IMP.ModelKey("PMI restraints")
58  _get_restraint_set_keys.rmf_rs_key = IMP.ModelKey("RMF restraints")
59  return (_get_restraint_set_keys.pmi_rs_key,
60  _get_restraint_set_keys.rmf_rs_key)
61 
62 def _add_restraint_sets(model, mk, mk_rmf):
63  rs = IMP.RestraintSet(model, "All PMI restraints")
64  rs_rmf = IMP.RestraintSet(model, "All PMI RMF restraints")
65  model.add_data(mk, rs)
66  model.add_data(mk_rmf, rs_rmf)
67  return rs, rs_rmf
68 
69 def add_restraint_to_model(model, restraint, add_to_rmf=False):
70  """Add a PMI restraint to the model.
71  Since Model.add_restraint() no longer exists (in modern IMP restraints
72  should be added to a ScoringFunction instead) store them instead in
73  a RestraintSet, and keep a reference to it in the Model.
74 
75  If `add_to_rmf` is True, also add the restraint to a separate list
76  of restraints that will be written out to RMF files (by default, most
77  PMI restraints are not)."""
78  mk, mk_rmf = _get_restraint_set_keys()
79  if model.get_has_data(mk):
80  rs = IMP.RestraintSet.get_from(model.get_data(mk))
81  rs_rmf = IMP.RestraintSet.get_from(model.get_data(mk_rmf))
82  else:
83  rs, rs_rmf = _add_restraint_sets(model, mk, mk_rmf)
84  rs.add_restraint(restraint)
85  if add_to_rmf:
86  rs_rmf.add_restraint(restraint)
87 
88 def get_restraint_set(model, rmf=False):
89  """Get a RestraintSet containing all PMI restraints added to the model.
90  If `rmf` is True, return only the subset of these restraints that
91  should be written out to RMF files."""
92  mk, mk_rmf = _get_restraint_set_keys()
93  if not model.get_has_data(mk):
94  print("WARNING: no restraints added to model yet")
95  _add_restraint_sets(model, mk, mk_rmf)
96  if rmf:
97  return IMP.RestraintSet.get_from(model.get_data(mk_rmf))
98  else:
99  return IMP.RestraintSet.get_from(model.get_data(mk))
100 
101 class Stopwatch(object):
102  """Collect timing information.
103  Add an instance of this class to outputobjects to get timing information
104  in a stat file."""
105 
106  def __init__(self, isdelta=True):
107  """Constructor.
108  @param isdelta if True (the default) then report the time since the
109  last use of this class; if False, report cumulative time."""
110  self.starttime = time.clock()
111  self.label = "None"
112  self.isdelta = isdelta
113 
114  def set_label(self, labelstr):
115  self.label = labelstr
116 
117  def get_output(self):
118  output = {}
119  if self.isdelta:
120  newtime = time.clock()
121  output["Stopwatch_" + self.label + "_delta_seconds"] \
122  = str(newtime - self.starttime)
123  self.starttime = newtime
124  else:
125  output["Stopwatch_" + self.label + "_elapsed_seconds"] \
126  = str(time.clock() - self.starttime)
127  return output
128 
129 
130 class SetupNuisance(object):
131 
132  def __init__(self, m, initialvalue, minvalue, maxvalue, isoptimized=True):
133 
134  nuisance = IMP.isd.Scale.setup_particle(IMP.Particle(m), initialvalue)
135  if minvalue:
136  nuisance.set_lower(minvalue)
137  if maxvalue:
138  nuisance.set_upper(maxvalue)
139 
140  # m.add_score_state(IMP.core.SingletonConstraint(IMP.isd.NuisanceRangeModifier(),None,nuisance))
141  nuisance.set_is_optimized(nuisance.get_nuisance_key(), isoptimized)
142  self.nuisance = nuisance
143 
144  def get_particle(self):
145  return self.nuisance
146 
147 
148 class SetupWeight(object):
149 
150  def __init__(self, m, isoptimized=True):
151  pw = IMP.Particle(m)
152  self.weight = IMP.isd.Weight.setup_particle(pw)
153  self.weight.set_weights_are_optimized(True)
154 
155  def get_particle(self):
156  return self.weight
157 
158 
159 class SetupSurface(object):
160 
161  def __init__(self, m, center, normal, isoptimized=True):
162  p = IMP.Particle(m)
163  self.surface = IMP.core.Surface.setup_particle(p, center, normal)
164  self.surface.set_coordinates_are_optimized(isoptimized)
165  self.surface.set_normal_is_optimized(isoptimized)
166 
167  def get_particle(self):
168  return self.surface
169 
170 
171 class ParticleToSampleList(object):
172 
173  def __init__(self, label="None"):
174 
175  self.dictionary_particle_type = {}
176  self.dictionary_particle_transformation = {}
177  self.dictionary_particle_name = {}
178  self.label = label
179 
180  def add_particle(
181  self,
182  particle,
183  particle_type,
184  particle_transformation,
185  name):
186  if not particle_type in ["Rigid_Bodies", "Floppy_Bodies", "Nuisances", "X_coord", "Weights", "Surfaces"]:
187  raise TypeError("not the right particle type")
188  else:
189  self.dictionary_particle_type[particle] = particle_type
190  if particle_type == "Rigid_Bodies":
191  if type(particle_transformation) == tuple and len(particle_transformation) == 2 and type(particle_transformation[0]) == float and type(particle_transformation[1]) == float:
192  self.dictionary_particle_transformation[
193  particle] = particle_transformation
194  self.dictionary_particle_name[particle] = name
195  else:
196  raise TypeError("ParticleToSampleList: not the right transformation format for Rigid_Bodies, should be a tuple of floats")
197  elif particle_type == "Surfaces":
198  if type(particle_transformation) == tuple and len(particle_transformation) == 3 and all(isinstance(x, float) for x in particle_transformation):
199  self.dictionary_particle_transformation[
200  particle] = particle_transformation
201  self.dictionary_particle_name[particle] = name
202  else:
203  raise TypeError("ParticleToSampleList: not the right transformation format for Surfaces, should be a tuple of floats")
204  else:
205  if type(particle_transformation) == float:
206  self.dictionary_particle_transformation[
207  particle] = particle_transformation
208  self.dictionary_particle_name[particle] = name
209  else:
210  raise TypeError("ParticleToSampleList: not the right transformation format, should be a float")
211 
212  def get_particles_to_sample(self):
213  ps = {}
214  for particle in self.dictionary_particle_type:
215  key = self.dictionary_particle_type[
216  particle] + "ParticleToSampleList_" + self.dictionary_particle_name[particle] + "_" + self.label
217  value = (
218  [particle],
219  self.dictionary_particle_transformation[particle])
220  ps[key] = value
221  return ps
222 
223 
225  "If you use this class please let the PMI developers know.")
226 class Variance(object):
227 
228  def __init__(self, model, tau, niter, prot, th_profile, write_data=False):
229 
230  self.model = model
231  self.write_data = write_data
232  self.tau = tau
233  self.niter = niter
234  #! select particles from the model
235  particles = IMP.atom.get_by_type(prot, IMP.atom.ATOM_TYPE)
236  self.particles = particles
237  # store reference coordinates and theoretical profile
238  self.refpos = [IMP.core.XYZ(p).get_coordinates() for p in particles]
239  self.model_profile = th_profile
240 
241  def perturb_particles(self, perturb=True):
242  for i, p in enumerate(self.particles):
243  newpos = array(self.refpos[i])
244  if perturb:
245  newpos += random.normal(0, self.tau, 3)
246  newpos = IMP.algebra.Vector3D(newpos)
247  IMP.core.XYZ(p).set_coordinates(newpos)
248 
249  def get_profile(self):
250  model_profile = self.model_profile
251  p = model_profile.calculate_profile(self.particles, IMP.saxs.CA_ATOMS)
252  return array([model_profile.get_intensity(i) for i in
253  range(model_profile.size())])
254 
255  def init_variances(self):
256  # create placeholders
257  N = self.model_profile.size()
258  a = self.profiles[0][:]
259  self.m = matrix(a).T # Nx1
260  self.V = self.m * self.m.T
261  self.normm = linalg.norm(self.m)
262  self.normV = linalg.norm(self.V)
263 
264  def update_variances(self):
265  a = matrix(self.profiles[-1]) # 1xN
266  n = float(len(self.profiles))
267  self.m = a.T / n + (n - 1) / n * self.m
268  self.V = a.T * a + self.V
269  self.oldnormm = self.normm
270  self.oldnormV = self.normV
271  self.normm = linalg.norm(self.m)
272  self.normV = linalg.norm(self.V)
273  self.diffm = (self.oldnormm - self.normm) / self.oldnormm
274  self.diffV = (self.oldnormV - self.normV) / self.oldnormV
275 
276  def get_direct_stats(self, a):
277  nq = len(a[0])
278  nprof = len(a)
279  m = [0] * nq
280  for prof in a:
281  for q, I in enumerate(prof):
282  m[q] += I
283  m = array(m) / nprof
284  V = matrix(a)
285  V = V.T * V
286  Sigma = (matrix(a - m))
287  Sigma = Sigma.T * Sigma / (nprof - 1)
288  mi = matrix(diag(1. / m))
289  Sigmarel = mi.T * Sigma * mi
290  return m, V, Sigma, Sigmarel
291 
292  def store_data(self):
293  if not os.path.isdir('data'):
294  os.mkdir('data')
295  profiles = matrix(self.profiles)
296  self.directm, self.directV, self.Sigma, self.Sigmarel = \
297  self.get_direct_stats(array(profiles))
298  directV = self.directV
299  # print "V comparison",(linalg.norm(directV-self.V)/self.normV)
300  save('data/profiles', profiles)
301  # absolute profile differences
302  fl = open('data/profiles.dat', 'w')
303  for i, l in enumerate(array(profiles).T):
304  self.model_profile.get_q(i)
305  fl.write('%s ' % i)
306  for k in l:
307  fl.write('%s ' % (k - self.directm[i]))
308  fl.write('\n')
309  # relative profile differences
310  fl = open('data/profiles_rel.dat', 'w')
311  for i, l in enumerate(array(profiles).T):
312  self.model_profile.get_q(i)
313  fl.write('%s ' % i)
314  for k in l:
315  fl.write('%s ' % ((k - self.directm[i]) / self.directm[i]))
316  fl.write('\n')
317  save('data/m', self.directm)
318  save('data/V', self.directV)
319  Sigma = self.Sigma
320  save('data/Sigma', Sigma)
321  # Sigma matrix
322  fl = open('data/Sigma.dat', 'w')
323  model_profile = self.model_profile
324  for i in range(model_profile.size()):
325  qi = model_profile.get_q(i)
326  for j in range(model_profile.size()):
327  qj = model_profile.get_q(j)
328  vij = self.Sigma[i, j]
329  fl.write('%s %s %s\n' % (qi, qj, vij))
330  fl.write('\n')
331  # Sigma eigenvalues
332  fl = open('data/eigenvals', 'w')
333  for i in linalg.eigvalsh(Sigma):
334  fl.write('%s\n' % i)
335  Sigmarel = self.Sigmarel
336  save('data/Sigmarel', Sigmarel)
337  # Sigmarel matrix
338  fl = open('data/Sigmarel.dat', 'w')
339  model_profile = self.model_profile
340  for i in range(model_profile.size()):
341  qi = model_profile.get_q(i)
342  for j in range(model_profile.size()):
343  qj = model_profile.get_q(j)
344  vij = self.Sigmarel[i, j]
345  fl.write('%s %s %s\n' % (qi, qj, vij))
346  fl.write('\n')
347  # Sigma eigenvalues
348  fl = open('data/eigenvals_rel', 'w')
349  for i in linalg.eigvalsh(Sigmarel):
350  fl.write('%s\n' % i)
351  # mean profile
352  fl = open('data/mean.dat', 'w')
353  for i in range(len(self.directm)):
354  qi = self.model_profile.get_q(i)
355  fl.write('%s ' % qi)
356  fl.write('%s ' % self.directm[i])
357  fl.write('%s ' % sqrt(self.Sigma[i, i]))
358  fl.write('\n')
359 
360  def try_chol(self, jitter):
361  Sigma = self.Sigma
362  try:
363  linalg.cholesky(Sigma + matrix(eye(len(Sigma))) * jitter)
364  except linalg.LinAlgError:
365  print("Decomposition failed with jitter =", jitter)
366  return
367  print("Successful decomposition with jitter =", jitter)
368 
369  def run(self):
370  self.profiles = [self.get_profile()]
371  # self.init_variances()
372  for n in range(self.niter):
373  self.perturb_particles()
374  self.profiles.append(self.get_profile())
375  # self.update_variances()
376  #profiles = matrix(self.profiles)
377  # print n,self.diffm,self.diffV
378  # print
379  #
380  if self.write_data:
381  self.store_data()
382  # self.try_chol(0.)
383  # for i in logspace(-7,0,num=8):
384  # self.try_chol(i)
385 
386  def get_cov(self, relative=True):
387  if not relative:
388  return self.Sigma
389  else:
390  return self.Sigmarel
391 
392  #-------------------------------
393 
394 def get_random_cross_link_dataset(representation,
395  resolution=1.0,
396  number_of_cross_links=10,
397  ambiguity_probability=0.1,
398  confidence_score_range=[0,100],
399  avoid_same_particles=False):
400  '''Return a random cross-link dataset as a string.
401  Every line is a residue pair, together with UniqueIdentifier
402  and XL score.'''
403 
404  residue_pairs=get_random_residue_pairs(representation, resolution, number_of_cross_links, avoid_same_particles=avoid_same_particles)
405 
406  unique_identifier=0
407  cmin=float(min(confidence_score_range))
408  cmax=float(max(confidence_score_range))
409 
410  dataset="#\n"
411 
412  for (name1, r1, name2, r2) in residue_pairs:
413  if random.random() > ambiguity_probability:
414  unique_identifier+=1
415  score=random.random()*(cmax-cmin)+cmin
416  dataset+=str(name1)+" "+str(name2)+" "+str(r1)+" "+str(r2)+" "+str(score)+" "+str(unique_identifier)+"\n"
417 
418  return dataset
419 
420 
421  #-------------------------------
422 
423 def get_cross_link_data(directory, filename, dist, omega, sigma,
424  don=None, doff=None, prior=0, type_of_profile="gofr"):
425 
426  (distmin, distmax, ndist) = dist
427  (omegamin, omegamax, nomega) = omega
428  (sigmamin, sigmamax, nsigma) = sigma
429 
430  filen = IMP.isd.get_data_path("CrossLinkPMFs.dict")
431  with open(filen) as xlpot:
432  dictionary = ast.literal_eval(xlpot.readline())
433 
434  xpot = dictionary[directory][filename]["distance"]
435  pot = dictionary[directory][filename][type_of_profile]
436 
437  dist_grid = get_grid(distmin, distmax, ndist, False)
438  omega_grid = get_log_grid(omegamin, omegamax, nomega)
439  sigma_grid = get_log_grid(sigmamin, sigmamax, nsigma)
440 
441  if not don is None and not doff is None:
442  xlmsdata = IMP.isd.CrossLinkData(
443  dist_grid,
444  omega_grid,
445  sigma_grid,
446  xpot,
447  pot,
448  don,
449  doff,
450  prior)
451  else:
452  xlmsdata = IMP.isd.CrossLinkData(
453  dist_grid,
454  omega_grid,
455  sigma_grid,
456  xpot,
457  pot)
458  return xlmsdata
459 
460  #-------------------------------
461 
462 
464  "If you use this function please let the PMI developers know.")
465 def get_cross_link_data_from_length(length, xxx_todo_changeme3, xxx_todo_changeme4, xxx_todo_changeme5):
466  (distmin, distmax, ndist) = xxx_todo_changeme3
467  (omegamin, omegamax, nomega) = xxx_todo_changeme4
468  (sigmamin, sigmamax, nsigma) = xxx_todo_changeme5
469 
470  dist_grid = get_grid(distmin, distmax, ndist, False)
471  omega_grid = get_log_grid(omegamin, omegamax, nomega)
472  sigma_grid = get_log_grid(sigmamin, sigmamax, nsigma)
473 
474  xlmsdata = IMP.isd.CrossLinkData(dist_grid, omega_grid, sigma_grid, length)
475  return xlmsdata
476 
477 
478 def get_grid(gmin, gmax, ngrid, boundaries):
479  grid = []
480  dx = (gmax - gmin) / float(ngrid)
481  for i in range(0, ngrid + 1):
482  if(not boundaries and i == 0):
483  continue
484  if(not boundaries and i == ngrid):
485  continue
486  grid.append(gmin + float(i) * dx)
487  return grid
488 
489  #-------------------------------
490 
491 
492 def get_log_grid(gmin, gmax, ngrid):
493  grid = []
494  for i in range(0, ngrid + 1):
495  grid.append(gmin * exp(float(i) / ngrid * log(gmax / gmin)))
496  return grid
497 
498  #-------------------------------
499 
500 
502  '''
503  example '"{ID_Score}" > 28 AND "{Sample}" ==
504  "%10_1%" OR ":Sample}" == "%10_2%" OR ":Sample}"
505  == "%10_3%" OR ":Sample}" == "%8_1%" OR ":Sample}" == "%8_2%"'
506  '''
507 
508  import pyparsing as pp
509 
510  operator = pp.Regex(">=|<=|!=|>|<|==|in").setName("operator")
511  value = pp.QuotedString(
512  '"') | pp.Regex(
513  r"[+-]?\d+(:?\.\d*)?(:?[eE][+-]?\d+)?")
514  identifier = pp.Word(pp.alphas, pp.alphanums + "_")
515  comparison_term = identifier | value
516  condition = pp.Group(comparison_term + operator + comparison_term)
517 
518  expr = pp.operatorPrecedence(condition, [
519  ("OR", 2, pp.opAssoc.LEFT, ),
520  ("AND", 2, pp.opAssoc.LEFT, ),
521  ])
522 
523  parsedstring = str(expr.parseString(inputstring)) \
524  .replace("[", "(") \
525  .replace("]", ")") \
526  .replace(",", " ") \
527  .replace("'", " ") \
528  .replace("%", "'") \
529  .replace("{", "float(entry['") \
530  .replace("}", "'])") \
531  .replace(":", "str(entry['") \
532  .replace("}", "'])") \
533  .replace("AND", "and") \
534  .replace("OR", "or")
535  return parsedstring
536 
537 
538 def open_file_or_inline_text(filename):
539  try:
540  fl = open(filename, "r")
541  except IOError:
542  fl = filename.split("\n")
543  return fl
544 
545 
547  "If you use this function please let the PMI developers know.")
548 def get_drmsd(prot0, prot1):
549  drmsd = 0.
550  npairs = 0.
551  for i in range(0, len(prot0) - 1):
552  for j in range(i + 1, len(prot0)):
553  dist0 = IMP.core.get_distance(prot0[i], prot0[j])
554  dist1 = IMP.core.get_distance(prot1[i], prot1[j])
555  drmsd += (dist0 - dist1) ** 2
556  npairs += 1.
557  return math.sqrt(drmsd / npairs)
558 
559  #-------------------------------
560 
561 
562 def get_ids_from_fasta_file(fastafile):
563  ids = []
564  with open(fastafile) as ff:
565  for l in ff:
566  if l[0] == ">":
567  ids.append(l[1:-1])
568  return ids
569 
570 
571 def get_closest_residue_position(hier, resindex, terminus="N"):
572  '''
573  this function works with plain hierarchies, as read from the pdb,
574  no multi-scale hierarchies
575  '''
576  p = []
577  niter = 0
578  while len(p) == 0:
579  niter += 1
580  sel = IMP.atom.Selection(hier, residue_index=resindex,
581  atom_type=IMP.atom.AT_CA)
582 
583  if terminus == "N":
584  resindex += 1
585  if terminus == "C":
586  resindex -= 1
587 
588  if niter >= 10000:
589  print("get_closest_residue_position: exiting while loop without result")
590  break
591  p = sel.get_selected_particles()
592 
593  if len(p) == 1:
594  return IMP.core.XYZ(p[0]).get_coordinates()
595  elif len(p) == 0:
596  print("get_closest_residue_position: got NO residues for hierarchy %s and residue %i" % (hier, resindex))
597  raise Exception("get_closest_residue_position: got NO residues for hierarchy %s and residue %i" % (
598  hier, resindex))
599  else:
600  raise ValueError("got multiple residues for hierarchy %s and residue %i; the list of particles is %s" % (hier, resindex, str([pp.get_name() for pp in p])))
601 
602 def get_terminal_residue(representation, hier, terminus="C", resolution=1):
603  '''
604  Get the particle of the terminal residue at the GIVEN resolution
605  (NOTE: not the closest resolution!).
606  To get the terminal residue at the closest resolution use:
607  particles=IMP.pmi.tools.select_by_tuple(representation,molecule_name)
608  particles[0] and particles[-1] will be the first and last particles
609  corresponding to the two termini.
610  It is needed for instance to determine the last residue of a pdb.
611  @param hier hierarchy containing the terminal residue
612  @param terminus either 'N' or 'C'
613  @param resolution resolution to use.
614  '''
615  termresidue = None
616  termparticle = None
617 
618  ps=select(representation,
619  resolution=resolution,
620  hierarchies=[hier])
621 
622  for p in ps:
623  if IMP.pmi.Resolution(p).get_resolution() == resolution:
625  if terminus == "C":
626  if termresidue is None:
627  termresidue = max(residues)
628  termparticle = p
629  elif max(residues) >= termresidue:
630  termresidue = max(residues)
631  termparticle = p
632  elif terminus == "N":
633  if termresidue is None:
634  termresidue = min(residues)
635  termparticle = p
636  elif min(residues) <= termresidue:
637  termresidue = min(residues)
638  termparticle = p
639  else:
640  raise ValueError("terminus argument should be either N or C")
641  return termparticle
642 
643 def get_terminal_residue_position(representation, hier, terminus="C",
644  resolution=1):
645  """Get XYZ coordinates of the terminal residue at the GIVEN resolution"""
646  p = get_terminal_residue(representation, hier, terminus, resolution)
647  return IMP.core.XYZ(p).get_coordinates()
648 
649 def get_residue_gaps_in_hierarchy(hierarchy, start, end):
650  '''
651  Return the residue index gaps and contiguous segments in the hierarchy.
652 
653  @param hierarchy hierarchy to examine
654  @param start first residue index
655  @param end last residue index
656 
657  @return A list of lists of the form
658  [[1,100,"cont"],[101,120,"gap"],[121,200,"cont"]]
659  '''
660  gaps = []
661  for n, rindex in enumerate(range(start, end + 1)):
662  sel = IMP.atom.Selection(hierarchy, residue_index=rindex,
663  atom_type=IMP.atom.AT_CA)
664 
665  if len(sel.get_selected_particles()) == 0:
666  if n == 0:
667  # set the initial condition
668  rindexgap = start
669  rindexcont = start - 1
670  if rindexgap == rindex - 1:
671  # residue is contiguous with the previously discovered gap
672  gaps[-1][1] += 1
673  else:
674  # residue is not contiguous with the previously discovered gap
675  # hence create a new gap tuple
676  gaps.append([rindex, rindex, "gap"])
677  # update the index of the last residue gap
678  rindexgap = rindex
679  else:
680  if n == 0:
681  # set the initial condition
682  rindexgap = start - 1
683  rindexcont = start
684  if rindexcont == rindex - 1:
685  # residue is contiguous with the previously discovered
686  # continuous part
687  gaps[-1][1] += 1
688  else:
689  # residue is not contiguous with the previously discovered continuous part
690  # hence create a new cont tuple
691  gaps.append([rindex, rindex, "cont"])
692  # update the index of the last residue gap
693  rindexcont = rindex
694  return gaps
695 
696 
697 class map(object):
698 
699  def __init__(self):
700  self.map = {}
701 
702  def set_map_element(self, xvalue, yvalue):
703  self.map[xvalue] = yvalue
704 
705  def get_map_element(self, invalue):
706  if type(invalue) == float:
707  n = 0
708  mindist = 1
709  for x in self.map:
710  dist = (invalue - x) * (invalue - x)
711 
712  if n == 0:
713  mindist = dist
714  minx = x
715  if dist < mindist:
716  mindist = dist
717  minx = x
718  n += 1
719  return self.map[minx]
720  elif type(invalue) == str:
721  return self.map[invalue]
722  else:
723  raise TypeError("wrong type for map")
724 
725 
726 def select(representation,
727  resolution=None,
728  hierarchies=None,
729  selection_arguments=None,
730  name=None,
731  name_is_ambiguous=False,
732  first_residue=None,
733  last_residue=None,
734  residue=None,
735  representation_type=None):
736  '''
737  this function uses representation=SimplifiedModel
738  it returns the corresponding selected particles
739  representation_type="Beads", "Res:X", "Densities", "Representation", "Molecule"
740  '''
741 
742  if resolution is None:
743  allparticles = IMP.atom.get_leaves(representation.prot)
744  resolution_particles = None
745  hierarchies_particles = None
746  names_particles = None
747  residue_range_particles = None
748  residue_particles = None
749  representation_type_particles = None
750 
751  if not resolution is None:
752  resolution_particles = []
753  hs = representation.get_hierarchies_at_given_resolution(resolution)
754  for h in hs:
755  resolution_particles += IMP.atom.get_leaves(h)
756 
757  if not hierarchies is None:
758  hierarchies_particles = []
759  for h in hierarchies:
760  hierarchies_particles += IMP.atom.get_leaves(h)
761 
762  if not name is None:
763  names_particles = []
764  if name_is_ambiguous:
765  for namekey in representation.hier_dict:
766  if name in namekey:
767  names_particles += IMP.atom.get_leaves(
768  representation.hier_dict[namekey])
769  elif name in representation.hier_dict:
770  names_particles += IMP.atom.get_leaves(representation.hier_dict[name])
771  else:
772  print("select: component %s is not there" % name)
773 
774  if not first_residue is None and not last_residue is None:
775  sel = IMP.atom.Selection(representation.prot,
776  residue_indexes=range(first_residue, last_residue + 1))
777  residue_range_particles = [IMP.atom.Hierarchy(p)
778  for p in sel.get_selected_particles()]
779 
780  if not residue is None:
781  sel = IMP.atom.Selection(representation.prot, residue_index=residue)
782  residue_particles = [IMP.atom.Hierarchy(p)
783  for p in sel.get_selected_particles()]
784 
785  if not representation_type is None:
786  representation_type_particles = []
787  if representation_type == "Molecule":
788  for name in representation.hier_representation:
789  for repr_type in representation.hier_representation[name]:
790  if repr_type == "Beads" or "Res:" in repr_type:
791  h = representation.hier_representation[name][repr_type]
792  representation_type_particles += IMP.atom.get_leaves(h)
793 
794  elif representation_type == "PDB":
795  for name in representation.hier_representation:
796  for repr_type in representation.hier_representation[name]:
797  if repr_type == "Res:" in repr_type:
798  h = representation.hier_representation[name][repr_type]
799  representation_type_particles += IMP.atom.get_leaves(h)
800 
801  else:
802  for name in representation.hier_representation:
803  h = representation.hier_representation[
804  name][
805  representation_type]
806  representation_type_particles += IMP.atom.get_leaves(h)
807 
808  selections = [hierarchies_particles, names_particles,
809  residue_range_particles, residue_particles, representation_type_particles]
810 
811  if resolution is None:
812  selected_particles = set(allparticles)
813  else:
814  selected_particles = set(resolution_particles)
815 
816  for s in selections:
817  if not s is None:
818  selected_particles = (set(s) & selected_particles)
819 
820  return list(selected_particles)
821 
822 
823 def select_by_tuple(
824  representation,
825  tupleselection,
826  resolution=None,
827  name_is_ambiguous=False):
828  if isinstance(tupleselection, tuple) and len(tupleselection) == 3:
829  particles = IMP.pmi.tools.select(representation, resolution=resolution,
830  name=tupleselection[2],
831  first_residue=tupleselection[0],
832  last_residue=tupleselection[1],
833  name_is_ambiguous=name_is_ambiguous)
834  elif isinstance(tupleselection, str):
835  particles = IMP.pmi.tools.select(representation, resolution=resolution,
836  name=tupleselection,
837  name_is_ambiguous=name_is_ambiguous)
838  else:
839  raise ValueError('you passed something bad to select_by_tuple()')
840  # now order the result by residue number
841  particles = IMP.pmi.tools.sort_by_residues(particles)
842 
843  return particles
844 
845 def select_by_tuple_2(hier,tuple_selection,resolution):
846  """New tuple format: molname OR (start,stop,molname,copynum,statenum)
847  Copy and state are optional. Can also use 'None' for them which will get all.
848  You can also pass -1 for stop which will go to the end.
849  Returns the particles
850  """
851  kwds = {} # going to accumulate keywords
852  kwds['resolution'] = resolution
853  if type(tuple_selection) is str:
854  kwds['molecule'] = tuple_selection
855  elif type(tuple_selection) is tuple:
856  rbegin = tuple_selection[0]
857  rend = tuple_selection[1]
858  kwds['molecule'] = tuple_selection[2]
859  try:
860  copynum = tuple_selection[3]
861  if copynum is not None:
862  kwds['copy_index'] = copynum
863  except:
864  pass
865  try:
866  statenum = tuple_selection[4]
867  if statenum is not None:
868  kwds['state_index'] = statenum
869  except:
870  pass
871  if rend==-1:
872  if rbegin>1:
873  s = IMP.atom.Selection(hier,**kwds)
874  s -= IMP.atom.Selection(hier,
875  residue_indexes=range(1,rbegin),
876  **kwds)
877  return s.get_selected_particles()
878  else:
879  kwds['residue_indexes'] = range(rbegin,rend+1)
880  s = IMP.atom.Selection(hier,**kwds)
881  return s.get_selected_particles()
882 
883 
884 
885 def get_db_from_csv(csvfilename):
886  import csv
887  outputlist = []
888  with open(csvfilename) as fh:
889  csvr = csv.DictReader(fh)
890  for l in csvr:
891  outputlist.append(l)
892  return outputlist
893 
894 
895 class HierarchyDatabase(object):
896  """Store the representations for a system."""
897 
898  def __init__(self):
899  self.db = {}
900  # this dictionary map a particle to its root hierarchy
901  self.root_hierarchy_dict = {}
902  self.preroot_fragment_hierarchy_dict = {}
903  self.particle_to_name = {}
904  self.model = None
905 
906  def add_name(self, name):
907  if name not in self.db:
908  self.db[name] = {}
909 
910  def add_residue_number(self, name, resn):
911  resn = int(resn)
912  self.add_name(name)
913  if resn not in self.db[name]:
914  self.db[name][resn] = {}
915 
916  def add_resolution(self, name, resn, resolution):
917  resn = int(resn)
918  resolution = float(resolution)
919  self.add_name(name)
920  self.add_residue_number(name, resn)
921  if resolution not in self.db[name][resn]:
922  self.db[name][resn][resolution] = []
923 
924  def add_particles(self, name, resn, resolution, particles):
925  resn = int(resn)
926  resolution = float(resolution)
927  self.add_name(name)
928  self.add_residue_number(name, resn)
929  self.add_resolution(name, resn, resolution)
930  self.db[name][resn][resolution] += particles
931  for p in particles:
932  (rh, prf) = self.get_root_hierarchy(p)
933  self.root_hierarchy_dict[p] = rh
934  self.preroot_fragment_hierarchy_dict[p] = prf
935  self.particle_to_name[p] = name
936  if self.model is None:
937  self.model = particles[0].get_model()
938 
939  def get_model(self):
940  return self.model
941 
942  def get_names(self):
943  names = list(self.db.keys())
944  names.sort()
945  return names
946 
947  def get_particles(self, name, resn, resolution):
948  resn = int(resn)
949  resolution = float(resolution)
950  return self.db[name][resn][resolution]
951 
952  def get_particles_at_closest_resolution(self, name, resn, resolution):
953  resn = int(resn)
954  resolution = float(resolution)
955  closestres = min(self.get_residue_resolutions(name, resn),
956  key=lambda x: abs(float(x) - float(resolution)))
957  return self.get_particles(name, resn, closestres)
958 
959  def get_residue_resolutions(self, name, resn):
960  resn = int(resn)
961  resolutions = list(self.db[name][resn].keys())
962  resolutions.sort()
963  return resolutions
964 
965  def get_molecule_resolutions(self, name):
966  resolutions = set()
967  for resn in self.db[name]:
968  resolutions.update(list(self.db[name][resn].keys()))
969  resolutions.sort()
970  return resolutions
971 
972  def get_residue_numbers(self, name):
973  residue_numbers = list(self.db[name].keys())
974  residue_numbers.sort()
975  return residue_numbers
976 
977  def get_particles_by_resolution(self, name, resolution):
978  resolution = float(resolution)
979  particles = []
980  for resn in self.get_residue_numbers(name):
981  result = self.get_particles_at_closest_resolution(
982  name,
983  resn,
984  resolution)
985  pstemp = [p for p in result if p not in particles]
986  particles += pstemp
987  return particles
988 
989  def get_all_particles_by_resolution(self, resolution):
990  resolution = float(resolution)
991  particles = []
992  for name in self.get_names():
993  particles += self.get_particles_by_resolution(name, resolution)
994  return particles
995 
996  def get_root_hierarchy(self, particle):
997  prerootfragment = particle
998  while IMP.atom.Atom.get_is_setup(particle) or \
999  IMP.atom.Residue.get_is_setup(particle) or \
1001  if IMP.atom.Atom.get_is_setup(particle):
1002  p = IMP.atom.Atom(particle).get_parent()
1003  elif IMP.atom.Residue.get_is_setup(particle):
1004  p = IMP.atom.Residue(particle).get_parent()
1005  elif IMP.atom.Fragment.get_is_setup(particle):
1006  p = IMP.atom.Fragment(particle).get_parent()
1007  prerootfragment = particle
1008  particle = p
1009  return (
1010  (IMP.atom.Hierarchy(particle), IMP.atom.Hierarchy(prerootfragment))
1011  )
1012 
1013  def get_all_root_hierarchies_by_resolution(self, resolution):
1014  hierarchies = []
1015  resolution = float(resolution)
1016  particles = self.get_all_particles_by_resolution(resolution)
1017  for p in particles:
1018  rh = self.root_hierarchy_dict[p]
1019  if rh not in hierarchies:
1020  hierarchies.append(IMP.atom.Hierarchy(rh))
1021  return hierarchies
1022 
1023  def get_preroot_fragments_by_resolution(self, name, resolution):
1024  fragments = []
1025  resolution = float(resolution)
1026  particles = self.get_particles_by_resolution(name, resolution)
1027  for p in particles:
1028  fr = self.preroot_fragment_hierarchy_dict[p]
1029  if fr not in fragments:
1030  fragments.append(fr)
1031  return fragments
1032 
1033  def show(self, name):
1034  print(name)
1035  for resn in self.get_residue_numbers(name):
1036  print(resn)
1037  for resolution in self.get_residue_resolutions(name, resn):
1038  print("----", resolution)
1039  for p in self.get_particles(name, resn, resolution):
1040  print("--------", p.get_name())
1041 
1042 
1043 def get_prot_name_from_particle(p, list_of_names):
1044  '''Return the component name provided a particle and a list of names'''
1045  root = p
1046  protname = root.get_name()
1047  is_a_bead = False
1048  while not protname in list_of_names:
1049  root0 = root.get_parent()
1050  if root0 == IMP.atom.Hierarchy():
1051  return (None, None)
1052  protname = root0.get_name()
1053 
1054  # check if that is a bead
1055  # this piece of code might be dangerous if
1056  # the hierarchy was called Bead :)
1057  if "Beads" in protname:
1058  is_a_bead = True
1059  root = root0
1060  return (protname, is_a_bead)
1061 
1062 
1064  '''
1065  Retrieve the residue indexes for the given particle.
1066 
1067  The particle must be an instance of Fragment,Residue, Atom or Molecule
1068  or else returns an empty list
1069  '''
1070  resind = []
1072  resind = IMP.atom.Fragment(hier).get_residue_indexes()
1073  elif IMP.atom.Residue.get_is_setup(hier):
1074  resind = [IMP.atom.Residue(hier).get_index()]
1075  elif IMP.atom.Atom.get_is_setup(hier):
1076  a = IMP.atom.Atom(hier)
1077  resind = [IMP.atom.Residue(a.get_parent()).get_index()]
1078  elif IMP.atom.Molecule.get_is_setup(hier):
1079  resind_tmp=IMP.pmi.tools.OrderedSet()
1080  for lv in IMP.atom.get_leaves(hier):
1081  if IMP.atom.Fragment.get_is_setup(lv) or \
1084  for ind in get_residue_indexes(lv): resind_tmp.add(ind)
1085  resind=list(resind_tmp)
1086  else:
1087  resind = []
1088  return resind
1089 
1090 
1091 def sort_by_residues(particles):
1092  particles_residues = [(p, IMP.pmi.tools.get_residue_indexes(p))
1093  for p in particles]
1094  sorted_particles_residues = sorted(
1095  particles_residues,
1096  key=lambda tup: tup[1])
1097  particles = [p[0] for p in sorted_particles_residues]
1098  return particles
1099 
1100 
1102  "If you use this function please let the PMI developers know.")
1103 def get_residue_to_particle_map(particles):
1104  # this function returns a dictionary that map particles to residue indexes
1105  particles = sort_by_residues(particles)
1106  particles_residues = [(p, IMP.pmi.tools.get_residue_indexes(p))
1107  for p in particles]
1108  return dict(zip(particles_residues, particles))
1109 
1110 #
1111 # Parallel Computation
1112 #
1113 
1114 
1116  """Synchronize data over a parallel run"""
1117  from mpi4py import MPI
1118  comm = MPI.COMM_WORLD
1119  rank = comm.Get_rank()
1120  number_of_processes = comm.size
1121  comm.Barrier()
1122  if rank != 0:
1123  comm.send(data, dest=0, tag=11)
1124 
1125  elif rank == 0:
1126  for i in range(1, number_of_processes):
1127  data_tmp = comm.recv(source=i, tag=11)
1128  if type(data) == list:
1129  data += data_tmp
1130  elif type(data) == dict:
1131  data.update(data_tmp)
1132  else:
1133  raise TypeError("data not supported, use list or dictionaries")
1134 
1135  for i in range(1, number_of_processes):
1136  comm.send(data, dest=i, tag=11)
1137 
1138  if rank != 0:
1139  data = comm.recv(source=0, tag=11)
1140  return data
1141 
1143  "If you use this function please let the PMI developers know.")
1144 def scatter_and_gather_dict_append(data):
1145  """Synchronize data over a parallel run"""
1146  from mpi4py import MPI
1147  comm = MPI.COMM_WORLD
1148  rank = comm.Get_rank()
1149  number_of_processes = comm.size
1150  comm.Barrier()
1151  if rank != 0:
1152  comm.send(data, dest=0, tag=11)
1153  elif rank == 0:
1154  for i in range(1, number_of_processes):
1155  data_tmp = comm.recv(source=i, tag=11)
1156  for k in data:
1157  data[k]+=data_tmp[k]
1158 
1159  for i in range(1, number_of_processes):
1160  comm.send(data, dest=i, tag=11)
1161 
1162  if rank != 0:
1163  data = comm.recv(source=0, tag=11)
1164  return data
1165 
1166 
1167 #
1168 ### Lists and iterators
1169 #
1170 
1171 def sublist_iterator(l, lmin=None, lmax=None):
1172  '''
1173  Yield all sublists of length >= lmin and <= lmax
1174  '''
1175  if lmin is None:
1176  lmin = 0
1177  if lmax is None:
1178  lmax = len(l)
1179  n = len(l) + 1
1180  for i in range(n):
1181  for j in range(i + 1, n):
1182  if len(l[i:j]) <= lmax and len(l[i:j]) >= lmin:
1183  yield l[i:j]
1184 
1185 
1186 def flatten_list(l):
1187  return [item for sublist in l for item in sublist]
1188 
1189 
1190 def list_chunks_iterator(list, length):
1191  """ Yield successive length-sized chunks from a list.
1192  """
1193  for i in range(0, len(list), length):
1194  yield list[i:i + length]
1195 
1196 
1197 def chunk_list_into_segments(seq, num):
1198  seq = list(seq)
1199  avg = len(seq) / float(num)
1200  out = []
1201  last = 0.0
1202 
1203  while last < len(seq):
1204  out.append(seq[int(last):int(last + avg)])
1205  last += avg
1206 
1207  return out
1208 
1209 
1210 class Segments(object):
1211 
1212  ''' This class stores integers
1213  in ordered compact lists eg:
1214  [[1,2,3],[6,7,8]]
1215  the methods help splitting and merging the internal lists
1216  Example:
1217  s=Segments([1,2,3]) is [[1,2,3]]
1218  s.add(4) is [[1,2,3,4]] (add right)
1219  s.add(3) is [[1,2,3,4]] (item already existing)
1220  s.add(7) is [[1,2,3,4],[7]] (new list)
1221  s.add([8,9]) is [[1,2,3,4],[7,8,9]] (add item right)
1222  s.add([5,6]) is [[1,2,3,4,5,6,7,8,9]] (merge)
1223  s.remove(3) is [[1,2],[4,5,6,7,8,9]] (split)
1224  etc.
1225  '''
1226 
1227  def __init__(self,index):
1228  '''index can be a integer or a list of integers '''
1229  if type(index) is int:
1230  self.segs=[[index]]
1231  elif type(index) is list:
1232  self.segs=[[index[0]]]
1233  for i in index[1:]:
1234  self.add(i)
1235  else:
1236  raise TypeError("index must be an int or list of ints")
1237 
1238  def add(self,index):
1239  '''index can be a integer or a list of integers '''
1240  if type(index) is int:
1241  mergeleft=None
1242  mergeright=None
1243  for n,s in enumerate(self.segs):
1244  if index in s:
1245  return 0
1246  else:
1247  if s[0]-index==1:
1248  mergeleft=n
1249  if index-s[-1]==1:
1250  mergeright=n
1251  if mergeright is None and mergeleft is None:
1252  self.segs.append([index])
1253  if not mergeright is None and mergeleft is None:
1254  self.segs[mergeright].append(index)
1255  if not mergeleft is None and mergeright is None:
1256  self.segs[mergeleft]=[index]+self.segs[mergeleft]
1257  if not mergeleft is None and not mergeright is None:
1258  self.segs[mergeright]=self.segs[mergeright]+[index]+self.segs[mergeleft]
1259  del self.segs[mergeleft]
1260 
1261  for n in range(len(self.segs)):
1262  self.segs[n].sort()
1263 
1264  self.segs.sort(key=lambda tup: tup[0])
1265 
1266  elif type(index) is list:
1267  for i in index:
1268  self.add(i)
1269  else:
1270  raise TypeError("index must be an int or list of ints")
1271 
1272  def remove(self,index):
1273  '''index can be a integer'''
1274  for n,s in enumerate(self.segs):
1275  if index in s:
1276  if s[0]==index:
1277  self.segs[n]=s[1:]
1278  elif s[-1]==index:
1279  self.segs[n]=s[:-1]
1280  else:
1281  i=self.segs[n].index(index)
1282  self.segs[n]=s[:i]
1283  self.segs.append(s[i+1:])
1284  for n in range(len(self.segs)):
1285  self.segs[n].sort()
1286  if len(self.segs[n])==0:
1287  del self.segs[n]
1288  self.segs.sort(key=lambda tup: tup[0])
1289 
1290  def get_flatten(self):
1291  ''' Returns a flatten list '''
1292  return [item for sublist in self.segs for item in sublist]
1293 
1294  def __repr__(self):
1295  ret_tmp="["
1296  for seg in self.segs:
1297  ret_tmp+=str(seg[0])+"-"+str(seg[-1])+","
1298  ret=ret_tmp[:-1]+"]"
1299  return ret
1300 
1301 
1302 
1303 #
1304 # COORDINATE MANIPULATION
1305 #
1306 
1307 
1309  "If you use this function please let the PMI developers know.")
1310 def translate_hierarchy(hierarchy, translation_vector):
1311  '''
1312  Apply a translation to a hierarchy along the input vector.
1313  '''
1314  rbs = set()
1315  xyzs = set()
1316  if type(translation_vector) == list:
1317  transformation = IMP.algebra.Transformation3D(
1318  IMP.algebra.Vector3D(translation_vector))
1319  else:
1320  transformation = IMP.algebra.Transformation3D(translation_vector)
1321  for p in IMP.atom.get_leaves(hierarchy):
1323  rbs.add(IMP.core.RigidBody(p))
1325  rb = IMP.core.RigidMember(p).get_rigid_body()
1326  rbs.add(rb)
1327  else:
1328  xyzs.add(p)
1329  for xyz in xyzs:
1330  IMP.core.transform(IMP.core.XYZ(xyz), transformation)
1331  for rb in rbs:
1332  IMP.core.transform(rb, transformation)
1333 
1334 
1336  "If you use this function please let the PMI developers know.")
1337 def translate_hierarchies(hierarchies, translation_vector):
1338  for h in hierarchies:
1339  IMP.pmi.tools.translate_hierarchy(h, translation_vector)
1340 
1341 
1343  "If you use this function please let the PMI developers know.")
1344 def translate_hierarchies_to_reference_frame(hierarchies):
1345  xc = 0
1346  yc = 0
1347  zc = 0
1348  nc = 0
1349  for h in hierarchies:
1350  for p in IMP.atom.get_leaves(h):
1351  coor = IMP.core.XYZ(p).get_coordinates()
1352  nc += 1
1353  xc += coor[0]
1354  yc += coor[1]
1355  zc += coor[2]
1356  xc = xc / nc
1357  yc = yc / nc
1358  zc = zc / nc
1359  IMP.pmi.tools.translate_hierarchies(hierarchies, (-xc, -yc, -zc))
1360 
1361 
1362 #
1363 # Tools to simulate data
1364 #
1365 
1366 def normal_density_function(expected_value, sigma, x):
1367  return (
1368  1 / math.sqrt(2 * math.pi) / sigma *
1369  math.exp(-(x - expected_value) ** 2 / 2 / sigma / sigma)
1370  )
1371 
1372 
1373 def log_normal_density_function(expected_value, sigma, x):
1374  return (
1375  1 / math.sqrt(2 * math.pi) / sigma / x *
1376  math.exp(-(math.log(x / expected_value) ** 2 / 2 / sigma / sigma))
1377  )
1378 
1379 
1380 def get_random_residue_pairs(representation, resolution,
1381  number,
1382  max_distance=None,
1383  avoid_same_particles=False,
1384  names=None):
1385 
1386  particles = []
1387  if names is None:
1388  names=list(representation.hier_dict.keys())
1389 
1390  for name in names:
1391  prot = representation.hier_dict[name]
1392  particles += select(representation,name=name,resolution=resolution)
1393  random_residue_pairs = []
1394  while len(random_residue_pairs)<=number:
1395  p1 = random.choice(particles)
1396  p2 = random.choice(particles)
1397  if max_distance is not None and \
1398  core.get_distance(core.XYZ(p1), core.XYZ(p2)) > max_distance:
1399  continue
1400  r1 = random.choice(IMP.pmi.tools.get_residue_indexes(p1))
1401  r2 = random.choice(IMP.pmi.tools.get_residue_indexes(p2))
1402  if r1==r2 and avoid_same_particles: continue
1403  name1 = representation.get_prot_name_from_particle(p1)
1404  name2 = representation.get_prot_name_from_particle(p2)
1405  random_residue_pairs.append((name1, r1, name2, r2))
1406 
1407  return random_residue_pairs
1408 
1409 
1411  "If you use this function please let the PMI developers know.")
1412 def get_random_data_point(
1413  expected_value,
1414  ntrials,
1415  sensitivity,
1416  sigma,
1417  outlierprob,
1418  begin_end_nbins_tuple,
1419  log=False,
1420  loggrid=False):
1421  import random
1422 
1423  begin = begin_end_nbins_tuple[0]
1424  end = begin_end_nbins_tuple[1]
1425  nbins = begin_end_nbins_tuple[2]
1426 
1427  if not loggrid:
1428  fmod_grid = get_grid(begin, end, nbins, True)
1429  else:
1430  fmod_grid = get_log_grid(begin, end, nbins)
1431 
1432  norm = 0
1433  cumul = []
1434  cumul.append(0)
1435 
1436  a = []
1437  for i in range(0, ntrials):
1438  a.append([random.random(), True])
1439 
1440  if sigma != 0.0:
1441  for j in range(1, len(fmod_grid)):
1442  fj = fmod_grid[j]
1443  fjm1 = fmod_grid[j - 1]
1444  df = fj - fjm1
1445 
1446  if not log:
1447  pj = normal_density_function(expected_value, sigma, fj)
1448  pjm1 = normal_density_function(expected_value, sigma, fjm1)
1449  else:
1450  pj = log_normal_density_function(expected_value, sigma, fj)
1451  pjm1 = log_normal_density_function(expected_value, sigma, fjm1)
1452 
1453  norm += (pj + pjm1) / 2.0 * df
1454  cumul.append(norm)
1455  # print fj, pj
1456 
1457  random_points = []
1458 
1459  for i in range(len(cumul)):
1460  # print i,a, cumul[i], norm
1461  for aa in a:
1462  if (aa[0] <= cumul[i] / norm and aa[1]):
1463  random_points.append(
1464  int(fmod_grid[i] / sensitivity) * sensitivity)
1465  aa[1] = False
1466 
1467  else:
1468  random_points = [expected_value] * ntrials
1469 
1470  for i in range(len(random_points)):
1471  if random.random() < outlierprob:
1472  a = random.uniform(begin, end)
1473  random_points[i] = int(a / sensitivity) * sensitivity
1474  print(random_points)
1475  '''
1476  for i in range(ntrials):
1477  if random.random() > OUTLIERPROB_:
1478  r=truncnorm.rvs(0.0,1.0,expected_value,BETA_)
1479  if r>1.0: print r,expected_value,BETA_
1480  else:
1481  r=random.random()
1482  random_points.append(int(r/sensitivity)*sensitivity)
1483  '''
1484 
1485  rmean = 0.
1486  rmean2 = 0.
1487  for r in random_points:
1488  rmean += r
1489  rmean2 += r * r
1490 
1491  rmean /= float(ntrials)
1492  rmean2 /= float(ntrials)
1493  stddev = math.sqrt(max(rmean2 - rmean * rmean, 0.))
1494  return rmean, stddev
1495 
1496 def print_multicolumn(list_of_strings, ncolumns=2, truncate=40):
1497 
1498  l = list_of_strings
1499 
1500  cols = ncolumns
1501  # add empty entries after l
1502  for i in range(len(l) % cols):
1503  l.append(" ")
1504 
1505  split = [l[i:i + len(l) / cols] for i in range(0, len(l), len(l) / cols)]
1506  for row in zip(*split):
1507  print("".join(str.ljust(i, truncate) for i in row))
1508 
1509 class ColorChange(object):
1510  '''Change color code to hexadecimal to rgb'''
1511  def __init__(self):
1512  self._NUMERALS = '0123456789abcdefABCDEF'
1513  self._HEXDEC = dict((v, int(v, 16)) for v in (x+y for x in self._NUMERALS for y in self._NUMERALS))
1514  self.LOWERCASE, self.UPPERCASE = 'x', 'X'
1515 
1516  def rgb(self,triplet):
1517  return float(self._HEXDEC[triplet[0:2]]), float(self._HEXDEC[triplet[2:4]]), float(self._HEXDEC[triplet[4:6]])
1518 
1519  def triplet(self,rgb, lettercase=None):
1520  if lettercase is None: lettercase=self.LOWERCASE
1521  return format(rgb[0]<<16 | rgb[1]<<8 | rgb[2], '06'+lettercase)
1522 
1523 # -------------- Collections --------------- #
1524 
1525 class OrderedSet(collections.MutableSet):
1526 
1527  def __init__(self, iterable=None):
1528  self.end = end = []
1529  end += [None, end, end] # sentinel node for doubly linked list
1530  self.map = {} # key --> [key, prev, next]
1531  if iterable is not None:
1532  self |= iterable
1533 
1534  def __len__(self):
1535  return len(self.map)
1536 
1537  def __contains__(self, key):
1538  return key in self.map
1539 
1540  def add(self, key):
1541  if key not in self.map:
1542  end = self.end
1543  curr = end[1]
1544  curr[2] = end[1] = self.map[key] = [key, curr, end]
1545 
1546  def discard(self, key):
1547  if key in self.map:
1548  key, prev, next = self.map.pop(key)
1549  prev[2] = next
1550  next[1] = prev
1551 
1552  def __iter__(self):
1553  end = self.end
1554  curr = end[2]
1555  while curr is not end:
1556  yield curr[0]
1557  curr = curr[2]
1558 
1559  def __reversed__(self):
1560  end = self.end
1561  curr = end[1]
1562  while curr is not end:
1563  yield curr[0]
1564  curr = curr[1]
1565 
1566  def pop(self, last=True):
1567  if not self:
1568  raise KeyError('set is empty')
1569  if last:
1570  key = self.end[1][0]
1571  else:
1572  key = self.end[2][0]
1573  self.discard(key)
1574  return key
1575 
1576  def __repr__(self):
1577  if not self:
1578  return '%s()' % (self.__class__.__name__,)
1579  return '%s(%r)' % (self.__class__.__name__, list(self))
1580 
1581  def __eq__(self, other):
1582  if isinstance(other, OrderedSet):
1583  return len(self) == len(other) and list(self) == list(other)
1584  return set(self) == set(other)
1585 
1586 
1587 class OrderedDefaultDict(OrderedDict):
1588  """Store objects in order they were added, but with default type.
1589  Source: http://stackoverflow.com/a/4127426/2608793
1590  """
1591  def __init__(self, *args, **kwargs):
1592  if not args:
1593  self.default_factory = None
1594  else:
1595  if not (args[0] is None or callable(args[0])):
1596  raise TypeError('first argument must be callable or None')
1597  self.default_factory = args[0]
1598  args = args[1:]
1599  super(OrderedDefaultDict, self).__init__(*args, **kwargs)
1600 
1601  def __missing__ (self, key):
1602  if self.default_factory is None:
1603  raise KeyError(key)
1604  self[key] = default = self.default_factory()
1605  return default
1606 
1607  def __reduce__(self): # optional, for pickle support
1608  args = (self.default_factory,) if self.default_factory else ()
1609  if sys.version_info[0] >= 3:
1610  return self.__class__, args, None, None, self.items()
1611  else:
1612  return self.__class__, args, None, None, self.iteritems()
1613 
1614 
1615 # -------------- PMI2 Tools --------------- #
1616 
1617 def set_coordinates_from_rmf(hier,rmf_fn,frame_num=0):
1618  """Extract frame from RMF file and fill coordinates. Must be identical topology.
1619  @param hier The (System) hierarchy to fill (e.g. after you've built it)
1620  @param rmf_fn The file to extract from
1621  @param frame_num The frame number to extract
1622  """
1623  rh = RMF.open_rmf_file_read_only(rmf_fn)
1624  IMP.rmf.link_hierarchies(rh,[hier])
1625  IMP.rmf.load_frame(rh, RMF.FrameID(frame_num))
1626  del rh
1627 
1628 def input_adaptor(stuff,
1629  pmi_resolution=0,
1630  flatten=False,
1631  selection_tuple=None,
1632  warn_about_slices=True):
1633  """Adapt things for PMI (degrees of freedom, restraints, ...)
1634  Returns list of list of hierarchies, separated into Molecules if possible.
1635  The input can be a list, or a list of lists (iterable of ^1 or iterable of ^2)
1636  (iterable of ^2) Hierarchy -> returns input as list of list of hierarchies,
1637  only one entry, not grouped by molecules.
1638  (iterable of ^2) PMI::System/State/Molecule/TempResidue ->
1639  returns residue hierarchies, grouped in molecules, at requested resolution
1640  @param stuff Can be one of the following inputs:
1641  IMP Hierarchy, PMI System/State/Molecule/TempResidue, or a list/set (of list/set) of them.
1642  Must be uniform input, however. No mixing object types.
1643  @param pmi_resolution For selecting, only does it if you pass PMI objects. Set it to "all"
1644  if you want all resolutions!
1645  @param flatten Set to True if you just want all hierarchies in one list.
1646  @param warn_about_slices Print a warning if you are requesting only part
1647  of a bead. Sometimes you just don't care!
1648  @note since this relies on IMP::atom::Selection, this will not return
1649  any objects if they weren't built! But there should be no problem
1650  if you request unbuilt residues - they should be ignored.
1651  """
1652 
1653  if stuff is None:
1654  return stuff
1655 
1656  if hasattr(stuff,'__iter__'):
1657  if len(stuff)==0:
1658  return stuff
1659  thelist=list(stuff)
1660 
1661  # iter of iter of should be ok
1662  if all(hasattr(el,'__iter__') for el in thelist):
1663  thelist = [i for sublist in thelist for i in sublist]
1664  elif any(hasattr(el,'__iter__') for el in thelist):
1665  raise Exception('input_adaptor: input_object must be a list or a list of lists')
1666 
1667  stuff = thelist
1668  else:
1669  stuff = [stuff]
1670 
1671  # check that it is a hierarchy homogenously:
1672  try:
1673  is_hierarchy=all(IMP.atom.Hierarchy.get_is_setup(s) for s in stuff)
1674  except NotImplementedError:
1675  is_hierarchy=False
1676  # get the other types homogenously
1677  is_system=all(isinstance(s, IMP.pmi.topology.System) for s in stuff)
1678  is_state=all(isinstance(s, IMP.pmi.topology.State) for s in stuff)
1679  is_molecule=all(isinstance(s, IMP.pmi.topology.Molecule) for s in stuff)
1680  is_temp_residue=all(isinstance(s, IMP.pmi.topology.TempResidue) for s in stuff)
1681 
1682  # now that things are ok, do selection if requested
1683  hier_list = []
1684  pmi_input = False
1685  if is_system or is_state or is_molecule or is_temp_residue:
1686  # if PMI, perform selection using gathered indexes
1687  pmi_input = True
1688  indexes_per_mol = OrderedDefaultDict(list) #key is Molecule object, value are residues
1689  if is_system:
1690  for system in stuff:
1691  for state in system.get_states():
1692  mdict = state.get_molecules()
1693  for molname in mdict:
1694  for copy in mdict[molname]:
1695  indexes_per_mol[copy] += [r.get_index() for r in copy.get_residues()]
1696  elif is_state:
1697  for state in stuff:
1698  mdict = state.get_molecules()
1699  for molname in mdict:
1700  for copy in mdict[molname]:
1701  indexes_per_mol[copy] += [r.get_index() for r in copy.get_residues()]
1702  elif is_molecule:
1703  for molecule in stuff:
1704  indexes_per_mol[molecule] += [r.get_index() for r in molecule.get_residues()]
1705  else: # is_temp_residue
1706  for tempres in stuff:
1707  indexes_per_mol[tempres.get_molecule()].append(tempres.get_index())
1708  for mol in indexes_per_mol:
1709  if pmi_resolution=='all':
1710  # because you select from the molecule,
1711  # this will start the search from the base resolution
1712  ps = select_at_all_resolutions(mol.get_hierarchy(),
1713  residue_indexes=indexes_per_mol[mol])
1714  else:
1715  sel = IMP.atom.Selection(mol.get_hierarchy(),
1716  resolution=pmi_resolution,
1717  residue_indexes=indexes_per_mol[mol])
1718  ps = sel.get_selected_particles()
1719 
1720  # check that you don't have any incomplete fragments!
1721  if warn_about_slices:
1722  rset = set(indexes_per_mol[mol])
1723  for p in ps:
1725  fset = set(IMP.atom.Fragment(p).get_residue_indexes())
1726  if not fset <= rset:
1727  minset = min(fset)
1728  maxset = max(fset)
1729  found = fset&rset
1730  minf = min(found)
1731  maxf = max(found)
1732  resbreak = maxf if minf==minset else minset-1
1733  print('WARNING: You are trying to select only part of the bead %s:%i-%i.\n'
1734  'The residues you requested are %i-%i. You can fix this by:\n'
1735  '1) requesting the whole bead/none of it or\n'
1736  '2) break the bead up by passing bead_extra_breaks=[\'%i\'] in '
1737  'molecule.add_representation()'
1738  %(mol.get_name(),minset,maxset,minf,maxf,resbreak))
1739  hier_list.append([IMP.atom.Hierarchy(p) for p in ps])
1740  elif is_hierarchy:
1741  #check
1742  ps=[]
1743  if pmi_resolution=='all':
1744  for h in stuff:
1745  ps+=select_at_all_resolutions(h)
1746  else:
1747  for h in stuff:
1748  ps+=IMP.atom.Selection(h,resolution=pmi_resolution).get_selected_particles()
1749  hier_list=[IMP.atom.Hierarchy(p) for p in ps]
1750  if not flatten:
1751  hier_list = [hier_list]
1752  else:
1753  raise Exception('input_adaptor: you passed something of wrong type or a list with mixed types')
1754 
1755  if flatten and pmi_input:
1756  return [h for sublist in hier_list for h in sublist]
1757  else:
1758  return hier_list
1759 
1760 
1761 def get_sorted_segments(mol):
1762  """Returns sequence-sorted segments array, each containing the first particle
1763  the last particle and the first residue index."""
1764 
1765  from operator import itemgetter
1766  hiers=IMP.pmi.tools.input_adaptor(mol)
1767  if len(hiers)>1:
1768  raise ValueError("only pass stuff from one Molecule, please")
1769  hiers = hiers[0]
1770  segs = []
1771  for h in hiers:
1772  try:
1773  start = IMP.atom.Hierarchy(h).get_children()[0]
1774  except IndexError:
1775  start = IMP.atom.Hierarchy(h)
1776 
1777  try:
1778  end = IMP.atom.Hierarchy(h).get_children()[-1]
1779  except IndexError:
1780  end = IMP.atom.Hierarchy(h)
1781 
1782  startres = IMP.pmi.tools.get_residue_indexes(start)[0]
1783  endres = IMP.pmi.tools.get_residue_indexes(end)[-1]
1784  segs.append((start, end, startres))
1785  return sorted(segs, key=itemgetter(2))
1786 
1787 def display_bonds(mol):
1788  """Decorate the sequence-consecutive particles from a PMI2 molecule with a bond,
1789  so that they appear connected in the rmf file"""
1790  SortedSegments=get_sorted_segments(mol)
1791  for x in range(len(SortedSegments) - 1):
1792 
1793  last = SortedSegments[x][1]
1794  first = SortedSegments[x + 1][0]
1795 
1796  p1 = last.get_particle()
1797  p2 = first.get_particle()
1798  if not IMP.atom.Bonded.get_is_setup(p1):
1800  if not IMP.atom.Bonded.get_is_setup(p2):
1802 
1805  IMP.atom.Bonded(p1),
1806  IMP.atom.Bonded(p2),1)
1807 
1808 
1809 class ThreeToOneConverter(defaultdict):
1810  """This class converts three to one letter codes, and return X for any unknown codes"""
1811  def __init__(self,is_nucleic=False):
1812 
1813  if not is_nucleic:
1814  threetoone = {'ALA': 'A', 'ARG': 'R', 'ASN': 'N', 'ASP': 'D',
1815  'CYS': 'C', 'GLU': 'E', 'GLN': 'Q', 'GLY': 'G',
1816  'HIS': 'H', 'ILE': 'I', 'LEU': 'L', 'LYS': 'K',
1817  'MET': 'M', 'PHE': 'F', 'PRO': 'P', 'SER': 'S',
1818  'THR': 'T', 'TRP': 'W', 'TYR': 'Y', 'VAL': 'V', 'UNK': 'X'}
1819  else:
1820  threetoone = {'ADE': 'A', 'URA': 'U', 'CYT': 'C', 'GUA': 'G',
1821  'THY': 'T', 'UNK': 'X'}
1822 
1823  defaultdict.__init__(self,lambda: "X", threetoone)
1824 
1825 
1826 
1827 
1828 def get_residue_type_from_one_letter_code(code,is_nucleic=False):
1829  threetoone=ThreeToOneConverter(is_nucleic)
1830  one_to_three={}
1831  for k in threetoone:
1832  one_to_three[threetoone[k]] = k
1833  return IMP.atom.ResidueType(one_to_three[code])
1834 
1835 
1836 def get_all_leaves(list_of_hs):
1837  """ Just get the leaves from a list of hierarchies """
1838  lvs = list(itertools.chain.from_iterable(IMP.atom.get_leaves(item) for item in list_of_hs))
1839  return lvs
1840 
1841 
1842 def select_at_all_resolutions(hier=None,
1843  hiers=None,
1844  **kwargs):
1845  """Perform selection using the usual keywords but return ALL resolutions (BEADS and GAUSSIANS).
1846  Returns in flat list!
1847  """
1848 
1849  if hiers is None:
1850  hiers = []
1851  if hier is not None:
1852  hiers.append(hier)
1853  if len(hiers)==0:
1854  print("WARNING: You passed nothing to select_at_all_resolutions()")
1855  return []
1856  ret = OrderedSet()
1857  for hsel in hiers:
1858  try:
1859  htest = IMP.atom.Hierarchy.get_is_setup(hsel)
1860  except:
1861  raise Exception('select_at_all_resolutions: you have to pass an IMP Hierarchy')
1862  if not htest:
1863  raise Exception('select_at_all_resolutions: you have to pass an IMP Hierarchy')
1864  if 'resolution' in kwargs or 'representation_type' in kwargs:
1865  raise Exception("don't pass resolution or representation_type to this function")
1866  selB = IMP.atom.Selection(hsel,resolution=IMP.atom.ALL_RESOLUTIONS,
1867  representation_type=IMP.atom.BALLS,**kwargs)
1868  selD = IMP.atom.Selection(hsel,resolution=IMP.atom.ALL_RESOLUTIONS,
1869  representation_type=IMP.atom.DENSITIES,**kwargs)
1870  ret |= OrderedSet(selB.get_selected_particles())
1871  ret |= OrderedSet(selD.get_selected_particles())
1872  return list(ret)
1873 
1874 
1875 def get_particles_within_zone(hier,
1876  target_ps,
1877  sel_zone,
1878  entire_residues,
1879  exclude_backbone):
1880  """Utility to retrieve particles from a hierarchy within a
1881  zone around a set of ps.
1882  @param hier The hierarchy in which to look for neighbors
1883  @param target_ps The particles for zoning
1884  @param sel_zone The maximum distance
1885  @param entire_residues If True, will grab entire residues
1886  @param exclude_backbone If True, will only return sidechain particles
1887  """
1888 
1889  test_sel = IMP.atom.Selection(hier)
1890  backbone_types=['C','N','CB','O']
1891  if exclude_backbone:
1892  test_sel -= IMP.atom.Selection(hier,atom_types=[IMP.atom.AtomType(n)
1893  for n in backbone_types])
1894  test_ps = test_sel.get_selected_particles()
1895  nn = IMP.algebra.NearestNeighbor3D([IMP.core.XYZ(p).get_coordinates()
1896  for p in test_ps])
1897  zone = set()
1898  for target in target_ps:
1899  zone|=set(nn.get_in_ball(IMP.core.XYZ(target).get_coordinates(),sel_zone))
1900  zone_ps = [test_ps[z] for z in zone]
1901  if entire_residues:
1902  final_ps = set()
1903  for z in zone_ps:
1904  final_ps|=set(IMP.atom.Hierarchy(z).get_parent().get_children())
1905  zone_ps = [h.get_particle() for h in final_ps]
1906  return zone_ps
1907 
1908 
1909 def get_rbs_and_beads(hiers):
1910  """Returns unique objects in original order"""
1911  rbs = set()
1912  beads = []
1913  rbs_ordered = []
1914  if not hasattr(hiers,'__iter__'):
1915  hiers = [hiers]
1916  for p in get_all_leaves(hiers):
1918  rb = IMP.core.RigidMember(p).get_rigid_body()
1919  if rb not in rbs:
1920  rbs.add(rb)
1921  rbs_ordered.append(rb)
1923  rb = IMP.core.NonRigidMember(p).get_rigid_body()
1924  if rb not in rbs:
1925  rbs.add(rb)
1926  rbs_ordered.append(rb)
1927  beads.append(p)
1928  else:
1929  beads.append(p)
1930  return rbs_ordered,beads
1931 
1932 def get_molecules(input_objects):
1933  "This function returns the parent molecule hierarchies of given objects"
1934  stuff=input_adaptor(input_objects, pmi_resolution='all',flatten=True)
1935  molecules=set()
1936  for h in stuff:
1937  is_root=False
1938  is_molecule=False
1939  while not (is_molecule or is_root):
1940  root=IMP.atom.get_root(h)
1941  if root == h:
1942  is_root=True
1943  is_molecule=IMP.atom.Molecule.get_is_setup(h)
1944  if is_molecule:
1945  molecules.add(IMP.atom.Molecule(h))
1946  h=h.get_parent()
1947  return list(molecules)
1948 
1949 def get_molecules_dictionary(input_objects):
1950  moldict=defaultdict(list)
1951  for mol in IMP.pmi.tools.get_molecules(input_objects):
1952  name=mol.get_name()
1953  moldict[name].append(mol)
1954 
1955  for mol in moldict:
1956  moldict[mol].sort(key=lambda x: IMP.atom.Copy(x).get_copy_index())
1957  return moldict
1958 
1959 def get_molecules_dictionary_by_copy(input_objects):
1960  moldict=defaultdict(dict)
1961  for mol in IMP.pmi.tools.get_molecules(input_objects):
1962  name=mol.get_name()
1963  c=IMP.atom.Copy(mol).get_copy_index()
1964  moldict[name][c]=mol
1965  return moldict
1966 
1967 def get_selections_dictionary(input_objects):
1968  moldict=IMP.pmi.tools.get_molecules_dictionary(input_objects)
1969  seldict=defaultdict(list)
1970  for name, mols in moldict.items():
1971  for m in mols:
1972  seldict[name].append(IMP.atom.Selection(m))
1973  return seldict
1974 
1975 def get_densities(input_objects):
1976  """Given a list of PMI objects, returns all density hierarchies within
1977  these objects. The output of this function can be inputted into
1978  things such as EM restraints. This function is intended to gather density particles
1979  appended to molecules (and not other hierarchies which might have been appended to the root node directly).
1980  """
1981  # Note that Densities can only be selected at the Root or Molecule level and not at the Leaves level.
1982  # we'll first get all molecule hierarchies corresponding to the leaves.
1983  molecules=get_molecules(input_objects)
1984  densities=[]
1985  for i in molecules:
1986  densities+=IMP.atom.Selection(i,representation_type=IMP.atom.DENSITIES).get_selected_particles()
1987  return densities
1988 
1989 def shuffle_configuration(objects,
1990  max_translation=300., max_rotation=2.0 * pi,
1991  avoidcollision_rb=True, avoidcollision_fb=False,
1992  cutoff=10.0, niterations=100,
1993  bounding_box=None,
1994  excluded_rigid_bodies=[],
1995  hierarchies_excluded_from_collision=[],
1996  hierarchies_included_in_collision=[],
1997  verbose=False,
1998  return_debug=False):
1999  """Shuffle particles. Used to restart the optimization.
2000  The configuration of the system is initialized by placing each
2001  rigid body and each bead randomly in a box. If `bounding_box` is
2002  specified, the particles are placed inside this box; otherwise, each
2003  particle is displaced by up to max_translation angstroms, and randomly
2004  rotated. Effort is made to place particles far enough from each other to
2005  prevent any steric clashes.
2006  @param objects Can be one of the following inputs:
2007  IMP Hierarchy, PMI System/State/Molecule/TempResidue, or a list/set of them
2008  @param max_translation Max translation (rbs and flexible beads)
2009  @param max_rotation Max rotation (rbs only)
2010  @param avoidcollision_rb check if the particle/rigid body was
2011  placed close to another particle; uses the optional
2012  arguments cutoff and niterations
2013  @param avoidcollision_fb Advanced. Generally you want this False because it's hard to shuffle beads.
2014  @param cutoff Distance less than this is a collision
2015  @param niterations How many times to try avoiding collision
2016  @param bounding_box Only shuffle particles within this box. Defined by ((x1,y1,z1),(x2,y2,z2)).
2017  @param excluded_rigid_bodies Don't shuffle these rigid body objects
2018  @param hierarchies_excluded_from_collision Don't count collision with these bodies
2019  @param hierarchies_included_in_collision Hierarchies that are not shuffled, but should be included in collision calculation (for fixed regions)
2020  @param verbose Give more output
2021  @note Best to only call this function after you've set up degrees of freedom
2022  For debugging purposes, returns: <shuffled indexes>, <collision avoided indexes>
2023  """
2024 
2025  ### checking input
2026  hierarchies = IMP.pmi.tools.input_adaptor(objects,
2027  pmi_resolution='all',
2028  flatten=True)
2029  rigid_bodies,flexible_beads = get_rbs_and_beads(hierarchies)
2030  if len(rigid_bodies)>0:
2031  mdl = rigid_bodies[0].get_model()
2032  elif len(flexible_beads)>0:
2033  mdl = flexible_beads[0].get_model()
2034  else:
2035  raise Exception("Could not find any particles in the hierarchy")
2036  if len(rigid_bodies) == 0:
2037  print("shuffle_configuration: rigid bodies were not intialized")
2038 
2039  ### gather all particles
2041  gcpf.set_distance(cutoff)
2042 
2043  # Add particles from excluded hierarchies to excluded list
2044  collision_excluded_hierarchies = IMP.pmi.tools.input_adaptor(hierarchies_excluded_from_collision,
2045  pmi_resolution='all',
2046  flatten=True)
2047 
2048  collision_included_hierarchies = IMP.pmi.tools.input_adaptor(hierarchies_included_in_collision,
2049  pmi_resolution='all',
2050  flatten=True)
2051 
2052  collision_excluded_idxs = set([l.get_particle().get_index() for h in collision_excluded_hierarchies \
2053  for l in IMP.core.get_leaves(h)])
2054 
2055  collision_included_idxs = set([l.get_particle().get_index() for h in collision_included_hierarchies \
2056  for l in IMP.core.get_leaves(h)])
2057 
2058  # Excluded collision with Gaussians
2059  all_idxs = [] #expand to representations?
2060  for p in IMP.pmi.tools.get_all_leaves(hierarchies):
2062  all_idxs.append(p.get_particle_index())
2064  collision_excluded_idxs.add(p.get_particle_index())
2065 
2066  if bounding_box is not None:
2067  ((x1, y1, z1), (x2, y2, z2)) = bounding_box
2068  ub = IMP.algebra.Vector3D(x1, y1, z1)
2069  lb = IMP.algebra.Vector3D(x2, y2, z2)
2070  bb = IMP.algebra.BoundingBox3D(ub, lb)
2071 
2072  all_idxs = set(all_idxs) | collision_included_idxs
2073  all_idxs = all_idxs - collision_excluded_idxs
2074  debug = []
2075  print('shuffling', len(rigid_bodies), 'rigid bodies')
2076  for rb in rigid_bodies:
2077  if rb not in excluded_rigid_bodies:
2078  # gather particles to avoid with this transform
2079  if avoidcollision_rb:
2080  rb_idxs = set(rb.get_member_particle_indexes()) - \
2081  collision_excluded_idxs
2082  other_idxs = all_idxs - rb_idxs
2083  if not other_idxs:
2084  continue
2085 
2086  # iterate, trying to avoid collisions
2087  niter = 0
2088  while niter < niterations:
2089  rbxyz = (rb.get_x(), rb.get_y(), rb.get_z())
2090 
2091  # local transform
2092  if bounding_box:
2093  translation = IMP.algebra.get_random_vector_in(bb)
2094  # First move to origin
2096  -IMP.core.XYZ(rb).get_coordinates())
2097  IMP.core.transform(rb, transformation_orig)
2098  old_coord = IMP.core.XYZ(rb).get_coordinates()
2100  transformation = IMP.algebra.Transformation3D(rotation,
2101  translation)
2102 
2103  else:
2105  rbxyz,
2106  max_translation,
2107  max_rotation)
2108 
2109  debug.append([rb, other_idxs if avoidcollision_rb else set()])
2110  IMP.core.transform(rb, transformation)
2111 
2112  # check collisions
2113  if avoidcollision_rb:
2114  mdl.update()
2115  npairs = len(gcpf.get_close_pairs(mdl,
2116  list(other_idxs),
2117  list(rb_idxs)))
2118  #print("NPAIRS:", npairs)
2119  if npairs==0:
2120  break
2121  else:
2122  niter += 1
2123  if verbose:
2124  print("shuffle_configuration: rigid body placed close to other %d particles, trying again..." % npairs)
2125  print("shuffle_configuration: rigid body name: " + rb.get_name())
2126  if niter == niterations:
2127  raise ValueError("tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
2128  else:
2129  break
2130 
2131  print('shuffling', len(flexible_beads), 'flexible beads')
2132  for fb in flexible_beads:
2133  # gather particles to avoid
2134  if avoidcollision_fb:
2135  fb_idxs = set(IMP.get_indexes([fb]))
2136  other_idxs = all_idxs - fb_idxs
2137  if not other_idxs:
2138  continue
2139 
2140  # iterate, trying to avoid collisions
2141  niter = 0
2142  while niter < niterations:
2143  if bounding_box:
2144  translation = IMP.algebra.get_random_vector_in(bb)
2145  transformation = IMP.algebra.Transformation3D(translation)
2146  else:
2147  fbxyz = IMP.core.XYZ(fb).get_coordinates()
2148  transformation = IMP.algebra.get_random_local_transformation(fbxyz,
2149  max_translation,
2150  max_rotation)
2151 
2152  # For gaussians, treat this fb as an rb
2154  memb = IMP.core.NonRigidMember(fb)
2155  xyz = memb.get_internal_coordinates()
2156  if bounding_box:
2157  # 'translation' is the new desired position in global
2158  # coordinates; we need to convert that to internal
2159  # coordinates first using the rigid body's ref frame
2160  rf = memb.get_rigid_body().get_reference_frame()
2161  glob_to_int = rf.get_transformation_from()
2162  memb.set_internal_coordinates(
2163  glob_to_int.get_transformed(translation))
2164  else:
2165  xyz_transformed=transformation.get_transformed(xyz)
2166  memb.set_internal_coordinates(xyz_transformed)
2167  debug.append([xyz,other_idxs if avoidcollision_fb else set()])
2168  else:
2169  d =IMP.core.XYZ(fb)
2170  if bounding_box:
2171  # Translate to origin first
2172  IMP.core.transform(d, -d.get_coordinates())
2173  d =IMP.core.XYZ(fb)
2174  debug.append([d,other_idxs if avoidcollision_fb else set()])
2175  if IMP.core.RigidBody.get_is_setup(fb.get_particle()):
2176  IMP.core.transform(IMP.core.RigidBody(fb.get_particle()), transformation)
2177  else:
2178  IMP.core.transform(d, transformation)
2179 
2180  if avoidcollision_fb:
2181  mdl.update()
2182  npairs = len(gcpf.get_close_pairs(mdl,
2183  list(other_idxs),
2184  list(fb_idxs)))
2185 
2186  if npairs==0:
2187  break
2188  else:
2189  niter += 1
2190  print("shuffle_configuration: floppy body placed close to other %d particles, trying again..." % npairs)
2191  if niter == niterations:
2192  raise ValueError("tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
2193  else:
2194  break
2195  if return_debug:
2196  return debug
2197 
2198 class ColorHierarchy(object):
2199 
2200  def __init__(self,hier):
2201  import matplotlib as mpl
2202  mpl.use('Agg')
2203  import matplotlib.pyplot as plt
2204  self.mpl = mpl
2205  self.plt = plt
2206 
2207  hier.ColorHierarchy=self
2208  self.hier=hier
2209  mols=IMP.pmi.tools.get_molecules(IMP.atom.get_leaves(self.hier))
2210  self.mols=[IMP.pmi.topology.PMIMoleculeHierarchy(mol) for mol in mols]
2211  self.method=self.nochange
2212  self.scheme=None
2213  self.first=None
2214  self.last=None
2215 
2216  def nochange(self):
2217  pass
2218 
2219  def get_color(self,fl):
2220  return IMP.display.Color(*self.scheme(fl)[0:3])
2221 
2222  def get_log_scale(self,fl):
2223  import math
2224  eps=1.0
2225  return math.log(fl+eps)
2226 
2227  def color_by_resid(self):
2228  self.method=self.color_by_resid
2229  self.scheme=self.mpl.cm.rainbow
2230  for mol in self.mols:
2231  self.first=1
2233  for p in IMP.atom.get_leaves(mol):
2235  ri=IMP.atom.Residue(p).get_index()
2236  c=self.get_color(float(ri)/self.last)
2237  IMP.display.Colored(p).set_color(c)
2240  avr=sum(ris)/len(ris)
2241  c=self.get_color(float(avr)/self.last)
2242  IMP.display.Colored(p).set_color(c)
2243 
2244  def color_by_uncertainty(self):
2245  self.method=self.color_by_uncertainty
2246  self.scheme=self.mpl.cm.jet
2247  ps=IMP.atom.get_leaves(self.hier)
2248  unc_dict={}
2249  for p in ps:
2251  u=IMP.pmi.Uncertainty(p).get_uncertainty()
2252  unc_dict[p]=u
2253  self.first=self.get_log_scale(1.0) #math.log(min(unc_dict.values())+eps)
2254  self.last=self.get_log_scale(100.0) #math.log(max(unc_dict.values())+eps)
2255  for p in unc_dict:
2256  value=self.get_log_scale(unc_dict[p])
2257  if value>=self.last: value=self.last
2258  if value<=self.first: value=self.first
2259  c=self.get_color((value-self.first)/(self.last-self.first))
2260  IMP.display.Colored(p).set_color(c)
2261 
2262  def get_color_bar(self,filename):
2263  import matplotlib as mpl
2264  mpl.use('Agg')
2265  import matplotlib.pyplot as plt
2266  import math
2267  plt.clf()
2268  fig = plt.figure(figsize=(8, 3))
2269  ax1 = fig.add_axes([0.05, 0.80, 0.9, 0.15])
2270 
2271  cmap = self.scheme
2272  norm = mpl.colors.Normalize(vmin=0.0, vmax=1.0)
2273 
2274  if self.method == self.color_by_uncertainty:
2275  angticks=[1.0,2.5,5.0,10.0,25.0,50.0,100.0]
2276  vvalues=[]
2277  marks=[]
2278  for at in angticks:
2279  vvalue=(self.get_log_scale(at)-self.first)/(self.last-self.first)
2280  if vvalue <= 1.0 and vvalue >= 0.0:
2281  vvalues.append(vvalue)
2282  marks.append(str(at))
2283  cb1 = mpl.colorbar.ColorbarBase(ax1, cmap=cmap,
2284  norm=norm,
2285  ticks=vvalues,
2286  orientation='horizontal')
2287  print(self.first,self.last,marks,vvalues)
2288  cb1.ax.set_xticklabels(marks)
2289  cb1.set_label('Angstorm')
2290  plt.savefig(filename, dpi=150, transparent=True)
2291  plt.show()
2292 
2293 
2294 
2295 
2296 def color2rgb(colorname):
2297  """Given a chimera color name, return RGB"""
2298  d = {'aquamarine': (0.4980392156862745, 1.0, 0.8313725490196079),
2299  'black': (0.0, 0.0, 0.0),
2300  'blue': (0.0, 0.0, 1.0),
2301  'brown': (0.6470588235294118, 0.16470588235294117, 0.16470588235294117),
2302  'chartreuse': (0.4980392156862745, 1.0, 0.0),
2303  'coral': (1.0, 0.4980392156862745, 0.3137254901960784),
2304  'cornflower blue': (0.39215686274509803, 0.5843137254901961, 0.9294117647058824),
2305  'cyan': (0.0, 1.0, 1.0),
2306  'dark cyan': (0.0, 0.5450980392156862, 0.5450980392156862),
2307  'dark gray': (0.6627450980392157, 0.6627450980392157, 0.6627450980392157),
2308  'dark green': (0.0, 0.39215686274509803, 0.0),
2309  'dark khaki': (0.7411764705882353, 0.7176470588235294, 0.4196078431372549),
2310  'dark magenta': (0.5450980392156862, 0.0, 0.5450980392156862),
2311  'dark olive green': (0.3333333333333333, 0.4196078431372549, 0.1843137254901961),
2312  'dark red': (0.5450980392156862, 0.0, 0.0),
2313  'dark slate blue': (0.2823529411764706, 0.23921568627450981, 0.5450980392156862),
2314  'dark slate gray': (0.1843137254901961, 0.30980392156862746, 0.30980392156862746),
2315  'deep pink': (1.0, 0.0784313725490196, 0.5764705882352941),
2316  'deep sky blue': (0.0, 0.7490196078431373, 1.0),
2317  'dim gray': (0.4117647058823529, 0.4117647058823529, 0.4117647058823529),
2318  'dodger blue': (0.11764705882352941, 0.5647058823529412, 1.0),
2319  'firebrick': (0.6980392156862745, 0.13333333333333333, 0.13333333333333333),
2320  'forest green': (0.13333333333333333, 0.5450980392156862, 0.13333333333333333),
2321  'gold': (1.0, 0.8431372549019608, 0.0),
2322  'goldenrod': (0.8549019607843137, 0.6470588235294118, 0.12549019607843137),
2323  'gray': (0.7450980392156863, 0.7450980392156863, 0.7450980392156863),
2324  'green': (0.0, 1.0, 0.0),
2325  'hot pink': (1.0, 0.4117647058823529, 0.7058823529411765),
2326  'khaki': (0.9411764705882353, 0.9019607843137255, 0.5490196078431373),
2327  'light blue': (0.6784313725490196, 0.8470588235294118, 0.9019607843137255),
2328  'light gray': (0.8274509803921568, 0.8274509803921568, 0.8274509803921568),
2329  'light green': (0.5647058823529412, 0.9333333333333333, 0.5647058823529412),
2330  'light sea green': (0.12549019607843137, 0.6980392156862745, 0.6666666666666666),
2331  'lime green': (0.19607843137254902, 0.803921568627451, 0.19607843137254902),
2332  'magenta': (1.0, 0.0, 1.0),
2333  'medium blue': (0.19607843137254902, 0.19607843137254902, 0.803921568627451),
2334  'medium purple': (0.5764705882352941, 0.4392156862745098, 0.8588235294117647),
2335  'navy blue': (0.0, 0.0, 0.5019607843137255),
2336  'olive drab': (0.4196078431372549, 0.5568627450980392, 0.13725490196078433),
2337  'orange red': (1.0, 0.27058823529411763, 0.0),
2338  'orange': (1.0, 0.4980392156862745, 0.0),
2339  'orchid': (0.8549019607843137, 0.4392156862745098, 0.8392156862745098),
2340  'pink': (1.0, 0.7529411764705882, 0.796078431372549),
2341  'plum': (0.8666666666666667, 0.6274509803921569, 0.8666666666666667),
2342  'purple': (0.6274509803921569, 0.12549019607843137, 0.9411764705882353),
2343  'red': (1.0, 0.0, 0.0),
2344  'rosy brown': (0.7372549019607844, 0.5607843137254902, 0.5607843137254902),
2345  'salmon': (0.9803921568627451, 0.5019607843137255, 0.4470588235294118),
2346  'sandy brown': (0.9568627450980393, 0.6431372549019608, 0.3764705882352941),
2347  'sea green': (0.1803921568627451, 0.5450980392156862, 0.3411764705882353),
2348  'sienna': (0.6274509803921569, 0.3215686274509804, 0.17647058823529413),
2349  'sky blue': (0.5294117647058824, 0.807843137254902, 0.9215686274509803),
2350  'slate gray': (0.4392156862745098, 0.5019607843137255, 0.5647058823529412),
2351  'spring green': (0.0, 1.0, 0.4980392156862745),
2352  'steel blue': (0.27450980392156865, 0.5098039215686274, 0.7058823529411765),
2353  'tan': (0.8235294117647058, 0.7058823529411765, 0.5490196078431373),
2354  'turquoise': (0.25098039215686274, 0.8784313725490196, 0.8156862745098039),
2355  'violet red': (0.8156862745098039, 0.12549019607843137, 0.5647058823529412),
2356  'white': (1.0, 1.0, 1.0),
2357  'yellow': (1.0, 1.0, 0.0)}
2358  return d[colorname]
2359 
2360 @IMP.deprecated_object(2.10,
2361  "If you use this class please let the PMI developers know.")
2362 class Colors(object):
2363  def __init__(self):
2364  self.colors={
2365  "reds":[("maroon","#800000",(128,0,0)),("dark red","#8B0000",(139,0,0)),
2366  ("brown","#A52A2A",(165,42,42)),("firebrick","#B22222",(178,34,34)),
2367  ("crimson","#DC143C",(220,20,60)),("red","#FF0000",(255,0,0)),
2368  ("tomato","#FF6347",(255,99,71)),("coral","#FF7F50",(255,127,80)),
2369  ("indian red","#CD5C5C",(205,92,92)),("light coral","#F08080",(240,128,128)),
2370  ("dark salmon","#E9967A",(233,150,122)),("salmon","#FA8072",(250,128,114)),
2371  ("light salmon","#FFA07A",(255,160,122)),("orange red","#FF4500",(255,69,0)),
2372  ("dark orange","#FF8C00",(255,140,0))],
2373  "yellows":[("orange","#FFA500",(255,165,0)),("gold","#FFD700",(255,215,0)),
2374  ("dark golden rod","#B8860B",(184,134,11)),("golden rod","#DAA520",(218,165,32)),
2375  ("pale golden rod","#EEE8AA",(238,232,170)),("dark khaki","#BDB76B",(189,183,107)),
2376  ("khaki","#F0E68C",(240,230,140)),("olive","#808000",(128,128,0)),
2377  ("yellow","#FFFF00",(255,255,0)),("antique white","#FAEBD7",(250,235,215)),
2378  ("beige","#F5F5DC",(245,245,220)),("bisque","#FFE4C4",(255,228,196)),
2379  ("blanched almond","#FFEBCD",(255,235,205)),("wheat","#F5DEB3",(245,222,179)),
2380  ("corn silk","#FFF8DC",(255,248,220)),("lemon chiffon","#FFFACD",(255,250,205)),
2381  ("light golden rod yellow","#FAFAD2",(250,250,210)),("light yellow","#FFFFE0",(255,255,224))],
2382  "greens":[("yellow green","#9ACD32",(154,205,50)),("dark olive green","#556B2F",(85,107,47)),
2383  ("olive drab","#6B8E23",(107,142,35)),("lawn green","#7CFC00",(124,252,0)),
2384  ("chart reuse","#7FFF00",(127,255,0)),("green yellow","#ADFF2F",(173,255,47)),
2385  ("dark green","#006400",(0,100,0)),("green","#008000",(0,128,0)),
2386  ("forest green","#228B22",(34,139,34)),("lime","#00FF00",(0,255,0)),
2387  ("lime green","#32CD32",(50,205,50)),("light green","#90EE90",(144,238,144)),
2388  ("pale green","#98FB98",(152,251,152)),("dark sea green","#8FBC8F",(143,188,143)),
2389  ("medium spring green","#00FA9A",(0,250,154)),("spring green","#00FF7F",(0,255,127)),
2390  ("sea green","#2E8B57",(46,139,87)),("medium aqua marine","#66CDAA",(102,205,170)),
2391  ("medium sea green","#3CB371",(60,179,113)),("light sea green","#20B2AA",(32,178,170)),
2392  ("dark slate gray","#2F4F4F",(47,79,79)),("teal","#008080",(0,128,128)),
2393  ("dark cyan","#008B8B",(0,139,139))],
2394  "blues":[("dark turquoise","#00CED1",(0,206,209)),
2395  ("turquoise","#40E0D0",(64,224,208)),("medium turquoise","#48D1CC",(72,209,204)),
2396  ("pale turquoise","#AFEEEE",(175,238,238)),("aqua marine","#7FFFD4",(127,255,212)),
2397  ("powder blue","#B0E0E6",(176,224,230)),("cadet blue","#5F9EA0",(95,158,160)),
2398  ("steel blue","#4682B4",(70,130,180)),("corn flower blue","#6495ED",(100,149,237)),
2399  ("deep sky blue","#00BFFF",(0,191,255)),("dodger blue","#1E90FF",(30,144,255)),
2400  ("light blue","#ADD8E6",(173,216,230)),("sky blue","#87CEEB",(135,206,235)),
2401  ("light sky blue","#87CEFA",(135,206,250)),("midnight blue","#191970",(25,25,112)),
2402  ("navy","#000080",(0,0,128)),("dark blue","#00008B",(0,0,139)),
2403  ("medium blue","#0000CD",(0,0,205)),("blue","#0000FF",(0,0,255)),("royal blue","#4169E1",(65,105,225)),
2404  ("aqua","#00FFFF",(0,255,255)),("cyan","#00FFFF",(0,255,255)),("light cyan","#E0FFFF",(224,255,255))],
2405  "violets":[("blue violet","#8A2BE2",(138,43,226)),("indigo","#4B0082",(75,0,130)),
2406  ("dark slate blue","#483D8B",(72,61,139)),("slate blue","#6A5ACD",(106,90,205)),
2407  ("medium slate blue","#7B68EE",(123,104,238)),("medium purple","#9370DB",(147,112,219)),
2408  ("dark magenta","#8B008B",(139,0,139)),("dark violet","#9400D3",(148,0,211)),
2409  ("dark orchid","#9932CC",(153,50,204)),("medium orchid","#BA55D3",(186,85,211)),
2410  ("purple","#800080",(128,0,128)),("thistle","#D8BFD8",(216,191,216)),
2411  ("plum","#DDA0DD",(221,160,221)),("violet","#EE82EE",(238,130,238)),
2412  ("magenta / fuchsia","#FF00FF",(255,0,255)),("orchid","#DA70D6",(218,112,214)),
2413  ("medium violet red","#C71585",(199,21,133)),("pale violet red","#DB7093",(219,112,147)),
2414  ("deep pink","#FF1493",(255,20,147)),("hot pink","#FF69B4",(255,105,180)),
2415  ("light pink","#FFB6C1",(255,182,193)),("pink","#FFC0CB",(255,192,203))],
2416  "browns":[("saddle brown","#8B4513",(139,69,19)),("sienna","#A0522D",(160,82,45)),
2417  ("chocolate","#D2691E",(210,105,30)),("peru","#CD853F",(205,133,63)),
2418  ("sandy brown","#F4A460",(244,164,96)),("burly wood","#DEB887",(222,184,135)),
2419  ("tan","#D2B48C",(210,180,140)),("rosy brown","#BC8F8F",(188,143,143)),
2420  ("moccasin","#FFE4B5",(255,228,181)),("navajo white","#FFDEAD",(255,222,173)),
2421  ("peach puff","#FFDAB9",(255,218,185)),("misty rose","#FFE4E1",(255,228,225)),
2422  ("lavender blush","#FFF0F5",(255,240,245)),("linen","#FAF0E6",(250,240,230)),
2423  ("old lace","#FDF5E6",(253,245,230)),("papaya whip","#FFEFD5",(255,239,213)),
2424  ("sea shell","#FFF5EE",(255,245,238))],
2425  "greys":[("black","#000000",(0,0,0)),("dim gray / dim grey","#696969",(105,105,105)),
2426  ("gray / grey","#808080",(128,128,128)),("dark gray / dark grey","#A9A9A9",(169,169,169)),
2427  ("silver","#C0C0C0",(192,192,192)),("light gray / light grey","#D3D3D3",(211,211,211)),
2428  ("gainsboro","#DCDCDC",(220,220,220)),("white smoke","#F5F5F5",(245,245,245)),
2429  ("white","#FFFFFF",(255,255,255))]}
2430 
2431  def assign_color_group(self,color_group,representation,component_names):
2432  for n,p in enumerate(component_names):
2433  s=IMP.atom.Selection(representation.prot,molecule=p)
2434  psel=s.get_selected_particles()
2435  ctuple=self.colors[color_group][n]
2436  print("Assigning "+p+" to color "+ctuple[0])
2437  c=ctuple[2]
2438  color=IMP.display.Color(float(c[0])/255,float(c[1])/255,float(c[2])/255)
2439  for part in psel:
2441  IMP.display.Colored(part).set_color(color)
2442  else:
2444 
2445  def get_list_distant_colors(self):
2446  cnames = ['#F0F8FF', '#FAEBD7', '#00FFFF', '#7FFFD4', '#F0FFFF', '#F5F5DC',
2447  '#FFE4C4', '#000000', '#FFEBCD', '#0000FF', '#8A2BE2', '#A52A2A', '#DEB887',
2448  '#5F9EA0', '#7FFF00', '#D2691E', '#FF7F50', '#6495ED', '#FFF8DC', '#DC143C',
2449  '#00FFFF', '#00008B', '#008B8B', '#B8860B', '#A9A9A9', '#006400', '#BDB76B',
2450  '#8B008B', '#556B2F', '#FF8C00', '#9932CC', '#8B0000', '#E9967A', '#8FBC8F',
2451  '#483D8B', '#2F4F4F', '#00CED1', '#9400D3', '#FF1493', '#00BFFF', '#696969',
2452  '#1E90FF', '#B22222', '#FFFAF0', '#228B22', '#FF00FF', '#DCDCDC', '#F8F8FF',
2453  '#FFD700', '#DAA520', '#808080', '#008000', '#ADFF2F', '#F0FFF0', '#FF69B4',
2454  '#CD5C5C', '#4B0082', '#FFFFF0', '#F0E68C', '#E6E6FA', '#FFF0F5', '#7CFC00',
2455  '#FFFACD', '#ADD8E6', '#F08080', '#E0FFFF', '#FAFAD2', '#90EE90', '#D3D3D3',
2456  '#FFB6C1', '#FFA07A', '#20B2AA', '#87CEFA', '#778899', '#B0C4DE', '#FFFFE0',
2457  '#00FF00', '#32CD32', '#FAF0E6', '#FF00FF', '#800000', '#66CDAA', '#0000CD',
2458  '#BA55D3', '#9370DB', '#3CB371', '#7B68EE', '#00FA9A', '#48D1CC', '#C71585',
2459  '#191970', '#F5FFFA', '#FFE4E1', '#FFE4B5', '#FFDEAD', '#000080', '#FDF5E6',
2460  '#808000', '#6B8E23', '#FFA500', '#FF4500', '#DA70D6', '#EEE8AA', '#98FB98',
2461  '#AFEEEE', '#DB7093', '#FFEFD5', '#FFDAB9', '#CD853F', '#FFC0CB', '#DDA0DD',
2462  '#B0E0E6', '#800080', '#FF0000', '#BC8F8F', '#4169E1', '#8B4513', '#FA8072',
2463  '#FAA460', '#2E8B57', '#FFF5EE', '#A0522D', '#C0C0C0', '#87CEEB', '#6A5ACD',
2464  '#708090', '#FFFAFA', '#00FF7F', '#4682B4', '#D2B48C', '#008080', '#D8BFD8',
2465  '#FF6347', '#40E0D0', '#EE82EE', '#F5DEB3', '#FFFFFF', '#F5F5F5', '#FFFF00',
2466  '#9ACD32']
2467  return cnames
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: Molecule.h:35
Simple 3D transformation class.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: Residue.h:155
Represent an RGB color.
Definition: Color.h:24
A decorator to associate a particle with a part of a protein/DNA/RNA.
Definition: Fragment.h:20
Extends the functionality of IMP.atom.Molecule.
def add_script_provenance
Tag the given particle with the current Python script.
def get_restraint_set
Get a RestraintSet containing all PMI restraints added to the model.
Definition: tools.py:88
Store the representations for a system.
Definition: tools.py:895
def get_random_cross_link_dataset
Return a random cross-link dataset as a string.
Definition: tools.py:394
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: atom/Atom.h:241
void add_particles(RMF::FileHandle fh, const ParticlesTemp &hs)
Set of Python classes to create a multi-state, multi-resolution IMP hierarchy.
def __init__
Constructor.
Definition: tools.py:106
static Weight setup_particle(Model *m, ParticleIndex pi)
Definition: Weight.h:31
A decorator for a particle which has bonds.
double get_drmsd(const Vector3DsOrXYZs0 &m0, const Vector3DsOrXYZs1 &m1)
Calculate distance the root mean square deviation between two sets of 3D points.
Definition: atom/distance.h:49
Rotation3D get_random_rotation_3d(const Rotation3D &center, double distance)
Pick a rotation at random near the provided one.
def deprecated_function
Python decorator to mark a function as deprecated.
Definition: __init__.py:9715
std::string get_module_version()
IMP::Vector< Color > Colors
Definition: Color.h:188
ParticlesTemp get_particles(Model *m, const ParticleIndexes &ps)
The type of an atom.
static Surface setup_particle(Model *m, ParticleIndex pi)
Definition: Surface.h:45
Add uncertainty to a particle.
Definition: Uncertainty.h:24
void add_particle(RMF::FileHandle fh, Particle *hs)
log
Definition: log.py:1
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:690
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: XYZ.h:49
GenericHierarchies get_leaves(Hierarchy mhd)
Get all the leaves of the bit of hierarchy.
def get_prot_name_from_particle
Return the component name provided a particle and a list of names.
Definition: tools.py:1043
This class initializes the root node of the global IMP.atom.Hierarchy.
double get_distance(XYZR a, XYZR b)
Compute the sphere distance between a and b.
Definition: XYZR.h:89
def get_residue_gaps_in_hierarchy
Return the residue index gaps and contiguous segments in the hierarchy.
Definition: tools.py:655
def add_imp_provenance
Tag the given particle as being created by the current version of IMP.
Vector3D get_random_vector_in(const Cylinder3D &c)
Generate a random vector in a cylinder with uniform density.
Add resolution to a particle.
Definition: Resolution.h:24
Bond create_bond(Bonded a, Bonded b, Bond o)
Connect the two wrapped particles by a custom bond.
Object used to hold a set of restraints.
Definition: RestraintSet.h:36
algebra::GridD< 3, algebra::DenseGridStorageD< 3, float >, float > get_grid(DensityMap *in_map)
Return a dense grid containing the voxels of the passed density map.
Stores a named protein chain.
def get_terminal_residue
Get the particle of the terminal residue at the GIVEN resolution (NOTE: not the closest resolution!)...
Definition: tools.py:602
static bool get_is_setup(Model *m, ParticleIndex pi)
Definition: Fragment.h:46
def add_software_provenance
Tag the given particle with the software used to create it.
A decorator for keeping track of copies of a molecule.
Definition: Copy.h:28
ParticleIndexPairs get_indexes(const ParticlePairsTemp &ps)
static bool get_is_setup(Model *m, ParticleIndex pi)
Definition: Colored.h:48
The standard decorator for manipulating molecular structures.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
def select_by_tuple_2
New tuple format: molname OR (start,stop,molname,copynum,statenum) Copy and state are optional...
Definition: tools.py:845
A decorator for a particle representing an atom.
Definition: atom/Atom.h:234
def scatter_and_gather
Synchronize data over a parallel run.
Definition: tools.py:1115
void transform(XYZ a, const algebra::Transformation3D &tr)
Apply a transformation to the particle.
def add_restraint_to_model
Add a PMI restraint to the model.
Definition: tools.py:69
static Bonded setup_particle(Model *m, ParticleIndex pi)
The type for a residue.
def get_terminal_residue_position
Get XYZ coordinates of the terminal residue at the GIVEN resolution.
Definition: tools.py:643
void load_frame(RMF::FileConstHandle file, RMF::FrameID frame)
Load the given RMF frame into the state of the linked objects.
A decorator for a particle with x,y,z coordinates.
Definition: XYZ.h:30
static Scale setup_particle(Model *m, ParticleIndex pi)
Definition: Scale.h:30
A base class for Keys.
Definition: Key.h:44
Collect timing information.
Definition: tools.py:101
static Colored setup_particle(Model *m, ParticleIndex pi, Color color)
Definition: Colored.h:62
A decorator for a particle that is part of a rigid body but not rigid.
Definition: rigid_bodies.h:709
def get_closest_residue_position
this function works with plain hierarchies, as read from the pdb, no multi-scale hierarchies ...
Definition: tools.py:571
std::ostream & show(Hierarchy h, std::ostream &out=std::cout)
Print the hierarchy using a given decorator to display each node.
Find all nearby pairs by testing all pairs.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:158
def deprecated_object
Python decorator to mark a class as deprecated.
Definition: __init__.py:9687
A decorator for a residue.
Definition: Residue.h:134
General purpose algebraic and geometric methods that are expected to be used by a wide variety of IMP...
static bool get_is_setup(const IMP::ParticleAdaptor &p)
static bool get_is_setup(Model *m, ParticleIndex p)
Check if the particle has the needed attributes for a cast to succeed.
The general base class for IMP exceptions.
Definition: exception.h:49
VectorD< 3 > Vector3D
Definition: VectorD.h:395
Rotation3D get_identity_rotation_3d()
Return a rotation that does not do anything.
Definition: Rotation3D.h:302
void link_hierarchies(RMF::FileConstHandle fh, const atom::Hierarchies &hs)
Class to handle individual particles of a Model object.
Definition: Particle.h:41
Bond get_bond(Bonded a, Bonded b)
Get the bond between two particles.
Stores a list of Molecules all with the same State index.
std::string get_data_path(std::string file_name)
Return the full path to one of this module's data files.
int get_copy_index(Hierarchy h)
Walk up the hierarchy to find the current copy index.
def select
this function uses representation=SimplifiedModel it returns the corresponding selected particles rep...
Definition: tools.py:726
Python classes to represent, score, sample and analyze models.
static bool get_is_setup(Model *m, ParticleIndex pi)
Definition: Uncertainty.h:30
double get_resolution(Model *m, ParticleIndex pi)
Estimate the resolution of the hierarchy as used by Representation.
A decorator for a rigid body.
Definition: rigid_bodies.h:82
def cross_link_db_filter_parser
example '"{ID_Score}" > 28 AND "{Sample}" == "%10_1%" OR ":Sample}" == "%10_2%" OR ":Sample...
Definition: tools.py:501
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: core/Gaussian.h:47
Hierarchies get_leaves(const Selection &h)
A decorator for a molecule.
Definition: Molecule.h:24
Select hierarchy particles identified by the biological name.
Definition: Selection.h:66
Support for the RMF file format for storing hierarchical molecular data and markup.
def get_residue_indexes
Retrieve the residue indexes for the given particle.
Definition: tools.py:1063
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:711
Transformation3D get_random_local_transformation(Vector3D origin, double max_translation=5., double max_angle_in_rad=0.26)
Get a local transformation.
Temporarily stores residue information, even without structure available.
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...
A particle with a color.
Definition: Colored.h:23