3 """@namespace IMP.pmi.tools
4 Miscellaneous utilities.
7 from __future__
import print_function
15 from math
import log,pi,sqrt,exp
20 from time
import process_time
22 from time
import clock
as process_time
25 from collections
import defaultdict
27 from collections
import OrderedDict
29 from IMP.pmi._compat_collections
import OrderedDict
32 def _get_system_for_hier(hier):
33 """Given a hierarchy, return the System that created it, or None"""
36 if hier
and not hasattr(hier,
'get_parent'):
43 if hasattr(hier,
'_pmi2_system'):
44 return hier._pmi2_system()
47 for ws
in IMP.pmi.topology.System._all_systems:
49 if s
and s.hier == hier:
52 hier = hier.get_parent()
54 def _all_protocol_outputs(representations, hier):
55 """Iterate over all (ProtocolOutput, State) pairs for the given
56 representations (PMI1) or hier (PMI2)"""
58 system = _get_system_for_hier(hier)
60 for state
in system.states:
61 for p
in state._protocol_output:
64 for p
in representations[0]._protocol_output:
67 def _add_pmi_provenance(p):
68 """Tag the given particle as being created by the current version of PMI."""
72 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)
82 def _add_restraint_sets(model, mk, mk_rmf):
85 model.add_data(mk, rs)
86 model.add_data(mk_rmf, rs_rmf)
90 """Add a PMI restraint to the model.
91 Since Model.add_restraint() no longer exists (in modern IMP restraints
92 should be added to a ScoringFunction instead) store them instead in
93 a RestraintSet, and keep a reference to it in the Model.
95 If `add_to_rmf` is True, also add the restraint to a separate list
96 of restraints that will be written out to RMF files (by default, most
97 PMI restraints are not)."""
98 mk, mk_rmf = _get_restraint_set_keys()
99 if model.get_has_data(mk):
100 rs = IMP.RestraintSet.get_from(model.get_data(mk))
101 rs_rmf = IMP.RestraintSet.get_from(model.get_data(mk_rmf))
103 rs, rs_rmf = _add_restraint_sets(model, mk, mk_rmf)
104 rs.add_restraint(restraint)
106 rs_rmf.add_restraint(restraint)
109 """Get a RestraintSet containing all PMI restraints added to the model.
110 If `rmf` is True, return only the subset of these restraints that
111 should be written out to RMF files."""
112 mk, mk_rmf = _get_restraint_set_keys()
113 if not model.get_has_data(mk):
114 warnings.warn(
"no restraints added to model yet",
116 _add_restraint_sets(model, mk, mk_rmf)
118 return IMP.RestraintSet.get_from(model.get_data(mk_rmf))
120 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):
157 nuisance.set_lower(minvalue)
159 nuisance.set_upper(maxvalue)
162 nuisance.set_is_optimized(nuisance.get_nuisance_key(), isoptimized)
163 self.nuisance = nuisance
165 def get_particle(self):
169 class SetupWeight(object):
171 def __init__(self, m, isoptimized=True, nweights_or_weights=None):
173 if isinstance(nweights_or_weights, int):
175 pw, nweights_or_weights
179 nweights_or_weights = list(nweights_or_weights)
181 pw, nweights_or_weights
185 self.weight.set_weights_are_optimized(isoptimized)
187 def get_particle(self):
191 class SetupSurface(object):
193 def __init__(self, m, center, normal, isoptimized=True):
196 self.surface.set_coordinates_are_optimized(isoptimized)
197 self.surface.set_normal_is_optimized(isoptimized)
199 def get_particle(self):
203 class ParticleToSampleList(object):
207 self.dictionary_particle_type = {}
208 self.dictionary_particle_transformation = {}
209 self.dictionary_particle_name = {}
216 particle_transformation,
218 if not particle_type
in [
"Rigid_Bodies",
"Floppy_Bodies",
"Nuisances",
"X_coord",
"Weights",
"Surfaces"]:
219 raise TypeError(
"not the right particle type")
221 self.dictionary_particle_type[particle] = particle_type
222 if particle_type ==
"Rigid_Bodies":
223 if (isinstance(particle_transformation, tuple)
224 and len(particle_transformation) == 2
225 and all(isinstance(x, float)
226 for x
in particle_transformation)):
227 self.dictionary_particle_transformation[
228 particle] = particle_transformation
229 self.dictionary_particle_name[particle] = name
231 raise TypeError(
"ParticleToSampleList: not the right transformation format for Rigid_Bodies, should be a tuple of floats")
232 elif particle_type ==
"Surfaces":
233 if (isinstance(particle_transformation, tuple)
234 and len(particle_transformation) == 3
235 and all(isinstance(x, float)
236 for x
in particle_transformation)):
237 self.dictionary_particle_transformation[
238 particle] = particle_transformation
239 self.dictionary_particle_name[particle] = name
241 raise TypeError(
"ParticleToSampleList: not the right transformation format for Surfaces, should be a tuple of floats")
243 if isinstance(particle_transformation, float):
244 self.dictionary_particle_transformation[
245 particle] = particle_transformation
246 self.dictionary_particle_name[particle] = name
248 raise TypeError(
"ParticleToSampleList: not the right transformation format, should be a float")
250 def get_particles_to_sample(self):
252 for particle
in self.dictionary_particle_type:
253 key = self.dictionary_particle_type[
254 particle] +
"ParticleToSampleList_" + self.dictionary_particle_name[particle] +
"_" + self.label
257 self.dictionary_particle_transformation[particle])
262 def get_cross_link_data(directory, filename, dist, omega, sigma,
263 don=
None, doff=
None, prior=0, type_of_profile=
"gofr"):
265 (distmin, distmax, ndist) = dist
266 (omegamin, omegamax, nomega) = omega
267 (sigmamin, sigmamax, nsigma) = sigma
270 with open(filen)
as xlpot:
271 dictionary = ast.literal_eval(xlpot.readline())
273 xpot = dictionary[directory][filename][
"distance"]
274 pot = dictionary[directory][filename][type_of_profile]
276 dist_grid =
get_grid(distmin, distmax, ndist,
False)
277 omega_grid = get_log_grid(omegamin, omegamax, nomega)
278 sigma_grid = get_log_grid(sigmamin, sigmamax, nsigma)
280 if not don
is None and not doff
is None:
299 def get_grid(gmin, gmax, ngrid, boundaries):
301 dx = (gmax - gmin) / float(ngrid)
302 for i
in range(0, ngrid + 1):
303 if(
not boundaries
and i == 0):
305 if(
not boundaries
and i == ngrid):
307 grid.append(gmin + float(i) * dx)
313 def get_log_grid(gmin, gmax, ngrid):
315 for i
in range(0, ngrid + 1):
316 grid.append(gmin * exp(float(i) / ngrid * log(gmax / gmin)))
324 example '"{ID_Score}" > 28 AND "{Sample}" ==
325 "%10_1%" OR ":Sample}" == "%10_2%" OR ":Sample}"
326 == "%10_3%" OR ":Sample}" == "%8_1%" OR ":Sample}" == "%8_2%"'
329 import pyparsing
as pp
331 operator = pp.Regex(
">=|<=|!=|>|<|==|in").setName(
"operator")
332 value = pp.QuotedString(
334 r"[+-]?\d+(:?\.\d*)?(:?[eE][+-]?\d+)?")
335 identifier = pp.Word(pp.alphas, pp.alphanums +
"_")
336 comparison_term = identifier | value
337 condition = pp.Group(comparison_term + operator + comparison_term)
339 expr = pp.operatorPrecedence(condition, [
340 (
"OR", 2, pp.opAssoc.LEFT, ),
341 (
"AND", 2, pp.opAssoc.LEFT, ),
344 parsedstring = str(expr.parseString(inputstring)) \
350 .replace(
"{",
"float(entry['") \
351 .replace(
"}",
"'])") \
352 .replace(
":",
"str(entry['") \
353 .replace(
"}",
"'])") \
354 .replace(
"AND",
"and") \
359 def open_file_or_inline_text(filename):
361 fl = open(filename,
"r")
363 fl = filename.split(
"\n")
366 def get_ids_from_fasta_file(fastafile):
368 with open(fastafile)
as ff:
377 this function works with plain hierarchies, as read from the pdb,
378 no multi-scale hierarchies
385 atom_type=IMP.atom.AT_CA)
393 print(
"get_closest_residue_position: exiting while loop without result")
395 p = sel.get_selected_particles()
400 print(
"get_closest_residue_position: got NO residues for hierarchy %s and residue %i" % (hier, resindex))
401 raise Exception(
"get_closest_residue_position: got NO residues for hierarchy %s and residue %i" % (
404 raise ValueError(
"got multiple residues for hierarchy %s and residue %i; the list of particles is %s" % (hier, resindex, str([pp.get_name()
for pp
in p])))
409 Return the residue index gaps and contiguous segments in the hierarchy.
411 @param hierarchy hierarchy to examine
412 @param start first residue index
413 @param end last residue index
415 @return A list of lists of the form
416 [[1,100,"cont"],[101,120,"gap"],[121,200,"cont"]]
419 for n, rindex
in enumerate(range(start, end + 1)):
421 atom_type=IMP.atom.AT_CA)
423 if len(sel.get_selected_particles()) == 0:
427 rindexcont = start - 1
428 if rindexgap == rindex - 1:
434 gaps.append([rindex, rindex,
"gap"])
440 rindexgap = start - 1
442 if rindexcont == rindex - 1:
449 gaps.append([rindex, rindex,
"cont"])
460 def set_map_element(self, xvalue, yvalue):
461 self.map[xvalue] = yvalue
463 def get_map_element(self, invalue):
464 if isinstance(invalue, float):
468 dist = (invalue - x) * (invalue - x)
477 return self.map[minx]
478 elif isinstance(invalue, str):
479 return self.map[invalue]
481 raise TypeError(
"wrong type for map")
485 """New tuple format: molname OR (start,stop,molname,copynum,statenum)
486 Copy and state are optional. Can also use 'None' for them which will get all.
487 You can also pass -1 for stop which will go to the end.
488 Returns the particles
491 kwds[
'resolution'] = resolution
492 if isinstance(tuple_selection, str):
493 kwds[
'molecule'] = tuple_selection
494 elif isinstance(tuple_selection, tuple):
495 rbegin = tuple_selection[0]
496 rend = tuple_selection[1]
497 kwds[
'molecule'] = tuple_selection[2]
499 copynum = tuple_selection[3]
500 if copynum
is not None:
501 kwds[
'copy_index'] = copynum
505 statenum = tuple_selection[4]
506 if statenum
is not None:
507 kwds[
'state_index'] = statenum
514 residue_indexes=range(1,rbegin),
516 return s.get_selected_particles()
518 kwds[
'residue_indexes'] = range(rbegin,rend+1)
520 return s.get_selected_particles()
524 def get_db_from_csv(csvfilename):
527 with open(csvfilename)
as fh:
528 csvr = csv.DictReader(fh)
535 '''Return the component name provided a particle and a list of names'''
537 protname = root.get_name()
539 while not protname
in list_of_names:
540 root0 = root.get_parent()
543 protname = root0.get_name()
548 if "Beads" in protname:
551 return (protname, is_a_bead)
556 Retrieve the residue indexes for the given particle.
558 The particle must be an instance of Fragment,Residue, Atom or Molecule
559 or else returns an empty list
570 resind_tmp=IMP.pmi.tools.OrderedSet()
576 resind=list(resind_tmp)
582 def sort_by_residues(particles):
585 sorted_particles_residues = sorted(
587 key=
lambda tup: tup[1])
588 particles = [p[0]
for p
in sorted_particles_residues]
597 """Synchronize data over a parallel run"""
598 from mpi4py
import MPI
599 comm = MPI.COMM_WORLD
600 rank = comm.Get_rank()
601 number_of_processes = comm.size
604 comm.send(data, dest=0, tag=11)
607 for i
in range(1, number_of_processes):
608 data_tmp = comm.recv(source=i, tag=11)
609 if isinstance(data, list):
611 elif isinstance(data, dict):
612 data.update(data_tmp)
614 raise TypeError(
"data not supported, use list or dictionaries")
616 for i
in range(1, number_of_processes):
617 comm.send(data, dest=i, tag=11)
620 data = comm.recv(source=0, tag=11)
629 Yield all sublists of length >= lmin and <= lmax
637 for j
in range(i + 1, n):
638 if len(l[i:j]) <= lmax
and len(l[i:j]) >= lmin:
643 return [item
for sublist
in l
for item
in sublist]
647 """ Yield successive length-sized chunks from a list.
649 for i
in range(0, len(list), length):
650 yield list[i:i + length]
653 def chunk_list_into_segments(seq, num):
655 avg = len(seq) / float(num)
659 while last < len(seq):
660 out.append(seq[int(last):int(last + avg)])
668 ''' This class stores integers
669 in ordered compact lists eg:
671 the methods help splitting and merging the internal lists
673 s=Segments([1,2,3]) is [[1,2,3]]
674 s.add(4) is [[1,2,3,4]] (add right)
675 s.add(3) is [[1,2,3,4]] (item already existing)
676 s.add(7) is [[1,2,3,4],[7]] (new list)
677 s.add([8,9]) is [[1,2,3,4],[7,8,9]] (add item right)
678 s.add([5,6]) is [[1,2,3,4,5,6,7,8,9]] (merge)
679 s.remove(3) is [[1,2],[4,5,6,7,8,9]] (split)
684 '''index can be a integer or a list of integers '''
685 if isinstance(index, int):
687 elif isinstance(index, list):
688 self.segs=[[index[0]]]
692 raise TypeError(
"index must be an int or list of ints")
695 '''index can be a integer or a list of integers '''
696 if isinstance(index, int):
699 for n,s
in enumerate(self.segs):
707 if mergeright
is None and mergeleft
is None:
708 self.segs.append([index])
709 if not mergeright
is None and mergeleft
is None:
710 self.segs[mergeright].append(index)
711 if not mergeleft
is None and mergeright
is None:
712 self.segs[mergeleft]=[index]+self.segs[mergeleft]
713 if not mergeleft
is None and not mergeright
is None:
714 self.segs[mergeright]=self.segs[mergeright]+[index]+self.segs[mergeleft]
715 del self.segs[mergeleft]
717 for n
in range(len(self.segs)):
720 self.segs.sort(key=
lambda tup: tup[0])
722 elif isinstance(index, list):
726 raise TypeError(
"index must be an int or list of ints")
729 '''index can be a integer'''
730 for n,s
in enumerate(self.segs):
737 i=self.segs[n].index(index)
739 self.segs.append(s[i+1:])
740 for n
in range(len(self.segs)):
742 if len(self.segs[n])==0:
744 self.segs.sort(key=
lambda tup: tup[0])
747 ''' Returns a flatten list '''
748 return [item
for sublist
in self.segs
for item
in sublist]
752 for seg
in self.segs:
753 ret_tmp+=str(seg[0])+
"-"+str(seg[-1])+
","
761 def normal_density_function(expected_value, sigma, x):
763 1 / math.sqrt(2 * math.pi) / sigma *
764 math.exp(-(x - expected_value) ** 2 / 2 / sigma / sigma)
768 def log_normal_density_function(expected_value, sigma, x):
770 1 / math.sqrt(2 * math.pi) / sigma / x *
771 math.exp(-(math.log(x / expected_value) ** 2 / 2 / sigma / sigma))
775 def print_multicolumn(list_of_strings, ncolumns=2, truncate=40):
781 for i
in range(len(l) % cols):
784 split = [l[i:i + len(l) / cols]
for i
in range(0, len(l), len(l) / cols)]
785 for row
in zip(*split):
786 print(
"".join(str.ljust(i, truncate)
for i
in row))
789 '''Change color code to hexadecimal to rgb'''
791 self._NUMERALS =
'0123456789abcdefABCDEF'
792 self._HEXDEC = dict((v, int(v, 16))
for v
in (x+y
for x
in self._NUMERALS
for y
in self._NUMERALS))
793 self.LOWERCASE, self.UPPERCASE =
'x',
'X'
795 def rgb(self,triplet):
796 return float(self._HEXDEC[triplet[0:2]]), float(self._HEXDEC[triplet[2:4]]), float(self._HEXDEC[triplet[4:6]])
798 def triplet(self,rgb, lettercase=None):
799 if lettercase
is None: lettercase=self.LOWERCASE
800 return format(rgb[0]<<16 | rgb[1]<<8 | rgb[2],
'06'+lettercase)
804 class OrderedSet(collections.MutableSet):
806 def __init__(self, iterable=None):
808 end += [
None, end, end]
810 if iterable
is not None:
816 def __contains__(self, key):
817 return key
in self.map
820 if key
not in self.map:
823 curr[2] = end[1] = self.map[key] = [key, curr, end]
825 def discard(self, key):
827 key, prev, next = self.map.pop(key)
834 while curr
is not end:
838 def __reversed__(self):
841 while curr
is not end:
845 def pop(self, last=True):
847 raise KeyError(
'set is empty')
857 return '%s()' % (self.__class__.__name__,)
858 return '%s(%r)' % (self.__class__.__name__, list(self))
860 def __eq__(self, other):
861 if isinstance(other, OrderedSet):
862 return len(self) == len(other)
and list(self) == list(other)
863 return set(self) == set(other)
867 """Store objects in order they were added, but with default type.
868 Source: http://stackoverflow.com/a/4127426/2608793
870 def __init__(self, *args, **kwargs):
872 self.default_factory =
None
874 if not (args[0]
is None or callable(args[0])):
875 raise TypeError(
'first argument must be callable or None')
876 self.default_factory = args[0]
878 super(OrderedDefaultDict, self).__init__(*args, **kwargs)
880 def __missing__ (self, key):
881 if self.default_factory
is None:
883 self[key] = default = self.default_factory()
886 def __reduce__(self):
887 args = (self.default_factory,)
if self.default_factory
else ()
888 if sys.version_info[0] >= 3:
889 return self.__class__, args,
None,
None, self.items()
891 return self.__class__, args,
None,
None, self.iteritems()
897 """Extract frame from RMF file and fill coordinates. Must be identical topology.
898 @param hier The (System) hierarchy to fill (e.g. after you've built it)
899 @param rmf_fn The file to extract from
900 @param frame_num The frame number to extract
902 rh = RMF.open_rmf_file_read_only(rmf_fn)
910 selection_tuple=
None,
911 warn_about_slices=
True):
912 """Adapt things for PMI (degrees of freedom, restraints, ...)
913 Returns list of list of hierarchies, separated into Molecules if possible.
914 The input can be a list, or a list of lists (iterable of ^1 or iterable of ^2)
915 (iterable of ^2) Hierarchy -> returns input as list of list of hierarchies,
916 only one entry, not grouped by molecules.
917 (iterable of ^2) PMI::System/State/Molecule/TempResidue ->
918 returns residue hierarchies, grouped in molecules, at requested resolution
919 @param stuff Can be one of the following inputs:
920 IMP Hierarchy, PMI System/State/Molecule/TempResidue, or a list/set (of list/set) of them.
921 Must be uniform input, however. No mixing object types.
922 @param pmi_resolution For selecting, only does it if you pass PMI objects. Set it to "all"
923 if you want all resolutions!
924 @param flatten Set to True if you just want all hierarchies in one list.
925 @param warn_about_slices Print a warning if you are requesting only part
926 of a bead. Sometimes you just don't care!
927 @note since this relies on IMP::atom::Selection, this will not return
928 any objects if they weren't built! But there should be no problem
929 if you request unbuilt residues - they should be ignored.
935 if hasattr(stuff,
'__iter__'):
941 if all(hasattr(el,
'__iter__')
for el
in thelist):
942 thelist = [i
for sublist
in thelist
for i
in sublist]
943 elif any(hasattr(el,
'__iter__')
for el
in thelist):
944 raise Exception(
'input_adaptor: input_object must be a list or a list of lists')
953 except (NotImplementedError, TypeError):
964 if is_system
or is_state
or is_molecule
or is_temp_residue:
970 for state
in system.get_states():
971 mdict = state.get_molecules()
972 for molname
in mdict:
973 for copy
in mdict[molname]:
974 indexes_per_mol[copy] += [r.get_index()
for r
in copy.get_residues()]
977 mdict = state.get_molecules()
978 for molname
in mdict:
979 for copy
in mdict[molname]:
980 indexes_per_mol[copy] += [r.get_index()
for r
in copy.get_residues()]
982 for molecule
in stuff:
983 indexes_per_mol[molecule] += [r.get_index()
for r
in molecule.get_residues()]
985 for tempres
in stuff:
986 indexes_per_mol[tempres.get_molecule()].append(tempres.get_index())
987 for mol
in indexes_per_mol:
988 if pmi_resolution==
'all':
992 residue_indexes=indexes_per_mol[mol])
995 resolution=pmi_resolution,
996 residue_indexes=indexes_per_mol[mol])
997 ps = sel.get_selected_particles()
1000 if warn_about_slices:
1001 rset = set(indexes_per_mol[mol])
1005 if not fset <= rset:
1011 resbreak = maxf
if minf==minset
else minset-1
1013 'You are trying to select only part of the '
1014 'bead %s:%i-%i. The residues you requested '
1015 'are %i-%i. You can fix this by: '
1016 '1) requesting the whole bead/none of it; or'
1017 '2) break the bead up by passing '
1018 'bead_extra_breaks=[\'%i\'] in '
1019 'molecule.add_representation()'
1020 %(mol.get_name(), minset, maxset, minf, maxf,
1026 if pmi_resolution==
'all':
1034 hier_list = [hier_list]
1036 raise Exception(
'input_adaptor: you passed something of wrong type or a list with mixed types')
1038 if flatten
and pmi_input:
1039 return [h
for sublist
in hier_list
for h
in sublist]
1045 """Returns sequence-sorted segments array, each containing the first particle
1046 the last particle and the first residue index."""
1048 from operator
import itemgetter
1051 raise ValueError(
"only pass stuff from one Molecule, please")
1067 segs.append((start, end, startres))
1068 return sorted(segs, key=itemgetter(2))
1071 """Decorate the sequence-consecutive particles from a PMI2 molecule with a bond,
1072 so that they appear connected in the rmf file"""
1074 for x
in range(len(SortedSegments) - 1):
1076 last = SortedSegments[x][1]
1077 first = SortedSegments[x + 1][0]
1079 p1 = last.get_particle()
1080 p2 = first.get_particle()
1093 """This class converts three to one letter codes, and return X for any unknown codes"""
1094 def __init__(self,is_nucleic=False):
1097 threetoone = {
'ALA':
'A',
'ARG':
'R', 'ASN': 'N', 'ASP': 'D',
1098 'CYS':
'C',
'GLU':
'E',
'GLN':
'Q',
'GLY':
'G',
1099 'HIS':
'H',
'ILE':
'I',
'LEU':
'L',
'LYS':
'K',
1100 'MET':
'M',
'PHE':
'F',
'PRO':
'P',
'SER':
'S',
1101 'THR':
'T',
'TRP':
'W',
'TYR':
'Y',
'VAL':
'V',
'UNK':
'X'}
1103 threetoone = {
'ADE':
'A',
'URA':
'U', 'CYT': 'C', 'GUA': 'G',
1104 'THY':
'T',
'UNK':
'X'}
1106 defaultdict.__init__(self,
lambda:
"X", threetoone)
1111 def get_residue_type_from_one_letter_code(code,is_nucleic=False):
1114 for k
in threetoone:
1115 one_to_three[threetoone[k]] = k
1120 """ Just get the leaves from a list of hierarchies """
1121 lvs = list(itertools.chain.from_iterable(
IMP.atom.get_leaves(item)
for item
in list_of_hs))
1128 """Perform selection using the usual keywords but return ALL resolutions (BEADS and GAUSSIANS).
1129 Returns in flat list!
1134 if hier
is not None:
1137 warnings.warn(
"You passed nothing to select_at_all_resolutions()",
1145 raise Exception(
'select_at_all_resolutions: you have to pass an IMP Hierarchy')
1147 raise Exception(
'select_at_all_resolutions: you have to pass an IMP Hierarchy')
1148 if 'resolution' in kwargs
or 'representation_type' in kwargs:
1149 raise Exception(
"don't pass resolution or representation_type to this function")
1151 representation_type=IMP.atom.BALLS,**kwargs)
1153 representation_type=IMP.atom.DENSITIES,**kwargs)
1154 ret |= OrderedSet(selB.get_selected_particles())
1155 ret |= OrderedSet(selD.get_selected_particles())
1164 """Utility to retrieve particles from a hierarchy within a
1165 zone around a set of ps.
1166 @param hier The hierarchy in which to look for neighbors
1167 @param target_ps The particles for zoning
1168 @param sel_zone The maximum distance
1169 @param entire_residues If True, will grab entire residues
1170 @param exclude_backbone If True, will only return sidechain particles
1174 backbone_types=[
'C',
'N',
'CB',
'O']
1175 if exclude_backbone:
1177 for n
in backbone_types])
1178 test_ps = test_sel.get_selected_particles()
1179 nn = IMP.algebra.NearestNeighbor3D([
IMP.core.XYZ(p).get_coordinates()
1182 for target
in target_ps:
1183 zone|=set(nn.get_in_ball(
IMP.core.XYZ(target).get_coordinates(),sel_zone))
1184 zone_ps = [test_ps[z]
for z
in zone]
1189 zone_ps = [h.get_particle()
for h
in final_ps]
1194 """Returns unique objects in original order"""
1198 if not hasattr(hiers,
'__iter__'):
1205 rbs_ordered.append(rb)
1210 rbs_ordered.append(rb)
1214 return rbs_ordered,beads
1217 "This function returns the parent molecule hierarchies of given objects"
1218 stuff=
input_adaptor(input_objects, pmi_resolution=
'all',flatten=
True)
1223 while not (is_molecule
or is_root):
1224 root=IMP.atom.get_root(h)
1231 return list(molecules)
1233 def get_molecules_dictionary(input_objects):
1234 moldict=defaultdict(list)
1237 moldict[name].append(mol)
1243 def get_molecules_dictionary_by_copy(input_objects):
1244 moldict=defaultdict(dict)
1248 moldict[name][c]=mol
1251 def get_selections_dictionary(input_objects):
1252 moldict=IMP.pmi.tools.get_molecules_dictionary(input_objects)
1253 seldict=defaultdict(list)
1254 for name, mols
in moldict.items():
1260 """Given a list of PMI objects, returns all density hierarchies within
1261 these objects. The output of this function can be inputted into
1262 things such as EM restraints. This function is intended to gather density particles
1263 appended to molecules (and not other hierarchies which might have been appended to the root node directly).
1270 densities+=
IMP.atom.Selection(i,representation_type=IMP.atom.DENSITIES).get_selected_particles()
1274 max_translation=300., max_rotation=2.0 * pi,
1275 avoidcollision_rb=
True, avoidcollision_fb=
False,
1276 cutoff=10.0, niterations=100,
1278 excluded_rigid_bodies=[],
1279 hierarchies_excluded_from_collision=[],
1280 hierarchies_included_in_collision=[],
1282 return_debug=
False):
1283 """Shuffle particles. Used to restart the optimization.
1284 The configuration of the system is initialized by placing each
1285 rigid body and each bead randomly in a box. If `bounding_box` is
1286 specified, the particles are placed inside this box; otherwise, each
1287 particle is displaced by up to max_translation angstroms, and randomly
1288 rotated. Effort is made to place particles far enough from each other to
1289 prevent any steric clashes.
1290 @param objects Can be one of the following inputs:
1291 IMP Hierarchy, PMI System/State/Molecule/TempResidue, or a list/set of them
1292 @param max_translation Max translation (rbs and flexible beads)
1293 @param max_rotation Max rotation (rbs only)
1294 @param avoidcollision_rb check if the particle/rigid body was
1295 placed close to another particle; uses the optional
1296 arguments cutoff and niterations
1297 @param avoidcollision_fb Advanced. Generally you want this False because it's hard to shuffle beads.
1298 @param cutoff Distance less than this is a collision
1299 @param niterations How many times to try avoiding collision
1300 @param bounding_box Only shuffle particles within this box. Defined by ((x1,y1,z1),(x2,y2,z2)).
1301 @param excluded_rigid_bodies Don't shuffle these rigid body objects
1302 @param hierarchies_excluded_from_collision Don't count collision with these bodies
1303 @param hierarchies_included_in_collision Hierarchies that are not shuffled, but should be included in collision calculation (for fixed regions)
1304 @param verbose Give more output
1305 @note Best to only call this function after you've set up degrees of freedom
1306 For debugging purposes, returns: <shuffled indexes>, <collision avoided indexes>
1311 pmi_resolution=
'all',
1314 if len(rigid_bodies)>0:
1315 mdl = rigid_bodies[0].get_model()
1316 elif len(flexible_beads)>0:
1317 mdl = flexible_beads[0].get_model()
1319 raise Exception(
"Could not find any particles in the hierarchy")
1320 if len(rigid_bodies) == 0:
1321 print(
"shuffle_configuration: rigid bodies were not initialized")
1325 gcpf.set_distance(cutoff)
1329 pmi_resolution=
'all',
1333 pmi_resolution=
'all',
1336 collision_excluded_idxs = set([l.get_particle().
get_index()
for h
in collision_excluded_hierarchies \
1339 collision_included_idxs = set([l.get_particle().
get_index()
for h
in collision_included_hierarchies \
1346 all_idxs.append(p.get_particle_index())
1348 collision_excluded_idxs.add(p.get_particle_index())
1350 if bounding_box
is not None:
1351 ((x1, y1, z1), (x2, y2, z2)) = bounding_box
1356 all_idxs = set(all_idxs) | collision_included_idxs
1357 all_idxs = all_idxs - collision_excluded_idxs
1359 print(
'shuffling', len(rigid_bodies),
'rigid bodies')
1360 for rb
in rigid_bodies:
1361 if rb
not in excluded_rigid_bodies:
1363 if avoidcollision_rb:
1364 rb_idxs = set(rb.get_member_particle_indexes()) - \
1365 collision_excluded_idxs
1366 other_idxs = all_idxs - rb_idxs
1372 while niter < niterations:
1373 rbxyz = (rb.get_x(), rb.get_y(), rb.get_z())
1393 debug.append([rb, other_idxs
if avoidcollision_rb
else set()])
1397 if avoidcollision_rb:
1399 npairs = len(gcpf.get_close_pairs(mdl,
1408 print(
"shuffle_configuration: rigid body placed close to other %d particles, trying again..." % npairs)
1409 print(
"shuffle_configuration: rigid body name: " + rb.get_name())
1410 if niter == niterations:
1411 raise ValueError(
"tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1415 print(
'shuffling', len(flexible_beads),
'flexible beads')
1416 for fb
in flexible_beads:
1418 if avoidcollision_fb:
1420 other_idxs = all_idxs - fb_idxs
1426 while niter < niterations:
1439 xyz = memb.get_internal_coordinates()
1444 rf = memb.get_rigid_body().get_reference_frame()
1445 glob_to_int = rf.get_transformation_from()
1446 memb.set_internal_coordinates(
1447 glob_to_int.get_transformed(translation))
1449 xyz_transformed=transformation.get_transformed(xyz)
1450 memb.set_internal_coordinates(xyz_transformed)
1451 debug.append([xyz,other_idxs
if avoidcollision_fb
else set()])
1458 debug.append([d,other_idxs
if avoidcollision_fb
else set()])
1464 if avoidcollision_fb:
1466 npairs = len(gcpf.get_close_pairs(mdl,
1474 print(
"shuffle_configuration: floppy body placed close to other %d particles, trying again..." % npairs)
1475 if niter == niterations:
1476 raise ValueError(
"tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
1482 class ColorHierarchy(object):
1484 def __init__(self,hier):
1485 import matplotlib
as mpl
1487 import matplotlib.pyplot
as plt
1491 hier.ColorHierarchy=self
1495 self.method=self.nochange
1503 def get_color(self,fl):
1506 def get_log_scale(self,fl):
1509 return math.log(fl+eps)
1511 def color_by_resid(self):
1512 self.method=self.color_by_resid
1513 self.scheme=self.mpl.cm.rainbow
1514 for mol
in self.mols:
1520 c=self.get_color(float(ri)/self.last)
1524 avr=sum(ris)/len(ris)
1525 c=self.get_color(float(avr)/self.last)
1528 def color_by_uncertainty(self):
1529 self.method=self.color_by_uncertainty
1530 self.scheme=self.mpl.cm.jet
1537 self.first=self.get_log_scale(1.0)
1538 self.last=self.get_log_scale(100.0)
1540 value=self.get_log_scale(unc_dict[p])
1541 if value>=self.last: value=self.last
1542 if value<=self.first: value=self.first
1543 c=self.get_color((value-self.first)/(self.last-self.first))
1546 def get_color_bar(self,filename):
1547 import matplotlib
as mpl
1549 import matplotlib.pyplot
as plt
1552 fig = plt.figure(figsize=(8, 3))
1553 ax1 = fig.add_axes([0.05, 0.80, 0.9, 0.15])
1556 norm = mpl.colors.Normalize(vmin=0.0, vmax=1.0)
1558 if self.method == self.color_by_uncertainty:
1559 angticks=[1.0,2.5,5.0,10.0,25.0,50.0,100.0]
1563 vvalue=(self.get_log_scale(at)-self.first)/(self.last-self.first)
1564 if vvalue <= 1.0
and vvalue >= 0.0:
1565 vvalues.append(vvalue)
1566 marks.append(str(at))
1567 cb1 = mpl.colorbar.ColorbarBase(ax1, cmap=cmap,
1570 orientation=
'horizontal')
1571 print(self.first,self.last,marks,vvalues)
1572 cb1.ax.set_xticklabels(marks)
1573 cb1.set_label(
'Angstorm')
1574 plt.savefig(filename, dpi=150, transparent=
True)
1581 """Given a Chimera color name or hex color value, return RGB"""
1582 d = {
'aquamarine': (0.4980392156862745, 1.0, 0.8313725490196079),
1583 'black': (0.0, 0.0, 0.0),
1584 'blue': (0.0, 0.0, 1.0),
1585 'brown': (0.6470588235294118, 0.16470588235294117, 0.16470588235294117),
1586 'chartreuse': (0.4980392156862745, 1.0, 0.0),
1587 'coral': (1.0, 0.4980392156862745, 0.3137254901960784),
1588 'cornflower blue': (0.39215686274509803, 0.5843137254901961, 0.9294117647058824),
1589 'cyan': (0.0, 1.0, 1.0),
1590 'dark cyan': (0.0, 0.5450980392156862, 0.5450980392156862),
1591 'dark gray': (0.6627450980392157, 0.6627450980392157, 0.6627450980392157),
1592 'dark green': (0.0, 0.39215686274509803, 0.0),
1593 'dark khaki': (0.7411764705882353, 0.7176470588235294, 0.4196078431372549),
1594 'dark magenta': (0.5450980392156862, 0.0, 0.5450980392156862),
1595 'dark olive green': (0.3333333333333333, 0.4196078431372549, 0.1843137254901961),
1596 'dark red': (0.5450980392156862, 0.0, 0.0),
1597 'dark slate blue': (0.2823529411764706, 0.23921568627450981, 0.5450980392156862),
1598 'dark slate gray': (0.1843137254901961, 0.30980392156862746, 0.30980392156862746),
1599 'deep pink': (1.0, 0.0784313725490196, 0.5764705882352941),
1600 'deep sky blue': (0.0, 0.7490196078431373, 1.0),
1601 'dim gray': (0.4117647058823529, 0.4117647058823529, 0.4117647058823529),
1602 'dodger blue': (0.11764705882352941, 0.5647058823529412, 1.0),
1603 'firebrick': (0.6980392156862745, 0.13333333333333333, 0.13333333333333333),
1604 'forest green': (0.13333333333333333, 0.5450980392156862, 0.13333333333333333),
1605 'gold': (1.0, 0.8431372549019608, 0.0),
1606 'goldenrod': (0.8549019607843137, 0.6470588235294118, 0.12549019607843137),
1607 'gray': (0.7450980392156863, 0.7450980392156863, 0.7450980392156863),
1608 'green': (0.0, 1.0, 0.0),
1609 'hot pink': (1.0, 0.4117647058823529, 0.7058823529411765),
1610 'khaki': (0.9411764705882353, 0.9019607843137255, 0.5490196078431373),
1611 'light blue': (0.6784313725490196, 0.8470588235294118, 0.9019607843137255),
1612 'light gray': (0.8274509803921568, 0.8274509803921568, 0.8274509803921568),
1613 'light green': (0.5647058823529412, 0.9333333333333333, 0.5647058823529412),
1614 'light sea green': (0.12549019607843137, 0.6980392156862745, 0.6666666666666666),
1615 'lime green': (0.19607843137254902, 0.803921568627451, 0.19607843137254902),
1616 'magenta': (1.0, 0.0, 1.0),
1617 'medium blue': (0.19607843137254902, 0.19607843137254902, 0.803921568627451),
1618 'medium purple': (0.5764705882352941, 0.4392156862745098, 0.8588235294117647),
1619 'navy blue': (0.0, 0.0, 0.5019607843137255),
1620 'olive drab': (0.4196078431372549, 0.5568627450980392, 0.13725490196078433),
1621 'orange red': (1.0, 0.27058823529411763, 0.0),
1622 'orange': (1.0, 0.4980392156862745, 0.0),
1623 'orchid': (0.8549019607843137, 0.4392156862745098, 0.8392156862745098),
1624 'pink': (1.0, 0.7529411764705882, 0.796078431372549),
1625 'plum': (0.8666666666666667, 0.6274509803921569, 0.8666666666666667),
1626 'purple': (0.6274509803921569, 0.12549019607843137, 0.9411764705882353),
1627 'red': (1.0, 0.0, 0.0),
1628 'rosy brown': (0.7372549019607844, 0.5607843137254902, 0.5607843137254902),
1629 'salmon': (0.9803921568627451, 0.5019607843137255, 0.4470588235294118),
1630 'sandy brown': (0.9568627450980393, 0.6431372549019608, 0.3764705882352941),
1631 'sea green': (0.1803921568627451, 0.5450980392156862, 0.3411764705882353),
1632 'sienna': (0.6274509803921569, 0.3215686274509804, 0.17647058823529413),
1633 'sky blue': (0.5294117647058824, 0.807843137254902, 0.9215686274509803),
1634 'slate gray': (0.4392156862745098, 0.5019607843137255, 0.5647058823529412),
1635 'spring green': (0.0, 1.0, 0.4980392156862745),
1636 'steel blue': (0.27450980392156865, 0.5098039215686274, 0.7058823529411765),
1637 'tan': (0.8235294117647058, 0.7058823529411765, 0.5490196078431373),
1638 'turquoise': (0.25098039215686274, 0.8784313725490196, 0.8156862745098039),
1639 'violet red': (0.8156862745098039, 0.12549019607843137, 0.5647058823529412),
1640 'white': (1.0, 1.0, 1.0),
1641 'yellow': (1.0, 1.0, 0.0)}
1642 if colorname.startswith(
'#'):
1643 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()
static Surface setup_particle(Model *m, ParticleIndex pi)
Add uncertainty to a particle.
void add_particle(RMF::FileHandle fh, Particle *hs)
static bool get_is_setup(const IMP::ParticleAdaptor &p)
static bool get_is_setup(const IMP::ParticleAdaptor &p)
GenericHierarchies get_leaves(Hierarchy mhd)
Get all the leaves of the bit of hierarchy.
This class initializes the root node of the global IMP.atom.Hierarchy.
def add_imp_provenance
Tag the given particle as being created by the current version of IMP.
Vector3D get_random_vector_in(const Cylinder3D &c)
Generate a random vector in a cylinder with uniform density.
Bond create_bond(Bonded a, Bonded b, Bond o)
Connect the two wrapped particles by a custom bond.
Object used to hold a set of restraints.
Stores a named protein chain.
static bool get_is_setup(Model *m, ParticleIndex pi)
def add_software_provenance
Tag the given particle with the software used to create it.
A decorator for keeping track of copies of a molecule.
ParticleIndexPairs get_indexes(const ParticlePairsTemp &ps)
The standard decorator for manipulating molecular structures.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
A decorator for a particle representing an atom.
void transform(XYZ a, const algebra::Transformation3D &tr)
Apply a transformation to the particle.
static Bonded setup_particle(Model *m, ParticleIndex pi)
void load_frame(RMF::FileConstHandle file, RMF::FrameID frame)
Load the given RMF frame into the state of the linked objects.
A decorator for a particle with x,y,z coordinates.
static Scale setup_particle(Model *m, ParticleIndex pi)
A decorator for a particle that is part of a rigid body but not rigid.
Find all nearby pairs by testing all pairs.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
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 ...