3 """@namespace IMP.pmi.tools
4 Miscellaneous utilities.
7 from __future__
import print_function, division
14 from collections.abc
import MutableSet
16 from collections
import MutableSet
22 from time
import process_time
24 from time
import clock
as process_time
27 from collections
import defaultdict, 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()
55 def _all_protocol_outputs(hier):
56 """Iterate over all (ProtocolOutput, State) pairs for the
58 system = _get_system_for_hier(hier)
60 for state
in system.states:
61 for p
in state._protocol_output:
65 def _add_pmi_provenance(p):
66 """Tag the given particle as being created by the current version
71 location=
"https://integrativemodeling.org")
75 def _get_restraint_set_keys():
76 if not hasattr(_get_restraint_set_keys,
'pmi_rs_key'):
77 _get_restraint_set_keys.pmi_rs_key =
IMP.ModelKey(
"PMI restraints")
78 _get_restraint_set_keys.rmf_rs_key =
IMP.ModelKey(
"RMF restraints")
79 return (_get_restraint_set_keys.pmi_rs_key,
80 _get_restraint_set_keys.rmf_rs_key)
83 def _add_restraint_sets(model, mk, mk_rmf):
86 model.add_data(mk, rs)
87 model.add_data(mk_rmf, rs_rmf)
92 """Add a PMI restraint to the model.
93 Since Model.add_restraint() no longer exists (in modern IMP restraints
94 should be added to a ScoringFunction instead) store them instead in
95 a RestraintSet, and keep a reference to it in the Model.
97 If `add_to_rmf` is True, also add the restraint to a separate list
98 of restraints that will be written out to RMF files (by default, most
99 PMI restraints are not)."""
100 mk, mk_rmf = _get_restraint_set_keys()
101 if model.get_has_data(mk):
102 rs = IMP.RestraintSet.get_from(model.get_data(mk))
103 rs_rmf = IMP.RestraintSet.get_from(model.get_data(mk_rmf))
105 rs, rs_rmf = _add_restraint_sets(model, mk, mk_rmf)
106 rs.add_restraint(restraint)
108 rs_rmf.add_restraint(restraint)
112 """Get a RestraintSet containing all PMI restraints added to the model.
113 If `rmf` is True, return only the subset of these restraints that
114 should be written out to RMF files."""
115 mk, mk_rmf = _get_restraint_set_keys()
116 if not model.get_has_data(mk):
117 warnings.warn(
"no restraints added to model yet",
119 _add_restraint_sets(model, mk, mk_rmf)
121 return IMP.RestraintSet.get_from(model.get_data(mk_rmf))
123 return IMP.RestraintSet.get_from(model.get_data(mk))
127 """Collect timing information.
128 Add an instance of this class to outputobjects to get timing information
133 @param isdelta if True (the default) then report the time since the
134 last use of this class; if False, report cumulative time."""
135 self.starttime = process_time()
137 self.isdelta = isdelta
139 def set_label(self, labelstr):
140 self.label = labelstr
142 def get_output(self):
145 newtime = process_time()
146 output[
"Stopwatch_" + self.label +
"_delta_seconds"] \
147 = str(newtime - self.starttime)
148 self.starttime = newtime
150 output[
"Stopwatch_" + self.label +
"_elapsed_seconds"] \
151 = str(process_time() - self.starttime)
155 class SetupNuisance(object):
157 def __init__(self, m, initialvalue, minvalue, maxvalue, isoptimized=True,
165 nuisance.set_lower(minvalue)
167 nuisance.set_upper(maxvalue)
170 nuisance.set_is_optimized(nuisance.get_nuisance_key(), isoptimized)
171 self.nuisance = nuisance
173 def get_particle(self):
177 class SetupWeight(object):
179 def __init__(self, m, isoptimized=True, nweights_or_weights=None):
181 if isinstance(nweights_or_weights, int):
183 pw, nweights_or_weights
187 nweights_or_weights = list(nweights_or_weights)
189 pw, nweights_or_weights
193 self.weight.set_weights_are_optimized(isoptimized)
195 def get_particle(self):
199 class SetupSurface(object):
201 def __init__(self, m, center, normal, isoptimized=True):
204 self.surface.set_coordinates_are_optimized(isoptimized)
205 self.surface.set_normal_is_optimized(isoptimized)
207 def get_particle(self):
211 def get_cross_link_data(directory, filename, dist, omega, sigma,
212 don=
None, doff=
None, prior=0, type_of_profile=
"gofr"):
214 (distmin, distmax, ndist) = dist
215 (omegamin, omegamax, nomega) = omega
216 (sigmamin, sigmamax, nsigma) = sigma
219 with open(filen)
as xlpot:
220 dictionary = ast.literal_eval(xlpot.readline())
222 xpot = dictionary[directory][filename][
"distance"]
223 pot = dictionary[directory][filename][type_of_profile]
225 dist_grid =
get_grid(distmin, distmax, ndist,
False)
226 omega_grid = get_log_grid(omegamin, omegamax, nomega)
227 sigma_grid = get_log_grid(sigmamin, sigmamax, nsigma)
229 if don
is not None and doff
is not None:
249 def get_grid(gmin, gmax, ngrid, boundaries):
251 dx = (gmax - gmin) / float(ngrid)
252 for i
in range(0, ngrid + 1):
253 if not boundaries
and i == 0:
255 if not boundaries
and i == ngrid:
257 grid.append(gmin + float(i) * dx)
261 def get_log_grid(gmin, gmax, ngrid):
263 for i
in range(0, ngrid + 1):
264 grid.append(gmin * math.exp(float(i) / ngrid * math.log(gmax / gmin)))
270 example '"{ID_Score}" > 28 AND "{Sample}" ==
271 "%10_1%" OR ":Sample}" == "%10_2%" OR ":Sample}"
272 == "%10_3%" OR ":Sample}" == "%8_1%" OR ":Sample}" == "%8_2%"'
275 import pyparsing
as pp
277 operator = pp.Regex(
">=|<=|!=|>|<|==|in").setName(
"operator")
278 value = pp.QuotedString(
280 r"[+-]?\d+(:?\.\d*)?(:?[eE][+-]?\d+)?")
281 identifier = pp.Word(pp.alphas, pp.alphanums +
"_")
282 comparison_term = identifier | value
283 condition = pp.Group(comparison_term + operator + comparison_term)
285 expr = pp.operatorPrecedence(condition, [
286 (
"OR", 2, pp.opAssoc.LEFT, ),
287 (
"AND", 2, pp.opAssoc.LEFT, ),
290 parsedstring = str(expr.parseString(inputstring)) \
296 .replace(
"{",
"float(entry['") \
297 .replace(
"}",
"'])") \
298 .replace(
":",
"str(entry['") \
299 .replace(
"}",
"'])") \
300 .replace(
"AND",
"and") \
305 def open_file_or_inline_text(filename):
307 fl = open(filename,
"r")
309 fl = filename.split(
"\n")
313 def get_ids_from_fasta_file(fastafile):
315 with open(fastafile)
as ff:
318 ids.append(line[1:-1])
324 this function works with plain hierarchies, as read from the pdb,
325 no multi-scale hierarchies
332 atom_type=IMP.atom.AT_CA)
340 print(
"get_closest_residue_position: exiting while loop "
343 p = sel.get_selected_particles()
348 print(
"get_closest_residue_position: got NO residues for hierarchy "
349 "%s and residue %i" % (hier, resindex))
351 "get_closest_residue_position: got NO residues for hierarchy "
352 "%s and residue %i" % (hier, resindex))
355 "got multiple residues for hierarchy %s and residue %i; the list "
357 % (hier, resindex, str([pp.get_name()
for pp
in p])))
362 Return the residue index gaps and contiguous segments in the hierarchy.
364 @param hierarchy hierarchy to examine
365 @param start first residue index
366 @param end last residue index
368 @return A list of lists of the form
369 [[1,100,"cont"],[101,120,"gap"],[121,200,"cont"]]
372 for n, rindex
in enumerate(range(start, end + 1)):
374 atom_type=IMP.atom.AT_CA)
376 if len(sel.get_selected_particles()) == 0:
380 rindexcont = start - 1
381 if rindexgap == rindex - 1:
387 gaps.append([rindex, rindex,
"gap"])
393 rindexgap = start - 1
395 if rindexcont == rindex - 1:
402 gaps.append([rindex, rindex,
"cont"])
413 def set_map_element(self, xvalue, yvalue):
414 self.map[xvalue] = yvalue
416 def get_map_element(self, invalue):
417 if isinstance(invalue, float):
421 dist = (invalue - x) * (invalue - x)
430 return self.map[minx]
431 elif isinstance(invalue, str):
432 return self.map[invalue]
434 raise TypeError(
"wrong type for map")
438 """New tuple format: molname OR (start,stop,molname,copynum,statenum)
439 Copy and state are optional. Can also use 'None' for them which will
440 get all. You can also pass -1 for stop which will go to the end.
441 Returns the particles
444 kwds[
'resolution'] = resolution
445 if isinstance(tuple_selection, str):
446 kwds[
'molecule'] = tuple_selection
447 elif isinstance(tuple_selection, tuple):
448 rbegin = tuple_selection[0]
449 rend = tuple_selection[1]
450 kwds[
'molecule'] = tuple_selection[2]
452 copynum = tuple_selection[3]
453 if copynum
is not None:
454 kwds[
'copy_index'] = copynum
458 statenum = tuple_selection[4]
459 if statenum
is not None:
460 kwds[
'state_index'] = statenum
467 residue_indexes=range(1, rbegin),
469 return s.get_selected_particles()
471 kwds[
'residue_indexes'] = range(rbegin, rend+1)
473 return s.get_selected_particles()
476 def get_db_from_csv(csvfilename, encoding=None):
477 if sys.version_info[0] == 2:
478 def open_with_encoding(fname, encoding):
481 open_with_encoding = open
484 with open_with_encoding(csvfilename, encoding=encoding)
as fh:
485 csvr = csv.DictReader(fh)
487 outputlist.append(ls)
492 '''Return the component name provided a particle and a list of names'''
494 protname = root.get_name()
496 while protname
not in list_of_names:
497 root0 = root.get_parent()
500 protname = root0.get_name()
505 if "Beads" in protname:
508 return (protname, is_a_bead)
513 Retrieve the residue indexes for the given particle.
515 The particle must be an instance of Fragment,Residue, Atom or Molecule
516 or else returns an empty list
527 resind_tmp = IMP.pmi.tools.OrderedSet()
534 resind = list(resind_tmp)
540 def sort_by_residues(particles):
543 sorted_particles_residues = sorted(
545 key=
lambda tup: tup[1])
546 particles = [p[0]
for p
in sorted_particles_residues]
555 """Synchronize data over a parallel run"""
556 from mpi4py
import MPI
557 comm = MPI.COMM_WORLD
558 rank = comm.Get_rank()
559 number_of_processes = comm.size
562 comm.send(data, dest=0, tag=11)
565 for i
in range(1, number_of_processes):
566 data_tmp = comm.recv(source=i, tag=11)
567 if isinstance(data, list):
569 elif isinstance(data, dict):
570 data.update(data_tmp)
572 raise TypeError(
"data not supported, use list or dictionaries")
574 for i
in range(1, number_of_processes):
575 comm.send(data, dest=i, tag=11)
578 data = comm.recv(source=0, tag=11)
588 Yield all sublists of length >= lmin and <= lmax
594 for j
in range(i + lmin, min(n + 1, i + 1 + lmax)):
598 def flatten_list(ls):
599 return [item
for sublist
in ls
for item
in sublist]
603 """ Yield successive length-sized chunks from a list.
605 for i
in range(0, len(list), length):
606 yield list[i:i + length]
609 def chunk_list_into_segments(seq, num):
611 avg = len(seq) / float(num)
615 while last < len(seq):
616 out.append(seq[int(last):int(last + avg)])
624 ''' This class stores integers
625 in ordered compact lists eg:
627 the methods help splitting and merging the internal lists
629 s=Segments([1,2,3]) is [[1,2,3]]
630 s.add(4) is [[1,2,3,4]] (add right)
631 s.add(3) is [[1,2,3,4]] (item already existing)
632 s.add(7) is [[1,2,3,4],[7]] (new list)
633 s.add([8,9]) is [[1,2,3,4],[7,8,9]] (add item right)
634 s.add([5,6]) is [[1,2,3,4,5,6,7,8,9]] (merge)
635 s.remove(3) is [[1,2],[4,5,6,7,8,9]] (split)
640 '''index can be a integer or a list of integers '''
641 if isinstance(index, int):
642 self.segs = [[index]]
643 elif isinstance(index, list):
644 self.segs = [[index[0]]]
648 raise TypeError(
"index must be an int or list of ints")
651 '''index can be a integer or a list of integers '''
652 if isinstance(index, (int, numpy.int32, numpy.int64)):
655 for n, s
in enumerate(self.segs):
663 if mergeright
is None and mergeleft
is None:
664 self.segs.append([index])
665 if mergeright
is not None and mergeleft
is None:
666 self.segs[mergeright].append(index)
667 if mergeleft
is not None and mergeright
is None:
668 self.segs[mergeleft] = [index]+self.segs[mergeleft]
669 if mergeleft
is not None and mergeright
is not None:
670 self.segs[mergeright] = \
671 self.segs[mergeright]+[index]+self.segs[mergeleft]
672 del self.segs[mergeleft]
674 for n
in range(len(self.segs)):
677 self.segs.sort(key=
lambda tup: tup[0])
679 elif isinstance(index, list):
683 raise TypeError(
"index must be an int or list of ints")
686 '''index can be a integer'''
687 for n, s
in enumerate(self.segs):
692 self.segs[n] = s[:-1]
694 i = self.segs[n].index(index)
696 self.segs.append(s[i+1:])
697 for n
in range(len(self.segs)):
699 if len(self.segs[n]) == 0:
701 self.segs.sort(key=
lambda tup: tup[0])
704 ''' Returns a flatten list '''
705 return [item
for sublist
in self.segs
for item
in sublist]
709 for seg
in self.segs:
710 ret_tmp += str(seg[0])+
"-"+str(seg[-1])+
","
711 ret = ret_tmp[:-1]+
"]"
719 def normal_density_function(expected_value, sigma, x):
721 1 / math.sqrt(2 * math.pi) / sigma *
722 math.exp(-(x - expected_value) ** 2 / 2 / sigma / sigma)
726 def log_normal_density_function(expected_value, sigma, x):
728 1 / math.sqrt(2 * math.pi) / sigma / x *
729 math.exp(-(math.log(x / expected_value) ** 2 / 2 / sigma / sigma))
733 def print_multicolumn(list_of_strings, ncolumns=2, truncate=40):
739 for i
in range(len(ls) % cols):
742 split = [ls[i:i + len(ls) // cols]
743 for i
in range(0, len(ls), len(ls) // cols)]
744 for row
in zip(*split):
745 print(
"".join(str.ljust(i, truncate)
for i
in row))
749 '''Change color code to hexadecimal to rgb'''
751 self._NUMERALS =
'0123456789abcdefABCDEF'
752 self._HEXDEC = dict((v, int(v, 16))
for v
in
753 (x+y
for x
in self._NUMERALS
754 for y
in self._NUMERALS))
755 self.LOWERCASE, self.UPPERCASE =
'x',
'X'
757 def rgb(self, triplet):
758 return (float(self._HEXDEC[triplet[0:2]]),
759 float(self._HEXDEC[triplet[2:4]]),
760 float(self._HEXDEC[triplet[4:6]]))
762 def triplet(self, rgb, lettercase=None):
763 if lettercase
is None:
764 lettercase = self.LOWERCASE
765 return format(rgb[0] << 16 | rgb[1] << 8 | rgb[2],
'06'+lettercase)
769 class OrderedSet(MutableSet):
771 def __init__(self, iterable=None):
773 end += [
None, end, end]
775 if iterable
is not None:
781 def __contains__(self, key):
782 return key
in self.map
785 if key
not in self.map:
788 curr[2] = end[1] = self.map[key] = [key, curr, end]
790 def discard(self, key):
792 key, prev, next = self.map.pop(key)
799 while curr
is not end:
803 def __reversed__(self):
806 while curr
is not end:
810 def pop(self, last=True):
812 raise KeyError(
'set is empty')
822 return '%s()' % (self.__class__.__name__,)
823 return '%s(%r)' % (self.__class__.__name__, list(self))
825 def __eq__(self, other):
826 if isinstance(other, OrderedSet):
827 return len(self) == len(other)
and list(self) == list(other)
828 return set(self) == set(other)
832 """Store objects in order they were added, but with default type.
833 Source: http://stackoverflow.com/a/4127426/2608793
835 def __init__(self, *args, **kwargs):
837 self.default_factory =
None
839 if not (args[0]
is None or callable(args[0])):
840 raise TypeError(
'first argument must be callable or None')
841 self.default_factory = args[0]
843 super(OrderedDefaultDict, self).__init__(*args, **kwargs)
845 def __missing__(self, key):
846 if self.default_factory
is None:
848 self[key] = default = self.default_factory()
851 def __reduce__(self):
852 args = (self.default_factory,)
if self.default_factory
else ()
853 if sys.version_info[0] >= 3:
854 return self.__class__, args,
None,
None, self.items()
856 return self.__class__, args,
None,
None, self.iteritems()
862 """Extract frame from RMF file and fill coordinates. Must be identical
865 @param hier The (System) hierarchy to fill (e.g. after you've built it)
866 @param rmf_fn The file to extract from
867 @param frame_num The frame number to extract
869 rh = RMF.open_rmf_file_read_only(rmf_fn)
875 def input_adaptor(stuff, pmi_resolution=0, flatten=False, selection_tuple=None,
876 warn_about_slices=
True):
877 """Adapt things for PMI (degrees of freedom, restraints, ...)
878 Returns list of list of hierarchies, separated into Molecules if possible.
879 The input can be a list, or a list of lists (iterable of ^1 or
881 (iterable of ^2) Hierarchy -> returns input as list of list of hierarchies,
882 only one entry, not grouped by molecules.
883 (iterable of ^2) PMI::System/State/Molecule/TempResidue ->
884 returns residue hierarchies, grouped in molecules, at requested
887 @param stuff Can be one of the following inputs:
888 IMP Hierarchy, PMI System/State/Molecule/TempResidue, or a
889 list/set (of list/set) of them.
890 Must be uniform input, however. No mixing object types.
891 @param pmi_resolution For selecting, only does it if you pass PMI
892 objects. Set it to "all" if you want all resolutions!
893 @param flatten Set to True if you just want all hierarchies in one list.
894 @param warn_about_slices Print a warning if you are requesting only part
895 of a bead. Sometimes you just don't care!
896 @note since this relies on IMP::atom::Selection, this will not return
897 any objects if they weren't built! But there should be no problem
898 if you request unbuilt residues - they should be ignored.
904 if hasattr(stuff,
'__iter__'):
907 thelist = list(stuff)
910 if all(hasattr(el,
'__iter__')
for el
in thelist):
911 thelist = [i
for sublist
in thelist
for i
in sublist]
912 elif any(hasattr(el,
'__iter__')
for el
in thelist):
913 raise Exception(
'input_adaptor: input_object must be a list '
914 'or a list of lists')
923 except (NotImplementedError, TypeError):
935 if is_system
or is_state
or is_molecule
or is_temp_residue:
942 for state
in system.get_states():
943 mdict = state.get_molecules()
944 for molname
in mdict:
945 for copy
in mdict[molname]:
946 indexes_per_mol[copy] += \
947 [r.get_index()
for r
in copy.get_residues()]
950 mdict = state.get_molecules()
951 for molname
in mdict:
952 for copy
in mdict[molname]:
953 indexes_per_mol[copy] += [r.get_index()
954 for r
in copy.get_residues()]
956 for molecule
in stuff:
957 indexes_per_mol[molecule] += [r.get_index()
958 for r
in molecule.get_residues()]
960 for tempres
in stuff:
961 indexes_per_mol[tempres.get_molecule()].append(
963 for mol
in indexes_per_mol:
964 if pmi_resolution ==
'all':
968 mol.get_hierarchy(), residue_indexes=indexes_per_mol[mol])
971 resolution=pmi_resolution,
972 residue_indexes=indexes_per_mol[mol])
973 ps = sel.get_selected_particles()
976 if warn_about_slices:
977 rset = set(indexes_per_mol[mol])
987 resbreak = maxf
if minf == minset
else minset-1
989 'You are trying to select only part of the '
990 'bead %s:%i-%i. The residues you requested '
991 'are %i-%i. You can fix this by: '
992 '1) requesting the whole bead/none of it; or'
993 '2) break the bead up by passing '
994 'bead_extra_breaks=[\'%i\'] in '
995 'molecule.add_representation()'
996 % (mol.get_name(), minset, maxset, minf, maxf,
1002 if pmi_resolution ==
'all':
1008 h, resolution=pmi_resolution).get_selected_particles()
1011 hier_list = [hier_list]
1013 raise Exception(
'input_adaptor: you passed something of wrong type '
1014 'or a list with mixed types')
1016 if flatten
and pmi_input:
1017 return [h
for sublist
in hier_list
for h
in sublist]
1023 """Returns sequence-sorted segments array, each containing the first
1024 particle the last particle and the first residue index."""
1026 from operator
import itemgetter
1029 raise ValueError(
"only pass stuff from one Molecule, please")
1044 segs.append((start, end, startres))
1045 return sorted(segs, key=itemgetter(2))
1049 """Decorate the sequence-consecutive particles from a PMI2 molecule
1050 with a bond, so that they appear connected in the rmf file"""
1052 for x
in range(len(SortedSegments) - 1):
1054 last = SortedSegments[x][1]
1055 first = SortedSegments[x + 1][0]
1057 p1 = last.get_particle()
1058 p2 = first.get_particle()
1071 """ Just get the leaves from a list of hierarchies """
1072 lvs = list(itertools.chain.from_iterable(
1078 """Perform selection using the usual keywords but return ALL
1079 resolutions (BEADS and GAUSSIANS).
1080 Returns in flat list!
1085 if hier
is not None:
1088 warnings.warn(
"You passed nothing to select_at_all_resolutions()",
1096 raise Exception(
'select_at_all_resolutions: you have to pass '
1099 raise Exception(
'select_at_all_resolutions: you have to pass '
1101 if 'resolution' in kwargs
or 'representation_type' in kwargs:
1102 raise Exception(
"don't pass resolution or representation_type "
1105 representation_type=IMP.atom.BALLS,
1108 representation_type=IMP.atom.DENSITIES,
1110 ret |= OrderedSet(selB.get_selected_particles())
1111 ret |= OrderedSet(selD.get_selected_particles())
1120 """Utility to retrieve particles from a hierarchy within a
1121 zone around a set of ps.
1122 @param hier The hierarchy in which to look for neighbors
1123 @param target_ps The particles for zoning
1124 @param sel_zone The maximum distance
1125 @param entire_residues If True, will grab entire residues
1126 @param exclude_backbone If True, will only return sidechain particles
1130 backbone_types = [
'C',
'N',
'CB',
'O']
1131 if exclude_backbone:
1134 test_ps = test_sel.get_selected_particles()
1135 nn = IMP.algebra.NearestNeighbor3D([
IMP.core.XYZ(p).get_coordinates()
1138 for target
in target_ps:
1139 zone |= set(nn.get_in_ball(
IMP.core.XYZ(target).get_coordinates(),
1141 zone_ps = [test_ps[z]
for z
in zone]
1146 zone_ps = [h.get_particle()
for h
in final_ps]
1151 """Returns unique objects in original order"""
1155 if not hasattr(hiers,
'__iter__'):
1162 rbs_ordered.append(rb)
1167 rbs_ordered.append(rb)
1171 return rbs_ordered, beads
1175 "This function returns the parent molecule hierarchies of given objects"
1176 stuff =
input_adaptor(input_objects, pmi_resolution=
'all', flatten=
True)
1181 while not (is_molecule
or is_root):
1182 root = IMP.atom.get_root(h)
1189 return list(molecules)
1192 def get_molecules_dictionary(input_objects):
1193 moldict = defaultdict(list)
1195 name = mol.get_name()
1196 moldict[name].append(mol)
1203 def get_molecules_dictionary_by_copy(input_objects):
1204 moldict = defaultdict(dict)
1206 name = mol.get_name()
1208 moldict[name][c] = mol
1212 def get_selections_dictionary(input_objects):
1213 moldict = IMP.pmi.tools.get_molecules_dictionary(input_objects)
1214 seldict = defaultdict(list)
1215 for name, mols
in moldict.items():
1222 """Given a list of PMI objects, returns all density hierarchies within
1223 these objects. The output of this function can be inputted into
1224 things such as EM restraints. This function is intended to gather
1225 density particles appended to molecules (and not other hierarchies
1226 which might have been appended to the root node directly).
1235 i, representation_type=IMP.atom.DENSITIES).get_selected_particles()
1240 max_translation=300., max_rotation=2.0 * math.pi,
1241 avoidcollision_rb=
True, avoidcollision_fb=
False,
1242 cutoff=10.0, niterations=100,
1244 excluded_rigid_bodies=[],
1245 hierarchies_excluded_from_collision=[],
1246 hierarchies_included_in_collision=[],
1248 return_debug=
False):
1249 """Shuffle particles. Used to restart the optimization.
1250 The configuration of the system is initialized by placing each
1251 rigid body and each bead randomly in a box. If `bounding_box` is
1252 specified, the particles are placed inside this box; otherwise, each
1253 particle is displaced by up to max_translation angstroms, and randomly
1254 rotated. Effort is made to place particles far enough from each other to
1255 prevent any steric clashes.
1256 @param objects Can be one of the following inputs:
1257 IMP Hierarchy, PMI System/State/Molecule/TempResidue, or
1259 @param max_translation Max translation (rbs and flexible beads)
1260 @param max_rotation Max rotation (rbs only)
1261 @param avoidcollision_rb check if the particle/rigid body was
1262 placed close to another particle; uses the optional
1263 arguments cutoff and niterations
1264 @param avoidcollision_fb Advanced. Generally you want this False because
1265 it's hard to shuffle beads.
1266 @param cutoff Distance less than this is a collision
1267 @param niterations How many times to try avoiding collision
1268 @param bounding_box Only shuffle particles within this box.
1269 Defined by ((x1,y1,z1),(x2,y2,z2)).
1270 @param excluded_rigid_bodies Don't shuffle these rigid body objects
1271 @param hierarchies_excluded_from_collision Don't count collision
1273 @param hierarchies_included_in_collision Hierarchies that are not
1274 shuffled, but should be included in collision calculation
1276 @param verbose Give more output
1277 @note Best to only call this function after you've set up degrees
1279 For debugging purposes, returns: <shuffled indexes>,
1280 <collision avoided indexes>
1285 pmi_resolution=
'all',
1288 if len(rigid_bodies) > 0:
1289 mdl = rigid_bodies[0].get_model()
1290 elif len(flexible_beads) > 0:
1291 mdl = flexible_beads[0].get_model()
1293 raise Exception(
"Could not find any particles in the hierarchy")
1294 if len(rigid_bodies) == 0:
1295 print(
"shuffle_configuration: rigid bodies were not initialized")
1299 gcpf.set_distance(cutoff)
1303 hierarchies_excluded_from_collision, pmi_resolution=
'all',
1307 hierarchies_included_in_collision, pmi_resolution=
'all', flatten=
True)
1309 collision_excluded_idxs = set(
1311 for h
in collision_excluded_hierarchies
1314 collision_included_idxs = set(
1316 for h
in collision_included_hierarchies
1323 all_idxs.append(p.get_particle_index())
1325 collision_excluded_idxs.add(p.get_particle_index())
1327 if bounding_box
is not None:
1328 ((x1, y1, z1), (x2, y2, z2)) = bounding_box
1333 all_idxs = set(all_idxs) | collision_included_idxs
1334 all_idxs = all_idxs - collision_excluded_idxs
1336 print(
'shuffling', len(rigid_bodies),
'rigid bodies')
1337 for rb
in rigid_bodies:
1338 if rb
not in excluded_rigid_bodies:
1340 if avoidcollision_rb:
1341 rb_idxs = set(rb.get_member_particle_indexes()) - \
1342 collision_excluded_idxs
1343 other_idxs = all_idxs - rb_idxs
1345 debug.append([rb, other_idxs
if avoidcollision_rb
else set()])
1348 while niter < niterations:
1349 rbxyz = (rb.get_x(), rb.get_y(), rb.get_z())
1366 rbxyz, max_translation, max_rotation)
1371 if avoidcollision_rb
and other_idxs:
1373 npairs = len(gcpf.get_close_pairs(mdl,
1381 print(
"shuffle_configuration: rigid body placed "
1382 "close to other %d particles, trying "
1383 "again..." % npairs)
1384 print(
"shuffle_configuration: rigid body name: "
1386 if niter == niterations:
1388 "tried the maximum number of iterations to "
1389 "avoid collisions, increase the distance "
1394 print(
'shuffling', len(flexible_beads),
'flexible beads')
1395 for fb
in flexible_beads:
1397 if avoidcollision_fb:
1399 other_idxs = all_idxs - fb_idxs
1405 while niter < niterations:
1412 fbxyz, max_translation, max_rotation)
1417 xyz = memb.get_internal_coordinates()
1422 rf = memb.get_rigid_body().get_reference_frame()
1423 glob_to_int = rf.get_transformation_from()
1424 memb.set_internal_coordinates(
1425 glob_to_int.get_transformed(translation))
1427 xyz_transformed = transformation.get_transformed(xyz)
1428 memb.set_internal_coordinates(xyz_transformed)
1431 [xyz, other_idxs
if avoidcollision_fb
else set()])
1439 -d.get_coordinates())
1445 [d, other_idxs
if avoidcollision_fb
else set()])
1452 if avoidcollision_fb:
1454 npairs = len(gcpf.get_close_pairs(mdl,
1462 print(
"shuffle_configuration: floppy body placed close "
1463 "to other %d particles, trying again..." % npairs)
1464 if niter == niterations:
1466 "tried the maximum number of iterations to avoid "
1467 "collisions, increase the distance cutoff")
1474 class ColorHierarchy(object):
1476 def __init__(self, hier):
1477 import matplotlib
as mpl
1479 import matplotlib.pyplot
as plt
1483 hier.ColorHierarchy = self
1488 self.method = self.nochange
1496 def get_color(self, fl):
1499 def get_log_scale(self, fl):
1502 return math.log(fl+eps)
1504 def color_by_resid(self):
1505 self.method = self.color_by_resid
1506 self.scheme = self.mpl.cm.rainbow
1507 for mol
in self.mols:
1514 c = self.get_color(float(ri)/self.last)
1518 avr = sum(ris)/len(ris)
1519 c = self.get_color(float(avr)/self.last)
1522 def color_by_uncertainty(self):
1523 self.method = self.color_by_uncertainty
1524 self.scheme = self.mpl.cm.jet
1531 self.first = self.get_log_scale(1.0)
1532 self.last = self.get_log_scale(100.0)
1534 value = self.get_log_scale(unc_dict[p])
1535 if value >= self.last:
1537 if value <= self.first:
1539 c = self.get_color((value-self.first) / (self.last-self.first))
1542 def get_color_bar(self, filename):
1543 import matplotlib
as mpl
1545 import matplotlib.pyplot
as plt
1547 fig = plt.figure(figsize=(8, 3))
1548 ax1 = fig.add_axes([0.05, 0.80, 0.9, 0.15])
1551 norm = mpl.colors.Normalize(vmin=0.0, vmax=1.0)
1553 if self.method == self.color_by_uncertainty:
1554 angticks = [1.0, 2.5, 5.0, 10.0, 25.0, 50.0, 100.0]
1558 vvalue = (self.get_log_scale(at)-self.first) \
1559 / (self.last-self.first)
1560 if vvalue <= 1.0
and vvalue >= 0.0:
1561 vvalues.append(vvalue)
1562 marks.append(str(at))
1563 cb1 = mpl.colorbar.ColorbarBase(
1564 ax1, cmap=cmap, norm=norm, ticks=vvalues,
1565 orientation=
'horizontal')
1566 print(self.first, self.last, marks, vvalues)
1567 cb1.ax.set_xticklabels(marks)
1568 cb1.set_label(
'Angstorm')
1569 plt.savefig(filename, dpi=150, transparent=
True)
1574 """Given a Chimera color name or hex color value, return RGB"""
1575 d = {
'aquamarine': (0.4980392156862745, 1.0, 0.8313725490196079),
1576 'black': (0.0, 0.0, 0.0),
1577 'blue': (0.0, 0.0, 1.0),
1578 'brown': (0.6470588235, 0.16470588235294117, 0.16470588235294117),
1579 'chartreuse': (0.4980392156862745, 1.0, 0.0),
1580 'coral': (1.0, 0.4980392156862745, 0.3137254901960784),
1581 'cornflower blue': (0.39215686, 0.58431372549, 0.9294117647058824),
1582 'cyan': (0.0, 1.0, 1.0),
1583 'dark cyan': (0.0, 0.5450980392156862, 0.5450980392156862),
1584 'dark gray': (0.6627450980, 0.6627450980392157, 0.6627450980392157),
1585 'dark green': (0.0, 0.39215686274509803, 0.0),
1586 'dark khaki': (0.74117647, 0.7176470588235294, 0.4196078431372549),
1587 'dark magenta': (0.5450980392156862, 0.0, 0.5450980392156862),
1588 'dark olive green': (0.333333333, 0.419607843, 0.1843137254901961),
1589 'dark red': (0.5450980392156862, 0.0, 0.0),
1590 'dark slate blue': (0.28235294, 0.239215686, 0.5450980392156862),
1591 'dark slate gray': (0.1843137, 0.30980392, 0.30980392156862746),
1592 'deep pink': (1.0, 0.0784313725490196, 0.5764705882352941),
1593 'deep sky blue': (0.0, 0.7490196078431373, 1.0),
1594 'dim gray': (0.41176470, 0.4117647058823529, 0.4117647058823529),
1595 'dodger blue': (0.11764705882352941, 0.5647058823529412, 1.0),
1596 'firebrick': (0.6980392, 0.13333333333333333, 0.13333333333333333),
1597 'forest green': (0.13333333, 0.5450980392156862, 0.13333333333333333),
1598 'gold': (1.0, 0.8431372549019608, 0.0),
1599 'goldenrod': (0.85490196, 0.6470588235294118, 0.12549019607843137),
1600 'gray': (0.7450980392156863, 0.7450980392156863, 0.7450980392156863),
1601 'green': (0.0, 1.0, 0.0),
1602 'hot pink': (1.0, 0.4117647058823529, 0.7058823529411765),
1603 'khaki': (0.9411764705882353, 0.9019607843137255, 0.5490196078431373),
1604 'light blue': (0.67843137, 0.8470588235294118, 0.9019607843137255),
1605 'light gray': (0.82745098, 0.8274509803921568, 0.8274509803921568),
1606 'light green': (0.56470588, 0.9333333333333333, 0.5647058823529412),
1607 'light sea green': (0.125490, 0.6980392156862745, 0.6666666666666666),
1608 'lime green': (0.1960784, 0.803921568627451, 0.19607843137254902),
1609 'magenta': (1.0, 0.0, 1.0),
1610 'medium blue': (0.1960784, 0.19607843137254902, 0.803921568627451),
1611 'medium purple': (0.57647, 0.4392156862745098, 0.8588235294117647),
1612 'navy blue': (0.0, 0.0, 0.5019607843137255),
1613 'olive drab': (0.4196078, 0.5568627450980392, 0.13725490196078433),
1614 'orange red': (1.0, 0.27058823529411763, 0.0),
1615 'orange': (1.0, 0.4980392156862745, 0.0),
1616 'orchid': (0.85490196, 0.4392156862745098, 0.8392156862745098),
1617 'pink': (1.0, 0.7529411764705882, 0.796078431372549),
1618 'plum': (0.8666666666666667, 0.6274509803921569, 0.8666666666666667),
1619 'purple': (0.62745098, 0.12549019607843137, 0.9411764705882353),
1620 'red': (1.0, 0.0, 0.0),
1621 'rosy brown': (0.7372549, 0.5607843137254902, 0.5607843137254902),
1622 'salmon': (0.980392, 0.5019607843137255, 0.4470588235294118),
1623 'sandy brown': (0.956862745, 0.6431372549019608, 0.3764705882352941),
1624 'sea green': (0.18039, 0.5450980392156862, 0.3411764705882353),
1625 'sienna': (0.6274509, 0.3215686274509804, 0.17647058823529413),
1626 'sky blue': (0.52941176, 0.807843137254902, 0.9215686274509803),
1627 'slate gray': (0.439215686, 0.50196078, 0.5647058823529412),
1628 'spring green': (0.0, 1.0, 0.4980392156862745),
1629 'steel blue': (0.2745098, 0.50980392, 0.70588235),
1630 'tan': (0.8235294117647058, 0.7058823529411765, 0.5490196078431373),
1631 'turquoise': (0.25098039, 0.87843137, 0.81568627),
1632 'violet red': (0.81568627, 0.125490196, 0.56470588235),
1633 'white': (1.0, 1.0, 1.0),
1634 'yellow': (1.0, 1.0, 0.0)}
1635 if colorname.startswith(
'#'):
1636 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.
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.
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.
Represent 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)
Get the indexes from a list of particle pairs.
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)
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 ...