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