3 """@namespace IMP.pmi.tools
4 Miscellaneous utilities.
7 from __future__
import print_function, division
15 from math
import log,pi,sqrt,exp
20 from time
import process_time
22 from time
import clock
as process_time
25 from collections
import defaultdict
27 from collections
import OrderedDict
29 from IMP.pmi._compat_collections
import OrderedDict
32 def _get_system_for_hier(hier):
33 """Given a hierarchy, return the System that created it, or None"""
36 if hier
and not hasattr(hier,
'get_parent'):
43 if hasattr(hier,
'_pmi2_system'):
44 return hier._pmi2_system()
47 for ws
in IMP.pmi.topology.System._all_systems:
49 if s
and s.hier == hier:
52 hier = hier.get_parent()
54 def _all_protocol_outputs(hier):
55 """Iterate over all (ProtocolOutput, State) pairs for the given hierarchy"""
56 system = _get_system_for_hier(hier)
58 for state
in system.states:
59 for p
in state._protocol_output:
62 def _add_pmi_provenance(p):
63 """Tag the given particle as being created by the current version of PMI."""
67 location=
"https://integrativemodeling.org")
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)
77 def _add_restraint_sets(model, mk, mk_rmf):
80 model.add_data(mk, rs)
81 model.add_data(mk_rmf, rs_rmf)
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.
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))
98 rs, rs_rmf = _add_restraint_sets(model, mk, mk_rmf)
99 rs.add_restraint(restraint)
101 rs_rmf.add_restraint(restraint)
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)
113 return IMP.RestraintSet.get_from(model.get_data(mk_rmf))
115 return IMP.RestraintSet.get_from(model.get_data(mk))
118 """Collect timing information.
119 Add an instance of this class to outputobjects to get timing information
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()
128 self.isdelta = isdelta
130 def set_label(self, labelstr):
131 self.label = labelstr
133 def get_output(self):
136 newtime = process_time()
137 output[
"Stopwatch_" + self.label +
"_delta_seconds"] \
138 = str(newtime - self.starttime)
139 self.starttime = newtime
141 output[
"Stopwatch_" + self.label +
"_elapsed_seconds"] \
142 = str(process_time() - self.starttime)
146 class SetupNuisance(object):
148 def __init__(self, m, initialvalue, minvalue, maxvalue, isoptimized=True,
156 nuisance.set_lower(minvalue)
158 nuisance.set_upper(maxvalue)
161 nuisance.set_is_optimized(nuisance.get_nuisance_key(), isoptimized)
162 self.nuisance = nuisance
164 def get_particle(self):
168 class SetupWeight(object):
170 def __init__(self, m, isoptimized=True, nweights_or_weights=None):
172 if isinstance(nweights_or_weights, int):
174 pw, nweights_or_weights
178 nweights_or_weights = list(nweights_or_weights)
180 pw, nweights_or_weights
184 self.weight.set_weights_are_optimized(isoptimized)
186 def get_particle(self):
190 class SetupSurface(object):
192 def __init__(self, m, center, normal, isoptimized=True):
195 self.surface.set_coordinates_are_optimized(isoptimized)
196 self.surface.set_normal_is_optimized(isoptimized)
198 def get_particle(self):
202 class ParticleToSampleList(object):
206 self.dictionary_particle_type = {}
207 self.dictionary_particle_transformation = {}
208 self.dictionary_particle_name = {}
215 particle_transformation,
217 if not particle_type
in [
"Rigid_Bodies",
"Floppy_Bodies",
"Nuisances",
"X_coord",
"Weights",
"Surfaces"]:
218 raise TypeError(
"not the right particle type")
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
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
240 raise TypeError(
"ParticleToSampleList: not the right transformation format for Surfaces, should be a tuple of floats")
242 if isinstance(particle_transformation, float):
243 self.dictionary_particle_transformation[
244 particle] = particle_transformation
245 self.dictionary_particle_name[particle] = name
247 raise TypeError(
"ParticleToSampleList: not the right transformation format, should be a float")
249 def get_particles_to_sample(self):
251 for particle
in self.dictionary_particle_type:
252 key = self.dictionary_particle_type[
253 particle] +
"ParticleToSampleList_" + self.dictionary_particle_name[particle] +
"_" + self.label
256 self.dictionary_particle_transformation[particle])
261 def get_cross_link_data(directory, filename, dist, omega, sigma,
262 don=
None, doff=
None, prior=0, type_of_profile=
"gofr"):
264 (distmin, distmax, ndist) = dist
265 (omegamin, omegamax, nomega) = omega
266 (sigmamin, sigmamax, nsigma) = sigma
269 with open(filen)
as xlpot:
270 dictionary = ast.literal_eval(xlpot.readline())
272 xpot = dictionary[directory][filename][
"distance"]
273 pot = dictionary[directory][filename][type_of_profile]
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)
279 if not don
is None and not doff
is None:
298 def get_grid(gmin, gmax, ngrid, boundaries):
300 dx = (gmax - gmin) / float(ngrid)
301 for i
in range(0, ngrid + 1):
302 if(
not boundaries
and i == 0):
304 if(
not boundaries
and i == ngrid):
306 grid.append(gmin + float(i) * dx)
312 def get_log_grid(gmin, gmax, ngrid):
314 for i
in range(0, ngrid + 1):
315 grid.append(gmin * exp(float(i) / ngrid * log(gmax / gmin)))
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%"'
328 import pyparsing
as pp
330 operator = pp.Regex(
">=|<=|!=|>|<|==|in").setName(
"operator")
331 value = pp.QuotedString(
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)
338 expr = pp.operatorPrecedence(condition, [
339 (
"OR", 2, pp.opAssoc.LEFT, ),
340 (
"AND", 2, pp.opAssoc.LEFT, ),
343 parsedstring = str(expr.parseString(inputstring)) \
349 .replace(
"{",
"float(entry['") \
350 .replace(
"}",
"'])") \
351 .replace(
":",
"str(entry['") \
352 .replace(
"}",
"'])") \
353 .replace(
"AND",
"and") \
358 def open_file_or_inline_text(filename):
360 fl = open(filename,
"r")
362 fl = filename.split(
"\n")
365 def get_ids_from_fasta_file(fastafile):
367 with open(fastafile)
as ff:
376 this function works with plain hierarchies, as read from the pdb,
377 no multi-scale hierarchies
384 atom_type=IMP.atom.AT_CA)
392 print(
"get_closest_residue_position: exiting while loop without result")
394 p = sel.get_selected_particles()
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" % (
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])))
408 Return the residue index gaps and contiguous segments in the hierarchy.
410 @param hierarchy hierarchy to examine
411 @param start first residue index
412 @param end last residue index
414 @return A list of lists of the form
415 [[1,100,"cont"],[101,120,"gap"],[121,200,"cont"]]
418 for n, rindex
in enumerate(range(start, end + 1)):
420 atom_type=IMP.atom.AT_CA)
422 if len(sel.get_selected_particles()) == 0:
426 rindexcont = start - 1
427 if rindexgap == rindex - 1:
433 gaps.append([rindex, rindex,
"gap"])
439 rindexgap = start - 1
441 if rindexcont == rindex - 1:
448 gaps.append([rindex, rindex,
"cont"])
459 def set_map_element(self, xvalue, yvalue):
460 self.map[xvalue] = yvalue
462 def get_map_element(self, invalue):
463 if isinstance(invalue, float):
467 dist = (invalue - x) * (invalue - x)
476 return self.map[minx]
477 elif isinstance(invalue, str):
478 return self.map[invalue]
480 raise TypeError(
"wrong type for map")
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
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]
498 copynum = tuple_selection[3]
499 if copynum
is not None:
500 kwds[
'copy_index'] = copynum
504 statenum = tuple_selection[4]
505 if statenum
is not None:
506 kwds[
'state_index'] = statenum
513 residue_indexes=range(1,rbegin),
515 return s.get_selected_particles()
517 kwds[
'residue_indexes'] = range(rbegin,rend+1)
519 return s.get_selected_particles()
523 def get_db_from_csv(csvfilename):
526 with open(csvfilename)
as fh:
527 csvr = csv.DictReader(fh)
534 '''Return the component name provided a particle and a list of names'''
536 protname = root.get_name()
538 while not protname
in list_of_names:
539 root0 = root.get_parent()
542 protname = root0.get_name()
547 if "Beads" in protname:
550 return (protname, is_a_bead)
555 Retrieve the residue indexes for the given particle.
557 The particle must be an instance of Fragment,Residue, Atom or Molecule
558 or else returns an empty list
569 resind_tmp=IMP.pmi.tools.OrderedSet()
575 resind=list(resind_tmp)
581 def sort_by_residues(particles):
584 sorted_particles_residues = sorted(
586 key=
lambda tup: tup[1])
587 particles = [p[0]
for p
in sorted_particles_residues]
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
603 comm.send(data, dest=0, tag=11)
606 for i
in range(1, number_of_processes):
607 data_tmp = comm.recv(source=i, tag=11)
608 if isinstance(data, list):
610 elif isinstance(data, dict):
611 data.update(data_tmp)
613 raise TypeError(
"data not supported, use list or dictionaries")
615 for i
in range(1, number_of_processes):
616 comm.send(data, dest=i, tag=11)
619 data = comm.recv(source=0, tag=11)
628 Yield all sublists of length >= lmin and <= lmax
636 for j
in range(i + 1, n):
637 if len(l[i:j]) <= lmax
and len(l[i:j]) >= lmin:
642 return [item
for sublist
in l
for item
in sublist]
646 """ Yield successive length-sized chunks from a list.
648 for i
in range(0, len(list), length):
649 yield list[i:i + length]
652 def chunk_list_into_segments(seq, num):
654 avg = len(seq) / float(num)
658 while last < len(seq):
659 out.append(seq[int(last):int(last + avg)])
667 ''' This class stores integers
668 in ordered compact lists eg:
670 the methods help splitting and merging the internal lists
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)
683 '''index can be a integer or a list of integers '''
684 if isinstance(index, int):
686 elif isinstance(index, list):
687 self.segs=[[index[0]]]
691 raise TypeError(
"index must be an int or list of ints")
694 '''index can be a integer or a list of integers '''
695 if isinstance(index, int):
698 for n,s
in enumerate(self.segs):
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]
716 for n
in range(len(self.segs)):
719 self.segs.sort(key=
lambda tup: tup[0])
721 elif isinstance(index, list):
725 raise TypeError(
"index must be an int or list of ints")
728 '''index can be a integer'''
729 for n,s
in enumerate(self.segs):
736 i=self.segs[n].index(index)
738 self.segs.append(s[i+1:])
739 for n
in range(len(self.segs)):
741 if len(self.segs[n])==0:
743 self.segs.sort(key=
lambda tup: tup[0])
746 ''' Returns a flatten list '''
747 return [item
for sublist
in self.segs
for item
in sublist]
751 for seg
in self.segs:
752 ret_tmp+=str(seg[0])+
"-"+str(seg[-1])+
","
760 def normal_density_function(expected_value, sigma, x):
762 1 / math.sqrt(2 * math.pi) / sigma *
763 math.exp(-(x - expected_value) ** 2 / 2 / sigma / sigma)
767 def log_normal_density_function(expected_value, sigma, x):
769 1 / math.sqrt(2 * math.pi) / sigma / x *
770 math.exp(-(math.log(x / expected_value) ** 2 / 2 / sigma / sigma))
774 def print_multicolumn(list_of_strings, ncolumns=2, truncate=40):
780 for i
in range(len(l) % cols):
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))
788 '''Change color code to hexadecimal to rgb'''
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'
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]])
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)
803 class OrderedSet(collections.MutableSet):
805 def __init__(self, iterable=None):
807 end += [
None, end, end]
809 if iterable
is not None:
815 def __contains__(self, key):
816 return key
in self.map
819 if key
not in self.map:
822 curr[2] = end[1] = self.map[key] = [key, curr, end]
824 def discard(self, key):
826 key, prev, next = self.map.pop(key)
833 while curr
is not end:
837 def __reversed__(self):
840 while curr
is not end:
844 def pop(self, last=True):
846 raise KeyError(
'set is empty')
856 return '%s()' % (self.__class__.__name__,)
857 return '%s(%r)' % (self.__class__.__name__, list(self))
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)
866 """Store objects in order they were added, but with default type.
867 Source: http://stackoverflow.com/a/4127426/2608793
869 def __init__(self, *args, **kwargs):
871 self.default_factory =
None
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]
877 super(OrderedDefaultDict, self).__init__(*args, **kwargs)
879 def __missing__ (self, key):
880 if self.default_factory
is None:
882 self[key] = default = self.default_factory()
885 def __reduce__(self):
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()
890 return self.__class__, args,
None,
None, self.iteritems()
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
901 rh = RMF.open_rmf_file_read_only(rmf_fn)
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.
934 if hasattr(stuff,
'__iter__'):
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')
952 except (NotImplementedError, TypeError):
963 if is_system
or is_state
or is_molecule
or is_temp_residue:
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()]
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()]
981 for molecule
in stuff:
982 indexes_per_mol[molecule] += [r.get_index()
for r
in molecule.get_residues()]
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':
991 residue_indexes=indexes_per_mol[mol])
994 resolution=pmi_resolution,
995 residue_indexes=indexes_per_mol[mol])
996 ps = sel.get_selected_particles()
999 if warn_about_slices:
1000 rset = set(indexes_per_mol[mol])
1004 if not fset <= rset:
1010 resbreak = maxf
if minf==minset
else minset-1
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,
1025 if pmi_resolution==
'all':
1033 hier_list = [hier_list]
1035 raise Exception(
'input_adaptor: you passed something of wrong type or a list with mixed types')
1037 if flatten
and pmi_input:
1038 return [h
for sublist
in hier_list
for h
in sublist]
1044 """Returns sequence-sorted segments array, each containing the first particle
1045 the last particle and the first residue index."""
1047 from operator
import itemgetter
1050 raise ValueError(
"only pass stuff from one Molecule, please")
1066 segs.append((start, end, startres))
1067 return sorted(segs, key=itemgetter(2))
1070 """Decorate the sequence-consecutive particles from a PMI2 molecule with a bond,
1071 so that they appear connected in the rmf file"""
1073 for x
in range(len(SortedSegments) - 1):
1075 last = SortedSegments[x][1]
1076 first = SortedSegments[x + 1][0]
1078 p1 = last.get_particle()
1079 p2 = first.get_particle()
1093 """This class converts three to one letter codes, and return X for
1094 any unknown codes"""
1095 def __init__(self, is_nucleic=False):
1097 threetoone = IMP.pmi.alphabets.rna.charmm_to_one
1099 threetoone = IMP.pmi.alphabets.amino_acid.charmm_to_one
1100 defaultdict.__init__(self,
lambda:
"X", threetoone)
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)
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))
1118 """Perform selection using the usual keywords but return ALL resolutions (BEADS and GAUSSIANS).
1119 Returns in flat list!
1124 if hier
is not None:
1127 warnings.warn(
"You passed nothing to select_at_all_resolutions()",
1135 raise Exception(
'select_at_all_resolutions: you have to pass an IMP Hierarchy')
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")
1141 representation_type=IMP.atom.BALLS,**kwargs)
1143 representation_type=IMP.atom.DENSITIES,**kwargs)
1144 ret |= OrderedSet(selB.get_selected_particles())
1145 ret |= OrderedSet(selD.get_selected_particles())
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
1164 backbone_types=[
'C',
'N',
'CB',
'O']
1165 if exclude_backbone:
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()
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]
1179 zone_ps = [h.get_particle()
for h
in final_ps]
1184 """Returns unique objects in original order"""
1188 if not hasattr(hiers,
'__iter__'):
1195 rbs_ordered.append(rb)
1200 rbs_ordered.append(rb)
1204 return rbs_ordered,beads
1207 "This function returns the parent molecule hierarchies of given objects"
1208 stuff=
input_adaptor(input_objects, pmi_resolution=
'all',flatten=
True)
1213 while not (is_molecule
or is_root):
1214 root=IMP.atom.get_root(h)
1221 return list(molecules)
1223 def get_molecules_dictionary(input_objects):
1224 moldict=defaultdict(list)
1227 moldict[name].append(mol)
1233 def get_molecules_dictionary_by_copy(input_objects):
1234 moldict=defaultdict(dict)
1238 moldict[name][c]=mol
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():
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).
1260 densities+=
IMP.atom.Selection(i,representation_type=IMP.atom.DENSITIES).get_selected_particles()
1264 max_translation=300., max_rotation=2.0 * pi,
1265 avoidcollision_rb=
True, avoidcollision_fb=
False,
1266 cutoff=10.0, niterations=100,
1268 excluded_rigid_bodies=[],
1269 hierarchies_excluded_from_collision=[],
1270 hierarchies_included_in_collision=[],
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>
1301 pmi_resolution=
'all',
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()
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")
1315 gcpf.set_distance(cutoff)
1319 pmi_resolution=
'all',
1323 pmi_resolution=
'all',
1326 collision_excluded_idxs = set([l.get_particle().
get_index()
for h
in collision_excluded_hierarchies \
1329 collision_included_idxs = set([l.get_particle().
get_index()
for h
in collision_included_hierarchies \
1336 all_idxs.append(p.get_particle_index())
1338 collision_excluded_idxs.add(p.get_particle_index())
1340 if bounding_box
is not None:
1341 ((x1, y1, z1), (x2, y2, z2)) = bounding_box
1346 all_idxs = set(all_idxs) | collision_included_idxs
1347 all_idxs = all_idxs - collision_excluded_idxs
1349 print(
'shuffling', len(rigid_bodies),
'rigid bodies')
1350 for rb
in rigid_bodies:
1351 if rb
not in excluded_rigid_bodies:
1353 if avoidcollision_rb:
1354 rb_idxs = set(rb.get_member_particle_indexes()) - \
1355 collision_excluded_idxs
1356 other_idxs = all_idxs - rb_idxs
1362 while niter < niterations:
1363 rbxyz = (rb.get_x(), rb.get_y(), rb.get_z())
1383 debug.append([rb, other_idxs
if avoidcollision_rb
else set()])
1387 if avoidcollision_rb:
1389 npairs = len(gcpf.get_close_pairs(mdl,
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")
1405 print(
'shuffling', len(flexible_beads),
'flexible beads')
1406 for fb
in flexible_beads:
1408 if avoidcollision_fb:
1410 other_idxs = all_idxs - fb_idxs
1416 while niter < niterations:
1429 xyz = memb.get_internal_coordinates()
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))
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()])
1448 debug.append([d,other_idxs
if avoidcollision_fb
else set()])
1454 if avoidcollision_fb:
1456 npairs = len(gcpf.get_close_pairs(mdl,
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")
1472 class ColorHierarchy(object):
1474 def __init__(self,hier):
1475 import matplotlib
as mpl
1477 import matplotlib.pyplot
as plt
1481 hier.ColorHierarchy=self
1485 self.method=self.nochange
1493 def get_color(self,fl):
1496 def get_log_scale(self,fl):
1499 return math.log(fl+eps)
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:
1510 c=self.get_color(float(ri)/self.last)
1514 avr=sum(ris)/len(ris)
1515 c=self.get_color(float(avr)/self.last)
1518 def color_by_uncertainty(self):
1519 self.method=self.color_by_uncertainty
1520 self.scheme=self.mpl.cm.jet
1527 self.first=self.get_log_scale(1.0)
1528 self.last=self.get_log_scale(100.0)
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))
1536 def get_color_bar(self,filename):
1537 import matplotlib
as mpl
1539 import matplotlib.pyplot
as plt
1542 fig = plt.figure(figsize=(8, 3))
1543 ax1 = fig.add_axes([0.05, 0.80, 0.9, 0.15])
1546 norm = mpl.colors.Normalize(vmin=0.0, vmax=1.0)
1548 if self.method == self.color_by_uncertainty:
1549 angticks=[1.0,2.5,5.0,10.0,25.0,50.0,100.0]
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,
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)
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))
static bool get_is_setup(const IMP::ParticleAdaptor &p)
static bool get_is_setup(const IMP::ParticleAdaptor &p)
A decorator to associate a particle with a part of a protein/DNA/RNA.
Extends the functionality of IMP.atom.Molecule.
def add_script_provenance
Tag the given particle with the current Python script.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Set of Python classes to create a multi-state, multi-resolution IMP hierarchy.
static Weight setup_particle(Model *m, ParticleIndex pi)
Set up an empty Weight.
A decorator for a particle which has bonds.
Rotation3D get_random_rotation_3d(const Rotation3D ¢er, double distance)
Pick a rotation at random near the provided one.
An exception for an invalid usage of IMP.
def deprecated_function
Python decorator to mark a function as deprecated.
std::string get_module_version()
Return the version of this module, as a string.
static Surface setup_particle(Model *m, ParticleIndex pi)
Add uncertainty to a particle.
void add_particle(RMF::FileHandle fh, Particle *hs)
static bool get_is_setup(const IMP::ParticleAdaptor &p)
static bool get_is_setup(const IMP::ParticleAdaptor &p)
GenericHierarchies get_leaves(Hierarchy mhd)
Get all the leaves of the bit of hierarchy.
This class initializes the root node of the global IMP.atom.Hierarchy.
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.
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.
Stores a named protein chain.
static bool get_is_setup(Model *m, ParticleIndex pi)
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.
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)
A decorator for a particle representing an atom.
void transform(XYZ a, const algebra::Transformation3D &tr)
Apply a transformation to the particle.
static Bonded setup_particle(Model *m, ParticleIndex pi)
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.
static Scale setup_particle(Model *m, ParticleIndex pi)
A decorator for a particle that is part of a rigid body but not rigid.
Find all nearby pairs by testing all pairs.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
def deprecated_object
Python decorator to mark a class as deprecated.
A decorator for a residue.
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.
Rotation3D get_identity_rotation_3d()
Return a rotation that does not do anything.
void link_hierarchies(RMF::FileConstHandle fh, const atom::Hierarchies &hs)
Class to handle individual particles of a Model object.
Bond get_bond(Bonded a, Bonded b)
Get the bond between two particles.
Stores a list of Molecules all with the same State index.
std::string get_data_path(std::string file_name)
Return the full path to one of this module's data files.
int get_copy_index(Hierarchy h)
Walk up the hierarchy to find the current copy index.
Python classes to represent, score, sample and analyze models.
static bool get_is_setup(Model *m, ParticleIndex pi)
A decorator for a rigid body.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
Hierarchies get_leaves(const Selection &h)
A decorator for a molecule.
Select hierarchy particles identified by the biological name.
Support for the RMF file format for storing hierarchical molecular data and markup.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
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.
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...