IMP logo
IMP Reference Guide  2.11.1
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 import warnings
28 
29 def _get_system_for_hier(hier):
30  """Given a hierarchy, return the System that created it, or None"""
31  # If we are given the raw particle, get the corresponding Hierarchy
32  # decorator if available
33  if hier and not hasattr(hier, 'get_parent'):
35  hier = IMP.atom.Hierarchy(hier)
36  else:
37  return None
38  while hier:
39  # See if we labeled the Python object directly with the System
40  if hasattr(hier, '_pmi2_system'):
41  return hier._pmi2_system()
42  # Otherwise (maybe we got a new Python wrapper around the same C++
43  # object), try all extant systems
44  for ws in IMP.pmi.topology.System._all_systems:
45  s = ws()
46  if s and s.hier == hier:
47  return s
48  # Try the next level up in the hierarchy
49  hier = hier.get_parent()
50 
51 def _all_protocol_outputs(representations, hier):
52  """Iterate over all (ProtocolOutput, State) pairs for the given
53  representations (PMI1) or hier (PMI2)"""
54  if hier:
55  system = _get_system_for_hier(hier)
56  if system:
57  for state in system.states:
58  for p in state._protocol_output:
59  yield p
60  else:
61  for p in representations[0]._protocol_output:
62  yield p
63 
64 def _add_pmi_provenance(p):
65  """Tag the given particle as being created by the current version of PMI."""
67  IMP.core.add_software_provenance(p, name="IMP PMI module",
69  location="https://integrativemodeling.org")
71 
72 def _get_restraint_set_keys():
73  if not hasattr(_get_restraint_set_keys, 'pmi_rs_key'):
74  _get_restraint_set_keys.pmi_rs_key = IMP.ModelKey("PMI restraints")
75  _get_restraint_set_keys.rmf_rs_key = IMP.ModelKey("RMF restraints")
76  return (_get_restraint_set_keys.pmi_rs_key,
77  _get_restraint_set_keys.rmf_rs_key)
78 
79 def _add_restraint_sets(model, mk, mk_rmf):
80  rs = IMP.RestraintSet(model, "All PMI restraints")
81  rs_rmf = IMP.RestraintSet(model, "All PMI RMF restraints")
82  model.add_data(mk, rs)
83  model.add_data(mk_rmf, rs_rmf)
84  return rs, rs_rmf
85 
86 def add_restraint_to_model(model, restraint, add_to_rmf=False):
87  """Add a PMI restraint to the model.
88  Since Model.add_restraint() no longer exists (in modern IMP restraints
89  should be added to a ScoringFunction instead) store them instead in
90  a RestraintSet, and keep a reference to it in the Model.
91 
92  If `add_to_rmf` is True, also add the restraint to a separate list
93  of restraints that will be written out to RMF files (by default, most
94  PMI restraints are not)."""
95  mk, mk_rmf = _get_restraint_set_keys()
96  if model.get_has_data(mk):
97  rs = IMP.RestraintSet.get_from(model.get_data(mk))
98  rs_rmf = IMP.RestraintSet.get_from(model.get_data(mk_rmf))
99  else:
100  rs, rs_rmf = _add_restraint_sets(model, mk, mk_rmf)
101  rs.add_restraint(restraint)
102  if add_to_rmf:
103  rs_rmf.add_restraint(restraint)
104 
105 def get_restraint_set(model, rmf=False):
106  """Get a RestraintSet containing all PMI restraints added to the model.
107  If `rmf` is True, return only the subset of these restraints that
108  should be written out to RMF files."""
109  mk, mk_rmf = _get_restraint_set_keys()
110  if not model.get_has_data(mk):
111  warnings.warn("no restraints added to model yet",
113  _add_restraint_sets(model, mk, mk_rmf)
114  if rmf:
115  return IMP.RestraintSet.get_from(model.get_data(mk_rmf))
116  else:
117  return IMP.RestraintSet.get_from(model.get_data(mk))
118 
119 class Stopwatch(object):
120  """Collect timing information.
121  Add an instance of this class to outputobjects to get timing information
122  in a stat file."""
123 
124  def __init__(self, isdelta=True):
125  """Constructor.
126  @param isdelta if True (the default) then report the time since the
127  last use of this class; if False, report cumulative time."""
128  self.starttime = time.clock()
129  self.label = "None"
130  self.isdelta = isdelta
131 
132  def set_label(self, labelstr):
133  self.label = labelstr
134 
135  def get_output(self):
136  output = {}
137  if self.isdelta:
138  newtime = time.clock()
139  output["Stopwatch_" + self.label + "_delta_seconds"] \
140  = str(newtime - self.starttime)
141  self.starttime = newtime
142  else:
143  output["Stopwatch_" + self.label + "_elapsed_seconds"] \
144  = str(time.clock() - self.starttime)
145  return output
146 
147 
148 class SetupNuisance(object):
149 
150  def __init__(self, m, initialvalue, minvalue, maxvalue, isoptimized=True):
151 
152  nuisance = IMP.isd.Scale.setup_particle(IMP.Particle(m), initialvalue)
153  if minvalue:
154  nuisance.set_lower(minvalue)
155  if maxvalue:
156  nuisance.set_upper(maxvalue)
157 
158  # m.add_score_state(IMP.core.SingletonConstraint(IMP.isd.NuisanceRangeModifier(),None,nuisance))
159  nuisance.set_is_optimized(nuisance.get_nuisance_key(), isoptimized)
160  self.nuisance = nuisance
161 
162  def get_particle(self):
163  return self.nuisance
164 
165 
166 class SetupWeight(object):
167 
168  def __init__(self, m, isoptimized=True):
169  pw = IMP.Particle(m)
170  self.weight = IMP.isd.Weight.setup_particle(pw)
171  self.weight.set_weights_are_optimized(True)
172 
173  def get_particle(self):
174  return self.weight
175 
176 
177 class SetupSurface(object):
178 
179  def __init__(self, m, center, normal, isoptimized=True):
180  p = IMP.Particle(m)
181  self.surface = IMP.core.Surface.setup_particle(p, center, normal)
182  self.surface.set_coordinates_are_optimized(isoptimized)
183  self.surface.set_normal_is_optimized(isoptimized)
184 
185  def get_particle(self):
186  return self.surface
187 
188 
189 class ParticleToSampleList(object):
190 
191  def __init__(self, label="None"):
192 
193  self.dictionary_particle_type = {}
194  self.dictionary_particle_transformation = {}
195  self.dictionary_particle_name = {}
196  self.label = label
197 
198  def add_particle(
199  self,
200  particle,
201  particle_type,
202  particle_transformation,
203  name):
204  if not particle_type in ["Rigid_Bodies", "Floppy_Bodies", "Nuisances", "X_coord", "Weights", "Surfaces"]:
205  raise TypeError("not the right particle type")
206  else:
207  self.dictionary_particle_type[particle] = particle_type
208  if particle_type == "Rigid_Bodies":
209  if (isinstance(particle_transformation, tuple)
210  and len(particle_transformation) == 2
211  and all(isinstance(x, float)
212  for x in particle_transformation)):
213  self.dictionary_particle_transformation[
214  particle] = particle_transformation
215  self.dictionary_particle_name[particle] = name
216  else:
217  raise TypeError("ParticleToSampleList: not the right transformation format for Rigid_Bodies, should be a tuple of floats")
218  elif particle_type == "Surfaces":
219  if (isinstance(particle_transformation, tuple)
220  and len(particle_transformation) == 3
221  and all(isinstance(x, float)
222  for x in particle_transformation)):
223  self.dictionary_particle_transformation[
224  particle] = particle_transformation
225  self.dictionary_particle_name[particle] = name
226  else:
227  raise TypeError("ParticleToSampleList: not the right transformation format for Surfaces, should be a tuple of floats")
228  else:
229  if isinstance(particle_transformation, float):
230  self.dictionary_particle_transformation[
231  particle] = particle_transformation
232  self.dictionary_particle_name[particle] = name
233  else:
234  raise TypeError("ParticleToSampleList: not the right transformation format, should be a float")
235 
236  def get_particles_to_sample(self):
237  ps = {}
238  for particle in self.dictionary_particle_type:
239  key = self.dictionary_particle_type[
240  particle] + "ParticleToSampleList_" + self.dictionary_particle_name[particle] + "_" + self.label
241  value = (
242  [particle],
243  self.dictionary_particle_transformation[particle])
244  ps[key] = value
245  return ps
246 
247 
248 def get_random_cross_link_dataset(representation,
249  resolution=1.0,
250  number_of_cross_links=10,
251  ambiguity_probability=0.1,
252  confidence_score_range=[0,100],
253  avoid_same_particles=False):
254  '''Return a random cross-link dataset as a string.
255  Every line is a residue pair, together with UniqueIdentifier
256  and XL score.'''
257 
258  residue_pairs=get_random_residue_pairs(representation, resolution, number_of_cross_links, avoid_same_particles=avoid_same_particles)
259 
260  unique_identifier=0
261  cmin=float(min(confidence_score_range))
262  cmax=float(max(confidence_score_range))
263 
264  dataset="#\n"
265 
266  for (name1, r1, name2, r2) in residue_pairs:
267  if random.random() > ambiguity_probability:
268  unique_identifier+=1
269  score=random.random()*(cmax-cmin)+cmin
270  dataset+=str(name1)+" "+str(name2)+" "+str(r1)+" "+str(r2)+" "+str(score)+" "+str(unique_identifier)+"\n"
271 
272  return dataset
273 
274 
275  #-------------------------------
276 
277 def get_cross_link_data(directory, filename, dist, omega, sigma,
278  don=None, doff=None, prior=0, type_of_profile="gofr"):
279 
280  (distmin, distmax, ndist) = dist
281  (omegamin, omegamax, nomega) = omega
282  (sigmamin, sigmamax, nsigma) = sigma
283 
284  filen = IMP.isd.get_data_path("CrossLinkPMFs.dict")
285  with open(filen) as xlpot:
286  dictionary = ast.literal_eval(xlpot.readline())
287 
288  xpot = dictionary[directory][filename]["distance"]
289  pot = dictionary[directory][filename][type_of_profile]
290 
291  dist_grid = get_grid(distmin, distmax, ndist, False)
292  omega_grid = get_log_grid(omegamin, omegamax, nomega)
293  sigma_grid = get_log_grid(sigmamin, sigmamax, nsigma)
294 
295  if not don is None and not doff is None:
296  xlmsdata = IMP.isd.CrossLinkData(
297  dist_grid,
298  omega_grid,
299  sigma_grid,
300  xpot,
301  pot,
302  don,
303  doff,
304  prior)
305  else:
306  xlmsdata = IMP.isd.CrossLinkData(
307  dist_grid,
308  omega_grid,
309  sigma_grid,
310  xpot,
311  pot)
312  return xlmsdata
313 
314 def get_grid(gmin, gmax, ngrid, boundaries):
315  grid = []
316  dx = (gmax - gmin) / float(ngrid)
317  for i in range(0, ngrid + 1):
318  if(not boundaries and i == 0):
319  continue
320  if(not boundaries and i == ngrid):
321  continue
322  grid.append(gmin + float(i) * dx)
323  return grid
324 
325  #-------------------------------
326 
327 
328 def get_log_grid(gmin, gmax, ngrid):
329  grid = []
330  for i in range(0, ngrid + 1):
331  grid.append(gmin * exp(float(i) / ngrid * log(gmax / gmin)))
332  return grid
333 
334  #-------------------------------
335 
336 
338  '''
339  example '"{ID_Score}" > 28 AND "{Sample}" ==
340  "%10_1%" OR ":Sample}" == "%10_2%" OR ":Sample}"
341  == "%10_3%" OR ":Sample}" == "%8_1%" OR ":Sample}" == "%8_2%"'
342  '''
343 
344  import pyparsing as pp
345 
346  operator = pp.Regex(">=|<=|!=|>|<|==|in").setName("operator")
347  value = pp.QuotedString(
348  '"') | pp.Regex(
349  r"[+-]?\d+(:?\.\d*)?(:?[eE][+-]?\d+)?")
350  identifier = pp.Word(pp.alphas, pp.alphanums + "_")
351  comparison_term = identifier | value
352  condition = pp.Group(comparison_term + operator + comparison_term)
353 
354  expr = pp.operatorPrecedence(condition, [
355  ("OR", 2, pp.opAssoc.LEFT, ),
356  ("AND", 2, pp.opAssoc.LEFT, ),
357  ])
358 
359  parsedstring = str(expr.parseString(inputstring)) \
360  .replace("[", "(") \
361  .replace("]", ")") \
362  .replace(",", " ") \
363  .replace("'", " ") \
364  .replace("%", "'") \
365  .replace("{", "float(entry['") \
366  .replace("}", "'])") \
367  .replace(":", "str(entry['") \
368  .replace("}", "'])") \
369  .replace("AND", "and") \
370  .replace("OR", "or")
371  return parsedstring
372 
373 
374 def open_file_or_inline_text(filename):
375  try:
376  fl = open(filename, "r")
377  except IOError:
378  fl = filename.split("\n")
379  return fl
380 
381 def get_ids_from_fasta_file(fastafile):
382  ids = []
383  with open(fastafile) as ff:
384  for l in ff:
385  if l[0] == ">":
386  ids.append(l[1:-1])
387  return ids
388 
389 
390 def get_closest_residue_position(hier, resindex, terminus="N"):
391  '''
392  this function works with plain hierarchies, as read from the pdb,
393  no multi-scale hierarchies
394  '''
395  p = []
396  niter = 0
397  while len(p) == 0:
398  niter += 1
399  sel = IMP.atom.Selection(hier, residue_index=resindex,
400  atom_type=IMP.atom.AT_CA)
401 
402  if terminus == "N":
403  resindex += 1
404  if terminus == "C":
405  resindex -= 1
406 
407  if niter >= 10000:
408  print("get_closest_residue_position: exiting while loop without result")
409  break
410  p = sel.get_selected_particles()
411 
412  if len(p) == 1:
413  return IMP.core.XYZ(p[0]).get_coordinates()
414  elif len(p) == 0:
415  print("get_closest_residue_position: got NO residues for hierarchy %s and residue %i" % (hier, resindex))
416  raise Exception("get_closest_residue_position: got NO residues for hierarchy %s and residue %i" % (
417  hier, resindex))
418  else:
419  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])))
420 
421 def get_terminal_residue(representation, hier, terminus="C", resolution=1):
422  '''
423  Get the particle of the terminal residue at the GIVEN resolution
424  (NOTE: not the closest resolution!).
425  To get the terminal residue at the closest resolution use:
426  particles=IMP.pmi.tools.select_by_tuple(representation,molecule_name)
427  particles[0] and particles[-1] will be the first and last particles
428  corresponding to the two termini.
429  It is needed for instance to determine the last residue of a pdb.
430  @param hier hierarchy containing the terminal residue
431  @param terminus either 'N' or 'C'
432  @param resolution resolution to use.
433  '''
434  termresidue = None
435  termparticle = None
436 
437  ps=select(representation,
438  resolution=resolution,
439  hierarchies=[hier])
440 
441  for p in ps:
442  if IMP.pmi.Resolution(p).get_resolution() == resolution:
444  if terminus == "C":
445  if termresidue is None:
446  termresidue = max(residues)
447  termparticle = p
448  elif max(residues) >= termresidue:
449  termresidue = max(residues)
450  termparticle = p
451  elif terminus == "N":
452  if termresidue is None:
453  termresidue = min(residues)
454  termparticle = p
455  elif min(residues) <= termresidue:
456  termresidue = min(residues)
457  termparticle = p
458  else:
459  raise ValueError("terminus argument should be either N or C")
460  return termparticle
461 
462 def get_terminal_residue_position(representation, hier, terminus="C",
463  resolution=1):
464  """Get XYZ coordinates of the terminal residue at the GIVEN resolution"""
465  p = get_terminal_residue(representation, hier, terminus, resolution)
466  return IMP.core.XYZ(p).get_coordinates()
467 
468 def get_residue_gaps_in_hierarchy(hierarchy, start, end):
469  '''
470  Return the residue index gaps and contiguous segments in the hierarchy.
471 
472  @param hierarchy hierarchy to examine
473  @param start first residue index
474  @param end last residue index
475 
476  @return A list of lists of the form
477  [[1,100,"cont"],[101,120,"gap"],[121,200,"cont"]]
478  '''
479  gaps = []
480  for n, rindex in enumerate(range(start, end + 1)):
481  sel = IMP.atom.Selection(hierarchy, residue_index=rindex,
482  atom_type=IMP.atom.AT_CA)
483 
484  if len(sel.get_selected_particles()) == 0:
485  if n == 0:
486  # set the initial condition
487  rindexgap = start
488  rindexcont = start - 1
489  if rindexgap == rindex - 1:
490  # residue is contiguous with the previously discovered gap
491  gaps[-1][1] += 1
492  else:
493  # residue is not contiguous with the previously discovered gap
494  # hence create a new gap tuple
495  gaps.append([rindex, rindex, "gap"])
496  # update the index of the last residue gap
497  rindexgap = rindex
498  else:
499  if n == 0:
500  # set the initial condition
501  rindexgap = start - 1
502  rindexcont = start
503  if rindexcont == rindex - 1:
504  # residue is contiguous with the previously discovered
505  # continuous part
506  gaps[-1][1] += 1
507  else:
508  # residue is not contiguous with the previously discovered continuous part
509  # hence create a new cont tuple
510  gaps.append([rindex, rindex, "cont"])
511  # update the index of the last residue gap
512  rindexcont = rindex
513  return gaps
514 
515 
516 class map(object):
517 
518  def __init__(self):
519  self.map = {}
520 
521  def set_map_element(self, xvalue, yvalue):
522  self.map[xvalue] = yvalue
523 
524  def get_map_element(self, invalue):
525  if isinstance(invalue, float):
526  n = 0
527  mindist = 1
528  for x in self.map:
529  dist = (invalue - x) * (invalue - x)
530 
531  if n == 0:
532  mindist = dist
533  minx = x
534  if dist < mindist:
535  mindist = dist
536  minx = x
537  n += 1
538  return self.map[minx]
539  elif isinstance(invalue, str):
540  return self.map[invalue]
541  else:
542  raise TypeError("wrong type for map")
543 
544 
545 def select(representation,
546  resolution=None,
547  hierarchies=None,
548  selection_arguments=None,
549  name=None,
550  name_is_ambiguous=False,
551  first_residue=None,
552  last_residue=None,
553  residue=None,
554  representation_type=None):
555  '''
556  this function uses representation=SimplifiedModel
557  it returns the corresponding selected particles
558  representation_type="Beads", "Res:X", "Densities", "Representation", "Molecule"
559  '''
560 
561  if resolution is None:
562  allparticles = IMP.atom.get_leaves(representation.prot)
563  resolution_particles = None
564  hierarchies_particles = None
565  names_particles = None
566  residue_range_particles = None
567  residue_particles = None
568  representation_type_particles = None
569 
570  if not resolution is None:
571  resolution_particles = []
572  hs = representation.get_hierarchies_at_given_resolution(resolution)
573  for h in hs:
574  resolution_particles += IMP.atom.get_leaves(h)
575 
576  if not hierarchies is None:
577  hierarchies_particles = []
578  for h in hierarchies:
579  hierarchies_particles += IMP.atom.get_leaves(h)
580 
581  if not name is None:
582  names_particles = []
583  if name_is_ambiguous:
584  for namekey in representation.hier_dict:
585  if name in namekey:
586  names_particles += IMP.atom.get_leaves(
587  representation.hier_dict[namekey])
588  elif name in representation.hier_dict:
589  names_particles += IMP.atom.get_leaves(representation.hier_dict[name])
590  else:
591  print("select: component %s is not there" % name)
592 
593  if not first_residue is None and not last_residue is None:
594  sel = IMP.atom.Selection(representation.prot,
595  residue_indexes=range(first_residue, last_residue + 1))
596  residue_range_particles = [IMP.atom.Hierarchy(p)
597  for p in sel.get_selected_particles()]
598 
599  if not residue is None:
600  sel = IMP.atom.Selection(representation.prot, residue_index=residue)
601  residue_particles = [IMP.atom.Hierarchy(p)
602  for p in sel.get_selected_particles()]
603 
604  if not representation_type is None:
605  representation_type_particles = []
606  if representation_type == "Molecule":
607  for name in representation.hier_representation:
608  for repr_type in representation.hier_representation[name]:
609  if repr_type == "Beads" or "Res:" in repr_type:
610  h = representation.hier_representation[name][repr_type]
611  representation_type_particles += IMP.atom.get_leaves(h)
612 
613  elif representation_type == "PDB":
614  for name in representation.hier_representation:
615  for repr_type in representation.hier_representation[name]:
616  if repr_type == "Res:" in repr_type:
617  h = representation.hier_representation[name][repr_type]
618  representation_type_particles += IMP.atom.get_leaves(h)
619 
620  else:
621  for name in representation.hier_representation:
622  h = representation.hier_representation[
623  name][
624  representation_type]
625  representation_type_particles += IMP.atom.get_leaves(h)
626 
627  selections = [hierarchies_particles, names_particles,
628  residue_range_particles, residue_particles, representation_type_particles]
629 
630  if resolution is None:
631  selected_particles = set(allparticles)
632  else:
633  selected_particles = set(resolution_particles)
634 
635  for s in selections:
636  if not s is None:
637  selected_particles = (set(s) & selected_particles)
638 
639  return list(selected_particles)
640 
641 
642 def select_by_tuple(
643  representation,
644  tupleselection,
645  resolution=None,
646  name_is_ambiguous=False):
647  if isinstance(tupleselection, tuple) and len(tupleselection) == 3:
648  particles = IMP.pmi.tools.select(representation, resolution=resolution,
649  name=tupleselection[2],
650  first_residue=tupleselection[0],
651  last_residue=tupleselection[1],
652  name_is_ambiguous=name_is_ambiguous)
653  elif isinstance(tupleselection, str):
654  particles = IMP.pmi.tools.select(representation, resolution=resolution,
655  name=tupleselection,
656  name_is_ambiguous=name_is_ambiguous)
657  else:
658  raise ValueError('you passed something bad to select_by_tuple()')
659  # now order the result by residue number
660  particles = IMP.pmi.tools.sort_by_residues(particles)
661 
662  return particles
663 
664 def select_by_tuple_2(hier,tuple_selection,resolution):
665  """New tuple format: molname OR (start,stop,molname,copynum,statenum)
666  Copy and state are optional. Can also use 'None' for them which will get all.
667  You can also pass -1 for stop which will go to the end.
668  Returns the particles
669  """
670  kwds = {} # going to accumulate keywords
671  kwds['resolution'] = resolution
672  if isinstance(tuple_selection, str):
673  kwds['molecule'] = tuple_selection
674  elif isinstance(tuple_selection, tuple):
675  rbegin = tuple_selection[0]
676  rend = tuple_selection[1]
677  kwds['molecule'] = tuple_selection[2]
678  try:
679  copynum = tuple_selection[3]
680  if copynum is not None:
681  kwds['copy_index'] = copynum
682  except:
683  pass
684  try:
685  statenum = tuple_selection[4]
686  if statenum is not None:
687  kwds['state_index'] = statenum
688  except:
689  pass
690  if rend==-1:
691  if rbegin>1:
692  s = IMP.atom.Selection(hier,**kwds)
693  s -= IMP.atom.Selection(hier,
694  residue_indexes=range(1,rbegin),
695  **kwds)
696  return s.get_selected_particles()
697  else:
698  kwds['residue_indexes'] = range(rbegin,rend+1)
699  s = IMP.atom.Selection(hier,**kwds)
700  return s.get_selected_particles()
701 
702 
703 
704 def get_db_from_csv(csvfilename):
705  import csv
706  outputlist = []
707  with open(csvfilename) as fh:
708  csvr = csv.DictReader(fh)
709  for l in csvr:
710  outputlist.append(l)
711  return outputlist
712 
713 
714 class HierarchyDatabase(object):
715  """Store the representations for a system."""
716 
717  def __init__(self):
718  self.db = {}
719  # this dictionary map a particle to its root hierarchy
720  self.root_hierarchy_dict = {}
721  self.preroot_fragment_hierarchy_dict = {}
722  self.particle_to_name = {}
723  self.model = None
724 
725  def add_name(self, name):
726  if name not in self.db:
727  self.db[name] = {}
728 
729  def add_residue_number(self, name, resn):
730  resn = int(resn)
731  self.add_name(name)
732  if resn not in self.db[name]:
733  self.db[name][resn] = {}
734 
735  def add_resolution(self, name, resn, resolution):
736  resn = int(resn)
737  resolution = float(resolution)
738  self.add_name(name)
739  self.add_residue_number(name, resn)
740  if resolution not in self.db[name][resn]:
741  self.db[name][resn][resolution] = []
742 
743  def add_particles(self, name, resn, resolution, particles):
744  resn = int(resn)
745  resolution = float(resolution)
746  self.add_name(name)
747  self.add_residue_number(name, resn)
748  self.add_resolution(name, resn, resolution)
749  self.db[name][resn][resolution] += particles
750  for p in particles:
751  (rh, prf) = self.get_root_hierarchy(p)
752  self.root_hierarchy_dict[p] = rh
753  self.preroot_fragment_hierarchy_dict[p] = prf
754  self.particle_to_name[p] = name
755  if self.model is None:
756  self.model = particles[0].get_model()
757 
758  def get_model(self):
759  return self.model
760 
761  def get_names(self):
762  names = list(self.db.keys())
763  names.sort()
764  return names
765 
766  def get_particles(self, name, resn, resolution):
767  resn = int(resn)
768  resolution = float(resolution)
769  return self.db[name][resn][resolution]
770 
771  def get_particles_at_closest_resolution(self, name, resn, resolution):
772  resn = int(resn)
773  resolution = float(resolution)
774  closestres = min(self.get_residue_resolutions(name, resn),
775  key=lambda x: abs(float(x) - float(resolution)))
776  return self.get_particles(name, resn, closestres)
777 
778  def get_residue_resolutions(self, name, resn):
779  resn = int(resn)
780  resolutions = list(self.db[name][resn].keys())
781  resolutions.sort()
782  return resolutions
783 
784  def get_molecule_resolutions(self, name):
785  resolutions = set()
786  for resn in self.db[name]:
787  resolutions.update(list(self.db[name][resn].keys()))
788  resolutions.sort()
789  return resolutions
790 
791  def get_residue_numbers(self, name):
792  residue_numbers = list(self.db[name].keys())
793  residue_numbers.sort()
794  return residue_numbers
795 
796  def get_particles_by_resolution(self, name, resolution):
797  resolution = float(resolution)
798  particles = []
799  for resn in self.get_residue_numbers(name):
800  result = self.get_particles_at_closest_resolution(
801  name,
802  resn,
803  resolution)
804  pstemp = [p for p in result if p not in particles]
805  particles += pstemp
806  return particles
807 
808  def get_all_particles_by_resolution(self, resolution):
809  resolution = float(resolution)
810  particles = []
811  for name in self.get_names():
812  particles += self.get_particles_by_resolution(name, resolution)
813  return particles
814 
815  def get_root_hierarchy(self, particle):
816  prerootfragment = particle
817  while IMP.atom.Atom.get_is_setup(particle) or \
818  IMP.atom.Residue.get_is_setup(particle) or \
820  if IMP.atom.Atom.get_is_setup(particle):
821  p = IMP.atom.Atom(particle).get_parent()
822  elif IMP.atom.Residue.get_is_setup(particle):
823  p = IMP.atom.Residue(particle).get_parent()
824  elif IMP.atom.Fragment.get_is_setup(particle):
825  p = IMP.atom.Fragment(particle).get_parent()
826  prerootfragment = particle
827  particle = p
828  return (
829  (IMP.atom.Hierarchy(particle), IMP.atom.Hierarchy(prerootfragment))
830  )
831 
832  def get_all_root_hierarchies_by_resolution(self, resolution):
833  hierarchies = []
834  resolution = float(resolution)
835  particles = self.get_all_particles_by_resolution(resolution)
836  for p in particles:
837  rh = self.root_hierarchy_dict[p]
838  if rh not in hierarchies:
839  hierarchies.append(IMP.atom.Hierarchy(rh))
840  return hierarchies
841 
842  def get_preroot_fragments_by_resolution(self, name, resolution):
843  fragments = []
844  resolution = float(resolution)
845  particles = self.get_particles_by_resolution(name, resolution)
846  for p in particles:
847  fr = self.preroot_fragment_hierarchy_dict[p]
848  if fr not in fragments:
849  fragments.append(fr)
850  return fragments
851 
852  def show(self, name):
853  print(name)
854  for resn in self.get_residue_numbers(name):
855  print(resn)
856  for resolution in self.get_residue_resolutions(name, resn):
857  print("----", resolution)
858  for p in self.get_particles(name, resn, resolution):
859  print("--------", p.get_name())
860 
861 
862 def get_prot_name_from_particle(p, list_of_names):
863  '''Return the component name provided a particle and a list of names'''
864  root = p
865  protname = root.get_name()
866  is_a_bead = False
867  while not protname in list_of_names:
868  root0 = root.get_parent()
869  if root0 == IMP.atom.Hierarchy():
870  return (None, None)
871  protname = root0.get_name()
872 
873  # check if that is a bead
874  # this piece of code might be dangerous if
875  # the hierarchy was called Bead :)
876  if "Beads" in protname:
877  is_a_bead = True
878  root = root0
879  return (protname, is_a_bead)
880 
881 
883  '''
884  Retrieve the residue indexes for the given particle.
885 
886  The particle must be an instance of Fragment,Residue, Atom or Molecule
887  or else returns an empty list
888  '''
889  resind = []
891  resind = IMP.atom.Fragment(hier).get_residue_indexes()
893  resind = [IMP.atom.Residue(hier).get_index()]
894  elif IMP.atom.Atom.get_is_setup(hier):
895  a = IMP.atom.Atom(hier)
896  resind = [IMP.atom.Residue(a.get_parent()).get_index()]
898  resind_tmp=IMP.pmi.tools.OrderedSet()
899  for lv in IMP.atom.get_leaves(hier):
903  for ind in get_residue_indexes(lv): resind_tmp.add(ind)
904  resind=list(resind_tmp)
905  else:
906  resind = []
907  return resind
908 
909 
910 def sort_by_residues(particles):
911  particles_residues = [(p, IMP.pmi.tools.get_residue_indexes(p))
912  for p in particles]
913  sorted_particles_residues = sorted(
914  particles_residues,
915  key=lambda tup: tup[1])
916  particles = [p[0] for p in sorted_particles_residues]
917  return particles
918 
919 #
920 # Parallel Computation
921 #
922 
923 
925  """Synchronize data over a parallel run"""
926  from mpi4py import MPI
927  comm = MPI.COMM_WORLD
928  rank = comm.Get_rank()
929  number_of_processes = comm.size
930  comm.Barrier()
931  if rank != 0:
932  comm.send(data, dest=0, tag=11)
933 
934  elif rank == 0:
935  for i in range(1, number_of_processes):
936  data_tmp = comm.recv(source=i, tag=11)
937  if isinstance(data, list):
938  data += data_tmp
939  elif isinstance(data, dict):
940  data.update(data_tmp)
941  else:
942  raise TypeError("data not supported, use list or dictionaries")
943 
944  for i in range(1, number_of_processes):
945  comm.send(data, dest=i, tag=11)
946 
947  if rank != 0:
948  data = comm.recv(source=0, tag=11)
949  return data
950 
951 #
952 ### Lists and iterators
953 #
954 
955 def sublist_iterator(l, lmin=None, lmax=None):
956  '''
957  Yield all sublists of length >= lmin and <= lmax
958  '''
959  if lmin is None:
960  lmin = 0
961  if lmax is None:
962  lmax = len(l)
963  n = len(l) + 1
964  for i in range(n):
965  for j in range(i + 1, n):
966  if len(l[i:j]) <= lmax and len(l[i:j]) >= lmin:
967  yield l[i:j]
968 
969 
970 def flatten_list(l):
971  return [item for sublist in l for item in sublist]
972 
973 
974 def list_chunks_iterator(list, length):
975  """ Yield successive length-sized chunks from a list.
976  """
977  for i in range(0, len(list), length):
978  yield list[i:i + length]
979 
980 
981 def chunk_list_into_segments(seq, num):
982  seq = list(seq)
983  avg = len(seq) / float(num)
984  out = []
985  last = 0.0
986 
987  while last < len(seq):
988  out.append(seq[int(last):int(last + avg)])
989  last += avg
990 
991  return out
992 
993 
994 class Segments(object):
995 
996  ''' This class stores integers
997  in ordered compact lists eg:
998  [[1,2,3],[6,7,8]]
999  the methods help splitting and merging the internal lists
1000  Example:
1001  s=Segments([1,2,3]) is [[1,2,3]]
1002  s.add(4) is [[1,2,3,4]] (add right)
1003  s.add(3) is [[1,2,3,4]] (item already existing)
1004  s.add(7) is [[1,2,3,4],[7]] (new list)
1005  s.add([8,9]) is [[1,2,3,4],[7,8,9]] (add item right)
1006  s.add([5,6]) is [[1,2,3,4,5,6,7,8,9]] (merge)
1007  s.remove(3) is [[1,2],[4,5,6,7,8,9]] (split)
1008  etc.
1009  '''
1010 
1011  def __init__(self,index):
1012  '''index can be a integer or a list of integers '''
1013  if isinstance(index, int):
1014  self.segs=[[index]]
1015  elif isinstance(index, list):
1016  self.segs=[[index[0]]]
1017  for i in index[1:]:
1018  self.add(i)
1019  else:
1020  raise TypeError("index must be an int or list of ints")
1021 
1022  def add(self,index):
1023  '''index can be a integer or a list of integers '''
1024  if isinstance(index, int):
1025  mergeleft=None
1026  mergeright=None
1027  for n,s in enumerate(self.segs):
1028  if index in s:
1029  return 0
1030  else:
1031  if s[0]-index==1:
1032  mergeleft=n
1033  if index-s[-1]==1:
1034  mergeright=n
1035  if mergeright is None and mergeleft is None:
1036  self.segs.append([index])
1037  if not mergeright is None and mergeleft is None:
1038  self.segs[mergeright].append(index)
1039  if not mergeleft is None and mergeright is None:
1040  self.segs[mergeleft]=[index]+self.segs[mergeleft]
1041  if not mergeleft is None and not mergeright is None:
1042  self.segs[mergeright]=self.segs[mergeright]+[index]+self.segs[mergeleft]
1043  del self.segs[mergeleft]
1044 
1045  for n in range(len(self.segs)):
1046  self.segs[n].sort()
1047 
1048  self.segs.sort(key=lambda tup: tup[0])
1049 
1050  elif isinstance(index, list):
1051  for i in index:
1052  self.add(i)
1053  else:
1054  raise TypeError("index must be an int or list of ints")
1055 
1056  def remove(self,index):
1057  '''index can be a integer'''
1058  for n,s in enumerate(self.segs):
1059  if index in s:
1060  if s[0]==index:
1061  self.segs[n]=s[1:]
1062  elif s[-1]==index:
1063  self.segs[n]=s[:-1]
1064  else:
1065  i=self.segs[n].index(index)
1066  self.segs[n]=s[:i]
1067  self.segs.append(s[i+1:])
1068  for n in range(len(self.segs)):
1069  self.segs[n].sort()
1070  if len(self.segs[n])==0:
1071  del self.segs[n]
1072  self.segs.sort(key=lambda tup: tup[0])
1073 
1074  def get_flatten(self):
1075  ''' Returns a flatten list '''
1076  return [item for sublist in self.segs for item in sublist]
1077 
1078  def __repr__(self):
1079  ret_tmp="["
1080  for seg in self.segs:
1081  ret_tmp+=str(seg[0])+"-"+str(seg[-1])+","
1082  ret=ret_tmp[:-1]+"]"
1083  return ret
1084 
1085 #
1086 # Tools to simulate data
1087 #
1088 
1089 def normal_density_function(expected_value, sigma, x):
1090  return (
1091  1 / math.sqrt(2 * math.pi) / sigma *
1092  math.exp(-(x - expected_value) ** 2 / 2 / sigma / sigma)
1093  )
1094 
1095 
1096 def log_normal_density_function(expected_value, sigma, x):
1097  return (
1098  1 / math.sqrt(2 * math.pi) / sigma / x *
1099  math.exp(-(math.log(x / expected_value) ** 2 / 2 / sigma / sigma))
1100  )
1101 
1102 
1103 def get_random_residue_pairs(representation, resolution,
1104  number,
1105  max_distance=None,
1106  avoid_same_particles=False,
1107  names=None):
1108 
1109  particles = []
1110  if names is None:
1111  names=list(representation.hier_dict.keys())
1112 
1113  for name in names:
1114  prot = representation.hier_dict[name]
1115  particles += select(representation,name=name,resolution=resolution)
1116  random_residue_pairs = []
1117  while len(random_residue_pairs)<=number:
1118  p1 = random.choice(particles)
1119  p2 = random.choice(particles)
1120  if max_distance is not None and \
1121  core.get_distance(core.XYZ(p1), core.XYZ(p2)) > max_distance:
1122  continue
1123  r1 = random.choice(IMP.pmi.tools.get_residue_indexes(p1))
1124  r2 = random.choice(IMP.pmi.tools.get_residue_indexes(p2))
1125  if r1==r2 and avoid_same_particles: continue
1126  name1 = representation.get_prot_name_from_particle(p1)
1127  name2 = representation.get_prot_name_from_particle(p2)
1128  random_residue_pairs.append((name1, r1, name2, r2))
1129 
1130  return random_residue_pairs
1131 
1132 def print_multicolumn(list_of_strings, ncolumns=2, truncate=40):
1133 
1134  l = list_of_strings
1135 
1136  cols = ncolumns
1137  # add empty entries after l
1138  for i in range(len(l) % cols):
1139  l.append(" ")
1140 
1141  split = [l[i:i + len(l) / cols] for i in range(0, len(l), len(l) / cols)]
1142  for row in zip(*split):
1143  print("".join(str.ljust(i, truncate) for i in row))
1144 
1145 class ColorChange(object):
1146  '''Change color code to hexadecimal to rgb'''
1147  def __init__(self):
1148  self._NUMERALS = '0123456789abcdefABCDEF'
1149  self._HEXDEC = dict((v, int(v, 16)) for v in (x+y for x in self._NUMERALS for y in self._NUMERALS))
1150  self.LOWERCASE, self.UPPERCASE = 'x', 'X'
1151 
1152  def rgb(self,triplet):
1153  return float(self._HEXDEC[triplet[0:2]]), float(self._HEXDEC[triplet[2:4]]), float(self._HEXDEC[triplet[4:6]])
1154 
1155  def triplet(self,rgb, lettercase=None):
1156  if lettercase is None: lettercase=self.LOWERCASE
1157  return format(rgb[0]<<16 | rgb[1]<<8 | rgb[2], '06'+lettercase)
1158 
1159 # -------------- Collections --------------- #
1160 
1161 class OrderedSet(collections.MutableSet):
1162 
1163  def __init__(self, iterable=None):
1164  self.end = end = []
1165  end += [None, end, end] # sentinel node for doubly linked list
1166  self.map = {} # key --> [key, prev, next]
1167  if iterable is not None:
1168  self |= iterable
1169 
1170  def __len__(self):
1171  return len(self.map)
1172 
1173  def __contains__(self, key):
1174  return key in self.map
1175 
1176  def add(self, key):
1177  if key not in self.map:
1178  end = self.end
1179  curr = end[1]
1180  curr[2] = end[1] = self.map[key] = [key, curr, end]
1181 
1182  def discard(self, key):
1183  if key in self.map:
1184  key, prev, next = self.map.pop(key)
1185  prev[2] = next
1186  next[1] = prev
1187 
1188  def __iter__(self):
1189  end = self.end
1190  curr = end[2]
1191  while curr is not end:
1192  yield curr[0]
1193  curr = curr[2]
1194 
1195  def __reversed__(self):
1196  end = self.end
1197  curr = end[1]
1198  while curr is not end:
1199  yield curr[0]
1200  curr = curr[1]
1201 
1202  def pop(self, last=True):
1203  if not self:
1204  raise KeyError('set is empty')
1205  if last:
1206  key = self.end[1][0]
1207  else:
1208  key = self.end[2][0]
1209  self.discard(key)
1210  return key
1211 
1212  def __repr__(self):
1213  if not self:
1214  return '%s()' % (self.__class__.__name__,)
1215  return '%s(%r)' % (self.__class__.__name__, list(self))
1216 
1217  def __eq__(self, other):
1218  if isinstance(other, OrderedSet):
1219  return len(self) == len(other) and list(self) == list(other)
1220  return set(self) == set(other)
1221 
1222 
1223 class OrderedDefaultDict(OrderedDict):
1224  """Store objects in order they were added, but with default type.
1225  Source: http://stackoverflow.com/a/4127426/2608793
1226  """
1227  def __init__(self, *args, **kwargs):
1228  if not args:
1229  self.default_factory = None
1230  else:
1231  if not (args[0] is None or callable(args[0])):
1232  raise TypeError('first argument must be callable or None')
1233  self.default_factory = args[0]
1234  args = args[1:]
1235  super(OrderedDefaultDict, self).__init__(*args, **kwargs)
1236 
1237  def __missing__ (self, key):
1238  if self.default_factory is None:
1239  raise KeyError(key)
1240  self[key] = default = self.default_factory()
1241  return default
1242 
1243  def __reduce__(self): # optional, for pickle support
1244  args = (self.default_factory,) if self.default_factory else ()
1245  if sys.version_info[0] >= 3:
1246  return self.__class__, args, None, None, self.items()
1247  else:
1248  return self.__class__, args, None, None, self.iteritems()
1249 
1250 
1251 # -------------- PMI2 Tools --------------- #
1252 
1253 def set_coordinates_from_rmf(hier,rmf_fn,frame_num=0):
1254  """Extract frame from RMF file and fill coordinates. Must be identical topology.
1255  @param hier The (System) hierarchy to fill (e.g. after you've built it)
1256  @param rmf_fn The file to extract from
1257  @param frame_num The frame number to extract
1258  """
1259  rh = RMF.open_rmf_file_read_only(rmf_fn)
1260  IMP.rmf.link_hierarchies(rh,[hier])
1261  IMP.rmf.load_frame(rh, RMF.FrameID(frame_num))
1262  del rh
1263 
1264 def input_adaptor(stuff,
1265  pmi_resolution=0,
1266  flatten=False,
1267  selection_tuple=None,
1268  warn_about_slices=True):
1269  """Adapt things for PMI (degrees of freedom, restraints, ...)
1270  Returns list of list of hierarchies, separated into Molecules if possible.
1271  The input can be a list, or a list of lists (iterable of ^1 or iterable of ^2)
1272  (iterable of ^2) Hierarchy -> returns input as list of list of hierarchies,
1273  only one entry, not grouped by molecules.
1274  (iterable of ^2) PMI::System/State/Molecule/TempResidue ->
1275  returns residue hierarchies, grouped in molecules, at requested resolution
1276  @param stuff Can be one of the following inputs:
1277  IMP Hierarchy, PMI System/State/Molecule/TempResidue, or a list/set (of list/set) of them.
1278  Must be uniform input, however. No mixing object types.
1279  @param pmi_resolution For selecting, only does it if you pass PMI objects. Set it to "all"
1280  if you want all resolutions!
1281  @param flatten Set to True if you just want all hierarchies in one list.
1282  @param warn_about_slices Print a warning if you are requesting only part
1283  of a bead. Sometimes you just don't care!
1284  @note since this relies on IMP::atom::Selection, this will not return
1285  any objects if they weren't built! But there should be no problem
1286  if you request unbuilt residues - they should be ignored.
1287  """
1288 
1289  if stuff is None:
1290  return stuff
1291 
1292  if hasattr(stuff,'__iter__'):
1293  if len(stuff)==0:
1294  return stuff
1295  thelist=list(stuff)
1296 
1297  # iter of iter of should be ok
1298  if all(hasattr(el,'__iter__') for el in thelist):
1299  thelist = [i for sublist in thelist for i in sublist]
1300  elif any(hasattr(el,'__iter__') for el in thelist):
1301  raise Exception('input_adaptor: input_object must be a list or a list of lists')
1302 
1303  stuff = thelist
1304  else:
1305  stuff = [stuff]
1306 
1307  # check that it is a hierarchy homogenously:
1308  try:
1309  is_hierarchy=all(IMP.atom.Hierarchy.get_is_setup(s) for s in stuff)
1310  except (NotImplementedError, TypeError):
1311  is_hierarchy=False
1312  # get the other types homogenously
1313  is_system=all(isinstance(s, IMP.pmi.topology.System) for s in stuff)
1314  is_state=all(isinstance(s, IMP.pmi.topology.State) for s in stuff)
1315  is_molecule=all(isinstance(s, IMP.pmi.topology.Molecule) for s in stuff)
1316  is_temp_residue=all(isinstance(s, IMP.pmi.topology.TempResidue) for s in stuff)
1317 
1318  # now that things are ok, do selection if requested
1319  hier_list = []
1320  pmi_input = False
1321  if is_system or is_state or is_molecule or is_temp_residue:
1322  # if PMI, perform selection using gathered indexes
1323  pmi_input = True
1324  indexes_per_mol = OrderedDefaultDict(list) #key is Molecule object, value are residues
1325  if is_system:
1326  for system in stuff:
1327  for state in system.get_states():
1328  mdict = state.get_molecules()
1329  for molname in mdict:
1330  for copy in mdict[molname]:
1331  indexes_per_mol[copy] += [r.get_index() for r in copy.get_residues()]
1332  elif is_state:
1333  for state in stuff:
1334  mdict = state.get_molecules()
1335  for molname in mdict:
1336  for copy in mdict[molname]:
1337  indexes_per_mol[copy] += [r.get_index() for r in copy.get_residues()]
1338  elif is_molecule:
1339  for molecule in stuff:
1340  indexes_per_mol[molecule] += [r.get_index() for r in molecule.get_residues()]
1341  else: # is_temp_residue
1342  for tempres in stuff:
1343  indexes_per_mol[tempres.get_molecule()].append(tempres.get_index())
1344  for mol in indexes_per_mol:
1345  if pmi_resolution=='all':
1346  # because you select from the molecule,
1347  # this will start the search from the base resolution
1348  ps = select_at_all_resolutions(mol.get_hierarchy(),
1349  residue_indexes=indexes_per_mol[mol])
1350  else:
1351  sel = IMP.atom.Selection(mol.get_hierarchy(),
1352  resolution=pmi_resolution,
1353  residue_indexes=indexes_per_mol[mol])
1354  ps = sel.get_selected_particles()
1355 
1356  # check that you don't have any incomplete fragments!
1357  if warn_about_slices:
1358  rset = set(indexes_per_mol[mol])
1359  for p in ps:
1361  fset = set(IMP.atom.Fragment(p).get_residue_indexes())
1362  if not fset <= rset:
1363  minset = min(fset)
1364  maxset = max(fset)
1365  found = fset&rset
1366  minf = min(found)
1367  maxf = max(found)
1368  resbreak = maxf if minf==minset else minset-1
1369  warnings.warn(
1370  'You are trying to select only part of the '
1371  'bead %s:%i-%i. The residues you requested '
1372  'are %i-%i. You can fix this by: '
1373  '1) requesting the whole bead/none of it; or'
1374  '2) break the bead up by passing '
1375  'bead_extra_breaks=[\'%i\'] in '
1376  'molecule.add_representation()'
1377  %(mol.get_name(), minset, maxset, minf, maxf,
1378  resbreak), IMP.pmi.ParameterWarning)
1379  hier_list.append([IMP.atom.Hierarchy(p) for p in ps])
1380  elif is_hierarchy:
1381  #check
1382  ps=[]
1383  if pmi_resolution=='all':
1384  for h in stuff:
1386  else:
1387  for h in stuff:
1388  ps+=IMP.atom.Selection(h,resolution=pmi_resolution).get_selected_particles()
1389  hier_list=[IMP.atom.Hierarchy(p) for p in ps]
1390  if not flatten:
1391  hier_list = [hier_list]
1392  else:
1393  raise Exception('input_adaptor: you passed something of wrong type or a list with mixed types')
1394 
1395  if flatten and pmi_input:
1396  return [h for sublist in hier_list for h in sublist]
1397  else:
1398  return hier_list
1399 
1400 
1402  """Returns sequence-sorted segments array, each containing the first particle
1403  the last particle and the first residue index."""
1404 
1405  from operator import itemgetter
1406  hiers=IMP.pmi.tools.input_adaptor(mol)
1407  if len(hiers)>1:
1408  raise ValueError("only pass stuff from one Molecule, please")
1409  hiers = hiers[0]
1410  segs = []
1411  for h in hiers:
1412  try:
1413  start = IMP.atom.Hierarchy(h).get_children()[0]
1414  except IndexError:
1415  start = IMP.atom.Hierarchy(h)
1416 
1417  try:
1418  end = IMP.atom.Hierarchy(h).get_children()[-1]
1419  except IndexError:
1420  end = IMP.atom.Hierarchy(h)
1421 
1422  startres = IMP.pmi.tools.get_residue_indexes(start)[0]
1423  endres = IMP.pmi.tools.get_residue_indexes(end)[-1]
1424  segs.append((start, end, startres))
1425  return sorted(segs, key=itemgetter(2))
1426 
1427 def display_bonds(mol):
1428  """Decorate the sequence-consecutive particles from a PMI2 molecule with a bond,
1429  so that they appear connected in the rmf file"""
1430  SortedSegments=get_sorted_segments(mol)
1431  for x in range(len(SortedSegments) - 1):
1432 
1433  last = SortedSegments[x][1]
1434  first = SortedSegments[x + 1][0]
1435 
1436  p1 = last.get_particle()
1437  p2 = first.get_particle()
1438  if not IMP.atom.Bonded.get_is_setup(p1):
1440  if not IMP.atom.Bonded.get_is_setup(p2):
1442 
1445  IMP.atom.Bonded(p1),
1446  IMP.atom.Bonded(p2),1)
1447 
1448 
1449 class ThreeToOneConverter(defaultdict):
1450  """This class converts three to one letter codes, and return X for any unknown codes"""
1451  def __init__(self,is_nucleic=False):
1452 
1453  if not is_nucleic:
1454  threetoone = {'ALA': 'A', 'ARG': 'R', 'ASN': 'N', 'ASP': 'D',
1455  'CYS': 'C', 'GLU': 'E', 'GLN': 'Q', 'GLY': 'G',
1456  'HIS': 'H', 'ILE': 'I', 'LEU': 'L', 'LYS': 'K',
1457  'MET': 'M', 'PHE': 'F', 'PRO': 'P', 'SER': 'S',
1458  'THR': 'T', 'TRP': 'W', 'TYR': 'Y', 'VAL': 'V', 'UNK': 'X'}
1459  else:
1460  threetoone = {'ADE': 'A', 'URA': 'U', 'CYT': 'C', 'GUA': 'G',
1461  'THY': 'T', 'UNK': 'X'}
1462 
1463  defaultdict.__init__(self,lambda: "X", threetoone)
1464 
1465 
1466 
1467 
1468 def get_residue_type_from_one_letter_code(code,is_nucleic=False):
1469  threetoone=ThreeToOneConverter(is_nucleic)
1470  one_to_three={}
1471  for k in threetoone:
1472  one_to_three[threetoone[k]] = k
1473  return IMP.atom.ResidueType(one_to_three[code])
1474 
1475 
1476 def get_all_leaves(list_of_hs):
1477  """ Just get the leaves from a list of hierarchies """
1478  lvs = list(itertools.chain.from_iterable(IMP.atom.get_leaves(item) for item in list_of_hs))
1479  return lvs
1480 
1481 
1483  hiers=None,
1484  **kwargs):
1485  """Perform selection using the usual keywords but return ALL resolutions (BEADS and GAUSSIANS).
1486  Returns in flat list!
1487  """
1488 
1489  if hiers is None:
1490  hiers = []
1491  if hier is not None:
1492  hiers.append(hier)
1493  if len(hiers)==0:
1494  warnings.warn("You passed nothing to select_at_all_resolutions()",
1496  return []
1497  ret = OrderedSet()
1498  for hsel in hiers:
1499  try:
1500  htest = IMP.atom.Hierarchy.get_is_setup(hsel)
1501  except:
1502  raise Exception('select_at_all_resolutions: you have to pass an IMP Hierarchy')
1503  if not htest:
1504  raise Exception('select_at_all_resolutions: you have to pass an IMP Hierarchy')
1505  if 'resolution' in kwargs or 'representation_type' in kwargs:
1506  raise Exception("don't pass resolution or representation_type to this function")
1507  selB = IMP.atom.Selection(hsel,resolution=IMP.atom.ALL_RESOLUTIONS,
1508  representation_type=IMP.atom.BALLS,**kwargs)
1509  selD = IMP.atom.Selection(hsel,resolution=IMP.atom.ALL_RESOLUTIONS,
1510  representation_type=IMP.atom.DENSITIES,**kwargs)
1511  ret |= OrderedSet(selB.get_selected_particles())
1512  ret |= OrderedSet(selD.get_selected_particles())
1513  return list(ret)
1514 
1515 
1517  target_ps,
1518  sel_zone,
1519  entire_residues,
1520  exclude_backbone):
1521  """Utility to retrieve particles from a hierarchy within a
1522  zone around a set of ps.
1523  @param hier The hierarchy in which to look for neighbors
1524  @param target_ps The particles for zoning
1525  @param sel_zone The maximum distance
1526  @param entire_residues If True, will grab entire residues
1527  @param exclude_backbone If True, will only return sidechain particles
1528  """
1529 
1530  test_sel = IMP.atom.Selection(hier)
1531  backbone_types=['C','N','CB','O']
1532  if exclude_backbone:
1533  test_sel -= IMP.atom.Selection(hier,atom_types=[IMP.atom.AtomType(n)
1534  for n in backbone_types])
1535  test_ps = test_sel.get_selected_particles()
1536  nn = IMP.algebra.NearestNeighbor3D([IMP.core.XYZ(p).get_coordinates()
1537  for p in test_ps])
1538  zone = set()
1539  for target in target_ps:
1540  zone|=set(nn.get_in_ball(IMP.core.XYZ(target).get_coordinates(),sel_zone))
1541  zone_ps = [test_ps[z] for z in zone]
1542  if entire_residues:
1543  final_ps = set()
1544  for z in zone_ps:
1545  final_ps|=set(IMP.atom.Hierarchy(z).get_parent().get_children())
1546  zone_ps = [h.get_particle() for h in final_ps]
1547  return zone_ps
1548 
1549 
1551  """Returns unique objects in original order"""
1552  rbs = set()
1553  beads = []
1554  rbs_ordered = []
1555  if not hasattr(hiers,'__iter__'):
1556  hiers = [hiers]
1557  for p in get_all_leaves(hiers):
1559  rb = IMP.core.RigidMember(p).get_rigid_body()
1560  if rb not in rbs:
1561  rbs.add(rb)
1562  rbs_ordered.append(rb)
1564  rb = IMP.core.NonRigidMember(p).get_rigid_body()
1565  if rb not in rbs:
1566  rbs.add(rb)
1567  rbs_ordered.append(rb)
1568  beads.append(p)
1569  else:
1570  beads.append(p)
1571  return rbs_ordered,beads
1572 
1573 def get_molecules(input_objects):
1574  "This function returns the parent molecule hierarchies of given objects"
1575  stuff=input_adaptor(input_objects, pmi_resolution='all',flatten=True)
1576  molecules=set()
1577  for h in stuff:
1578  is_root=False
1579  is_molecule=False
1580  while not (is_molecule or is_root):
1581  root=IMP.atom.get_root(h)
1582  if root == h:
1583  is_root=True
1584  is_molecule=IMP.atom.Molecule.get_is_setup(h)
1585  if is_molecule:
1586  molecules.add(IMP.atom.Molecule(h))
1587  h=h.get_parent()
1588  return list(molecules)
1589 
1590 def get_molecules_dictionary(input_objects):
1591  moldict=defaultdict(list)
1592  for mol in IMP.pmi.tools.get_molecules(input_objects):
1593  name=mol.get_name()
1594  moldict[name].append(mol)
1595 
1596  for mol in moldict:
1597  moldict[mol].sort(key=lambda x: IMP.atom.Copy(x).get_copy_index())
1598  return moldict
1599 
1600 def get_molecules_dictionary_by_copy(input_objects):
1601  moldict=defaultdict(dict)
1602  for mol in IMP.pmi.tools.get_molecules(input_objects):
1603  name=mol.get_name()
1604  c=IMP.atom.Copy(mol).get_copy_index()
1605  moldict[name][c]=mol
1606  return moldict
1607 
1608 def get_selections_dictionary(input_objects):
1609  moldict=IMP.pmi.tools.get_molecules_dictionary(input_objects)
1610  seldict=defaultdict(list)
1611  for name, mols in moldict.items():
1612  for m in mols:
1613  seldict[name].append(IMP.atom.Selection(m))
1614  return seldict
1615 
1616 def get_densities(input_objects):
1617  """Given a list of PMI objects, returns all density hierarchies within
1618  these objects. The output of this function can be inputted into
1619  things such as EM restraints. This function is intended to gather density particles
1620  appended to molecules (and not other hierarchies which might have been appended to the root node directly).
1621  """
1622  # Note that Densities can only be selected at the Root or Molecule level and not at the Leaves level.
1623  # we'll first get all molecule hierarchies corresponding to the leaves.
1624  molecules=get_molecules(input_objects)
1625  densities=[]
1626  for i in molecules:
1627  densities+=IMP.atom.Selection(i,representation_type=IMP.atom.DENSITIES).get_selected_particles()
1628  return densities
1629 
1630 def shuffle_configuration(objects,
1631  max_translation=300., max_rotation=2.0 * pi,
1632  avoidcollision_rb=True, avoidcollision_fb=False,
1633  cutoff=10.0, niterations=100,
1634  bounding_box=None,
1635  excluded_rigid_bodies=[],
1636  hierarchies_excluded_from_collision=[],
1637  hierarchies_included_in_collision=[],
1638  verbose=False,
1639  return_debug=False):
1640  """Shuffle particles. Used to restart the optimization.
1641  The configuration of the system is initialized by placing each
1642  rigid body and each bead randomly in a box. If `bounding_box` is
1643  specified, the particles are placed inside this box; otherwise, each
1644  particle is displaced by up to max_translation angstroms, and randomly
1645  rotated. Effort is made to place particles far enough from each other to
1646  prevent any steric clashes.
1647  @param objects Can be one of the following inputs:
1648  IMP Hierarchy, PMI System/State/Molecule/TempResidue, or a list/set of them
1649  @param max_translation Max translation (rbs and flexible beads)
1650  @param max_rotation Max rotation (rbs only)
1651  @param avoidcollision_rb check if the particle/rigid body was
1652  placed close to another particle; uses the optional
1653  arguments cutoff and niterations
1654  @param avoidcollision_fb Advanced. Generally you want this False because it's hard to shuffle beads.
1655  @param cutoff Distance less than this is a collision
1656  @param niterations How many times to try avoiding collision
1657  @param bounding_box Only shuffle particles within this box. Defined by ((x1,y1,z1),(x2,y2,z2)).
1658  @param excluded_rigid_bodies Don't shuffle these rigid body objects
1659  @param hierarchies_excluded_from_collision Don't count collision with these bodies
1660  @param hierarchies_included_in_collision Hierarchies that are not shuffled, but should be included in collision calculation (for fixed regions)
1661  @param verbose Give more output
1662  @note Best to only call this function after you've set up degrees of freedom
1663  For debugging purposes, returns: <shuffled indexes>, <collision avoided indexes>
1664  """
1665 
1666  ### checking input
1667  hierarchies = IMP.pmi.tools.input_adaptor(objects,
1668  pmi_resolution='all',
1669  flatten=True)
1670  rigid_bodies,flexible_beads = get_rbs_and_beads(hierarchies)
1671  if len(rigid_bodies)>0:
1672  mdl = rigid_bodies[0].get_model()
1673  elif len(flexible_beads)>0:
1674  mdl = flexible_beads[0].get_model()
1675  else:
1676  raise Exception("Could not find any particles in the hierarchy")
1677  if len(rigid_bodies) == 0:
1678  print("shuffle_configuration: rigid bodies were not initialized")
1679 
1680  ### gather all particles
1682  gcpf.set_distance(cutoff)
1683 
1684  # Add particles from excluded hierarchies to excluded list
1685  collision_excluded_hierarchies = IMP.pmi.tools.input_adaptor(hierarchies_excluded_from_collision,
1686  pmi_resolution='all',
1687  flatten=True)
1688 
1689  collision_included_hierarchies = IMP.pmi.tools.input_adaptor(hierarchies_included_in_collision,
1690  pmi_resolution='all',
1691  flatten=True)
1692 
1693  collision_excluded_idxs = set([l.get_particle().get_index() for h in collision_excluded_hierarchies \
1694  for l in IMP.core.get_leaves(h)])
1695 
1696  collision_included_idxs = set([l.get_particle().get_index() for h in collision_included_hierarchies \
1697  for l in IMP.core.get_leaves(h)])
1698 
1699  # Excluded collision with Gaussians
1700  all_idxs = [] #expand to representations?
1701  for p in IMP.pmi.tools.get_all_leaves(hierarchies):
1703  all_idxs.append(p.get_particle_index())
1705  collision_excluded_idxs.add(p.get_particle_index())
1706 
1707  if bounding_box is not None:
1708  ((x1, y1, z1), (x2, y2, z2)) = bounding_box
1709  ub = IMP.algebra.Vector3D(x1, y1, z1)
1710  lb = IMP.algebra.Vector3D(x2, y2, z2)
1711  bb = IMP.algebra.BoundingBox3D(ub, lb)
1712 
1713  all_idxs = set(all_idxs) | collision_included_idxs
1714  all_idxs = all_idxs - collision_excluded_idxs
1715  debug = []
1716  print('shuffling', len(rigid_bodies), 'rigid bodies')
1717  for rb in rigid_bodies:
1718  if rb not in excluded_rigid_bodies:
1719  # gather particles to avoid with this transform
1720  if avoidcollision_rb:
1721  rb_idxs = set(rb.get_member_particle_indexes()) - \
1722  collision_excluded_idxs
1723  other_idxs = all_idxs - rb_idxs
1724  if not other_idxs:
1725  continue
1726 
1727  # iterate, trying to avoid collisions
1728  niter = 0
1729  while niter < niterations:
1730  rbxyz = (rb.get_x(), rb.get_y(), rb.get_z())
1731 
1732  # local transform
1733  if bounding_box:
1734  translation = IMP.algebra.get_random_vector_in(bb)
1735  # First move to origin
1737  -IMP.core.XYZ(rb).get_coordinates())
1738  IMP.core.transform(rb, transformation_orig)
1739  old_coord = IMP.core.XYZ(rb).get_coordinates()
1741  transformation = IMP.algebra.Transformation3D(rotation,
1742  translation)
1743 
1744  else:
1746  rbxyz,
1747  max_translation,
1748  max_rotation)
1749 
1750  debug.append([rb, other_idxs if avoidcollision_rb else set()])
1751  IMP.core.transform(rb, transformation)
1752 
1753  # check collisions
1754  if avoidcollision_rb:
1755  mdl.update()
1756  npairs = len(gcpf.get_close_pairs(mdl,
1757  list(other_idxs),
1758  list(rb_idxs)))
1759  #print("NPAIRS:", npairs)
1760  if npairs==0:
1761  break
1762  else:
1763  niter += 1
1764  if verbose:
1765  print("shuffle_configuration: rigid body placed close to other %d particles, trying again..." % npairs)
1766  print("shuffle_configuration: rigid body name: " + rb.get_name())
1767  if niter == niterations:
1768  raise ValueError("tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1769  else:
1770  break
1771 
1772  print('shuffling', len(flexible_beads), 'flexible beads')
1773  for fb in flexible_beads:
1774  # gather particles to avoid
1775  if avoidcollision_fb:
1776  fb_idxs = set(IMP.get_indexes([fb]))
1777  other_idxs = all_idxs - fb_idxs
1778  if not other_idxs:
1779  continue
1780 
1781  # iterate, trying to avoid collisions
1782  niter = 0
1783  while niter < niterations:
1784  if bounding_box:
1785  translation = IMP.algebra.get_random_vector_in(bb)
1786  transformation = IMP.algebra.Transformation3D(translation)
1787  else:
1788  fbxyz = IMP.core.XYZ(fb).get_coordinates()
1789  transformation = IMP.algebra.get_random_local_transformation(fbxyz,
1790  max_translation,
1791  max_rotation)
1792 
1793  # For gaussians, treat this fb as an rb
1795  memb = IMP.core.NonRigidMember(fb)
1796  xyz = memb.get_internal_coordinates()
1797  if bounding_box:
1798  # 'translation' is the new desired position in global
1799  # coordinates; we need to convert that to internal
1800  # coordinates first using the rigid body's ref frame
1801  rf = memb.get_rigid_body().get_reference_frame()
1802  glob_to_int = rf.get_transformation_from()
1803  memb.set_internal_coordinates(
1804  glob_to_int.get_transformed(translation))
1805  else:
1806  xyz_transformed=transformation.get_transformed(xyz)
1807  memb.set_internal_coordinates(xyz_transformed)
1808  debug.append([xyz,other_idxs if avoidcollision_fb else set()])
1809  else:
1810  d =IMP.core.XYZ(fb)
1811  if bounding_box:
1812  # Translate to origin first
1813  IMP.core.transform(d, -d.get_coordinates())
1814  d =IMP.core.XYZ(fb)
1815  debug.append([d,other_idxs if avoidcollision_fb else set()])
1816  if IMP.core.RigidBody.get_is_setup(fb.get_particle()):
1817  IMP.core.transform(IMP.core.RigidBody(fb.get_particle()), transformation)
1818  else:
1819  IMP.core.transform(d, transformation)
1820 
1821  if avoidcollision_fb:
1822  mdl.update()
1823  npairs = len(gcpf.get_close_pairs(mdl,
1824  list(other_idxs),
1825  list(fb_idxs)))
1826 
1827  if npairs==0:
1828  break
1829  else:
1830  niter += 1
1831  print("shuffle_configuration: floppy body placed close to other %d particles, trying again..." % npairs)
1832  if niter == niterations:
1833  raise ValueError("tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1834  else:
1835  break
1836  if return_debug:
1837  return debug
1838 
1839 class ColorHierarchy(object):
1840 
1841  def __init__(self,hier):
1842  import matplotlib as mpl
1843  mpl.use('Agg')
1844  import matplotlib.pyplot as plt
1845  self.mpl = mpl
1846  self.plt = plt
1847 
1848  hier.ColorHierarchy=self
1849  self.hier=hier
1851  self.mols=[IMP.pmi.topology.PMIMoleculeHierarchy(mol) for mol in mols]
1852  self.method=self.nochange
1853  self.scheme=None
1854  self.first=None
1855  self.last=None
1856 
1857  def nochange(self):
1858  pass
1859 
1860  def get_color(self,fl):
1861  return IMP.display.Color(*self.scheme(fl)[0:3])
1862 
1863  def get_log_scale(self,fl):
1864  import math
1865  eps=1.0
1866  return math.log(fl+eps)
1867 
1868  def color_by_resid(self):
1869  self.method=self.color_by_resid
1870  self.scheme=self.mpl.cm.rainbow
1871  for mol in self.mols:
1872  self.first=1
1874  for p in IMP.atom.get_leaves(mol):
1876  ri=IMP.atom.Residue(p).get_index()
1877  c=self.get_color(float(ri)/self.last)
1878  IMP.display.Colored(p).set_color(c)
1881  avr=sum(ris)/len(ris)
1882  c=self.get_color(float(avr)/self.last)
1883  IMP.display.Colored(p).set_color(c)
1884 
1885  def color_by_uncertainty(self):
1886  self.method=self.color_by_uncertainty
1887  self.scheme=self.mpl.cm.jet
1888  ps=IMP.atom.get_leaves(self.hier)
1889  unc_dict={}
1890  for p in ps:
1892  u=IMP.pmi.Uncertainty(p).get_uncertainty()
1893  unc_dict[p]=u
1894  self.first=self.get_log_scale(1.0) #math.log(min(unc_dict.values())+eps)
1895  self.last=self.get_log_scale(100.0) #math.log(max(unc_dict.values())+eps)
1896  for p in unc_dict:
1897  value=self.get_log_scale(unc_dict[p])
1898  if value>=self.last: value=self.last
1899  if value<=self.first: value=self.first
1900  c=self.get_color((value-self.first)/(self.last-self.first))
1901  IMP.display.Colored(p).set_color(c)
1902 
1903  def get_color_bar(self,filename):
1904  import matplotlib as mpl
1905  mpl.use('Agg')
1906  import matplotlib.pyplot as plt
1907  import math
1908  plt.clf()
1909  fig = plt.figure(figsize=(8, 3))
1910  ax1 = fig.add_axes([0.05, 0.80, 0.9, 0.15])
1911 
1912  cmap = self.scheme
1913  norm = mpl.colors.Normalize(vmin=0.0, vmax=1.0)
1914 
1915  if self.method == self.color_by_uncertainty:
1916  angticks=[1.0,2.5,5.0,10.0,25.0,50.0,100.0]
1917  vvalues=[]
1918  marks=[]
1919  for at in angticks:
1920  vvalue=(self.get_log_scale(at)-self.first)/(self.last-self.first)
1921  if vvalue <= 1.0 and vvalue >= 0.0:
1922  vvalues.append(vvalue)
1923  marks.append(str(at))
1924  cb1 = mpl.colorbar.ColorbarBase(ax1, cmap=cmap,
1925  norm=norm,
1926  ticks=vvalues,
1927  orientation='horizontal')
1928  print(self.first,self.last,marks,vvalues)
1929  cb1.ax.set_xticklabels(marks)
1930  cb1.set_label('Angstorm')
1931  plt.savefig(filename, dpi=150, transparent=True)
1932  plt.show()
1933 
1934 
1935 
1936 
1937 def color2rgb(colorname):
1938  """Given a Chimera color name or hex color value, return RGB"""
1939  d = {'aquamarine': (0.4980392156862745, 1.0, 0.8313725490196079),
1940  'black': (0.0, 0.0, 0.0),
1941  'blue': (0.0, 0.0, 1.0),
1942  'brown': (0.6470588235294118, 0.16470588235294117, 0.16470588235294117),
1943  'chartreuse': (0.4980392156862745, 1.0, 0.0),
1944  'coral': (1.0, 0.4980392156862745, 0.3137254901960784),
1945  'cornflower blue': (0.39215686274509803, 0.5843137254901961, 0.9294117647058824),
1946  'cyan': (0.0, 1.0, 1.0),
1947  'dark cyan': (0.0, 0.5450980392156862, 0.5450980392156862),
1948  'dark gray': (0.6627450980392157, 0.6627450980392157, 0.6627450980392157),
1949  'dark green': (0.0, 0.39215686274509803, 0.0),
1950  'dark khaki': (0.7411764705882353, 0.7176470588235294, 0.4196078431372549),
1951  'dark magenta': (0.5450980392156862, 0.0, 0.5450980392156862),
1952  'dark olive green': (0.3333333333333333, 0.4196078431372549, 0.1843137254901961),
1953  'dark red': (0.5450980392156862, 0.0, 0.0),
1954  'dark slate blue': (0.2823529411764706, 0.23921568627450981, 0.5450980392156862),
1955  'dark slate gray': (0.1843137254901961, 0.30980392156862746, 0.30980392156862746),
1956  'deep pink': (1.0, 0.0784313725490196, 0.5764705882352941),
1957  'deep sky blue': (0.0, 0.7490196078431373, 1.0),
1958  'dim gray': (0.4117647058823529, 0.4117647058823529, 0.4117647058823529),
1959  'dodger blue': (0.11764705882352941, 0.5647058823529412, 1.0),
1960  'firebrick': (0.6980392156862745, 0.13333333333333333, 0.13333333333333333),
1961  'forest green': (0.13333333333333333, 0.5450980392156862, 0.13333333333333333),
1962  'gold': (1.0, 0.8431372549019608, 0.0),
1963  'goldenrod': (0.8549019607843137, 0.6470588235294118, 0.12549019607843137),
1964  'gray': (0.7450980392156863, 0.7450980392156863, 0.7450980392156863),
1965  'green': (0.0, 1.0, 0.0),
1966  'hot pink': (1.0, 0.4117647058823529, 0.7058823529411765),
1967  'khaki': (0.9411764705882353, 0.9019607843137255, 0.5490196078431373),
1968  'light blue': (0.6784313725490196, 0.8470588235294118, 0.9019607843137255),
1969  'light gray': (0.8274509803921568, 0.8274509803921568, 0.8274509803921568),
1970  'light green': (0.5647058823529412, 0.9333333333333333, 0.5647058823529412),
1971  'light sea green': (0.12549019607843137, 0.6980392156862745, 0.6666666666666666),
1972  'lime green': (0.19607843137254902, 0.803921568627451, 0.19607843137254902),
1973  'magenta': (1.0, 0.0, 1.0),
1974  'medium blue': (0.19607843137254902, 0.19607843137254902, 0.803921568627451),
1975  'medium purple': (0.5764705882352941, 0.4392156862745098, 0.8588235294117647),
1976  'navy blue': (0.0, 0.0, 0.5019607843137255),
1977  'olive drab': (0.4196078431372549, 0.5568627450980392, 0.13725490196078433),
1978  'orange red': (1.0, 0.27058823529411763, 0.0),
1979  'orange': (1.0, 0.4980392156862745, 0.0),
1980  'orchid': (0.8549019607843137, 0.4392156862745098, 0.8392156862745098),
1981  'pink': (1.0, 0.7529411764705882, 0.796078431372549),
1982  'plum': (0.8666666666666667, 0.6274509803921569, 0.8666666666666667),
1983  'purple': (0.6274509803921569, 0.12549019607843137, 0.9411764705882353),
1984  'red': (1.0, 0.0, 0.0),
1985  'rosy brown': (0.7372549019607844, 0.5607843137254902, 0.5607843137254902),
1986  'salmon': (0.9803921568627451, 0.5019607843137255, 0.4470588235294118),
1987  'sandy brown': (0.9568627450980393, 0.6431372549019608, 0.3764705882352941),
1988  'sea green': (0.1803921568627451, 0.5450980392156862, 0.3411764705882353),
1989  'sienna': (0.6274509803921569, 0.3215686274509804, 0.17647058823529413),
1990  'sky blue': (0.5294117647058824, 0.807843137254902, 0.9215686274509803),
1991  'slate gray': (0.4392156862745098, 0.5019607843137255, 0.5647058823529412),
1992  'spring green': (0.0, 1.0, 0.4980392156862745),
1993  'steel blue': (0.27450980392156865, 0.5098039215686274, 0.7058823529411765),
1994  'tan': (0.8235294117647058, 0.7058823529411765, 0.5490196078431373),
1995  'turquoise': (0.25098039215686274, 0.8784313725490196, 0.8156862745098039),
1996  'violet red': (0.8156862745098039, 0.12549019607843137, 0.5647058823529412),
1997  'white': (1.0, 1.0, 1.0),
1998  'yellow': (1.0, 1.0, 0.0)}
1999  if colorname.startswith('#'):
2000  return tuple(int(colorname[i:i+2], 16) / 255. for i in (1, 3, 5))
2001  else:
2002  return d[colorname]
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: Molecule.h:35
def list_chunks_iterator
Yield successive length-sized chunks from a list.
Definition: tools.py:974
Simple 3D transformation class.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: Residue.h:155
def select_at_all_resolutions
Perform selection using the usual keywords but return ALL resolutions (BEADS and GAUSSIANS).
Definition: tools.py:1482
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 display_bonds
Decorate the sequence-consecutive particles from a PMI2 molecule with a bond, so that they appear con...
Definition: tools.py:1427
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:105
Store the representations for a system.
Definition: tools.py:714
def remove
index can be a integer
Definition: tools.py:1056
def get_random_cross_link_dataset
Return a random cross-link dataset as a string.
Definition: tools.py:248
def shuffle_configuration
Shuffle particles.
Definition: tools.py:1644
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: atom/Atom.h:241
void add_particles(RMF::FileHandle fh, const ParticlesTemp &hs)
def get_particles_within_zone
Utility to retrieve particles from a hierarchy within a zone around a set of ps.
Definition: tools.py:1516
Set of Python classes to create a multi-state, multi-resolution IMP hierarchy.
def __init__
Constructor.
Definition: tools.py:124
static Weight setup_particle(Model *m, ParticleIndex pi)
Definition: Weight.h:31
A decorator for a particle which has bonds.
Rotation3D get_random_rotation_3d(const Rotation3D &center, double distance)
Pick a rotation at random near the provided one.
def get_molecules
This function returns the parent molecule hierarchies of given objects.
Definition: tools.py:1573
std::string get_module_version()
Change color code to hexadecimal to rgb.
Definition: tools.py:1145
ParticlesTemp get_particles(Model *m, const ParticleIndexes &ps)
This class stores integers in ordered compact lists eg: [[1,2,3],[6,7,8]] the methods help splitting ...
Definition: tools.py:994
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:862
This class initializes the root node of the global IMP.atom.Hierarchy.
def get_residue_gaps_in_hierarchy
Return the residue index gaps and contiguous segments in the hierarchy.
Definition: tools.py:474
def add_imp_provenance
Tag the given particle as being created by the current version of IMP.
def add
index can be a integer or a list of integers
Definition: tools.py:1022
Vector3D get_random_vector_in(const Cylinder3D &c)
Generate a random vector in a cylinder with uniform density.
def get_all_leaves
Just get the leaves from a list of hierarchies.
Definition: tools.py:1476
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 color2rgb
Given a Chimera color name or hex color value, return RGB.
Definition: tools.py:1937
def input_adaptor
Adapt things for PMI (degrees of freedom, restraints, ...) Returns list of list of hierarchies...
Definition: tools.py:1275
def get_terminal_residue
Get the particle of the terminal residue at the GIVEN resolution (NOTE: not the closest resolution!)...
Definition: tools.py:421
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)
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:664
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:924
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:86
This class converts three to one letter codes, and return X for any unknown codes.
Definition: tools.py:1449
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:462
def __init__
index can be a integer or a list of integers
Definition: tools.py:1011
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
def get_flatten
Returns a flatten list.
Definition: tools.py:1074
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:119
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:390
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
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:309
def get_densities
Given a list of PMI objects, returns all density hierarchies within these objects.
Definition: tools.py:1616
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.
def set_coordinates_from_rmf
Extract frame from RMF file and fill coordinates.
Definition: tools.py:1253
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:545
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:337
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: core/Gaussian.h:47
def get_sorted_segments
Returns sequence-sorted segments array, each containing the first particle the last particle and the ...
Definition: tools.py:1401
def get_rbs_and_beads
Returns unique objects in original order.
Definition: tools.py:1550
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:882
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Definition: rigid_bodies.h:711
Warning for probably incorrect input parameters.
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.
Store objects in order they were added, but with default type.
Definition: tools.py:1223
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...
def sublist_iterator
Yield all sublists of length >= lmin and <= lmax.
Definition: tools.py:955
A particle with a color.
Definition: Colored.h:23