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 h = hier._pmi2_system()
49 for s
in IMP.pmi.topology.System._all_systems:
53 hier = hier.get_parent()
56 def _all_protocol_outputs(hier):
57 """Iterate over all (ProtocolOutput, State) pairs for the
59 system = _get_system_for_hier(hier)
61 for state
in system.states:
62 for p
in state._protocol_output:
66 def _add_pmi_provenance(p):
67 """Tag the given particle as being created by the current version
72 location=
"https://integrativemodeling.org")
76 def _get_restraint_set_keys():
77 if not hasattr(_get_restraint_set_keys,
'pmi_rs_key'):
78 _get_restraint_set_keys.pmi_rs_key =
IMP.ModelKey(
"PMI restraints")
79 _get_restraint_set_keys.rmf_rs_key =
IMP.ModelKey(
"RMF restraints")
80 return (_get_restraint_set_keys.pmi_rs_key,
81 _get_restraint_set_keys.rmf_rs_key)
84 def _add_restraint_sets(model, mk, mk_rmf):
87 model.add_data(mk, rs)
88 model.add_data(mk_rmf, rs_rmf)
93 """Add a PMI restraint to the model.
94 Since Model.add_restraint() no longer exists (in modern IMP restraints
95 should be added to a ScoringFunction instead) store them instead in
96 a RestraintSet, and keep a reference to it in the Model.
98 If `add_to_rmf` is True, also add the restraint to a separate list
99 of restraints that will be written out to RMF files (by default, most
100 PMI restraints are not)."""
101 mk, mk_rmf = _get_restraint_set_keys()
102 if model.get_has_data(mk):
103 rs = IMP.RestraintSet.get_from(model.get_data(mk))
104 rs_rmf = IMP.RestraintSet.get_from(model.get_data(mk_rmf))
106 rs, rs_rmf = _add_restraint_sets(model, mk, mk_rmf)
107 rs.add_restraint(restraint)
109 rs_rmf.add_restraint(restraint)
113 """Get a RestraintSet containing all PMI restraints added to the model.
114 If `rmf` is True, return only the subset of these restraints that
115 should be written out to RMF files."""
116 mk, mk_rmf = _get_restraint_set_keys()
117 if not model.get_has_data(mk):
118 warnings.warn(
"no restraints added to model yet",
120 _add_restraint_sets(model, mk, mk_rmf)
122 return IMP.RestraintSet.get_from(model.get_data(mk_rmf))
124 return IMP.RestraintSet.get_from(model.get_data(mk))
128 """Collect timing information.
129 Add an instance of this class to outputobjects to get timing information
134 @param isdelta if True (the default) then report the time since the
135 last use of this class; if False, report cumulative time."""
136 self.starttime = process_time()
138 self.isdelta = isdelta
140 def set_label(self, labelstr):
141 self.label = labelstr
143 def get_output(self):
146 newtime = process_time()
147 output[
"Stopwatch_" + self.label +
"_delta_seconds"] \
148 = str(newtime - self.starttime)
149 self.starttime = newtime
151 output[
"Stopwatch_" + self.label +
"_elapsed_seconds"] \
152 = str(process_time() - self.starttime)
156 class SetupNuisance(object):
158 def __init__(self, m, initialvalue, minvalue, maxvalue, isoptimized=True,
166 nuisance.set_lower(minvalue)
168 nuisance.set_upper(maxvalue)
171 nuisance.set_is_optimized(nuisance.get_nuisance_key(), isoptimized)
172 self.nuisance = nuisance
174 def get_particle(self):
178 class SetupWeight(object):
180 def __init__(self, m, isoptimized=True, nweights_or_weights=None):
182 if isinstance(nweights_or_weights, int):
184 pw, nweights_or_weights
188 nweights_or_weights = list(nweights_or_weights)
190 pw, nweights_or_weights
194 self.weight.set_weights_are_optimized(isoptimized)
196 def get_particle(self):
200 class SetupSurface(object):
202 def __init__(self, m, center, normal, isoptimized=True):
205 self.surface.set_coordinates_are_optimized(isoptimized)
206 self.surface.set_normal_is_optimized(isoptimized)
208 def get_particle(self):
212 def get_cross_link_data(directory, filename, dist, omega, sigma,
213 don=
None, doff=
None, prior=0, type_of_profile=
"gofr"):
215 (distmin, distmax, ndist) = dist
216 (omegamin, omegamax, nomega) = omega
217 (sigmamin, sigmamax, nsigma) = sigma
220 with open(filen)
as xlpot:
221 dictionary = ast.literal_eval(xlpot.readline())
223 xpot = dictionary[directory][filename][
"distance"]
224 pot = dictionary[directory][filename][type_of_profile]
226 dist_grid =
get_grid(distmin, distmax, ndist,
False)
227 omega_grid = get_log_grid(omegamin, omegamax, nomega)
228 sigma_grid = get_log_grid(sigmamin, sigmamax, nsigma)
230 if don
is not None and doff
is not None:
250 def get_grid(gmin, gmax, ngrid, boundaries):
252 dx = (gmax - gmin) / float(ngrid)
253 for i
in range(0, ngrid + 1):
254 if not boundaries
and i == 0:
256 if not boundaries
and i == ngrid:
258 grid.append(gmin + float(i) * dx)
262 def get_log_grid(gmin, gmax, ngrid):
264 for i
in range(0, ngrid + 1):
265 grid.append(gmin * math.exp(float(i) / ngrid * math.log(gmax / gmin)))
271 example '"{ID_Score}" > 28 AND "{Sample}" ==
272 "%10_1%" OR ":Sample}" == "%10_2%" OR ":Sample}"
273 == "%10_3%" OR ":Sample}" == "%8_1%" OR ":Sample}" == "%8_2%"'
276 import pyparsing
as pp
278 operator = pp.Regex(
">=|<=|!=|>|<|==|in").setName(
"operator")
279 value = pp.QuotedString(
281 r"[+-]?\d+(:?\.\d*)?(:?[eE][+-]?\d+)?")
282 identifier = pp.Word(pp.alphas, pp.alphanums +
"_")
283 comparison_term = identifier | value
284 condition = pp.Group(comparison_term + operator + comparison_term)
286 expr = pp.operatorPrecedence(condition, [
287 (
"OR", 2, pp.opAssoc.LEFT, ),
288 (
"AND", 2, pp.opAssoc.LEFT, ),
291 parsedstring = str(expr.parseString(inputstring)) \
297 .replace(
"{",
"float(entry['") \
298 .replace(
"}",
"'])") \
299 .replace(
":",
"str(entry['") \
300 .replace(
"}",
"'])") \
301 .replace(
"AND",
"and") \
306 def open_file_or_inline_text(filename):
308 fl = open(filename,
"r")
310 fl = filename.split(
"\n")
314 def get_ids_from_fasta_file(fastafile):
316 with open(fastafile)
as ff:
319 ids.append(line[1:-1])
325 this function works with plain hierarchies, as read from the pdb,
326 no multi-scale hierarchies
333 atom_type=IMP.atom.AT_CA)
341 print(
"get_closest_residue_position: exiting while loop "
344 p = sel.get_selected_particles()
349 print(
"get_closest_residue_position: got NO residues for hierarchy "
350 "%s and residue %i" % (hier, resindex))
352 "get_closest_residue_position: got NO residues for hierarchy "
353 "%s and residue %i" % (hier, resindex))
356 "got multiple residues for hierarchy %s and residue %i; the list "
358 % (hier, resindex, str([pp.get_name()
for pp
in p])))
363 Return the residue index gaps and contiguous segments in the hierarchy.
365 @param hierarchy hierarchy to examine
366 @param start first residue index
367 @param end last residue index
369 @return A list of lists of the form
370 [[1,100,"cont"],[101,120,"gap"],[121,200,"cont"]]
373 for n, rindex
in enumerate(range(start, end + 1)):
375 atom_type=IMP.atom.AT_CA)
377 if len(sel.get_selected_particles()) == 0:
381 rindexcont = start - 1
382 if rindexgap == rindex - 1:
388 gaps.append([rindex, rindex,
"gap"])
394 rindexgap = start - 1
396 if rindexcont == rindex - 1:
403 gaps.append([rindex, rindex,
"cont"])
414 def set_map_element(self, xvalue, yvalue):
415 self.map[xvalue] = yvalue
417 def get_map_element(self, invalue):
418 if isinstance(invalue, float):
422 dist = (invalue - x) * (invalue - x)
431 return self.map[minx]
432 elif isinstance(invalue, str):
433 return self.map[invalue]
435 raise TypeError(
"wrong type for map")
439 """New tuple format: molname OR (start,stop,molname,copynum,statenum)
440 Copy and state are optional. Can also use 'None' for them which will
441 get all. You can also pass -1 for stop which will go to the end.
442 Returns the particles
445 kwds[
'resolution'] = resolution
446 if isinstance(tuple_selection, str):
447 kwds[
'molecule'] = tuple_selection
448 elif isinstance(tuple_selection, tuple):
449 rbegin = tuple_selection[0]
450 rend = tuple_selection[1]
451 kwds[
'molecule'] = tuple_selection[2]
453 copynum = tuple_selection[3]
454 if copynum
is not None:
455 kwds[
'copy_index'] = copynum
459 statenum = tuple_selection[4]
460 if statenum
is not None:
461 kwds[
'state_index'] = statenum
468 residue_indexes=range(1, rbegin),
470 return s.get_selected_particles()
472 kwds[
'residue_indexes'] = range(rbegin, rend+1)
474 return s.get_selected_particles()
477 def get_db_from_csv(csvfilename, encoding=None):
478 if sys.version_info[0] == 2:
479 def open_with_encoding(fname, encoding):
482 open_with_encoding = open
485 with open_with_encoding(csvfilename, encoding=encoding)
as fh:
486 csvr = csv.DictReader(fh)
488 outputlist.append(ls)
493 '''Return the component name provided a particle and a list of names'''
495 protname = root.get_name()
497 while protname
not in list_of_names:
498 root0 = root.get_parent()
501 protname = root0.get_name()
506 if "Beads" in protname:
509 return (protname, is_a_bead)
514 Retrieve the residue indexes for the given particle.
516 The particle must be an instance of Fragment,Residue, Atom or Molecule
517 or else returns an empty list
528 resind_tmp = IMP.pmi.tools.OrderedSet()
535 resind = list(resind_tmp)
541 def sort_by_residues(particles):
544 sorted_particles_residues = sorted(
546 key=
lambda tup: tup[1])
547 particles = [p[0]
for p
in sorted_particles_residues]
556 """Synchronize data over a parallel run"""
557 from mpi4py
import MPI
558 comm = MPI.COMM_WORLD
559 rank = comm.Get_rank()
560 number_of_processes = comm.size
563 comm.send(data, dest=0, tag=11)
566 for i
in range(1, number_of_processes):
567 data_tmp = comm.recv(source=i, tag=11)
568 if isinstance(data, list):
570 elif isinstance(data, dict):
571 data.update(data_tmp)
573 raise TypeError(
"data not supported, use list or dictionaries")
575 for i
in range(1, number_of_processes):
576 comm.send(data, dest=i, tag=11)
579 data = comm.recv(source=0, tag=11)
589 Yield all sublists of length >= lmin and <= lmax
595 for j
in range(i + lmin, min(n + 1, i + 1 + lmax)):
599 def flatten_list(ls):
600 return [item
for sublist
in ls
for item
in sublist]
604 """ Yield successive length-sized chunks from a list.
606 for i
in range(0, len(list), length):
607 yield list[i:i + length]
610 def chunk_list_into_segments(seq, num):
612 avg = len(seq) / float(num)
616 while last < len(seq):
617 out.append(seq[int(last):int(last + avg)])
625 ''' This class stores integers
626 in ordered compact lists eg:
628 the methods help splitting and merging the internal lists
630 s=Segments([1,2,3]) is [[1,2,3]]
631 s.add(4) is [[1,2,3,4]] (add right)
632 s.add(3) is [[1,2,3,4]] (item already existing)
633 s.add(7) is [[1,2,3,4],[7]] (new list)
634 s.add([8,9]) is [[1,2,3,4],[7,8,9]] (add item right)
635 s.add([5,6]) is [[1,2,3,4,5,6,7,8,9]] (merge)
636 s.remove(3) is [[1,2],[4,5,6,7,8,9]] (split)
641 '''index can be a integer or a list of integers '''
642 if isinstance(index, int):
643 self.segs = [[index]]
644 elif isinstance(index, list):
645 self.segs = [[index[0]]]
649 raise TypeError(
"index must be an int or list of ints")
652 '''index can be a integer or a list of integers '''
653 if isinstance(index, (int, numpy.int32, numpy.int64)):
656 for n, s
in enumerate(self.segs):
664 if mergeright
is None and mergeleft
is None:
665 self.segs.append([index])
666 if mergeright
is not None and mergeleft
is None:
667 self.segs[mergeright].append(index)
668 if mergeleft
is not None and mergeright
is None:
669 self.segs[mergeleft] = [index]+self.segs[mergeleft]
670 if mergeleft
is not None and mergeright
is not None:
671 self.segs[mergeright] = \
672 self.segs[mergeright]+[index]+self.segs[mergeleft]
673 del self.segs[mergeleft]
675 for n
in range(len(self.segs)):
678 self.segs.sort(key=
lambda tup: tup[0])
680 elif isinstance(index, list):
684 raise TypeError(
"index must be an int or list of ints")
687 '''index can be a integer'''
688 for n, s
in enumerate(self.segs):
693 self.segs[n] = s[:-1]
695 i = self.segs[n].index(index)
697 self.segs.append(s[i+1:])
698 for n
in range(len(self.segs)):
700 if len(self.segs[n]) == 0:
702 self.segs.sort(key=
lambda tup: tup[0])
705 ''' Returns a flatten list '''
706 return [item
for sublist
in self.segs
for item
in sublist]
710 for seg
in self.segs:
711 ret_tmp += str(seg[0])+
"-"+str(seg[-1])+
","
712 ret = ret_tmp[:-1]+
"]"
720 def normal_density_function(expected_value, sigma, x):
722 1 / math.sqrt(2 * math.pi) / sigma *
723 math.exp(-(x - expected_value) ** 2 / 2 / sigma / sigma)
727 def log_normal_density_function(expected_value, sigma, x):
729 1 / math.sqrt(2 * math.pi) / sigma / x *
730 math.exp(-(math.log(x / expected_value) ** 2 / 2 / sigma / sigma))
734 def print_multicolumn(list_of_strings, ncolumns=2, truncate=40):
740 for i
in range(len(ls) % cols):
743 split = [ls[i:i + len(ls) // cols]
744 for i
in range(0, len(ls), len(ls) // cols)]
745 for row
in zip(*split):
746 print(
"".join(str.ljust(i, truncate)
for i
in row))
750 '''Change color code to hexadecimal to rgb'''
752 self._NUMERALS =
'0123456789abcdefABCDEF'
753 self._HEXDEC = dict((v, int(v, 16))
for v
in
754 (x+y
for x
in self._NUMERALS
755 for y
in self._NUMERALS))
756 self.LOWERCASE, self.UPPERCASE =
'x',
'X'
758 def rgb(self, triplet):
759 return (float(self._HEXDEC[triplet[0:2]]),
760 float(self._HEXDEC[triplet[2:4]]),
761 float(self._HEXDEC[triplet[4:6]]))
763 def triplet(self, rgb, lettercase=None):
764 if lettercase
is None:
765 lettercase = self.LOWERCASE
766 return format(rgb[0] << 16 | rgb[1] << 8 | rgb[2],
'06'+lettercase)
770 class OrderedSet(MutableSet):
772 def __init__(self, iterable=None):
774 end += [
None, end, end]
776 if iterable
is not None:
782 def __contains__(self, key):
783 return key
in self.map
786 if key
not in self.map:
789 curr[2] = end[1] = self.map[key] = [key, curr, end]
791 def discard(self, key):
793 key, prev, next = self.map.pop(key)
800 while curr
is not end:
804 def __reversed__(self):
807 while curr
is not end:
811 def pop(self, last=True):
813 raise KeyError(
'set is empty')
823 return '%s()' % (self.__class__.__name__,)
824 return '%s(%r)' % (self.__class__.__name__, list(self))
826 def __eq__(self, other):
827 if isinstance(other, OrderedSet):
828 return len(self) == len(other)
and list(self) == list(other)
829 return set(self) == set(other)
833 """Store objects in order they were added, but with default type.
834 Source: http://stackoverflow.com/a/4127426/2608793
836 def __init__(self, *args, **kwargs):
838 self.default_factory =
None
840 if not (args[0]
is None or callable(args[0])):
841 raise TypeError(
'first argument must be callable or None')
842 self.default_factory = args[0]
844 super(OrderedDefaultDict, self).__init__(*args, **kwargs)
846 def __missing__(self, key):
847 if self.default_factory
is None:
849 self[key] = default = self.default_factory()
852 def __reduce__(self):
853 args = (self.default_factory,)
if self.default_factory
else ()
854 if sys.version_info[0] >= 3:
855 return self.__class__, args,
None,
None, self.items()
857 return self.__class__, args,
None,
None, self.iteritems()
863 """Extract frame from RMF file and fill coordinates. Must be identical
866 @param hier The (System) hierarchy to fill (e.g. after you've built it)
867 @param rmf_fn The file to extract from
868 @param frame_num The frame number to extract
870 rh = RMF.open_rmf_file_read_only(rmf_fn)
876 def input_adaptor(stuff, pmi_resolution=0, flatten=False, selection_tuple=None,
877 warn_about_slices=
True):
878 """Adapt things for PMI (degrees of freedom, restraints, ...)
879 Returns list of list of hierarchies, separated into Molecules if possible.
880 The input can be a list, or a list of lists (iterable of ^1 or
882 (iterable of ^2) Hierarchy -> returns input as list of list of hierarchies,
883 only one entry, not grouped by molecules.
884 (iterable of ^2) PMI::System/State/Molecule/TempResidue ->
885 returns residue hierarchies, grouped in molecules, at requested
888 @param stuff Can be one of the following inputs:
889 IMP Hierarchy, PMI System/State/Molecule/TempResidue, or a
890 list/set (of list/set) of them.
891 Must be uniform input, however. No mixing object types.
892 @param pmi_resolution For selecting, only does it if you pass PMI
893 objects. Set it to "all" if you want all resolutions!
894 @param flatten Set to True if you just want all hierarchies in one list.
895 @param warn_about_slices Print a warning if you are requesting only part
896 of a bead. Sometimes you just don't care!
897 @note since this relies on IMP::atom::Selection, this will not return
898 any objects if they weren't built! But there should be no problem
899 if you request unbuilt residues - they should be ignored.
905 if hasattr(stuff,
'__iter__'):
908 thelist = list(stuff)
911 if all(hasattr(el,
'__iter__')
for el
in thelist):
912 thelist = [i
for sublist
in thelist
for i
in sublist]
913 elif any(hasattr(el,
'__iter__')
for el
in thelist):
914 raise Exception(
'input_adaptor: input_object must be a list '
915 'or a list of lists')
924 except (NotImplementedError, TypeError):
936 if is_system
or is_state
or is_molecule
or is_temp_residue:
943 for state
in system.get_states():
944 mdict = state.get_molecules()
945 for molname
in mdict:
946 for copy
in mdict[molname]:
947 indexes_per_mol[copy] += \
948 [r.get_index()
for r
in copy.get_residues()]
951 mdict = state.get_molecules()
952 for molname
in mdict:
953 for copy
in mdict[molname]:
954 indexes_per_mol[copy] += [r.get_index()
955 for r
in copy.get_residues()]
957 for molecule
in stuff:
958 indexes_per_mol[molecule] += [r.get_index()
959 for r
in molecule.get_residues()]
961 for tempres
in stuff:
962 indexes_per_mol[tempres.get_molecule()].append(
964 for mol
in indexes_per_mol:
965 if pmi_resolution ==
'all':
969 mol.get_hierarchy(), residue_indexes=indexes_per_mol[mol])
972 resolution=pmi_resolution,
973 residue_indexes=indexes_per_mol[mol])
974 ps = sel.get_selected_particles()
977 if warn_about_slices:
978 rset = set(indexes_per_mol[mol])
988 resbreak = maxf
if minf == minset
else minset-1
990 'You are trying to select only part of the '
991 'bead %s:%i-%i. The residues you requested '
992 'are %i-%i. You can fix this by: '
993 '1) requesting the whole bead/none of it; or'
994 '2) break the bead up by passing '
995 'bead_extra_breaks=[\'%i\'] in '
996 'molecule.add_representation()'
997 % (mol.get_name(), minset, maxset, minf, maxf,
1003 if pmi_resolution ==
'all':
1009 h, resolution=pmi_resolution).get_selected_particles()
1012 hier_list = [hier_list]
1014 raise Exception(
'input_adaptor: you passed something of wrong type '
1015 'or a list with mixed types')
1017 if flatten
and pmi_input:
1018 return [h
for sublist
in hier_list
for h
in sublist]
1024 """Returns sequence-sorted segments array, each containing the first
1025 particle the last particle and the first residue index."""
1027 from operator
import itemgetter
1030 raise ValueError(
"only pass stuff from one Molecule, please")
1045 segs.append((start, end, startres))
1046 return sorted(segs, key=itemgetter(2))
1050 """Decorate the sequence-consecutive particles from a PMI2 molecule
1051 with a bond, so that they appear connected in the rmf file"""
1053 for x
in range(len(SortedSegments) - 1):
1055 last = SortedSegments[x][1]
1056 first = SortedSegments[x + 1][0]
1058 p1 = last.get_particle()
1059 p2 = first.get_particle()
1072 """ Just get the leaves from a list of hierarchies """
1073 lvs = list(itertools.chain.from_iterable(
1079 """Perform selection using the usual keywords but return ALL
1080 resolutions (BEADS and GAUSSIANS).
1081 Returns in flat list!
1086 if hier
is not None:
1089 warnings.warn(
"You passed nothing to select_at_all_resolutions()",
1097 raise Exception(
'select_at_all_resolutions: you have to pass '
1100 raise Exception(
'select_at_all_resolutions: you have to pass '
1102 if 'resolution' in kwargs
or 'representation_type' in kwargs:
1103 raise Exception(
"don't pass resolution or representation_type "
1106 representation_type=IMP.atom.BALLS,
1109 representation_type=IMP.atom.DENSITIES,
1111 ret |= OrderedSet(selB.get_selected_particles())
1112 ret |= OrderedSet(selD.get_selected_particles())
1121 """Utility to retrieve particles from a hierarchy within a
1122 zone around a set of ps.
1123 @param hier The hierarchy in which to look for neighbors
1124 @param target_ps The particles for zoning
1125 @param sel_zone The maximum distance
1126 @param entire_residues If True, will grab entire residues
1127 @param exclude_backbone If True, will only return sidechain particles
1131 backbone_types = [
'C',
'N',
'CB',
'O']
1132 if exclude_backbone:
1135 test_ps = test_sel.get_selected_particles()
1136 nn = IMP.algebra.NearestNeighbor3D([
IMP.core.XYZ(p).get_coordinates()
1139 for target
in target_ps:
1140 zone |= set(nn.get_in_ball(
IMP.core.XYZ(target).get_coordinates(),
1142 zone_ps = [test_ps[z]
for z
in zone]
1147 zone_ps = [h.get_particle()
for h
in final_ps]
1152 """Returns unique objects in original order"""
1156 if not hasattr(hiers,
'__iter__'):
1163 rbs_ordered.append(rb)
1168 rbs_ordered.append(rb)
1172 return rbs_ordered, beads
1176 "This function returns the parent molecule hierarchies of given objects"
1177 stuff =
input_adaptor(input_objects, pmi_resolution=
'all', flatten=
True)
1182 while not (is_molecule
or is_root):
1183 root = IMP.atom.get_root(h)
1190 return list(molecules)
1193 def get_molecules_dictionary(input_objects):
1194 moldict = defaultdict(list)
1196 name = mol.get_name()
1197 moldict[name].append(mol)
1204 def get_molecules_dictionary_by_copy(input_objects):
1205 moldict = defaultdict(dict)
1207 name = mol.get_name()
1209 moldict[name][c] = mol
1213 def get_selections_dictionary(input_objects):
1214 moldict = IMP.pmi.tools.get_molecules_dictionary(input_objects)
1215 seldict = defaultdict(list)
1216 for name, mols
in moldict.items():
1223 """Given a list of PMI objects, returns all density hierarchies within
1224 these objects. The output of this function can be inputted into
1225 things such as EM restraints. This function is intended to gather
1226 density particles appended to molecules (and not other hierarchies
1227 which might have been appended to the root node directly).
1236 i, representation_type=IMP.atom.DENSITIES).get_selected_particles()
1241 max_translation=300., max_rotation=2.0 * math.pi,
1242 avoidcollision_rb=
True, avoidcollision_fb=
False,
1243 cutoff=10.0, niterations=100,
1245 excluded_rigid_bodies=[],
1246 hierarchies_excluded_from_collision=[],
1247 hierarchies_included_in_collision=[],
1249 return_debug=
False):
1250 """Shuffle particles. Used to restart the optimization.
1251 The configuration of the system is initialized by placing each
1252 rigid body and each bead randomly in a box. If `bounding_box` is
1253 specified, the particles are placed inside this box; otherwise, each
1254 particle is displaced by up to max_translation angstroms, and randomly
1255 rotated. Effort is made to place particles far enough from each other to
1256 prevent any steric clashes.
1257 @param objects Can be one of the following inputs:
1258 IMP Hierarchy, PMI System/State/Molecule/TempResidue, or
1260 @param max_translation Max translation (rbs and flexible beads)
1261 @param max_rotation Max rotation (rbs only)
1262 @param avoidcollision_rb check if the particle/rigid body was
1263 placed close to another particle; uses the optional
1264 arguments cutoff and niterations
1265 @param avoidcollision_fb Advanced. Generally you want this False because
1266 it's hard to shuffle beads.
1267 @param cutoff Distance less than this is a collision
1268 @param niterations How many times to try avoiding collision
1269 @param bounding_box Only shuffle particles within this box.
1270 Defined by ((x1,y1,z1),(x2,y2,z2)).
1271 @param excluded_rigid_bodies Don't shuffle these rigid body objects
1272 @param hierarchies_excluded_from_collision Don't count collision
1274 @param hierarchies_included_in_collision Hierarchies that are not
1275 shuffled, but should be included in collision calculation
1277 @param verbose Give more output
1278 @note Best to only call this function after you've set up degrees
1280 For debugging purposes, returns: <shuffled indexes>,
1281 <collision avoided indexes>
1286 pmi_resolution=
'all',
1289 if len(rigid_bodies) > 0:
1290 mdl = rigid_bodies[0].get_model()
1291 elif len(flexible_beads) > 0:
1292 mdl = flexible_beads[0].get_model()
1294 raise Exception(
"Could not find any particles in the hierarchy")
1295 if len(rigid_bodies) == 0:
1296 print(
"shuffle_configuration: rigid bodies were not initialized")
1300 gcpf.set_distance(cutoff)
1304 hierarchies_excluded_from_collision, pmi_resolution=
'all',
1308 hierarchies_included_in_collision, pmi_resolution=
'all', flatten=
True)
1310 collision_excluded_idxs = set(
1312 for h
in collision_excluded_hierarchies
1315 collision_included_idxs = set(
1317 for h
in collision_included_hierarchies
1324 all_idxs.append(p.get_particle_index())
1326 collision_excluded_idxs.add(p.get_particle_index())
1328 if bounding_box
is not None:
1329 ((x1, y1, z1), (x2, y2, z2)) = bounding_box
1334 all_idxs = set(all_idxs) | collision_included_idxs
1335 all_idxs = all_idxs - collision_excluded_idxs
1337 print(
'shuffling', len(rigid_bodies),
'rigid bodies')
1338 for rb
in rigid_bodies:
1339 if rb
not in excluded_rigid_bodies:
1341 if avoidcollision_rb:
1342 rb_idxs = set(rb.get_member_particle_indexes()) - \
1343 collision_excluded_idxs
1344 other_idxs = all_idxs - rb_idxs
1346 debug.append([rb, other_idxs
if avoidcollision_rb
else set()])
1349 while niter < niterations:
1350 rbxyz = (rb.get_x(), rb.get_y(), rb.get_z())
1367 rbxyz, max_translation, max_rotation)
1372 if avoidcollision_rb
and other_idxs:
1374 npairs = len(gcpf.get_close_pairs(mdl,
1382 print(
"shuffle_configuration: rigid body placed "
1383 "close to other %d particles, trying "
1384 "again..." % npairs)
1385 print(
"shuffle_configuration: rigid body name: "
1387 if niter == niterations:
1389 "tried the maximum number of iterations to "
1390 "avoid collisions, increase the distance "
1395 print(
'shuffling', len(flexible_beads),
'flexible beads')
1396 for fb
in flexible_beads:
1398 if avoidcollision_fb:
1400 other_idxs = all_idxs - fb_idxs
1406 while niter < niterations:
1413 fbxyz, max_translation, max_rotation)
1418 xyz = memb.get_internal_coordinates()
1423 rf = memb.get_rigid_body().get_reference_frame()
1424 glob_to_int = rf.get_transformation_from()
1425 memb.set_internal_coordinates(
1426 glob_to_int.get_transformed(translation))
1428 xyz_transformed = transformation.get_transformed(xyz)
1429 memb.set_internal_coordinates(xyz_transformed)
1432 [xyz, other_idxs
if avoidcollision_fb
else set()])
1440 -d.get_coordinates())
1446 [d, other_idxs
if avoidcollision_fb
else set()])
1453 if avoidcollision_fb:
1455 npairs = len(gcpf.get_close_pairs(mdl,
1463 print(
"shuffle_configuration: floppy body placed close "
1464 "to other %d particles, trying again..." % npairs)
1465 if niter == niterations:
1467 "tried the maximum number of iterations to avoid "
1468 "collisions, increase the distance cutoff")
1475 class ColorHierarchy(object):
1477 def __init__(self, hier):
1478 import matplotlib
as mpl
1480 import matplotlib.pyplot
as plt
1484 hier.ColorHierarchy = self
1489 self.method = self.nochange
1497 def get_color(self, fl):
1500 def get_log_scale(self, fl):
1503 return math.log(fl+eps)
1505 def color_by_resid(self):
1506 self.method = self.color_by_resid
1507 self.scheme = self.mpl.cm.rainbow
1508 for mol
in self.mols:
1515 c = self.get_color(float(ri)/self.last)
1519 avr = sum(ris)/len(ris)
1520 c = self.get_color(float(avr)/self.last)
1523 def color_by_uncertainty(self):
1524 self.method = self.color_by_uncertainty
1525 self.scheme = self.mpl.cm.jet
1532 self.first = self.get_log_scale(1.0)
1533 self.last = self.get_log_scale(100.0)
1535 value = self.get_log_scale(unc_dict[p])
1536 if value >= self.last:
1538 if value <= self.first:
1540 c = self.get_color((value-self.first) / (self.last-self.first))
1543 def get_color_bar(self, filename):
1544 import matplotlib
as mpl
1546 import matplotlib.pyplot
as plt
1548 fig = plt.figure(figsize=(8, 3))
1549 ax1 = fig.add_axes([0.05, 0.80, 0.9, 0.15])
1552 norm = mpl.colors.Normalize(vmin=0.0, vmax=1.0)
1554 if self.method == self.color_by_uncertainty:
1555 angticks = [1.0, 2.5, 5.0, 10.0, 25.0, 50.0, 100.0]
1559 vvalue = (self.get_log_scale(at)-self.first) \
1560 / (self.last-self.first)
1561 if vvalue <= 1.0
and vvalue >= 0.0:
1562 vvalues.append(vvalue)
1563 marks.append(str(at))
1564 cb1 = mpl.colorbar.ColorbarBase(
1565 ax1, cmap=cmap, norm=norm, ticks=vvalues,
1566 orientation=
'horizontal')
1567 print(self.first, self.last, marks, vvalues)
1568 cb1.ax.set_xticklabels(marks)
1569 cb1.set_label(
'Angstorm')
1570 plt.savefig(filename, dpi=150, transparent=
True)
1575 """Given a Chimera color name or hex color value, return RGB"""
1576 d = {
'aquamarine': (0.4980392156862745, 1.0, 0.8313725490196079),
1577 'black': (0.0, 0.0, 0.0),
1578 'blue': (0.0, 0.0, 1.0),
1579 'brown': (0.6470588235, 0.16470588235294117, 0.16470588235294117),
1580 'chartreuse': (0.4980392156862745, 1.0, 0.0),
1581 'coral': (1.0, 0.4980392156862745, 0.3137254901960784),
1582 'cornflower blue': (0.39215686, 0.58431372549, 0.9294117647058824),
1583 'cyan': (0.0, 1.0, 1.0),
1584 'dark cyan': (0.0, 0.5450980392156862, 0.5450980392156862),
1585 'dark gray': (0.6627450980, 0.6627450980392157, 0.6627450980392157),
1586 'dark green': (0.0, 0.39215686274509803, 0.0),
1587 'dark khaki': (0.74117647, 0.7176470588235294, 0.4196078431372549),
1588 'dark magenta': (0.5450980392156862, 0.0, 0.5450980392156862),
1589 'dark olive green': (0.333333333, 0.419607843, 0.1843137254901961),
1590 'dark red': (0.5450980392156862, 0.0, 0.0),
1591 'dark slate blue': (0.28235294, 0.239215686, 0.5450980392156862),
1592 'dark slate gray': (0.1843137, 0.30980392, 0.30980392156862746),
1593 'deep pink': (1.0, 0.0784313725490196, 0.5764705882352941),
1594 'deep sky blue': (0.0, 0.7490196078431373, 1.0),
1595 'dim gray': (0.41176470, 0.4117647058823529, 0.4117647058823529),
1596 'dodger blue': (0.11764705882352941, 0.5647058823529412, 1.0),
1597 'firebrick': (0.6980392, 0.13333333333333333, 0.13333333333333333),
1598 'forest green': (0.13333333, 0.5450980392156862, 0.13333333333333333),
1599 'gold': (1.0, 0.8431372549019608, 0.0),
1600 'goldenrod': (0.85490196, 0.6470588235294118, 0.12549019607843137),
1601 'gray': (0.7450980392156863, 0.7450980392156863, 0.7450980392156863),
1602 'green': (0.0, 1.0, 0.0),
1603 'hot pink': (1.0, 0.4117647058823529, 0.7058823529411765),
1604 'khaki': (0.9411764705882353, 0.9019607843137255, 0.5490196078431373),
1605 'light blue': (0.67843137, 0.8470588235294118, 0.9019607843137255),
1606 'light gray': (0.82745098, 0.8274509803921568, 0.8274509803921568),
1607 'light green': (0.56470588, 0.9333333333333333, 0.5647058823529412),
1608 'light sea green': (0.125490, 0.6980392156862745, 0.6666666666666666),
1609 'lime green': (0.1960784, 0.803921568627451, 0.19607843137254902),
1610 'magenta': (1.0, 0.0, 1.0),
1611 'medium blue': (0.1960784, 0.19607843137254902, 0.803921568627451),
1612 'medium purple': (0.57647, 0.4392156862745098, 0.8588235294117647),
1613 'navy blue': (0.0, 0.0, 0.5019607843137255),
1614 'olive drab': (0.4196078, 0.5568627450980392, 0.13725490196078433),
1615 'orange red': (1.0, 0.27058823529411763, 0.0),
1616 'orange': (1.0, 0.4980392156862745, 0.0),
1617 'orchid': (0.85490196, 0.4392156862745098, 0.8392156862745098),
1618 'pink': (1.0, 0.7529411764705882, 0.796078431372549),
1619 'plum': (0.8666666666666667, 0.6274509803921569, 0.8666666666666667),
1620 'purple': (0.62745098, 0.12549019607843137, 0.9411764705882353),
1621 'red': (1.0, 0.0, 0.0),
1622 'rosy brown': (0.7372549, 0.5607843137254902, 0.5607843137254902),
1623 'salmon': (0.980392, 0.5019607843137255, 0.4470588235294118),
1624 'sandy brown': (0.956862745, 0.6431372549019608, 0.3764705882352941),
1625 'sea green': (0.18039, 0.5450980392156862, 0.3411764705882353),
1626 'sienna': (0.6274509, 0.3215686274509804, 0.17647058823529413),
1627 'sky blue': (0.52941176, 0.807843137254902, 0.9215686274509803),
1628 'slate gray': (0.439215686, 0.50196078, 0.5647058823529412),
1629 'spring green': (0.0, 1.0, 0.4980392156862745),
1630 'steel blue': (0.2745098, 0.50980392, 0.70588235),
1631 'tan': (0.8235294117647058, 0.7058823529411765, 0.5490196078431373),
1632 'turquoise': (0.25098039, 0.87843137, 0.81568627),
1633 'violet red': (0.81568627, 0.125490196, 0.56470588235),
1634 'white': (1.0, 1.0, 1.0),
1635 'yellow': (1.0, 1.0, 0.0)}
1636 if colorname.startswith(
'#'):
1637 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 ...