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 class ParticleToSampleList(object):
215 self.dictionary_particle_type = {}
216 self.dictionary_particle_transformation = {}
217 self.dictionary_particle_name = {}
220 def add_particle(self, particle, particle_type, particle_transformation,
222 if particle_type
not in [
"Rigid_Bodies",
"Floppy_Bodies",
"Nuisances",
223 "X_coord",
"Weights",
"Surfaces"]:
224 raise TypeError(
"not the right particle type")
226 self.dictionary_particle_type[particle] = particle_type
227 if particle_type ==
"Rigid_Bodies":
228 if (isinstance(particle_transformation, tuple)
229 and len(particle_transformation) == 2
230 and all(isinstance(x, float)
231 for x
in particle_transformation)):
232 self.dictionary_particle_transformation[
233 particle] = particle_transformation
234 self.dictionary_particle_name[particle] = name
237 "ParticleToSampleList: not the right transformation "
238 "format for Rigid_Bodies, should be a tuple of floats")
239 elif particle_type ==
"Surfaces":
240 if (isinstance(particle_transformation, tuple)
241 and len(particle_transformation) == 3
242 and all(isinstance(x, float)
243 for x
in particle_transformation)):
244 self.dictionary_particle_transformation[
245 particle] = particle_transformation
246 self.dictionary_particle_name[particle] = name
249 "ParticleToSampleList: not the right transformation "
250 "format for Surfaces, should be a tuple of floats")
252 if isinstance(particle_transformation, float):
253 self.dictionary_particle_transformation[
254 particle] = particle_transformation
255 self.dictionary_particle_name[particle] = name
258 "ParticleToSampleList: not the right transformation "
259 "format, should be a float")
261 def get_particles_to_sample(self):
263 for particle
in self.dictionary_particle_type:
264 key = self.dictionary_particle_type[particle] + \
265 "ParticleToSampleList_" + \
266 self.dictionary_particle_name[particle] +
"_" + self.label
269 self.dictionary_particle_transformation[particle])
274 def get_cross_link_data(directory, filename, dist, omega, sigma,
275 don=
None, doff=
None, prior=0, type_of_profile=
"gofr"):
277 (distmin, distmax, ndist) = dist
278 (omegamin, omegamax, nomega) = omega
279 (sigmamin, sigmamax, nsigma) = sigma
282 with open(filen)
as xlpot:
283 dictionary = ast.literal_eval(xlpot.readline())
285 xpot = dictionary[directory][filename][
"distance"]
286 pot = dictionary[directory][filename][type_of_profile]
288 dist_grid =
get_grid(distmin, distmax, ndist,
False)
289 omega_grid = get_log_grid(omegamin, omegamax, nomega)
290 sigma_grid = get_log_grid(sigmamin, sigmamax, nsigma)
292 if don
is not None and doff
is not None:
312 def get_grid(gmin, gmax, ngrid, boundaries):
314 dx = (gmax - gmin) / float(ngrid)
315 for i
in range(0, ngrid + 1):
316 if(
not boundaries
and i == 0):
318 if(
not boundaries
and i == ngrid):
320 grid.append(gmin + float(i) * dx)
324 def get_log_grid(gmin, gmax, ngrid):
326 for i
in range(0, ngrid + 1):
327 grid.append(gmin * math.exp(float(i) / ngrid * math.log(gmax / gmin)))
333 example '"{ID_Score}" > 28 AND "{Sample}" ==
334 "%10_1%" OR ":Sample}" == "%10_2%" OR ":Sample}"
335 == "%10_3%" OR ":Sample}" == "%8_1%" OR ":Sample}" == "%8_2%"'
338 import pyparsing
as pp
340 operator = pp.Regex(
">=|<=|!=|>|<|==|in").setName(
"operator")
341 value = pp.QuotedString(
343 r"[+-]?\d+(:?\.\d*)?(:?[eE][+-]?\d+)?")
344 identifier = pp.Word(pp.alphas, pp.alphanums +
"_")
345 comparison_term = identifier | value
346 condition = pp.Group(comparison_term + operator + comparison_term)
348 expr = pp.operatorPrecedence(condition, [
349 (
"OR", 2, pp.opAssoc.LEFT, ),
350 (
"AND", 2, pp.opAssoc.LEFT, ),
353 parsedstring = str(expr.parseString(inputstring)) \
359 .replace(
"{",
"float(entry['") \
360 .replace(
"}",
"'])") \
361 .replace(
":",
"str(entry['") \
362 .replace(
"}",
"'])") \
363 .replace(
"AND",
"and") \
368 def open_file_or_inline_text(filename):
370 fl = open(filename,
"r")
372 fl = filename.split(
"\n")
376 def get_ids_from_fasta_file(fastafile):
378 with open(fastafile)
as ff:
381 ids.append(line[1:-1])
387 this function works with plain hierarchies, as read from the pdb,
388 no multi-scale hierarchies
395 atom_type=IMP.atom.AT_CA)
403 print(
"get_closest_residue_position: exiting while loop "
406 p = sel.get_selected_particles()
411 print(
"get_closest_residue_position: got NO residues for hierarchy "
412 "%s and residue %i" % (hier, resindex))
414 "get_closest_residue_position: got NO residues for hierarchy "
415 "%s and residue %i" % (hier, resindex))
418 "got multiple residues for hierarchy %s and residue %i; the list "
420 % (hier, resindex, str([pp.get_name()
for pp
in p])))
425 Return the residue index gaps and contiguous segments in the hierarchy.
427 @param hierarchy hierarchy to examine
428 @param start first residue index
429 @param end last residue index
431 @return A list of lists of the form
432 [[1,100,"cont"],[101,120,"gap"],[121,200,"cont"]]
435 for n, rindex
in enumerate(range(start, end + 1)):
437 atom_type=IMP.atom.AT_CA)
439 if len(sel.get_selected_particles()) == 0:
443 rindexcont = start - 1
444 if rindexgap == rindex - 1:
450 gaps.append([rindex, rindex,
"gap"])
456 rindexgap = start - 1
458 if rindexcont == rindex - 1:
465 gaps.append([rindex, rindex,
"cont"])
476 def set_map_element(self, xvalue, yvalue):
477 self.map[xvalue] = yvalue
479 def get_map_element(self, invalue):
480 if isinstance(invalue, float):
484 dist = (invalue - x) * (invalue - x)
493 return self.map[minx]
494 elif isinstance(invalue, str):
495 return self.map[invalue]
497 raise TypeError(
"wrong type for map")
501 """New tuple format: molname OR (start,stop,molname,copynum,statenum)
502 Copy and state are optional. Can also use 'None' for them which will
503 get all. You can also pass -1 for stop which will go to the end.
504 Returns the particles
507 kwds[
'resolution'] = resolution
508 if isinstance(tuple_selection, str):
509 kwds[
'molecule'] = tuple_selection
510 elif isinstance(tuple_selection, tuple):
511 rbegin = tuple_selection[0]
512 rend = tuple_selection[1]
513 kwds[
'molecule'] = tuple_selection[2]
515 copynum = tuple_selection[3]
516 if copynum
is not None:
517 kwds[
'copy_index'] = copynum
521 statenum = tuple_selection[4]
522 if statenum
is not None:
523 kwds[
'state_index'] = statenum
530 residue_indexes=range(1, rbegin),
532 return s.get_selected_particles()
534 kwds[
'residue_indexes'] = range(rbegin, rend+1)
536 return s.get_selected_particles()
539 def get_db_from_csv(csvfilename):
542 with open(csvfilename)
as fh:
543 csvr = csv.DictReader(fh)
545 outputlist.append(ls)
550 '''Return the component name provided a particle and a list of names'''
552 protname = root.get_name()
554 while protname
not in list_of_names:
555 root0 = root.get_parent()
558 protname = root0.get_name()
563 if "Beads" in protname:
566 return (protname, is_a_bead)
571 Retrieve the residue indexes for the given particle.
573 The particle must be an instance of Fragment,Residue, Atom or Molecule
574 or else returns an empty list
585 resind_tmp = IMP.pmi.tools.OrderedSet()
592 resind = list(resind_tmp)
598 def sort_by_residues(particles):
601 sorted_particles_residues = sorted(
603 key=
lambda tup: tup[1])
604 particles = [p[0]
for p
in sorted_particles_residues]
613 """Synchronize data over a parallel run"""
614 from mpi4py
import MPI
615 comm = MPI.COMM_WORLD
616 rank = comm.Get_rank()
617 number_of_processes = comm.size
620 comm.send(data, dest=0, tag=11)
623 for i
in range(1, number_of_processes):
624 data_tmp = comm.recv(source=i, tag=11)
625 if isinstance(data, list):
627 elif isinstance(data, dict):
628 data.update(data_tmp)
630 raise TypeError(
"data not supported, use list or dictionaries")
632 for i
in range(1, number_of_processes):
633 comm.send(data, dest=i, tag=11)
636 data = comm.recv(source=0, tag=11)
646 Yield all sublists of length >= lmin and <= lmax
652 for j
in range(i + lmin, min(n + 1, i + 1 + lmax)):
656 def flatten_list(ls):
657 return [item
for sublist
in ls
for item
in sublist]
661 """ Yield successive length-sized chunks from a list.
663 for i
in range(0, len(list), length):
664 yield list[i:i + length]
667 def chunk_list_into_segments(seq, num):
669 avg = len(seq) / float(num)
673 while last < len(seq):
674 out.append(seq[int(last):int(last + avg)])
682 ''' This class stores integers
683 in ordered compact lists eg:
685 the methods help splitting and merging the internal lists
687 s=Segments([1,2,3]) is [[1,2,3]]
688 s.add(4) is [[1,2,3,4]] (add right)
689 s.add(3) is [[1,2,3,4]] (item already existing)
690 s.add(7) is [[1,2,3,4],[7]] (new list)
691 s.add([8,9]) is [[1,2,3,4],[7,8,9]] (add item right)
692 s.add([5,6]) is [[1,2,3,4,5,6,7,8,9]] (merge)
693 s.remove(3) is [[1,2],[4,5,6,7,8,9]] (split)
698 '''index can be a integer or a list of integers '''
699 if isinstance(index, int):
700 self.segs = [[index]]
701 elif isinstance(index, list):
702 self.segs = [[index[0]]]
706 raise TypeError(
"index must be an int or list of ints")
709 '''index can be a integer or a list of integers '''
710 if isinstance(index, (int, numpy.int32, numpy.int64)):
713 for n, s
in enumerate(self.segs):
721 if mergeright
is None and mergeleft
is None:
722 self.segs.append([index])
723 if mergeright
is not None and mergeleft
is None:
724 self.segs[mergeright].append(index)
725 if mergeleft
is not None and mergeright
is None:
726 self.segs[mergeleft] = [index]+self.segs[mergeleft]
727 if mergeleft
is not None and mergeright
is not None:
728 self.segs[mergeright] = \
729 self.segs[mergeright]+[index]+self.segs[mergeleft]
730 del self.segs[mergeleft]
732 for n
in range(len(self.segs)):
735 self.segs.sort(key=
lambda tup: tup[0])
737 elif isinstance(index, list):
741 raise TypeError(
"index must be an int or list of ints")
744 '''index can be a integer'''
745 for n, s
in enumerate(self.segs):
750 self.segs[n] = s[:-1]
752 i = self.segs[n].index(index)
754 self.segs.append(s[i+1:])
755 for n
in range(len(self.segs)):
757 if len(self.segs[n]) == 0:
759 self.segs.sort(key=
lambda tup: tup[0])
762 ''' Returns a flatten list '''
763 return [item
for sublist
in self.segs
for item
in sublist]
767 for seg
in self.segs:
768 ret_tmp += str(seg[0])+
"-"+str(seg[-1])+
","
769 ret = ret_tmp[:-1]+
"]"
777 def normal_density_function(expected_value, sigma, x):
779 1 / math.sqrt(2 * math.pi) / sigma *
780 math.exp(-(x - expected_value) ** 2 / 2 / sigma / sigma)
784 def log_normal_density_function(expected_value, sigma, x):
786 1 / math.sqrt(2 * math.pi) / sigma / x *
787 math.exp(-(math.log(x / expected_value) ** 2 / 2 / sigma / sigma))
791 def print_multicolumn(list_of_strings, ncolumns=2, truncate=40):
797 for i
in range(len(ls) % cols):
800 split = [ls[i:i + len(ls) // cols]
801 for i
in range(0, len(ls), len(ls) // cols)]
802 for row
in zip(*split):
803 print(
"".join(str.ljust(i, truncate)
for i
in row))
807 '''Change color code to hexadecimal to rgb'''
809 self._NUMERALS =
'0123456789abcdefABCDEF'
810 self._HEXDEC = dict((v, int(v, 16))
for v
in
811 (x+y
for x
in self._NUMERALS
812 for y
in self._NUMERALS))
813 self.LOWERCASE, self.UPPERCASE =
'x',
'X'
815 def rgb(self, triplet):
816 return (float(self._HEXDEC[triplet[0:2]]),
817 float(self._HEXDEC[triplet[2:4]]),
818 float(self._HEXDEC[triplet[4:6]]))
820 def triplet(self, rgb, lettercase=None):
821 if lettercase
is None:
822 lettercase = self.LOWERCASE
823 return format(rgb[0] << 16 | rgb[1] << 8 | rgb[2],
'06'+lettercase)
827 class OrderedSet(MutableSet):
829 def __init__(self, iterable=None):
831 end += [
None, end, end]
833 if iterable
is not None:
839 def __contains__(self, key):
840 return key
in self.map
843 if key
not in self.map:
846 curr[2] = end[1] = self.map[key] = [key, curr, end]
848 def discard(self, key):
850 key, prev, next = self.map.pop(key)
857 while curr
is not end:
861 def __reversed__(self):
864 while curr
is not end:
868 def pop(self, last=True):
870 raise KeyError(
'set is empty')
880 return '%s()' % (self.__class__.__name__,)
881 return '%s(%r)' % (self.__class__.__name__, list(self))
883 def __eq__(self, other):
884 if isinstance(other, OrderedSet):
885 return len(self) == len(other)
and list(self) == list(other)
886 return set(self) == set(other)
890 """Store objects in order they were added, but with default type.
891 Source: http://stackoverflow.com/a/4127426/2608793
893 def __init__(self, *args, **kwargs):
895 self.default_factory =
None
897 if not (args[0]
is None or callable(args[0])):
898 raise TypeError(
'first argument must be callable or None')
899 self.default_factory = args[0]
901 super(OrderedDefaultDict, self).__init__(*args, **kwargs)
903 def __missing__(self, key):
904 if self.default_factory
is None:
906 self[key] = default = self.default_factory()
909 def __reduce__(self):
910 args = (self.default_factory,)
if self.default_factory
else ()
911 if sys.version_info[0] >= 3:
912 return self.__class__, args,
None,
None, self.items()
914 return self.__class__, args,
None,
None, self.iteritems()
920 """Extract frame from RMF file and fill coordinates. Must be identical
923 @param hier The (System) hierarchy to fill (e.g. after you've built it)
924 @param rmf_fn The file to extract from
925 @param frame_num The frame number to extract
927 rh = RMF.open_rmf_file_read_only(rmf_fn)
933 def input_adaptor(stuff, pmi_resolution=0, flatten=False, selection_tuple=None,
934 warn_about_slices=
True):
935 """Adapt things for PMI (degrees of freedom, restraints, ...)
936 Returns list of list of hierarchies, separated into Molecules if possible.
937 The input can be a list, or a list of lists (iterable of ^1 or
939 (iterable of ^2) Hierarchy -> returns input as list of list of hierarchies,
940 only one entry, not grouped by molecules.
941 (iterable of ^2) PMI::System/State/Molecule/TempResidue ->
942 returns residue hierarchies, grouped in molecules, at requested
945 @param stuff Can be one of the following inputs:
946 IMP Hierarchy, PMI System/State/Molecule/TempResidue, or a
947 list/set (of list/set) of them.
948 Must be uniform input, however. No mixing object types.
949 @param pmi_resolution For selecting, only does it if you pass PMI
950 objects. Set it to "all" if you want all resolutions!
951 @param flatten Set to True if you just want all hierarchies in one list.
952 @param warn_about_slices Print a warning if you are requesting only part
953 of a bead. Sometimes you just don't care!
954 @note since this relies on IMP::atom::Selection, this will not return
955 any objects if they weren't built! But there should be no problem
956 if you request unbuilt residues - they should be ignored.
962 if hasattr(stuff,
'__iter__'):
965 thelist = list(stuff)
968 if all(hasattr(el,
'__iter__')
for el
in thelist):
969 thelist = [i
for sublist
in thelist
for i
in sublist]
970 elif any(hasattr(el,
'__iter__')
for el
in thelist):
971 raise Exception(
'input_adaptor: input_object must be a list '
972 'or a list of lists')
981 except (NotImplementedError, TypeError):
993 if is_system
or is_state
or is_molecule
or is_temp_residue:
1000 for state
in system.get_states():
1001 mdict = state.get_molecules()
1002 for molname
in mdict:
1003 for copy
in mdict[molname]:
1004 indexes_per_mol[copy] += \
1005 [r.get_index()
for r
in copy.get_residues()]
1008 mdict = state.get_molecules()
1009 for molname
in mdict:
1010 for copy
in mdict[molname]:
1011 indexes_per_mol[copy] += [r.get_index()
1012 for r
in copy.get_residues()]
1014 for molecule
in stuff:
1015 indexes_per_mol[molecule] += [r.get_index()
1016 for r
in molecule.get_residues()]
1018 for tempres
in stuff:
1019 indexes_per_mol[tempres.get_molecule()].append(
1020 tempres.get_index())
1021 for mol
in indexes_per_mol:
1022 if pmi_resolution ==
'all':
1026 mol.get_hierarchy(), residue_indexes=indexes_per_mol[mol])
1029 resolution=pmi_resolution,
1030 residue_indexes=indexes_per_mol[mol])
1031 ps = sel.get_selected_particles()
1034 if warn_about_slices:
1035 rset = set(indexes_per_mol[mol])
1039 if not fset <= rset:
1045 resbreak = maxf
if minf == minset
else minset-1
1047 'You are trying to select only part of the '
1048 'bead %s:%i-%i. The residues you requested '
1049 'are %i-%i. You can fix this by: '
1050 '1) requesting the whole bead/none of it; or'
1051 '2) break the bead up by passing '
1052 'bead_extra_breaks=[\'%i\'] in '
1053 'molecule.add_representation()'
1054 % (mol.get_name(), minset, maxset, minf, maxf,
1060 if pmi_resolution ==
'all':
1066 h, resolution=pmi_resolution).get_selected_particles()
1069 hier_list = [hier_list]
1071 raise Exception(
'input_adaptor: you passed something of wrong type '
1072 'or a list with mixed types')
1074 if flatten
and pmi_input:
1075 return [h
for sublist
in hier_list
for h
in sublist]
1081 """Returns sequence-sorted segments array, each containing the first
1082 particle the last particle and the first residue index."""
1084 from operator
import itemgetter
1087 raise ValueError(
"only pass stuff from one Molecule, please")
1102 segs.append((start, end, startres))
1103 return sorted(segs, key=itemgetter(2))
1107 """Decorate the sequence-consecutive particles from a PMI2 molecule
1108 with a bond, so that they appear connected in the rmf file"""
1110 for x
in range(len(SortedSegments) - 1):
1112 last = SortedSegments[x][1]
1113 first = SortedSegments[x + 1][0]
1115 p1 = last.get_particle()
1116 p2 = first.get_particle()
1129 """ Just get the leaves from a list of hierarchies """
1130 lvs = list(itertools.chain.from_iterable(
1136 """Perform selection using the usual keywords but return ALL
1137 resolutions (BEADS and GAUSSIANS).
1138 Returns in flat list!
1143 if hier
is not None:
1146 warnings.warn(
"You passed nothing to select_at_all_resolutions()",
1154 raise Exception(
'select_at_all_resolutions: you have to pass '
1157 raise Exception(
'select_at_all_resolutions: you have to pass '
1159 if 'resolution' in kwargs
or 'representation_type' in kwargs:
1160 raise Exception(
"don't pass resolution or representation_type "
1163 representation_type=IMP.atom.BALLS,
1166 representation_type=IMP.atom.DENSITIES,
1168 ret |= OrderedSet(selB.get_selected_particles())
1169 ret |= OrderedSet(selD.get_selected_particles())
1178 """Utility to retrieve particles from a hierarchy within a
1179 zone around a set of ps.
1180 @param hier The hierarchy in which to look for neighbors
1181 @param target_ps The particles for zoning
1182 @param sel_zone The maximum distance
1183 @param entire_residues If True, will grab entire residues
1184 @param exclude_backbone If True, will only return sidechain particles
1188 backbone_types = [
'C',
'N',
'CB',
'O']
1189 if exclude_backbone:
1192 test_ps = test_sel.get_selected_particles()
1193 nn = IMP.algebra.NearestNeighbor3D([
IMP.core.XYZ(p).get_coordinates()
1196 for target
in target_ps:
1197 zone |= set(nn.get_in_ball(
IMP.core.XYZ(target).get_coordinates(),
1199 zone_ps = [test_ps[z]
for z
in zone]
1204 zone_ps = [h.get_particle()
for h
in final_ps]
1209 """Returns unique objects in original order"""
1213 if not hasattr(hiers,
'__iter__'):
1220 rbs_ordered.append(rb)
1225 rbs_ordered.append(rb)
1229 return rbs_ordered, beads
1233 "This function returns the parent molecule hierarchies of given objects"
1234 stuff =
input_adaptor(input_objects, pmi_resolution=
'all', flatten=
True)
1239 while not (is_molecule
or is_root):
1240 root = IMP.atom.get_root(h)
1247 return list(molecules)
1250 def get_molecules_dictionary(input_objects):
1251 moldict = defaultdict(list)
1253 name = mol.get_name()
1254 moldict[name].append(mol)
1261 def get_molecules_dictionary_by_copy(input_objects):
1262 moldict = defaultdict(dict)
1264 name = mol.get_name()
1266 moldict[name][c] = mol
1270 def get_selections_dictionary(input_objects):
1271 moldict = IMP.pmi.tools.get_molecules_dictionary(input_objects)
1272 seldict = defaultdict(list)
1273 for name, mols
in moldict.items():
1280 """Given a list of PMI objects, returns all density hierarchies within
1281 these objects. The output of this function can be inputted into
1282 things such as EM restraints. This function is intended to gather
1283 density particles appended to molecules (and not other hierarchies
1284 which might have been appended to the root node directly).
1293 i, representation_type=IMP.atom.DENSITIES).get_selected_particles()
1298 max_translation=300., max_rotation=2.0 * math.pi,
1299 avoidcollision_rb=
True, avoidcollision_fb=
False,
1300 cutoff=10.0, niterations=100,
1302 excluded_rigid_bodies=[],
1303 hierarchies_excluded_from_collision=[],
1304 hierarchies_included_in_collision=[],
1306 return_debug=
False):
1307 """Shuffle particles. Used to restart the optimization.
1308 The configuration of the system is initialized by placing each
1309 rigid body and each bead randomly in a box. If `bounding_box` is
1310 specified, the particles are placed inside this box; otherwise, each
1311 particle is displaced by up to max_translation angstroms, and randomly
1312 rotated. Effort is made to place particles far enough from each other to
1313 prevent any steric clashes.
1314 @param objects Can be one of the following inputs:
1315 IMP Hierarchy, PMI System/State/Molecule/TempResidue, or
1317 @param max_translation Max translation (rbs and flexible beads)
1318 @param max_rotation Max rotation (rbs only)
1319 @param avoidcollision_rb check if the particle/rigid body was
1320 placed close to another particle; uses the optional
1321 arguments cutoff and niterations
1322 @param avoidcollision_fb Advanced. Generally you want this False because
1323 it's hard to shuffle beads.
1324 @param cutoff Distance less than this is a collision
1325 @param niterations How many times to try avoiding collision
1326 @param bounding_box Only shuffle particles within this box.
1327 Defined by ((x1,y1,z1),(x2,y2,z2)).
1328 @param excluded_rigid_bodies Don't shuffle these rigid body objects
1329 @param hierarchies_excluded_from_collision Don't count collision
1331 @param hierarchies_included_in_collision Hierarchies that are not
1332 shuffled, but should be included in collision calculation
1334 @param verbose Give more output
1335 @note Best to only call this function after you've set up degrees
1337 For debugging purposes, returns: <shuffled indexes>,
1338 <collision avoided indexes>
1343 pmi_resolution=
'all',
1346 if len(rigid_bodies) > 0:
1347 mdl = rigid_bodies[0].get_model()
1348 elif len(flexible_beads) > 0:
1349 mdl = flexible_beads[0].get_model()
1351 raise Exception(
"Could not find any particles in the hierarchy")
1352 if len(rigid_bodies) == 0:
1353 print(
"shuffle_configuration: rigid bodies were not initialized")
1357 gcpf.set_distance(cutoff)
1361 hierarchies_excluded_from_collision, pmi_resolution=
'all',
1365 hierarchies_included_in_collision, pmi_resolution=
'all', flatten=
True)
1367 collision_excluded_idxs = set(
1369 for h
in collision_excluded_hierarchies
1372 collision_included_idxs = set(
1374 for h
in collision_included_hierarchies
1381 all_idxs.append(p.get_particle_index())
1383 collision_excluded_idxs.add(p.get_particle_index())
1385 if bounding_box
is not None:
1386 ((x1, y1, z1), (x2, y2, z2)) = bounding_box
1391 all_idxs = set(all_idxs) | collision_included_idxs
1392 all_idxs = all_idxs - collision_excluded_idxs
1394 print(
'shuffling', len(rigid_bodies),
'rigid bodies')
1395 for rb
in rigid_bodies:
1396 if rb
not in excluded_rigid_bodies:
1398 if avoidcollision_rb:
1399 rb_idxs = set(rb.get_member_particle_indexes()) - \
1400 collision_excluded_idxs
1401 other_idxs = all_idxs - rb_idxs
1405 while niter < niterations:
1406 rbxyz = (rb.get_x(), rb.get_y(), rb.get_z())
1423 rbxyz, max_translation, max_rotation)
1425 debug.append([rb, other_idxs
if avoidcollision_rb
else set()])
1429 if avoidcollision_rb
and other_idxs:
1431 npairs = len(gcpf.get_close_pairs(mdl,
1439 print(
"shuffle_configuration: rigid body placed "
1440 "close to other %d particles, trying "
1441 "again..." % npairs)
1442 print(
"shuffle_configuration: rigid body name: "
1444 if niter == niterations:
1446 "tried the maximum number of iterations to "
1447 "avoid collisions, increase the distance "
1452 print(
'shuffling', len(flexible_beads),
'flexible beads')
1453 for fb
in flexible_beads:
1455 if avoidcollision_fb:
1457 other_idxs = all_idxs - fb_idxs
1463 while niter < niterations:
1470 fbxyz, max_translation, max_rotation)
1475 xyz = memb.get_internal_coordinates()
1480 rf = memb.get_rigid_body().get_reference_frame()
1481 glob_to_int = rf.get_transformation_from()
1482 memb.set_internal_coordinates(
1483 glob_to_int.get_transformed(translation))
1485 xyz_transformed = transformation.get_transformed(xyz)
1486 memb.set_internal_coordinates(xyz_transformed)
1487 debug.append([xyz, other_idxs
if avoidcollision_fb
else set()])
1495 -d.get_coordinates())
1499 debug.append([d, other_idxs
if avoidcollision_fb
else set()])
1506 if avoidcollision_fb:
1508 npairs = len(gcpf.get_close_pairs(mdl,
1516 print(
"shuffle_configuration: floppy body placed close "
1517 "to other %d particles, trying again..." % npairs)
1518 if niter == niterations:
1520 "tried the maximum number of iterations to avoid "
1521 "collisions, increase the distance cutoff")
1528 class ColorHierarchy(object):
1530 def __init__(self, hier):
1531 import matplotlib
as mpl
1533 import matplotlib.pyplot
as plt
1537 hier.ColorHierarchy = self
1542 self.method = self.nochange
1550 def get_color(self, fl):
1553 def get_log_scale(self, fl):
1556 return math.log(fl+eps)
1558 def color_by_resid(self):
1559 self.method = self.color_by_resid
1560 self.scheme = self.mpl.cm.rainbow
1561 for mol
in self.mols:
1568 c = self.get_color(float(ri)/self.last)
1572 avr = sum(ris)/len(ris)
1573 c = self.get_color(float(avr)/self.last)
1576 def color_by_uncertainty(self):
1577 self.method = self.color_by_uncertainty
1578 self.scheme = self.mpl.cm.jet
1585 self.first = self.get_log_scale(1.0)
1586 self.last = self.get_log_scale(100.0)
1588 value = self.get_log_scale(unc_dict[p])
1589 if value >= self.last:
1591 if value <= self.first:
1593 c = self.get_color((value-self.first) / (self.last-self.first))
1596 def get_color_bar(self, filename):
1597 import matplotlib
as mpl
1599 import matplotlib.pyplot
as plt
1601 fig = plt.figure(figsize=(8, 3))
1602 ax1 = fig.add_axes([0.05, 0.80, 0.9, 0.15])
1605 norm = mpl.colors.Normalize(vmin=0.0, vmax=1.0)
1607 if self.method == self.color_by_uncertainty:
1608 angticks = [1.0, 2.5, 5.0, 10.0, 25.0, 50.0, 100.0]
1612 vvalue = (self.get_log_scale(at)-self.first) \
1613 / (self.last-self.first)
1614 if vvalue <= 1.0
and vvalue >= 0.0:
1615 vvalues.append(vvalue)
1616 marks.append(str(at))
1617 cb1 = mpl.colorbar.ColorbarBase(
1618 ax1, cmap=cmap, norm=norm, ticks=vvalues,
1619 orientation=
'horizontal')
1620 print(self.first, self.last, marks, vvalues)
1621 cb1.ax.set_xticklabels(marks)
1622 cb1.set_label(
'Angstorm')
1623 plt.savefig(filename, dpi=150, transparent=
True)
1628 """Given a Chimera color name or hex color value, return RGB"""
1629 d = {
'aquamarine': (0.4980392156862745, 1.0, 0.8313725490196079),
1630 'black': (0.0, 0.0, 0.0),
1631 'blue': (0.0, 0.0, 1.0),
1632 'brown': (0.6470588235, 0.16470588235294117, 0.16470588235294117),
1633 'chartreuse': (0.4980392156862745, 1.0, 0.0),
1634 'coral': (1.0, 0.4980392156862745, 0.3137254901960784),
1635 'cornflower blue': (0.39215686, 0.58431372549, 0.9294117647058824),
1636 'cyan': (0.0, 1.0, 1.0),
1637 'dark cyan': (0.0, 0.5450980392156862, 0.5450980392156862),
1638 'dark gray': (0.6627450980, 0.6627450980392157, 0.6627450980392157),
1639 'dark green': (0.0, 0.39215686274509803, 0.0),
1640 'dark khaki': (0.74117647, 0.7176470588235294, 0.4196078431372549),
1641 'dark magenta': (0.5450980392156862, 0.0, 0.5450980392156862),
1642 'dark olive green': (0.333333333, 0.419607843, 0.1843137254901961),
1643 'dark red': (0.5450980392156862, 0.0, 0.0),
1644 'dark slate blue': (0.28235294, 0.239215686, 0.5450980392156862),
1645 'dark slate gray': (0.1843137, 0.30980392, 0.30980392156862746),
1646 'deep pink': (1.0, 0.0784313725490196, 0.5764705882352941),
1647 'deep sky blue': (0.0, 0.7490196078431373, 1.0),
1648 'dim gray': (0.41176470, 0.4117647058823529, 0.4117647058823529),
1649 'dodger blue': (0.11764705882352941, 0.5647058823529412, 1.0),
1650 'firebrick': (0.6980392, 0.13333333333333333, 0.13333333333333333),
1651 'forest green': (0.13333333, 0.5450980392156862, 0.13333333333333333),
1652 'gold': (1.0, 0.8431372549019608, 0.0),
1653 'goldenrod': (0.85490196, 0.6470588235294118, 0.12549019607843137),
1654 'gray': (0.7450980392156863, 0.7450980392156863, 0.7450980392156863),
1655 'green': (0.0, 1.0, 0.0),
1656 'hot pink': (1.0, 0.4117647058823529, 0.7058823529411765),
1657 'khaki': (0.9411764705882353, 0.9019607843137255, 0.5490196078431373),
1658 'light blue': (0.67843137, 0.8470588235294118, 0.9019607843137255),
1659 'light gray': (0.82745098, 0.8274509803921568, 0.8274509803921568),
1660 'light green': (0.56470588, 0.9333333333333333, 0.5647058823529412),
1661 'light sea green': (0.125490, 0.6980392156862745, 0.6666666666666666),
1662 'lime green': (0.1960784, 0.803921568627451, 0.19607843137254902),
1663 'magenta': (1.0, 0.0, 1.0),
1664 'medium blue': (0.1960784, 0.19607843137254902, 0.803921568627451),
1665 'medium purple': (0.57647, 0.4392156862745098, 0.8588235294117647),
1666 'navy blue': (0.0, 0.0, 0.5019607843137255),
1667 'olive drab': (0.4196078, 0.5568627450980392, 0.13725490196078433),
1668 'orange red': (1.0, 0.27058823529411763, 0.0),
1669 'orange': (1.0, 0.4980392156862745, 0.0),
1670 'orchid': (0.85490196, 0.4392156862745098, 0.8392156862745098),
1671 'pink': (1.0, 0.7529411764705882, 0.796078431372549),
1672 'plum': (0.8666666666666667, 0.6274509803921569, 0.8666666666666667),
1673 'purple': (0.62745098, 0.12549019607843137, 0.9411764705882353),
1674 'red': (1.0, 0.0, 0.0),
1675 'rosy brown': (0.7372549, 0.5607843137254902, 0.5607843137254902),
1676 'salmon': (0.980392, 0.5019607843137255, 0.4470588235294118),
1677 'sandy brown': (0.956862745, 0.6431372549019608, 0.3764705882352941),
1678 'sea green': (0.18039, 0.5450980392156862, 0.3411764705882353),
1679 'sienna': (0.6274509, 0.3215686274509804, 0.17647058823529413),
1680 'sky blue': (0.52941176, 0.807843137254902, 0.9215686274509803),
1681 'slate gray': (0.439215686, 0.50196078, 0.5647058823529412),
1682 'spring green': (0.0, 1.0, 0.4980392156862745),
1683 'steel blue': (0.2745098, 0.50980392, 0.70588235),
1684 'tan': (0.8235294117647058, 0.7058823529411765, 0.5490196078431373),
1685 'turquoise': (0.25098039, 0.87843137, 0.81568627),
1686 'violet red': (0.81568627, 0.125490196, 0.56470588235),
1687 'white': (1.0, 1.0, 1.0),
1688 'yellow': (1.0, 1.0, 0.0)}
1689 if colorname.startswith(
'#'):
1690 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.
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.
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 ...