1 """@namespace IMP.pmi.output
2 Classes for writing output files and processing them.
5 from __future__
import print_function, division
24 import cPickle
as pickle
29 class _ChainIDs(object):
30 """Map indices to multi-character chain IDs.
31 We label the first 26 chains A-Z, then we move to two-letter
32 chain IDs: AA through AZ, then BA through BZ, through to ZZ.
33 This continues with longer chain IDs."""
34 def __getitem__(self, ind):
35 chars = string.ascii_uppercase
39 ids.append(chars[ind % lc])
41 ids.append(chars[ind])
42 return "".join(reversed(ids))
46 """Base class for capturing a modeling protocol.
47 Unlike simple output of model coordinates, a complete
48 protocol includes the input data used, details on the restraints,
49 sampling, and clustering, as well as output models.
50 Use via IMP.pmi.topology.System.add_protocol_output().
52 @see IMP.pmi.mmcif.ProtocolOutput for a concrete subclass that outputs
60 if isinstance(elt, (tuple, list)):
61 for elt2
in _flatten(elt):
67 def _disambiguate_chain(chid, seen_chains):
68 """Make sure that the chain ID is unique; warn and correct if it isn't"""
73 if chid
in seen_chains:
74 warnings.warn(
"Duplicate chain ID '%s' encountered" % chid,
77 for suffix
in itertools.count(1):
78 new_chid = chid +
"%d" % suffix
79 if new_chid
not in seen_chains:
80 seen_chains.add(new_chid)
86 def _write_pdb_internal(flpdb, particle_infos_for_pdb, geometric_center,
87 write_all_residues_per_bead):
88 for n, tupl
in enumerate(particle_infos_for_pdb):
89 (xyz, atom_type, residue_type,
90 chain_id, residue_index, all_indexes, radius) = tupl
92 atom_type = IMP.atom.AT_CA
93 if write_all_residues_per_bead
and all_indexes
is not None:
94 for residue_number
in all_indexes:
96 IMP.atom.get_pdb_string((xyz[0] - geometric_center[0],
97 xyz[1] - geometric_center[1],
98 xyz[2] - geometric_center[2]),
99 n+1, atom_type, residue_type,
100 chain_id[:1], residue_number,
' ',
104 IMP.atom.get_pdb_string((xyz[0] - geometric_center[0],
105 xyz[1] - geometric_center[1],
106 xyz[2] - geometric_center[2]),
107 n+1, atom_type, residue_type,
108 chain_id[:1], residue_index,
' ',
110 flpdb.write(
"ENDMDL\n")
113 _Entity = collections.namedtuple(
'_Entity', (
'id',
'seq'))
114 _ChainInfo = collections.namedtuple(
'_ChainInfo', (
'entity',
'name'))
117 def _get_chain_info(chains, root_hier):
121 for mol
in IMP.atom.get_by_type(root_hier, IMP.atom.MOLECULE_TYPE):
123 chain_id = chains[molname]
125 seq = chain.get_sequence()
126 if seq
not in entities:
127 entities[seq] = e = _Entity(id=len(entities)+1, seq=seq)
128 all_entities.append(e)
129 entity = entities[seq]
130 info = _ChainInfo(entity=entity, name=molname)
131 chain_info[chain_id] = info
132 return chain_info, all_entities
135 def _write_mmcif_internal(flpdb, particle_infos_for_pdb, geometric_center,
136 write_all_residues_per_bead, chains, root_hier):
138 chain_info, entities = _get_chain_info(chains, root_hier)
140 writer = ihm.format.CifWriter(flpdb)
141 writer.start_block(
'model')
142 with writer.category(
"_entry")
as lp:
145 with writer.loop(
"_entity", [
"id",
"type"])
as lp:
147 lp.write(id=e.id, type=
"polymer")
149 with writer.loop(
"_entity_poly",
150 [
"entity_id",
"pdbx_seq_one_letter_code"])
as lp:
152 lp.write(entity_id=e.id, pdbx_seq_one_letter_code=e.seq)
154 with writer.loop(
"_struct_asym", [
"id",
"entity_id",
"details"])
as lp:
156 for chid
in sorted(chains.values(), key=
lambda x: (len(x.strip()), x)):
157 ci = chain_info[chid]
158 lp.write(id=chid, entity_id=ci.entity.id, details=ci.name)
160 with writer.loop(
"_atom_site",
161 [
"group_PDB",
"type_symbol",
"label_atom_id",
162 "label_comp_id",
"label_asym_id",
"label_seq_id",
164 "Cartn_x",
"Cartn_y",
"Cartn_z",
"label_entity_id",
165 "pdbx_pdb_model_num",
168 for n, tupl
in enumerate(particle_infos_for_pdb):
169 (xyz, atom_type, residue_type,
170 chain_id, residue_index, all_indexes, radius) = tupl
171 ci = chain_info[chain_id]
172 if atom_type
is None:
173 atom_type = IMP.atom.AT_CA
174 c = (xyz[0] - geometric_center[0],
175 xyz[1] - geometric_center[1],
176 xyz[2] - geometric_center[2])
177 if write_all_residues_per_bead
and all_indexes
is not None:
178 for residue_number
in all_indexes:
179 lp.write(group_PDB=
'ATOM',
181 label_atom_id=atom_type.get_string(),
182 label_comp_id=residue_type.get_string(),
183 label_asym_id=chain_id,
184 label_seq_id=residue_index,
185 auth_seq_id=residue_index, Cartn_x=c[0],
186 Cartn_y=c[1], Cartn_z=c[2], id=ordinal,
187 pdbx_pdb_model_num=1,
188 label_entity_id=ci.entity.id)
191 lp.write(group_PDB=
'ATOM', type_symbol=
'C',
192 label_atom_id=atom_type.get_string(),
193 label_comp_id=residue_type.get_string(),
194 label_asym_id=chain_id,
195 label_seq_id=residue_index,
196 auth_seq_id=residue_index, Cartn_x=c[0],
197 Cartn_y=c[1], Cartn_z=c[2], id=ordinal,
198 pdbx_pdb_model_num=1,
199 label_entity_id=ci.entity.id)
204 """Class for easy writing of PDBs, RMFs, and stat files
206 @note Model should be updated prior to writing outputs.
208 def __init__(self, ascii=True, atomistic=False):
209 self.dictionary_pdbs = {}
211 self.dictionary_rmfs = {}
212 self.dictionary_stats = {}
213 self.dictionary_stats2 = {}
214 self.best_score_list =
None
215 self.nbestscoring =
None
217 self.replica_exchange =
False
222 self.chainids =
"ABCDEFGHIJKLMNOPQRSTUVWXYZ" \
223 "abcdefghijklmnopqrstuvwxyz0123456789"
225 self.multi_chainids = _ChainIDs()
227 self.particle_infos_for_pdb = {}
228 self.atomistic = atomistic
229 self.use_pmi2 =
False
232 """Get a list of all PDB files being output by this instance"""
233 return list(self.dictionary_pdbs.keys())
235 def get_rmf_names(self):
236 return list(self.dictionary_rmfs.keys())
238 def get_stat_names(self):
239 return list(self.dictionary_stats.keys())
241 def _init_dictchain(self, name, prot, multichar_chain=False, mmcif=False):
242 self.dictchain[name] = {}
243 self.use_pmi2 =
False
249 self.atomistic =
True
250 for n, mol
in enumerate(IMP.atom.get_by_type(
251 prot, IMP.atom.MOLECULE_TYPE)):
253 if not mmcif
and len(chid) > 1:
255 "The system contains at least one chain ID (%s) that "
256 "is more than 1 character long; this cannot be "
257 "represented in PDB. Either write mmCIF files "
258 "instead, or assign 1-character IDs to all chains "
259 "(this can be done with the `chain_ids` argument to "
260 "BuildSystem.add_state())." % chid)
261 chid = _disambiguate_chain(chid, seen_chains)
263 self.dictchain[name][molname] = chid
265 chainids = self.multi_chainids
if multichar_chain \
267 for n, i
in enumerate(self.dictionary_pdbs[name].get_children()):
268 self.dictchain[name][i.get_name()] = chainids[n]
272 @param name The PDB filename
273 @param prot The hierarchy to write to this pdb file
274 @param mmcif If True, write PDBs in mmCIF format
275 @note if the PDB name is 'System' then will use Selection
278 flpdb = open(name,
'w')
280 self.dictionary_pdbs[name] = prot
281 self._pdb_mmcif[name] = mmcif
282 self._init_dictchain(name, prot, mmcif=mmcif)
284 def write_psf(self, filename, name):
285 flpsf = open(filename,
'w')
286 flpsf.write(
"PSF CMAP CHEQ" +
"\n")
287 index_residue_pair_list = {}
288 (particle_infos_for_pdb, geometric_center) = \
289 self.get_particle_infos_for_pdb_writing(name)
290 nparticles = len(particle_infos_for_pdb)
291 flpsf.write(str(nparticles) +
" !NATOM" +
"\n")
292 for n, p
in enumerate(particle_infos_for_pdb):
297 flpsf.write(
'{0:8d}{1:1s}{2:4s}{3:1s}{4:4s}{5:1s}{6:4s}{7:1s}'
298 '{8:4s}{9:1s}{10:4s}{11:14.6f}{12:14.6f}{13:8d}'
299 '{14:14.6f}{15:14.6f}'.format(
300 atom_index,
" ", chain,
" ", str(resid),
" ",
301 '"'+residue_type.get_string()+
'"',
" ",
"C",
302 " ",
"C", 1.0, 0.0, 0, 0.0, 0.0))
304 if chain
not in index_residue_pair_list:
305 index_residue_pair_list[chain] = [(atom_index, resid)]
307 index_residue_pair_list[chain].append((atom_index, resid))
311 for chain
in sorted(index_residue_pair_list.keys()):
313 ls = index_residue_pair_list[chain]
315 ls = sorted(ls, key=
lambda tup: tup[1])
317 indexes = [x[0]
for x
in ls]
320 indexes, lmin=2, lmax=2))
321 nbonds = len(indexes_pairs)
322 flpsf.write(str(nbonds)+
" !NBOND: bonds"+
"\n")
325 for i
in range(0, len(indexes_pairs), 4):
326 for bond
in indexes_pairs[i:i+4]:
327 flpsf.write(
'{0:8d}{1:8d}'.format(*bond))
330 del particle_infos_for_pdb
333 def write_pdb(self, name, appendmode=True,
334 translate_to_geometric_center=
False,
335 write_all_residues_per_bead=
False):
337 (particle_infos_for_pdb,
338 geometric_center) = self.get_particle_infos_for_pdb_writing(name)
340 if not translate_to_geometric_center:
341 geometric_center = (0, 0, 0)
343 filemode =
'a' if appendmode
else 'w'
344 with open(name, filemode)
as flpdb:
345 if self._pdb_mmcif[name]:
346 _write_mmcif_internal(flpdb, particle_infos_for_pdb,
348 write_all_residues_per_bead,
349 self.dictchain[name],
350 self.dictionary_pdbs[name])
352 _write_pdb_internal(flpdb, particle_infos_for_pdb,
354 write_all_residues_per_bead)
357 """Get the protein name from the particle.
358 This is done by traversing the hierarchy."""
363 p, self.dictchain[name])
365 def get_particle_infos_for_pdb_writing(self, name):
375 particle_infos_for_pdb = []
377 geometric_center = [0, 0, 0]
383 ps = sel.get_selected_particles()
387 for n, p
in enumerate(ps):
390 if protname
not in resindexes_dict:
391 resindexes_dict[protname] = []
395 rt = residue.get_residue_type()
396 resind = residue.get_index()
400 geometric_center[0] += xyz[0]
401 geometric_center[1] += xyz[1]
402 geometric_center[2] += xyz[2]
404 particle_infos_for_pdb.append(
405 (xyz, atomtype, rt, self.dictchain[name][protname],
406 resind,
None, radius))
407 resindexes_dict[protname].append(resind)
412 resind = residue.get_index()
415 if resind
in resindexes_dict[protname]:
418 resindexes_dict[protname].append(resind)
419 rt = residue.get_residue_type()
422 geometric_center[0] += xyz[0]
423 geometric_center[1] += xyz[1]
424 geometric_center[2] += xyz[2]
426 particle_infos_for_pdb.append(
427 (xyz,
None, rt, self.dictchain[name][protname], resind,
432 resind = resindexes[len(resindexes) // 2]
433 if resind
in resindexes_dict[protname]:
436 resindexes_dict[protname].append(resind)
440 geometric_center[0] += xyz[0]
441 geometric_center[1] += xyz[1]
442 geometric_center[2] += xyz[2]
444 particle_infos_for_pdb.append(
445 (xyz,
None, rt, self.dictchain[name][protname], resind,
452 if len(resindexes) > 0:
453 resind = resindexes[len(resindexes) // 2]
456 geometric_center[0] += xyz[0]
457 geometric_center[1] += xyz[1]
458 geometric_center[2] += xyz[2]
460 particle_infos_for_pdb.append(
461 (xyz,
None, rt, self.dictchain[name][protname],
462 resind, resindexes, radius))
465 geometric_center = (geometric_center[0] / atom_count,
466 geometric_center[1] / atom_count,
467 geometric_center[2] / atom_count)
471 particle_infos_for_pdb = sorted(particle_infos_for_pdb,
472 key=
lambda x: (len(x[3]), x[3], x[4]))
474 return (particle_infos_for_pdb, geometric_center)
476 def write_pdbs(self, appendmode=True, mmcif=False):
477 for pdb
in self.dictionary_pdbs.keys():
478 self.write_pdb(pdb, appendmode)
481 replica_exchange=
False, mmcif=
False,
482 best_score_file=
'best.scores.rex.py'):
483 """Prepare for writing best-scoring PDBs (or mmCIFs) for a
486 @param prefix Initial part of each PDB filename (e.g. 'model').
487 @param prot The top-level Hierarchy to output.
488 @param nbestscoring The number of best-scoring files to output.
489 @param replica_exchange Whether to combine best scores from a
490 replica exchange run.
491 @param mmcif If True, output models in mmCIF format. If False
492 (the default) output in legacy PDB format.
493 @param best_score_file The filename to use for replica
497 self._pdb_best_scoring_mmcif = mmcif
498 fileext =
'.cif' if mmcif
else '.pdb'
499 self.prefixes.append(prefix)
500 self.replica_exchange = replica_exchange
501 if not self.replica_exchange:
505 self.best_score_list = []
509 self.best_score_file_name = best_score_file
510 self.best_score_list = []
511 with open(self.best_score_file_name,
"w")
as best_score_file:
512 best_score_file.write(
513 "self.best_score_list=" + str(self.best_score_list) +
"\n")
515 self.nbestscoring = nbestscoring
516 for i
in range(self.nbestscoring):
517 name = prefix +
"." + str(i) + fileext
518 flpdb = open(name,
'w')
520 self.dictionary_pdbs[name] = prot
521 self._pdb_mmcif[name] = mmcif
522 self._init_dictchain(name, prot, mmcif=mmcif)
524 def write_pdb_best_scoring(self, score):
525 if self.nbestscoring
is None:
526 print(
"Output.write_pdb_best_scoring: init_pdb_best_scoring "
529 mmcif = self._pdb_best_scoring_mmcif
530 fileext =
'.cif' if mmcif
else '.pdb'
532 if self.replica_exchange:
534 with open(self.best_score_file_name)
as fh:
535 self.best_score_list = ast.literal_eval(
536 fh.read().split(
'=')[1])
538 if len(self.best_score_list) < self.nbestscoring:
539 self.best_score_list.append(score)
540 self.best_score_list.sort()
541 index = self.best_score_list.index(score)
542 for prefix
in self.prefixes:
543 for i
in range(len(self.best_score_list) - 2, index - 1, -1):
544 oldname = prefix +
"." + str(i) + fileext
545 newname = prefix +
"." + str(i + 1) + fileext
547 if os.path.exists(newname):
549 os.rename(oldname, newname)
550 filetoadd = prefix +
"." + str(index) + fileext
551 self.write_pdb(filetoadd, appendmode=
False)
554 if score < self.best_score_list[-1]:
555 self.best_score_list.append(score)
556 self.best_score_list.sort()
557 self.best_score_list.pop(-1)
558 index = self.best_score_list.index(score)
559 for prefix
in self.prefixes:
560 for i
in range(len(self.best_score_list) - 1,
562 oldname = prefix +
"." + str(i) + fileext
563 newname = prefix +
"." + str(i + 1) + fileext
564 os.rename(oldname, newname)
565 filenametoremove = prefix + \
566 "." + str(self.nbestscoring) + fileext
567 os.remove(filenametoremove)
568 filetoadd = prefix +
"." + str(index) + fileext
569 self.write_pdb(filetoadd, appendmode=
False)
571 if self.replica_exchange:
573 with open(self.best_score_file_name,
"w")
as best_score_file:
574 best_score_file.write(
575 "self.best_score_list=" + str(self.best_score_list) +
'\n')
577 def init_rmf(self, name, hierarchies, rs=None, geometries=None,
580 Initialize an RMF file
582 @param name the name of the RMF file
583 @param hierarchies the hierarchies to be included (it is a list)
584 @param rs optional, the restraint sets (it is a list)
585 @param geometries optional, the geometries (it is a list)
586 @param listofobjects optional, the list of objects for the stat
589 rh = RMF.create_rmf_file(name)
592 outputkey_rmfkey =
None
596 if geometries
is not None:
598 if listofobjects
is not None:
599 cat = rh.get_category(
"stat")
600 outputkey_rmfkey = {}
601 for o
in listofobjects:
602 if "get_output" not in dir(o):
604 "Output: object %s doesn't have get_output() method"
606 output = o.get_output()
607 for outputkey
in output:
608 rmftag = RMF.string_tag
609 if isinstance(output[outputkey], float):
610 rmftag = RMF.float_tag
611 elif isinstance(output[outputkey], int):
613 elif isinstance(output[outputkey], str):
614 rmftag = RMF.string_tag
616 rmftag = RMF.string_tag
617 rmfkey = rh.get_key(cat, outputkey, rmftag)
618 outputkey_rmfkey[outputkey] = rmfkey
619 outputkey_rmfkey[
"rmf_file"] = \
620 rh.get_key(cat,
"rmf_file", RMF.string_tag)
621 outputkey_rmfkey[
"rmf_frame_index"] = \
622 rh.get_key(cat,
"rmf_frame_index", RMF.int_tag)
624 self.dictionary_rmfs[name] = (rh, cat, outputkey_rmfkey, listofobjects)
626 def add_restraints_to_rmf(self, name, objectlist):
627 for o
in _flatten(objectlist):
629 rs = o.get_restraint_for_rmf()
630 if not isinstance(rs, (list, tuple)):
633 rs = [o.get_restraint()]
635 self.dictionary_rmfs[name][0], rs)
637 def add_geometries_to_rmf(self, name, objectlist):
639 geos = o.get_geometries()
642 def add_particle_pair_from_restraints_to_rmf(self, name, objectlist):
645 pps = o.get_particle_pairs()
648 self.dictionary_rmfs[name][0],
651 def write_rmf(self, name):
653 if self.dictionary_rmfs[name][1]
is not None:
654 outputkey_rmfkey = self.dictionary_rmfs[name][2]
655 listofobjects = self.dictionary_rmfs[name][3]
656 for o
in listofobjects:
657 output = o.get_output()
658 for outputkey
in output:
659 rmfkey = outputkey_rmfkey[outputkey]
661 n = self.dictionary_rmfs[name][0].get_root_node()
662 n.set_value(rmfkey, output[outputkey])
663 except NotImplementedError:
665 rmfkey = outputkey_rmfkey[
"rmf_file"]
666 self.dictionary_rmfs[name][0].get_root_node().set_value(
668 rmfkey = outputkey_rmfkey[
"rmf_frame_index"]
670 self.dictionary_rmfs[name][0].get_root_node().set_value(
672 self.dictionary_rmfs[name][0].flush()
674 def close_rmf(self, name):
675 rh = self.dictionary_rmfs[name][0]
676 del self.dictionary_rmfs[name]
679 def write_rmfs(self):
680 for rmfinfo
in self.dictionary_rmfs.keys():
681 self.write_rmf(rmfinfo[0])
683 def init_stat(self, name, listofobjects):
685 flstat = open(name,
'w')
688 flstat = open(name,
'wb')
692 for o
in listofobjects:
693 if "get_output" not in dir(o):
695 "Output: object %s doesn't have get_output() method"
697 self.dictionary_stats[name] = listofobjects
699 def set_output_entry(self, key, value):
700 self.initoutput.update({key: value})
702 def write_stat(self, name, appendmode=True):
703 output = self.initoutput
704 for obj
in self.dictionary_stats[name]:
707 dfiltered = dict((k, v)
for k, v
in d.items()
if k[0] !=
"_")
708 output.update(dfiltered)
716 flstat = open(name, writeflag)
717 flstat.write(
"%s \n" % output)
720 flstat = open(name, writeflag +
'b')
721 pickle.dump(output, flstat, 2)
724 def write_stats(self):
725 for stat
in self.dictionary_stats.keys():
726 self.write_stat(stat)
728 def get_stat(self, name):
730 for obj
in self.dictionary_stats[name]:
731 output.update(obj.get_output())
734 def write_test(self, name, listofobjects):
735 flstat = open(name,
'w')
736 output = self.initoutput
737 for o
in listofobjects:
738 if "get_test_output" not in dir(o)
and "get_output" not in dir(o):
740 "Output: object %s doesn't have get_output() or "
741 "get_test_output() method" % str(o))
742 self.dictionary_stats[name] = listofobjects
744 for obj
in self.dictionary_stats[name]:
746 d = obj.get_test_output()
750 dfiltered = dict((k, v)
for k, v
in d.items()
if k[0] !=
"_")
751 output.update(dfiltered)
752 flstat.write(
"%s \n" % output)
755 def test(self, name, listofobjects, tolerance=1e-5):
756 output = self.initoutput
757 for o
in listofobjects:
758 if "get_test_output" not in dir(o)
and "get_output" not in dir(o):
760 "Output: object %s doesn't have get_output() or "
761 "get_test_output() method" % str(o))
762 for obj
in listofobjects:
764 output.update(obj.get_test_output())
766 output.update(obj.get_output())
768 flstat = open(name,
'r')
772 test_dict = ast.literal_eval(fl)
775 old_value = str(test_dict[k])
776 new_value = str(output[k])
784 fold = float(old_value)
785 fnew = float(new_value)
786 diff = abs(fold - fnew)
788 print(
"%s: test failed, old value: %s new value %s; "
789 "diff %f > %f" % (str(k), str(old_value),
790 str(new_value), diff,
791 tolerance), file=sys.stderr)
793 elif test_dict[k] != output[k]:
794 if len(old_value) < 50
and len(new_value) < 50:
795 print(
"%s: test failed, old value: %s new value %s"
796 % (str(k), old_value, new_value),
800 print(
"%s: test failed, omitting results (too long)"
801 % str(k), file=sys.stderr)
805 print(
"%s from old objects (file %s) not in new objects"
806 % (str(k), str(name)), file=sys.stderr)
810 def get_environment_variables(self):
812 return str(os.environ)
814 def get_versions_of_relevant_modules(self):
821 versions[
"ISD2_VERSION"] = IMP.isd2.get_module_version()
826 versions[
"ISD_EMXL_VERSION"] = IMP.isd_emxl.get_module_version()
831 def init_stat2(self, name, listofobjects, extralabels=None,
832 listofsummedobjects=
None):
839 if listofsummedobjects
is None:
840 listofsummedobjects = []
841 if extralabels
is None:
843 flstat = open(name,
'w')
845 stat2_keywords = {
"STAT2HEADER":
"STAT2HEADER"}
846 stat2_keywords.update(
847 {
"STAT2HEADER_ENVIRON": str(self.get_environment_variables())})
848 stat2_keywords.update(
849 {
"STAT2HEADER_IMP_VERSIONS":
850 str(self.get_versions_of_relevant_modules())})
853 for obj
in listofobjects:
854 if "get_output" not in dir(obj):
856 "Output: object %s doesn't have get_output() method"
861 dfiltered = dict((k, v)
862 for k, v
in d.items()
if k[0] !=
"_")
863 output.update(dfiltered)
866 for obj
in listofsummedobjects:
868 if "get_output" not in dir(t):
870 "Output: object %s doesn't have get_output() method"
873 if "_TotalScore" not in t.get_output():
875 "Output: object %s doesn't have _TotalScore "
876 "entry to be summed" % str(t))
878 output.update({obj[1]: 0.0})
880 for k
in extralabels:
881 output.update({k: 0.0})
883 for n, k
in enumerate(output):
884 stat2_keywords.update({n: k})
885 stat2_inverse.update({k: n})
887 flstat.write(
"%s \n" % stat2_keywords)
889 self.dictionary_stats2[name] = (
895 def write_stat2(self, name, appendmode=True):
897 (listofobjects, stat2_inverse, listofsummedobjects,
898 extralabels) = self.dictionary_stats2[name]
901 for obj
in listofobjects:
902 od = obj.get_output()
903 dfiltered = dict((k, v)
for k, v
in od.items()
if k[0] !=
"_")
905 output.update({stat2_inverse[k]: od[k]})
908 for so
in listofsummedobjects:
912 partial_score += float(d[
"_TotalScore"])
913 output.update({stat2_inverse[so[1]]: str(partial_score)})
916 for k
in extralabels:
917 if k
in self.initoutput:
918 output.update({stat2_inverse[k]: self.initoutput[k]})
920 output.update({stat2_inverse[k]:
"None"})
927 flstat = open(name, writeflag)
928 flstat.write(
"%s \n" % output)
931 def write_stats2(self):
932 for stat
in self.dictionary_stats2.keys():
933 self.write_stat2(stat)
937 """Collect statistics from ProcessOutput.get_fields().
938 Counters of the total number of frames read, plus the models that
939 passed the various filters used in get_fields(), are provided."""
942 self.passed_get_every = 0
943 self.passed_filterout = 0
944 self.passed_filtertuple = 0
948 """A class for reading stat files (either rmf or ascii v1 and v2)"""
949 def __init__(self, filename):
950 self.filename = filename
955 if self.filename
is None:
956 raise ValueError(
"No file name provided. Use -h for help")
960 rh = RMF.open_rmf_file_read_only(self.filename)
962 cat = rh.get_category(
'stat')
963 rmf_klist = rh.get_keys(cat)
964 self.rmf_names_keys = dict([(rh.get_name(k), k)
969 f = open(self.filename,
"r")
972 for line
in f.readlines():
973 d = ast.literal_eval(line)
974 self.klist = list(d.keys())
976 if "STAT2HEADER" in self.klist:
979 if "STAT2HEADER" in str(k):
985 for k
in sorted(stat2_dict.items(),
986 key=operator.itemgetter(1))]
988 for k
in sorted(stat2_dict.items(),
989 key=operator.itemgetter(1))]
990 self.invstat2_dict = {}
992 self.invstat2_dict.update({stat2_dict[k]: k})
995 "statfile v1 is deprecated. "
996 "Please convert to statfile v2.\n")
1005 return sorted(self.rmf_names_keys.keys())
1009 def show_keys(self, ncolumns=2, truncate=65):
1010 IMP.pmi.tools.print_multicolumn(self.get_keys(), ncolumns, truncate)
1012 def get_fields(self, fields, filtertuple=None, filterout=None, get_every=1,
1015 Get the desired field names, and return a dictionary.
1016 Namely, "fields" are the queried keys in the stat file
1017 (eg. ["Total_Score",...])
1018 The returned data structure is a dictionary, where each key is
1019 a field and the value is the time series (ie, frame ordered series)
1020 of that field (ie, {"Total_Score":[Score_0,Score_1,Score_2,,...],....})
1022 @param fields (list of strings) queried keys in the stat file
1023 (eg. "Total_Score"....)
1024 @param filterout specify if you want to "grep" out something from
1025 the file, so that it is faster
1026 @param filtertuple a tuple that contains
1027 ("TheKeyToBeFiltered",relationship,value)
1028 where relationship = "<", "==", or ">"
1029 @param get_every only read every Nth line from the file
1030 @param statistics if provided, accumulate statistics in an
1031 OutputStatistics object
1034 if statistics
is None:
1037 for field
in fields:
1042 rh = RMF.open_rmf_file_read_only(self.filename)
1043 nframes = rh.get_number_of_frames()
1044 for i
in range(nframes):
1045 statistics.total += 1
1047 statistics.passed_get_every += 1
1048 statistics.passed_filterout += 1
1050 if filtertuple
is not None:
1051 keytobefiltered = filtertuple[0]
1052 relationship = filtertuple[1]
1053 value = filtertuple[2]
1054 datavalue = rh.get_root_node().get_value(
1055 self.rmf_names_keys[keytobefiltered])
1056 if self.isfiltered(datavalue, relationship, value):
1059 statistics.passed_filtertuple += 1
1060 for field
in fields:
1061 outdict[field].append(rh.get_root_node().get_value(
1062 self.rmf_names_keys[field]))
1065 f = open(self.filename,
"r")
1068 for line
in f.readlines():
1069 statistics.total += 1
1070 if filterout
is not None:
1071 if filterout
in line:
1073 statistics.passed_filterout += 1
1076 if line_number % get_every != 0:
1077 if line_number == 1
and self.isstat2:
1078 statistics.total -= 1
1079 statistics.passed_filterout -= 1
1081 statistics.passed_get_every += 1
1083 d = ast.literal_eval(line)
1085 print(
"# Warning: skipped line number " + str(line_number)
1086 +
" not a valid line")
1091 if filtertuple
is not None:
1092 keytobefiltered = filtertuple[0]
1093 relationship = filtertuple[1]
1094 value = filtertuple[2]
1095 datavalue = d[keytobefiltered]
1096 if self.isfiltered(datavalue, relationship, value):
1099 statistics.passed_filtertuple += 1
1100 [outdict[field].append(d[field])
for field
in fields]
1103 if line_number == 1:
1104 statistics.total -= 1
1105 statistics.passed_filterout -= 1
1106 statistics.passed_get_every -= 1
1109 if filtertuple
is not None:
1110 keytobefiltered = filtertuple[0]
1111 relationship = filtertuple[1]
1112 value = filtertuple[2]
1113 datavalue = d[self.invstat2_dict[keytobefiltered]]
1114 if self.isfiltered(datavalue, relationship, value):
1117 statistics.passed_filtertuple += 1
1118 [outdict[field].append(d[self.invstat2_dict[field]])
1119 for field
in fields]
1125 def isfiltered(self, datavalue, relationship, refvalue):
1128 _ = float(datavalue)
1130 raise ValueError(
"ProcessOutput.filter: datavalue cannot be "
1131 "converted into a float")
1133 if relationship ==
"<":
1134 if float(datavalue) >= refvalue:
1136 if relationship ==
">":
1137 if float(datavalue) <= refvalue:
1139 if relationship ==
"==":
1140 if float(datavalue) != refvalue:
1146 """ class to allow more advanced handling of RMF files.
1147 It is both a container and a IMP.atom.Hierarchy.
1148 - it is iterable (while loading the corresponding frame)
1149 - Item brackets [] load the corresponding frame
1150 - slice create an iterator
1151 - can relink to another RMF file
1155 @param model: the IMP.Model()
1156 @param rmf_file_name: str, path of the rmf file
1160 self.rh_ref = RMF.open_rmf_file_read_only(rmf_file_name)
1162 raise TypeError(
"Wrong rmf file name or type: %s"
1163 % str(rmf_file_name))
1166 self.root_hier_ref = hs[0]
1167 IMP.atom.Hierarchy.__init__(self, self.root_hier_ref)
1169 self.ColorHierarchy =
None
1173 Link to another RMF file
1175 self.rh_ref = RMF.open_rmf_file_read_only(rmf_file_name)
1177 if self.ColorHierarchy:
1178 self.ColorHierarchy.method()
1179 RMFHierarchyHandler.set_frame(self, 0)
1181 def set_frame(self, index):
1185 print(
"skipping frame %s:%d\n" % (self.current_rmf, index))
1189 return self.rh_ref.get_number_of_frames()
1191 def __getitem__(self, int_slice_adaptor):
1192 if isinstance(int_slice_adaptor, int):
1193 self.set_frame(int_slice_adaptor)
1194 return int_slice_adaptor
1195 elif isinstance(int_slice_adaptor, slice):
1196 return self.__iter__(int_slice_adaptor)
1198 raise TypeError(
"Unknown Type")
1201 return self.get_number_of_frames()
1203 def __iter__(self, slice_key=None):
1204 if slice_key
is None:
1205 for nframe
in range(len(self)):
1208 for nframe
in list(range(len(self)))[slice_key]:
1212 class CacheHierarchyCoordinates(object):
1213 def __init__(self, StatHierarchyHandler):
1220 self.current_index =
None
1221 self.rmfh = StatHierarchyHandler
1223 self.model = self.rmfh.get_model()
1228 self.nrms.append(nrm)
1231 self.xyzs.append(fb)
1233 def do_store(self, index):
1234 self.rb_trans[index] = {}
1235 self.nrm_coors[index] = {}
1236 self.xyz_coors[index] = {}
1238 self.rb_trans[index][rb] = rb.get_reference_frame()
1239 for nrm
in self.nrms:
1240 self.nrm_coors[index][nrm] = nrm.get_internal_coordinates()
1241 for xyz
in self.xyzs:
1242 self.xyz_coors[index][xyz] = xyz.get_coordinates()
1243 self.current_index = index
1245 def do_update(self, index):
1246 if self.current_index != index:
1248 rb.set_reference_frame(self.rb_trans[index][rb])
1249 for nrm
in self.nrms:
1250 nrm.set_internal_coordinates(self.nrm_coors[index][nrm])
1251 for xyz
in self.xyzs:
1252 xyz.set_coordinates(self.xyz_coors[index][xyz])
1253 self.current_index = index
1257 return len(self.rb_trans.keys())
1259 def __getitem__(self, index):
1260 if isinstance(index, int):
1261 return index
in self.rb_trans.keys()
1263 raise TypeError(
"Unknown Type")
1266 return self.get_number_of_frames()
1270 """ class to link stat files to several rmf files """
1271 def __init__(self, model=None, stat_file=None,
1272 number_best_scoring_models=
None, score_key=
None,
1273 StatHierarchyHandler=
None, cache=
None):
1276 @param model: IMP.Model()
1277 @param stat_file: either 1) a list or 2) a single stat file names
1278 (either rmfs or ascii, or pickled data or pickled cluster),
1279 3) a dictionary containing an rmf/ascii
1280 stat file name as key and a list of frames as values
1281 @param number_best_scoring_models:
1282 @param StatHierarchyHandler: copy constructor input object
1283 @param cache: cache coordinates and rigid body transformations.
1286 if StatHierarchyHandler
is not None:
1290 self.model = StatHierarchyHandler.model
1291 self.data = StatHierarchyHandler.data
1292 self.number_best_scoring_models = \
1293 StatHierarchyHandler.number_best_scoring_models
1294 self.is_setup =
True
1295 self.current_rmf = StatHierarchyHandler.current_rmf
1296 self.current_frame =
None
1297 self.current_index =
None
1298 self.score_threshold = StatHierarchyHandler.score_threshold
1299 self.score_key = StatHierarchyHandler.score_key
1300 self.cache = StatHierarchyHandler.cache
1301 RMFHierarchyHandler.__init__(self, self.model,
1304 self.cache = CacheHierarchyCoordinates(self)
1313 self.number_best_scoring_models = number_best_scoring_models
1316 if score_key
is None:
1317 self.score_key =
"Total_Score"
1319 self.score_key = score_key
1320 self.is_setup =
None
1321 self.current_rmf =
None
1322 self.current_frame =
None
1323 self.current_index =
None
1324 self.score_threshold =
None
1326 if isinstance(stat_file, str):
1327 self.add_stat_file(stat_file)
1328 elif isinstance(stat_file, list):
1330 self.add_stat_file(f)
1332 def add_stat_file(self, stat_file):
1334 '''check that it is not a pickle file with saved data
1335 from a previous calculation'''
1336 self.load_data(stat_file)
1338 if self.number_best_scoring_models:
1339 scores = self.get_scores()
1340 max_score = sorted(scores)[
1341 0:min(len(self), self.number_best_scoring_models)][-1]
1342 self.do_filter_by_score(max_score)
1344 except pickle.UnpicklingError:
1345 '''alternatively read the ascii stat files'''
1347 scores, rmf_files, rmf_frame_indexes, features = \
1348 self.get_info_from_stat_file(stat_file,
1349 self.score_threshold)
1350 except (KeyError, SyntaxError):
1355 rh = RMF.open_rmf_file_read_only(stat_file)
1356 nframes = rh.get_number_of_frames()
1357 scores = [0.0]*nframes
1358 rmf_files = [stat_file]*nframes
1359 rmf_frame_indexes = range(nframes)
1364 if len(set(rmf_files)) > 1:
1365 raise (
"Multiple RMF files found")
1368 print(
"StatHierarchyHandler: Error: Trying to set none as "
1369 "rmf_file (probably empty stat file), aborting")
1372 for n, index
in enumerate(rmf_frame_indexes):
1373 featn_dict = dict([(k, features[k][n])
for k
in features])
1375 stat_file, rmf_files[n], index, scores[n], featn_dict))
1377 if self.number_best_scoring_models:
1378 scores = self.get_scores()
1379 max_score = sorted(scores)[
1380 0:min(len(self), self.number_best_scoring_models)][-1]
1381 self.do_filter_by_score(max_score)
1383 if not self.is_setup:
1384 RMFHierarchyHandler.__init__(
1385 self, self.model, self.get_rmf_names()[0])
1387 self.cache = CacheHierarchyCoordinates(self)
1390 self.is_setup =
True
1391 self.current_rmf = self.get_rmf_names()[0]
1395 def save_data(self, filename='data.pkl'):
1396 with open(filename,
'wb')
as fl:
1397 pickle.dump(self.data, fl)
1399 def load_data(self, filename='data.pkl'):
1400 with open(filename,
'rb')
as fl:
1401 data_structure = pickle.load(fl)
1403 if not isinstance(data_structure, list):
1405 "%filename should contain a list of IMP.pmi.output.DataEntry "
1406 "or IMP.pmi.output.Cluster" % filename)
1409 for item
in data_structure):
1410 self.data = data_structure
1412 for item
in data_structure):
1414 for cluster
in data_structure:
1415 nmodels += len(cluster)
1416 self.data = [
None]*nmodels
1417 for cluster
in data_structure:
1418 for n, data
in enumerate(cluster):
1419 index = cluster.members[n]
1420 self.data[index] = data
1423 "%filename should contain a list of IMP.pmi.output.DataEntry "
1424 "or IMP.pmi.output.Cluster" % filename)
1426 def set_frame(self, index):
1427 if self.cache
is not None and self.cache[index]:
1428 self.cache.do_update(index)
1430 nm = self.data[index].rmf_name
1431 fidx = self.data[index].rmf_index
1432 if nm != self.current_rmf:
1434 self.current_rmf = nm
1435 self.current_frame = -1
1436 if fidx != self.current_frame:
1437 RMFHierarchyHandler.set_frame(self, fidx)
1438 self.current_frame = fidx
1439 if self.cache
is not None:
1440 self.cache.do_store(index)
1442 self.current_index = index
1444 def __getitem__(self, int_slice_adaptor):
1445 if isinstance(int_slice_adaptor, int):
1446 self.set_frame(int_slice_adaptor)
1447 return self.data[int_slice_adaptor]
1448 elif isinstance(int_slice_adaptor, slice):
1449 return self.__iter__(int_slice_adaptor)
1451 raise TypeError(
"Unknown Type")
1454 return len(self.data)
1456 def __iter__(self, slice_key=None):
1457 if slice_key
is None:
1458 for i
in range(len(self)):
1461 for i
in range(len(self))[slice_key]:
1464 def do_filter_by_score(self, maximum_score):
1465 self.data = [d
for d
in self.data
if d.score <= maximum_score]
1467 def get_scores(self):
1468 return [d.score
for d
in self.data]
1470 def get_feature_series(self, feature_name):
1471 return [d.features[feature_name]
for d
in self.data]
1473 def get_feature_names(self):
1474 return self.data[0].features.keys()
1476 def get_rmf_names(self):
1477 return [d.rmf_name
for d
in self.data]
1479 def get_stat_files_names(self):
1480 return [d.stat_file
for d
in self.data]
1482 def get_rmf_indexes(self):
1483 return [d.rmf_index
for d
in self.data]
1485 def get_info_from_stat_file(self, stat_file, score_threshold=None):
1489 [stat_file], score_key=self.score_key, feature_keys=fs,
1490 rmf_file_key=
"rmf_file", rmf_file_frame_key=
"rmf_frame_index",
1491 prefiltervalue=score_threshold, get_every=1)
1493 scores = [float(y)
for y
in models[2]]
1494 rmf_files = models[0]
1495 rmf_frame_indexes = models[1]
1496 features = models[3]
1497 return scores, rmf_files, rmf_frame_indexes, features
1502 A class to store data associated to a model
1504 def __init__(self, stat_file=None, rmf_name=None, rmf_index=None,
1505 score=
None, features=
None):
1506 self.rmf_name = rmf_name
1507 self.rmf_index = rmf_index
1509 self.features = features
1510 self.stat_file = stat_file
1513 s =
"IMP.pmi.output.DataEntry\n"
1514 s +=
"---- stat file %s \n" % (self.stat_file)
1515 s +=
"---- rmf file %s \n" % (self.rmf_name)
1516 s +=
"---- rmf index %s \n" % (str(self.rmf_index))
1517 s +=
"---- score %s \n" % (str(self.score))
1518 s +=
"---- number of features %s \n" % (str(len(self.features.keys())))
1524 A container for models organized into clusters
1526 def __init__(self, cid=None):
1527 self.cluster_id = cid
1529 self.precision =
None
1530 self.center_index =
None
1531 self.members_data = {}
1533 def add_member(self, index, data=None):
1534 self.members.append(index)
1535 self.members_data[index] = data
1536 self.average_score = self.compute_score()
1538 def compute_score(self):
1540 score = sum([d.score
for d
in self])/len(self)
1541 except AttributeError:
1546 s =
"IMP.pmi.output.Cluster\n"
1547 s +=
"---- cluster_id %s \n" % str(self.cluster_id)
1548 s +=
"---- precision %s \n" % str(self.precision)
1549 s +=
"---- average score %s \n" % str(self.average_score)
1550 s +=
"---- number of members %s \n" % str(len(self.members))
1551 s +=
"---- center index %s \n" % str(self.center_index)
1554 def __getitem__(self, int_slice_adaptor):
1555 if isinstance(int_slice_adaptor, int):
1556 index = self.members[int_slice_adaptor]
1557 return self.members_data[index]
1558 elif isinstance(int_slice_adaptor, slice):
1559 return self.__iter__(int_slice_adaptor)
1561 raise TypeError(
"Unknown Type")
1564 return len(self.members)
1566 def __iter__(self, slice_key=None):
1567 if slice_key
is None:
1568 for i
in range(len(self)):
1571 for i
in range(len(self))[slice_key]:
1574 def __add__(self, other):
1575 self.members += other.members
1576 self.members_data.update(other.members_data)
1577 self.average_score = self.compute_score()
1578 self.precision =
None
1579 self.center_index =
None
1583 def plot_clusters_populations(clusters):
1586 for cluster
in clusters:
1587 indexes.append(cluster.cluster_id)
1588 populations.append(len(cluster))
1590 import matplotlib.pyplot
as plt
1591 fig, ax = plt.subplots()
1592 ax.bar(indexes, populations, 0.5, color=
'r')
1593 ax.set_ylabel('Population')
1594 ax.set_xlabel((
'Cluster index'))
1598 def plot_clusters_precisions(clusters):
1601 for cluster
in clusters:
1602 indexes.append(cluster.cluster_id)
1604 prec = cluster.precision
1605 print(cluster.cluster_id, prec)
1608 precisions.append(prec)
1610 import matplotlib.pyplot
as plt
1611 fig, ax = plt.subplots()
1612 ax.bar(indexes, precisions, 0.5, color=
'r')
1613 ax.set_ylabel('Precision [A]')
1614 ax.set_xlabel((
'Cluster index'))
1618 def plot_clusters_scores(clusters):
1621 for cluster
in clusters:
1622 indexes.append(cluster.cluster_id)
1624 for data
in cluster:
1625 values[-1].append(data.score)
1628 valuename=
"Scores", positionname=
"Cluster index",
1629 xlabels=
None, scale_plot_length=1.0)
1632 class CrossLinkIdentifierDatabase(object):
1636 def check_key(self, key):
1637 if key
not in self.clidb:
1638 self.clidb[key] = {}
1640 def set_unique_id(self, key, value):
1642 self.clidb[key][
"XLUniqueID"] = str(value)
1644 def set_protein1(self, key, value):
1646 self.clidb[key][
"Protein1"] = str(value)
1648 def set_protein2(self, key, value):
1650 self.clidb[key][
"Protein2"] = str(value)
1652 def set_residue1(self, key, value):
1654 self.clidb[key][
"Residue1"] = int(value)
1656 def set_residue2(self, key, value):
1658 self.clidb[key][
"Residue2"] = int(value)
1660 def set_idscore(self, key, value):
1662 self.clidb[key][
"IDScore"] = float(value)
1664 def set_state(self, key, value):
1666 self.clidb[key][
"State"] = int(value)
1668 def set_sigma1(self, key, value):
1670 self.clidb[key][
"Sigma1"] = str(value)
1672 def set_sigma2(self, key, value):
1674 self.clidb[key][
"Sigma2"] = str(value)
1676 def set_psi(self, key, value):
1678 self.clidb[key][
"Psi"] = str(value)
1680 def get_unique_id(self, key):
1681 return self.clidb[key][
"XLUniqueID"]
1683 def get_protein1(self, key):
1684 return self.clidb[key][
"Protein1"]
1686 def get_protein2(self, key):
1687 return self.clidb[key][
"Protein2"]
1689 def get_residue1(self, key):
1690 return self.clidb[key][
"Residue1"]
1692 def get_residue2(self, key):
1693 return self.clidb[key][
"Residue2"]
1695 def get_idscore(self, key):
1696 return self.clidb[key][
"IDScore"]
1698 def get_state(self, key):
1699 return self.clidb[key][
"State"]
1701 def get_sigma1(self, key):
1702 return self.clidb[key][
"Sigma1"]
1704 def get_sigma2(self, key):
1705 return self.clidb[key][
"Sigma2"]
1707 def get_psi(self, key):
1708 return self.clidb[key][
"Psi"]
1710 def set_float_feature(self, key, value, feature_name):
1712 self.clidb[key][feature_name] = float(value)
1714 def set_int_feature(self, key, value, feature_name):
1716 self.clidb[key][feature_name] = int(value)
1718 def set_string_feature(self, key, value, feature_name):
1720 self.clidb[key][feature_name] = str(value)
1722 def get_feature(self, key, feature_name):
1723 return self.clidb[key][feature_name]
1725 def write(self, filename):
1726 with open(filename,
'wb')
as handle:
1727 pickle.dump(self.clidb, handle)
1729 def load(self, filename):
1730 with open(filename,
'rb')
as handle:
1731 self.clidb = pickle.load(handle)
1735 """Plot the given fields and save a figure as `output`.
1736 The fields generally are extracted from a stat file
1737 using ProcessOutput.get_fields()."""
1738 import matplotlib
as mpl
1740 import matplotlib.pyplot
as plt
1742 plt.rc(
'lines', linewidth=4)
1743 fig, axs = plt.subplots(nrows=len(fields))
1744 fig.set_size_inches(10.5, 5.5 * len(fields))
1749 if framemin
is None:
1751 if framemax
is None:
1752 framemax = len(fields[key])
1753 x = list(range(framemin, framemax))
1754 y = [float(y)
for y
in fields[key][framemin:framemax]]
1757 axs[n].set_title(key, size=
"xx-large")
1758 axs[n].tick_params(labelsize=18, pad=10)
1761 axs.set_title(key, size=
"xx-large")
1762 axs.tick_params(labelsize=18, pad=10)
1766 plt.subplots_adjust(hspace=0.3)
1771 colors=
None, format=
"png", reference_xline=
None,
1772 yplotrange=
None, xplotrange=
None, normalized=
True,
1774 '''Plot a list of histograms from a value list.
1775 @param name the name of the plot
1776 @param value_lists the list of list of values eg: [[...],[...],[...]]
1777 @param valuename the y-label
1778 @param bins the number of bins
1779 @param colors If None, will use rainbow. Else will use specific list
1780 @param format output format
1781 @param reference_xline plot a reference line parallel to the y-axis
1782 @param yplotrange the range for the y-axis
1783 @param normalized whether the histogram is normalized or not
1784 @param leg_names names for the legend
1787 import matplotlib
as mpl
1789 import matplotlib.pyplot
as plt
1790 import matplotlib.cm
as cm
1791 plt.figure(figsize=(18.0, 9.0))
1794 colors = cm.rainbow(np.linspace(0, 1, len(values_lists)))
1795 for nv, values
in enumerate(values_lists):
1797 if leg_names
is not None:
1798 label = leg_names[nv]
1803 [float(y)
for y
in values], bins=bins, color=col,
1804 density=normalized, histtype=
'step', lw=4, label=label)
1805 except AttributeError:
1807 [float(y)
for y
in values], bins=bins, color=col,
1808 normed=normalized, histtype=
'step', lw=4, label=label)
1811 plt.tick_params(labelsize=12, pad=10)
1812 if valuename
is None:
1813 plt.xlabel(name, size=
"xx-large")
1815 plt.xlabel(valuename, size=
"xx-large")
1816 plt.ylabel(
"Frequency", size=
"xx-large")
1818 if yplotrange
is not None:
1820 if xplotrange
is not None:
1821 plt.xlim(xplotrange)
1825 if reference_xline
is not None:
1832 plt.savefig(name +
"." + format, dpi=150, transparent=
True)
1836 valuename=
"None", positionname=
"None",
1837 xlabels=
None, scale_plot_length=1.0):
1839 Plot time series as boxplots.
1840 fields is a list of time series, positions are the x-values
1841 valuename is the y-label, positionname is the x-label
1844 import matplotlib
as mpl
1846 import matplotlib.pyplot
as plt
1849 fig = plt.figure(figsize=(float(len(positions))*scale_plot_length, 5.0))
1850 fig.canvas.manager.set_window_title(name)
1852 ax1 = fig.add_subplot(111)
1854 plt.subplots_adjust(left=0.1, right=0.990, top=0.95, bottom=0.4)
1856 bps.append(plt.boxplot(values, notch=0, sym=
'', vert=1,
1857 whis=1.5, positions=positions))
1859 plt.setp(bps[-1][
'boxes'], color=
'black', lw=1.5)
1860 plt.setp(bps[-1][
'whiskers'], color=
'black', ls=
":", lw=1.5)
1862 if frequencies
is not None:
1863 for n, v
in enumerate(values):
1864 plist = [positions[n]]*len(v)
1865 ax1.plot(plist, v,
'gx', alpha=0.7, markersize=7)
1868 if xlabels
is not None:
1869 ax1.set_xticklabels(xlabels)
1870 plt.xticks(rotation=90)
1871 plt.xlabel(positionname)
1872 plt.ylabel(valuename)
1874 plt.savefig(name +
".pdf", dpi=150)
1878 def plot_xy_data(x, y, title=None, out_fn=None, display=True,
1879 set_plot_yaxis_range=
None, xlabel=
None, ylabel=
None):
1880 import matplotlib
as mpl
1882 import matplotlib.pyplot
as plt
1883 plt.rc(
'lines', linewidth=2)
1885 fig, ax = plt.subplots(nrows=1)
1886 fig.set_size_inches(8, 4.5)
1887 if title
is not None:
1888 fig.canvas.manager.set_window_title(title)
1890 ax.plot(x, y, color=
'r')
1891 if set_plot_yaxis_range
is not None:
1892 x1, x2, y1, y2 = plt.axis()
1893 y1 = set_plot_yaxis_range[0]
1894 y2 = set_plot_yaxis_range[1]
1895 plt.axis((x1, x2, y1, y2))
1896 if title
is not None:
1898 if xlabel
is not None:
1899 ax.set_xlabel(xlabel)
1900 if ylabel
is not None:
1901 ax.set_ylabel(ylabel)
1902 if out_fn
is not None:
1903 plt.savefig(out_fn +
".pdf")
1909 def plot_scatter_xy_data(x, y, labelx="None", labely="None",
1910 xmin=
None, xmax=
None, ymin=
None, ymax=
None,
1911 savefile=
False, filename=
"None.eps", alpha=0.75):
1913 import matplotlib
as mpl
1915 import matplotlib.pyplot
as plt
1916 from matplotlib
import rc
1917 rc(
'font', **{
'family':
'sans-serif',
'sans-serif': [
'Helvetica']})
1919 fig, axs = plt.subplots(1)
1923 axs0.set_xlabel(labelx, size=
"xx-large")
1924 axs0.set_ylabel(labely, size=
"xx-large")
1925 axs0.tick_params(labelsize=18, pad=10)
1929 plot2.append(axs0.plot(x, y,
'o', color=
'k', lw=2, ms=0.1, alpha=alpha,
1939 fig.set_size_inches(8.0, 8.0)
1940 fig.subplots_adjust(left=0.161, right=0.850, top=0.95, bottom=0.11)
1941 if (ymin
is not None)
and (ymax
is not None):
1942 axs0.set_ylim(ymin, ymax)
1943 if (xmin
is not None)
and (xmax
is not None):
1944 axs0.set_xlim(xmin, xmax)
1947 fig.savefig(filename, dpi=300)
1950 def get_graph_from_hierarchy(hier):
1954 (graph, depth, depth_dict) = recursive_graph(
1955 hier, graph, depth, depth_dict)
1958 node_labels_dict = {}
1959 for key
in depth_dict:
1960 if depth_dict[key] < 3:
1961 node_labels_dict[key] = key
1963 node_labels_dict[key] =
""
1964 draw_graph(graph, labels_dict=node_labels_dict)
1967 def recursive_graph(hier, graph, depth, depth_dict):
1970 index = str(hier.get_particle().
get_index())
1971 name1 = nameh +
"|#" + index
1972 depth_dict[name1] = depth
1976 if len(children) == 1
or children
is None:
1978 return (graph, depth, depth_dict)
1982 (graph, depth, depth_dict) = recursive_graph(
1983 c, graph, depth, depth_dict)
1985 index = str(c.get_particle().
get_index())
1986 namec = nameh +
"|#" + index
1987 graph.append((name1, namec))
1990 return (graph, depth, depth_dict)
1993 def draw_graph(graph, labels_dict=None, graph_layout='spring',
1994 node_size=5, node_color=
None, node_alpha=0.3,
1995 node_text_size=11, fixed=
None, pos=
None,
1996 edge_color=
'blue', edge_alpha=0.3, edge_thickness=1,
1998 validation_edges=
None,
1999 text_font=
'sans-serif',
2002 import matplotlib
as mpl
2004 import networkx
as nx
2005 import matplotlib.pyplot
as plt
2006 from math
import sqrt, pi
2012 if isinstance(edge_thickness, list):
2013 for edge, weight
in zip(graph, edge_thickness):
2014 G.add_edge(edge[0], edge[1], weight=weight)
2017 G.add_edge(edge[0], edge[1])
2019 if node_color
is None:
2020 node_color_rgb = (0, 0, 0)
2021 node_color_hex =
"000000"
2026 for node
in G.nodes():
2027 cctuple = cc.rgb(node_color[node])
2028 tmpcolor_rgb.append((cctuple[0]/255,
2031 tmpcolor_hex.append(node_color[node])
2032 node_color_rgb = tmpcolor_rgb
2033 node_color_hex = tmpcolor_hex
2036 if isinstance(node_size, dict):
2038 for node
in G.nodes():
2039 size = sqrt(node_size[node])/pi*10.0
2040 tmpsize.append(size)
2043 for n, node
in enumerate(G.nodes()):
2044 color = node_color_hex[n]
2046 nx.set_node_attributes(
2048 {node: {
'type':
'ellipse',
'w': size,
'h': size,
2049 'fill':
'#' + color,
'label': node}})
2050 nx.set_node_attributes(
2052 {node: {
'type':
'text',
'text': node,
'color':
'#000000',
2053 'visible':
'true'}})
2055 for edge
in G.edges():
2056 nx.set_edge_attributes(
2058 {edge: {
'width': 1,
'fill':
'#000000'}})
2060 for ve
in validation_edges:
2062 if (ve[0], ve[1])
in G.edges():
2063 print(
"found forward")
2064 nx.set_edge_attributes(
2066 {ve: {
'width': 1,
'fill':
'#00FF00'}})
2067 elif (ve[1], ve[0])
in G.edges():
2068 print(
"found backward")
2069 nx.set_edge_attributes(
2071 {(ve[1], ve[0]): {
'width': 1,
'fill':
'#00FF00'}})
2073 G.add_edge(ve[0], ve[1])
2075 nx.set_edge_attributes(
2077 {ve: {
'width': 1,
'fill':
'#FF0000'}})
2081 if graph_layout ==
'spring':
2083 graph_pos = nx.spring_layout(G, k=1.0/8.0, fixed=fixed, pos=pos)
2084 elif graph_layout ==
'spectral':
2085 graph_pos = nx.spectral_layout(G)
2086 elif graph_layout ==
'random':
2087 graph_pos = nx.random_layout(G)
2089 graph_pos = nx.shell_layout(G)
2092 nx.draw_networkx_nodes(G, graph_pos, node_size=node_size,
2093 alpha=node_alpha, node_color=node_color_rgb,
2095 nx.draw_networkx_edges(G, graph_pos, width=edge_thickness,
2096 alpha=edge_alpha, edge_color=edge_color)
2097 nx.draw_networkx_labels(
2098 G, graph_pos, labels=labels_dict, font_size=node_text_size,
2099 font_family=text_font)
2101 plt.savefig(out_filename)
2102 nx.write_gml(G,
'out.gml')
2110 from ipyD3
import d3object
2111 from IPython.display
import display
2113 d3 = d3object(width=800,
2118 title=
'Example table with d3js',
2119 desc=
'An example table created created with d3js with '
2120 'data generated with Python.')
2121 data = [[1277.0, 654.0, 288.0, 1976.0, 3281.0, 3089.0, 10336.0, 4650.0,
2122 4441.0, 4670.0, 944.0, 110.0],
2123 [1318.0, 664.0, 418.0, 1952.0, 3581.0, 4574.0, 11457.0, 6139.0,
2124 7078.0, 6561.0, 2354.0, 710.0],
2125 [1783.0, 774.0, 564.0, 1470.0, 3571.0, 3103.0, 9392.0, 5532.0,
2126 5661.0, 4991.0, 2032.0, 680.0],
2127 [1301.0, 604.0, 286.0, 2152.0, 3282.0, 3369.0, 10490.0, 5406.0,
2128 4727.0, 3428.0, 1559.0, 620.0],
2129 [1537.0, 1714.0, 724.0, 4824.0, 5551.0, 8096.0, 16589.0, 13650.0,
2130 9552.0, 13709.0, 2460.0, 720.0],
2131 [5691.0, 2995.0, 1680.0, 11741.0, 16232.0, 14731.0, 43522.0,
2132 32794.0, 26634.0, 31400.0, 7350.0, 3010.0],
2133 [1650.0, 2096.0, 60.0, 50.0, 1180.0, 5602.0, 15728.0, 6874.0,
2134 5115.0, 3510.0, 1390.0, 170.0],
2135 [72.0, 60.0, 60.0, 10.0, 120.0, 172.0, 1092.0, 675.0, 408.0,
2136 360.0, 156.0, 100.0]]
2137 data = [list(i)
for i
in zip(*data)]
2138 sRows = [[
'January',
2150 sColumns = [[
'Prod {0}'.format(i)
for i
in range(1, 9)],
2151 [
None,
'',
None,
None,
'Group 1',
None,
None,
'Group 2']]
2152 d3.addSimpleTable(data,
2153 fontSizeCells=[12, ],
2156 sRowsMargins=[5, 50, 0],
2157 sColsMargins=[5, 20, 10],
2160 addOutsideBorders=-1,
2164 html = d3.render(mode=[
'html',
'show'])
static bool get_is_setup(const IMP::ParticleAdaptor &p)
A container for models organized into clusters.
A class for reading stat files (either rmf or ascii v1 and v2)
atom::Hierarchies create_hierarchies(RMF::FileConstHandle fh, Model *m)
RMF::FrameID save_frame(RMF::FileHandle file, std::string name="")
Save the current state of the linked objects as a new RMF frame.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
def plot_field_histogram
Plot a list of histograms from a value list.
def plot_fields_box_plots
Plot time series as boxplots.
Utility classes and functions for reading and storing PMI files.
def get_best_models
Given a list of stat files, read them all and find the best models.
A class to store data associated to a model.
void handle_use_deprecated(std::string message)
Break in this method in gdb to find deprecated uses at runtime.
std::string get_module_version()
Return the version of this module, as a string.
void write_pdb(const Selection &mhd, TextOutput out, unsigned int model=1)
Collect statistics from ProcessOutput.get_fields().
def get_fields
Get the desired field names, and return a dictionary.
Warning related to handling of structures.
static bool get_is_setup(Model *m, ParticleIndex pi)
def link_to_rmf
Link to another RMF file.
std::string get_molecule_name_and_copy(atom::Hierarchy h)
Walk up a PMI2 hierarchy/representations and get the "molname.copynum".
The standard decorator for manipulating molecular structures.
Ints get_index(const ParticlesTemp &particles, const Subset &subset, const Subsets &excluded)
def init_pdb
Init PDB Writing.
int get_number_of_frames(const ::npctransport_proto::Assignment &config, double time_step)
A decorator for a particle representing an atom.
Base class for capturing a modeling protocol.
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.
void add_hierarchies(RMF::NodeHandle fh, const atom::Hierarchies &hs)
Class for easy writing of PDBs, RMFs, and stat files.
void add_geometries(RMF::NodeHandle parent, const display::GeometriesTemp &r)
Add geometries to a given parent node.
void add_restraints(RMF::NodeHandle fh, const Restraints &hs)
A decorator for a particle that is part of a rigid body but not rigid.
bool get_is_canonical(atom::Hierarchy h)
Walk up a PMI2 hierarchy/representations and check if the root is named System.
Display a segment connecting a pair of particles.
A decorator for a residue.
Basic functionality that is expected to be used by a wide variety of IMP users.
def get_pdb_names
Get a list of all PDB files being output by this instance.
def get_prot_name_from_particle
Get the protein name from the particle.
class to link stat files to several rmf files
class to allow more advanced handling of RMF files.
void link_hierarchies(RMF::FileConstHandle fh, const atom::Hierarchies &hs)
def plot_fields
Plot the given fields and save a figure as output.
void add_geometry(RMF::FileHandle file, display::Geometry *r)
Add a single geometry to the file.
Store info for a chain of a protein.
Python classes to represent, score, sample and analyze models.
def init_pdb_best_scoring
Prepare for writing best-scoring PDBs (or mmCIFs) for a sampling run.
Functionality for loading, creating, manipulating and scoring atomic structures.
Hierarchies get_leaves(const Selection &h)
Select hierarchy particles identified by the biological name.
def init_rmf
Initialize an RMF file.
static bool get_is_setup(const IMP::ParticleAdaptor &p)
std::string get_module_version()
Return the version of this module, as a string.
A decorator for a particle with x,y,z coordinates and a radius.