3 """@namespace IMP.pmi.tools
4 Miscellaneous utilities.
7 from __future__
import print_function, division
19 from time
import process_time
21 from time
import clock
as process_time
24 from collections
import defaultdict, OrderedDict
28 def _get_system_for_hier(hier):
29 """Given a hierarchy, return the System that created it, or None"""
32 if hier
and not hasattr(hier,
'get_parent'):
39 if hasattr(hier,
'_pmi2_system'):
40 return hier._pmi2_system()
43 for ws
in IMP.pmi.topology.System._all_systems:
45 if s
and s.hier == hier:
48 hier = hier.get_parent()
51 def _all_protocol_outputs(hier):
52 """Iterate over all (ProtocolOutput, State) pairs for the
54 system = _get_system_for_hier(hier)
56 for state
in system.states:
57 for p
in state._protocol_output:
61 def _add_pmi_provenance(p):
62 """Tag the given particle as being created by the current version
67 location=
"https://integrativemodeling.org")
71 def _get_restraint_set_keys():
72 if not hasattr(_get_restraint_set_keys,
'pmi_rs_key'):
73 _get_restraint_set_keys.pmi_rs_key =
IMP.ModelKey(
"PMI restraints")
74 _get_restraint_set_keys.rmf_rs_key =
IMP.ModelKey(
"RMF restraints")
75 return (_get_restraint_set_keys.pmi_rs_key,
76 _get_restraint_set_keys.rmf_rs_key)
79 def _add_restraint_sets(model, mk, mk_rmf):
82 model.add_data(mk, rs)
83 model.add_data(mk_rmf, rs_rmf)
88 """Add a PMI restraint to the model.
89 Since Model.add_restraint() no longer exists (in modern IMP restraints
90 should be added to a ScoringFunction instead) store them instead in
91 a RestraintSet, and keep a reference to it in the Model.
93 If `add_to_rmf` is True, also add the restraint to a separate list
94 of restraints that will be written out to RMF files (by default, most
95 PMI restraints are not)."""
96 mk, mk_rmf = _get_restraint_set_keys()
97 if model.get_has_data(mk):
98 rs = IMP.RestraintSet.get_from(model.get_data(mk))
99 rs_rmf = IMP.RestraintSet.get_from(model.get_data(mk_rmf))
101 rs, rs_rmf = _add_restraint_sets(model, mk, mk_rmf)
102 rs.add_restraint(restraint)
104 rs_rmf.add_restraint(restraint)
108 """Get a RestraintSet containing all PMI restraints added to the model.
109 If `rmf` is True, return only the subset of these restraints that
110 should be written out to RMF files."""
111 mk, mk_rmf = _get_restraint_set_keys()
112 if not model.get_has_data(mk):
113 warnings.warn(
"no restraints added to model yet",
115 _add_restraint_sets(model, mk, mk_rmf)
117 return IMP.RestraintSet.get_from(model.get_data(mk_rmf))
119 return IMP.RestraintSet.get_from(model.get_data(mk))
123 """Collect timing information.
124 Add an instance of this class to outputobjects to get timing information
129 @param isdelta if True (the default) then report the time since the
130 last use of this class; if False, report cumulative time."""
131 self.starttime = process_time()
133 self.isdelta = isdelta
135 def set_label(self, labelstr):
136 self.label = labelstr
138 def get_output(self):
141 newtime = process_time()
142 output[
"Stopwatch_" + self.label +
"_delta_seconds"] \
143 = str(newtime - self.starttime)
144 self.starttime = newtime
146 output[
"Stopwatch_" + self.label +
"_elapsed_seconds"] \
147 = str(process_time() - self.starttime)
151 class SetupNuisance(object):
153 def __init__(self, m, initialvalue, minvalue, maxvalue, isoptimized=True,
161 nuisance.set_lower(minvalue)
163 nuisance.set_upper(maxvalue)
166 nuisance.set_is_optimized(nuisance.get_nuisance_key(), isoptimized)
167 self.nuisance = nuisance
169 def get_particle(self):
173 class SetupWeight(object):
175 def __init__(self, m, isoptimized=True, nweights_or_weights=None):
177 if isinstance(nweights_or_weights, int):
179 pw, nweights_or_weights
183 nweights_or_weights = list(nweights_or_weights)
185 pw, nweights_or_weights
189 self.weight.set_weights_are_optimized(isoptimized)
191 def get_particle(self):
195 class SetupSurface(object):
197 def __init__(self, m, center, normal, isoptimized=True):
200 self.surface.set_coordinates_are_optimized(isoptimized)
201 self.surface.set_normal_is_optimized(isoptimized)
203 def get_particle(self):
207 class ParticleToSampleList(object):
211 self.dictionary_particle_type = {}
212 self.dictionary_particle_transformation = {}
213 self.dictionary_particle_name = {}
216 def add_particle(self, particle, particle_type, particle_transformation,
218 if particle_type
not in [
"Rigid_Bodies",
"Floppy_Bodies",
"Nuisances",
219 "X_coord",
"Weights",
"Surfaces"]:
220 raise TypeError(
"not the right particle type")
222 self.dictionary_particle_type[particle] = particle_type
223 if particle_type ==
"Rigid_Bodies":
224 if (isinstance(particle_transformation, tuple)
225 and len(particle_transformation) == 2
226 and all(isinstance(x, float)
227 for x
in particle_transformation)):
228 self.dictionary_particle_transformation[
229 particle] = particle_transformation
230 self.dictionary_particle_name[particle] = name
233 "ParticleToSampleList: not the right transformation "
234 "format for Rigid_Bodies, should be a tuple of floats")
235 elif particle_type ==
"Surfaces":
236 if (isinstance(particle_transformation, tuple)
237 and len(particle_transformation) == 3
238 and all(isinstance(x, float)
239 for x
in particle_transformation)):
240 self.dictionary_particle_transformation[
241 particle] = particle_transformation
242 self.dictionary_particle_name[particle] = name
245 "ParticleToSampleList: not the right transformation "
246 "format for Surfaces, should be a tuple of floats")
248 if isinstance(particle_transformation, float):
249 self.dictionary_particle_transformation[
250 particle] = particle_transformation
251 self.dictionary_particle_name[particle] = name
254 "ParticleToSampleList: not the right transformation "
255 "format, should be a float")
257 def get_particles_to_sample(self):
259 for particle
in self.dictionary_particle_type:
260 key = self.dictionary_particle_type[particle] + \
261 "ParticleToSampleList_" + \
262 self.dictionary_particle_name[particle] +
"_" + self.label
265 self.dictionary_particle_transformation[particle])
270 def get_cross_link_data(directory, filename, dist, omega, sigma,
271 don=
None, doff=
None, prior=0, type_of_profile=
"gofr"):
273 (distmin, distmax, ndist) = dist
274 (omegamin, omegamax, nomega) = omega
275 (sigmamin, sigmamax, nsigma) = sigma
278 with open(filen)
as xlpot:
279 dictionary = ast.literal_eval(xlpot.readline())
281 xpot = dictionary[directory][filename][
"distance"]
282 pot = dictionary[directory][filename][type_of_profile]
284 dist_grid =
get_grid(distmin, distmax, ndist,
False)
285 omega_grid = get_log_grid(omegamin, omegamax, nomega)
286 sigma_grid = get_log_grid(sigmamin, sigmamax, nsigma)
288 if don
is not None and doff
is not None:
308 def get_grid(gmin, gmax, ngrid, boundaries):
310 dx = (gmax - gmin) / float(ngrid)
311 for i
in range(0, ngrid + 1):
312 if(
not boundaries
and i == 0):
314 if(
not boundaries
and i == ngrid):
316 grid.append(gmin + float(i) * dx)
320 def get_log_grid(gmin, gmax, ngrid):
322 for i
in range(0, ngrid + 1):
323 grid.append(gmin * math.exp(float(i) / ngrid * math.log(gmax / gmin)))
329 example '"{ID_Score}" > 28 AND "{Sample}" ==
330 "%10_1%" OR ":Sample}" == "%10_2%" OR ":Sample}"
331 == "%10_3%" OR ":Sample}" == "%8_1%" OR ":Sample}" == "%8_2%"'
334 import pyparsing
as pp
336 operator = pp.Regex(
">=|<=|!=|>|<|==|in").setName(
"operator")
337 value = pp.QuotedString(
339 r"[+-]?\d+(:?\.\d*)?(:?[eE][+-]?\d+)?")
340 identifier = pp.Word(pp.alphas, pp.alphanums +
"_")
341 comparison_term = identifier | value
342 condition = pp.Group(comparison_term + operator + comparison_term)
344 expr = pp.operatorPrecedence(condition, [
345 (
"OR", 2, pp.opAssoc.LEFT, ),
346 (
"AND", 2, pp.opAssoc.LEFT, ),
349 parsedstring = str(expr.parseString(inputstring)) \
355 .replace(
"{",
"float(entry['") \
356 .replace(
"}",
"'])") \
357 .replace(
":",
"str(entry['") \
358 .replace(
"}",
"'])") \
359 .replace(
"AND",
"and") \
364 def open_file_or_inline_text(filename):
366 fl = open(filename,
"r")
368 fl = filename.split(
"\n")
372 def get_ids_from_fasta_file(fastafile):
374 with open(fastafile)
as ff:
377 ids.append(line[1:-1])
383 this function works with plain hierarchies, as read from the pdb,
384 no multi-scale hierarchies
391 atom_type=IMP.atom.AT_CA)
399 print(
"get_closest_residue_position: exiting while loop "
402 p = sel.get_selected_particles()
407 print(
"get_closest_residue_position: got NO residues for hierarchy "
408 "%s and residue %i" % (hier, resindex))
410 "get_closest_residue_position: got NO residues for hierarchy "
411 "%s and residue %i" % (hier, resindex))
414 "got multiple residues for hierarchy %s and residue %i; the list "
416 % (hier, resindex, str([pp.get_name()
for pp
in p])))
421 Return the residue index gaps and contiguous segments in the hierarchy.
423 @param hierarchy hierarchy to examine
424 @param start first residue index
425 @param end last residue index
427 @return A list of lists of the form
428 [[1,100,"cont"],[101,120,"gap"],[121,200,"cont"]]
431 for n, rindex
in enumerate(range(start, end + 1)):
433 atom_type=IMP.atom.AT_CA)
435 if len(sel.get_selected_particles()) == 0:
439 rindexcont = start - 1
440 if rindexgap == rindex - 1:
446 gaps.append([rindex, rindex,
"gap"])
452 rindexgap = start - 1
454 if rindexcont == rindex - 1:
461 gaps.append([rindex, rindex,
"cont"])
472 def set_map_element(self, xvalue, yvalue):
473 self.map[xvalue] = yvalue
475 def get_map_element(self, invalue):
476 if isinstance(invalue, float):
480 dist = (invalue - x) * (invalue - x)
489 return self.map[minx]
490 elif isinstance(invalue, str):
491 return self.map[invalue]
493 raise TypeError(
"wrong type for map")
497 """New tuple format: molname OR (start,stop,molname,copynum,statenum)
498 Copy and state are optional. Can also use 'None' for them which will
499 get all. You can also pass -1 for stop which will go to the end.
500 Returns the particles
503 kwds[
'resolution'] = resolution
504 if isinstance(tuple_selection, str):
505 kwds[
'molecule'] = tuple_selection
506 elif isinstance(tuple_selection, tuple):
507 rbegin = tuple_selection[0]
508 rend = tuple_selection[1]
509 kwds[
'molecule'] = tuple_selection[2]
511 copynum = tuple_selection[3]
512 if copynum
is not None:
513 kwds[
'copy_index'] = copynum
517 statenum = tuple_selection[4]
518 if statenum
is not None:
519 kwds[
'state_index'] = statenum
526 residue_indexes=range(1, rbegin),
528 return s.get_selected_particles()
530 kwds[
'residue_indexes'] = range(rbegin, rend+1)
532 return s.get_selected_particles()
535 def get_db_from_csv(csvfilename):
538 with open(csvfilename)
as fh:
539 csvr = csv.DictReader(fh)
541 outputlist.append(ls)
546 '''Return the component name provided a particle and a list of names'''
548 protname = root.get_name()
550 while protname
not in list_of_names:
551 root0 = root.get_parent()
554 protname = root0.get_name()
559 if "Beads" in protname:
562 return (protname, is_a_bead)
567 Retrieve the residue indexes for the given particle.
569 The particle must be an instance of Fragment,Residue, Atom or Molecule
570 or else returns an empty list
581 resind_tmp = IMP.pmi.tools.OrderedSet()
588 resind = list(resind_tmp)
594 def sort_by_residues(particles):
597 sorted_particles_residues = sorted(
599 key=
lambda tup: tup[1])
600 particles = [p[0]
for p
in sorted_particles_residues]
609 """Synchronize data over a parallel run"""
610 from mpi4py
import MPI
611 comm = MPI.COMM_WORLD
612 rank = comm.Get_rank()
613 number_of_processes = comm.size
616 comm.send(data, dest=0, tag=11)
619 for i
in range(1, number_of_processes):
620 data_tmp = comm.recv(source=i, tag=11)
621 if isinstance(data, list):
623 elif isinstance(data, dict):
624 data.update(data_tmp)
626 raise TypeError(
"data not supported, use list or dictionaries")
628 for i
in range(1, number_of_processes):
629 comm.send(data, dest=i, tag=11)
632 data = comm.recv(source=0, tag=11)
642 Yield all sublists of length >= lmin and <= lmax
650 for j
in range(i + 1, n):
651 if len(ls[i:j]) <= lmax
and len(ls[i:j]) >= lmin:
655 def flatten_list(ls):
656 return [item
for sublist
in ls
for item
in sublist]
660 """ Yield successive length-sized chunks from a list.
662 for i
in range(0, len(list), length):
663 yield list[i:i + length]
666 def chunk_list_into_segments(seq, num):
668 avg = len(seq) / float(num)
672 while last < len(seq):
673 out.append(seq[int(last):int(last + avg)])
681 ''' This class stores integers
682 in ordered compact lists eg:
684 the methods help splitting and merging the internal lists
686 s=Segments([1,2,3]) is [[1,2,3]]
687 s.add(4) is [[1,2,3,4]] (add right)
688 s.add(3) is [[1,2,3,4]] (item already existing)
689 s.add(7) is [[1,2,3,4],[7]] (new list)
690 s.add([8,9]) is [[1,2,3,4],[7,8,9]] (add item right)
691 s.add([5,6]) is [[1,2,3,4,5,6,7,8,9]] (merge)
692 s.remove(3) is [[1,2],[4,5,6,7,8,9]] (split)
697 '''index can be a integer or a list of integers '''
698 if isinstance(index, int):
699 self.segs = [[index]]
700 elif isinstance(index, list):
701 self.segs = [[index[0]]]
705 raise TypeError(
"index must be an int or list of ints")
708 '''index can be a integer or a list of integers '''
709 if isinstance(index, int):
712 for n, s
in enumerate(self.segs):
720 if mergeright
is None and mergeleft
is None:
721 self.segs.append([index])
722 if mergeright
is not None and mergeleft
is None:
723 self.segs[mergeright].append(index)
724 if mergeleft
is not None and mergeright
is None:
725 self.segs[mergeleft] = [index]+self.segs[mergeleft]
726 if mergeleft
is not None and mergeright
is not None:
727 self.segs[mergeright] = \
728 self.segs[mergeright]+[index]+self.segs[mergeleft]
729 del self.segs[mergeleft]
731 for n
in range(len(self.segs)):
734 self.segs.sort(key=
lambda tup: tup[0])
736 elif isinstance(index, list):
740 raise TypeError(
"index must be an int or list of ints")
743 '''index can be a integer'''
744 for n, s
in enumerate(self.segs):
749 self.segs[n] = s[:-1]
751 i = self.segs[n].index(index)
753 self.segs.append(s[i+1:])
754 for n
in range(len(self.segs)):
756 if len(self.segs[n]) == 0:
758 self.segs.sort(key=
lambda tup: tup[0])
761 ''' Returns a flatten list '''
762 return [item
for sublist
in self.segs
for item
in sublist]
766 for seg
in self.segs:
767 ret_tmp += str(seg[0])+
"-"+str(seg[-1])+
","
768 ret = ret_tmp[:-1]+
"]"
776 def normal_density_function(expected_value, sigma, x):
778 1 / math.sqrt(2 * math.pi) / sigma *
779 math.exp(-(x - expected_value) ** 2 / 2 / sigma / sigma)
783 def log_normal_density_function(expected_value, sigma, x):
785 1 / math.sqrt(2 * math.pi) / sigma / x *
786 math.exp(-(math.log(x / expected_value) ** 2 / 2 / sigma / sigma))
790 def print_multicolumn(list_of_strings, ncolumns=2, truncate=40):
796 for i
in range(len(ls) % cols):
799 split = [ls[i:i + len(ls) // cols]
800 for i
in range(0, len(ls), len(ls) // cols)]
801 for row
in zip(*split):
802 print(
"".join(str.ljust(i, truncate)
for i
in row))
806 '''Change color code to hexadecimal to rgb'''
808 self._NUMERALS =
'0123456789abcdefABCDEF'
809 self._HEXDEC = dict((v, int(v, 16))
for v
in
810 (x+y
for x
in self._NUMERALS
811 for y
in self._NUMERALS))
812 self.LOWERCASE, self.UPPERCASE =
'x',
'X'
814 def rgb(self, triplet):
815 return (float(self._HEXDEC[triplet[0:2]]),
816 float(self._HEXDEC[triplet[2:4]]),
817 float(self._HEXDEC[triplet[4:6]]))
819 def triplet(self, rgb, lettercase=None):
820 if lettercase
is None:
821 lettercase = self.LOWERCASE
822 return format(rgb[0] << 16 | rgb[1] << 8 | rgb[2],
'06'+lettercase)
826 class OrderedSet(collections.MutableSet):
828 def __init__(self, iterable=None):
830 end += [
None, end, end]
832 if iterable
is not None:
838 def __contains__(self, key):
839 return key
in self.map
842 if key
not in self.map:
845 curr[2] = end[1] = self.map[key] = [key, curr, end]
847 def discard(self, key):
849 key, prev, next = self.map.pop(key)
856 while curr
is not end:
860 def __reversed__(self):
863 while curr
is not end:
867 def pop(self, last=True):
869 raise KeyError(
'set is empty')
879 return '%s()' % (self.__class__.__name__,)
880 return '%s(%r)' % (self.__class__.__name__, list(self))
882 def __eq__(self, other):
883 if isinstance(other, OrderedSet):
884 return len(self) == len(other)
and list(self) == list(other)
885 return set(self) == set(other)
889 """Store objects in order they were added, but with default type.
890 Source: http://stackoverflow.com/a/4127426/2608793
892 def __init__(self, *args, **kwargs):
894 self.default_factory =
None
896 if not (args[0]
is None or callable(args[0])):
897 raise TypeError(
'first argument must be callable or None')
898 self.default_factory = args[0]
900 super(OrderedDefaultDict, self).__init__(*args, **kwargs)
902 def __missing__(self, key):
903 if self.default_factory
is None:
905 self[key] = default = self.default_factory()
908 def __reduce__(self):
909 args = (self.default_factory,)
if self.default_factory
else ()
910 if sys.version_info[0] >= 3:
911 return self.__class__, args,
None,
None, self.items()
913 return self.__class__, args,
None,
None, self.iteritems()
919 """Extract frame from RMF file and fill coordinates. Must be identical
922 @param hier The (System) hierarchy to fill (e.g. after you've built it)
923 @param rmf_fn The file to extract from
924 @param frame_num The frame number to extract
926 rh = RMF.open_rmf_file_read_only(rmf_fn)
932 def input_adaptor(stuff, pmi_resolution=0, flatten=False, selection_tuple=None,
933 warn_about_slices=
True):
934 """Adapt things for PMI (degrees of freedom, restraints, ...)
935 Returns list of list of hierarchies, separated into Molecules if possible.
936 The input can be a list, or a list of lists (iterable of ^1 or
938 (iterable of ^2) Hierarchy -> returns input as list of list of hierarchies,
939 only one entry, not grouped by molecules.
940 (iterable of ^2) PMI::System/State/Molecule/TempResidue ->
941 returns residue hierarchies, grouped in molecules, at requested
944 @param stuff Can be one of the following inputs:
945 IMP Hierarchy, PMI System/State/Molecule/TempResidue, or a
946 list/set (of list/set) of them.
947 Must be uniform input, however. No mixing object types.
948 @param pmi_resolution For selecting, only does it if you pass PMI
949 objects. Set it to "all" if you want all resolutions!
950 @param flatten Set to True if you just want all hierarchies in one list.
951 @param warn_about_slices Print a warning if you are requesting only part
952 of a bead. Sometimes you just don't care!
953 @note since this relies on IMP::atom::Selection, this will not return
954 any objects if they weren't built! But there should be no problem
955 if you request unbuilt residues - they should be ignored.
961 if hasattr(stuff,
'__iter__'):
964 thelist = list(stuff)
967 if all(hasattr(el,
'__iter__')
for el
in thelist):
968 thelist = [i
for sublist
in thelist
for i
in sublist]
969 elif any(hasattr(el,
'__iter__')
for el
in thelist):
970 raise Exception(
'input_adaptor: input_object must be a list '
971 'or a list of lists')
980 except (NotImplementedError, TypeError):
992 if is_system
or is_state
or is_molecule
or is_temp_residue:
999 for state
in system.get_states():
1000 mdict = state.get_molecules()
1001 for molname
in mdict:
1002 for copy
in mdict[molname]:
1003 indexes_per_mol[copy] += \
1004 [r.get_index()
for r
in copy.get_residues()]
1007 mdict = state.get_molecules()
1008 for molname
in mdict:
1009 for copy
in mdict[molname]:
1010 indexes_per_mol[copy] += [r.get_index()
1011 for r
in copy.get_residues()]
1013 for molecule
in stuff:
1014 indexes_per_mol[molecule] += [r.get_index()
1015 for r
in molecule.get_residues()]
1017 for tempres
in stuff:
1018 indexes_per_mol[tempres.get_molecule()].append(
1019 tempres.get_index())
1020 for mol
in indexes_per_mol:
1021 if pmi_resolution ==
'all':
1025 mol.get_hierarchy(), residue_indexes=indexes_per_mol[mol])
1028 resolution=pmi_resolution,
1029 residue_indexes=indexes_per_mol[mol])
1030 ps = sel.get_selected_particles()
1033 if warn_about_slices:
1034 rset = set(indexes_per_mol[mol])
1038 if not fset <= rset:
1044 resbreak = maxf
if minf == minset
else minset-1
1046 'You are trying to select only part of the '
1047 'bead %s:%i-%i. The residues you requested '
1048 'are %i-%i. You can fix this by: '
1049 '1) requesting the whole bead/none of it; or'
1050 '2) break the bead up by passing '
1051 'bead_extra_breaks=[\'%i\'] in '
1052 'molecule.add_representation()'
1053 % (mol.get_name(), minset, maxset, minf, maxf,
1059 if pmi_resolution ==
'all':
1065 h, resolution=pmi_resolution).get_selected_particles()
1068 hier_list = [hier_list]
1070 raise Exception(
'input_adaptor: you passed something of wrong type '
1071 'or a list with mixed types')
1073 if flatten
and pmi_input:
1074 return [h
for sublist
in hier_list
for h
in sublist]
1080 """Returns sequence-sorted segments array, each containing the first
1081 particle the last particle and the first residue index."""
1083 from operator
import itemgetter
1086 raise ValueError(
"only pass stuff from one Molecule, please")
1101 segs.append((start, end, startres))
1102 return sorted(segs, key=itemgetter(2))
1106 """Decorate the sequence-consecutive particles from a PMI2 molecule
1107 with a bond, so that they appear connected in the rmf file"""
1109 for x
in range(len(SortedSegments) - 1):
1111 last = SortedSegments[x][1]
1112 first = SortedSegments[x + 1][0]
1114 p1 = last.get_particle()
1115 p2 = first.get_particle()
1128 """ Just get the leaves from a list of hierarchies """
1129 lvs = list(itertools.chain.from_iterable(
1135 """Perform selection using the usual keywords but return ALL
1136 resolutions (BEADS and GAUSSIANS).
1137 Returns in flat list!
1142 if hier
is not None:
1145 warnings.warn(
"You passed nothing to select_at_all_resolutions()",
1153 raise Exception(
'select_at_all_resolutions: you have to pass '
1156 raise Exception(
'select_at_all_resolutions: you have to pass '
1158 if 'resolution' in kwargs
or 'representation_type' in kwargs:
1159 raise Exception(
"don't pass resolution or representation_type "
1162 representation_type=IMP.atom.BALLS,
1165 representation_type=IMP.atom.DENSITIES,
1167 ret |= OrderedSet(selB.get_selected_particles())
1168 ret |= OrderedSet(selD.get_selected_particles())
1177 """Utility to retrieve particles from a hierarchy within a
1178 zone around a set of ps.
1179 @param hier The hierarchy in which to look for neighbors
1180 @param target_ps The particles for zoning
1181 @param sel_zone The maximum distance
1182 @param entire_residues If True, will grab entire residues
1183 @param exclude_backbone If True, will only return sidechain particles
1187 backbone_types = [
'C',
'N',
'CB',
'O']
1188 if exclude_backbone:
1191 test_ps = test_sel.get_selected_particles()
1192 nn = IMP.algebra.NearestNeighbor3D([
IMP.core.XYZ(p).get_coordinates()
1195 for target
in target_ps:
1196 zone |= set(nn.get_in_ball(
IMP.core.XYZ(target).get_coordinates(),
1198 zone_ps = [test_ps[z]
for z
in zone]
1203 zone_ps = [h.get_particle()
for h
in final_ps]
1208 """Returns unique objects in original order"""
1212 if not hasattr(hiers,
'__iter__'):
1219 rbs_ordered.append(rb)
1224 rbs_ordered.append(rb)
1228 return rbs_ordered, beads
1232 "This function returns the parent molecule hierarchies of given objects"
1233 stuff =
input_adaptor(input_objects, pmi_resolution=
'all', flatten=
True)
1238 while not (is_molecule
or is_root):
1239 root = IMP.atom.get_root(h)
1246 return list(molecules)
1249 def get_molecules_dictionary(input_objects):
1250 moldict = defaultdict(list)
1252 name = mol.get_name()
1253 moldict[name].append(mol)
1260 def get_molecules_dictionary_by_copy(input_objects):
1261 moldict = defaultdict(dict)
1263 name = mol.get_name()
1265 moldict[name][c] = mol
1269 def get_selections_dictionary(input_objects):
1270 moldict = IMP.pmi.tools.get_molecules_dictionary(input_objects)
1271 seldict = defaultdict(list)
1272 for name, mols
in moldict.items():
1279 """Given a list of PMI objects, returns all density hierarchies within
1280 these objects. The output of this function can be inputted into
1281 things such as EM restraints. This function is intended to gather
1282 density particles appended to molecules (and not other hierarchies
1283 which might have been appended to the root node directly).
1292 i, representation_type=IMP.atom.DENSITIES).get_selected_particles()
1297 max_translation=300., max_rotation=2.0 * math.pi,
1298 avoidcollision_rb=
True, avoidcollision_fb=
False,
1299 cutoff=10.0, niterations=100,
1301 excluded_rigid_bodies=[],
1302 hierarchies_excluded_from_collision=[],
1303 hierarchies_included_in_collision=[],
1305 return_debug=
False):
1306 """Shuffle particles. Used to restart the optimization.
1307 The configuration of the system is initialized by placing each
1308 rigid body and each bead randomly in a box. If `bounding_box` is
1309 specified, the particles are placed inside this box; otherwise, each
1310 particle is displaced by up to max_translation angstroms, and randomly
1311 rotated. Effort is made to place particles far enough from each other to
1312 prevent any steric clashes.
1313 @param objects Can be one of the following inputs:
1314 IMP Hierarchy, PMI System/State/Molecule/TempResidue, or
1316 @param max_translation Max translation (rbs and flexible beads)
1317 @param max_rotation Max rotation (rbs only)
1318 @param avoidcollision_rb check if the particle/rigid body was
1319 placed close to another particle; uses the optional
1320 arguments cutoff and niterations
1321 @param avoidcollision_fb Advanced. Generally you want this False because
1322 it's hard to shuffle beads.
1323 @param cutoff Distance less than this is a collision
1324 @param niterations How many times to try avoiding collision
1325 @param bounding_box Only shuffle particles within this box.
1326 Defined by ((x1,y1,z1),(x2,y2,z2)).
1327 @param excluded_rigid_bodies Don't shuffle these rigid body objects
1328 @param hierarchies_excluded_from_collision Don't count collision
1330 @param hierarchies_included_in_collision Hierarchies that are not
1331 shuffled, but should be included in collision calculation
1333 @param verbose Give more output
1334 @note Best to only call this function after you've set up degrees
1336 For debugging purposes, returns: <shuffled indexes>,
1337 <collision avoided indexes>
1342 pmi_resolution=
'all',
1345 if len(rigid_bodies) > 0:
1346 mdl = rigid_bodies[0].get_model()
1347 elif len(flexible_beads) > 0:
1348 mdl = flexible_beads[0].get_model()
1350 raise Exception(
"Could not find any particles in the hierarchy")
1351 if len(rigid_bodies) == 0:
1352 print(
"shuffle_configuration: rigid bodies were not initialized")
1356 gcpf.set_distance(cutoff)
1360 hierarchies_excluded_from_collision, pmi_resolution=
'all',
1364 hierarchies_included_in_collision, pmi_resolution=
'all', flatten=
True)
1366 collision_excluded_idxs = set(
1368 for h
in collision_excluded_hierarchies
1371 collision_included_idxs = set(
1373 for h
in collision_included_hierarchies
1380 all_idxs.append(p.get_particle_index())
1382 collision_excluded_idxs.add(p.get_particle_index())
1384 if bounding_box
is not None:
1385 ((x1, y1, z1), (x2, y2, z2)) = bounding_box
1390 all_idxs = set(all_idxs) | collision_included_idxs
1391 all_idxs = all_idxs - collision_excluded_idxs
1393 print(
'shuffling', len(rigid_bodies),
'rigid bodies')
1394 for rb
in rigid_bodies:
1395 if rb
not in excluded_rigid_bodies:
1397 if avoidcollision_rb:
1398 rb_idxs = set(rb.get_member_particle_indexes()) - \
1399 collision_excluded_idxs
1400 other_idxs = all_idxs - rb_idxs
1404 while niter < niterations:
1405 rbxyz = (rb.get_x(), rb.get_y(), rb.get_z())
1422 rbxyz, max_translation, max_rotation)
1424 debug.append([rb, other_idxs
if avoidcollision_rb
else set()])
1428 if avoidcollision_rb
and other_idxs:
1430 npairs = len(gcpf.get_close_pairs(mdl,
1438 print(
"shuffle_configuration: rigid body placed "
1439 "close to other %d particles, trying "
1440 "again..." % npairs)
1441 print(
"shuffle_configuration: rigid body name: "
1443 if niter == niterations:
1445 "tried the maximum number of iterations to "
1446 "avoid collisions, increase the distance "
1451 print(
'shuffling', len(flexible_beads),
'flexible beads')
1452 for fb
in flexible_beads:
1454 if avoidcollision_fb:
1456 other_idxs = all_idxs - fb_idxs
1462 while niter < niterations:
1469 fbxyz, max_translation, max_rotation)
1474 xyz = memb.get_internal_coordinates()
1479 rf = memb.get_rigid_body().get_reference_frame()
1480 glob_to_int = rf.get_transformation_from()
1481 memb.set_internal_coordinates(
1482 glob_to_int.get_transformed(translation))
1484 xyz_transformed = transformation.get_transformed(xyz)
1485 memb.set_internal_coordinates(xyz_transformed)
1486 debug.append([xyz, other_idxs
if avoidcollision_fb
else set()])
1494 -d.get_coordinates())
1498 debug.append([d, other_idxs
if avoidcollision_fb
else set()])
1505 if avoidcollision_fb:
1507 npairs = len(gcpf.get_close_pairs(mdl,
1515 print(
"shuffle_configuration: floppy body placed close "
1516 "to other %d particles, trying again..." % npairs)
1517 if niter == niterations:
1519 "tried the maximum number of iterations to avoid "
1520 "collisions, increase the distance cutoff")
1527 class ColorHierarchy(object):
1529 def __init__(self, hier):
1530 import matplotlib
as mpl
1532 import matplotlib.pyplot
as plt
1536 hier.ColorHierarchy = self
1541 self.method = self.nochange
1549 def get_color(self, fl):
1552 def get_log_scale(self, fl):
1555 return math.log(fl+eps)
1557 def color_by_resid(self):
1558 self.method = self.color_by_resid
1559 self.scheme = self.mpl.cm.rainbow
1560 for mol
in self.mols:
1567 c = self.get_color(float(ri)/self.last)
1571 avr = sum(ris)/len(ris)
1572 c = self.get_color(float(avr)/self.last)
1575 def color_by_uncertainty(self):
1576 self.method = self.color_by_uncertainty
1577 self.scheme = self.mpl.cm.jet
1584 self.first = self.get_log_scale(1.0)
1585 self.last = self.get_log_scale(100.0)
1587 value = self.get_log_scale(unc_dict[p])
1588 if value >= self.last:
1590 if value <= self.first:
1592 c = self.get_color((value-self.first) / (self.last-self.first))
1595 def get_color_bar(self, filename):
1596 import matplotlib
as mpl
1598 import matplotlib.pyplot
as plt
1600 fig = plt.figure(figsize=(8, 3))
1601 ax1 = fig.add_axes([0.05, 0.80, 0.9, 0.15])
1604 norm = mpl.colors.Normalize(vmin=0.0, vmax=1.0)
1606 if self.method == self.color_by_uncertainty:
1607 angticks = [1.0, 2.5, 5.0, 10.0, 25.0, 50.0, 100.0]
1611 vvalue = (self.get_log_scale(at)-self.first) \
1612 / (self.last-self.first)
1613 if vvalue <= 1.0
and vvalue >= 0.0:
1614 vvalues.append(vvalue)
1615 marks.append(str(at))
1616 cb1 = mpl.colorbar.ColorbarBase(
1617 ax1, cmap=cmap, norm=norm, ticks=vvalues,
1618 orientation=
'horizontal')
1619 print(self.first, self.last, marks, vvalues)
1620 cb1.ax.set_xticklabels(marks)
1621 cb1.set_label(
'Angstorm')
1622 plt.savefig(filename, dpi=150, transparent=
True)
1627 """Given a Chimera color name or hex color value, return RGB"""
1628 d = {
'aquamarine': (0.4980392156862745, 1.0, 0.8313725490196079),
1629 'black': (0.0, 0.0, 0.0),
1630 'blue': (0.0, 0.0, 1.0),
1631 'brown': (0.6470588235, 0.16470588235294117, 0.16470588235294117),
1632 'chartreuse': (0.4980392156862745, 1.0, 0.0),
1633 'coral': (1.0, 0.4980392156862745, 0.3137254901960784),
1634 'cornflower blue': (0.39215686, 0.58431372549, 0.9294117647058824),
1635 'cyan': (0.0, 1.0, 1.0),
1636 'dark cyan': (0.0, 0.5450980392156862, 0.5450980392156862),
1637 'dark gray': (0.6627450980, 0.6627450980392157, 0.6627450980392157),
1638 'dark green': (0.0, 0.39215686274509803, 0.0),
1639 'dark khaki': (0.74117647, 0.7176470588235294, 0.4196078431372549),
1640 'dark magenta': (0.5450980392156862, 0.0, 0.5450980392156862),
1641 'dark olive green': (0.333333333, 0.419607843, 0.1843137254901961),
1642 'dark red': (0.5450980392156862, 0.0, 0.0),
1643 'dark slate blue': (0.28235294, 0.239215686, 0.5450980392156862),
1644 'dark slate gray': (0.1843137, 0.30980392, 0.30980392156862746),
1645 'deep pink': (1.0, 0.0784313725490196, 0.5764705882352941),
1646 'deep sky blue': (0.0, 0.7490196078431373, 1.0),
1647 'dim gray': (0.41176470, 0.4117647058823529, 0.4117647058823529),
1648 'dodger blue': (0.11764705882352941, 0.5647058823529412, 1.0),
1649 'firebrick': (0.6980392, 0.13333333333333333, 0.13333333333333333),
1650 'forest green': (0.13333333, 0.5450980392156862, 0.13333333333333333),
1651 'gold': (1.0, 0.8431372549019608, 0.0),
1652 'goldenrod': (0.85490196, 0.6470588235294118, 0.12549019607843137),
1653 'gray': (0.7450980392156863, 0.7450980392156863, 0.7450980392156863),
1654 'green': (0.0, 1.0, 0.0),
1655 'hot pink': (1.0, 0.4117647058823529, 0.7058823529411765),
1656 'khaki': (0.9411764705882353, 0.9019607843137255, 0.5490196078431373),
1657 'light blue': (0.67843137, 0.8470588235294118, 0.9019607843137255),
1658 'light gray': (0.82745098, 0.8274509803921568, 0.8274509803921568),
1659 'light green': (0.56470588, 0.9333333333333333, 0.5647058823529412),
1660 'light sea green': (0.125490, 0.6980392156862745, 0.6666666666666666),
1661 'lime green': (0.1960784, 0.803921568627451, 0.19607843137254902),
1662 'magenta': (1.0, 0.0, 1.0),
1663 'medium blue': (0.1960784, 0.19607843137254902, 0.803921568627451),
1664 'medium purple': (0.57647, 0.4392156862745098, 0.8588235294117647),
1665 'navy blue': (0.0, 0.0, 0.5019607843137255),
1666 'olive drab': (0.4196078, 0.5568627450980392, 0.13725490196078433),
1667 'orange red': (1.0, 0.27058823529411763, 0.0),
1668 'orange': (1.0, 0.4980392156862745, 0.0),
1669 'orchid': (0.85490196, 0.4392156862745098, 0.8392156862745098),
1670 'pink': (1.0, 0.7529411764705882, 0.796078431372549),
1671 'plum': (0.8666666666666667, 0.6274509803921569, 0.8666666666666667),
1672 'purple': (0.62745098, 0.12549019607843137, 0.9411764705882353),
1673 'red': (1.0, 0.0, 0.0),
1674 'rosy brown': (0.7372549, 0.5607843137254902, 0.5607843137254902),
1675 'salmon': (0.980392, 0.5019607843137255, 0.4470588235294118),
1676 'sandy brown': (0.956862745, 0.6431372549019608, 0.3764705882352941),
1677 'sea green': (0.18039, 0.5450980392156862, 0.3411764705882353),
1678 'sienna': (0.6274509, 0.3215686274509804, 0.17647058823529413),
1679 'sky blue': (0.52941176, 0.807843137254902, 0.9215686274509803),
1680 'slate gray': (0.439215686, 0.50196078, 0.5647058823529412),
1681 'spring green': (0.0, 1.0, 0.4980392156862745),
1682 'steel blue': (0.2745098, 0.50980392, 0.70588235),
1683 'tan': (0.8235294117647058, 0.7058823529411765, 0.5490196078431373),
1684 'turquoise': (0.25098039, 0.87843137, 0.81568627),
1685 'violet red': (0.81568627, 0.125490196, 0.56470588235),
1686 'white': (1.0, 1.0, 1.0),
1687 'yellow': (1.0, 1.0, 0.0)}
1688 if colorname.startswith(
'#'):
1689 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 ...