diff --git a/modules/qa/pymod/CMakeLists.txt b/modules/qa/pymod/CMakeLists.txt
index c4fff8b16c9a6f8fc579717f9d2ccd988c59002c..a63e9c2d9bc05f8a671c8bb172aac63cf7c8ccd5 100644
--- a/modules/qa/pymod/CMakeLists.txt
+++ b/modules/qa/pymod/CMakeLists.txt
@@ -2,6 +2,7 @@ set(OST_QA_PYMOD_SOURCES
export_interaction.cc
export_torsion.cc
export_packing.cc
+ export_reduced.cc
wrap_qa.cc
export_clash.cc
)
diff --git a/modules/qa/pymod/export_reduced.cc b/modules/qa/pymod/export_reduced.cc
new file mode 100644
index 0000000000000000000000000000000000000000..0c230fb7d74aa88be2531ff122cadae94aeacd83
--- /dev/null
+++ b/modules/qa/pymod/export_reduced.cc
@@ -0,0 +1,63 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2010 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 boost::python;
+
+uint64_t (ReducedStatistics::*get_count)(AminoAcid,AminoAcid,int,int)=&ReducedStatistics::GetCount;
+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", &ReducedStatistics::Extract)
+ .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", &ReducedPotential::GetTotalEnergy,
+ (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
index a02e46f6c67d00ba90d9e65e353426affd948ae9..ed430c871e4369042b7b77127afbd66ff7656df7 100644
--- a/modules/qa/pymod/export_torsion.cc
+++ b/modules/qa/pymod/export_torsion.cc
@@ -84,6 +84,7 @@ void export_Torsion()
.value("His", His)
.value("Phe", Phe)
.value("Trp", Trp)
+ .export_values()
;
class_<AminoAcidSet>("AminoAcidSet", init<>())
diff --git a/modules/qa/pymod/wrap_qa.cc b/modules/qa/pymod/wrap_qa.cc
index bb6519155ef12476d31c1e1fac7752069ca182d2..8020ed29e5c1714d65911f289a7bf31cb998e912 100644
--- a/modules/qa/pymod/wrap_qa.cc
+++ b/modules/qa/pymod/wrap_qa.cc
@@ -22,10 +22,12 @@ void export_Torsion();
void export_Interaction();
void export_Packing();
void export_Clash();
+void export_Reduced();
BOOST_PYTHON_MODULE(_qa)
{
export_Torsion();
export_Interaction();
export_Packing();
export_Clash();
+ export_Reduced();
}
diff --git a/modules/qa/src/CMakeLists.txt b/modules/qa/src/CMakeLists.txt
index cc75ea3b85e8c3a6de3f98aeb75732ccf14f0849..b4c7baf97796e995ac5d3aca326dea2ec45154b1 100644
--- a/modules/qa/src/CMakeLists.txt
+++ b/modules/qa/src/CMakeLists.txt
@@ -13,10 +13,14 @@ 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
@@ -28,8 +32,10 @@ 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 ${OST_QA_HEADERS} DEPENDS_ON io)
+ HEADERS reduced_impl.hh IN_DIR impl ${OST_QA_HEADERS} DEPENDS_ON io)
diff --git a/modules/qa/src/impl/reduced_impl.cc b/modules/qa/src/impl/reduced_impl.cc
new file mode 100644
index 0000000000000000000000000000000000000000..0efa531b3e3befc326b22489b39e84f90a3af8ca
--- /dev/null
+++ b/modules/qa/src/impl/reduced_impl.cc
@@ -0,0 +1,91 @@
+#include <ost/log.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;
+ }
+
+ 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) {
+ mol::ResidueHandle res_two=i->GetResidue();
+ if (res_two.GetIndex()-index<opts_.sequence_sep) {
+ continue;
+ }
+ if (i->GetName()!="CA" || !res_two.IsPeptideLinking()) {
+ continue;
+ }
+
+
+ AminoAcid aa_two=OneLetterCodeToAminoAcid(res_two.GetOneLetterCode());
+ if (aa_two==Xxx) {
+ continue;
+ }
+ geom::Vec3 ca_pos_two;
+ geom::Vec3 cb_pos_two;
+ if (!this->GetCAlphaCBetaPos(res_two, ca_pos_two, cb_pos_two)) {
+ continue;
+ }
+
+ Real dist=geom::Length(ca_pos_one-ca_pos_two);
+ if (dist<opts_.lower_cutoff) {
+ continue;
+ }
+ Real angle=geom::Angle(cb_pos_one-ca_pos_one, cb_pos_two-ca_pos_two);
+
+ this->OnInteraction(aa_one, aa_two, dist, angle);
+ }
+ return false;
+}
+
+bool ReducedPotentialImpl::GetCAlphaCBetaPos(const mol::ResidueHandle& res,
+ geom::Vec3& ca_pos,
+ geom::Vec3& cb_pos)
+ {
+ const static Real bond_length=1.5;
+ 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;
+ }
+ geom::Vec3 v1=geom::Normalize(ca.GetPos()-n.GetPos());
+ geom::Vec3 v2=geom::Normalize(ca.GetPos()-c.GetPos());
+ geom::Vec3 in_plane_v=geom::Normalize(v1+v2);
+ geom::Plane p(ca.GetPos() ,n.GetPos(), c.GetPos());
+ // rotate around vector perpendicular to p and in_plane_v
+ geom::Vec3 axis=geom::Normalize(geom::Cross(p.GetNormal(), in_plane_v));
+ geom::Mat3 rot_mat=geom::AxisRotation(axis, (-54/180.0)*M_PI);
+ cb_pos=ca.GetPos()+rot_mat*in_plane_v*bond_length;
+ return true;
+ }
+}}}
+
diff --git a/modules/qa/src/impl/reduced_impl.hh b/modules/qa/src/impl/reduced_impl.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e10ef2c380ac9eb5aee1093f9fd57a80256483b4
--- /dev/null
+++ b/modules/qa/src/impl/reduced_impl.hh
@@ -0,0 +1,53 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2010 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/mol.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)
+ { }
+
+ virtual bool VisitResidue(const mol::ResidueHandle& res);
+
+ 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:
+ ReducedStatOptions opts_;
+ mol::EntityHandle ent_;
+};
+
+}}}
+
+#endif
diff --git a/modules/qa/src/reduced_potential.cc b/modules/qa/src/reduced_potential.cc
new file mode 100644
index 0000000000000000000000000000000000000000..e2bbc2bbd3cf06aff0e945a9fb229c0545120ed6
--- /dev/null
+++ b/modules/qa/src/reduced_potential.cc
@@ -0,0 +1,135 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2010 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)
+ { }
+
+ 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_;
+};
+
+
+}
+
+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();
+}
+
+
+}}
diff --git a/modules/qa/src/reduced_potential.hh b/modules/qa/src/reduced_potential.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a85698ffe57753bfa30f0ca246a9a8269eabc815
--- /dev/null
+++ b/modules/qa/src/reduced_potential.hh
@@ -0,0 +1,73 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2010 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 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
new file mode 100644
index 0000000000000000000000000000000000000000..9c20a9cd20d3aa9f0c026417eae46479c4b4cdb5
--- /dev/null
+++ b/modules/qa/src/reduced_statistics.cc
@@ -0,0 +1,133 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2010 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)
+ { }
+
+ 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);
+}
+
+}}
\ No newline at end of file
diff --git a/modules/qa/src/reduced_statistics.hh b/modules/qa/src/reduced_statistics.hh
new file mode 100644
index 0000000000000000000000000000000000000000..419ac550f4b00c6ac216def800ba9860da295dbe
--- /dev/null
+++ b/modules/qa/src/reduced_statistics.hh
@@ -0,0 +1,120 @@
+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2010 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
+ void Extract(mol::EntityHandle 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