3 """@namespace IMP.pmi1.tools
4 Miscellaneous utilities.
7 from __future__
import print_function
14 from collections.abc
import MutableSet
16 from collections
import MutableSet
18 from math
import log,pi,sqrt,exp
23 from time
import process_time
25 from time
import clock
as process_time
28 from collections
import defaultdict, OrderedDict
31 def _add_pmi_provenance(p):
32 """Tag the given particle as being created by the current version of PMI."""
36 location=
"https://integrativemodeling.org")
39 def _get_restraint_set_keys():
40 if not hasattr(_get_restraint_set_keys,
'pmi_rs_key'):
41 _get_restraint_set_keys.pmi_rs_key =
IMP.ModelKey(
"PMI restraints")
42 _get_restraint_set_keys.rmf_rs_key =
IMP.ModelKey(
"RMF restraints")
43 return (_get_restraint_set_keys.pmi_rs_key,
44 _get_restraint_set_keys.rmf_rs_key)
46 def _add_restraint_sets(model, mk, mk_rmf):
49 model.add_data(mk, rs)
50 model.add_data(mk_rmf, rs_rmf)
54 """Add a PMI restraint to the model.
55 Since Model.add_restraint() no longer exists (in modern IMP restraints
56 should be added to a ScoringFunction instead) store them instead in
57 a RestraintSet, and keep a reference to it in the Model.
59 If `add_to_rmf` is True, also add the restraint to a separate list
60 of restraints that will be written out to RMF files (by default, most
61 PMI restraints are not)."""
62 mk, mk_rmf = _get_restraint_set_keys()
63 if model.get_has_data(mk):
64 rs = IMP.RestraintSet.get_from(model.get_data(mk))
65 rs_rmf = IMP.RestraintSet.get_from(model.get_data(mk_rmf))
67 rs, rs_rmf = _add_restraint_sets(model, mk, mk_rmf)
68 rs.add_restraint(restraint)
70 rs_rmf.add_restraint(restraint)
73 """Get a RestraintSet containing all PMI restraints added to the model.
74 If `rmf` is True, return only the subset of these restraints that
75 should be written out to RMF files."""
76 mk, mk_rmf = _get_restraint_set_keys()
77 if not model.get_has_data(mk):
78 print(
"WARNING: no restraints added to model yet")
79 _add_restraint_sets(model, mk, mk_rmf)
81 return IMP.RestraintSet.get_from(model.get_data(mk_rmf))
83 return IMP.RestraintSet.get_from(model.get_data(mk))
86 """Collect timing information.
87 Add an instance of this class to outputobjects to get timing information
92 @param isdelta if True (the default) then report the time since the
93 last use of this class; if False, report cumulative time."""
94 self.starttime = process_time()
96 self.isdelta = isdelta
98 def set_label(self, labelstr):
101 def get_output(self):
104 newtime = process_time()
105 output[
"Stopwatch_" + self.label +
"_delta_seconds"] \
106 = str(newtime - self.starttime)
107 self.starttime = newtime
109 output[
"Stopwatch_" + self.label +
"_elapsed_seconds"] \
110 = str(process_time() - self.starttime)
114 class SetupNuisance(object):
116 def __init__(self, m, initialvalue, minvalue, maxvalue, isoptimized=True):
120 nuisance.set_lower(minvalue)
122 nuisance.set_upper(maxvalue)
125 nuisance.set_is_optimized(nuisance.get_nuisance_key(), isoptimized)
126 self.nuisance = nuisance
128 def get_particle(self):
132 class SetupWeight(object):
134 def __init__(self, m, isoptimized=True):
137 self.weight.set_weights_are_optimized(isoptimized)
139 def get_particle(self):
143 class SetupSurface(object):
145 def __init__(self, m, center, normal, isoptimized=True):
148 self.surface.set_coordinates_are_optimized(isoptimized)
149 self.surface.set_normal_is_optimized(isoptimized)
151 def get_particle(self):
156 "If you use this class please let the PMI developers know.")
157 class ParticleToSampleFilter(object):
158 def __init__(self, sampled_objects):
159 self.sampled_objects=sampled_objects
162 def add_filter(self,filter_string):
163 self.filters.append(filter_string)
165 def get_particles_to_sample(self):
166 particles_to_sample={}
167 for so
in self.sampled_objects:
168 ps_dict=so.get_particles_to_sample()
170 for f
in self.filters:
172 if key
not in particles_to_sample:
173 particles_to_sample[key]=ps_dict[key]
175 particles_to_sample[key]+=ps_dict[key]
176 return particles_to_sample
178 class ParticleToSampleList(object):
182 self.dictionary_particle_type = {}
183 self.dictionary_particle_transformation = {}
184 self.dictionary_particle_name = {}
191 particle_transformation,
193 if not particle_type
in [
"Rigid_Bodies",
"Floppy_Bodies",
"Nuisances",
"X_coord",
"Weights",
"Surfaces"]:
194 raise TypeError(
"not the right particle type")
196 self.dictionary_particle_type[particle] = particle_type
197 if particle_type ==
"Rigid_Bodies":
198 if type(particle_transformation) == tuple
and len(particle_transformation) == 2
and type(particle_transformation[0]) == float
and type(particle_transformation[1]) == float:
199 self.dictionary_particle_transformation[
200 particle] = particle_transformation
201 self.dictionary_particle_name[particle] = name
203 raise TypeError(
"ParticleToSampleList: not the right transformation format for Rigid_Bodies, should be a tuple of floats")
204 elif particle_type ==
"Surfaces":
205 if type(particle_transformation) == tuple
and len(particle_transformation) == 3
and all(isinstance(x, float)
for x
in particle_transformation):
206 self.dictionary_particle_transformation[
207 particle] = particle_transformation
208 self.dictionary_particle_name[particle] = name
210 raise TypeError(
"ParticleToSampleList: not the right transformation format for Surfaces, should be a tuple of floats")
212 if type(particle_transformation) == float:
213 self.dictionary_particle_transformation[
214 particle] = particle_transformation
215 self.dictionary_particle_name[particle] = name
217 raise TypeError(
"ParticleToSampleList: not the right transformation format, should be a float")
219 def get_particles_to_sample(self):
221 for particle
in self.dictionary_particle_type:
222 key = self.dictionary_particle_type[
223 particle] +
"ParticleToSampleList_" + self.dictionary_particle_name[particle] +
"_" + self.label
226 self.dictionary_particle_transformation[particle])
231 class Variance(object):
233 def __init__(self, model, tau, niter, prot, th_profile, write_data=False):
236 self.write_data = write_data
240 particles = IMP.atom.get_by_type(prot, IMP.atom.ATOM_TYPE)
241 self.particles = particles
243 self.refpos = [
IMP.core.XYZ(p).get_coordinates()
for p
in particles]
244 self.model_profile = th_profile
246 def perturb_particles(self, perturb=True):
247 for i, p
in enumerate(self.particles):
248 newpos = array(self.refpos[i])
250 newpos += random.normal(0, self.tau, 3)
254 def get_profile(self):
255 model_profile = self.model_profile
256 p = model_profile.calculate_profile(self.particles, IMP.saxs.CA_ATOMS)
257 return array([model_profile.get_intensity(i)
for i
in
258 range(model_profile.size())])
260 def init_variances(self):
262 N = self.model_profile.size()
263 a = self.profiles[0][:]
265 self.V = self.m * self.m.T
266 self.normm = linalg.norm(self.m)
267 self.normV = linalg.norm(self.V)
269 def update_variances(self):
270 a = matrix(self.profiles[-1])
271 n = float(len(self.profiles))
272 self.m = a.T / n + (n - 1) / n * self.m
273 self.V = a.T * a + self.V
274 self.oldnormm = self.normm
275 self.oldnormV = self.normV
276 self.normm = linalg.norm(self.m)
277 self.normV = linalg.norm(self.V)
278 self.diffm = (self.oldnormm - self.normm) / self.oldnormm
279 self.diffV = (self.oldnormV - self.normV) / self.oldnormV
281 def get_direct_stats(self, a):
286 for q, I
in enumerate(prof):
291 Sigma = (matrix(a - m))
292 Sigma = Sigma.T * Sigma / (nprof - 1)
293 mi = matrix(diag(1. / m))
294 Sigmarel = mi.T * Sigma * mi
295 return m, V, Sigma, Sigmarel
297 def store_data(self):
298 if not os.path.isdir(
'data'):
300 profiles = matrix(self.profiles)
301 self.directm, self.directV, self.Sigma, self.Sigmarel = \
302 self.get_direct_stats(array(profiles))
303 directV = self.directV
305 save(
'data/profiles', profiles)
307 fl = open(
'data/profiles.dat',
'w')
308 for i, l
in enumerate(array(profiles).T):
309 self.model_profile.get_q(i)
312 fl.write(
'%s ' % (k - self.directm[i]))
315 fl = open(
'data/profiles_rel.dat',
'w')
316 for i, l
in enumerate(array(profiles).T):
317 self.model_profile.get_q(i)
320 fl.write(
'%s ' % ((k - self.directm[i]) / self.directm[i]))
322 save(
'data/m', self.directm)
323 save(
'data/V', self.directV)
325 save(
'data/Sigma', Sigma)
327 fl = open(
'data/Sigma.dat',
'w')
328 model_profile = self.model_profile
329 for i
in range(model_profile.size()):
330 qi = model_profile.get_q(i)
331 for j
in range(model_profile.size()):
332 qj = model_profile.get_q(j)
333 vij = self.Sigma[i, j]
334 fl.write(
'%s %s %s\n' % (qi, qj, vij))
337 fl = open(
'data/eigenvals',
'w')
338 for i
in linalg.eigvalsh(Sigma):
340 Sigmarel = self.Sigmarel
341 save(
'data/Sigmarel', Sigmarel)
343 fl = open(
'data/Sigmarel.dat',
'w')
344 model_profile = self.model_profile
345 for i
in range(model_profile.size()):
346 qi = model_profile.get_q(i)
347 for j
in range(model_profile.size()):
348 qj = model_profile.get_q(j)
349 vij = self.Sigmarel[i, j]
350 fl.write(
'%s %s %s\n' % (qi, qj, vij))
353 fl = open(
'data/eigenvals_rel',
'w')
354 for i
in linalg.eigvalsh(Sigmarel):
357 fl = open(
'data/mean.dat',
'w')
358 for i
in range(len(self.directm)):
359 qi = self.model_profile.get_q(i)
361 fl.write(
'%s ' % self.directm[i])
362 fl.write(
'%s ' % sqrt(self.Sigma[i, i]))
365 def try_chol(self, jitter):
368 linalg.cholesky(Sigma + matrix(eye(len(Sigma))) * jitter)
369 except linalg.LinAlgError:
370 print(
"Decomposition failed with jitter =", jitter)
372 print(
"Successful decomposition with jitter =", jitter)
375 self.profiles = [self.get_profile()]
377 for n
in range(self.niter):
378 self.perturb_particles()
379 self.profiles.append(self.get_profile())
391 def get_cov(self, relative=True):
401 number_of_cross_links=10,
402 ambiguity_probability=0.1,
403 confidence_score_range=[0,100],
404 avoid_same_particles=
False):
405 '''Return a random cross-link dataset as a string.
406 Every line is a residue pair, together with UniqueIdentifier
409 residue_pairs=get_random_residue_pairs(representation, resolution, number_of_cross_links, avoid_same_particles=avoid_same_particles)
412 cmin=float(min(confidence_score_range))
413 cmax=float(max(confidence_score_range))
417 for (name1, r1, name2, r2)
in residue_pairs:
418 if random.random() > ambiguity_probability:
420 score=random.random()*(cmax-cmin)+cmin
421 dataset+=str(name1)+
" "+str(name2)+
" "+str(r1)+
" "+str(r2)+
" "+str(score)+
" "+str(unique_identifier)+
"\n"
428 def get_cross_link_data(directory, filename, dist, omega, sigma,
429 don=
None, doff=
None, prior=0, type_of_profile=
"gofr"):
431 (distmin, distmax, ndist) = dist
432 (omegamin, omegamax, nomega) = omega
433 (sigmamin, sigmamax, nsigma) = sigma
436 with open(filen)
as xlpot:
437 dictionary = ast.literal_eval(xlpot.readline())
439 xpot = dictionary[directory][filename][
"distance"]
440 pot = dictionary[directory][filename][type_of_profile]
442 dist_grid =
get_grid(distmin, distmax, ndist,
False)
443 omega_grid = get_log_grid(omegamin, omegamax, nomega)
444 sigma_grid = get_log_grid(sigmamin, sigmamax, nsigma)
446 if not don
is None and not doff
is None:
468 def get_cross_link_data_from_length(length, xxx_todo_changeme3, xxx_todo_changeme4, xxx_todo_changeme5):
469 (distmin, distmax, ndist) = xxx_todo_changeme3
470 (omegamin, omegamax, nomega) = xxx_todo_changeme4
471 (sigmamin, sigmamax, nsigma) = xxx_todo_changeme5
473 dist_grid =
get_grid(distmin, distmax, ndist,
False)
474 omega_grid = get_log_grid(omegamin, omegamax, nomega)
475 sigma_grid = get_log_grid(sigmamin, sigmamax, nsigma)
481 def get_grid(gmin, gmax, ngrid, boundaries):
483 dx = (gmax - gmin) / float(ngrid)
484 for i
in range(0, ngrid + 1):
485 if(
not boundaries
and i == 0):
487 if(
not boundaries
and i == ngrid):
489 grid.append(gmin + float(i) * dx)
495 def get_log_grid(gmin, gmax, ngrid):
497 for i
in range(0, ngrid + 1):
498 grid.append(gmin * exp(float(i) / ngrid * log(gmax / gmin)))
506 example '"{ID_Score}" > 28 AND "{Sample}" ==
507 "%10_1%" OR ":Sample}" == "%10_2%" OR ":Sample}"
508 == "%10_3%" OR ":Sample}" == "%8_1%" OR ":Sample}" == "%8_2%"'
511 import pyparsing
as pp
513 operator = pp.Regex(
">=|<=|!=|>|<|==|in").setName(
"operator")
514 value = pp.QuotedString(
516 r"[+-]?\d+(:?\.\d*)?(:?[eE][+-]?\d+)?")
517 identifier = pp.Word(pp.alphas, pp.alphanums +
"_")
518 comparison_term = identifier | value
519 condition = pp.Group(comparison_term + operator + comparison_term)
521 expr = pp.operatorPrecedence(condition, [
522 (
"OR", 2, pp.opAssoc.LEFT, ),
523 (
"AND", 2, pp.opAssoc.LEFT, ),
526 parsedstring = str(expr.parseString(inputstring)) \
532 .replace(
"{",
"float(entry['") \
533 .replace(
"}",
"'])") \
534 .replace(
":",
"str(entry['") \
535 .replace(
"}",
"'])") \
536 .replace(
"AND",
"and") \
541 def open_file_or_inline_text(filename):
543 fl = open(filename,
"r")
545 fl = filename.split(
"\n")
552 for i
in range(0, len(prot0) - 1):
553 for j
in range(i + 1, len(prot0)):
556 drmsd += (dist0 - dist1) ** 2
558 return math.sqrt(drmsd / npairs)
563 def get_ids_from_fasta_file(fastafile):
565 with open(fastafile)
as ff:
574 this function works with plain hierarchies, as read from the pdb,
575 no multi-scale hierarchies
582 atom_type=IMP.atom.AT_CA)
590 print(
"get_closest_residue_position: exiting while loop without result")
592 p = sel.get_selected_particles()
597 print(
"get_closest_residue_position: got NO residues for hierarchy %s and residue %i" % (hier, resindex))
598 raise Exception(
"get_closest_residue_position: got NO residues for hierarchy %s and residue %i" % (
601 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])))
605 Get the particle of the terminal residue at the GIVEN resolution
606 (NOTE: not the closest resolution!).
607 To get the terminal residue at the closest resolution use:
608 particles=IMP.pmi1.tools.select_by_tuple(representation,molecule_name)
609 particles[0] and particles[-1] will be the first and last particles
610 corresponding to the two termini.
611 It is needed for instance to determine the last residue of a pdb.
612 @param hier hierarchy containing the terminal residue
613 @param terminus either 'N' or 'C'
614 @param resolution resolution to use.
620 resolution=resolution,
627 if termresidue
is None:
628 termresidue = max(residues)
630 elif max(residues) >= termresidue:
631 termresidue = max(residues)
633 elif terminus ==
"N":
634 if termresidue
is None:
635 termresidue = min(residues)
637 elif min(residues) <= termresidue:
638 termresidue = min(residues)
641 raise ValueError(
"terminus argument should be either N or C")
646 """Get XYZ coordinates of the terminal residue at the GIVEN resolution"""
652 Return the residue index gaps and contiguous segments in the hierarchy.
654 @param hierarchy hierarchy to examine
655 @param start first residue index
656 @param end last residue index
658 @return A list of lists of the form
659 [[1,100,"cont"],[101,120,"gap"],[121,200,"cont"]]
662 for n, rindex
in enumerate(range(start, end + 1)):
664 atom_type=IMP.atom.AT_CA)
666 if len(sel.get_selected_particles()) == 0:
670 rindexcont = start - 1
671 if rindexgap == rindex - 1:
677 gaps.append([rindex, rindex,
"gap"])
683 rindexgap = start - 1
685 if rindexcont == rindex - 1:
692 gaps.append([rindex, rindex,
"cont"])
703 def set_map_element(self, xvalue, yvalue):
704 self.map[xvalue] = yvalue
706 def get_map_element(self, invalue):
707 if type(invalue) == float:
711 dist = (invalue - x) * (invalue - x)
720 return self.map[minx]
721 elif type(invalue) == str:
722 return self.map[invalue]
724 raise TypeError(
"wrong type for map")
730 selection_arguments=
None,
732 name_is_ambiguous=
False,
736 representation_type=
None):
738 this function uses representation=SimplifiedModel
739 it returns the corresponding selected particles
740 representation_type="Beads", "Res:X", "Densities", "Representation", "Molecule"
743 if resolution
is None:
745 resolution_particles =
None
746 hierarchies_particles =
None
747 names_particles =
None
748 residue_range_particles =
None
749 residue_particles =
None
750 representation_type_particles =
None
752 if not resolution
is None:
753 resolution_particles = []
754 hs = representation.get_hierarchies_at_given_resolution(resolution)
758 if not hierarchies
is None:
759 hierarchies_particles = []
760 for h
in hierarchies:
765 if name_is_ambiguous:
766 for namekey
in representation.hier_dict:
769 representation.hier_dict[namekey])
770 elif name
in representation.hier_dict:
773 print(
"select: component %s is not there" % name)
775 if not first_residue
is None and not last_residue
is None:
777 residue_indexes=range(first_residue, last_residue + 1))
779 for p
in sel.get_selected_particles()]
781 if not residue
is None:
784 for p
in sel.get_selected_particles()]
786 if not representation_type
is None:
787 representation_type_particles = []
788 if representation_type ==
"Molecule":
789 for name
in representation.hier_representation:
790 for repr_type
in representation.hier_representation[name]:
791 if repr_type ==
"Beads" or "Res:" in repr_type:
792 h = representation.hier_representation[name][repr_type]
795 elif representation_type ==
"PDB":
796 for name
in representation.hier_representation:
797 for repr_type
in representation.hier_representation[name]:
798 if repr_type ==
"Res:" in repr_type:
799 h = representation.hier_representation[name][repr_type]
803 for name
in representation.hier_representation:
804 h = representation.hier_representation[
809 selections = [hierarchies_particles, names_particles,
810 residue_range_particles, residue_particles, representation_type_particles]
812 if resolution
is None:
813 selected_particles = set(allparticles)
815 selected_particles = set(resolution_particles)
819 selected_particles = (set(s) & selected_particles)
821 return list(selected_particles)
828 name_is_ambiguous=
False):
829 if isinstance(tupleselection, tuple)
and len(tupleselection) == 3:
831 name=tupleselection[2],
832 first_residue=tupleselection[0],
833 last_residue=tupleselection[1],
834 name_is_ambiguous=name_is_ambiguous)
835 elif isinstance(tupleselection, str):
838 name_is_ambiguous=name_is_ambiguous)
840 raise ValueError(
'you passed something bad to select_by_tuple()')
842 particles = IMP.pmi1.tools.sort_by_residues(particles)
847 """New tuple format: molname OR (start,stop,molname,copynum,statenum)
848 Copy and state are optional. Can also use 'None' for them which will get all.
849 You can also pass -1 for stop which will go to the end.
850 Returns the particles
853 kwds[
'resolution'] = resolution
854 if type(tuple_selection)
is str:
855 kwds[
'molecule'] = tuple_selection
856 elif type(tuple_selection)
is tuple:
857 rbegin = tuple_selection[0]
858 rend = tuple_selection[1]
859 kwds[
'molecule'] = tuple_selection[2]
861 copynum = tuple_selection[3]
862 if copynum
is not None:
863 kwds[
'copy_index'] = copynum
867 statenum = tuple_selection[4]
868 if statenum
is not None:
869 kwds[
'state_index'] = statenum
876 residue_indexes=range(1,rbegin),
878 return s.get_selected_particles()
880 kwds[
'residue_indexes'] = range(rbegin,rend+1)
882 return s.get_selected_particles()
886 def get_db_from_csv(csvfilename):
889 with open(csvfilename)
as fh:
890 csvr = csv.DictReader(fh)
897 """Store the representations for a system."""
902 self.root_hierarchy_dict = {}
903 self.preroot_fragment_hierarchy_dict = {}
904 self.particle_to_name = {}
907 def add_name(self, name):
908 if name
not in self.db:
911 def add_residue_number(self, name, resn):
914 if resn
not in self.db[name]:
915 self.db[name][resn] = {}
917 def add_resolution(self, name, resn, resolution):
919 resolution = float(resolution)
921 self.add_residue_number(name, resn)
922 if resolution
not in self.db[name][resn]:
923 self.db[name][resn][resolution] = []
927 resolution = float(resolution)
929 self.add_residue_number(name, resn)
930 self.add_resolution(name, resn, resolution)
931 self.db[name][resn][resolution] += particles
933 (rh, prf) = self.get_root_hierarchy(p)
934 self.root_hierarchy_dict[p] = rh
935 self.preroot_fragment_hierarchy_dict[p] = prf
936 self.particle_to_name[p] = name
937 if self.model
is None:
938 self.model = particles[0].get_model()
944 names = list(self.db.keys())
950 resolution = float(resolution)
951 return self.db[name][resn][resolution]
953 def get_particles_at_closest_resolution(self, name, resn, resolution):
955 resolution = float(resolution)
956 closestres = min(self.get_residue_resolutions(name, resn),
957 key=
lambda x: abs(float(x) - float(resolution)))
958 return self.get_particles(name, resn, closestres)
960 def get_residue_resolutions(self, name, resn):
962 resolutions = list(self.db[name][resn].keys())
966 def get_molecule_resolutions(self, name):
968 for resn
in self.db[name]:
969 resolutions.update(list(self.db[name][resn].keys()))
973 def get_residue_numbers(self, name):
974 residue_numbers = list(self.db[name].keys())
975 residue_numbers.sort()
976 return residue_numbers
978 def get_particles_by_resolution(self, name, resolution):
979 resolution = float(resolution)
981 for resn
in self.get_residue_numbers(name):
982 result = self.get_particles_at_closest_resolution(
986 pstemp = [p
for p
in result
if p
not in particles]
990 def get_all_particles_by_resolution(self, resolution):
991 resolution = float(resolution)
993 for name
in self.get_names():
994 particles += self.get_particles_by_resolution(name, resolution)
997 def get_root_hierarchy(self, particle):
998 prerootfragment = particle
1008 prerootfragment = particle
1014 def get_all_root_hierarchies_by_resolution(self, resolution):
1016 resolution = float(resolution)
1017 particles = self.get_all_particles_by_resolution(resolution)
1019 rh = self.root_hierarchy_dict[p]
1020 if rh
not in hierarchies:
1024 def get_preroot_fragments_by_resolution(self, name, resolution):
1026 resolution = float(resolution)
1027 particles = self.get_particles_by_resolution(name, resolution)
1029 fr = self.preroot_fragment_hierarchy_dict[p]
1030 if fr
not in fragments:
1031 fragments.append(fr)
1034 def show(self, name):
1036 for resn
in self.get_residue_numbers(name):
1038 for resolution
in self.get_residue_resolutions(name, resn):
1039 print(
"----", resolution)
1040 for p
in self.get_particles(name, resn, resolution):
1041 print(
"--------", p.get_name())
1045 '''Return the component name provided a particle and a list of names'''
1047 protname = root.get_name()
1049 while not protname
in list_of_names:
1050 root0 = root.get_parent()
1053 protname = root0.get_name()
1058 if "Beads" in protname:
1061 return (protname, is_a_bead)
1066 Retrieve the residue indexes for the given particle.
1068 The particle must be an instance of Fragment,Residue, Atom or Molecule
1069 or else returns an empty list
1080 resind_tmp=IMP.pmi1.tools.OrderedSet()
1086 resind=list(resind_tmp)
1092 def sort_by_residues(particles):
1095 sorted_particles_residues = sorted(
1097 key=
lambda tup: tup[1])
1098 particles = [p[0]
for p
in sorted_particles_residues]
1102 def get_residue_to_particle_map(particles):
1104 particles = sort_by_residues(particles)
1107 return dict(zip(particles_residues, particles))
1115 """Synchronize data over a parallel run"""
1116 from mpi4py
import MPI
1117 comm = MPI.COMM_WORLD
1118 rank = comm.Get_rank()
1119 number_of_processes = comm.size
1122 comm.send(data, dest=0, tag=11)
1125 for i
in range(1, number_of_processes):
1126 data_tmp = comm.recv(source=i, tag=11)
1127 if type(data) == list:
1129 elif type(data) == dict:
1130 data.update(data_tmp)
1132 raise TypeError(
"data not supported, use list or dictionaries")
1134 for i
in range(1, number_of_processes):
1135 comm.send(data, dest=i, tag=11)
1138 data = comm.recv(source=0, tag=11)
1142 """Synchronize data over a parallel run"""
1143 from mpi4py
import MPI
1144 comm = MPI.COMM_WORLD
1145 rank = comm.Get_rank()
1146 number_of_processes = comm.size
1149 comm.send(data, dest=0, tag=11)
1151 for i
in range(1, number_of_processes):
1152 data_tmp = comm.recv(source=i, tag=11)
1154 data[k]+=data_tmp[k]
1156 for i
in range(1, number_of_processes):
1157 comm.send(data, dest=i, tag=11)
1160 data = comm.recv(source=0, tag=11)
1170 Yield all sublists of length >= lmin and <= lmax
1178 for j
in range(i + 1, n):
1179 if len(l[i:j]) <= lmax
and len(l[i:j]) >= lmin:
1183 def flatten_list(l):
1184 return [item
for sublist
in l
for item
in sublist]
1188 """ Yield successive length-sized chunks from a list.
1190 for i
in range(0, len(list), length):
1191 yield list[i:i + length]
1194 def chunk_list_into_segments(seq, num):
1196 avg = len(seq) / float(num)
1200 while last < len(seq):
1201 out.append(seq[int(last):int(last + avg)])
1209 ''' This class stores integers
1210 in ordered compact lists eg:
1212 the methods help splitting and merging the internal lists
1214 s=Segments([1,2,3]) is [[1,2,3]]
1215 s.add(4) is [[1,2,3,4]] (add right)
1216 s.add(3) is [[1,2,3,4]] (item already existing)
1217 s.add(7) is [[1,2,3,4],[7]] (new list)
1218 s.add([8,9]) is [[1,2,3,4],[7,8,9]] (add item right)
1219 s.add([5,6]) is [[1,2,3,4,5,6,7,8,9]] (merge)
1220 s.remove(3) is [[1,2],[4,5,6,7,8,9]] (split)
1225 '''index can be a integer or a list of integers '''
1226 if type(index)
is int:
1228 elif type(index)
is list:
1229 self.segs=[[index[0]]]
1234 '''index can be a integer or a list of integers '''
1235 if type(index)
is int:
1238 for n,s
in enumerate(self.segs):
1246 if mergeright
is None and mergeleft
is None:
1247 self.segs.append([index])
1248 if not mergeright
is None and mergeleft
is None:
1249 self.segs[mergeright].append(index)
1250 if not mergeleft
is None and mergeright
is None:
1251 self.segs[mergeleft]=[index]+self.segs[mergeleft]
1252 if not mergeleft
is None and not mergeright
is None:
1253 self.segs[mergeright]=self.segs[mergeright]+[index]+self.segs[mergeleft]
1254 del self.segs[mergeleft]
1256 for n
in range(len(self.segs)):
1259 self.segs.sort(key=
lambda tup: tup[0])
1261 elif type(index)
is list:
1266 '''index can be a integer'''
1267 for n,s
in enumerate(self.segs):
1274 i=self.segs[n].index(index)
1276 self.segs.append(s[i+1:])
1277 for n
in range(len(self.segs)):
1279 if len(self.segs[n])==0:
1281 self.segs.sort(key=
lambda tup: tup[0])
1284 ''' Returns a flatten list '''
1285 return [item
for sublist
in self.segs
for item
in sublist]
1289 for seg
in self.segs:
1290 ret_tmp+=str(seg[0])+
"-"+str(seg[-1])+
","
1291 ret=ret_tmp[:-1]+
"]"
1303 Apply a translation to a hierarchy along the input vector.
1307 if type(translation_vector) == list:
1326 def translate_hierarchies(hierarchies, translation_vector):
1327 for h
in hierarchies:
1331 def translate_hierarchies_to_reference_frame(hierarchies):
1336 for h
in hierarchies:
1346 IMP.pmi1.tools.translate_hierarchies(hierarchies, (-xc, -yc, -zc))
1353 def normal_density_function(expected_value, sigma, x):
1355 1 / math.sqrt(2 * math.pi) / sigma *
1356 math.exp(-(x - expected_value) ** 2 / 2 / sigma / sigma)
1360 def log_normal_density_function(expected_value, sigma, x):
1362 1 / math.sqrt(2 * math.pi) / sigma / x *
1363 math.exp(-(math.log(x / expected_value) ** 2 / 2 / sigma / sigma))
1367 def get_random_residue_pairs(representation, resolution,
1370 avoid_same_particles=
False,
1375 names=list(representation.hier_dict.keys())
1378 prot = representation.hier_dict[name]
1379 particles +=
select(representation,name=name,resolution=resolution)
1380 random_residue_pairs = []
1381 while len(random_residue_pairs)<=number:
1382 p1 = random.choice(particles)
1383 p2 = random.choice(particles)
1384 if max_distance
is not None and \
1389 if r1==r2
and avoid_same_particles:
continue
1390 name1 = representation.get_prot_name_from_particle(p1)
1391 name2 = representation.get_prot_name_from_particle(p2)
1392 random_residue_pairs.append((name1, r1, name2, r2))
1394 return random_residue_pairs
1397 def get_random_data_point(
1403 begin_end_nbins_tuple,
1408 begin = begin_end_nbins_tuple[0]
1409 end = begin_end_nbins_tuple[1]
1410 nbins = begin_end_nbins_tuple[2]
1413 fmod_grid =
get_grid(begin, end, nbins,
True)
1415 fmod_grid = get_log_grid(begin, end, nbins)
1422 for i
in range(0, ntrials):
1423 a.append([random.random(),
True])
1426 for j
in range(1, len(fmod_grid)):
1428 fjm1 = fmod_grid[j - 1]
1432 pj = normal_density_function(expected_value, sigma, fj)
1433 pjm1 = normal_density_function(expected_value, sigma, fjm1)
1435 pj = log_normal_density_function(expected_value, sigma, fj)
1436 pjm1 = log_normal_density_function(expected_value, sigma, fjm1)
1438 norm += (pj + pjm1) / 2.0 * df
1444 for i
in range(len(cumul)):
1447 if (aa[0] <= cumul[i] / norm
and aa[1]):
1448 random_points.append(
1449 int(fmod_grid[i] / sensitivity) * sensitivity)
1453 random_points = [expected_value] * ntrials
1455 for i
in range(len(random_points)):
1456 if random.random() < outlierprob:
1457 a = random.uniform(begin, end)
1458 random_points[i] = int(a / sensitivity) * sensitivity
1459 print(random_points)
1461 for i in range(ntrials):
1462 if random.random() > OUTLIERPROB_:
1463 r=truncnorm.rvs(0.0,1.0,expected_value,BETA_)
1464 if r>1.0: print r,expected_value,BETA_
1467 random_points.append(int(r/sensitivity)*sensitivity)
1472 for r
in random_points:
1476 rmean /= float(ntrials)
1477 rmean2 /= float(ntrials)
1478 stddev = math.sqrt(max(rmean2 - rmean * rmean, 0.))
1479 return rmean, stddev
1481 def print_multicolumn(list_of_strings, ncolumns=2, truncate=40):
1487 for i
in range(len(l) % cols):
1490 split = [l[i:i + len(l) / cols]
for i
in range(0, len(l), len(l) / cols)]
1491 for row
in zip(*split):
1492 print(
"".join(str.ljust(i, truncate)
for i
in row))
1495 '''Change color code to hexadecimal to rgb'''
1497 self._NUMERALS =
'0123456789abcdefABCDEF'
1498 self._HEXDEC = dict((v, int(v, 16))
for v
in (x+y
for x
in self._NUMERALS
for y
in self._NUMERALS))
1499 self.LOWERCASE, self.UPPERCASE =
'x',
'X'
1501 def rgb(self,triplet):
1502 return float(self._HEXDEC[triplet[0:2]]), float(self._HEXDEC[triplet[2:4]]), float(self._HEXDEC[triplet[4:6]])
1504 def triplet(self,rgb, lettercase=None):
1505 if lettercase
is None: lettercase=self.LOWERCASE
1506 return format(rgb[0]<<16 | rgb[1]<<8 | rgb[2],
'06'+lettercase)
1510 class OrderedSet(MutableSet):
1512 def __init__(self, iterable=None):
1514 end += [
None, end, end]
1516 if iterable
is not None:
1520 return len(self.map)
1522 def __contains__(self, key):
1523 return key
in self.map
1526 if key
not in self.map:
1529 curr[2] = end[1] = self.map[key] = [key, curr, end]
1531 def discard(self, key):
1533 key, prev, next = self.map.pop(key)
1540 while curr
is not end:
1544 def __reversed__(self):
1547 while curr
is not end:
1551 def pop(self, last=True):
1553 raise KeyError(
'set is empty')
1555 key = self.end[1][0]
1557 key = self.end[2][0]
1563 return '%s()' % (self.__class__.__name__,)
1564 return '%s(%r)' % (self.__class__.__name__, list(self))
1566 def __eq__(self, other):
1567 if isinstance(other, OrderedSet):
1568 return len(self) == len(other)
and list(self) == list(other)
1569 return set(self) == set(other)
1573 """Store objects in order they were added, but with default type.
1574 Source: http://stackoverflow.com/a/4127426/2608793
1576 def __init__(self, *args, **kwargs):
1578 self.default_factory =
None
1580 if not (args[0]
is None or callable(args[0])):
1581 raise TypeError(
'first argument must be callable or None')
1582 self.default_factory = args[0]
1584 super(OrderedDefaultDict, self).__init__(*args, **kwargs)
1586 def __missing__ (self, key):
1587 if self.default_factory
is None:
1589 self[key] = default = self.default_factory()
1592 def __reduce__(self):
1593 args = (self.default_factory,)
if self.default_factory
else ()
1594 return self.__class__, args,
None,
None, self.iteritems()
1599 """Extract frame from RMF file and fill coordinates. Must be identical topology.
1600 @param hier The (System) hierarchy to fill (e.g. after you've built it)
1601 @param rmf_fn The file to extract from
1602 @param frame_num The frame number to extract
1604 rh = RMF.open_rmf_file_read_only(rmf_fn)
1612 selection_tuple=
None,
1613 warn_about_slices=
True):
1614 """Adapt things for PMI (degrees of freedom, restraints, ...)
1615 Returns list of list of hierarchies, separated into Molecules if possible.
1616 The input can be a list, or a list of lists (iterable of ^1 or iterable of ^2)
1617 (iterable of ^2) Hierarchy -> returns input as list of list of hierarchies,
1618 only one entry, not grouped by molecules.
1619 (iterable of ^2) PMI::System/State/Molecule/TempResidue ->
1620 returns residue hierarchies, grouped in molecules, at requested resolution
1621 @param stuff Can be one of the following inputs:
1622 IMP Hierarchy, PMI System/State/Molecule/TempResidue, or a list/set (of list/set) of them.
1623 Must be uniform input, however. No mixing object types.
1624 @param pmi_resolution For selecting, only does it if you pass PMI objects. Set it to "all"
1625 if you want all resolutions!
1626 @param flatten Set to True if you just want all hierarchies in one list.
1627 @param warn_about_slices Print a warning if you are requesting only part of a bead.
1628 Sometimes you just don't care!
1629 \note since this relies on IMP::atom::Selection, this will not return any objects if they weren't built!
1630 But there should be no problem if you request unbuilt residues, they should be ignored.
1636 if hasattr(stuff,
'__iter__'):
1642 if all(hasattr(el,
'__iter__')
for el
in thelist):
1643 thelist = [i
for sublist
in thelist
for i
in sublist]
1644 elif any(hasattr(el,
'__iter__')
for el
in thelist):
1645 raise Exception(
'input_adaptor: input_object must be a list or a list of lists')
1654 except NotImplementedError:
1660 is_temp_residue=
False
1665 if is_system
or is_state
or is_molecule
or is_temp_residue:
1670 for system
in stuff:
1671 for state
in system.get_states():
1672 mdict = state.get_molecules()
1673 for molname
in mdict:
1674 for copy
in mdict[molname]:
1675 indexes_per_mol[copy] += [r.get_index()
for r
in copy.get_residues()]
1678 mdict = state.get_molecules()
1679 for molname
in mdict:
1680 for copy
in mdict[molname]:
1681 indexes_per_mol[copy] += [r.get_index()
for r
in copy.get_residues()]
1683 for molecule
in stuff:
1684 indexes_per_mol[molecule] += [r.get_index()
for r
in molecule.get_residues()]
1685 elif is_temp_residue:
1686 for tempres
in stuff:
1687 indexes_per_mol[tempres.get_molecule()].append(tempres.get_index())
1688 for mol
in indexes_per_mol:
1689 if pmi_resolution==
'all':
1693 residue_indexes=indexes_per_mol[mol])
1696 resolution=pmi_resolution,
1697 residue_indexes=indexes_per_mol[mol])
1698 ps = sel.get_selected_particles()
1701 if warn_about_slices:
1702 rset = set(indexes_per_mol[mol])
1706 if not fset <= rset:
1713 resbreak = maxf
if minf==minset
else minset-1
1714 print(
'WARNING: You are trying to select only part of the bead %s:%i-%i.\n'
1715 'The residues you requested are %i-%i. You can fix this by:\n'
1716 '1) requesting the whole bead/none of it or\n'
1717 '2) break the bead up by passing bead_extra_breaks=[\'%i\'] in '
1718 'molecule.add_representation()'
1719 %(mol.get_name(),minset,maxset,minf,maxf,resbreak))
1724 if pmi_resolution==
'all':
1732 hier_list = [hier_list]
1734 raise Exception(
'input_adaptor: you passed something of wrong type or a list with mixed types')
1736 if flatten
and pmi_input:
1737 return [h
for sublist
in hier_list
for h
in sublist]
1743 """Returns sequence-sorted segments array, each containing the first particle
1744 the last particle and the first residue index."""
1746 from operator
import itemgetter
1749 raise Exception(
"IMP.pmi1.tools.get_sorted_segments: only pass stuff from one Molecule, please")
1765 SortedSegments.append((start, end, startres))
1766 SortedSegments = sorted(SortedSegments, key=itemgetter(2))
1767 return SortedSegments
1770 """Decorate the sequence-consecutive particles from a PMI2 molecule with a bond,
1771 so that they appear connected in the rmf file"""
1773 for x
in range(len(SortedSegments) - 1):
1775 last = SortedSegments[x][1]
1776 first = SortedSegments[x + 1][0]
1778 p1 = last.get_particle()
1779 p2 = first.get_particle()
1792 """This class converts three to one letter codes, and return X for any unknown codes"""
1793 def __init__(self,is_nucleic=False):
1796 threetoone = {
'ALA':
'A',
'ARG':
'R', 'ASN': 'N', 'ASP': 'D',
1797 'CYS':
'C',
'GLU':
'E',
'GLN':
'Q',
'GLY':
'G',
1798 'HIS':
'H',
'ILE':
'I',
'LEU':
'L',
'LYS':
'K',
1799 'MET':
'M',
'PHE':
'F',
'PRO':
'P',
'SER':
'S',
1800 'THR':
'T',
'TRP':
'W',
'TYR':
'Y',
'VAL':
'V',
'UNK':
'X'}
1802 threetoone = {
'ADE':
'A',
'URA':
'U', 'CYT': 'C', 'GUA': 'G',
1803 'THY':
'T',
'UNK':
'X'}
1805 defaultdict.__init__(self,
lambda:
"X", threetoone)
1810 def get_residue_type_from_one_letter_code(code,is_nucleic=False):
1813 for k
in threetoone:
1814 one_to_three[threetoone[k]] = k
1819 """ Just get the leaves from a list of hierarchies """
1820 lvs = list(itertools.chain.from_iterable(
IMP.atom.get_leaves(item)
for item
in list_of_hs))
1827 """Perform selection using the usual keywords but return ALL resolutions (BEADS and GAUSSIANS).
1828 Returns in flat list!
1833 if hier
is not None:
1836 print(
"WARNING: You passed nothing to select_at_all_resolutions()")
1843 raise Exception(
'select_at_all_resolutions: you have to pass an IMP Hierarchy')
1845 raise Exception(
'select_at_all_resolutions: you have to pass an IMP Hierarchy')
1846 if 'resolution' in kwargs
or 'representation_type' in kwargs:
1847 raise Exception(
"don't pass resolution or representation_type to this function")
1849 representation_type=IMP.atom.BALLS,**kwargs)
1851 representation_type=IMP.atom.DENSITIES,**kwargs)
1852 ret |= OrderedSet(selB.get_selected_particles())
1853 ret |= OrderedSet(selD.get_selected_particles())
1862 """Utility to retrieve particles from a hierarchy within a
1863 zone around a set of ps.
1864 @param hier The hierarchy in which to look for neighbors
1865 @param target_ps The particles for zoning
1866 @param sel_zone The maximum distance
1867 @param entire_residues If True, will grab entire residues
1868 @param exclude_backbone If True, will only return sidechain particles
1872 backbone_types=[
'C',
'N',
'CB',
'O']
1873 if exclude_backbone:
1875 for n
in backbone_types])
1876 test_ps = test_sel.get_selected_particles()
1877 nn = IMP.algebra.NearestNeighbor3D([
IMP.core.XYZ(p).get_coordinates()
1880 for target
in target_ps:
1881 zone|=set(nn.get_in_ball(
IMP.core.XYZ(target).get_coordinates(),sel_zone))
1882 zone_ps = [test_ps[z]
for z
in zone]
1887 zone_ps = [h.get_particle()
for h
in final_ps]
1892 """Returns unique objects in original order"""
1896 if not hasattr(hiers,
'__iter__'):
1903 rbs_ordered.append(rb)
1908 rbs_ordered.append(rb)
1912 return rbs_ordered,beads
1915 "This function returns the parent molecule hierarchies of given objects"
1916 stuff=
input_adaptor(input_objects, pmi_resolution=
'all',flatten=
True)
1921 while not (is_molecule
or is_root):
1922 root=IMP.atom.get_root(h)
1929 return list(molecules)
1931 def get_molecules_dictionary(input_objects):
1932 moldict=defaultdict(list)
1935 moldict[name].append(mol)
1941 def get_molecules_dictionary_by_copy(input_objects):
1942 moldict=defaultdict(dict)
1946 moldict[name][c]=mol
1949 def get_selections_dictionary(input_objects):
1950 moldict=IMP.pmi1.tools.get_molecules_dictionary(input_objects)
1951 seldict=defaultdict(list)
1952 for name, mols
in moldict.items():
1958 """Given a list of PMI objects, returns all density hierarchies within
1959 these objects. The output of this function can be inputted into
1960 things such as EM restraints. This function is intended to gather density particles
1961 appended to molecules (and not other hierarchies which might have been appended to the root node directly).
1968 densities+=
IMP.atom.Selection(i,representation_type=IMP.atom.DENSITIES).get_selected_particles()
1972 max_translation=300., max_rotation=2.0 * pi,
1973 avoidcollision_rb=
True, avoidcollision_fb=
False,
1974 cutoff=10.0, niterations=100,
1976 excluded_rigid_bodies=[],
1977 hierarchies_excluded_from_collision=[],
1978 hierarchies_included_in_collision=[],
1980 return_debug=
False):
1981 """Shuffle particles. Used to restart the optimization.
1982 The configuration of the system is initialized by placing each
1983 rigid body and each bead randomly in a box. If `bounding_box` is
1984 specified, the particles are placed inside this box; otherwise, each
1985 particle is displaced by up to max_translation angstroms, and randomly
1986 rotated. Effort is made to place particles far enough from each other to
1987 prevent any steric clashes.
1988 @param objects Can be one of the following inputs:
1989 IMP Hierarchy, PMI System/State/Molecule/TempResidue, or a list/set of them
1990 @param max_translation Max translation (rbs and flexible beads)
1991 @param max_rotation Max rotation (rbs only)
1992 @param avoidcollision_rb check if the particle/rigid body was
1993 placed close to another particle; uses the optional
1994 arguments cutoff and niterations
1995 @param avoidcollision_fb Advanced. Generally you want this False because it's hard to shuffle beads.
1996 @param cutoff Distance less than this is a collision
1997 @param niterations How many times to try avoiding collision
1998 @param bounding_box Only shuffle particles within this box. Defined by ((x1,y1,z1),(x2,y2,z2)).
1999 @param excluded_rigid_bodies Don't shuffle these rigid body objects
2000 @param hierarchies_excluded_from_collision Don't count collision with these bodies
2001 @param hierarchies_included_in_collision Hierarchies that are not shuffled, but should be included in collision calculation (for fixed regions)
2002 @param verbose Give more output
2003 \note Best to only call this function after you've set up degrees of freedom
2004 For debugging purposes, returns: <shuffled indexes>, <collision avoided indexes>
2009 pmi_resolution=
'all',
2012 if len(rigid_bodies)>0:
2013 mdl = rigid_bodies[0].get_model()
2014 elif len(flexible_beads)>0:
2015 mdl = flexible_beads[0].get_model()
2017 raise Exception(
"Could not find any particles in the hierarchy")
2018 if len(rigid_bodies) == 0:
2019 print(
"shuffle_configuration: rigid bodies were not initialized")
2023 gcpf.set_distance(cutoff)
2027 pmi_resolution=
'all',
2031 pmi_resolution=
'all',
2034 collision_excluded_idxs = set([l.get_particle().
get_index()
for h
in collision_excluded_hierarchies \
2037 collision_included_idxs = set([l.get_particle().
get_index()
for h
in collision_included_hierarchies \
2044 all_idxs.append(p.get_particle_index())
2046 collision_excluded_idxs.add(p.get_particle_index())
2048 if bounding_box
is not None:
2049 ((x1, y1, z1), (x2, y2, z2)) = bounding_box
2054 all_idxs = set(all_idxs) | collision_included_idxs
2055 all_idxs = all_idxs - collision_excluded_idxs
2057 print(
'shuffling', len(rigid_bodies),
'rigid bodies')
2058 for rb
in rigid_bodies:
2059 if rb
not in excluded_rigid_bodies:
2061 if avoidcollision_rb:
2062 rb_idxs = set(rb.get_member_particle_indexes()) - \
2063 collision_excluded_idxs
2064 other_idxs = all_idxs - rb_idxs
2070 while niter < niterations:
2071 rbxyz = (rb.get_x(), rb.get_y(), rb.get_z())
2091 debug.append([rb, other_idxs
if avoidcollision_rb
else set()])
2095 if avoidcollision_rb:
2097 npairs = len(gcpf.get_close_pairs(mdl,
2106 print(
"shuffle_configuration: rigid body placed close to other %d particles, trying again..." % npairs)
2107 print(
"shuffle_configuration: rigid body name: " + rb.get_name())
2108 if niter == niterations:
2109 raise ValueError(
"tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
2113 print(
'shuffling', len(flexible_beads),
'flexible beads')
2114 for fb
in flexible_beads:
2116 if avoidcollision_fb:
2118 other_idxs = all_idxs - fb_idxs
2124 while niter < niterations:
2137 xyz = memb.get_internal_coordinates()
2142 rf = memb.get_rigid_body().get_reference_frame()
2143 glob_to_int = rf.get_transformation_from()
2144 memb.set_internal_coordinates(
2145 glob_to_int.get_transformed(translation))
2147 xyz_transformed=transformation.get_transformed(xyz)
2148 memb.set_internal_coordinates(xyz_transformed)
2149 debug.append([xyz,other_idxs
if avoidcollision_fb
else set()])
2156 debug.append([d,other_idxs
if avoidcollision_fb
else set()])
2162 if avoidcollision_fb:
2164 npairs = len(gcpf.get_close_pairs(mdl,
2172 print(
"shuffle_configuration: floppy body placed close to other %d particles, trying again..." % npairs)
2173 if niter == niterations:
2174 raise ValueError(
"tried the maximum number of iterations to avoid collisions, increase the distance cutoff")
2180 class ColorHierarchy(object):
2182 def __init__(self,hier):
2183 import matplotlib
as mpl
2185 import matplotlib.pyplot
as plt
2189 hier.ColorHierarchy=self
2193 self.method=self.nochange
2201 def get_color(self,fl):
2204 def get_log_scale(self,fl):
2207 return math.log(fl+eps)
2209 def color_by_resid(self):
2210 self.method=self.color_by_resid
2211 self.scheme=self.mpl.cm.rainbow
2212 for mol
in self.mols:
2218 c=self.get_color(float(ri)/self.last)
2222 avr=sum(ris)/len(ris)
2223 c=self.get_color(float(avr)/self.last)
2226 def color_by_uncertainty(self):
2227 self.method=self.color_by_uncertainty
2228 self.scheme=self.mpl.cm.jet
2235 self.first=self.get_log_scale(1.0)
2236 self.last=self.get_log_scale(100.0)
2238 value=self.get_log_scale(unc_dict[p])
2239 if value>=self.last: value=self.last
2240 if value<=self.first: value=self.first
2241 c=self.get_color((value-self.first)/(self.last-self.first))
2244 def get_color_bar(self,filename):
2245 import matplotlib
as mpl
2247 import matplotlib.pyplot
as plt
2250 fig = plt.figure(figsize=(8, 3))
2251 ax1 = fig.add_axes([0.05, 0.80, 0.9, 0.15])
2254 norm = mpl.colors.Normalize(vmin=0.0, vmax=1.0)
2256 if self.method == self.color_by_uncertainty:
2257 angticks=[1.0,2.5,5.0,10.0,25.0,50.0,100.0]
2261 vvalue=(self.get_log_scale(at)-self.first)/(self.last-self.first)
2262 if vvalue <= 1.0
and vvalue >= 0.0:
2263 vvalues.append(vvalue)
2264 marks.append(str(at))
2265 cb1 = mpl.colorbar.ColorbarBase(ax1, cmap=cmap,
2268 orientation=
'horizontal')
2269 print(self.first,self.last,marks,vvalues)
2270 cb1.ax.set_xticklabels(marks)
2271 cb1.set_label(
'Angstorm')
2272 plt.savefig(filename, dpi=150, transparent=
True)
2279 """Given a chimera color name, return RGB"""
2280 d = {
'aquamarine': (0.4980392156862745, 1.0, 0.8313725490196079),
2281 'black': (0.0, 0.0, 0.0),
2282 'blue': (0.0, 0.0, 1.0),
2283 'brown': (0.6470588235294118, 0.16470588235294117, 0.16470588235294117),
2284 'chartreuse': (0.4980392156862745, 1.0, 0.0),
2285 'coral': (1.0, 0.4980392156862745, 0.3137254901960784),
2286 'cornflower blue': (0.39215686274509803, 0.5843137254901961, 0.9294117647058824),
2287 'cyan': (0.0, 1.0, 1.0),
2288 'dark cyan': (0.0, 0.5450980392156862, 0.5450980392156862),
2289 'dark gray': (0.6627450980392157, 0.6627450980392157, 0.6627450980392157),
2290 'dark green': (0.0, 0.39215686274509803, 0.0),
2291 'dark khaki': (0.7411764705882353, 0.7176470588235294, 0.4196078431372549),
2292 'dark magenta': (0.5450980392156862, 0.0, 0.5450980392156862),
2293 'dark olive green': (0.3333333333333333, 0.4196078431372549, 0.1843137254901961),
2294 'dark red': (0.5450980392156862, 0.0, 0.0),
2295 'dark slate blue': (0.2823529411764706, 0.23921568627450981, 0.5450980392156862),
2296 'dark slate gray': (0.1843137254901961, 0.30980392156862746, 0.30980392156862746),
2297 'deep pink': (1.0, 0.0784313725490196, 0.5764705882352941),
2298 'deep sky blue': (0.0, 0.7490196078431373, 1.0),
2299 'dim gray': (0.4117647058823529, 0.4117647058823529, 0.4117647058823529),
2300 'dodger blue': (0.11764705882352941, 0.5647058823529412, 1.0),
2301 'firebrick': (0.6980392156862745, 0.13333333333333333, 0.13333333333333333),
2302 'forest green': (0.13333333333333333, 0.5450980392156862, 0.13333333333333333),
2303 'gold': (1.0, 0.8431372549019608, 0.0),
2304 'goldenrod': (0.8549019607843137, 0.6470588235294118, 0.12549019607843137),
2305 'gray': (0.7450980392156863, 0.7450980392156863, 0.7450980392156863),
2306 'green': (0.0, 1.0, 0.0),
2307 'hot pink': (1.0, 0.4117647058823529, 0.7058823529411765),
2308 'khaki': (0.9411764705882353, 0.9019607843137255, 0.5490196078431373),
2309 'light blue': (0.6784313725490196, 0.8470588235294118, 0.9019607843137255),
2310 'light gray': (0.8274509803921568, 0.8274509803921568, 0.8274509803921568),
2311 'light green': (0.5647058823529412, 0.9333333333333333, 0.5647058823529412),
2312 'light sea green': (0.12549019607843137, 0.6980392156862745, 0.6666666666666666),
2313 'lime green': (0.19607843137254902, 0.803921568627451, 0.19607843137254902),
2314 'magenta': (1.0, 0.0, 1.0),
2315 'medium blue': (0.19607843137254902, 0.19607843137254902, 0.803921568627451),
2316 'medium purple': (0.5764705882352941, 0.4392156862745098, 0.8588235294117647),
2317 'navy blue': (0.0, 0.0, 0.5019607843137255),
2318 'olive drab': (0.4196078431372549, 0.5568627450980392, 0.13725490196078433),
2319 'orange red': (1.0, 0.27058823529411763, 0.0),
2320 'orange': (1.0, 0.4980392156862745, 0.0),
2321 'orchid': (0.8549019607843137, 0.4392156862745098, 0.8392156862745098),
2322 'pink': (1.0, 0.7529411764705882, 0.796078431372549),
2323 'plum': (0.8666666666666667, 0.6274509803921569, 0.8666666666666667),
2324 'purple': (0.6274509803921569, 0.12549019607843137, 0.9411764705882353),
2325 'red': (1.0, 0.0, 0.0),
2326 'rosy brown': (0.7372549019607844, 0.5607843137254902, 0.5607843137254902),
2327 'salmon': (0.9803921568627451, 0.5019607843137255, 0.4470588235294118),
2328 'sandy brown': (0.9568627450980393, 0.6431372549019608, 0.3764705882352941),
2329 'sea green': (0.1803921568627451, 0.5450980392156862, 0.3411764705882353),
2330 'sienna': (0.6274509803921569, 0.3215686274509804, 0.17647058823529413),
2331 'sky blue': (0.5294117647058824, 0.807843137254902, 0.9215686274509803),
2332 'slate gray': (0.4392156862745098, 0.5019607843137255, 0.5647058823529412),
2333 'spring green': (0.0, 1.0, 0.4980392156862745),
2334 'steel blue': (0.27450980392156865, 0.5098039215686274, 0.7058823529411765),
2335 'tan': (0.8235294117647058, 0.7058823529411765, 0.5490196078431373),
2336 'turquoise': (0.25098039215686274, 0.8784313725490196, 0.8156862745098039),
2337 'violet red': (0.8156862745098039, 0.12549019607843137, 0.5647058823529412),
2338 'white': (1.0, 1.0, 1.0),
2339 'yellow': (1.0, 1.0, 0.0)}
2345 "reds":[(
"maroon",
"#800000",(128,0,0)),(
"dark red",
"#8B0000",(139,0,0)),
2346 (
"brown",
"#A52A2A",(165,42,42)),(
"firebrick",
"#B22222",(178,34,34)),
2347 (
"crimson",
"#DC143C",(220,20,60)),(
"red",
"#FF0000",(255,0,0)),
2348 (
"tomato",
"#FF6347",(255,99,71)),(
"coral",
"#FF7F50",(255,127,80)),
2349 (
"indian red",
"#CD5C5C",(205,92,92)),(
"light coral",
"#F08080",(240,128,128)),
2350 (
"dark salmon",
"#E9967A",(233,150,122)),(
"salmon",
"#FA8072",(250,128,114)),
2351 (
"light salmon",
"#FFA07A",(255,160,122)),(
"orange red",
"#FF4500",(255,69,0)),
2352 (
"dark orange",
"#FF8C00",(255,140,0))],
2353 "yellows":[(
"orange",
"#FFA500",(255,165,0)),(
"gold",
"#FFD700",(255,215,0)),
2354 (
"dark golden rod",
"#B8860B",(184,134,11)),(
"golden rod",
"#DAA520",(218,165,32)),
2355 (
"pale golden rod",
"#EEE8AA",(238,232,170)),(
"dark khaki",
"#BDB76B",(189,183,107)),
2356 (
"khaki",
"#F0E68C",(240,230,140)),(
"olive",
"#808000",(128,128,0)),
2357 (
"yellow",
"#FFFF00",(255,255,0)),(
"antique white",
"#FAEBD7",(250,235,215)),
2358 (
"beige",
"#F5F5DC",(245,245,220)),(
"bisque",
"#FFE4C4",(255,228,196)),
2359 (
"blanched almond",
"#FFEBCD",(255,235,205)),(
"wheat",
"#F5DEB3",(245,222,179)),
2360 (
"corn silk",
"#FFF8DC",(255,248,220)),(
"lemon chiffon",
"#FFFACD",(255,250,205)),
2361 (
"light golden rod yellow",
"#FAFAD2",(250,250,210)),(
"light yellow",
"#FFFFE0",(255,255,224))],
2362 "greens":[(
"yellow green",
"#9ACD32",(154,205,50)),(
"dark olive green",
"#556B2F",(85,107,47)),
2363 (
"olive drab",
"#6B8E23",(107,142,35)),(
"lawn green",
"#7CFC00",(124,252,0)),
2364 (
"chart reuse",
"#7FFF00",(127,255,0)),(
"green yellow",
"#ADFF2F",(173,255,47)),
2365 (
"dark green",
"#006400",(0,100,0)),(
"green",
"#008000",(0,128,0)),
2366 (
"forest green",
"#228B22",(34,139,34)),(
"lime",
"#00FF00",(0,255,0)),
2367 (
"lime green",
"#32CD32",(50,205,50)),(
"light green",
"#90EE90",(144,238,144)),
2368 (
"pale green",
"#98FB98",(152,251,152)),(
"dark sea green",
"#8FBC8F",(143,188,143)),
2369 (
"medium spring green",
"#00FA9A",(0,250,154)),(
"spring green",
"#00FF7F",(0,255,127)),
2370 (
"sea green",
"#2E8B57",(46,139,87)),(
"medium aqua marine",
"#66CDAA",(102,205,170)),
2371 (
"medium sea green",
"#3CB371",(60,179,113)),(
"light sea green",
"#20B2AA",(32,178,170)),
2372 (
"dark slate gray",
"#2F4F4F",(47,79,79)),(
"teal",
"#008080",(0,128,128)),
2373 (
"dark cyan",
"#008B8B",(0,139,139))],
2374 "blues":[(
"dark turquoise",
"#00CED1",(0,206,209)),
2375 (
"turquoise",
"#40E0D0",(64,224,208)),(
"medium turquoise",
"#48D1CC",(72,209,204)),
2376 (
"pale turquoise",
"#AFEEEE",(175,238,238)),(
"aqua marine",
"#7FFFD4",(127,255,212)),
2377 (
"powder blue",
"#B0E0E6",(176,224,230)),(
"cadet blue",
"#5F9EA0",(95,158,160)),
2378 (
"steel blue",
"#4682B4",(70,130,180)),(
"corn flower blue",
"#6495ED",(100,149,237)),
2379 (
"deep sky blue",
"#00BFFF",(0,191,255)),(
"dodger blue",
"#1E90FF",(30,144,255)),
2380 (
"light blue",
"#ADD8E6",(173,216,230)),(
"sky blue",
"#87CEEB",(135,206,235)),
2381 (
"light sky blue",
"#87CEFA",(135,206,250)),(
"midnight blue",
"#191970",(25,25,112)),
2382 (
"navy",
"#000080",(0,0,128)),(
"dark blue",
"#00008B",(0,0,139)),
2383 (
"medium blue",
"#0000CD",(0,0,205)),(
"blue",
"#0000FF",(0,0,255)),(
"royal blue",
"#4169E1",(65,105,225)),
2384 (
"aqua",
"#00FFFF",(0,255,255)),(
"cyan",
"#00FFFF",(0,255,255)),(
"light cyan",
"#E0FFFF",(224,255,255))],
2385 "violets":[(
"blue violet",
"#8A2BE2",(138,43,226)),(
"indigo",
"#4B0082",(75,0,130)),
2386 (
"dark slate blue",
"#483D8B",(72,61,139)),(
"slate blue",
"#6A5ACD",(106,90,205)),
2387 (
"medium slate blue",
"#7B68EE",(123,104,238)),(
"medium purple",
"#9370DB",(147,112,219)),
2388 (
"dark magenta",
"#8B008B",(139,0,139)),(
"dark violet",
"#9400D3",(148,0,211)),
2389 (
"dark orchid",
"#9932CC",(153,50,204)),(
"medium orchid",
"#BA55D3",(186,85,211)),
2390 (
"purple",
"#800080",(128,0,128)),(
"thistle",
"#D8BFD8",(216,191,216)),
2391 (
"plum",
"#DDA0DD",(221,160,221)),(
"violet",
"#EE82EE",(238,130,238)),
2392 (
"magenta / fuchsia",
"#FF00FF",(255,0,255)),(
"orchid",
"#DA70D6",(218,112,214)),
2393 (
"medium violet red",
"#C71585",(199,21,133)),(
"pale violet red",
"#DB7093",(219,112,147)),
2394 (
"deep pink",
"#FF1493",(255,20,147)),(
"hot pink",
"#FF69B4",(255,105,180)),
2395 (
"light pink",
"#FFB6C1",(255,182,193)),(
"pink",
"#FFC0CB",(255,192,203))],
2396 "browns":[(
"saddle brown",
"#8B4513",(139,69,19)),(
"sienna",
"#A0522D",(160,82,45)),
2397 (
"chocolate",
"#D2691E",(210,105,30)),(
"peru",
"#CD853F",(205,133,63)),
2398 (
"sandy brown",
"#F4A460",(244,164,96)),(
"burly wood",
"#DEB887",(222,184,135)),
2399 (
"tan",
"#D2B48C",(210,180,140)),(
"rosy brown",
"#BC8F8F",(188,143,143)),
2400 (
"moccasin",
"#FFE4B5",(255,228,181)),(
"navajo white",
"#FFDEAD",(255,222,173)),
2401 (
"peach puff",
"#FFDAB9",(255,218,185)),(
"misty rose",
"#FFE4E1",(255,228,225)),
2402 (
"lavender blush",
"#FFF0F5",(255,240,245)),(
"linen",
"#FAF0E6",(250,240,230)),
2403 (
"old lace",
"#FDF5E6",(253,245,230)),(
"papaya whip",
"#FFEFD5",(255,239,213)),
2404 (
"sea shell",
"#FFF5EE",(255,245,238))],
2405 "greys":[(
"black",
"#000000",(0,0,0)),(
"dim gray / dim grey",
"#696969",(105,105,105)),
2406 (
"gray / grey",
"#808080",(128,128,128)),(
"dark gray / dark grey",
"#A9A9A9",(169,169,169)),
2407 (
"silver",
"#C0C0C0",(192,192,192)),(
"light gray / light grey",
"#D3D3D3",(211,211,211)),
2408 (
"gainsboro",
"#DCDCDC",(220,220,220)),(
"white smoke",
"#F5F5F5",(245,245,245)),
2409 (
"white",
"#FFFFFF",(255,255,255))]}
2411 def assign_color_group(self,color_group,representation,component_names):
2412 for n,p
in enumerate(component_names):
2414 psel=s.get_selected_particles()
2415 ctuple=self.colors[color_group][n]
2416 print(
"Assigning "+p+
" to color "+ctuple[0])
2425 def get_list_distant_colors(self):
2426 cnames = [
'#F0F8FF',
'#FAEBD7',
'#00FFFF',
'#7FFFD4',
'#F0FFFF',
'#F5F5DC',
2427 '#FFE4C4',
'#000000',
'#FFEBCD',
'#0000FF',
'#8A2BE2',
'#A52A2A',
'#DEB887',
2428 '#5F9EA0',
'#7FFF00',
'#D2691E',
'#FF7F50',
'#6495ED',
'#FFF8DC',
'#DC143C',
2429 '#00FFFF',
'#00008B',
'#008B8B',
'#B8860B',
'#A9A9A9',
'#006400',
'#BDB76B',
2430 '#8B008B',
'#556B2F',
'#FF8C00',
'#9932CC',
'#8B0000',
'#E9967A',
'#8FBC8F',
2431 '#483D8B',
'#2F4F4F',
'#00CED1',
'#9400D3',
'#FF1493',
'#00BFFF',
'#696969',
2432 '#1E90FF',
'#B22222',
'#FFFAF0',
'#228B22',
'#FF00FF',
'#DCDCDC',
'#F8F8FF',
2433 '#FFD700',
'#DAA520',
'#808080',
'#008000',
'#ADFF2F',
'#F0FFF0',
'#FF69B4',
2434 '#CD5C5C',
'#4B0082',
'#FFFFF0',
'#F0E68C',
'#E6E6FA',
'#FFF0F5',
'#7CFC00',
2435 '#FFFACD',
'#ADD8E6',
'#F08080',
'#E0FFFF',
'#FAFAD2',
'#90EE90',
'#D3D3D3',
2436 '#FFB6C1',
'#FFA07A',
'#20B2AA',
'#87CEFA',
'#778899',
'#B0C4DE',
'#FFFFE0',
2437 '#00FF00',
'#32CD32',
'#FAF0E6',
'#FF00FF',
'#800000',
'#66CDAA',
'#0000CD',
2438 '#BA55D3',
'#9370DB',
'#3CB371',
'#7B68EE',
'#00FA9A',
'#48D1CC',
'#C71585',
2439 '#191970',
'#F5FFFA',
'#FFE4E1',
'#FFE4B5',
'#FFDEAD',
'#000080',
'#FDF5E6',
2440 '#808000',
'#6B8E23',
'#FFA500',
'#FF4500',
'#DA70D6',
'#EEE8AA',
'#98FB98',
2441 '#AFEEEE',
'#DB7093',
'#FFEFD5',
'#FFDAB9',
'#CD853F',
'#FFC0CB',
'#DDA0DD',
2442 '#B0E0E6',
'#800080',
'#FF0000',
'#BC8F8F',
'#4169E1',
'#8B4513',
'#FA8072',
2443 '#FAA460',
'#2E8B57',
'#FFF5EE',
'#A0522D',
'#C0C0C0',
'#87CEEB',
'#6A5ACD',
2444 '#708090',
'#FFFAFA',
'#00FF7F',
'#4682B4',
'#D2B48C',
'#008080',
'#D8BFD8',
2445 '#FF6347',
'#40E0D0',
'#EE82EE',
'#F5DEB3',
'#FFFFFF',
'#F5F5F5',
'#FFFF00',
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.
def add_script_provenance
Tag the given particle with the current Python script.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
void add_particles(RMF::FileHandle fh, const ParticlesTemp &hs)
Add resolution to a particle.
static Weight setup_particle(Model *m, ParticleIndex pi)
Set up an empty Weight.
A decorator for a particle which has bonds.
double get_drmsd(const Vector3DsOrXYZs0 &m0, const Vector3DsOrXYZs1 &m1)
Calculate distance the root mean square deviation between two sets of 3D points.
Rotation3D get_random_rotation_3d(const Rotation3D ¢er, double distance)
Pick a rotation at random near the provided one.
Legacy PMI1 module to represent, score, sample and analyze models.
Extends the functionality of IMP.atom.Molecule.
IMP::Vector< Color > Colors
Set of python classes to create a multi-state, multi-resolution IMP hierarchy.
ParticlesTemp get_particles(Model *m, const ParticleIndexes &ps)
Get the particles from a list of indexes.
static Surface setup_particle(Model *m, ParticleIndex pi)
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.
double get_distance(XYZR a, XYZR b)
Compute the sphere distance between a and b.
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.
Add uncertainty to a particle.
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.
static bool get_is_setup(Model *m, ParticleIndex pi)
ParticleIndexPairs get_indexes(const ParticlePairsTemp &ps)
Get the indexes from a list of particle pairs.
std::string get_module_version()
Return the version of this module, as a string.
static bool get_is_setup(Model *m, ParticleIndex pi)
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)
static Colored setup_particle(Model *m, ParticleIndex pi, Color color)
A decorator for a particle that is part of a rigid body but not rigid.
std::ostream & show(Hierarchy h, std::ostream &out=std::cout)
Print the hierarchy using a given decorator to display each node.
Find all nearby pairs by testing all pairs.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
def deprecated_object
Python decorator to mark a class as deprecated.
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.
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.
double get_resolution(Model *m, ParticleIndex pi)
Estimate the resolution of the hierarchy as used by Representation.
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)
Transformation3D get_random_local_transformation(Vector3D origin, double max_translation=5., double max_angle_in_rad=0.26)
Get a local transformation.
Inferential scoring building on methods developed as part of the Inferential Structure Determination ...