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