diff --git a/examples/demos/the_hammer.py b/examples/demos/the_hammer.py index 48970cefb16a8b9bee9a70395b2d5555f86e82cb..4f6266d287430788fa99d3b56647bba48c9d9588 100644 --- a/examples/demos/the_hammer.py +++ b/examples/demos/the_hammer.py @@ -1,7 +1,5 @@ from PyQt4 import QtCore import math -from ost import qa - # remove all objects from scene, just in case scene.RemoveAll() @@ -26,7 +24,7 @@ class Anim(QtCore.QTimer): self.edi.SetTransform(geom.Mat4(rot)) self.edi.UpdateICS() for a in self.b.view.atoms: - score=qa.ClashScore(a.handle, self.a.view) + score=mol.alg.ClashScore(a.handle, self.a.view) a.SetFloatProp('clash', score) self.a.UpdatePositions() self.b.ReapplyColorOps() diff --git a/modules/conop/pymod/CMakeLists.txt b/modules/conop/pymod/CMakeLists.txt index 57eeb82c986e9b2220fa1e1c188b65d07ca6886b..5affa4b457d780b7368b34a097bec832d2c3f36c 100644 --- a/modules/conop/pymod/CMakeLists.txt +++ b/modules/conop/pymod/CMakeLists.txt @@ -2,6 +2,7 @@ set(OST_CONOP_PYMOD_SOURCES wrap_conop.cc export_builder.cc export_compound.cc + export_amino_acids.cc export_conop.cc export_ring_finder.cc ) diff --git a/modules/conop/pymod/export_amino_acids.cc b/modules/conop/pymod/export_amino_acids.cc new file mode 100644 index 0000000000000000000000000000000000000000..cf06f17324f53667243573edc8b26f6011f9b3c7 --- /dev/null +++ b/modules/conop/pymod/export_amino_acids.cc @@ -0,0 +1,84 @@ +//------------------------------------------------------------------------------ +// This file is part of the OpenStructure project <www.openstructure.org> +// +// Copyright (C) 2008-2011 by the OpenStructure authors +// +// This library is free software; you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License as published by the Free +// Software Foundation; either version 3.0 of the License, or (at your option) +// any later version. +// This library is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with this library; if not, write to the Free Software Foundation, Inc., +// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA +//------------------------------------------------------------------------------ +#include <boost/python.hpp> +#include <boost/python/suite/indexing/vector_indexing_suite.hpp> +#include <ost/conop/amino_acids.hh> + +using namespace boost::python; + +using namespace ost::conop; + +void export_AminoAcids() +{ + enum_<AminoAcid>("AminoAcid") + .value("ALA", ALA) + .value("ARG", ARG) + .value("ASN", ASN) + .value("ASP", ASP) + .value("GLN", GLN) + .value("GLU", GLU) + .value("LYS", LYS) + .value("SER", SER) + .value("CYS", CYS) + .value("TYR", TYR) + .value("THR", THR) + .value("VAL", VAL) + .value("ILE", ILE) + .value("LEU", LEU) + .value("GLY", GLY) + .value("PRO", PRO) + .value("MET", MET) + .value("HIS", HIS) + .value("PHE", PHE) + .value("TRP", TRP) + .export_values() + ; + + class_<AminoAcidSet>("AminoAcidSet", init<>()) + .def("Add", &AminoAcidSet::Add) + .def("Remove", &AminoAcidSet::Remove) + .def("Contains", &AminoAcidSet::Contains) + .def("Empty", &AminoAcidSet::Empty) + .def("CreatePolarSet", + &AminoAcidSet::CreatePolarSet).staticmethod("CreatePolarSet") + .def("CreateApolarSet", + &AminoAcidSet::CreateApolarSet).staticmethod("CreateApolarSet") + .def("CreateAromaticSet", + &AminoAcidSet::CreateAromaticSet).staticmethod("CreateAromaticSet") + .def("CreateThreeStateSet", + &AminoAcidSet::CreateThreeStateSet).staticmethod("CreateThreeStateSet") + .def("CreatePseudoSet", + &AminoAcidSet::CreatePseudoSet).staticmethod("CreatePseudoSet") + .def("CreateCompleteSet", + &AminoAcidSet::CreateCompleteSet).staticmethod("CreateCompleteSet") + .def("CreateSet", &AminoAcidSet::CreateSet).staticmethod("CreateSet") + .def(self_ns::str(self)) + ; + + class_<AminoAcidAlphabet>("AminoAcidAlphabet", init<>()) + .def(vector_indexing_suite<AminoAcidAlphabet>()) + ; + + def("ResidueToAminoAcid",&ResidueToAminoAcid); + def("AminoAcidToResidueName",&AminoAcidToResidueName); + def("OneLetterCodeToResidueName",&OneLetterCodeToResidueName); + def("ResidueNameToOneLetterCode",&ResidueNameToOneLetterCode); + def("OneLetterCodeToAminoAcid",&OneLetterCodeToAminoAcid); +} + diff --git a/modules/conop/src/CMakeLists.txt b/modules/conop/src/CMakeLists.txt index e69c45689648a82ad4d086eb9e1657fa639d1131..65201bea130c811a12bd0a23b9c181e4d7fa795d 100644 --- a/modules/conop/src/CMakeLists.txt +++ b/modules/conop/src/CMakeLists.txt @@ -3,6 +3,7 @@ builder.hh builder_fw.hh conop.hh heuristic_builder.hh +amino_acids.hh compound.hh compound_lib.hh module_config.hh @@ -12,6 +13,7 @@ ring_finder.hh set(OST_CONOP_SOURCES builder.cc +amino_acids.cc conop.cc heuristic_builder.cc compound.cc diff --git a/modules/qa/src/amino_acids.cc b/modules/conop/src/amino_acids.cc similarity index 68% rename from modules/qa/src/amino_acids.cc rename to modules/conop/src/amino_acids.cc index 7111844e3e4035c5fb3d255c23c0c4174c89cb76..7665e97b34240a78c886a06896fe48ef00026ed2 100644 --- a/modules/qa/src/amino_acids.cc +++ b/modules/conop/src/amino_acids.cc @@ -33,7 +33,7 @@ /* Author: Marco Biasini */ -namespace ost { namespace qa { +namespace ost { namespace conop { using namespace boost::spirit; @@ -41,26 +41,26 @@ namespace { struct AminoAcidKeys : public symbols<AminoAcid> { AminoAcidKeys() { - add("ALA", Ala) - ("ARG", Arg) - ("ASN", Asn) - ("ASP", Asp) - ("GLN", Gln) - ("GLU", Glu) - ("LYS", Lys) - ("SER", Ser) - ("CYS", Cys) - ("TYR", Tyr) - ("THR", Thr) - ("VAL", Val) - ("ILE", Ile) - ("LEU", Leu) - ("PRO", Pro) - ("GLY", Gly) - ("MET", Met) - ("HIS", His) - ("TRP", Trp) - ("PHE", Phe); + add("ALA", ALA) + ("ARG", ARG) + ("ASN", ASN) + ("ASP", ASP) + ("GLN", GLN) + ("GLU", GLU) + ("LYS", LYS) + ("SER", SER) + ("CYS", CYS) + ("TYR", TYR) + ("THR", THR) + ("VAL", VAL) + ("ILE", ILE) + ("LEU", LEU) + ("PRO", PRO) + ("GLY", GLY) + ("MET", MET) + ("HIS", HIS) + ("TRP", TRP) + ("PHE", PHE); } }; } // anon namespace @@ -74,51 +74,51 @@ AminoAcid ResidueToAminoAcid(const mol::ResidueHandle& r) { if (aa) return *aa; - return Xxx; + return XXX; } String AminoAcidToResidueName(AminoAcid aa) { switch(aa) { - case Ala: + case ALA: return "ALA"; - case Arg: + case ARG: return "ARG"; - case Asn: + case ASN: return "ASN"; - case Asp: + case ASP: return "ASP"; - case Gln: + case GLN: return "GLN"; - case Glu: + case GLU: return "GLU"; - case Lys: + case LYS: return "LYS"; - case Ser: + case SER: return "SER"; - case Cys: + case CYS: return "CYS"; - case Tyr: + case TYR: return "TYR"; - case Trp: + case TRP: return "TRP"; - case Thr: + case THR: return "THR"; - case Val: + case VAL: return "VAL"; - case Ile: + case ILE: return "ILE"; - case Met: + case MET: return "MET"; - case Leu: + case LEU: return "LEU"; - case Gly: + case GLY: return "GLY"; - case Pro: + case PRO: return "PRO"; - case His: + case HIS: return "HIS"; - case Phe: + case PHE: return "PHE"; default: return "UNK"; @@ -246,85 +246,85 @@ AminoAcid OneLetterCodeToAminoAcid(char olc) char upper_olc=toupper(olc); switch(upper_olc) { case 'A': - return Ala; + return ALA; case 'R': - return Arg; + return ARG; case 'N': - return Asn; + return ASN; case 'D': - return Asp; + return ASP; case 'Q': - return Gln; + return GLN; case 'E': - return Glu; + return GLU; case 'K': - return Lys; + return LYS; case 'S': - return Ser; + return SER; case 'C': - return Cys; + return CYS; case 'Y': - return Tyr; + return TYR; case 'W': - return Trp; + return TRP; case 'T': - return Thr; + return THR; case 'V': - return Val; + return VAL; case 'I': - return Ile; + return ILE; case 'M': - return Met; + return MET; case 'L': - return Leu; + return LEU; case 'G': - return Gly; + return GLY; case 'P': - return Pro; + return PRO; case 'H': - return His; + return HIS; case 'F': - return Phe; + return PHE; default: - return Xxx; + return XXX; } } AminoAcidSet AminoAcidSet::CreatePolarSet() { AminoAcidSet polar; - polar.Add(Arg); - polar.Add(Asn); - polar.Add(Asp); - polar.Add(Gln); - polar.Add(Glu); - polar.Add(Ser); - polar.Add(His); - polar.Add(Lys); - polar.Add(Thr); + polar.Add(ARG); + polar.Add(ASN); + polar.Add(ASP); + polar.Add(GLN); + polar.Add(GLU); + polar.Add(SER); + polar.Add(HIS); + polar.Add(LYS); + polar.Add(THR); return polar; } AminoAcidSet AminoAcidSet::CreateAromaticSet() { AminoAcidSet aromatic; - aromatic.Add(Trp); - aromatic.Add(Tyr); - aromatic.Add(Phe); + aromatic.Add(TRP); + aromatic.Add(TYR); + aromatic.Add(PHE); return aromatic; } AminoAcidSet AminoAcidSet::CreateApolarSet() { AminoAcidSet apolar; - apolar.Add(Ala); - apolar.Add(Val); - apolar.Add(Leu); - apolar.Add(Ile); - apolar.Add(Gly); - apolar.Add(Pro); - apolar.Add(Cys); - apolar.Add(Met); + apolar.Add(ALA); + apolar.Add(VAL); + apolar.Add(LEU); + apolar.Add(ILE); + apolar.Add(GLY); + apolar.Add(PRO); + apolar.Add(CYS); + apolar.Add(MET); return apolar; } @@ -338,26 +338,26 @@ AminoAcidSet AminoAcidSet::CreateSet(AminoAcid aa) std::vector<AminoAcidSet> AminoAcidSet::CreateCompleteSet() { std::vector<AminoAcidSet> alphabet; - alphabet.push_back(AminoAcidSet::CreateSet(Ala)); - alphabet.push_back(AminoAcidSet::CreateSet(Arg)); - alphabet.push_back(AminoAcidSet::CreateSet(Asn)); - alphabet.push_back(AminoAcidSet::CreateSet(Asp)); - alphabet.push_back(AminoAcidSet::CreateSet(Gln)); - alphabet.push_back(AminoAcidSet::CreateSet(Glu)); - alphabet.push_back(AminoAcidSet::CreateSet(Lys)); - alphabet.push_back(AminoAcidSet::CreateSet(Ser)); - alphabet.push_back(AminoAcidSet::CreateSet(Cys)); - alphabet.push_back(AminoAcidSet::CreateSet(Met)); - alphabet.push_back(AminoAcidSet::CreateSet(Trp)); - alphabet.push_back(AminoAcidSet::CreateSet(Tyr)); - alphabet.push_back(AminoAcidSet::CreateSet(Thr)); - alphabet.push_back(AminoAcidSet::CreateSet(Val)); - alphabet.push_back(AminoAcidSet::CreateSet(Ile)); - alphabet.push_back(AminoAcidSet::CreateSet(Leu)); - alphabet.push_back(AminoAcidSet::CreateSet(Gly)); - alphabet.push_back(AminoAcidSet::CreateSet(Pro)); - alphabet.push_back(AminoAcidSet::CreateSet(His)); - alphabet.push_back(AminoAcidSet::CreateSet(Phe)); + alphabet.push_back(AminoAcidSet::CreateSet(ALA)); + alphabet.push_back(AminoAcidSet::CreateSet(ARG)); + alphabet.push_back(AminoAcidSet::CreateSet(ASN)); + alphabet.push_back(AminoAcidSet::CreateSet(ASP)); + alphabet.push_back(AminoAcidSet::CreateSet(GLN)); + alphabet.push_back(AminoAcidSet::CreateSet(GLU)); + alphabet.push_back(AminoAcidSet::CreateSet(LYS)); + alphabet.push_back(AminoAcidSet::CreateSet(SER)); + alphabet.push_back(AminoAcidSet::CreateSet(CYS)); + alphabet.push_back(AminoAcidSet::CreateSet(MET)); + alphabet.push_back(AminoAcidSet::CreateSet(TRP)); + alphabet.push_back(AminoAcidSet::CreateSet(TYR)); + alphabet.push_back(AminoAcidSet::CreateSet(THR)); + alphabet.push_back(AminoAcidSet::CreateSet(VAL)); + alphabet.push_back(AminoAcidSet::CreateSet(ILE)); + alphabet.push_back(AminoAcidSet::CreateSet(LEU)); + alphabet.push_back(AminoAcidSet::CreateSet(GLY)); + alphabet.push_back(AminoAcidSet::CreateSet(PRO)); + alphabet.push_back(AminoAcidSet::CreateSet(HIS)); + alphabet.push_back(AminoAcidSet::CreateSet(PHE)); return alphabet; } @@ -375,26 +375,26 @@ std::vector<AminoAcidSet> AminoAcidSet::CreatePseudoSet() { std::vector<AminoAcidSet> alphabet; AminoAcidSet full_alphabet; - full_alphabet.Add(Ala); - full_alphabet.Add(Val); - full_alphabet.Add(Leu); - full_alphabet.Add(Ile); - full_alphabet.Add(Gly); - full_alphabet.Add(Pro); - full_alphabet.Add(Cys); - full_alphabet.Add(Met); - full_alphabet.Add(Trp); - full_alphabet.Add(Tyr); - full_alphabet.Add(Phe); - full_alphabet.Add(Arg); - full_alphabet.Add(Asn); - full_alphabet.Add(Asp); - full_alphabet.Add(Gln); - full_alphabet.Add(Glu); - full_alphabet.Add(Ser); - full_alphabet.Add(His); - full_alphabet.Add(Lys); - full_alphabet.Add(Thr); + full_alphabet.Add(ALA); + full_alphabet.Add(VAL); + full_alphabet.Add(LEU); + full_alphabet.Add(ILE); + full_alphabet.Add(GLY); + full_alphabet.Add(PRO); + full_alphabet.Add(CYS); + full_alphabet.Add(MET); + full_alphabet.Add(TRP); + full_alphabet.Add(TYR); + full_alphabet.Add(PHE); + full_alphabet.Add(ARG); + full_alphabet.Add(ASN); + full_alphabet.Add(ASP); + full_alphabet.Add(GLN); + full_alphabet.Add(GLU); + full_alphabet.Add(SER); + full_alphabet.Add(HIS); + full_alphabet.Add(LYS); + full_alphabet.Add(THR); alphabet.push_back(full_alphabet); return alphabet; @@ -431,7 +431,7 @@ AminoAcidSet::Iterator AminoAcidSet::Begin() const { // find first bit that is set int start=0; - while (start<=Xxx+1 && !(bits_ & (1 << start))) { + while (start<=XXX+1 && !(bits_ & (1 << start))) { ++start; } return AminoAcidSetIterator(bits_, start); @@ -439,7 +439,7 @@ AminoAcidSet::Iterator AminoAcidSet::Begin() const AminoAcidSet::Iterator AminoAcidSet::End() const { - return AminoAcidSetIterator(bits_, Xxx+1); + return AminoAcidSetIterator(bits_, XXX+1); } bool AminoAcidSet::operator==(const AminoAcidSet& rhs) const diff --git a/modules/qa/src/amino_acids.hh b/modules/conop/src/amino_acids.hh similarity index 81% rename from modules/qa/src/amino_acids.hh rename to modules/conop/src/amino_acids.hh index 03521fa40eb3fcd3c9fa3770d0b73693f39cc2ba..aa778e3c0090f20a56752e9dbf87b2f210e8dc53 100644 --- a/modules/qa/src/amino_acids.hh +++ b/modules/conop/src/amino_acids.hh @@ -16,8 +16,8 @@ // along with this library; if not, write to the Free Software Foundation, Inc., // 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA //------------------------------------------------------------------------------ -#ifndef OST_QA_AMINO_ACIDS_HH -#define OST_QA_AMINO_ACIDS_HH +#ifndef OST_CONOP_AMINO_ACIDS_HH +#define OST_CONOP_AMINO_ACIDS_HH /* Author: Marco Biasini */ @@ -26,35 +26,35 @@ #include <ost/mol/residue_handle.hh> #include <ost/mol/torsion_handle.hh> -#include <ost/qa/module_config.hh> -namespace ost { namespace qa { +#include <ost/conop/module_config.hh> + +namespace ost { namespace conop { // convenience enum for standard set of 20 amino acids // Xxx signifies an unknown amino acid. typedef enum { - Ala, Arg, Asn, - Asp, Gln, Glu, - Lys, Ser, Cys, - Met, Trp, Tyr, - Thr, Val, Ile, - Leu, Gly, Pro, - His, Phe, Xxx - + ALA, ARG, ASN, + ASP, GLN, GLU, + LYS, SER, CYS, + MET, TRP, TYR, + THR, VAL, ILE, + LEU, GLY, PRO, + HIS, PHE, XXX } AminoAcid; /// \brief from residue name to amino acid. /// /// for non standard amino acids Xxx is returned. -DLLEXPORT_OST_QA AminoAcid ResidueToAminoAcid(const mol::ResidueHandle& r); +DLLEXPORT_OST_CONOP AminoAcid ResidueToAminoAcid(const mol::ResidueHandle& r); // \brief from amino acid to residue name -DLLEXPORT_OST_QA String AminoAcidToResidueName(AminoAcid aa); +DLLEXPORT_OST_CONOP String AminoAcidToResidueName(AminoAcid aa); -DLLEXPORT_OST_QA String OneLetterCodeToResidueName(char olc); +DLLEXPORT_OST_CONOP String OneLetterCodeToResidueName(char olc); -DLLEXPORT_OST_QA AminoAcid OneLetterCodeToAminoAcid(char olc); +DLLEXPORT_OST_CONOP AminoAcid OneLetterCodeToAminoAcid(char olc); -char DLLEXPORT_OST_QA ResidueNameToOneLetterCode(String rn); +char DLLEXPORT_OST_CONOP ResidueNameToOneLetterCode(String rn); class AminoAcidSetIterator : public std::iterator<std::forward_iterator_tag, AminoAcid> { @@ -86,16 +86,14 @@ private: void Advance() { ++curr_; - while (curr_<=Xxx && !(bits_ & (1 << curr_))) { ++curr_; } + while (curr_<=XXX && !(bits_ & (1 << curr_))) { ++curr_; } } unsigned int bits_; int curr_; }; - - /// \brief Amino acid bit set -class DLLEXPORT_OST_QA AminoAcidSet { +class DLLEXPORT_OST_CONOP AminoAcidSet { public: typedef AminoAcidSetIterator Iterator; @@ -152,9 +150,11 @@ private: unsigned int bits_; }; -DLLEXPORT_OST_QA std::ostream& operator<<(std::ostream& os, +DLLEXPORT_OST_CONOP std::ostream& operator<<(std::ostream& os, const AminoAcidSet& aa_set); +typedef std::vector<AminoAcidSet> AminoAcidAlphabet; + }} #endif diff --git a/modules/mol/alg/pymod/CMakeLists.txt b/modules/mol/alg/pymod/CMakeLists.txt index 30fae2eb199a4ff4f4c461eefcd5bfc4f00d12c8..a1880895f2096619147137d66cfa94c0480ce7f2 100644 --- a/modules/mol/alg/pymod/CMakeLists.txt +++ b/modules/mol/alg/pymod/CMakeLists.txt @@ -1,6 +1,7 @@ set(OST_MOL_ALG_PYMOD_SOURCES wrap_mol_alg.cc export_svd_superpose.cc + export_clash.cc ) set(OST_MOL_ALG_PYMOD_MODULES diff --git a/modules/qa/pymod/export_clash.cc b/modules/mol/alg/pymod/export_clash.cc similarity index 95% rename from modules/qa/pymod/export_clash.cc rename to modules/mol/alg/pymod/export_clash.cc index 82be0e395f0954cdbdbc227fb6e727ab68de2dcb..0720999be1cf7e65eb42ddb256768f3072f9f2f1 100644 --- a/modules/qa/pymod/export_clash.cc +++ b/modules/mol/alg/pymod/export_clash.cc @@ -17,11 +17,12 @@ // 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA //------------------------------------------------------------------------------ #include <boost/python.hpp> + using namespace boost::python; -#include <ost/qa/clash_score.hh> +#include <ost/mol/alg/clash_score.hh> -using namespace ost::qa; +using namespace ost::mol::alg; using namespace ost; namespace { diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc index c41707f1aea60ab41c726a473ee1840835244506..2d409899673e6feeb6d46caf76d6a1d6b2bbc35b 100644 --- a/modules/mol/alg/pymod/wrap_mol_alg.cc +++ b/modules/mol/alg/pymod/wrap_mol_alg.cc @@ -27,7 +27,7 @@ using namespace boost::python; using namespace ost; void export_svdSuperPose(); - +void export_Clash(); #if OST_IMG_ENABLED void export_entity_to_density(); #endif @@ -56,4 +56,6 @@ BOOST_PYTHON_MODULE(_mol_alg) arg("end")=-1, arg("ref")=-1)); // def("ConstructCBetas", one_arg, args("entity_handle")); def("ConstructCBetas", &ost::mol::alg::ConstructCBetas, (arg("entity_handle"), arg("include_gly")=false)); + + export_Clash(); } diff --git a/modules/mol/alg/src/CMakeLists.txt b/modules/mol/alg/src/CMakeLists.txt index 9ba04a5ee60442edc6ab9522b23222b46a8331f0..f2a23bdd6fbde06023483b25fe210329336a52ee 100644 --- a/modules/mol/alg/src/CMakeLists.txt +++ b/modules/mol/alg/src/CMakeLists.txt @@ -6,10 +6,12 @@ set(OST_MOL_ALG_HEADERS superpose_frames.hh filter_clashes.hh construct_cbeta.hh + clash_score.hh ) set(OST_MOL_ALG_SOURCES svd_superpose.cc + clash_score.cc sec_structure_segments.cc local_dist_test.cc superpose_frames.cc diff --git a/modules/qa/src/clash_score.cc b/modules/mol/alg/src/clash_score.cc similarity index 76% rename from modules/qa/src/clash_score.cc rename to modules/mol/alg/src/clash_score.cc index e942da7e0fd636f219c4e2a43d3fa9787e4c7360..81e97c6322b9bade135b26cd7e33a74a79240eaa 100644 --- a/modules/qa/src/clash_score.cc +++ b/modules/mol/alg/src/clash_score.cc @@ -21,20 +21,20 @@ #include <ost/mol/iterator.hh> #include "clash_score.hh" -namespace ost { namespace qa { +namespace ost { namespace mol { namespace alg { namespace { template <typename T, typename I> -Real do_clash_score(const T& ent_a, const mol::EntityView& ent_b) +Real do_clash_score(const T& ent_a, const EntityView& ent_b) { Real energy=0.0; for (I i=ent_a.AtomsBegin(), e=ent_a.AtomsEnd(); i!=e; ++i) { - mol::AtomViewList clashees=ent_b.FindWithin((*i).GetPos(), + AtomViewList clashees=ent_b.FindWithin((*i).GetPos(), (*i).GetRadius()+1.7); - for (mol::AtomViewList::iterator j=clashees.begin(), + for (AtomViewList::iterator j=clashees.begin(), e2=clashees.end(); j!=e2; ++j) { energy+=StericEnergy((*j).GetPos(), (*j).GetRadius()-0.25, (*i).GetPos(), (*i).GetRadius()-0.25); @@ -62,22 +62,22 @@ Real StericEnergy(const geom::Vec3& pos1, Real r1, return 57.273*(1.0-sqrt(distance_sqr)/rr); } -Real ClashScore(const mol::EntityView& ent_a, const mol::EntityView& ent_b) +Real ClashScore(const EntityView& ent_a, const EntityView& ent_b) { - return do_clash_score<mol::EntityView, mol::AtomViewIter>(ent_a, ent_b); + return do_clash_score<EntityView, AtomViewIter>(ent_a, ent_b); } -Real ClashScore(const mol::EntityHandle& ent_a, const mol::EntityView& ent_b) +Real ClashScore(const EntityHandle& ent_a, const EntityView& ent_b) { - return do_clash_score<mol::EntityHandle, mol::AtomHandleIter>(ent_a, ent_b); + return do_clash_score<EntityHandle, AtomHandleIter>(ent_a, ent_b); } -Real ClashScore(const mol::AtomHandle& atom, const mol::EntityView& ent_b) +Real ClashScore(const AtomHandle& atom, const EntityView& ent_b) { Real energy=0.0; - mol::AtomViewList clashees=ent_b.FindWithin(atom.GetPos(), + AtomViewList clashees=ent_b.FindWithin(atom.GetPos(), atom.GetRadius()+2.0); - for (mol::AtomViewList::iterator j=clashees.begin(), + for (AtomViewList::iterator j=clashees.begin(), e2=clashees.end(); j!=e2; ++j) { energy+=StericEnergy((*j).GetPos(), (*j).GetRadius(), atom.GetPos(), atom.GetRadius()); @@ -85,4 +85,4 @@ Real ClashScore(const mol::AtomHandle& atom, const mol::EntityView& ent_b) return energy; } -}} +}}} diff --git a/modules/qa/src/clash_score.hh b/modules/mol/alg/src/clash_score.hh similarity index 83% rename from modules/qa/src/clash_score.hh rename to modules/mol/alg/src/clash_score.hh index fea13bbeaab037354a044a3ba537b29f95e4d581..ba1700295b341ded0caa6ae727f1ff2c86c5c5e6 100644 --- a/modules/qa/src/clash_score.hh +++ b/modules/mol/alg/src/clash_score.hh @@ -16,8 +16,8 @@ // along with this library; if not, write to the Free Software Foundation, Inc., // 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA //------------------------------------------------------------------------------ -#ifndef OST_QA_CLASH_SCORE_HH -#define OST_QA_CLASH_SCORE_HH +#ifndef OST_MOL_ALG_CLASH_SCORE_HH +#define OST_MOL_ALG_CLASH_SCORE_HH /* Author: Marco Biasini @@ -25,9 +25,9 @@ #include <ost/mol/entity_view.hh> #include <ost/mol/entity_handle.hh> -#include <ost/qa/module_config.hh> +#include <ost/mol/alg/module_config.hh> -namespace ost { namespace qa { +namespace ost { namespace mol { namespace alg { /// \defgroup Clash Steric Clash Score Calculation /// @@ -41,20 +41,20 @@ namespace ost { namespace qa { /// For each atom of ent_a the interaction with atoms of ent_b calculated. /// \return 0.0 if there are no clashes and a positive clash score otherwise. /// \sa \ref the_hammer.py "The Hammer Example" -Real DLLEXPORT_OST_QA ClashScore(const mol::EntityView& ent_a, const mol::EntityView& ent_b); +Real DLLEXPORT_OST_MOL_ALG ClashScore(const mol::EntityView& ent_a, const mol::EntityView& ent_b); /// \brief calculate clash score between full entity and view -Real DLLEXPORT_OST_QA ClashScore(const mol::EntityHandle& ent_a, +Real DLLEXPORT_OST_MOL_ALG ClashScore(const mol::EntityHandle& ent_a, const mol::EntityView& ent_b); //// \brief calculate clash score of one single atom /// /// \return floating point between 0 and 10 /// \sa \ref the_hammer.py "The Hammer Example" -Real DLLEXPORT_OST_QA ClashScore(const mol::AtomHandle& atom, +Real DLLEXPORT_OST_MOL_ALG ClashScore(const mol::AtomHandle& atom, const mol::EntityView& ent_b); /// \brief calculate steric energy of two atoms -Real DLLEXPORT_OST_QA StericEnergy(const geom::Vec3& pos1, Real r1, +Real DLLEXPORT_OST_MOL_ALG StericEnergy(const geom::Vec3& pos1, Real r1, const geom::Vec3& pos2, Real r2); //@} @@ -62,6 +62,6 @@ Real DLLEXPORT_OST_QA StericEnergy(const geom::Vec3& pos1, Real r1, /// /// Dynamic recalculation of clash score for a moving object. The real-valued /// clash score is then color-mapped onto the objects. -}} +}}} #endif diff --git a/modules/qa/CMakeLists.txt b/modules/qa/CMakeLists.txt deleted file mode 100644 index 6955f55f9b8d0237224a2d4aab25aa3b5310003d..0000000000000000000000000000000000000000 --- a/modules/qa/CMakeLists.txt +++ /dev/null @@ -1,2 +0,0 @@ -add_subdirectory(src) -add_subdirectory(pymod) \ No newline at end of file diff --git a/modules/qa/pymod/CMakeLists.txt b/modules/qa/pymod/CMakeLists.txt deleted file mode 100644 index c9681924294521ffbe8119f585f24bb9e7a3f120..0000000000000000000000000000000000000000 --- a/modules/qa/pymod/CMakeLists.txt +++ /dev/null @@ -1,13 +0,0 @@ -set(OST_QA_PYMOD_SOURCES - export_interaction.cc - export_torsion.cc - export_packing.cc - export_reduced.cc - export_utilities.cc - wrap_qa.cc - export_clash.cc -) -set(OST_QA_PYMOD_MODULES -__init__.py -) -pymod(NAME qa CPP ${OST_QA_PYMOD_SOURCES} PY ${OST_QA_PYMOD_MODULES}) diff --git a/modules/qa/pymod/__init__.py b/modules/qa/pymod/__init__.py deleted file mode 100644 index 6e144e870165d695cb921ac8493706179b9bbb25..0000000000000000000000000000000000000000 --- a/modules/qa/pymod/__init__.py +++ /dev/null @@ -1,19 +0,0 @@ -#------------------------------------------------------------------------------ -# This file is part of the OpenStructure project <www.openstructure.org> -# -# Copyright (C) 2008-2011 by the OpenStructure authors -# -# This library is free software; you can redistribute it and/or modify it under -# the terms of the GNU Lesser General Public License as published by the Free -# Software Foundation; either version 3.0 of the License, or (at your option) -# any later version. -# This library is distributed in the hope that it will be useful, but WITHOUT -# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -# details. -# -# You should have received a copy of the GNU Lesser General Public License -# along with this library; if not, write to the Free Software Foundation, Inc., -# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -#------------------------------------------------------------------------------ -from _qa import * \ No newline at end of file diff --git a/modules/qa/pymod/export_interaction.cc b/modules/qa/pymod/export_interaction.cc deleted file mode 100644 index 47f8a1c9ba1f8e543bda4772ec9527b2bd949fc5..0000000000000000000000000000000000000000 --- a/modules/qa/pymod/export_interaction.cc +++ /dev/null @@ -1,269 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#include <boost/python.hpp> -using namespace boost::python; - -#include <ost/qa/interaction_statistics.hh> -#include <ost/qa/all_atom_potential.hh> -#include <boost/python/register_ptr_to_python.hpp> -using namespace ost::qa; -using namespace ost; - -namespace { -typedef uint64_t(InteractionStatistics::*F1)(Real) const; -typedef uint64_t(InteractionStatistics::*F2)(atom::ChemType, - atom::ChemType, Real) const; -F1 f1=&InteractionStatistics::GetCount; -F2 f2=&InteractionStatistics::GetCount; - - - -typedef float (AllAtomPotential::*oneArg)(mol::EntityView, - std::string); -typedef float (AllAtomPotential::*twoArgs)(mol::EntityView, - mol::EntityView, - std::string); - -typedef float (AllAtomPotential::*threeArgs)(atom::ChemType, - atom::ChemType, float); - -oneArg one_arg =&AllAtomPotential::GetTotalEnergy; -twoArgs two_args =&AllAtomPotential::GetTotalEnergy; -threeArgs three_args =&AllAtomPotential::GetEnergy; -} - -void export_Interaction() -{ - enum_<atom::ChemType>("ChemType") - .value("O_GLY", atom::O_GLY) - .value("N_GLY", atom::N_GLY) - .value("C_GLY", atom::C_GLY) - .value("C_GLY_A", atom::C_GLY_A) - - .value("O_ALA", atom::O_ALA) - .value("N_ALA", atom::N_ALA) - .value("C_ALA", atom::C_ALA) - .value("C_ALA_A", atom::C_ALA_A) - .value("C_ALA_B", atom::C_ALA_B) - - .value("O_VAL", atom::O_VAL) - .value("N_VAL", atom::N_VAL) - .value("C_VAL", atom::C_VAL) - .value("C_VAL_A", atom::C_VAL_A) - .value("C_VAL_B", atom::C_VAL_B) - .value("C_VAL_G", atom::C_VAL_G) - - .value("O_LEU", atom::O_LEU) - .value("N_LEU", atom::N_LEU) - .value("C_LEU", atom::C_LEU) - .value("C_LEU_A", atom::C_LEU_A) - .value("C_LEU_B", atom::C_LEU_B) - .value("C_LEU_G", atom::C_LEU_G) - .value("C_LEU_D", atom::C_LEU_D) - - .value("O_ILE", atom::O_ILE) - .value("N_ILE", atom::N_ILE) - .value("C_ILE", atom::C_ILE) - .value("C_ILE_A", atom::C_ILE_A) - .value("C_ILE_B", atom::C_ILE_B) - .value("C_ILE_G1", atom::C_ILE_G1) - .value("C_ILE_G2", atom::C_ILE_G2) - .value("C_ILE_D1", atom::C_ILE_D1) - - .value("O_THR", atom::O_THR) - .value("N_THR", atom::N_THR) - .value("C_THR", atom::C_THR) - .value("C_THR_A", atom::C_THR_A) - .value("C_THR_B", atom::C_THR_B) - .value("O_THR_G1", atom::O_THR_G1) - .value("C_THR_G2", atom::C_THR_G2) - - .value("O_SER", atom::O_SER) - .value("N_SER", atom::N_SER) - .value("C_SER", atom::C_SER) - .value("C_SER_A", atom::C_SER_A) - .value("C_SER_B", atom::C_SER_B) - .value("O_SER_G", atom::O_SER_G) - - .value("O_CYS", atom::O_CYS) - .value("N_CYS", atom::N_CYS) - .value("C_CYS", atom::C_CYS) - .value("C_CYS_A", atom::C_CYS_A) - .value("C_CYS_B", atom::C_CYS_B) - .value("S_CYS_G", atom::S_CYS_G) - - .value("O_MET", atom::O_MET) - .value("N_MET", atom::N_MET) - .value("C_MET", atom::C_MET) - .value("C_MET_A", atom::C_MET_A) - .value("C_MET_B", atom::C_MET_B) - .value("C_MET_G", atom::C_MET_G) - .value("S_MET_D", atom::S_MET_D) - .value("C_MET_E", atom::C_MET_E) - - .value("O_ASP", atom::O_ASP) - .value("N_ASP", atom::N_ASP) - .value("C_ASP", atom::C_ASP) - .value("C_ASP_A", atom::C_ASP_A) - .value("C_ASP_B", atom::C_ASP_B) - .value("C_ASP_G", atom::C_ASP_G) - .value("O_ASP_D", atom::O_ASP_D) - - .value("O_GLU", atom::O_GLU) - .value("N_GLU", atom::N_GLU) - .value("C_GLU", atom::C_GLU) - .value("C_GLU_A", atom::C_GLU_A) - .value("C_GLU_B", atom::C_GLU_B) - .value("C_GLU_G", atom::C_GLU_G) - .value("C_GLU_D", atom::C_GLU_D) - .value("O_GLU_E", atom::O_GLU_E) - - .value("O_ASN", atom::O_ASN) - .value("N_ASN", atom::N_ASN) - .value("C_ASN", atom::C_ASN) - .value("C_ASN_A", atom::C_ASN_A) - .value("C_ASN_B", atom::C_ASN_B) - .value("C_ASN_G", atom::C_ASN_G) - .value("O_ASN_D", atom::O_ASN_D) - .value("N_ASN_D", atom::N_ASN_D) - - .value("O_GLN", atom::O_GLN) - .value("N_GLN", atom::N_GLN) - .value("C_GLN", atom::C_GLN) - .value("C_GLN_A", atom::C_GLN_A) - .value("C_GLN_B", atom::C_GLN_B) - .value("C_GLN_G", atom::C_GLN_G) - .value("C_GLN_D", atom::C_GLN_D) - .value("O_GLN_E", atom::O_GLN_E) - .value("N_GLN_E", atom::N_GLN_E) - - .value("O_LYS", atom::O_LYS) - .value("N_LYS", atom::N_LYS) - .value("C_LYS", atom::C_LYS) - .value("C_LYS_A", atom::C_LYS_A) - .value("C_LYS_B", atom::C_LYS_B) - .value("C_LYS_G", atom::C_LYS_G) - .value("C_LYS_D", atom::C_LYS_D) - .value("C_LYS_E", atom::C_LYS_E) - .value("N_LYS_Z", atom::N_LYS_Z) - - .value("O_ARG", atom::O_ARG) - .value("N_ARG", atom::N_ARG) - .value("C_ARG", atom::C_ARG) - .value("C_ARG_A", atom::C_ARG_A) - .value("C_ARG_B", atom::C_ARG_B) - .value("C_ARG_G", atom::C_ARG_G) - .value("C_ARG_D", atom::C_ARG_D) - .value("N_ARG_E", atom::N_ARG_E) - .value("C_ARG_Z", atom::C_ARG_Z) - .value("N_ARG_H", atom::N_ARG_H) - - .value("O_TYR", atom::O_TYR) - .value("N_TYR", atom::N_TYR) - .value("C_TYR", atom::C_TYR) - .value("C_TYR_A", atom::C_TYR_A) - .value("C_TYR_B", atom::C_TYR_B) - .value("C_TYR_G", atom::C_TYR_G) - .value("C_TYR_D", atom::C_TYR_D) - .value("C_TYR_E", atom::C_TYR_E) - .value("C_TYR_Z", atom::C_TYR_Z) - .value("O_TYR_H", atom::O_TYR_H) - - .value("O_PHE", atom::O_PHE) - .value("N_PHE", atom::N_PHE) - .value("C_PHE", atom::C_PHE) - .value("C_PHE_A", atom::C_PHE_A) - .value("C_PHE_B", atom::C_PHE_B) - .value("C_PHE_G", atom::C_PHE_G) - .value("C_PHE_D", atom::C_PHE_D) - .value("C_PHE_E", atom::C_PHE_E) - .value("C_PHE_Z", atom::C_PHE_Z) - - .value("O_HIS", atom::O_HIS) - .value("N_HIS", atom::N_HIS) - .value("C_HIS", atom::C_HIS) - .value("C_HIS_A", atom::C_HIS_A) - .value("C_HIS_B", atom::C_HIS_B) - .value("C_HIS_G", atom::C_HIS_G) - .value("N_HIS_D1", atom::N_HIS_D1) - .value("C_HIS_D2", atom::C_HIS_D2) - .value("C_HIS_E1", atom::C_HIS_E1) - .value("N_HIS_E2", atom::N_HIS_E2) - - .value("O_TRP", atom::O_TRP) - .value("N_TRP", atom::N_TRP) - .value("C_TRP", atom::C_TRP) - .value("C_TRP_A", atom::C_TRP_A) - .value("C_TRP_B", atom::C_TRP_B) - .value("C_TRP_G", atom::C_TRP_G) - .value("C_TRP_D1", atom::C_TRP_D1) - .value("C_TRP_D2", atom::C_TRP_D2) - .value("N_TRP_E1", atom::N_TRP_E1) - .value("C_TRP_E2", atom::C_TRP_E2) - .value("C_TRP_E3", atom::C_TRP_E3) - .value("C_TRP_Z2", atom::C_TRP_Z2) - .value("C_TRP_Z3", atom::C_TRP_Z3) - .value("C_TRP_H2", atom::C_TRP_H2) - - .value("O_PRO", atom::O_PRO) - .value("N_PRO", atom::N_PRO) - .value("C_PRO", atom::C_PRO) - .value("C_PRO_A", atom::C_PRO_A) - .value("C_PRO_B", atom::C_PRO_B) - .value("C_PRO_G", atom::C_PRO_G) - .value("C_PRO_D", atom::C_PRO_D) - ; - - class_<InteractionStatistics>("InteractionStatistics", - init<Real, Real, Real, int>(args("lower_cutoff","upper_cutoff", - "bucket_size", "sequence_sep"))) - .def("Extract", &InteractionStatistics::Extract, args("view")) - .def("GetCount", f1, args("distance")) - .def("GetCount", f2, args("partner1", "partner2", "distance")) - .def("SaveToFile", &InteractionStatistics::SaveToFile, args("file_name")) - .def("LoadFromFile", &InteractionStatistics::LoadFromFile).staticmethod("LoadFromFile") - ; - - class_<AllAtomPotentialOpts>("AllAtomPotentialOpts", init<>()) - .def(init<float,float,float,int,float>((arg("lower_cutoff"), - arg("upper_cutoff"), arg("dist_bin_size"), - arg("sequence_sep"), arg("sigma")=0.02))) - .def_readwrite("LowerCutoff", &AllAtomPotentialOpts::lower_cutoff) - .def_readwrite("UpperCutoff", &AllAtomPotentialOpts::upper_cutoff) - .def_readwrite("DistanceBucketSize", - &AllAtomPotentialOpts::distance_bucket_size) - .def_readwrite("SequenceSeparation", &AllAtomPotentialOpts::sequence_sep) - .def_readwrite("Sigma", &AllAtomPotentialOpts::sigma) - ; - class_<AllAtomPotential>("AllAtomPotential", no_init) - .def("Create", &AllAtomPotential::Create).staticmethod("Create") - .def("Repair", &AllAtomPotential::Repair).staticmethod("Repair") - .def("LoadFromFile", - &AllAtomPotential::LoadFromFile).staticmethod("LoadFromFile") - .def("GetTotalEnergy", one_arg, args("view")) - .def("GetTotalEnergy", two_args, args("view", "target_view")) - .def("GetEnergyCounts", &AllAtomPotential::GetEnergyCounts) - .def("GetEnergy", three_args, args("atom_type_1", "atom_type_1", "distance")) - .def("SaveToFile", &AllAtomPotential::SaveToFile) - .add_property("options", make_function(&AllAtomPotential::GetOptions, return_value_policy<copy_const_reference>())) - .def("SetSequenceSeparation", &AllAtomPotential::SetSequenceSeparation) - ; - register_ptr_to_python<AllAtomPotentialPtr>(); - register_ptr_to_python<InteractionStatisticsPtr>(); -} diff --git a/modules/qa/pymod/export_packing.cc b/modules/qa/pymod/export_packing.cc deleted file mode 100644 index 607d0618250f66575756a469d13b7ebcc75f7cb7..0000000000000000000000000000000000000000 --- a/modules/qa/pymod/export_packing.cc +++ /dev/null @@ -1,63 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#include <boost/python.hpp> -using namespace boost::python; -#include <boost/python/register_ptr_to_python.hpp> -#include <ost/qa/packing_statistics.hh> -#include <ost/qa/packing_potential.hh> -using namespace ost::qa; -using namespace ost; - -namespace { - typedef uint64_t (PackingStatistics::*F1)(AminoAcid) const; - typedef uint64_t (PackingStatistics::*F2)(AminoAcid, int) const; - F1 f1=&PackingStatistics::GetCount; - F2 f2=&PackingStatistics::GetCount; -} - -void export_Packing() -{ - class_<PackingStatistics>("PackingStatistics", - init<Real,int,int>()) - .def(init<>()) - .def("Extract", &PackingStatistics::Extract, arg("model")) - .def("GetCount", f1, - arg("amino_acid")) - .def("GetCount", f2, - arg("amino_acid"), arg("count")) - .def("LoadFromFile", &PackingStatistics::LoadFromFile).staticmethod("LoadFromFile") - .def("SaveToFile", &PackingStatistics::SaveToFile) - ; - register_ptr_to_python<PackingStatisticsPtr>(); - - class_<PackingPotentialOpts>("PackingPotentialOpts", init<>()) - .def_readwrite("MaxCounts", &PackingPotentialOpts::max_counts) - .def_readwrite("Sigma", &PackingPotentialOpts::sigma) - .def_readwrite("Cutoff", &PackingPotentialOpts::cutoff) - .def_readwrite("BucketSize", &PackingPotentialOpts::bucket_size) - ; - class_<PackingPotential>("PackingPotential", no_init) - .def("Create", &PackingPotential::Create).staticmethod("Create") - .def("LoadFromFile", &PackingPotential::LoadFromFile).staticmethod("LoadFromFile") - .def("GetTotalEnergy", &PackingPotential::GetTotalEnergy) - .def("GetEnergyCounts", &PackingPotential::GetEnergyCounts) - .def("SaveToFile", &PackingPotential::SaveToFile) - ; - register_ptr_to_python<PackingPotentialPtr>(); -} diff --git a/modules/qa/pymod/export_reduced.cc b/modules/qa/pymod/export_reduced.cc deleted file mode 100644 index 14f9d546532e6456cee32cd23ffb60e26cfeb8f7..0000000000000000000000000000000000000000 --- a/modules/qa/pymod/export_reduced.cc +++ /dev/null @@ -1,73 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ - -#include <boost/python.hpp> -#include <ost/mol/mol.hh> -#include <ost/qa/reduced_statistics.hh> -#include <ost/qa/reduced_potential.hh> - -using namespace ost::qa; -using namespace ost::mol; -using namespace boost::python; - -uint64_t (ReducedStatistics::*get_count)(AminoAcid, - AminoAcid, - int,int)=&ReducedStatistics::GetCount; -void (ReducedStatistics::*ex_h)(EntityHandle)=&ReducedStatistics::Extract; -void (ReducedStatistics::*ex_v)(EntityView)=&ReducedStatistics::Extract; - -Real (ReducedPotential::*get_h)(EntityHandle,bool)=&ReducedPotential::GetTotalEnergy; -Real (ReducedPotential::*get_v)(EntityView,bool)=&ReducedPotential::GetTotalEnergy; -void export_Reduced() -{ - class_<ReducedStatOptions>("ReducedStatOptions", init<>()) - .def(init<Real,Real,uint,uint,uint>((arg("l_cutoff"), arg("d_cutoff"), - arg("num_ang_bins"), arg("num_dist_bins"), arg("sequence_sep")))) - .def_readwrite("lower_cutoff", &ReducedStatOptions::lower_cutoff) - .def_readwrite("upper_cutoff", &ReducedStatOptions::upper_cutoff) - .def_readwrite("num_angular_bins", &ReducedStatOptions::num_angular_bins) - .def_readwrite("num_dist_bins", &ReducedStatOptions::num_dist_bins) - .def_readwrite("sequence_sep", &ReducedStatOptions::sequence_sep) - ; - - - class_<ReducedStatistics, ReducedStatisticsPtr>("ReducedStatistics", no_init) - .def(init<Real,Real,uint,uint,uint>((arg("l_cutoff"), arg("u_cutoff"), - arg("num_ang_bins"), arg("num_dist_bins"), arg("sequence_sep")))) - .def("Extract", ex_h, arg("ent")) - .def("Extract", ex_v, arg("ent")) - .add_property("options", make_function(&ReducedStatistics::GetOptions, - return_value_policy<copy_const_reference>())) - .def("Save", &ReducedStatistics::Save) - .def("Load", &ReducedStatistics::Load).staticmethod("Load") - .def("GetCount", get_count) - ; - class_<ReducedPotential, ReducedPotentialPtr>("ReducedPotential", no_init) - .def("Load", &ReducedPotential::Load).staticmethod("Load") - .def("Create", &ReducedPotential::Create).staticmethod("Create") - .def("Save", &ReducedPotential::Save) - .def("GetEnergy", &ReducedPotential::GetEnergy) - .def("GetTotalEnergy", get_h, - (arg("ent"), arg("normalize")=true)) - .def("GetTotalEnergy", get_v, - (arg("ent"), arg("normalize")=true)) - .add_property("options", make_function(&ReducedPotential::GetOptions, - return_value_policy<copy_const_reference>())) - ; -} \ No newline at end of file diff --git a/modules/qa/pymod/export_torsion.cc b/modules/qa/pymod/export_torsion.cc deleted file mode 100644 index 77d76a0d5c860848c470b8f8d47f2b4a9fdf7c59..0000000000000000000000000000000000000000 --- a/modules/qa/pymod/export_torsion.cc +++ /dev/null @@ -1,182 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#include <boost/python.hpp> -using namespace boost::python; - -#include <ost/qa/solis_torsion_statistics.hh> -#include <ost/qa/torsion_statistics.hh> -#include <ost/mol/mol.hh> -#include <boost/python/suite/indexing/vector_indexing_suite.hpp> -#include <boost/python/register_ptr_to_python.hpp> -#include <ost/qa/solis_torsion_potential.hh> -#include <ost/qa/torsion_potential.hh> - -using namespace ost::qa; -using namespace ost; - - -namespace { - typedef Real (TorsionStatisticsSolis::*ADDMethod)(AminoAcid, - Real, Real) const; - typedef Real (TorsionStatisticsSolis::*DDMethod)(Real, Real) const; - typedef void (TorsionStatisticsSolis::*CMethod)(mol::EntityView); - typedef void (TorsionStatisticsSolis::*EMethod)(mol::EntityHandle); - ADDMethod GetF1 = &TorsionStatisticsSolis::GetFrequency; - DDMethod GetF2 = &TorsionStatisticsSolis::GetFrequency; - CMethod Extract1 = &TorsionStatisticsSolis::Extract; - EMethod Extract2 = &TorsionStatisticsSolis::Extract; - -Real (TorsionPotentialSolis::*total_en_view_solis)(mol::EntityView)= - &TorsionPotentialSolis::GetTotalEnergy; -Real (TorsionPotentialSolis::*total_en_handl_solis)(mol::EntityHandle)= - &TorsionPotentialSolis::GetTotalEnergy; - - - - typedef void (TorsionStatistics::*XMethod)(mol::EntityView); - typedef void (TorsionStatistics::*YMethod)(mol::EntityHandle); - XMethod Extract3 = &TorsionStatistics::Extract; - YMethod Extract4 = &TorsionStatistics::Extract; - -Real (TorsionPotential::*total_en_view)(mol::EntityView)= - &TorsionPotential::GetTotalEnergy; -Real (TorsionPotential::*total_en_handl)(mol::EntityHandle)= - &TorsionPotential::GetTotalEnergy; -} - -void export_Torsion() -{ - - enum_<AminoAcid>("AminoAcid") - .value("Ala", Ala) - .value("Arg", Arg) - .value("Asn", Asn) - .value("Asp", Asp) - .value("Gln", Gln) - .value("Glu", Glu) - .value("Lys", Lys) - .value("Ser", Ser) - .value("Cys", Cys) - .value("Tyr", Tyr) - .value("Thr", Thr) - .value("Val", Val) - .value("Ile", Ile) - .value("Leu", Leu) - .value("Gly", Gly) - .value("Pro", Pro) - .value("Met", Met) - .value("His", His) - .value("Phe", Phe) - .value("Trp", Trp) - .export_values() - ; - - class_<AminoAcidSet>("AminoAcidSet", init<>()) - .def("Add", &AminoAcidSet::Add) - .def("Remove", &AminoAcidSet::Remove) - .def("Contains", &AminoAcidSet::Contains) - .def("Empty", &AminoAcidSet::Empty) - .def("CreatePolarSet", - &AminoAcidSet::CreatePolarSet).staticmethod("CreatePolarSet") - .def("CreateApolarSet", - &AminoAcidSet::CreateApolarSet).staticmethod("CreateApolarSet") - .def("CreateAromaticSet", - &AminoAcidSet::CreateAromaticSet).staticmethod("CreateAromaticSet") - .def("CreateThreeStateSet", - &AminoAcidSet::CreateThreeStateSet).staticmethod("CreateThreeStateSet") - .def("CreatePseudoSet", - &AminoAcidSet::CreatePseudoSet).staticmethod("CreatePseudoSet") - .def("CreateCompleteSet", - &AminoAcidSet::CreateCompleteSet).staticmethod("CreateCompleteSet") - .def("CreateSet", &AminoAcidSet::CreateSet).staticmethod("CreateSet") - .def(self_ns::str(self)) - ; - - class_<AminoAcidAlphabet>("AminoAcidAlphabet", init<>()) - .def(vector_indexing_suite<AminoAcidAlphabet>()) - ; - - - class_<TorsionPotentialOptsSolis>("TorsionPotentialOptsSolis", init<>()) - .def_readwrite("AngularBucketSize", - &TorsionPotentialOptsSolis::angular_bucket_size) - .def_readwrite("OuterAlphabet", &TorsionPotentialOptsSolis::outer_alphabet) - .def_readwrite("InnerAlphabet", &TorsionPotentialOptsSolis::inner_alphabet) - ; - - class_<TorsionPotentialOpts>("TorsionPotentialOpts", init<>()) - .def_readwrite("PrevAngularBucketSize", - &TorsionPotentialOpts::prev_angular_bucket_size) - .def_readwrite("CentralAngularBucketSize", - &TorsionPotentialOpts::central_angular_bucket_size) - .def_readwrite("NextAngularBucketSize", - &TorsionPotentialOpts::next_angular_bucket_size) - .def_readwrite("Alphabet", &TorsionPotentialOpts::alphabet) - ; - - - class_<TorsionPotentialSolis, boost::noncopyable>("TorsionPotentialSolis", no_init) - .def("Create", &TorsionPotentialSolis::Create).staticmethod("Create") - .def("LoadFromFile", &TorsionPotentialSolis::LoadFromFile).staticmethod("LoadFromFile") - .def("SaveToFile", &TorsionPotentialSolis::SaveToFile) - .def("GetEnergyCounts", &TorsionPotentialSolis::GetEnergyCounts) - .def("GetTotalEnergy", total_en_view_solis) - .def("GetTotalEnergy", total_en_handl_solis) - ; - - class_<TorsionPotential, boost::noncopyable>("TorsionPotential", no_init) - .def("Create", &TorsionPotential::Create).staticmethod("Create") - .def("LoadFromFile", &TorsionPotential::LoadFromFile).staticmethod("LoadFromFile") - .def("SaveToFile", &TorsionPotential::SaveToFile) - .def("GetEnergyCounts", &TorsionPotential::GetEnergyCounts) - .def("GetTotalEnergy", total_en_view) - .def("GetTotalEnergy", total_en_handl) - ; - - - register_ptr_to_python<TorsionPotentialSolisPtr>(); - register_ptr_to_python<TorsionStatisticsSolisPtr>(); - class_<TorsionStatisticsSolis>("TorsionStatisticsSolis", init<>()) - .def(init<int>(args("torsion_bucket_size"))) - .def("Extract", Extract2, args("model")) - .def("Extract", Extract1, args("chain")) - .def("GetFrequency", GetF1, - arg("amino_acid"), arg("phi_angle"), arg("psi_angle")) - .def("GetFrequency", GetF2, - arg("phi_angle"), arg("psi_angle")) - .def("SaveToFile", &TorsionStatisticsSolis::SaveToFile, arg("file_name")) - .def("LoadFromFile", - &TorsionStatisticsSolis::LoadFromFile, - arg("file_name")).staticmethod("LoadFromFile") - ; - - register_ptr_to_python<TorsionPotentialPtr>(); - register_ptr_to_python<TorsionStatisticsPtr>(); - class_<TorsionStatistics>("TorsionStatistics", init<>()) - .def(init<int,int,int>(args("prev_torsion_bucket_size", - "central_torsion_bucket_size", - "next_torsion_bucket_size"))) - .def("Extract", Extract4, args("model")) - .def("Extract", Extract3, args("chain")) - .def("SaveToFile", &TorsionStatistics::SaveToFile, arg("file_name")) - .def("LoadFromFile", - &TorsionStatistics::LoadFromFile, - arg("file_name")).staticmethod("LoadFromFile") - ; -} diff --git a/modules/qa/pymod/export_utilities.cc b/modules/qa/pymod/export_utilities.cc deleted file mode 100644 index 31a6f2ab8d5bd954b75c9f2829a8753b2840fbc5..0000000000000000000000000000000000000000 --- a/modules/qa/pymod/export_utilities.cc +++ /dev/null @@ -1,34 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#include <boost/python.hpp> -#include <ost/qa/amino_acids.hh> -using namespace boost::python; -using namespace ost::qa; - - -void export_Utilties() -{ - - def ("ResidueToAminoAcid",&ResidueToAminoAcid); - def ("AminoAcidToResidueName",&AminoAcidToResidueName); - def ("OneLetterCodeToResidueName",&OneLetterCodeToResidueName); - def ("ResidueNameToOneLetterCode",&ResidueNameToOneLetterCode); - def ("OneLetterCodeToAminoAcid",&OneLetterCodeToAminoAcid); - -} \ No newline at end of file diff --git a/modules/qa/pymod/pdbtools.py b/modules/qa/pymod/pdbtools.py deleted file mode 100644 index d53d213a74df50c5382971cc5945930a6baec13d..0000000000000000000000000000000000000000 --- a/modules/qa/pymod/pdbtools.py +++ /dev/null @@ -1,59 +0,0 @@ -#------------------------------------------------------------------------------ -# This file is part of the OpenStructure project <www.openstructure.org> -# -# Copyright (C) 2008-2011 by the OpenStructure authors -# -# This library is free software; you can redistribute it and/or modify it under -# the terms of the GNU Lesser General Public License as published by the Free -# Software Foundation; either version 3.0 of the License, or (at your option) -# any later version. -# This library is distributed in the hope that it will be useful, but WITHOUT -# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -# details. -# -# You should have received a copy of the GNU Lesser General Public License -# along with this library; if not, write to the Free Software Foundation, Inc., -# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -#------------------------------------------------------------------------------ -import httplib -import os.path -import types -def GetModelsFromPISCES(pisces_list, skip_header=True): - """ - generator function to iterate over models and chains listed in a PISCES - file. - """ - pisces=open(pisces_list) - if skip_header: - pisces.readline() - for line in pisces: - if line=='': - continue - pdb_id_and_chain=line.split(' ')[0] - yield pdb_id_and_chain[:4], pdb_id_and_chain[4:] - pisces.close() - -def RetrieveModels(models, output_dir, file_pattern='pdb%s.ent'): - """ - Retrieve one or more pdb models from the pdb.org server. - - models must be a pdb id or list of pdb ids - """ - is_tuple=isinstance(models, types.TupleType) - is_list=isinstance(models, types.ListType) - is_generator=isinstance(models, types.GeneratorType) - if not (is_list or is_tuple or is_generator): - models=[models] - conn=httplib.HTTPConnection('www.pdb.org') - for model_id, chain in models: - print 'getting file for model with id=%s...' % model_id.lower(), - conn.request('GET', '/pdb/files/%s.pdb' % model_id ) - response=conn.getresponse() - print response.reason - if response.status==200: - data=response.read() - f=open(os.path.join(output_dir, file_pattern % model_id), 'w+') - f.write(data) - f.close() - conn.close() \ No newline at end of file diff --git a/modules/qa/pymod/wrap_qa.cc b/modules/qa/pymod/wrap_qa.cc deleted file mode 100644 index 14c66d90e3b2ebce817671682f49dc93bea624b0..0000000000000000000000000000000000000000 --- a/modules/qa/pymod/wrap_qa.cc +++ /dev/null @@ -1,36 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#include <boost/python.hpp> - - -void export_Torsion(); -void export_Interaction(); -void export_Packing(); -void export_Clash(); -void export_Reduced(); -void export_Utilties(); -BOOST_PYTHON_MODULE(_qa) -{ - export_Torsion(); - export_Interaction(); - export_Packing(); - export_Clash(); - export_Reduced(); - export_Utilties(); -} diff --git a/modules/qa/src/CMakeLists.txt b/modules/qa/src/CMakeLists.txt deleted file mode 100644 index fee79e13a73d1c7a582c68ce7df81a4c5563f1a6..0000000000000000000000000000000000000000 --- a/modules/qa/src/CMakeLists.txt +++ /dev/null @@ -1,42 +0,0 @@ -set(OST_QA_HEADERS -histogram.hh -index.hh -multi_classifier.hh -amino_acids.hh -torsion_statistics.hh -torsion_potential.hh -solis_torsion_statistics.hh -solis_torsion_potential.hh -interaction_statistics.hh -all_atom_potential.hh -packing_statistics.hh -atom_types.hh -packing_potential.hh -clash_score.hh -reduced_statistics.hh -reduced_potential.hh -module_config.hh -) - -set(OST_QA_SOURCES -reduced_statistics.cc -reduced_potential.cc -torsion_statistics.cc -torsion_potential.cc -solis_torsion_statistics.cc -solis_torsion_potential.cc -interaction_statistics.cc -packing_statistics.cc -amino_acids.cc -atom_types.cc -all_atom_potential.cc -packing_potential.cc -clash_score.cc -impl/reduced_impl.cc -) - - -module(NAME qa SOURCES ${OST_QA_SOURCES} - HEADERS reduced_impl.hh IN_DIR impl ${OST_QA_HEADERS} - DEPENDS_ON io mol_alg) - diff --git a/modules/qa/src/all_atom_potential.cc b/modules/qa/src/all_atom_potential.cc deleted file mode 100644 index 0de014227eee5a568a3cff0826e24b965f485000..0000000000000000000000000000000000000000 --- a/modules/qa/src/all_atom_potential.cc +++ /dev/null @@ -1,334 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#include "all_atom_potential.hh" -#include <ost/mol/mol.hh> -#include <ost/io/binary_data_source.hh> -#include <ost/io/binary_data_sink.hh> -#include <ost/io/container_serialization.hh> -#include <boost/filesystem/convenience.hpp> -#include <ost/message.hh> -#include <ost/io/io_exception.hh> - -/* - Author: Marco Biasini, Pascal Benkert - */ - -namespace ost { namespace qa { - -namespace { - -class AllAtomPotentialCalculator : public mol::EntityVisitor { -public: - AllAtomPotentialCalculator(AllAtomPotential::AllAtomEnergies& energies, - AllAtomPotentialOpts& opts, - mol::EntityView target_view, - int s): - energies_(energies), - options_(opts), - energy_(0.0), - target_view_(target_view), - interaction_counts_(0), - interaction_counts_vector_(std::vector<int>(s, 0)), - score_vector_(std::vector<float>(s, 0.0)), - pos_counter(0) - { - } - - virtual bool VisitResidue(const mol::ResidueHandle& residue) - { - if(curr_aa_!=Xxx) pos_counter++; - curr_aa_=ResidueToAminoAcid(residue); -// std::cout << residue_energy_ << " " << interaction_counts_ << std::endl; - return curr_aa_!=Xxx; - } - - virtual bool VisitAtom(const mol::AtomHandle& atom_handle) - { - atom::ChemType type_a=GetAtomTypeByName(curr_aa_, atom_handle.GetName()); - if (type_a!=atom::UNKNOWN) { - mol::AtomViewList nearby_atoms=target_view_.FindWithin(atom_handle.GetPos(), - options_.upper_cutoff); - for (mol::AtomViewList::iterator i=nearby_atoms.begin(), - e=nearby_atoms.end(); i!=e; ++i) { - mol::AtomView a=*i; - AminoAcid aa_b=ResidueToAminoAcid(a.GetHandle().GetResidue()); - if (aa_b==Xxx) - continue; - atom::ChemType type_b=GetAtomTypeByName(aa_b, a.GetName()); - if (type_b!=atom::UNKNOWN) { - - float d=geom::Distance(atom_handle.GetPos(), - a.GetPos()); - - // avoid rounding problems: distances very close to the cut-off may - // be rounded up the to cut-off and lead to an overflow in container - // energies_ - d=d-0.0001; //used to prevent overflow if distance == upper_cutoff - if (d<0) d=0.0; - - if (d>options_.lower_cutoff) { - // GetHandle() necessary to retriev original index of the view (as - // compared to index in the view): - int abs_sep=abs(atom_handle.GetResidue().GetIndex()- - a.GetResidue().GetHandle().GetIndex()); - if (abs_sep<options_.sequence_sep) { - continue; - } - energy_+=energies_.Get(type_a, type_b, d); - interaction_counts_vector_[pos_counter-1]++; - score_vector_[pos_counter-1]+=energies_.Get(type_a, type_b, d); - interaction_counts_++; - } - } - } - } - return false; - } - - float GetEnergy() const - { - return energy_; - } - - int GetEnergyCounts() const - { - return interaction_counts_; - } - - std::vector<float> GetResidueEnergyVector() const - { - return score_vector_; - } - - std::vector<int> GetResidueCountsVector() const - { - return interaction_counts_vector_; - ; - } - -private: - AllAtomPotential::AllAtomEnergies& energies_; - AllAtomPotentialOpts options_; - AminoAcid curr_aa_; - float energy_; - mol::EntityView target_view_; - int interaction_counts_; - std::vector<int> interaction_counts_vector_; - std::vector<float> score_vector_; - int pos_counter; -}; - -} - -AllAtomPotentialPtr AllAtomPotential::Create(const InteractionStatisticsPtr& s, - const AllAtomPotentialOpts & o) -{ - - AllAtomPotentialPtr p(new AllAtomPotential); - p->options_=o; - int num=int((o.upper_cutoff-o.lower_cutoff)/o.distance_bucket_size); - p->energies_=AllAtomEnergies(0.0, IntegralClassifier(atom::UNKNOWN, 0), - IntegralClassifier(atom::UNKNOWN, 0), - ContinuousClassifier(num, o.lower_cutoff, - o.upper_cutoff)); - p->Fill(s); - return p; -} - - -void AllAtomPotential::SaveToFile(const String& filename) -{ - std::ofstream stream(filename.c_str(), std::ios_base::binary); - io::BinaryDataSink ds(stream); - ds << *this; -} - -AllAtomPotentialPtr AllAtomPotential::LoadFromFile(const String& filename) -{ - if (!boost::filesystem::exists(filename)) { - std::stringstream ss; - ss << "Could not open interaction potential. File '" - << filename << "' does not exist"; - throw io::IOException(ss.str()); - } - - - std::ifstream stream(filename.c_str(), std::ios_base::binary); - io::BinaryDataSource ds(stream); - AllAtomPotentialPtr all_atom_p(new AllAtomPotential); - ds >> *all_atom_p.get(); - return all_atom_p; -} - -float AllAtomPotential::GetTotalEnergy(mol::EntityView view, - mol::EntityView target_view, std::string property_identifier="") -{ - mol::EntityHandle e=view.GetHandle(); - int res_count=view.GetResidueList().size(); - AllAtomPotentialCalculator c(energies_, options_, target_view, res_count); - view.Apply(c); - interaction_counts_=c.GetEnergyCounts(); - - std::vector<float> rev=c.GetResidueEnergyVector(); - std::vector<int> rcv=c.GetResidueCountsVector(); - - if(property_identifier!="") { - for(int i=0;i<res_count;++i) { - view.GetResidueList()[i].SetFloatProp(property_identifier, rev[i]); - view.GetResidueList()[i].SetIntProp(property_identifier+"_counts", rcv[i]); - } - } - return c.GetEnergy(); -} - -float AllAtomPotential::GetTotalEnergy(mol::EntityView view, std::string property_identifier="") -{ - mol::EntityHandle e=view.GetHandle(); - int res_count=view.GetResidueList().size(); - AllAtomPotentialCalculator c(energies_, options_, view, res_count); - view.Apply(c); - interaction_counts_=c.GetEnergyCounts(); - std::string a=""; - std::vector<float> rev=c.GetResidueEnergyVector(); - std::vector<int> rcv=c.GetResidueCountsVector(); - - if(property_identifier!="") { - for(int i=0;i<res_count;++i) { - view.GetResidueList()[i].SetFloatProp(property_identifier, rev[i]); - view.GetResidueList()[i].SetIntProp(property_identifier+"_counts", rcv[i]); - } - } - return c.GetEnergy(); -} - -void AllAtomPotential::SetSequenceSeparation(int seq_sep) { - options_.sequence_sep=seq_sep; -} - - -AllAtomPotentialOpts::AllAtomPotentialOpts(): - sigma(0.02), lower_cutoff(3.0), upper_cutoff(20.0), distance_bucket_size(1.0), - sequence_sep(8) -{ -} - -AllAtomPotentialOpts::AllAtomPotentialOpts(float lc, float uc, - float dbs, int sp, float s): - sigma(s), lower_cutoff(lc), upper_cutoff(uc), distance_bucket_size(dbs), - sequence_sep(sp) -{ - -} - -template <typename DS> -void AllAtomPotential::Serialize(DS& ds) -{ - ds & options_; - ds & energies_; -} - -template <typename DS> -void AllAtomPotentialOpts::Serialize(DS& ds) -{ - ds & lower_cutoff; - ds & upper_cutoff; - ds & sigma; - ds & sequence_sep; - ds & distance_bucket_size; -} - -void AllAtomPotential::Repair() -{ - // for C-apha (backbone only) models fill Calpha interactions with corresponding Cbeta values - typedef AllAtomPotential::AllAtomEnergies::IndexType Index; - int num=int((options_.upper_cutoff- - options_.lower_cutoff)/options_.distance_bucket_size); - - for (int i=0; i<atom::UNKNOWN-1; ++i) { - for (int j=0; j<atom::UNKNOWN-1; ++j) { - - for (int k=0; k<num; ++k) { - float e=energies_.Get(Index(i, j, k)); - - //check if Cbeta (has counts) and not Glycin-Calpha - if(i==0 or j==0){ //mysteriously needed for one case of j=0 (TODO: check later, perhaps in RepairCbetaStatistics()) - continue; - } - if (e!= 0 and i==3 and energies_.Get(Index(i, j+1, k))==0) { - energies_.Set(Index(i, j-1, k), e); -// std::cout << e << " "; - } - else if (e!= 0 and j==3 and energies_.Get(Index(i+1, j, k))==0) { - energies_.Set(Index(i-1, j, k), e); -// std::cout << e << " "; - } - else { - if(e!= 0 and energies_.Get(Index(i, j+1, k))==0 and energies_.Get(Index(i+1, j, k))==0) { -// std::cout << e << " "; - - energies_.Set(Index(i-1, j, k), e); - energies_.Set(Index(i, j-1, k), e); - energies_.Set(Index(i-1, j-1, k), e); - } - } - - } - } - - } -} - -void AllAtomPotential::Fill(const InteractionStatisticsPtr& stats) -{ - typedef AllAtomPotential::AllAtomEnergies::IndexType Index; - int num=int((options_.upper_cutoff- - options_.lower_cutoff)/options_.distance_bucket_size); - - uint64_t total_counts = 0; - for (int k=0; k<num; ++k) { - total_counts = total_counts+stats->GetCount(options_.distance_bucket_size*k+ - options_.lower_cutoff+ - options_.distance_bucket_size*0.5); - } - - for (int i=0; i<atom::UNKNOWN; ++i) { - for (int j=0; j<atom::UNKNOWN; ++j) { - uint64_t t2=stats->GetCount(atom::ChemType(i), atom::ChemType(j)); - - for (int k=0; k<num; ++k) { - uint64_t t3=stats->GetCount(atom::ChemType(i), atom::ChemType(j), - options_.distance_bucket_size*k+ - options_.lower_cutoff+ - Real(options_.distance_bucket_size*0.5)); - // could be moved outside of the loop and cached... - uint64_t t4=stats->GetCount(options_.distance_bucket_size*k+ - options_.lower_cutoff+ - options_.distance_bucket_size*0.5); - float f1=t2 > 0 ?(float(t3)/float(t2)) : 0.0; - float f2=total_counts > 0 ? (float(t4)/float(total_counts)) : 0.0; - float d = f2>0.000000001 ? f1/f2: 0.0; - // prop = (Nd_xy / Nxy) / (Nd / Ntot) - float e=0.582*log(1+options_.sigma*t2)-0.582*log(1+options_.sigma*t2*d); - energies_.Set(Index(i, j, k), e); - } - } - } -} - -}} diff --git a/modules/qa/src/all_atom_potential.hh b/modules/qa/src/all_atom_potential.hh deleted file mode 100644 index 9381e702474f7b36c6ea85d487bb9d67c604450a..0000000000000000000000000000000000000000 --- a/modules/qa/src/all_atom_potential.hh +++ /dev/null @@ -1,116 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#ifndef OST_QA_ALL_ATOM_POTENTIAL_HH -#define OST_QA_ALL_ATOM_POTENTIAL_HH -/* - Author: Marco Biasini - */ - -#include <boost/shared_ptr.hpp> -#include <ost/qa/module_config.hh> -#include <ost/qa/interaction_statistics.hh> - - -namespace ost { namespace qa { - -class AllAtomPotential; -typedef boost::shared_ptr<AllAtomPotential> AllAtomPotentialPtr; - -/// \brief interaction statistics options -struct DLLEXPORT_OST_QA AllAtomPotentialOpts { - - AllAtomPotentialOpts(); - - AllAtomPotentialOpts(float lower_cutoff, float upper_cutoff, - float distance_bucket_size, int sequence_sep, - float sigma=0.02); - - public: - /// \brief weight factor - float sigma; - /// \brief atoms that are closer than the lower cutoff are not considered - float lower_cutoff; - /// \brief atoms that are further apart than the upper cutoff are ignored - float upper_cutoff; - - /// \brief distance bucket size - float distance_bucket_size; - /// \brief sequence separation - int sequence_sep; - - template <typename DS> - void Serialize(DS& ds); -}; - - - -class DLLEXPORT_OST_QA AllAtomPotential { -public: - /// \brief calculate new statistical potential from the given statistics - /// and options - static AllAtomPotentialPtr Create(const InteractionStatisticsPtr& s, - const AllAtomPotentialOpts& o); - - /// \brief load interaction potential from file - static AllAtomPotentialPtr LoadFromFile(const String& filename); - - /// \brief save interaction potential to file - void SaveToFile(const String& filename); - - /// \brief calculate all-atom interaction score for whole entity - float GetTotalEnergy(mol::EntityView view, std::string property_identifier); - - /// \brief extract energy of a specific interaction - /// (for plotting pseudo Lennard-Jones potential). - float GetEnergy(atom::ChemType type_a, atom::ChemType type_b, - float distance) - { - return energies_.Get(type_a, type_b, distance); - } - /// \brief calculate all-atom interaction between two entities. - /// Two entities need to be provided: - /// the atoms for which the energy should be derived and - /// the atoms with respect to which the energy should be calculted. - float GetTotalEnergy(mol::EntityView view, mol::EntityView target_view, std::string property_identifier); - - /// \brief retrieve total number of interactions (for normalisation) - int GetEnergyCounts() const { return interaction_counts_; } - - /// \brief set different seqeunce separation than used for training - void SetSequenceSeparation(int seq_sep); - - const AllAtomPotentialOpts& GetOptions() const { return options_; } - template <typename DS> - void Serialize(DS& ds); -public: - void Repair(); - void Fill(const InteractionStatisticsPtr& stats); - - /// parameters: atom type one, atom type two, distance - typedef MultiClassifier<float, int, int, float> AllAtomEnergies; -private: - AllAtomPotentialOpts options_; - AllAtomEnergies energies_; - mol::EntityView target_view_; - int interaction_counts_; -}; - -}} - -#endif diff --git a/modules/qa/src/atom_types.cc b/modules/qa/src/atom_types.cc deleted file mode 100644 index 7933e93d6de590e9fcdc3fffb6af42fba56bde13..0000000000000000000000000000000000000000 --- a/modules/qa/src/atom_types.cc +++ /dev/null @@ -1,241 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#include <boost/version.hpp> -#if BOOST_VERSION>=103800 -#define BOOST_SPIRIT_USE_OLD_NAMESPACE -#include <boost/spirit/include/classic_symbols.hpp> -#else -#include <boost/spirit/symbols.hpp> -#endif - -#include <map> -#include "atom_types.hh" - -namespace ost { namespace qa { - -using namespace boost::spirit; - - -// let's have a little bit of syntactic sugar to define the atom name to -// atom type mapping. -class ResidueMap : public std::map<AminoAcid, atom::ChemType> { -public: - ResidueMap(AminoAcid amino_acid, atom::ChemType type) { - this->operator()(amino_acid, type); - } - ResidueMap& operator()(AminoAcid amino_acid, atom::ChemType type) { - this->insert(std::make_pair(amino_acid, type)); - return *this; - } - atom::ChemType Get(AminoAcid amino_acid) { - assert(amino_acid!=Xxx); - std::map<AminoAcid, atom::ChemType>::const_iterator i; - i=this->find(amino_acid); - if (i!=this->end()) - return i->second; - return atom::UNKNOWN; - } -}; -typedef ResidueMap R; - - -struct AtomNames : public symbols<R> { - AtomNames() { - add - ("C", R(Ala, atom::C_ALA) - (Arg, atom::C_ARG) - (Asn, atom::C_ASN) - (Asp, atom::C_ASP) - (Gln, atom::C_GLN) - (Glu, atom::C_GLU) - (Lys, atom::C_LYS) - (Ser, atom::C_SER) - (Cys, atom::C_CYS) - (Met, atom::C_MET) - (Trp, atom::C_TRP) - (Tyr, atom::C_TYR) - (Thr, atom::C_THR) - (Val, atom::C_VAL) - (Ile, atom::C_ILE) - (Leu, atom::C_LEU) - (Gly, atom::C_GLY) - (Pro, atom::C_PRO) - (His, atom::C_HIS) - (Phe, atom::C_PHE)) - ("N", R(Ala, atom::N_ALA) - (Arg, atom::N_ARG) - (Asn, atom::N_ASN) - (Asp, atom::N_ASP) - (Gln, atom::N_GLN) - (Glu, atom::N_GLU) - (Lys, atom::N_LYS) - (Ser, atom::N_SER) - (Cys, atom::N_CYS) - (Met, atom::N_MET) - (Trp, atom::N_TRP) - (Tyr, atom::N_TYR) - (Thr, atom::N_THR) - (Val, atom::N_VAL) - (Ile, atom::N_ILE) - (Leu, atom::N_LEU) - (Gly, atom::N_GLY) - (Pro, atom::N_PRO) - (His, atom::N_HIS) - (Phe, atom::N_PHE)) - ("O", R(Ala, atom::O_ALA) - (Arg, atom::O_ARG) - (Asn, atom::O_ASN) - (Asp, atom::O_ASP) - (Gln, atom::O_GLN) - (Glu, atom::O_GLU) - (Lys, atom::O_LYS) - (Ser, atom::O_SER) - (Cys, atom::O_CYS) - (Met, atom::O_MET) - (Trp, atom::O_TRP) - (Tyr, atom::O_TYR) - (Thr, atom::O_THR) - (Val, atom::O_VAL) - (Ile, atom::O_ILE) - (Leu, atom::O_LEU) - (Gly, atom::O_GLY) - (Pro, atom::O_PRO) - (His, atom::O_HIS) - (Phe, atom::O_PHE)) - ("CA", R(Ala, atom::C_ALA_A) - (Arg, atom::C_ARG_A) - (Asn, atom::C_ASN_A) - (Asp, atom::C_ASP_A) - (Gln, atom::C_GLN_A) - (Glu, atom::C_GLU_A) - (Lys, atom::C_LYS_A) - (Ser, atom::C_SER_A) - (Cys, atom::C_CYS_A) - (Met, atom::C_MET_A) - (Trp, atom::C_TRP_A) - (Tyr, atom::C_TYR_A) - (Thr, atom::C_THR_A) - (Val, atom::C_VAL_A) - (Gly, atom::C_GLY_A) - (Ile, atom::C_ILE_A) - (Leu, atom::C_LEU_A) - (Pro, atom::C_PRO_A) - (His, atom::C_HIS_A) - (Phe, atom::C_PHE_A)) - ("CB", R(Ala, atom::C_ALA_B) - (Arg, atom::C_ARG_B) - (Asn, atom::C_ASN_B) - (Asp, atom::C_ASP_B) - (Gln, atom::C_GLN_B) - (Glu, atom::C_GLU_B) - (Lys, atom::C_LYS_B) - (Ser, atom::C_SER_B) - (Cys, atom::C_CYS_B) - (Met, atom::C_MET_B) - (Trp, atom::C_TRP_B) - (Tyr, atom::C_TYR_B) - (Thr, atom::C_THR_B) - (Val, atom::C_VAL_B) - (Ile, atom::C_ILE_B) - (Leu, atom::C_LEU_B) - (Pro, atom::C_PRO_B) - (His, atom::C_HIS_B) - (Phe, atom::C_PHE_B)) - ("CG", R(Arg, atom::C_ARG_G) - (Asn, atom::C_ASN_G) - (Asp, atom::C_ASP_G) - (Gln, atom::C_GLN_G) - (Glu, atom::C_GLU_G) - (Lys, atom::C_LYS_G) - (Met, atom::C_MET_G) - (Trp, atom::C_TRP_G) - (Leu, atom::C_LEU_G) - (Pro, atom::C_PRO_G) - (Phe, atom::C_PHE_G) - (Tyr, atom::C_TYR_G) - (His, atom::C_HIS_G)) - ("CG1", R(Ile, atom::C_ILE_G1) - (Val, atom::C_VAL_G)) - ("CG2", R(Ile, atom::C_ILE_G2) - (Thr, atom::C_THR_G2) - (Val, atom::C_VAL_G)) - ("CD", R(Glu, atom::C_GLU_D) - (Ile, atom::C_GLN_D) - (Lys, atom::C_LYS_D) - (Pro, atom::C_PRO_D) - (Gln, atom::C_GLN_D) - (Glu, atom::C_GLU_D) - (Arg, atom::C_ARG_D)) - ("CD1", R(Trp, atom::C_TRP_D1) - (Tyr, atom::C_TYR_D) - (Phe, atom::C_PHE_D) - (Leu, atom::C_LEU_D) - (Ile, atom::C_ILE_D1)) - ("CD2", R(Trp, atom::C_TRP_D2) - (Tyr, atom::C_TYR_D) - (Phe, atom::C_PHE_D) - (Leu, atom::C_LEU_D) - (His, atom::C_HIS_D2)) - ("CE", R(Lys, atom::C_LYS_E)) - ("CE1", R(Phe, atom::C_PHE_E) - (Tyr, atom::C_TYR_E) - (His, atom::C_HIS_E1)) - ("CE2", R(Phe, atom::C_PHE_E) - (Tyr, atom::C_TYR_E)) - ("CE3", R(Trp, atom::C_TRP_E3)) - ("CZ", R(Arg, atom::C_ARG_Z) - (Phe, atom::C_PHE_Z) - (Tyr, atom::C_TYR_Z)) - ("CZ3", R(Trp, atom::C_TRP_Z3)) - ("OG", R(Ser, atom::O_SER_G)) - ("OG1", R(Thr, atom::O_THR_G1)) - ("SG", R(Cys, atom::S_CYS_G)) - ("SD", R(Met, atom::S_MET_D)) - ("OD1", R(Asp, atom::O_ASP_D) - (Asn, atom::O_ASN_D)) - ("OD2", R(Asp, atom::O_ASP_D)) - ("OE1", R(Glu, atom::O_GLU_E) - (Gln, atom::O_GLN_E)) - ("OE2", R(Glu, atom::O_GLU_E)) - ("ND1", R(His, atom::N_HIS_D1)) - ("ND2", R(Asn, atom::N_ASN_D)) - ("NE", R(Arg, atom::N_ARG_E)) - ("NE1", R(Arg, atom::N_ARG_E)) - ("NE2", R(Gln, atom::N_GLN_E) - (His, atom::N_HIS_E2)) - ("NH1", R(Arg, atom::N_ARG_H)) - ("NH2", R(Arg, atom::N_ARG_H)) - ("NZ", R(Lys, atom::N_LYS_Z)) - ("OH", R(Tyr, atom::O_TYR_H)) - ; - } -}; - -atom::ChemType GetAtomTypeByName(AminoAcid amino_acid, const String& aname) { - - static AtomNames aa_name_symbols; - R* m=find(aa_name_symbols, aname.c_str()); - if (m) { - atom::ChemType t=m->Get(amino_acid); - return t; - } - return atom::UNKNOWN; - -} -}} diff --git a/modules/qa/src/atom_types.hh b/modules/qa/src/atom_types.hh deleted file mode 100644 index 733edb00449008e3aa8b39e922a773eac8aeb207..0000000000000000000000000000000000000000 --- a/modules/qa/src/atom_types.hh +++ /dev/null @@ -1,217 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#ifndef OST_QA_ATOM_TYPES_HH -#define OST_QA_ATOM_TYPES_HH - -#include <ost/qa/amino_acids.hh> - -namespace ost { namespace qa { - -namespace atom { - - // list of chemically distinguishable atoms. - typedef enum { - O_GLY = 0, - N_GLY, - C_GLY, - C_GLY_A, - - O_ALA, - N_ALA, - C_ALA, - C_ALA_A, - C_ALA_B, - - O_VAL, - N_VAL, - C_VAL, - C_VAL_A, - C_VAL_B, - C_VAL_G, - - O_LEU, - N_LEU, - C_LEU, - C_LEU_A, - C_LEU_B, - C_LEU_G, - C_LEU_D, - - O_ILE, - N_ILE, - C_ILE, - C_ILE_A, - C_ILE_B, - C_ILE_G1, - C_ILE_G2, - C_ILE_D1, - - O_THR, - N_THR, - C_THR, - C_THR_A, - C_THR_B, - O_THR_G1, - C_THR_G2, - - O_SER, - N_SER, - C_SER, - C_SER_A, - C_SER_B, - O_SER_G, - - O_CYS, - N_CYS, - C_CYS, - C_CYS_A, - C_CYS_B, - S_CYS_G, - - O_MET, - N_MET, - C_MET, - C_MET_A, - C_MET_B, - C_MET_G, - S_MET_D, - C_MET_E, - - O_ASP, - N_ASP, - C_ASP, - C_ASP_A, - C_ASP_B, - C_ASP_G, - O_ASP_D, - - O_GLU, - N_GLU, - C_GLU, - C_GLU_A, - C_GLU_B, - C_GLU_G, - C_GLU_D, - O_GLU_E, - - O_ASN, - N_ASN, - C_ASN, - C_ASN_A, - C_ASN_B, - C_ASN_G, - O_ASN_D, - N_ASN_D, - - O_GLN, - N_GLN, - C_GLN, - C_GLN_A, - C_GLN_B, - C_GLN_G, - C_GLN_D, - O_GLN_E, - N_GLN_E, - - O_LYS, - N_LYS, - C_LYS, - C_LYS_A, - C_LYS_B, - C_LYS_G, - C_LYS_D, - C_LYS_E, - N_LYS_Z, - - O_ARG, - N_ARG, - C_ARG, - C_ARG_A, - C_ARG_B, - C_ARG_G, - C_ARG_D, - N_ARG_E, - C_ARG_Z, - N_ARG_H, - - O_TYR, - N_TYR, - C_TYR, - C_TYR_A, - C_TYR_B, - C_TYR_G, - C_TYR_D, //D1 AND D2 - C_TYR_E, //E1 and E2 - C_TYR_Z, - O_TYR_H, - - O_PHE, - N_PHE, - C_PHE, - C_PHE_A, - C_PHE_B, - C_PHE_G, - C_PHE_D, //D1 AND D2 - C_PHE_E, //E1 AND E2 - C_PHE_Z, - - O_HIS, - N_HIS, - C_HIS, - C_HIS_A, - C_HIS_B, - C_HIS_G, - N_HIS_D1, - C_HIS_D2, - C_HIS_E1, - N_HIS_E2, - - O_TRP, - N_TRP, - C_TRP, - C_TRP_A, - C_TRP_B, - C_TRP_G, - C_TRP_D1, - C_TRP_D2, - N_TRP_E1, - C_TRP_E2, - C_TRP_E3, - C_TRP_Z2, - C_TRP_Z3, - C_TRP_H2, - - O_PRO, - N_PRO, - C_PRO, - C_PRO_A, - C_PRO_B, - C_PRO_G, - C_PRO_D, - UNKNOWN - } ChemType; -} - - -DLLEXPORT_OST_QA atom::ChemType GetAtomTypeByName(AminoAcid amino_acid, - const String& aname); - -}} - -#endif diff --git a/modules/qa/src/histogram.hh b/modules/qa/src/histogram.hh deleted file mode 100644 index e8abb8ff9d894960d36ea5290b2d51be034794cd..0000000000000000000000000000000000000000 --- a/modules/qa/src/histogram.hh +++ /dev/null @@ -1,78 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#ifndef OST_QA_HISTOGRAM_HH -#define OST_QA_HISTOGRAM_HH -/* - Author: Marco Biasini - */ -#include <ost/stdint.hh> - -#include <ost/qa/module_config.hh> -#include <ost/qa/multi_classifier.hh> - -namespace ost { namespace qa { - -/// \brief histogram specialization for multi classifier -template <typename T1, - typename T2=impl::NullType, - typename T3=impl::NullType, - typename T4=impl::NullType, - typename T5=impl::NullType, - typename T6=impl::NullType, - typename T7=impl::NullType> -class Histogram : public MultiClassifier<uint32_t, T1, T2, T3, T4, T5, T6, T7> -{ -public: - typedef Classifier<T1> C1; - typedef Classifier<T2> C2; - typedef Classifier<T3> C3; - typedef Classifier<T4> C4; - typedef Classifier<T5> C5; - typedef Classifier<T6> C6; - typedef Classifier<T7> C7; -#ifdef _MSC_VER - Histogram(typename C1::ConstRefType c1, - typename C2::ConstRefType c2=C2::Type(), - typename C3::ConstRefType c3=C3::Type(), - typename C4::ConstRefType c4=C4::Type(), - typename C5::ConstRefType c5=C5::Type(), - typename C6::ConstRefType c6=C6::Type(), - typename C7::ConstRefType c7=C7::Type()) -#else - Histogram(typename C1::ConstRefType c1, - typename C2::ConstRefType c2=typename C2::Type(), - typename C3::ConstRefType c3=typename C3::Type(), - typename C4::ConstRefType c4=typename C4::Type(), - typename C5::ConstRefType c5=typename C5::Type(), - typename C6::ConstRefType c6=typename C6::Type(), - typename C7::ConstRefType c7=typename C7::Type()) -#endif - : MultiClassifier<uint32_t, T1, T2, T3, T4, - T5, T6, T7>(0, c1, c2, c3, c4, c5, c6, c7) { - } - - Histogram(): - MultiClassifier<uint32_t, T1, T2, T3, T4, T5, T6, T7>() - { - } -}; - -}} // ost::qa - -#endif diff --git a/modules/qa/src/impl/reduced_impl.cc b/modules/qa/src/impl/reduced_impl.cc deleted file mode 100644 index 59e8c02095a9e1bce00c5afc019b39f9b219b15a..0000000000000000000000000000000000000000 --- a/modules/qa/src/impl/reduced_impl.cc +++ /dev/null @@ -1,113 +0,0 @@ -#include <ost/log.hh> -#include <ost/mol/residue_handle.hh> -#include <ost/mol/atom_handle.hh> -#include <ost/mol/atom_view.hh> -#include <ost/mol/alg/construct_cbeta.hh> -#include "reduced_impl.hh" - - -namespace ost { namespace qa { namespace impl { - -bool ReducedPotentialImpl::VisitResidue(const mol::ResidueHandle& res) -{ - if (!res.IsPeptideLinking()) { - return false; - } - AminoAcid aa_one=OneLetterCodeToAminoAcid(res.GetOneLetterCode()); - if (aa_one==Xxx) { - return false; - } - geom::Vec3 ca_pos_one; - geom::Vec3 cb_pos_one; - uint index=res.GetIndex(); - if (!this->GetCAlphaCBetaPos(res, ca_pos_one, cb_pos_one)) { - return false; - } - - if (ent_) { - // we got a full entity handle. - mol::AtomHandleList within=ent_.FindWithin(ca_pos_one, - opts_.upper_cutoff-0.00001); - for (mol::AtomHandleList::const_iterator - i=within.begin(), e=within.end(); i!=e; ++i) { - if (i->GetName()=="CA") { - this->HandleResidue(aa_one, ca_pos_one, cb_pos_one, index, *i); - } - - } - return false; - } - mol::AtomViewList within=view_.FindWithin(ca_pos_one, - opts_.upper_cutoff-0.00001); - for (mol::AtomViewList::const_iterator - i=within.begin(), e=within.end(); i!=e; ++i) { - if (i->GetName()=="CA") { - this->HandleResidue(aa_one, ca_pos_one, cb_pos_one, index, i->GetHandle()); - } - - } - return false; -} - -void ReducedPotentialImpl::HandleResidue(AminoAcid aa_one, - const geom::Vec3& ca_pos_one, - const geom::Vec3& cb_pos_one, - uint index_one, - const mol::AtomHandle& ca_two) -{ - - mol::ResidueHandle res_two=ca_two.GetResidue(); - if (res_two.GetIndex()-index_one<opts_.sequence_sep) { - return; - } - if (!res_two.IsPeptideLinking()) { - return; - } - - - AminoAcid aa_two=OneLetterCodeToAminoAcid(res_two.GetOneLetterCode()); - if (aa_two==Xxx) { - return; - } - geom::Vec3 ca_pos_two; - geom::Vec3 cb_pos_two; - if (!this->GetCAlphaCBetaPos(res_two, ca_pos_two, cb_pos_two)) { - return; - } - - Real dist=geom::Length(ca_pos_one-ca_pos_two); - if (dist<opts_.lower_cutoff) { - return; - } - Real angle=geom::Angle(cb_pos_one-ca_pos_one, cb_pos_two-ca_pos_two); - - this->OnInteraction(aa_one, aa_two, dist, angle); -} - -bool ReducedPotentialImpl::GetCAlphaCBetaPos(const mol::ResidueHandle& res, - geom::Vec3& ca_pos, - geom::Vec3& cb_pos) -{ - mol::AtomHandle ca=res.FindAtom("CA"); - if (!ca.IsValid()) { - return false; - } - ca_pos=ca.GetPos(); - mol::AtomHandle cb=res.FindAtom("CB"); - if (cb.IsValid()) { - cb_pos=cb.GetPos(); - return true; - } - mol::AtomHandle n=res.FindAtom("N"); - mol::AtomHandle c=res.FindAtom("C"); - if (!(ca.IsValid() && c.IsValid() && n.IsValid())) { - LOG_WARNING("residue " << res.GetQualifiedName() - << " doesn't have enough atoms to reconstruct Cbeta position"); - return false; - } - cb_pos=mol::alg::CBetaPosition(n.GetPos(), ca.GetPos(), c.GetPos()); - return true; -} - -}}} - diff --git a/modules/qa/src/impl/reduced_impl.hh b/modules/qa/src/impl/reduced_impl.hh deleted file mode 100644 index 1d1157228733a2157681e3d8bcb46a368d1bd085..0000000000000000000000000000000000000000 --- a/modules/qa/src/impl/reduced_impl.hh +++ /dev/null @@ -1,63 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ - -#ifndef OST_QA_IMPL_REDUCED_HH -#define OST_QA_IMPL_REDUCED_HH - - -#include <ost/mol/entity_visitor.hh> -#include <ost/mol/entity_handle.hh> -#include <ost/mol/entity_view.hh> -#include <ost/qa/module_config.hh> -#include <ost/qa/amino_acids.hh> -#include <ost/qa/reduced_statistics.hh> - -namespace ost { namespace qa { namespace impl { - -class DLLEXPORT_OST_QA ReducedPotentialImpl : public mol::EntityVisitor { -public: - - ReducedPotentialImpl(const ReducedStatOptions& opts, mol::EntityHandle ent): - opts_(opts), ent_(ent) - { } - - ReducedPotentialImpl(const ReducedStatOptions& opts, mol::EntityView ent): - opts_(opts), view_(ent) - { } - - virtual bool VisitResidue(const mol::ResidueHandle& res); - - static bool GetCAlphaCBetaPos(const mol::ResidueHandle& res, - geom::Vec3& ca_pos, - geom::Vec3& cb_pos); - - virtual void OnInteraction(AminoAcid aa_one, AminoAcid aa_two, - Real dist, Real angle)=0; -private: - void HandleResidue(AminoAcid aa_one, const geom::Vec3& ca_pos_one, - const geom::Vec3& cb_pos_one, - uint index_one, const mol::AtomHandle& ca_two); - ReducedStatOptions opts_; - mol::EntityHandle ent_; - mol::EntityView view_; -}; - -}}} - -#endif diff --git a/modules/qa/src/index.hh b/modules/qa/src/index.hh deleted file mode 100644 index 8ae741851125060775b8afbd9a8b6eee8db51163..0000000000000000000000000000000000000000 --- a/modules/qa/src/index.hh +++ /dev/null @@ -1,185 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#ifndef OST_QA_INDEX_HH -#define OST_QA_INDEX_HH - -/* - Author: Marco Biasini - */ - -#include <ost/stdint.hh> - -#include <cstring> -#include <boost/shared_ptr.hpp> -#include <ost/qa/module_config.hh> - -namespace ost { namespace qa { - -namespace impl { -template <uint32_t D> -class IndexBase { -public: - enum { Dimension = D }; - IndexBase(const IndexBase& rhs) { - memcpy(data_, rhs.data_, sizeof(uint32_t[D])); - } - IndexBase() { - memset(data_, 0, sizeof(uint32_t[D])); - } - IndexBase& operator=(const IndexBase& rhs) { - memcpy(data_, rhs.data_, sizeof(uint32_t[D])); - return *this; - } - uint32_t operator[](uint32_t idx) const { - assert(idx <= D); - return data_[idx]; - } - uint32_t& operator[](uint32_t idx) { - assert(idx <= D); - return data_[idx]; - } -private: - uint32_t data_[D]; -}; - -} // namespace impl - -template <uint32_t D> -class Index; - -template <> -class Index<1> : public impl::IndexBase<1> { -public: - Index() : impl::IndexBase<1>() {} - Index(uint32_t a) { - (*this)[0]=a; - } -}; -template <> -class Index<2> : public impl::IndexBase<2> { -public: - Index() : impl::IndexBase<2>() {} - Index(uint32_t a, uint32_t b) { - (*this)[0]=a; - (*this)[1]=b; - } -}; - -template <> -class Index<3> : public impl::IndexBase<3> { -public: - Index() : impl::IndexBase<3>() {} - Index(uint32_t a, uint32_t b, uint32_t c) { - (*this)[0]=a; - (*this)[1]=b; - (*this)[2]=c; - } -}; - -template <> -class Index<4> : public impl::IndexBase<4> { -public: - Index() : impl::IndexBase<4>() {} - Index(uint32_t a, uint32_t b, uint32_t c, uint32_t d) { - (*this)[0]=a; - (*this)[1]=b; - (*this)[2]=c; - (*this)[3]=d; - } -}; - -template <> -class Index<5> : public impl::IndexBase<5> { -public: - Index() : impl::IndexBase<5>() {} - Index(uint32_t a, uint32_t b, uint32_t c, uint32_t d, uint32_t e) { - (*this)[0]=a; - (*this)[1]=b; - (*this)[2]=c; - (*this)[3]=d; - (*this)[4]=e; - } -}; - -template <> -class Index<6> : public impl::IndexBase<6> { -public: - Index() : impl::IndexBase<6>() {} - Index(uint32_t a, uint32_t b, uint32_t c, uint32_t d, uint32_t e, uint32_t f) { - (*this)[0]=a; - (*this)[1]=b; - (*this)[2]=c; - (*this)[3]=d; - (*this)[4]=e; - (*this)[5]=f; - } -}; - -template <> -class Index<7> : public impl::IndexBase<7> { -public: - Index() : impl::IndexBase<7>() {} - Index(uint32_t a, uint32_t b, uint32_t c, uint32_t d, uint32_t e, uint32_t f, uint32_t g) { - (*this)[0]=a; - (*this)[1]=b; - (*this)[2]=c; - (*this)[3]=d; - (*this)[4]=e; - (*this)[5]=f; - (*this)[6]=g; - } -}; - -template<uint32_t D> -class IndexIterator { -public: - typedef Index<D> IndexType; - IndexIterator(const IndexType& s, const IndexType& e) - : start_(s), end_(e), current_(s) { - - } - - IndexIterator<D>& operator++() { - uint32_t current_it=0; - while (++current_[current_it] > end_[current_it]) { - current_it++; - if (current_it < D) { - current_[current_it-1] = start_[current_it-1]; - } else { - break; - } - } - return *this; - } - const IndexType& operator *() const { - return current_; - } - bool AtEnd() { - return current_[D-1] > end_[D-1]; - } -private: - IndexType start_; - IndexType end_; - IndexType current_; - -}; - -}} - -#endif diff --git a/modules/qa/src/interaction_statistics.cc b/modules/qa/src/interaction_statistics.cc deleted file mode 100644 index c30fcf132999e9df487a133a8ba4ddd47128ba41..0000000000000000000000000000000000000000 --- a/modules/qa/src/interaction_statistics.cc +++ /dev/null @@ -1,194 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#include "interaction_statistics.hh" -#include <ost/mol/atom_view.hh> -#include <ost/mol/residue_view.hh> -#include <ost/io/binary_data_source.hh> -#include <ost/io/binary_data_sink.hh> -#include <ost/io/container_serialization.hh> - -using namespace ost::mol::impl; - -namespace ost { namespace qa { - -InteractionStatistics::InteractionStatistics() -{ - -} - -void InteractionStatistics::Extract(mol::EntityView view_a, mol::EntityView view_b) { - view_a_=view_a; - view_b_=view_b; - view_a.Apply(*this); -} - -bool InteractionStatistics::VisitResidue(const mol::ResidueHandle& r) { - amino_acid_=ResidueToAminoAcid(r); - return amino_acid_!=Xxx; -} - -bool InteractionStatistics::VisitAtom(const mol::AtomHandle& atom) { - atom::ChemType atype=GetAtomTypeByName(amino_acid_, atom.GetName()); - mol::ResidueHandle my_res=atom.GetResidue(); - if (atype!=atom::UNKNOWN) { - mol::AtomViewList atoms_in_reach=view_b_.FindWithin(atom.GetPos(), - upper_cutoff_); - for (mol::AtomViewList::iterator i=atoms_in_reach.begin(), - e=atoms_in_reach.end(); i!=e; ++i) { - mol::AtomView other=*i; - // avoid rounding problems: distances very close to the cut-off may - // be rounded up the to cut-off and lead to an overflow in histogram_ - Real dist_sqr=geom::Length2(other.GetPos()-atom.GetPos())-0.0001; - if (dist_sqr<lower_sqr_) { - continue; - } - mol::ResidueView other_res=other.GetResidue(); - AminoAcid other_aa=ResidueToAminoAcid(other_res.GetHandle()); - if (other_aa==Xxx) { - continue; - } - if (abs(other_res.GetNumber().GetNum()-my_res.GetNumber().GetNum())<sequence_sep_) { - continue; - } - atom::ChemType btype=GetAtomTypeByName(other_aa, other.GetName()); - if (btype==atom::UNKNOWN) { - continue; - } - histogram_.Add(1, atype, btype, sqrt(dist_sqr)); - } - } - return true; -} - -InteractionStatistics::InteractionStatistics(Real lower_cutoff, - Real upper_cutoff, - Real bucket_size, - int sequence_sep) - : lower_cutoff_(lower_cutoff), upper_cutoff_(upper_cutoff), - bucket_size_(bucket_size), - sequence_sep_(sequence_sep), - histogram_(IntegralClassifier(atom::UNKNOWN, 0), - IntegralClassifier(atom::UNKNOWN, 0), - ContinuousClassifier(int((upper_cutoff_-lower_cutoff_)/bucket_size), - lower_cutoff_, upper_cutoff_)) -{ - upper_sqr_=upper_cutoff_*upper_cutoff_; - lower_sqr_=lower_cutoff_*lower_cutoff_; -} - -int InteractionStatistics::GetSequenceSeparation() const -{ - return sequence_sep_; -} - -Real InteractionStatistics::GetUpperCutoff() const -{ - return upper_cutoff_; -} - -Real InteractionStatistics::GetLowerCutoff() const -{ - return lower_cutoff_; -} - -void InteractionStatistics::SaveToFile(const String& file_name) const -{ - std::ofstream stream(file_name.c_str(), std::ios_base::binary); - io::BinaryDataSink ds(stream); - ds << *this; -} - -Real InteractionStatistics::GetDistanceBucketSize() const -{ - return bucket_size_; -} - -uint64_t InteractionStatistics::GetCount(Real distance) const -{ - uint64_t total=0; - for (uint64_t i=0; i<atom::UNKNOWN; ++i) { - for (uint64_t j=0; j<atom::UNKNOWN; ++j) { - total+=histogram_.Get(i, j, distance); - } - } - return total; -} - -uint64_t InteractionStatistics::GetCount(atom::ChemType a, - atom::ChemType b, - Real distance) const -{ - return histogram_.Get(a, b, distance); -} - - -uint64_t InteractionStatistics::GetCount(atom::ChemType a, - atom::ChemType b, - int distance_bin) const -{ - typedef InteractionHistogram::IndexType Index; - - return histogram_.Get(Index(a, b, distance_bin)); -} - - - -uint64_t InteractionStatistics::GetCount(atom::ChemType a, - atom::ChemType b) const -{ - typedef InteractionHistogram::IndexType Index; - uint64_t num=uint64_t((upper_cutoff_-lower_cutoff_)/bucket_size_); - uint64_t total=0; - for (uint64_t k=0; k<num; ++k) { - total+=histogram_.Get(Index(a, b, k)); - } - return total; -} - -void InteractionStatistics::Set(atom::ChemType a, atom::ChemType b, - int distance_bin, int counts) { - typedef InteractionHistogram::IndexType Index; - histogram_.Set(Index(a, b, distance_bin), counts); -} - - -template <typename DS> -void InteractionStatistics::Serialize(DS& ds) -{ - ds & lower_cutoff_; - ds & upper_cutoff_; - ds & bucket_size_; - ds & sequence_sep_; - if (ds.IsSource()) { - upper_sqr_=upper_cutoff_*upper_cutoff_; - lower_sqr_=lower_cutoff_*lower_cutoff_; - } - ds & histogram_; -} - -InteractionStatisticsPtr InteractionStatistics::LoadFromFile(const String& fn) -{ - std::ifstream stream(fn.c_str(), std::ios_base::binary); - io::BinaryDataSource ds(stream); - InteractionStatisticsPtr p(new InteractionStatistics); - ds >> *p.get(); - return p; -} - -}} diff --git a/modules/qa/src/interaction_statistics.hh b/modules/qa/src/interaction_statistics.hh deleted file mode 100644 index 046a24cc4fe02f5df5c14e0fb4b9910de1138f2b..0000000000000000000000000000000000000000 --- a/modules/qa/src/interaction_statistics.hh +++ /dev/null @@ -1,122 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#ifndef OST_QA_INTERACTION_STATISTICS_HH -#define OST_QA_INTERACTION_STATISTICS_HH -/* - Author: Marco Biasini - */ -#include <ost/mol/entity_view.hh> -#include <ost/mol/entity_visitor.hh> -#include <ost/qa/module_config.hh> -#include <ost/qa/amino_acids.hh> -#include <ost/qa/histogram.hh> -#include <ost/qa/atom_types.hh> - - - -namespace ost { namespace qa { - -/// \brief histogram for interaction histogram -/// -/// the meaning of the dimensions is: type of first interaction partner, type -/// of second interaction partner, distance between the two paterns. -typedef Histogram<int, int, Real> InteractionHistogram; -class InteractionStatistics; -typedef boost::shared_ptr<InteractionStatistics> InteractionStatisticsPtr; - -class DLLEXPORT_OST_QA InteractionStatistics : public mol::EntityVisitor { -public: - /// \brief construct new interaction statistics. - /// - /// \param lower_cutoff only interaction parters that are further apart than - /// this cutoff are taken into account. - /// \param upper_cutoff only interaction partners below the cutoff are taken - /// into account - /// \param bucket_size size of distance bins in Angstroem. - /// \param sequence_sep residues that are closer than sequence_sep in - /// sequence are not considered. - InteractionStatistics(Real lower_cutoff, - Real upper_cutoff, - Real bucket_size, - int sequence_sep); - - /// \brief load interaction statistics from file - static InteractionStatisticsPtr LoadFromFile(const String& file_name); - - - /// Extract the interaction potential from the given entity view. - void Extract(mol::EntityView a, mol::EntityView b); - - void SaveToFile(const String& file_name) const; - - - /// \brief Get distance bucket size - Real GetDistanceBucketSize() const; - - int GetSequenceSeparation() const; - - Real GetUpperCutoff() const; - - Real GetLowerCutoff() const; - - /// \brief Get number of atoms at the given distance regardless of their atom - /// type. - uint64_t GetCount(Real distance) const; - - /// \brief get number of atoms at given distance by taking their type into - /// account. - uint64_t GetCount(atom::ChemType a, atom::ChemType b, Real distance) const; - - /// \brief get number of atoms in given distance-bin by taking their type into - /// account. - uint64_t GetCount(atom::ChemType a, atom::ChemType b, int distance_bin) const; - - /// \brief get number of atoms at given distance by taking their type into - /// account. - uint64_t GetCount(atom::ChemType a, atom::ChemType b) const; - - /// \brief get number of atoms at given distance by taking their type into - /// account. - void Set(atom::ChemType a, atom::ChemType b, int distance_bin, int counts); - - /// \internal - template <typename DS> - void Serialize(DS& ds); -public: - virtual bool VisitResidue(const mol::ResidueHandle& r); - virtual bool VisitAtom(const mol::AtomHandle& a); - -private: - InteractionStatistics(); - - Real lower_cutoff_; - Real lower_sqr_; - Real upper_cutoff_; - Real upper_sqr_; - - Real bucket_size_; - int sequence_sep_; - mol::EntityView view_a_; - mol::EntityView view_b_; - AminoAcid amino_acid_; - InteractionHistogram histogram_; -}; - -}} -#endif diff --git a/modules/qa/src/module_config.hh b/modules/qa/src/module_config.hh deleted file mode 100644 index 1d07f3dd7c95819cb2e4c4315c2eae5302e59722..0000000000000000000000000000000000000000 --- a/modules/qa/src/module_config.hh +++ /dev/null @@ -1,30 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#ifndef OST_QA_MODULE_CONFIG_HH -#define OST_QA_MODULE_CONFIG_HH - -#include <ost/qa/module_config.hh> - -#if defined(OST_MODULE_OST_QA) -# define DLLEXPORT_OST_QA DLLEXPORT -#else -# define DLLEXPORT_OST_QA DLLIMPORT -#endif - -#endif diff --git a/modules/qa/src/multi_classifier.hh b/modules/qa/src/multi_classifier.hh deleted file mode 100644 index 846fdddecbadd99b92093482bcd1024923d6797a..0000000000000000000000000000000000000000 --- a/modules/qa/src/multi_classifier.hh +++ /dev/null @@ -1,455 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#ifndef OST_QA_MULTI_CLASSIFIER_HH -#define OST_QA_MULTI_CLASSIFIER_HH - -#include <ost/stdint.hh> -#include <vector> -#include <cassert> -#include <fstream> -#include <ost/message.hh> -#include <iostream> - -#include <boost/shared_ptr.hpp> - -#include "index.hh" -#include <ost/config.hh> - -namespace ost { namespace qa { - -namespace impl { - // Determine length of a typelist, i.e. the number of template arguments - // different from NullType. The length can be accessed over the Value - // member. - // usage example: - // int length=LengthOf<int, int>::Value; - struct NullType { - template <typename DS> - void Serialize(DS&) {} - }; - template <typename T1, - typename T2, - typename T3, - typename T4, - typename T5, - typename T6, - typename T7> - struct LengthOf; - - template <typename T1, - typename T2, - typename T3, - typename T4, - typename T5, - typename T6, - typename T7> - struct LengthOf { - enum { Value = 1+LengthOf<T2, T3, T4, T5, T6, T7, NullType>::Value }; - }; - template <> - struct LengthOf<NullType, NullType, NullType, NullType, NullType, NullType, - NullType> - { - enum { Value = 0 }; - }; - - //! compile time if statement. If<C, T, F>::Type is set to T if C is true - // and set to F otherwise. - template <bool C, - typename T, - typename F> - struct If; - - template <typename T, - typename F> - struct If<true, T, F> { - typedef T Type; - }; - - template <typename T, - typename F> - struct If<false, T, F> { - typedef F Type; - }; - - //! Compile-time check for identity of two types. IsEqual<T1,T2>::Value - // is true if T1 and T2 are identical, and false otherwise. - template <typename T1, - typename T2> - struct IsEqual; - template <typename T1, - typename T2> - struct IsEqual { - enum { Value = false }; - }; - template <typename T> - struct IsEqual<T,T> { - enum { Value = true }; - }; - - // Compile-time selector statement to execute conditional statement based - // on the identity of \c C. If C is identical to the NullType - // then IfNull<...>::Type is a typedef of T, otherwise, IfNull<...>::Type is a - // typedef for F. - template <typename C, typename T, typename F> - struct IfNull; - - template <typename C, - typename T, - typename F> - struct IfNull { - typedef typename If<IsEqual<NullType, C>::Value, T, F>::Type Type; - }; - -} - -//! Base class for classifiers. -class DLLEXPORT_OST_QA ClassifierBase { -public: - ClassifierBase(uint32_t number_of_classes) - : number_of_classes_(number_of_classes) { - } - ClassifierBase() - : number_of_classes_(0) {} - virtual ~ClassifierBase() {} - uint32_t GetNumberOfClasses() const { - return number_of_classes_; - } -protected: - uint32_t number_of_classes_; -}; - -//! Classifier for integral classes. -class DLLEXPORT_OST_QA IntegralClassifier : public ClassifierBase { -public: - IntegralClassifier(uint32_t number_of_classes, - int lower_bound) - : ClassifierBase(number_of_classes), - lower_bound_(lower_bound) { - } - uint32_t GetIndexOf(int value) const { - uint32_t idx=(value-lower_bound_); - assert(this->GetNumberOfClasses()>idx); - return idx; - } - IntegralClassifier() - : ClassifierBase(0), - lower_bound_(0) { - } - - template <typename DS> - void Serialize(DS& ds) - { - ds & number_of_classes_; - ds & lower_bound_; - } -private: - int32_t lower_bound_; -}; - -//! Classifier for real valued classes. -class DLLEXPORT_OST_QA ContinuousClassifier : public ClassifierBase { -public: - ContinuousClassifier(uint32_t number_of_classes, - Real lower_bound, - Real upper_bound) - : ClassifierBase(number_of_classes), - lower_bound_(lower_bound), - upper_bound_(upper_bound) { - } - uint32_t GetIndexOf(Real value) const { - Real factor=(value-lower_bound_)/(upper_bound_-lower_bound_); - uint32_t idx=uint32_t(floor(this->GetNumberOfClasses()*factor)); -// std::cout << value << " " << factor << std::endl; - assert(this->GetNumberOfClasses()>idx); - return idx; - } - ContinuousClassifier() - : ClassifierBase(1), - lower_bound_(0), upper_bound_(1) { - } - template <typename DS> - void Serialize(DS& ds) - { - ds & number_of_classes_; - ds & lower_bound_; - ds & upper_bound_; - } -private: - Real lower_bound_; - Real upper_bound_; -}; - - -template <typename T> -struct Classifier; - -template <> -struct DLLEXPORT_OST_QA Classifier<int> { - typedef IntegralClassifier Type; - typedef const IntegralClassifier& ConstRefType; - typedef IntegralClassifier& RefType; -}; -template <> -struct DLLEXPORT_OST_QA Classifier<Real> { - typedef ContinuousClassifier Type; - typedef const ContinuousClassifier& ConstRefType; - typedef ContinuousClassifier& RefType; -}; -#if OST_DOUBLE_PRECISION -template <> -struct DLLEXPORT_OST_QA Classifier<float> { - typedef ContinuousClassifier Type; - typedef const ContinuousClassifier& ConstRefType; - typedef ContinuousClassifier& RefType; -}; -#endif -template <> -struct DLLEXPORT_OST_QA Classifier<impl::NullType> { - typedef impl::NullType Type; - typedef const impl::NullType& ConstRefType; - typedef impl::NullType& RefType; -}; - -template <typename I> -struct DLLEXPORT_OST_QA NullFind { - NullFind(const ClassifierBase&,uint32_t,const impl::NullType&,I&) {}; -}; -template <typename C, typename T, typename I> -struct IndexFind; - -template <typename C, typename I> -struct DLLEXPORT_OST_QA IndexFind<C,impl::NullType,I> { - IndexFind(const C&, - uint32_t, - const impl::NullType&, I&) { - } -}; - -template <typename C, typename T, typename I> -struct DLLEXPORT_OST_QA IndexFind { - IndexFind(const C& classifier, uint32_t i, const T& value, I& index) { - index[i]=classifier.GetIndexOf(value); - } -}; -template <typename T> -struct NumberOfClasses; - - -template <> -struct DLLEXPORT_OST_QA NumberOfClasses<impl::NullType> { - uint32_t operator ()(const impl::NullType& t) { - return 1; - } -}; - -template <typename T> -struct DLLEXPORT_OST_QA NumberOfClasses { - uint32_t operator ()(const T& t) { - return t.GetNumberOfClasses(); - } -}; - -//! \brief generic n-dimensional classifier -template <typename V, typename T1, - typename T2=impl::NullType, - typename T3=impl::NullType, - typename T4=impl::NullType, - typename T5=impl::NullType, - typename T6=impl::NullType, - typename T7=impl::NullType> -class DLLEXPORT_OST_QA MultiClassifier { -public: - enum { Dimension = impl::LengthOf<T1, T2, T3, T4, T5, T6, T7>::Value }; - typedef V ValueType; - typedef Index<MultiClassifier::Dimension> IndexType; - typedef IndexIterator<Dimension> Iterator; - typedef Classifier<T1> C1; - typedef Classifier<T2> C2; - typedef Classifier<T3> C3; - typedef Classifier<T4> C4; - typedef Classifier<T5> C5; - typedef Classifier<T6> C6; - typedef Classifier<T7> C7; -#if WIN32 - MultiClassifier(const V& initial_value, - typename C1::ConstRefType c1, - typename C2::ConstRefType c2=C2::Type(), - typename C3::ConstRefType c3=C3::Type(), - typename C4::ConstRefType c4=C4::Type(), - typename C5::ConstRefType c5=C5::Type(), - typename C6::ConstRefType c6=C6::Type(), - typename C7::ConstRefType c7=C7::Type()) -#else - MultiClassifier(const V& initial_value, - typename C1::ConstRefType c1, - typename C2::ConstRefType c2=typename C2::Type(), - typename C3::ConstRefType c3=typename C3::Type(), - typename C4::ConstRefType c4=typename C4::Type(), - typename C5::ConstRefType c5=typename C5::Type(), - typename C6::ConstRefType c6=typename C6::Type(), - typename C7::ConstRefType c7=typename C7::Type()) -#endif - : classifier1_(c1), classifier2_(c2), classifier3_(c3), - classifier4_(c4), classifier5_(c5), classifier6_(c6), - classifier7_(c7) { - this->ExtractNumberOfClasses(); - // allocate enough memory for all the buckets - uint32_t total=this->CalculateNumberOfBuckets(); - buckets_.resize(total, initial_value); - } - - MultiClassifier() - { - memset(number_of_classes_, 0, sizeof(number_of_classes_)); - } - - template <typename DS> - void Serialize(DS& ds) - { - ds & classifier1_; - ds & classifier2_; - ds & classifier3_; - ds & classifier4_; - ds & classifier5_; - ds & classifier6_; - ds & classifier7_; - if (ds.IsSource()) { - this->ExtractNumberOfClasses(); - } - ds & buckets_; - } - - MultiClassifier(const MultiClassifier& rhs) - : classifier1_(rhs.classifier1_), classifier2_(rhs.classifier2_), - classifier3_(rhs.classifier3_), classifier4_(rhs.classifier4_), - classifier5_(rhs.classifier5_), classifier6_(rhs.classifier6_), - classifier7_(rhs.classifier7_) { - this->ExtractNumberOfClasses(); - uint32_t total=this->CalculateNumberOfBuckets(); - buckets_.resize(total); - memcpy(&buckets_.front(), &rhs.buckets_.front(), sizeof(V)*total); - } - - uint32_t GetBucketCount() const { - return static_cast<uint32_t>(buckets_.size()); - } - - void Add(const ValueType& value, - T1 x1=T1(), T2 x2=T2(), - T3 x3=T3(), T4 x4=T4(), - T5 x5=T5(), T6 x6=T6(), - T7 x7=T7()) { - IndexType index=this->FindBucket(x1, x2, x3, x4, x5, x6, x7); - uint32_t linear_index=this->LinearizeBucketIndex(index); - buckets_[linear_index]+=value; - } - - const ValueType& Get(T1 x1=T1(), T2 x2=T2(), - T3 x3=T3(), T4 x4=T4(), - T5 x5=T5(), T6 x6=T6(), T7 x7=T7()) const { - IndexType index=this->FindBucket(x1, x2, x3, x4, x5, x6, x7); - uint32_t linear_index=this->LinearizeBucketIndex(index); - return buckets_[linear_index]; - } - - const ValueType& Get(const IndexType& index) const - { - return buckets_[this->LinearizeBucketIndex(index)]; - } - - void Set(const IndexType& index, const ValueType& value) - { - buckets_[this->LinearizeBucketIndex(index)]=value; - } - //TODO Make sure that FindBucket is called with the correct number of - // arguments. Can be done over a wrapper type around T1..T6 that only - // provides a default c'tor when T is equal to NullType. - IndexType FindBucket(T1 x1=T1(), T2 x2=T2(), T3 x3=T3(), - T4 x4=T4(), T5 x5=T5(), T6 x6=T6(), - T7 x7=T7()) const { - // determine indices for parameters whose type is not equal to - // NullType - IndexType index; - IndexFind<typename C1::Type, T1, - IndexType> find_index_1(classifier1_, 0, x1, index); - IndexFind<typename C2::Type, T2, - IndexType> find_index_2(classifier2_, 1, x2, index); - IndexFind<typename C3::Type, T3, - IndexType> find_index_3(classifier3_, 2, x3, index); - IndexFind<typename C4::Type, T4, - IndexType> find_index_4(classifier4_, 3, x4, index); - IndexFind<typename C5::Type, T5, - IndexType> find_index_5(classifier5_, 4, x5, index); - IndexFind<typename C6::Type, T6, - IndexType> find_index_6(classifier6_, 5, x6, index); - IndexFind<typename C7::Type, T7, - IndexType> find_index_7(classifier7_, 6, x7, index); - return index; - } - - void Add(const ValueType& value, const IndexType& index) - { - buckets_[this->LinearizeBucketIndex(index)]+=value; - } -private: - void ExtractNumberOfClasses() - { - number_of_classes_[0]=NumberOfClasses<typename C1::Type>()(classifier1_); - number_of_classes_[1]=NumberOfClasses<typename C2::Type>()(classifier2_); - number_of_classes_[2]=NumberOfClasses<typename C3::Type>()(classifier3_); - number_of_classes_[3]=NumberOfClasses<typename C4::Type>()(classifier4_); - number_of_classes_[4]=NumberOfClasses<typename C5::Type>()(classifier5_); - number_of_classes_[5]=NumberOfClasses<typename C6::Type>()(classifier6_); - number_of_classes_[6]=NumberOfClasses<typename C7::Type>()(classifier7_); - } - - uint32_t LinearizeBucketIndex(const IndexType& index) const - { - uint32_t factor=1; - uint32_t linear_index=0; - for (uint32_t i=0; i<MultiClassifier::Dimension; ++i) { - linear_index+=factor*index[i]; - factor*=number_of_classes_[i]; - } - return linear_index; - } - - uint32_t CalculateNumberOfBuckets() const - { - uint32_t total=1; - for (uint32_t i=0; i<MultiClassifier::Dimension; ++i) { - total*=number_of_classes_[i]; - } - return total; - } - typename C1::Type classifier1_; - typename C2::Type classifier2_; - typename C3::Type classifier3_; - typename C4::Type classifier4_; - typename C5::Type classifier5_; - typename C6::Type classifier6_; - typename C7::Type classifier7_; - uint32_t number_of_classes_[7]; - std::vector<ValueType> buckets_; -}; - -}} - -#endif diff --git a/modules/qa/src/packing_potential.cc b/modules/qa/src/packing_potential.cc deleted file mode 100644 index b64ab96d0c05742cc6736be5fca87edf882962fb..0000000000000000000000000000000000000000 --- a/modules/qa/src/packing_potential.cc +++ /dev/null @@ -1,200 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#include "packing_potential.hh" -#include <ost/io/binary_data_source.hh> -#include <ost/io/binary_data_sink.hh> -#include <ost/io/container_serialization.hh> -#include <boost/filesystem/convenience.hpp> -#include <ost/message.hh> -#include <sstream> -#include <ost/mol/atom_view.hh> -#include <ost/io/io_exception.hh> -#include <ost/integrity_error.hh> - -/* - Authors: Marco Biasini, Pascal Benkert - */ - -namespace ost { namespace qa { - - -PackingPotentialOpts::PackingPotentialOpts(): - sigma(0.02), cutoff(10.0), bucket_size(1), max_counts(25) -{ - -} - - -PackingPotentialOpts::PackingPotentialOpts(int b, int m, Real c, Real s): - sigma(s), cutoff(c), bucket_size(b), max_counts(m) -{ - -} - -template <typename DS> -void PackingPotentialOpts::Serialize(DS& ds) -{ - ds & sigma; - ds & cutoff; - ds & bucket_size; - ds & max_counts; -} - - -PackingPotential::PackingPotential() -{ - -} - -void PackingPotential::Fill(const PackingStatisticsPtr& stat, bool calculate_average_energy_flag) -{ - typedef PackingEnergies::IndexType Index; - int buckets=(options_.max_counts/options_.bucket_size)+1; // +1 since zero counts are in the first bin - uint64_t M=stat->GetTotalCount(); - for (int i=0; i < Xxx; ++i) { - uint64_t Mi=stat->GetCount(AminoAcid(i)); - - if (calculate_average_energy_flag==1) std::cout << AminoAcid(i) << " " << float(Mi)/M << " "; - Real per_aa_energy = 0; // for average energy - - for (int j=0; j<buckets; ++j) { - uint64_t fxi=stat->GetCount(AminoAcid(i), j*options_.bucket_size); - // could be moved out of the loop and cached for more speed. at the - // moment I just don't care... - uint64_t fx=stat->GetCount(j*options_.bucket_size); - - // propensity = (fxi/Mi)/(fx/M) - Real propensity=0.0; - //avoid division by zero: - if (fx != 0 && Mi != 0 && M != 0) { - propensity=Real(Real(fxi)/Mi)/(Real(fx)/M); - } - - Real e=log(1+options_.sigma*Mi)- - log(1+options_.sigma*Mi*propensity); - - per_aa_energy=per_aa_energy+fxi*e; // for average energy - - energies_.Set(Index(i, j), e); - } - if (calculate_average_energy_flag==1) std::cout << per_aa_energy/Mi << std::endl; - } -} - -PackingPotentialPtr PackingPotential::Create(const PackingStatisticsPtr& stat, - const PackingPotentialOpts& opts, - bool calculate_average_energy_flag) -{ -// stat->bucket_size_=opts.bucket_size; -// stat->max_counts_=opts.max_counts; -// stat->cutoff_=opts.cutoff; - - if (stat->GetBucketSize() != opts.bucket_size) { - std::stringstream error_message; - error_message << "BucketSize limit specified in opts does not agree with those of the statistics file: " << opts.bucket_size << " vs " << stat->GetBucketSize(); - throw IntegrityError(error_message.str()); - } - - if (stat->GetMaxCounts() != opts.max_counts) { - std::stringstream error_message; - error_message << "MaxCounts limit specified in opts does not agree with those of the statistics file: " << opts.max_counts << " vs " << stat->GetMaxCounts(); - throw IntegrityError(error_message.str()); - } - - if (stat->GetCutoffDistance() != opts.cutoff) { - std::stringstream error_message; - error_message << "CutoffDistance limit specified in opts does not agree with those of the statistics file: " << opts.cutoff << " vs " << stat->GetCutoffDistance(); - throw IntegrityError(error_message.str()); - } - - PackingPotentialPtr p(new PackingPotential); - p->options_=opts; - - p->options_.cutoff=opts.cutoff; - p->options_.bucket_size=opts.bucket_size; - p->options_.max_counts=opts.max_counts; - IntegralClassifier counts((opts.max_counts/opts.bucket_size)+1, 0); // +1 since zero counts are in the first bin - p->energies_=PackingEnergies(0.0, IntegralClassifier(Xxx, 0), counts); - p->Fill(stat, calculate_average_energy_flag); - return p; -} - -PackingPotentialPtr PackingPotential::LoadFromFile(const String& filename) -{ - if(!boost::filesystem::exists(filename)) - throw io::IOException("Could not open packing potential data file.\nFile does not exist at: "+filename); - - PackingPotentialPtr p(new PackingPotential); - std::ifstream stream(filename.c_str(), std::ios_base::binary); - io::BinaryDataSource ds(stream); - ds >> *p.get(); - return p; -} - -void PackingPotential::SaveToFile(const String& filename) -{ - std::ofstream stream(filename.c_str(), std::ios_base::binary); - io::BinaryDataSink ds(stream); - ds << *this; -} - -template <typename DS> -void PackingPotential::Serialize(DS& ds) -{ - ds & options_; - ds & energies_; -} - - -Real PackingPotential::GetTotalEnergy(const mol::EntityView& view, - const mol::EntityViewList& views) -{ - view_=view; - views_=views; - energy_=0.0; - energy_counts_=0; - view_.Apply(*this); - return energy_; -} - -int PackingPotential::GetEnergyCounts() -{ - return energy_counts_; -} - -bool PackingPotential::VisitAtom(const mol::AtomHandle& atom) -{ - AminoAcid aa=ResidueToAminoAcid(atom.GetResidue()); - if (aa==Xxx) { - atom.GetResidue().SetFloatProp("solvation_energy",0); - return false; - } - int count=0; - for (mol::EntityViewList::iterator i=views_.begin(), - e=views_.end(); i!=e; ++i) { - count+=i->FindWithin(atom.GetPos(), options_.cutoff).size(); - } - float local_energy=this->GetPackingEnergy(aa, count); - atom.GetResidue().SetFloatProp("solvation_energy",local_energy); - energy_+=local_energy; - energy_counts_++; - return false; -} - -}} diff --git a/modules/qa/src/packing_potential.hh b/modules/qa/src/packing_potential.hh deleted file mode 100644 index 00c3229f6ea55b5e39d6c2b495baee7796bbe87b..0000000000000000000000000000000000000000 --- a/modules/qa/src/packing_potential.hh +++ /dev/null @@ -1,126 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#ifndef OST_QA_PACKING_POTENTIAL_HH -#define OST_QA_PACKING_POTENTIAL_HH - -/* - Author: Marco Biasini - */ - - -#include <boost/shared_ptr.hpp> - -#include <ost/mol/entity_visitor.hh> - -#include "packing_statistics.hh" -#include "multi_classifier.hh" - -namespace ost { namespace qa { - -struct DLLEXPORT_OST_QA PackingPotentialOpts { - /// \brief initialise with sigma 0.02 and bucket_size 1 distance cutoff - /// 2.0 and max_counts 100. - PackingPotentialOpts(); - - /// \brief initialise with specified bucket size, max count and sigma - PackingPotentialOpts(int bucket_size, int max_count, Real cutoff, - Real sigma=0.02); - /// \brief sigma for smoothing, after Sippl, 1990. Usually set to 0.02 - Real sigma; - /// \brief upper distance cutoff - Real cutoff; - /// \brief bin size for counts. - int bucket_size; - /// \brief maximal number of atoms - int max_counts; - - template <typename DS> - void Serialize(DS& ds); -}; - -class PackingPotential; -typedef boost::shared_ptr<PackingPotential> PackingPotentialPtr; - - -/// \brief packing (solvent accessibility) potential -/// -/// Not thread safe -class DLLEXPORT_OST_QA PackingPotential : public mol::EntityVisitor { -public: - typedef MultiClassifier<float, int, int> PackingEnergies; -public: - PackingPotential(); - - /// \brief create new solvation potential with the given solvation statistics - /// and options - /// If the flag ist set the average energy per is calculated - /// and shown on stdout. May be used to correct compositional bias in - /// scoring functions (e.g. Pro and Gly have on average lower torsion energies). - static PackingPotentialPtr Create(const PackingStatisticsPtr& stat, - const PackingPotentialOpts& opts, - bool calculate_average_energy_flag=0); - - static PackingPotentialPtr LoadFromFile(const String& filename); - - void SaveToFile(const String& filename); - - /// \brief extract packing statistics of entity view. - /// - /// See documentation for PackingStatistics::Extract() - Real GetTotalEnergy(const mol::EntityView& view, - const mol::EntityViewList& views); - - int GetEnergyCounts(); - - - Real GetPackingEnergy(AminoAcid aa, int count) const - { - count=count>options_.max_counts ? options_.max_counts : count; - return energies_.Get(aa, count/options_.bucket_size); - } - const PackingPotentialOpts& GetOptions() const { return options_; } - template <typename DS> - void Serialize(DS& ds); - -public: - /// \internal - virtual bool VisitAtom(const mol::AtomHandle& atom); -private: - /// \brief potenials are calculated based on the derived statistics - /// - /// If the flag ist set the average energy per is calculated - /// and shown on stdout. May be used to correct compositional bias in - /// scoring functions (e.g. Pro and Gly have on average lower torsion energies). - void Fill(const PackingStatisticsPtr& stat, - bool calculate_average_energy_flag=0); - - PackingPotentialOpts options_; - PackingEnergies energies_; - - // used to calculate total energy... - mol::EntityView view_; - mol::EntityViewList views_; - Real energy_; - int energy_counts_; -}; - -}} - -#endif - diff --git a/modules/qa/src/packing_statistics.cc b/modules/qa/src/packing_statistics.cc deleted file mode 100644 index d557931ca9ec23af1b7cd81a5cf792a71230eabc..0000000000000000000000000000000000000000 --- a/modules/qa/src/packing_statistics.cc +++ /dev/null @@ -1,143 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#include "packing_statistics.hh" - -#include <ost/log.hh> -#include <ost/geom/geom.hh> -#include <ost/mol/atom_handle.hh> -#include <ost/mol/residue_handle.hh> -#include <ost/mol/atom_view.hh> -#include <ost/io/binary_data_source.hh> -#include <ost/io/binary_data_sink.hh> -#include <ost/io/container_serialization.hh> - -namespace ost { namespace qa { - -PackingStatistics::PackingStatistics(Real cutoff, - int max_counts, - int bucket_size) - : cutoff_(cutoff), - max_counts_(max_counts), - bucket_size_(bucket_size), - histogram_(IntegralClassifier(Xxx, 0), // Xxx is the last object in enum and therefore corresponds to the total number of atom types - IntegralClassifier((max_counts/bucket_size)+1, 0)) { // since zero counts is possible, we have n+1 bins - if (max_counts%bucket_size!=0) { - LOG_INFO("initialising packing statistics with max_count=" - << max_counts << "which is not divisible by bucket_size=" - << bucket_size << "without remainder. Will continue anyway."); - } -} - - -PackingStatistics::PackingStatistics(): - cutoff_(0.0), - max_counts_(0), - bucket_size_(0), - histogram_() -{ -} - - -void PackingStatistics::Extract(mol::EntityView view, const mol::EntityViewList& views) -{ - view_=view; - views_=views; - view.Apply(*this); -} - -bool PackingStatistics::VisitAtom(const mol::AtomHandle& atom) { - AminoAcid aa=ResidueToAminoAcid(atom.GetResidue()); - if (aa==Xxx) - return false; - int count=0; - for (mol::EntityViewList::iterator i=views_.begin(), - e=views_.end(); i!=e; ++i) { - count+=i->FindWithin(atom.GetPos(), cutoff_).size(); - } - count=count>=max_counts_ ? max_counts_ : count; // if more counts observed than maximum expected -> set to maximum - histogram_.Add(1, aa, count/bucket_size_); - return false; -} - -int PackingStatistics::GetBucketSize() const { - return bucket_size_; -} - -int PackingStatistics::GetMaxCounts() const { - return max_counts_; -} - -//! \brief Get cutoff radius of sphere around CB. -Real PackingStatistics::GetCutoffDistance() const { - return cutoff_; -} - -uint64_t PackingStatistics::GetCount(AminoAcid aa, int counts) const { - return histogram_.Get(aa, counts/bucket_size_); -} - -uint64_t PackingStatistics::GetCount(AminoAcid aa) const { - uint64_t b=(max_counts_/bucket_size_)+1; - uint64_t total=0; - for (uint64_t i=0; i<b;++i) { - total+=histogram_.Get(aa, i); - } - return total; -} - -PackingStatisticsPtr PackingStatistics::LoadFromFile(const String& filename) -{ - std::ifstream stream(filename.c_str(), std::ios_base::binary); - io::BinaryDataSource ds(stream); - PackingStatisticsPtr ptr(new PackingStatistics); - ds >> *ptr.get(); - return ptr; -} - -void PackingStatistics::SaveToFile(const String& filename) -{ - std::ofstream stream(filename.c_str(), std::ios_base::binary); - io::BinaryDataSink ds(stream); - ds << *this; -} - -uint64_t PackingStatistics::GetTotalCount() const -{ - uint64_t b=(max_counts_/bucket_size_)+1; - uint64_t total=0; - for (uint64_t i=0; i<Xxx; ++i) { - for (uint64_t j=0; j<b;++j) { - total+=histogram_.Get(i, j); - } - } - return total; -} - -uint64_t PackingStatistics::GetCount(int count) const -{ - uint64_t total=0; - for (uint64_t i=0; i<Xxx; ++i) { - total+=histogram_.Get(i, count/bucket_size_); - } - return total; -} - - - -}} diff --git a/modules/qa/src/packing_statistics.hh b/modules/qa/src/packing_statistics.hh deleted file mode 100644 index 68e4e616c5d37d35f1b5ab4c9b8798772276861e..0000000000000000000000000000000000000000 --- a/modules/qa/src/packing_statistics.hh +++ /dev/null @@ -1,90 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#ifndef OST_QA_PACKING_STATISTICS_HH -#define OST_QA_PACKING_STATISTICS_HH - -#include <ost/mol/entity_visitor.hh> -#include <ost/qa/histogram.hh> -#include <ost/qa/amino_acids.hh> -#include <ost/mol/entity_view.hh> - -namespace ost { namespace qa { - -/// \brief two dimensional histogram -/// the meaning of the histogram axes is: amino acid type, number of -/// atoms inside of cutoff radius. -typedef Histogram<int, int> PackingHistogram; -class PackingStatistics; -typedef boost::shared_ptr<PackingStatistics> PackingStatisticsPtr; -class DLLEXPORT_OST_QA PackingStatistics : public mol::EntityVisitor { -public: - PackingStatistics(Real cutoff, int max_counts, int bucket_size); - PackingStatistics(); - - /// \brief Get bucket size for counts, i.e. 2 indicates that a residue with - /// 6 and 7 seven CB within the cutoff radius have the same packing. - int GetBucketSize() const; - - /// \brief Get maximal count of CB atoms inside cutoff distance. larger - /// values are truncated. - int GetMaxCounts() const; - - /// \brief Get cutoff radius of sphere around CB. - Real GetCutoffDistance() const; - - /// \brief extract packing statistics of entity view. - /// - /// For each peptide residue in view, the total number of CB (CA for glycine) - /// in views are counted. - void Extract(mol::EntityView view, const mol::EntityViewList& views); - - static PackingStatisticsPtr LoadFromFile(const String& filename); - - void SaveToFile(const String& filename); - - virtual bool VisitAtom(const mol::AtomHandle& atom); - - uint64_t GetCount(AminoAcid aa) const; - uint64_t GetCount(AminoAcid aa, int count) const; - - uint64_t GetTotalCount() const; - uint64_t GetCount(int count) const; - - template <typename DS> - void Serialize(DS& ds) - { - ds & cutoff_; - ds & max_counts_; - ds & bucket_size_; - ds & histogram_; - } - - Real cutoff_; - int max_counts_; - int bucket_size_; - -private: - mol::EntityView view_; - mol::EntityViewList views_; - PackingHistogram histogram_; -}; - -}} - -#endif diff --git a/modules/qa/src/rapdf.cc b/modules/qa/src/rapdf.cc deleted file mode 100644 index d530c27b12d23c991233f324cdf96487d0618e04..0000000000000000000000000000000000000000 --- a/modules/qa/src/rapdf.cc +++ /dev/null @@ -1,95 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#include "rapdf.hh" -#include <ost/mol/mol.hh> -#include <ost/geom/geom.hh> -namespace ost { namespace qa { - - -RAPDF::RAPDF(const String& file) { - energies_.reset(new InteractionEnergies(file)); -} - -RAPDF::RAPDF(InteractionStatisticsPtr statistics, Real sigma) { - typedef InteractionEnergies::IndexType Index; - Real upper=statistics->GetUpperCutoff(); - Real lower=statistics->GetLowerCutoff(); - Real bucket_size=statistics->GetDistanceBucketSize(); - uint32_t distance_buckets=uint32_t((upper-lower)/bucket_size); - ContinuousClassifier distance_classifier(distance_buckets, lower, upper); - energies_.reset(new InteractionEnergies(0.0, - IntegralClassifier(atom::UNKNOWN, 0), - IntegralClassifier(atom::UNKNOWN, 0), - distance_classifier)); - // loop over all distance buckets and calculate energies for all atom pairs. - for (uint32_t i=0; i<distance_buckets;++i) { - Index start(0, 0, i), end(atom::UNKNOWN-1, atom::UNKNOWN-1, i); - Real radius=i*bucket_size+bucket_size/2; - uint32_t total=statistics->GetCount(radius); - for (uint32_t j=0; j<atom::UNKNOWN; ++j) { - for (uint32_t k=0; k<atom::UNKNOWN;++k) { - uint32_t lc=statistics->GetCount(atom::ChemType(j), - atom::ChemType(k), radius); - Real frequency=lc/Real(total); - Real energy=log((1+lc*sigma)/(1+lc*sigma*frequency)); - energies_->Set(Index(j, k, i), energy); - } - } - } -} - -Real RAPDF::GetTotalEnergy(mol::EntityView entity) { - tmp_energy_=0.0; - entity.Apply(*this); - return tmp_energy_; -} - -RAPDF::RAPDF(const RAPDF& rhs) { - energies_.reset(new InteractionEnergies(*rhs.energies_.get())); -} - -bool RAPDF::VisitAtom(const mol::AtomHandle& atom) { - mol::EntityHandle e=atom.GetResidue().GetChain().GetEntity(); - atom::ChemType atype=GetAtomTypeByName(amino_acid_, atom.GetName()); - //FIXME: do not hard-code the distance cutoff here. - mol::AtomHandleList atoms=e.FindWithin(atom.GetPos(), 20.0); - mol::AtomHandleList::iterator i; - for(i=atoms.begin(); i!=atoms.end(); ++i) { - mol::AtomHandle& a=*i; - if (abs(a.GetResidue().GetNumber().GetNum()-rnumber_)<3) - continue; - AminoAcid b=ResidueToAminoAcid(a.GetResidue()); - if (b==Xxx) - continue; - atom::ChemType btype=GetAtomTypeByName(b, a.GetName()); - if (btype==atom::UNKNOWN) - continue; - Real d=Length(a.GetPos()-atom.GetPos()); - tmp_energy_+=energies_->Get(atype, btype, d); - } - return true; -} - -bool RAPDF::VisitResidue(const mol::ResidueHandle& residue) { - amino_acid_=ResidueToAminoAcid(residue); - rnumber_=residue.GetNumber().GetNum(); - return amino_acid_!=Xxx; -} - -}} diff --git a/modules/qa/src/rapdf.hh b/modules/qa/src/rapdf.hh deleted file mode 100644 index ae69b054e33210992f1ca532960a013ef7989565..0000000000000000000000000000000000000000 --- a/modules/qa/src/rapdf.hh +++ /dev/null @@ -1,64 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#ifndef OST_QA_RAPDF_HH -#define OST_QA_RAPDF_HH - -#include <ost/mol/entity_view.hh> -#include <boost/shared_ptr.hpp> -#include <boost/scoped_ptr.hpp> -#include <ost/mol/entity_visitor.hh> -#include <ost/qa/interaction_statistics.hh> - -namespace ost { namespace qa { - -typedef MultiClassifier<Real, int, int, Real> InteractionEnergies; -class RAPDF; - -typedef boost::shared_ptr<RAPDF> RAPDFPtr; - -/// \brief RAPDF after (Samudrala & Moult, 1998) -class DLLEXPORT_OST_QA RAPDF : public EntityVisitor { -public: - /// load RAPDF potential from file. - RAPDF(const String& file); - RAPDF(const RAPDF& rhs); - /// create new interaction potential from interaction statistics - /// \param statistics interaction statistics - /// \param sigma mixing parameter (see Sippl et al., 1990) - RAPDF(InteractionStatisticsPtr statistics, Real sigma=50.0); - - /// calculate all-atom interaction score for entity. all atoms of the - /// entity are taken into account, but only energies from atoms in the - /// view contribute to the energy. - /// not thread-safe (just factor out visitor to make thread-safe) - Real GetTotalEnergy(mol::EntityView entity); - - virtual bool VisitAtom(const mol::AtomHandle& atom); - virtual bool VisitResidue(const mol::ResidueHandle& residue); - -private: - boost::scoped_ptr<InteractionEnergies> energies_; - AminoAcid amino_acid_; - int rnumber_; - Real tmp_energy_; -}; - -}} - -#endif diff --git a/modules/qa/src/reduced_potential.cc b/modules/qa/src/reduced_potential.cc deleted file mode 100644 index 3963fc132f0880ae94f715057af0beda3dfaf2b9..0000000000000000000000000000000000000000 --- a/modules/qa/src/reduced_potential.cc +++ /dev/null @@ -1,153 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ - -#include <ost/io/binary_data_source.hh> -#include <ost/io/binary_data_sink.hh> -#include <ost/io/container_serialization.hh> -#include "reduced_potential.hh" -#include "impl/reduced_impl.hh" -namespace ost { namespace qa { - - -namespace { - -class ReducedEnergiesCalc : public impl::ReducedPotentialImpl { -public: - - ReducedEnergiesCalc(const ReducedStatOptions& opts, ReducedEnergies& energies, - mol::EntityHandle ent, bool norm): - impl::ReducedPotentialImpl(opts, ent), - energies_(energies), energy_(0.0), norm_(norm), count_(0) - { } - ReducedEnergiesCalc(const ReducedStatOptions& opts, ReducedEnergies& energies, - mol::EntityView ent, bool norm): - impl::ReducedPotentialImpl(opts, ent), - energies_(energies), energy_(0.0), norm_(norm), count_(0) - { } - virtual void OnInteraction(AminoAcid aa_one, AminoAcid aa_two, - Real dist, Real angle) - { - //std::cout << aa_one << " " << aa_two << " " << dist << " " << angle << std::endl; - energy_+=energies_.Get(aa_one, aa_two, dist, angle); - count_+=1; - } - - Real GetEnergy() const - { - if (norm_) { - return count_ > 0 ? energy_ /count_ : 0.0; - } - return energy_; - } -private: - ReducedEnergies& energies_; - Real energy_; - bool norm_; - uint64_t count_; -}; - - -} - -bool ReducedPotential::GetCAlphaCBetaPos(const mol::ResidueHandle& res, - geom::Vec3& ca_pos, - geom::Vec3& cb_pos) -{ - return impl::ReducedPotentialImpl::GetCAlphaCBetaPos(res, ca_pos, cb_pos); -} - -ReducedPotentialPtr ReducedPotential::Load(const String& filename) -{ - std::ifstream stream(filename.c_str(), std::ios_base::binary); - io::BinaryDataSource ds(stream); - ReducedPotentialPtr p(new ReducedPotential); - ds >> *p.get(); - return p; -} - - -ReducedPotentialPtr ReducedPotential::Create(const ReducedStatisticsPtr& stat) -{ - ReducedPotentialPtr pot(new ReducedPotential); - pot->opts_=stat->GetOptions(); - pot->energies_=ReducedEnergies(0.0, IntegralClassifier(20, 0), - IntegralClassifier(20, 0), - ContinuousClassifier(pot->opts_.num_dist_bins, - pot->opts_.lower_cutoff, - pot->opts_.upper_cutoff), - ContinuousClassifier(pot->opts_.num_angular_bins, - 0.0, M_PI)); - pot->InitFromStats(stat); - return pot; -} - -void ReducedPotential::InitFromStats(const ReducedStatisticsPtr& stats) -{ - typedef ReducedEnergies::IndexType Index; - uint64_t total_count=stats->GetTotalCount(); - for (size_t i=0; i<Xxx; ++i) { - for (size_t j=0; j<Xxx; ++j) { - uint64_t t2=stats->GetCount(AminoAcid(i), AminoAcid(j)); - for (size_t k=0; k<opts_.num_dist_bins; ++k) { - for (size_t l=0; l<opts_.num_angular_bins; ++l) { - uint64_t t3=stats->GetCount(AminoAcid(i), AminoAcid(j), k, l); - // could be moved outside of the loop and cached... - uint64_t t4=stats->GetCount(k, l); - float f1=t2 > 0 ?(float(t3)/float(t2)) : 0.0; - float f2=total_count > 0 ? (float(t4)/float(total_count)) : 0.0; - float d = f2>1e-10 ? f1/f2: 0.0; - float e=0.582*(log(1+0.02*t2)-log(1+0.02*t2*d)); - energies_.Set(Index(i, j, k, l), e); - } - } - } - } -} - -void ReducedPotential::Save(const String& filename) -{ - std::ofstream stream(filename.c_str(), std::ios_base::binary); - io::BinaryDataSink ds(stream); - ds << *this; -} - - -template <typename DS> -void ReducedPotential::Serialize(DS& ds) -{ - ds & opts_; - ds & energies_; -} - -Real ReducedPotential::GetTotalEnergy(ost::mol::EntityHandle ent, bool norm) -{ - ReducedEnergiesCalc calc(opts_, energies_, ent, norm); - ent.Apply(calc); - return calc.GetEnergy(); -} - -Real ReducedPotential::GetTotalEnergy(ost::mol::EntityView ent, bool norm) -{ - ReducedEnergiesCalc calc(opts_, energies_, ent.GetHandle(), norm); - ent.Apply(calc); - return calc.GetEnergy(); -} - - -}} diff --git a/modules/qa/src/reduced_potential.hh b/modules/qa/src/reduced_potential.hh deleted file mode 100644 index e2d600b7f1f01d58197225ada94ed3ec298f4b22..0000000000000000000000000000000000000000 --- a/modules/qa/src/reduced_potential.hh +++ /dev/null @@ -1,78 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#ifndef OST_QA_REDUCED_POTENTIAL_HH -#define OST_QA_REDUCED_POTENTIAL_HH - -#include <ost/qa/reduced_statistics.hh> - - -namespace ost { namespace qa { - - -class ReducedPotential; - -typedef boost::shared_ptr<ReducedPotential> ReducedPotentialPtr; - -typedef MultiClassifier<float, int, int, float, float> ReducedEnergies; -/// \brief reduced model statistical potential operating on the residue level -class DLLEXPORT_OST_QA ReducedPotential { -public: - - - - static ReducedPotentialPtr Create(const ReducedStatisticsPtr& stats); - static ReducedPotentialPtr Load(const String& filename); - - void Save(const String& filename); - - - template <typename DS> - void Serialize(DS& ds); - - const ReducedStatOptions& GetOptions() const { return opts_; } - - Real GetTotalEnergy(ost::mol::EntityHandle ent, bool normalize=true); - - Real GetTotalEnergy(ost::mol::EntityView ent, bool normalize=true); - static bool GetCAlphaCBetaPos(const mol::ResidueHandle& res, - geom::Vec3& ca_pos, - geom::Vec3& cb_pos); - - Real GetEnergy(AminoAcid aa_one, AminoAcid aa_two, - Real ca_dist, Real angle) const - { - if (aa_one==Xxx || aa_two==Xxx) { return 0.0; } - if (ca_dist<opts_.lower_cutoff || ca_dist>=opts_.upper_cutoff) { - return 0.0; - } - if (angle<0.0 || angle>M_PI) { return 0.0; } - return energies_.Get(aa_one, aa_two, ca_dist, angle); - } -private: - void InitFromStats(const ReducedStatisticsPtr& stats); - ReducedPotential() { } - - ReducedStatOptions opts_; - ReducedEnergies energies_; -}; - - -}} - -#endif diff --git a/modules/qa/src/reduced_statistics.cc b/modules/qa/src/reduced_statistics.cc deleted file mode 100644 index 403153e8b295055ee4b1e8fea606dc8f30b14451..0000000000000000000000000000000000000000 --- a/modules/qa/src/reduced_statistics.cc +++ /dev/null @@ -1,152 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ - -#include <ost/log.hh> -#include <ost/mol/mol.hh> -#include "amino_acids.hh" -#include "reduced_statistics.hh" -#include "impl/reduced_impl.hh" - -#include <ost/io/binary_data_source.hh> -#include <ost/io/binary_data_sink.hh> -#include <ost/io/container_serialization.hh> - -namespace ost { namespace qa { - - -namespace { - -class ReducedStatExtractor : public impl::ReducedPotentialImpl { -public: - - ReducedStatExtractor(const ReducedStatOptions& opts, ReducedHistogram& histo, - mol::EntityHandle ent): - impl::ReducedPotentialImpl(opts, ent), - histo_(histo) - { } - - ReducedStatExtractor(const ReducedStatOptions& opts, ReducedHistogram& histo, - mol::EntityView ent): - impl::ReducedPotentialImpl(opts, ent), - histo_(histo) - { } - - virtual void OnInteraction(AminoAcid aa_one, AminoAcid aa_two, - Real dist, Real angle) - { - histo_.Add(1, aa_one, aa_two, dist, angle); - } -private: - ReducedHistogram& histo_; -}; - - -} - -uint64_t ReducedStatistics::GetTotalCount() const -{ - typedef ReducedHistogram::IndexType Index; - uint64_t count=0; - for (size_t i=0; i<Xxx; ++i) { - for (size_t j=0; j<Xxx; ++j) { - for (size_t k=0; k<opts_.num_dist_bins; ++k) { - for (size_t l=0; l<opts_.num_angular_bins; ++l) { - count+=histo_.Get(Index(i, j, k, l)); - } - } - } - } - return count; -} - -uint64_t ReducedStatistics::GetCount(AminoAcid aa_one, AminoAcid aa_two) const -{ - typedef ReducedHistogram::IndexType Index; - uint64_t count=0; - for (size_t k=0; k<opts_.num_dist_bins; ++k) { - for (size_t l=0; l<opts_.num_angular_bins; ++l) { - count+=histo_.Get(Index(aa_one, aa_two, k, l)); - } - } - return count; -} - -uint64_t ReducedStatistics::GetCount(int dist_bin, int ang_bin) const -{ - typedef ReducedHistogram::IndexType Index; - uint64_t count=0; - for (size_t i=0; i<Xxx; ++i) { - for (size_t j=0; j<Xxx; ++j) { - count+=histo_.Get(Index(i, j, dist_bin, ang_bin)); - } - } - return count; -} - -void ReducedStatistics::Save(const String& filename) -{ - std::ofstream stream(filename.c_str(), std::ios_base::binary); - io::BinaryDataSink ds(stream); - ds << *this; -} - - -ReducedStatisticsPtr ReducedStatistics::Load(const String& filename) -{ - std::ifstream stream(filename.c_str(), std::ios_base::binary); - io::BinaryDataSource ds(stream); - ReducedStatisticsPtr p(new ReducedStatistics); - ds >> *p.get(); - return p; -} - -template <typename DS> -void ReducedStatistics::Serialize(DS& ds) -{ - ds & opts_; - ds & histo_; -} - - -void ReducedStatistics::Extract(mol::EntityHandle ent) -{ - if (ent.GetChainCount()!=1) { - std::stringstream ss; - ss << "Expected exactly one chain, but entity has " - << ent.GetChainCount() << " chains"; - throw std::runtime_error(ss.str()); - } - ReducedStatExtractor extractor(opts_, histo_, ent); - ent.Apply(extractor); -} - -void ReducedStatistics::Extract(mol::EntityView ent) -{ - if (ent.GetChainCount()!=1) { - std::stringstream ss; - ss << "Expected exactly one chain, but entity has " - << ent.GetChainCount() << " chains"; - throw std::runtime_error(ss.str()); - } - ReducedStatExtractor extractor(opts_, histo_, ent.GetHandle()); - ent.Apply(extractor); -} - - -}} \ No newline at end of file diff --git a/modules/qa/src/reduced_statistics.hh b/modules/qa/src/reduced_statistics.hh deleted file mode 100644 index 1f1ab2426abf419e1b2e05a02998cb390320567a..0000000000000000000000000000000000000000 --- a/modules/qa/src/reduced_statistics.hh +++ /dev/null @@ -1,123 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#ifndef OST_QA_REDUCED_STATISTICS_HH -#define OST_QA_REDUCED_STATISTICS_HH - -#include <ost/qa/histogram.hh> -#include <ost/qa/amino_acids.hh> -namespace ost { - -namespace mol { - -class EntityHandle; - -} -namespace qa { - -/// \brief aggregates the option for the reduced model statistics -struct DLLEXPORT_OST_QA ReducedStatOptions { - ReducedStatOptions(): lower_cutoff(0), upper_cutoff(0), - num_angular_bins(0), num_dist_bins(0), - sequence_sep(0) - { } - - ReducedStatOptions(Real l_cutoff, Real u_cutoff, uint n_ang_bins, - uint n_dist_bins, uint ssep): - lower_cutoff(l_cutoff), upper_cutoff(u_cutoff), - num_angular_bins(n_ang_bins), num_dist_bins(n_dist_bins), - sequence_sep(ssep) - { } - Real lower_cutoff; - Real upper_cutoff; - uint num_angular_bins; - uint num_dist_bins; - uint sequence_sep; - - - template <typename DS> - void Serialize(DS& ds) - { - ds & lower_cutoff; - ds & upper_cutoff; - ds & num_angular_bins; - ds & num_dist_bins; - ds & sequence_sep; - } -}; - - -class ReducedStatistics; -typedef boost::shared_ptr<ReducedStatistics> ReducedStatisticsPtr; - - -// parametrized as first amino acid, second amino acid, -// calpha-calpha distance, angle -typedef Histogram<int, int, float, float> ReducedHistogram; - -class DLLEXPORT_OST_QA ReducedStatistics { -public: - ReducedStatistics(Real l_cutoff, Real u_cutoff, uint num_ang_bins, - uint num_dist_bins, uint ssep): - opts_(l_cutoff, u_cutoff, num_ang_bins, num_dist_bins, ssep), - histo_(IntegralClassifier(20, 0), IntegralClassifier(20, 0), - ContinuousClassifier(num_dist_bins, l_cutoff, u_cutoff), - ContinuousClassifier(num_ang_bins, 0.0, M_PI)) - { } - - const ReducedStatOptions& GetOptions() const { return opts_; } - - /// \brief extract the statistics from the given entity (handle) - void Extract(mol::EntityHandle ent); - - /// \brief extract statistics from given entity (view) - void Extract(mol::EntityView ent); - - void Save(const String& filename); - - - static ReducedStatisticsPtr Load(const String& filename); - /// \internal - template <typename DS> - void Serialize(DS& ds); - - uint64_t GetTotalCount() const; - - uint64_t GetCount(AminoAcid aa_one, AminoAcid aa_two) const; - - - uint64_t GetCount(AminoAcid aa_one, AminoAcid aa_two, int dist_bin, - int ang_bin) - { - return histo_.Get(ReducedHistogram::IndexType(aa_one, aa_two, dist_bin, - ang_bin)); - } - - uint64_t GetCount(int dist_bin, int ang_bin) const; -private: - ReducedStatistics(): opts_(), histo_() {} - ReducedStatOptions opts_; - ReducedHistogram histo_; -}; - - - -}} - - -#endif diff --git a/modules/qa/src/solis_torsion_potential.cc b/modules/qa/src/solis_torsion_potential.cc deleted file mode 100644 index 3c463ac00da97ca91994ae6cbd1b1dd7fb901dc0..0000000000000000000000000000000000000000 --- a/modules/qa/src/solis_torsion_potential.cc +++ /dev/null @@ -1,266 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#include "solis_torsion_potential.hh" -#include "amino_acids.hh" - -#include <ost/mol/mol.hh> -#include <ost/log.hh> -#include <ost/io/binary_data_source.hh> -#include <ost/io/binary_data_sink.hh> -#include <ost/io/container_serialization.hh> -#include <ost/integrity_error.hh> -#include <boost/filesystem/convenience.hpp> -#include <ost/message.hh> -#include <ost/io/io_exception.hh> - -namespace ost { namespace qa { - -namespace { - -class SolisTorsionEnergyCalc : public mol::EntityVisitor { -public: - SolisTorsionEnergyCalc(TorsionPotentialSolis::TorsionEnergies& energies, - TorsionPotentialOptsSolis opts): - energies_(energies), prev_(Xxx), center_(Xxx), energy_(0.0), - opts_(opts), num_torsions_(0) - { - } - - virtual bool VisitResidue(const mol::ResidueHandle& rh) - { - AminoAcid c=ResidueToAminoAcid(rh); - - if (prev_==Xxx || center_==Xxx || c==Xxx) { - prev_=center_; - center_=c; - cr_=rh; - return false; - } - mol::TorsionHandle phit=cr_.GetPhiTorsion(); - mol::TorsionHandle psit=cr_.GetPsiTorsion(); - if (!phit || !psit) { - prev_=center_; - center_=c; - cr_=rh; - return false; - } - Real phi=phit.GetAngle()*180/M_PI - 0.000000001; - Real psi=psit.GetAngle()*180/M_PI - 0.000000001; - int ibefore=this->GetOuterAAIndex(prev_); - int iafter=this->GetOuterAAIndex(c); - int icenter=this->GetInnerAAIndex(center_); - - - if (iafter>-1 && icenter >-1 && ibefore>-1) { - energy_+=energies_.Get(ibefore, icenter, iafter, phi, psi); - num_torsions_++; - } - else { - LOG_INFO("Amino acid not found in alphabets..."); - } - prev_=center_; - center_=c; - cr_=rh; - return false; - } - - Real GetEnergy() const { - return energy_; - //return num_torsions_==0?0:energy_/num_torsions_; - } - int GetEnergyCounts() const { - return num_torsions_; - } -private: - int GetOuterAAIndex(AminoAcid aa) - { - return this->GetIndex(aa, opts_.outer_alphabet); - } - - int GetInnerAAIndex(AminoAcid aa) - { - return this->GetIndex(aa, opts_.inner_alphabet); - } - - int GetIndex(AminoAcid aa, AminoAcidAlphabet& alpha) { - AminoAcidAlphabet::iterator i=alpha.begin(); - for (int index=0; i!=alpha.end(); ++i, ++index) { - if ((*i).Contains(aa)) { - return index; - } - - } - return -1; - } -private: - TorsionPotentialSolis::TorsionEnergies& energies_; - AminoAcid prev_; - AminoAcid center_; - mol::ResidueHandle cr_; - Real energy_; - TorsionPotentialOptsSolis opts_; - int num_torsions_; -}; - -} - -template <typename DS> -void TorsionPotentialOptsSolis::Serialize(DS& ds) -{ - ds & angular_bucket_size; - ds & outer_alphabet; - ds & inner_alphabet; - ds & sigma; -} - -template <typename DS> -void TorsionPotentialSolis::Serialize(DS& ds) -{ - ds & options_; - ds & energies_; -} - -TorsionPotentialOptsSolis::TorsionPotentialOptsSolis(): - angular_bucket_size(30), sigma(0.02) -{ - -} - - -TorsionPotentialSolisPtr TorsionPotentialSolis::Create( - TorsionStatisticsSolisPtr statistics, - const TorsionPotentialOptsSolis& opts) -{ - typedef TorsionPotentialSolis::TorsionEnergies TorsionEnergies; - int bucket_size=statistics->GetTorsionBucketSize(); - - if ((opts.angular_bucket_size % bucket_size)!=0) { - throw IntegrityError("Angular bucket size for torsion potential is not an " - "integral multiple of the statistics bucket size"); - } - int buckets=360/opts.angular_bucket_size; - TorsionPotentialSolisPtr ptr(new TorsionPotentialSolis); - - ptr->options_.angular_bucket_size = opts.angular_bucket_size; - ptr->options_.outer_alphabet = opts.outer_alphabet; - ptr->options_.inner_alphabet = opts.inner_alphabet; - - uint32_t outer_s=ptr->options_.outer_alphabet.size(); - uint32_t inner_s=ptr->options_.inner_alphabet.size(); - ptr->energies_=TorsionEnergies(0.0, IntegralClassifier(outer_s, 0), - IntegralClassifier(inner_s, 0), - IntegralClassifier(outer_s, 0), - ContinuousClassifier(buckets,-180,180), - ContinuousClassifier(buckets,-180,180)); - ptr->FillAll(statistics); - return ptr; -} - -TorsionPotentialSolisPtr TorsionPotentialSolis::LoadFromFile(const String& path) -{ - if(!boost::filesystem::exists(path)) - throw io::IOException("Could not open torsion potential data file.\nFile does not exist at: "+path); - - std::ifstream stream(path.c_str(), std::ios_base::binary); - io::BinaryDataSource ds(stream); - TorsionPotentialSolisPtr ptr(new TorsionPotentialSolis); - ds >> *ptr.get(); - return ptr; -} - -void TorsionPotentialSolis::SaveToFile(const String& path) -{ - std::ofstream stream(path.c_str(), std::ios_base::binary); - io::BinaryDataSink ds(stream); - ds << *this; -} - -Real TorsionPotentialSolis::GetTotalEnergy(mol::EntityHandle entity) -{ - SolisTorsionEnergyCalc c(energies_, options_); - entity.Apply(c); - num_torsions_ = c.GetEnergyCounts(); - return c.GetEnergy(); -} - -Real TorsionPotentialSolis::GetTotalEnergy(mol::EntityView entity) -{ - SolisTorsionEnergyCalc c(energies_, options_); - entity.Apply(c); - num_torsions_ = c.GetEnergyCounts(); - return c.GetEnergy(); -} - -int TorsionPotentialSolis::GetEnergyCounts() const { - return num_torsions_; -} - - -void TorsionPotentialSolis::FillAll(const TorsionStatisticsSolisPtr& stat) -{ - uint32_t total=stat->GetCount(); - int num_b=360/(options_.angular_bucket_size); - int half_b=options_.angular_bucket_size/2; - for (int phi=0; phi<num_b; ++phi) { - for (int psi=0; psi<num_b; ++psi) { - int phi_a=-180+options_.angular_bucket_size*phi+half_b; - int psi_a=-180+options_.angular_bucket_size*psi+half_b; - this->FillPhiPsi(stat, phi_a, psi_a, total); - } - } -} - -void TorsionPotentialSolis::FillPhiPsi(const TorsionStatisticsSolisPtr& stat, - int phi, int psi, uint32_t M) -{ - // collect total counts for all amino acids first - uint32_t fx=stat->GetCount(phi, psi); - typedef AminoAcidAlphabet::const_iterator AAA; - for (AAA i=options_.outer_alphabet.begin(), - e1=options_.outer_alphabet.end(); i!=e1; ++i) { - for (AAA j=options_.inner_alphabet.begin(), - e2=options_.inner_alphabet.end(); j!=e2; ++j) { - for (AAA k=options_.outer_alphabet.begin(), - e3=options_.outer_alphabet.end(); k!=e3; ++k) { - - uint32_t fxi=stat->GetCount(*i, *j, *k, phi, psi); - uint32_t Mi=stat->GetCount(*i, *j, *k); - - // propensity = (fxi/Mi)/(fx/M) - Real propensity=0.0; - //avoid division by zero: - if (fx != 0 && Mi != 0 && M != 0) { - propensity=Real(Real(fxi)/Mi)/(Real(fx)/M); - } - - Real energy=log(1+options_.sigma*Mi)- - log(1+options_.sigma*Mi*propensity); - - energies_.Add(energy, - i-options_.outer_alphabet.begin(), - j-options_.inner_alphabet.begin(), - k-options_.outer_alphabet.begin(), phi, psi); - } - } - } -} - - - -}} diff --git a/modules/qa/src/solis_torsion_potential.hh b/modules/qa/src/solis_torsion_potential.hh deleted file mode 100644 index acfb2aed108e7d46f1673b37b68f1be5b9b20314..0000000000000000000000000000000000000000 --- a/modules/qa/src/solis_torsion_potential.hh +++ /dev/null @@ -1,121 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#ifndef OST_QA_SOLIS_TORSION_POTENTIAL_HH -#define OST_QA_SOLIS_TORSION_POTENTIAL_HH -/* - - - Author: Marco Biasini, Pascal Benkert - */ -#include <ost/qa/solis_torsion_statistics.hh> -#include <boost/shared_ptr.hpp> -#include <vector> -#include <boost/scoped_ptr.hpp> - -namespace ost { namespace qa { - -typedef std::vector<AminoAcidSet> AminoAcidAlphabet; - -/// \brief torsion potential options -/// -/// These options may be used to configure the different torsion potentials -/// based on a sliding window of 3 amino acids. -struct DLLEXPORT_OST_QA TorsionPotentialOptsSolis { - /// \brief initialize torsion potential with angular bucket size of 30 and - /// sigma of 0.02 - /// - /// The amino acid alphabets are empty. - TorsionPotentialOptsSolis(); - - /// angular bucket size, 360 % angular_bucket_size must be zero - int angular_bucket_size; - /// alphabet of residues before and after the central amino acids in the - /// window. The amino acids in the set must be mutually exclusive - AminoAcidAlphabet outer_alphabet; - /// alphabet of the central amino acid. The amino acid sets must be - /// mutually exclusive - AminoAcidAlphabet inner_alphabet; - /// weighting factor after (Sippl et al., 1990), usually 0.02 - Real sigma; - /// \internal - template <typename DS> - void Serialize(DS& ds); -}; - -class TorsionPotentialSolis; -typedef boost::shared_ptr<TorsionPotentialSolis> TorsionPotentialSolisPtr; - -/// \brief Torsion potential -/// -/// The torsion potential class is parametrisable by TorsionPotentialOptsSolis. -/// -/// This implementation of the torsion angle potential combines 1 torsion angle -/// pair with up to 3 residues in analogy to the potential described by: -/// Armando D. Solis, S. Rackovsky -/// Optimally informative backbone structural propensities in proteins Proteins, -/// 2002 The ordinary trosion angle potential (1 angle-pair, 1 residue) is -/// covered as well. -class DLLEXPORT_OST_QA TorsionPotentialSolis { -public: - /// \brief create new torsion potential with the given torsion statistics - /// and options - static TorsionPotentialSolisPtr Create(TorsionStatisticsSolisPtr statistics, - const TorsionPotentialOptsSolis& opts); - /// \brief Load potential from file. - static TorsionPotentialSolisPtr LoadFromFile(const String& path); - - /// \brief Calculate energy for entity. - Real GetTotalEnergy(mol::EntityHandle entity); - - /// \brief Calculate energy for entity view. Only residues in the view are - /// considered. - Real GetTotalEnergy(mol::EntityView view); - - /// \brief retrieve total number of energy local (i.e. valid residues) - int GetEnergyCounts() const; - - /// \brief save torsion potential - /// - /// The output is in a non-portable binary form. - void SaveToFile(const String& path); - - /// \internal - template <typename DS> - void Serialize(DS& ds); - - -public: - /// Meaning of parameters: AA before AA center AA after, - /// phi torsion angle, psi torsion angle - typedef MultiClassifier<float, int, int, - int, Real, Real> TorsionEnergies; -private: - void FillAll(const TorsionStatisticsSolisPtr& stat); - - void FillPhiPsi(const TorsionStatisticsSolisPtr& stat, int phi, int psi, - uint32_t total); - TorsionPotentialOptsSolis options_; - TorsionEnergies energies_; - int num_torsions_; -}; - - -}} - -#endif diff --git a/modules/qa/src/solis_torsion_statistics.cc b/modules/qa/src/solis_torsion_statistics.cc deleted file mode 100644 index ed2c6ac1683e15d9bcfaa5d76d261d2425543aec..0000000000000000000000000000000000000000 --- a/modules/qa/src/solis_torsion_statistics.cc +++ /dev/null @@ -1,261 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#include "solis_torsion_statistics.hh" -#include "amino_acids.hh" -#include <ost/log.hh> -#include <ost/mol/mol.hh> -#include <fstream> -#include <ost/message.hh> -#include <ost/io/binary_data_source.hh> -#include <ost/io/binary_data_sink.hh> -#include <ost/io/container_serialization.hh> - -/* - Author: Marco Biasini, Pascal Benkert - */ -namespace ost { namespace qa { - - -namespace { - -class TorsionStatSolis : public mol::EntityVisitor { -public: - TorsionStatSolis(TorsionStatisticsSolis::TorsionHistogram& histo): - histo_(histo) - {} - - virtual bool VisitResidue(const mol::ResidueHandle& residue) { - - mol::ResidueHandle prev=residue.GetPrev(); - mol::ResidueHandle next=residue.GetNext(); - if (!(next && prev && next.IsPeptideLinking() && - prev.IsPeptideLinking() && residue.IsPeptideLinking())) { - return false; - } - - AminoAcid pa=ResidueToAminoAcid(prev); - AminoAcid na=ResidueToAminoAcid(next); - AminoAcid ca=ResidueToAminoAcid(residue); - - if (pa==Xxx || na==Xxx || ca==Xxx) - return false; - mol::TorsionHandle phit=residue.GetPhiTorsion(); - mol::TorsionHandle psit=residue.GetPsiTorsion(); - if (!phit || !psit) - return false; - Real phi=phit.GetAngle()*180/M_PI - 0.000000001; - Real psi=psit.GetAngle()*180/M_PI - 0.000000001; - histo_.Add(1, pa, ca, na, phi, psi); - return false; - } -private: - TorsionStatisticsSolis::TorsionHistogram& histo_; -}; - -} - -// initialise with 0 or 1 (1 better for combined potentials) -TorsionStatisticsSolis::TorsionStatisticsSolis(int torsion_bucket_size) - : histogram_(1,IntegralClassifier(20, 0), - IntegralClassifier(20, 0), - IntegralClassifier(20, 0), - ContinuousClassifier(int(360/torsion_bucket_size), -180, 180), - ContinuousClassifier(int(360/torsion_bucket_size), -180, 180)), - torsion_bucket_size_(torsion_bucket_size) -{ - -} - -uint32_t TorsionStatisticsSolis::GetTorsionBucketCount() const -{ - return uint32_t(360/torsion_bucket_size_); -} - -int TorsionStatisticsSolis::GetTorsionBucketSize() const -{ - return torsion_bucket_size_; -} - -void TorsionStatisticsSolis::Extract(mol::EntityView view) -{ - TorsionStatSolis stat(histogram_); - view.Apply(stat); -} - -void TorsionStatisticsSolis::Extract(mol::EntityHandle entity) -{ - TorsionStatSolis stat(histogram_); - entity.Apply(stat); -} - -TorsionStatisticsSolisPtr TorsionStatisticsSolis::LoadFromFile(const String& file_name) -{ - std::ifstream stream(file_name.c_str(), std::ios_base::binary); - TorsionStatisticsSolisPtr ptr(new TorsionStatisticsSolis); - io::BinaryDataSource ds(stream); - ds >> *ptr.get(); - return ptr; -} - -void TorsionStatisticsSolis::SaveToFile(const String& file_name) const -{ - std::ofstream stream(file_name.c_str(), std::ios_base::binary); - io::BinaryDataSink ds(stream); - ds << *this; -} - -uint32_t TorsionStatisticsSolis::GetCount(const AminoAcidSet& prev_aa, - const AminoAcidSet& central_aa, - const AminoAcidSet& next_aa, - Real phi_angle, Real psi_angle) const -{ - uint32_t count=0; - uint32_t phi_bucket=this->IndexForAngle(phi_angle); - uint32_t psi_bucket=this->IndexForAngle(psi_angle); - for (AminoAcidSet::Iterator i=prev_aa.Begin(), e1=prev_aa.End(); i!=e1; ++i) { - for (AminoAcidSet::Iterator j=central_aa.Begin(), - e2=central_aa.End(); j!=e2; ++j) { - for (AminoAcidSet::Iterator k=next_aa.Begin(), - e3=next_aa.End(); k!=e3; ++k) { - TorsionHistogram::IndexType start(*i, *j, *k, phi_bucket, psi_bucket), - end(*i, *j, *k, phi_bucket, psi_bucket); - count+=this->Accumulate(start, end); - } - } - } - return count; -} - -uint32_t TorsionStatisticsSolis::GetCount(const AminoAcidSet& prev_aa, - const AminoAcidSet& central_aa, - const AminoAcidSet& next_aa) const -{ - uint32_t count=0; - uint32_t torsion_buckets=uint32_t(360.0/torsion_bucket_size_); - for (AminoAcidSet::Iterator i=prev_aa.Begin(), e1=prev_aa.End(); i!=e1; ++i) { - for (AminoAcidSet::Iterator j=central_aa.Begin(), - e2=central_aa.End(); j!=e2; ++j) { - for (AminoAcidSet::Iterator k=next_aa.Begin(), - e3=next_aa.End(); k!=e3; ++k) { - TorsionHistogram::IndexType start(*i, *j, *k, 0, 0), - end(*i, *j, *k, torsion_buckets-1, - torsion_buckets-1); - count+=this->Accumulate(start, end); - } - } - } - return count; -} - - -uint32_t TorsionStatisticsSolis::GetCount(Real phi_angle, Real psi_angle) const -{ - uint32_t phi_bucket=this->IndexForAngle(phi_angle); - uint32_t psi_bucket=this->IndexForAngle(psi_angle); - TorsionHistogram::IndexType start(0,0,0, phi_bucket, psi_bucket), - end(Xxx-1, Xxx-1, Xxx-1, phi_bucket, psi_bucket); - return this->Accumulate(start, end); -} - -uint32_t TorsionStatisticsSolis::Accumulate(const IndexType& start, - const IndexType& end) const -{ - uint32_t total=0; - for (TorsionHistogram::Iterator i(start, end); !i.AtEnd(); ++i) { - total+=histogram_.Get(*i); - } - return total; -} - - -uint32_t TorsionStatisticsSolis::IndexForAngle(Real angle) const -{ - return (histogram_.FindBucket(0,0,0,angle,0))[3]; -} - -uint32_t TorsionStatisticsSolis::GetCount(AminoAcid central_aa, - Real phi_angle, Real psi_angle) const -{ - uint32_t phi_bucket=this->IndexForAngle(phi_angle); - uint32_t psi_bucket=this->IndexForAngle(psi_angle); - TorsionHistogram::IndexType start(0,central_aa,0, phi_bucket, psi_bucket), - end(Xxx-1, central_aa, Xxx-1, - phi_bucket, psi_bucket); - return this->Accumulate(start, end); -} - - -Real TorsionStatisticsSolis::GetFrequency(AminoAcid aa, Real phi_angle, - Real psi_angle) const -{ - uint32_t phi_bucket=this->IndexForAngle(phi_angle); - uint32_t psi_bucket=this->IndexForAngle(psi_angle); - uint32_t local_count=0, total_count=0; - uint32_t torsion_buckets=uint32_t(360.0/torsion_bucket_size_); - for (uint32_t phi=0; phi<torsion_buckets; ++phi) { - for (uint32_t psi=0; psi<torsion_buckets; ++psi) { - TorsionHistogram::IndexType start(0, aa, 0, phi, psi), - end(Xxx-1, aa, Xxx-1, phi, psi); - uint32_t phi_psi_count=this->Accumulate(start, end); - if (phi==phi_bucket && psi==psi_bucket) { - local_count+=phi_psi_count; - } - total_count+=phi_psi_count; - } - } - return Real(local_count)/Real(total_count); -} - -Real TorsionStatisticsSolis::GetFrequency(Real phi_angle, - Real psi_angle) const -{ - uint32_t phi_bucket=this->IndexForAngle(phi_angle); - uint32_t psi_bucket=this->IndexForAngle(psi_angle); - uint32_t local_count=0, total_count=0; - uint32_t torsion_buckets=uint32_t(360.0/torsion_bucket_size_); - for (uint32_t phi=0; phi<torsion_buckets; ++phi) { - for (uint32_t psi=0; psi<torsion_buckets; ++psi) { - TorsionHistogram::IndexType start(0, 0, 0, phi, psi), - end(Xxx-1, Xxx-1, Xxx-1, phi, psi); - uint32_t phi_psi_count=this->Accumulate(start, end); - if (phi==phi_bucket && psi==psi_bucket) { - local_count+=phi_psi_count; - } - total_count+=phi_psi_count; - } - } - return Real(local_count)/Real(total_count); -} - -int TorsionStatisticsSolis::GetCount() const -{ - uint32_t torsion_buckets=uint32_t(360/torsion_bucket_size_); - TorsionHistogram::IndexType start(0, 0, 0, 0, 0), - end(Xxx-1, Xxx-1, Xxx-1, torsion_buckets-1, - torsion_buckets-1); - return this->Accumulate(start, end); -} - - -TorsionStatisticsSolis::TorsionStatisticsSolis(): - torsion_bucket_size_(0) -{ -} - -}} diff --git a/modules/qa/src/solis_torsion_statistics.hh b/modules/qa/src/solis_torsion_statistics.hh deleted file mode 100644 index 5c189595b2aa88c2e6e76d89752fc945881d1302..0000000000000000000000000000000000000000 --- a/modules/qa/src/solis_torsion_statistics.hh +++ /dev/null @@ -1,127 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#ifndef OST_QA_SOLIS_TORSION_STATISTICS_HH -#define OST_QA_SOLIS_TORSION_STATISTICS_HH -/* - This implementation of the torsion angle potential combines 1 - torsion angle pair with up to 3 residues in analogy to the - potential described by: - Armando D. Solis, S. Rackovsky - Optimally informative backbone structural propensities in proteins - Protei - The ordinary trosion angle potential (1 angle-pair, 1 residue) is covered - as well. - - Author: Marco Biasini, Pascal Benkert - */ -#include <ost/mol/entity_handle.hh> -#include <ost/qa/module_config.hh> -#include <ost/qa/histogram.hh> -#include <ost/qa/amino_acids.hh> -#include <boost/shared_ptr.hpp> - -namespace ost { namespace qa { - -class TorsionStatisticsSolis; -typedef boost::shared_ptr<TorsionStatisticsSolis> TorsionStatisticsSolisPtr; - -/// \brief sequence dependent torsion statistics -class DLLEXPORT_OST_QA TorsionStatisticsSolis { -public: - /// construct new torsion statistics. In order to get something useful - /// you need to either extract statistics from a set of PDBs or the - /// statistics from file. - /// - /// \param torsion_bucket_size is the histogram bucket size for torsion - /// angles. - TorsionStatisticsSolis(int torsion_bucket_size); - TorsionStatisticsSolis(); - /// \brief extract torsion counts from whole entity. - void Extract(mol::EntityHandle entity); - - /// \brief extract torsion counts from (filtered) entity - void Extract(mol::EntityView view); - - /// \brief load torsion statistics from file - static TorsionStatisticsSolisPtr LoadFromFile(const String& file_name); - - void SaveToFile(const String& file_name) const; - - /// \brief Get count of phi/psi torsions with the specified torsion - /// bucket size away conditional on the central, previous - /// and next amino acid. - uint32_t GetCount(const AminoAcidSet& prev_aa, - const AminoAcidSet& central_aa, - const AminoAcidSet& next_aa, - Real phi_angle, Real psi_angle) const; - - uint32_t GetCount(const AminoAcidSet& prev_aa, - const AminoAcidSet& central_aa, - const AminoAcidSet& next_aa) const; - - /// \brief Get count of amino acids with the given phi/psi angle. - uint32_t GetCount(Real phi_angle, Real psi_angle) const; - - /// \brief Get count of amino acids with the given phi/psi angle - /// conditional on the central residue. - uint32_t GetCount(AminoAcid central_aa, - Real phi_angle, Real psi_angle) const; - - /// \brief Get frequency of phi/psi angle in bucket phi_angle/psi_angle when - /// the amino acid is aa. - /// \return frequency in the range of [0, 1] - Real GetFrequency(AminoAcid aa, Real phi_angle, Real psi_angle) const; - - /// \brief Get frequency of phi/psi angle in bucket phi_angle/psi_angle - /// regardless of the amino acid type. - /// \return frequency in the range of [0, 1] - Real GetFrequency(Real phi_angle, Real psi_angle) const; - - /// \brief get the number of torsion buckets for every of the torsion - /// angles. - uint32_t GetTorsionBucketCount() const; - int GetTorsionBucketSize() const; - - int GetCount() const; - /// \internal - template <typename DS> - void Serialize(DS& ds) - { - ds & torsion_bucket_size_; - ds & histogram_; - } -public: - typedef MultiClassifier<int, int, int, int, Real, Real> TorsionHistogram; - typedef TorsionHistogram::IndexType IndexType; -private: - /// \brief Convenience function to add counts in the inclusive range - /// [start, end] - uint32_t Accumulate(const IndexType& start, const IndexType& end) const; - - uint32_t IndexForAngle(Real angle) const; - - - TorsionHistogram histogram_; - int torsion_bucket_size_; -}; -}} - -#endif - - diff --git a/modules/qa/src/torsion_potential.cc b/modules/qa/src/torsion_potential.cc deleted file mode 100644 index d3fd7bcd81f5aea49b5bc6d6aa4c1a071a996cb8..0000000000000000000000000000000000000000 --- a/modules/qa/src/torsion_potential.cc +++ /dev/null @@ -1,321 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ - -#include "torsion_potential.hh" -#include "amino_acids.hh" - -#include <ost/log.hh> -#include <ost/io/binary_data_source.hh> -#include <ost/io/binary_data_sink.hh> -#include <ost/io/container_serialization.hh> -#include <ost/integrity_error.hh> -#include <ost/mol/mol.hh> -#include <boost/filesystem/convenience.hpp> -#include <ost/message.hh> -#include <ost/io/io_exception.hh> - -/* - Authors: Marco Biasini, Pascal Benkert - */ - -namespace ost { namespace qa { - -namespace { - -class TorsionEnergyCalc : public mol::EntityVisitor { -public: - TorsionEnergyCalc(const TorsionPotentialPtr& pot, - TorsionPotentialOpts opts, int length): - pot_(pot), energy_(0.0), num_torsions_(0), score_vector_(std::vector<float>(length, 0.0)), - pos_counter_(0) - { - } - - virtual bool VisitResidue(const mol::ResidueHandle& rh) - { - pos_counter_++; - mol::ResidueHandle prev=rh.GetPrev(); - mol::ResidueHandle next=rh.GetNext(); - if (!(next && prev && next.IsPeptideLinking() && - prev.IsPeptideLinking() && rh.IsPeptideLinking())) { - return false; - } - - AminoAcid ca=ResidueToAminoAcid(rh); - if (ca==Xxx) - return false; - - mol::TorsionHandle prev_phit=prev.GetPhiTorsion(); - mol::TorsionHandle prev_psit=prev.GetPsiTorsion(); - mol::TorsionHandle central_phit=rh.GetPhiTorsion(); - mol::TorsionHandle central_psit=rh.GetPsiTorsion(); - mol::TorsionHandle next_phit=next.GetPhiTorsion(); - mol::TorsionHandle next_psit=next.GetPsiTorsion(); - if (!prev_phit || !prev_psit || !central_phit || !central_psit || - !next_phit || !next_psit) - return false; - Real prev_phi=prev_phit.GetAngle()*180/M_PI - 0.001; if (prev_phi<-180) prev_phi=-180; - Real prev_psi=prev_psit.GetAngle()*180/M_PI - 0.001; if (prev_psi<-180) prev_psi=-180; - Real central_phi=central_phit.GetAngle()*180/M_PI - 0.001; if (central_phi<-180) central_phi=-180; - Real central_psi=central_psit.GetAngle()*180/M_PI - 0.001; if (central_psi<-180) central_psi=-180; - Real next_phi=next_phit.GetAngle()*180/M_PI - 0.001; if (next_phi<-180) next_phi=-180; - Real next_psi=next_psit.GetAngle()*180/M_PI - 0.001; if (next_psi<-180) next_psi=-180; - - - // calculate position of the amino acid in the alphabet - float local_energy=pot_->GetTorsionEnergy(ca, prev_phi, prev_psi, central_phi, - central_psi, next_phi, next_psi); - energy_+=local_energy; - - score_vector_[pos_counter_-1]=local_energy; - ++num_torsions_; - return false; - } - - Real GetEnergy() const { - return energy_; - } - - int GetEnergyCounts() const { - return num_torsions_; - } - - std::vector<float> GetResidueEnergyVector() const { - return score_vector_; - } - -private: - TorsionPotentialPtr pot_; - AminoAcid prev_; - AminoAcid center_; - mol::ResidueHandle cr_; - Real energy_; - int num_torsions_; - std::vector<float> score_vector_; - int pos_counter_; -}; - -} - -template <typename DS> -void TorsionPotentialOpts::Serialize(DS& ds) -{ - ds & prev_angular_bucket_size; - ds & central_angular_bucket_size; - ds & next_angular_bucket_size; - ds & alphabet; - ds & sigma; -} - -template <typename DS> -void TorsionPotential::Serialize(DS& ds) -{ - ds & options_; - ds & energies_; -} - -TorsionPotentialOpts::TorsionPotentialOpts(): - prev_angular_bucket_size(90), central_angular_bucket_size(45), - next_angular_bucket_size(90), sigma(0.0005) -{ - -} - - -int TorsionPotential::GetAAIndex(AminoAcid aa) const -{ - AminoAcidAlphabet::const_iterator i=options_.alphabet.begin(); - for (int index=0; i!=options_.alphabet.end(); ++i, ++index) { - if ((*i).Contains(aa)) { - return index; - } - } - return -1; -} - - -TorsionPotentialPtr TorsionPotential::Create(TorsionStatisticsPtr statistics, - const TorsionPotentialOpts& opts, - bool calculate_average_energy_flag) -{ - - typedef TorsionPotential::TorsionEnergies TorsionEnergies; - int prev_bucket_size=statistics->GetTorsionBucketSizePrev(); - int central_bucket_size=statistics->GetTorsionBucketSizeCentral(); - int next_bucket_size=statistics->GetTorsionBucketSizeNext(); - - //std::cout << prev_bucket_size << " " << central_bucket_size << " " << next_bucket_size << std::endl; - - if (prev_bucket_size == 0 || central_bucket_size == 0 || next_bucket_size == 0) { - throw IntegrityError("Torsion angle statistics does not exist. Re-extract counts."); - } - - if ((opts.prev_angular_bucket_size % prev_bucket_size)!=0) { - throw IntegrityError("Previous angular bucket size for torsion potential is not an " - "integral multiple of the statistics bucket size"); - } - if ((opts.central_angular_bucket_size % central_bucket_size)!=0) { - throw IntegrityError("Central angular bucket size for torsion potential is not an " - "integral multiple of the statistics bucket size"); - } - if ((opts.next_angular_bucket_size % next_bucket_size)!=0) { - throw IntegrityError("Next angular bucket size for torsion potential is not an " - "integral multiple of the statistics bucket size"); - } - int prev_buckets=360/opts.prev_angular_bucket_size; - int central_buckets=360/opts.central_angular_bucket_size; - int next_buckets=360/opts.next_angular_bucket_size; - TorsionPotentialPtr ptr(new TorsionPotential); - - ptr->options_.prev_angular_bucket_size = opts.prev_angular_bucket_size; - ptr->options_.central_angular_bucket_size = opts.central_angular_bucket_size; - ptr->options_.next_angular_bucket_size = opts.next_angular_bucket_size; - - ptr->options_.alphabet = opts.alphabet; - uint64_t alphabet_size=ptr->options_.alphabet.size(); - - ptr->energies_=TorsionEnergies(0.0, IntegralClassifier(alphabet_size, 0), - ContinuousClassifier(prev_buckets,-180,180), - ContinuousClassifier(prev_buckets,-180,180), - ContinuousClassifier(central_buckets,-180,180), - ContinuousClassifier(central_buckets,-180,180), - ContinuousClassifier(next_buckets,-180,180), - ContinuousClassifier(next_buckets,-180,180)); - ptr->Fill(statistics, calculate_average_energy_flag); - return ptr; -} - -TorsionPotentialPtr TorsionPotential::LoadFromFile(const String& path) -{ - if(!boost::filesystem::exists(path)) - throw io::IOException("Could not open torsion potential data file.\nFile does not exist at: "+path); - - std::ifstream stream(path.c_str(), std::ios_base::binary); - io::BinaryDataSource ds(stream); - TorsionPotentialPtr ptr(new TorsionPotential); - ds >> *ptr.get(); - return ptr; -} - -void TorsionPotential::SaveToFile(const String& path) -{ - std::ofstream stream(path.c_str(), std::ios_base::binary); - io::BinaryDataSink ds(stream); - ds << *this; -} - -Real TorsionPotential::GetTotalEnergy(mol::EntityHandle entity) -{ - int res_count=entity.GetResidueList().size(); - TorsionEnergyCalc c(this->shared_from_this(), options_, res_count); - entity.Apply(c); - num_torsions_ = c.GetEnergyCounts(); - std::vector<float> rev=c.GetResidueEnergyVector(); - for(int i=0;i<res_count;++i) { - entity.GetResidueList()[i].SetFloatProp("torsion_energy", rev[i]); - } - return c.GetEnergy(); -} - -Real TorsionPotential::GetTotalEnergy(mol::EntityView entity) -{ - int res_count=entity.GetResidueList().size(); - TorsionEnergyCalc c(this->shared_from_this(), options_, res_count); - entity.Apply(c); - num_torsions_ = c.GetEnergyCounts(); - std::vector<float> rev=c.GetResidueEnergyVector(); - for(int i=0;i<res_count;++i) { - entity.GetResidueList()[i].SetFloatProp("torsion_energy", rev[i]); - } - return c.GetEnergy(); -} - -int TorsionPotential::GetEnergyCounts() const { - return num_torsions_; -} - -void TorsionPotential::Fill(const TorsionStatisticsPtr& stat, bool calculate_average_energy_flag) -{ - int prev_num_b=360/(options_.prev_angular_bucket_size); - int central_num_b=360/(options_.central_angular_bucket_size); - int next_num_b=360/(options_.next_angular_bucket_size); - - int prev_half_b=options_.prev_angular_bucket_size/2; - int central_half_b=options_.central_angular_bucket_size/2; - int next_half_b=options_.next_angular_bucket_size/2; - - uint64_t M=stat->GetCount(); - //std::cout << "total counts: " << M << std::endl; - - typedef AminoAcidAlphabet::const_iterator AAA; - for (AAA i=options_.alphabet.begin(), - e1=options_.alphabet.end(); i!=e1; ++i) { - - uint64_t Mi=stat->GetCount(*i); - if (calculate_average_energy_flag==1) std::cout << *i << " " << float(Mi)/M << " "; - Real per_aa_energy = 0; // for average energy - - - for (int prev_phi=0; prev_phi<prev_num_b; ++prev_phi) { - for (int prev_psi=0; prev_psi<prev_num_b; ++prev_psi) { - int prev_phi_a=-180+options_.prev_angular_bucket_size*prev_phi+prev_half_b; - int prev_psi_a=-180+options_.prev_angular_bucket_size*prev_psi+prev_half_b; - - for (int central_phi=0; central_phi<central_num_b; ++central_phi) { - for (int central_psi=0; central_psi<central_num_b; ++central_psi) { - int central_phi_a=-180+options_.central_angular_bucket_size*central_phi+central_half_b; - int central_psi_a=-180+options_.central_angular_bucket_size*central_psi+central_half_b; - - for (int next_phi=0; next_phi<next_num_b; ++next_phi) { - for (int next_psi=0; next_psi<next_num_b; ++next_psi) { - int next_phi_a=-180+options_.next_angular_bucket_size*next_phi+next_half_b; - int next_psi_a=-180+options_.next_angular_bucket_size*next_psi+next_half_b; - - uint64_t fxi=stat->GetCount(*i, prev_phi_a, prev_psi_a, central_phi_a, central_psi_a, next_phi_a, next_psi_a); - uint64_t fx=stat->GetCount(prev_phi_a, prev_psi_a, central_phi_a, central_psi_a, next_phi_a, next_psi_a); - - // propensity = (fxi/Mi)/(fx/M) - Real propensity=0.0; - Real energy=0.0; - //avoid division by zero: - if (fx != 0 && Mi != 0 && M != 0) { - propensity=Real(Real(fxi)/Mi)/(Real(fx)/M); - } - - energy=log(1+options_.sigma*Mi)- - log(1+options_.sigma*Mi*propensity); - - per_aa_energy=per_aa_energy+fxi*energy; // for average energy - - - energies_.Add(energy, - i-options_.alphabet.begin(), - prev_phi_a, prev_psi_a, central_phi_a, central_psi_a, next_phi_a, next_psi_a); - } - } - - } - } - - } - } - if (calculate_average_energy_flag==1) std::cout << per_aa_energy/Mi << std::endl; -} -} -}} diff --git a/modules/qa/src/torsion_potential.hh b/modules/qa/src/torsion_potential.hh deleted file mode 100644 index 91b5407db4a9877460c033e31f78b36b0db82aa5..0000000000000000000000000000000000000000 --- a/modules/qa/src/torsion_potential.hh +++ /dev/null @@ -1,137 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#ifndef OST_QA_TORSION_POTENTIAL_HH -#define OST_QA_TORSION_POTENTIAL_HH -/* - This implementation of the torsion angle potential combines 3 - torsion angle pairs with 1 residues in analogy to the combined - potential described by: - Pascal Benkert, Silvio C. E. Tosatto, Dietmar Schomburg - QMEAN: A comprehensive scoring function for model quality assessment - Proteins, 2008 - - Author: Pascal, Benkert, Marco Biasini -*/ -#include <ost/qa/torsion_statistics.hh> -#include <boost/shared_ptr.hpp> -#include <boost/enable_shared_from_this.hpp> -#include <vector> -#include <boost/scoped_ptr.hpp> - -namespace ost { namespace qa { - -typedef std::vector<AminoAcidSet> AminoAcidAlphabet; - -/// \brief torsion potential options -/// -/// These options may be used to configure the different torsion potentials -/// based on a sliding window of 3 amino acids. -struct DLLEXPORT_OST_QA TorsionPotentialOpts { - /// \brief initialize torsion potential with angular bucket size of 30 and - /// sigma of 0.02 - /// - /// The amino acid alphabets are empty. - TorsionPotentialOpts(); - - /// angular bucket size, 360 % angular_bucket_size must be zero - int prev_angular_bucket_size; - int central_angular_bucket_size; - int next_angular_bucket_size; - /// alphabet of the central. The amino acids in the set must be mutually excl. - AminoAcidAlphabet alphabet; - /// weighting factor after (Sippl et al., 1990), usually 0.02 - Real sigma; - /// \internal - template <typename DS> - void Serialize(DS& ds); -}; - -class TorsionPotential; -typedef boost::shared_ptr<TorsionPotential> TorsionPotentialPtr; - -/// \brief Torsion potential -/// -/// The torsion potential class is parametrisable by TorsionPotentialOpts. -class DLLEXPORT_OST_QA TorsionPotential : - public boost::enable_shared_from_this<TorsionPotential> { -public: - /// \brief create new torsion potential with the given torsion statistics - /// and options - /// If the flag ist set the average energy per is calculated - /// and shown on stdout. May be used to correct compositional bias in - /// scoring functions (e.g. Pro and Gly have on average lower torsion energies). - static TorsionPotentialPtr Create(TorsionStatisticsPtr statistics, - const TorsionPotentialOpts& opts, - bool calculate_average_energy_flag); - /// \brief Load potential from file. - static TorsionPotentialPtr LoadFromFile(const String& path); - - /// \brief Calculate energy for entity. - Real GetTotalEnergy(mol::EntityHandle entity); - - /// \brief Calculate energy for entity view. Only residues in the view are - /// considered. - Real GetTotalEnergy(mol::EntityView view); - - /// \brief retrieve total number of energy local (i.e. valid residues) - int GetEnergyCounts() const; - - - - Real GetTorsionEnergy(AminoAcid central_aa, Real prev_phi, Real prev_psi, - Real central_phi, Real central_psi, - Real next_phi, Real next_psi) const { - int icenter=this->GetAAIndex(central_aa); - return energies_.Get(icenter, prev_phi, prev_psi, - central_phi, central_psi, - next_phi, next_psi); - } - - - /// \brief save torsion potential - /// - /// The output is in a non-portable binary form. - void SaveToFile(const String& path); - - /// \internal - template <typename DS> - void Serialize(DS& ds); - - -public: - /// Meaning of parameters: AA before AA center AA after, - /// phi torsion angle, psi torsion angle - typedef MultiClassifier<float, int, Real, Real, - Real, Real, Real, Real> TorsionEnergies; -private: - - int GetAAIndex(AminoAcid aa) const; - - void Fill(const TorsionStatisticsPtr& stat, - bool calculate_average_energy_flag); - - TorsionPotentialOpts options_; - TorsionEnergies energies_; - int num_torsions_; -}; - - -}} - -#endif diff --git a/modules/qa/src/torsion_statistics.cc b/modules/qa/src/torsion_statistics.cc deleted file mode 100644 index fd67b7f864e1f99e92d58f84201573f3e50e3564..0000000000000000000000000000000000000000 --- a/modules/qa/src/torsion_statistics.cc +++ /dev/null @@ -1,271 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#include "torsion_statistics.hh" -#include "amino_acids.hh" -#include <ost/log.hh> -#include <ost/mol/mol.hh> -#include <fstream> -#include <ost/message.hh> -#include <ost/io/binary_data_source.hh> -#include <ost/io/binary_data_sink.hh> -#include <ost/io/container_serialization.hh> - -/* - Author: Marco Biasini, Pascal Benkert - */ -namespace ost { namespace qa { - - -namespace { - -class TorsionStat : public mol::EntityVisitor { -public: - TorsionStat(TorsionStatistics::TorsionHistogram& histo): - histo_(histo) - {} - - virtual bool VisitResidue(const mol::ResidueHandle& residue) { - - mol::ResidueHandle prev=residue.GetPrev(); - mol::ResidueHandle next=residue.GetNext(); - if (!(next && prev && next.IsPeptideLinking() && - prev.IsPeptideLinking() && residue.IsPeptideLinking())) { - return false; - } - - AminoAcid ca=ResidueToAminoAcid(residue); - if (ca==Xxx) - return false; - - mol::TorsionHandle prev_phit=prev.GetPhiTorsion(); - mol::TorsionHandle prev_psit=prev.GetPsiTorsion(); - mol::TorsionHandle central_phit=residue.GetPhiTorsion(); - mol::TorsionHandle central_psit=residue.GetPsiTorsion(); - mol::TorsionHandle next_phit=next.GetPhiTorsion(); - mol::TorsionHandle next_psit=next.GetPsiTorsion(); - if (!prev_phit || !prev_psit || !central_phit || !central_psit || - !next_phit || !next_psit) - return false; - Real prev_phi=prev_phit.GetAngle()*180/M_PI - 0.001; - Real prev_psi=prev_psit.GetAngle()*180/M_PI - 0.001; - Real central_phi=central_phit.GetAngle()*180/M_PI - 0.001; - Real central_psi=central_psit.GetAngle()*180/M_PI - 0.001; - Real next_phi=next_phit.GetAngle()*180/M_PI - 0.001; - Real next_psi=next_psit.GetAngle()*180/M_PI - 0.001; - - histo_.Add(1, ca, prev_phi, prev_psi, central_phi, central_psi, next_phi, next_psi); - return false; - } -private: - TorsionStatistics::TorsionHistogram& histo_; -}; - -} - -TorsionStatistics::TorsionStatistics(int prev_torsion_bucket_size, - int central_torsion_bucket_size, - int next_torsion_bucket_size) - : histogram_(1, IntegralClassifier(20, 0), - ContinuousClassifier(int(360/prev_torsion_bucket_size), -180, 180), - ContinuousClassifier(int(360/prev_torsion_bucket_size), -180, 180), - ContinuousClassifier(int(360/central_torsion_bucket_size), -180, 180), - ContinuousClassifier(int(360/central_torsion_bucket_size), -180, 180), - ContinuousClassifier(int(360/next_torsion_bucket_size), -180, 180), - ContinuousClassifier(int(360/next_torsion_bucket_size), -180, 180)), - prev_torsion_bucket_size_(prev_torsion_bucket_size), - central_torsion_bucket_size_(central_torsion_bucket_size), - next_torsion_bucket_size_(next_torsion_bucket_size) -{ - -} - -uint32_t TorsionStatistics::GetTorsionBucketCountPrev() const -{ - return uint32_t(360/prev_torsion_bucket_size_); -} - -uint32_t TorsionStatistics::GetTorsionBucketCountCentral() const -{ - return uint32_t(360/central_torsion_bucket_size_); -} - -uint32_t TorsionStatistics::GetTorsionBucketCountNext() const -{ - return uint32_t(360/next_torsion_bucket_size_); -} - - -int TorsionStatistics::GetTorsionBucketSizePrev() const -{ - return prev_torsion_bucket_size_; -} - -int TorsionStatistics::GetTorsionBucketSizeCentral() const -{ - return central_torsion_bucket_size_; -} - -int TorsionStatistics::GetTorsionBucketSizeNext() const -{ - return next_torsion_bucket_size_; -} - - -void TorsionStatistics::Extract(mol::EntityView view) -{ - TorsionStat stat(histogram_); - view.Apply(stat); -} - -void TorsionStatistics::Extract(mol::EntityHandle entity) -{ - TorsionStat stat(histogram_); - entity.Apply(stat); -} - -TorsionStatisticsPtr TorsionStatistics::LoadFromFile(const String& file_name) -{ - std::ifstream stream(file_name.c_str(), std::ios_base::binary); - TorsionStatisticsPtr ptr(new TorsionStatistics); - io::BinaryDataSource ds(stream); - ds >> *ptr.get(); - return ptr; -} - -void TorsionStatistics::SaveToFile(const String& file_name) const -{ - std::ofstream stream(file_name.c_str(), std::ios_base::binary); - io::BinaryDataSink ds(stream); - ds << *this; -} - - -uint64_t TorsionStatistics::GetCount(const AminoAcidSet& central_aa, - Real prev_phi_angle, Real prev_psi_angle, - Real central_phi_angle, Real central_psi_angle, - Real next_phi_angle, Real next_psi_angle) const -{ - uint64_t count=0; - uint32_t prev_phi_bucket=this->IndexForAnglePrev(prev_phi_angle); - uint32_t prev_psi_bucket=this->IndexForAnglePrev(prev_psi_angle); - uint32_t central_phi_bucket=this->IndexForAngleCentral(central_phi_angle); - uint32_t central_psi_bucket=this->IndexForAngleCentral(central_psi_angle); - uint32_t next_phi_bucket=this->IndexForAngleNext(next_phi_angle); - uint32_t next_psi_bucket=this->IndexForAngleNext(next_psi_angle); - for (AminoAcidSet::Iterator i=central_aa.Begin(), e1=central_aa.End(); i!=e1; ++i) { - TorsionHistogram::IndexType start(*i, prev_phi_bucket, prev_psi_bucket, - central_phi_bucket, central_psi_bucket, - next_phi_bucket, next_psi_bucket), - end(*i, prev_phi_bucket, prev_psi_bucket, - central_phi_bucket, central_psi_bucket, - next_phi_bucket, next_psi_bucket); - count+=this->Accumulate(start, end); - } - return count; -} - - -uint64_t TorsionStatistics::GetCount(Real prev_phi_angle, Real prev_psi_angle, - Real central_phi_angle, Real central_psi_angle, - Real next_phi_angle, Real next_psi_angle) const -{ - uint32_t prev_phi_bucket=this->IndexForAnglePrev(prev_phi_angle); - uint32_t prev_psi_bucket=this->IndexForAnglePrev(prev_psi_angle); - uint32_t central_phi_bucket=this->IndexForAngleCentral(central_phi_angle); - uint32_t central_psi_bucket=this->IndexForAngleCentral(central_psi_angle); - uint32_t next_phi_bucket=this->IndexForAngleNext(next_phi_angle); - uint32_t next_psi_bucket=this->IndexForAngleNext(next_psi_angle); - TorsionHistogram::IndexType start(0, prev_phi_bucket, prev_psi_bucket, - central_phi_bucket, central_psi_bucket, - next_phi_bucket, next_psi_bucket), - end(Xxx-1, prev_phi_bucket, prev_psi_bucket, - central_phi_bucket, central_psi_bucket, - next_phi_bucket, next_psi_bucket); - return this->Accumulate(start, end); -} - - -uint64_t TorsionStatistics::GetCount(const AminoAcidSet& central_aa) const -{ - uint64_t count=0; - uint32_t prev_buckets=uint32_t(360.0/prev_torsion_bucket_size_); - uint32_t central_buckets=uint32_t(360.0/central_torsion_bucket_size_); - uint32_t next_buckets=uint32_t(360.0/next_torsion_bucket_size_); - - for (AminoAcidSet::Iterator i=central_aa.Begin(), e1=central_aa.End(); i!=e1; ++i) { - TorsionHistogram::IndexType start(*i, 0, 0, 0, 0, 0, 0), - end(*i, prev_buckets-1, prev_buckets-1, - central_buckets-1, central_buckets-1, - next_buckets-1, next_buckets-1); - count+=this->Accumulate(start, end); - } - return count; -} - - -uint64_t TorsionStatistics::GetCount() const -{ - uint32_t prev_buckets=uint32_t(360.0/prev_torsion_bucket_size_); - uint32_t central_buckets=uint32_t(360.0/central_torsion_bucket_size_); - uint32_t next_buckets=uint32_t(360.0/next_torsion_bucket_size_); - - TorsionHistogram::IndexType start(0, 0, 0, 0, 0, 0, 0), - end(Xxx-1, prev_buckets-1, prev_buckets-1, - central_buckets-1, central_buckets-1, - next_buckets-1, next_buckets-1); - return this->Accumulate(start, end); -} - - -uint64_t TorsionStatistics::Accumulate(const IndexType& start, - const IndexType& end) const -{ - uint64_t total=0; - for (TorsionHistogram::Iterator i(start, end); !i.AtEnd(); ++i) { - total+=histogram_.Get(*i); - } - return total; -} - - -uint32_t TorsionStatistics::IndexForAnglePrev(Real angle) const -{ - return (histogram_.FindBucket(0,angle,0,0,0,0,0))[1]; -} - -uint32_t TorsionStatistics::IndexForAngleCentral(Real angle) const -{ - return (histogram_.FindBucket(0,0,0,angle,0,0,0))[3]; -} - -uint32_t TorsionStatistics::IndexForAngleNext(Real angle) const -{ - return (histogram_.FindBucket(0,0,0,0,0,angle,0))[5]; -} - - - -TorsionStatistics::TorsionStatistics(): - prev_torsion_bucket_size_(0), - central_torsion_bucket_size_(0), - next_torsion_bucket_size_(0) -{ -} - -} } diff --git a/modules/qa/src/torsion_statistics.hh b/modules/qa/src/torsion_statistics.hh deleted file mode 100644 index 73105e6e5c0d7f2fd2e16740e7f6150f04191ebb..0000000000000000000000000000000000000000 --- a/modules/qa/src/torsion_statistics.hh +++ /dev/null @@ -1,129 +0,0 @@ -//------------------------------------------------------------------------------ -// This file is part of the OpenStructure project <www.openstructure.org> -// -// Copyright (C) 2008-2011 by the OpenStructure authors -// -// This library is free software; you can redistribute it and/or modify it under -// the terms of the GNU Lesser General Public License as published by the Free -// Software Foundation; either version 3.0 of the License, or (at your option) -// any later version. -// This library is distributed in the hope that it will be useful, but WITHOUT -// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more -// details. -// -// You should have received a copy of the GNU Lesser General Public License -// along with this library; if not, write to the Free Software Foundation, Inc., -// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA -//------------------------------------------------------------------------------ -#ifndef OST_QA_TORSION_STATISTICS_HH -#define OST_QA_TORSION_STATISTICS_HH -/* - Author: Pascal, Benkert, Marco Biasini -*/ -#include <ost/mol/entity_handle.hh> -#include <ost/qa/histogram.hh> -#include <ost/qa/amino_acids.hh> -#include <boost/shared_ptr.hpp> - -namespace ost { namespace qa { - -class TorsionStatistics; -typedef boost::shared_ptr<TorsionStatistics> TorsionStatisticsPtr; - -/// \brief sequence dependent torsion statistics -/// -/// This implementation of the torsion angle potential combines 3 torsion angle -/// pairs with 1 residues in analogy to the combined potential described by: -/// Pascal Benkert, Silvio C. E. Tosatto, Dietmar Schomburg QMEAN: A -/// comprehensive scoring function for model quality assessment Proteins, 2008 -class DLLEXPORT_OST_QA TorsionStatistics { -public: - /// construct new torsion statistics. In order to get something useful - /// you need to either extract statistics from a set of PDBs or the - /// statistics from file. - /// - /// \param prev_torsion_bucket_size is torsion bucket size for the first - /// torsion - /// \param central_torsion_bucket_size is the bucket size for the second - /// torsion - /// \param next_torsion_bucket_size is the bucket size for the third torsion - TorsionStatistics(int prev_torsion_bucket_size, int central_torsion_bucket_size, - int next_torsion_bucket_size); - TorsionStatistics(); - /// \brief extract torsion counts from whole entity. - void Extract(mol::EntityHandle entity); - - /// \brief extract torsion counts from (filtered) entity - void Extract(mol::EntityView view); - - /// \brief load torsion statistics from file - static TorsionStatisticsPtr LoadFromFile(const String& file_name); - - void SaveToFile(const String& file_name) const; - - /// \brief Get count of phi/psi torsions with the specified torsion - /// bucket size away conditional on the central amino acid. - uint64_t GetCount(const AminoAcidSet& central_aa, - Real prev_phi_angle, Real prev_psi_angle, - Real central_phi_angle, Real central_psi_angle, - Real next_phi_angle, Real next_psi_angle) const; - - uint64_t GetCount(Real prev_phi_angle, Real prev_psi_angle, - Real central_phi_angle, Real central_psi_angle, - Real next_phi_angle, Real next_psi_angle) const; - - uint64_t GetCount(const AminoAcidSet& central_aa) const; - - uint64_t GetCount() const; - -// /// \brief Get frequency of the amino acid aa. -// /// \return frequency in the range of [0, 1] -// Real GetFrequency(AminoAcid aa) const; - - - /// \brief get the number of torsion buckets for every of the torsion - /// angles. - uint32_t GetTorsionBucketCountPrev() const; - int GetTorsionBucketSizePrev() const; - uint32_t GetTorsionBucketCountCentral() const; - int GetTorsionBucketSizeCentral() const; - uint32_t GetTorsionBucketCountNext() const; - int GetTorsionBucketSizeNext() const; - - - /// \internal - template <typename DS> - void Serialize(DS& ds) - { - ds & prev_torsion_bucket_size_; - ds & central_torsion_bucket_size_; - ds & next_torsion_bucket_size_; - ds & histogram_; - } -public: - - typedef MultiClassifier<uint32_t, int, Real, Real, - Real, Real, Real, Real> TorsionHistogram; - typedef TorsionHistogram::IndexType IndexType; - - -private: - /// \brief Convenience function to add counts in the inclusive range - /// [start, end] - uint64_t Accumulate(const IndexType& start, const IndexType& end) const; - - uint32_t IndexForAnglePrev(Real angle) const; - uint32_t IndexForAngleCentral(Real angle) const; - uint32_t IndexForAngleNext(Real angle) const; - - TorsionHistogram histogram_; - int prev_torsion_bucket_size_; - int central_torsion_bucket_size_; - int next_torsion_bucket_size_; -}; -}} - -#endif - -