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