diff --git a/sidechain/pymod/CMakeLists.txt b/sidechain/pymod/CMakeLists.txt
index 5abd39150d025d662c9681cffab7f3fe6bb1a5a0..d08fe422596d6e81cdf570715e973360b5c05c0f 100644
--- a/sidechain/pymod/CMakeLists.txt
+++ b/sidechain/pymod/CMakeLists.txt
@@ -11,6 +11,7 @@ export_connector.cc
export_graph.cc
export_disulfid.cc
export_rotamer_density.cc
+export_rotamer_cruncher.cc
wrap_sidechain.cc
)
diff --git a/sidechain/pymod/export_rotamer_cruncher.cc b/sidechain/pymod/export_rotamer_cruncher.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c4c8e9a91a9ca303e4cae620f24ac829a4a5354a
--- /dev/null
+++ b/sidechain/pymod/export_rotamer_cruncher.cc
@@ -0,0 +1,41 @@
+#include <boost/python.hpp>
+#include <boost/python/register_ptr_to_python.hpp>
+#include <promod3/sidechain/rotamer_cruncher.hh>
+
+
+using namespace boost::python;
+using namespace promod3::sidechain;
+
+
+namespace{
+
+ void wrap_Initialize(RotamerCruncherPtr cruncher, const String& seqres,
+ const boost::python::list& weights){
+ std::vector<Real> v_weights;
+ for(int i = 0; i < boost::python::len(weights); ++i){
+ v_weights.push_back(boost::python::extract<Real>(weights[i]));
+ }
+ cruncher->Initialize(seqres,v_weights);
+ }
+ void wrap_UpdateRotamerProbabilities(RotamerCruncherPtr cruncher, boost::python::list rotamers,
+ uint seqres_pos, Real sigma){
+ std::vector<RotamerLibEntryPtr> v_rotamers;
+ for(int i = 0; i < boost::python::len(rotamers); ++i){
+ v_rotamers.push_back(boost::python::extract<RotamerLibEntryPtr>(rotamers[i]));
+ }
+ cruncher->UpdateRotamerProbabilities(v_rotamers,seqres_pos,sigma);
+ }
+}
+
+void export_RotamerCruncher()
+{
+ class_<RotamerCruncher>("RotamerCruncher", init<>())
+ .def("Initialize", &wrap_Initialize)
+ .def("AddStructure", &RotamerCruncher::AddStructure)
+ .def("MergeDensities", &RotamerCruncher::MergeDensities)
+ .def("GetExtractedProbability",&RotamerCruncher::GetExtractedProbability,(arg("seqres_pos"),arg("chi1"),arg("chi2")=std::numeric_limits<Real>::quiet_NaN()))
+ .def("UpdateRotamerProbabilities", &wrap_UpdateRotamerProbabilities)
+ ;
+
+ register_ptr_to_python<RotamerCruncherPtr>();
+}
diff --git a/sidechain/pymod/export_rotamer_density.cc b/sidechain/pymod/export_rotamer_density.cc
index 65cda779d623fad7aa7ce02c57b9a99e7f8c2026..86f692070a003c95598b0f8886e6af2e80147c50 100644
--- a/sidechain/pymod/export_rotamer_density.cc
+++ b/sidechain/pymod/export_rotamer_density.cc
@@ -43,6 +43,16 @@ namespace{
return RotamerDensity2DPtr(new RotamerDensity2D(num_bins));
}
+ Real wrap_GetPOne(RotamerDensity1DPtr p, Real chi1){
+ Real chi_angles[] = {chi1};
+ return p->GetP(chi_angles);
+ }
+
+ Real wrap_GetPTwo(RotamerDensity2DPtr p, Real chi1, Real chi2){
+ Real chi_angles[] = {chi1, chi2};
+ return p->GetP(chi_angles);
+ }
+
}
@@ -51,26 +61,39 @@ namespace{
void export_RotamerDensity()
{
- class_<RotamerDensity1D>("RotamerDensity1D", no_init)
+ class_<RotamerDensity, boost::noncopyable>("RotamerDensity", no_init);
+ register_ptr_to_python<RotamerDensityPtr>();
+
+ class_<RotamerDensity1D,bases<RotamerDensity> >("RotamerDensity1D", no_init)
.def("__init__",make_constructor(&WrapRotamerDensity1DInit))
.def("__init__",make_constructor(&WrapRotamerDensity1DListInit))
.def("Normalize",&RotamerDensity1D::Normalize)
.def("AddRotamer",&RotamerDensity1D::AddRotamer,(arg("rot")))
.def("GetNumRotamers",&RotamerDensity1D::GetNumRotamers)
.def("BinsPerDimension",&RotamerDensity1D::BinsPerDimension)
- .def("GetP",&RotamerDensity1D::GetP,(arg("chi1")))
+ .def("Bins",&RotamerDensity1D::Bins)
+ .def("ApplyFactor",&RotamerDensity1D::ApplyFactor)
+ .def("GetP",&wrap_GetPOne,(arg("chi1")))
+ .def("MergeProbabilities",&RotamerDensity1D::MergeProbabilities)
+ .def("Clear",&RotamerDensity1D::Clear)
+ //.def("Data",&RotamerDensity1D::Data)
;
register_ptr_to_python<RotamerDensity1DPtr>();
- class_<RotamerDensity2D>("RotamerDensity2D", no_init)
+ class_<RotamerDensity2D,bases<RotamerDensity> >("RotamerDensity2D", no_init)
.def("__init__",make_constructor(&WrapRotamerDensity2DInit))
.def("__init__",make_constructor(&WrapRotamerDensity2DListInit))
.def("Normalize",&RotamerDensity2D::Normalize)
.def("AddRotamer",&RotamerDensity2D::AddRotamer,(arg("rot")))
.def("GetNumRotamers",&RotamerDensity2D::GetNumRotamers)
.def("BinsPerDimension",&RotamerDensity2D::BinsPerDimension)
- .def("GetP",&RotamerDensity2D::GetP,(arg("chi1"),arg("chi2")))
+ .def("Bins",&RotamerDensity2D::Bins)
+ .def("ApplyFactor",&RotamerDensity2D::ApplyFactor)
+ .def("GetP",&wrap_GetPTwo,(arg("chi1"),arg("chi2")))
+ .def("MergeProbabilities",&RotamerDensity2D::MergeProbabilities)
+ .def("Clear",&RotamerDensity2D::Clear)
+ //.def("Data",&RotamerDensity2D::Data)
;
register_ptr_to_python<RotamerDensity2DPtr>();
diff --git a/sidechain/pymod/wrap_sidechain.cc b/sidechain/pymod/wrap_sidechain.cc
index eef99c233225fb6d169ee50af192ce56e07561f7..1b6263248a8cba799ed59189d72b6e2a785131c9 100644
--- a/sidechain/pymod/wrap_sidechain.cc
+++ b/sidechain/pymod/wrap_sidechain.cc
@@ -13,6 +13,7 @@ void export_Connector();
void export_Graph();
void export_Disulfid();
void export_RotamerDensity();
+void export_RotamerCruncher();
using namespace boost::python;
@@ -30,4 +31,6 @@ BOOST_PYTHON_MODULE(_sidechain)
export_Graph();
export_Disulfid();
export_RotamerDensity();
+ export_RotamerCruncher();
}
+
diff --git a/sidechain/src/CMakeLists.txt b/sidechain/src/CMakeLists.txt
index 68f5a732305f4d21432fbc55ac888e0e27bb7c4a..5a23a8382b123e3d1706594e2c8205f53eea578b 100644
--- a/sidechain/src/CMakeLists.txt
+++ b/sidechain/src/CMakeLists.txt
@@ -20,6 +20,7 @@ sidechain_connector.hh
settings.hh
disulfid.hh
rotamer_density.hh
+rotamer_cruncher.hh
)
set(SIDECHAIN_SOURCES
@@ -42,6 +43,7 @@ sidechain_connector.cc
disulfid.cc
rotamer_id.cc
rotamer_density.cc
+rotamer_cruncher.cc
)
module(NAME sidechain HEADERS ${SIDECHAIN_HEADERS} SOURCES ${SIDECHAIN_SOURCES} DEPENDS_ON
diff --git a/sidechain/src/rotamer_cruncher.cc b/sidechain/src/rotamer_cruncher.cc
new file mode 100644
index 0000000000000000000000000000000000000000..2e1947cd83c60072324b9f0b4b70d632e9ea8d46
--- /dev/null
+++ b/sidechain/src/rotamer_cruncher.cc
@@ -0,0 +1,416 @@
+
+#include <promod3/sidechain/rotamer_cruncher.hh>
+
+namespace{
+
+ Real _bin_size = 0.0872664625997;
+ Real _margin = 0.349065850399;
+ int _num_sampling_points = 9;
+
+}
+
+namespace promod3{ namespace sidechain{
+
+RotamerCruncher::RotamerCruncher():readonly_(false) { std::cout<<"construct rotamer cruncher "<<42<<std::endl;}
+
+void RotamerCruncher::Initialize(const String& seqres, const std::vector<Real>& c_weights){
+ if(readonly_){
+ throw promod3::Error("Cannot reinitialize readonly RotamerCruncher!");
+ }
+ std::cout<<"start with initialization"<<std::endl;
+ seqres_ = seqres;
+ c_weights_ = c_weights;
+ dimensionality_.clear();
+ extracted_rotamers_.clear();
+ extracted_rotamers_.resize(seqres_.size());
+ for(uint i = 0; i < seqres_.size(); ++i){
+ extracted_rotamers_[i].resize(c_weights_.size());
+ }
+
+ for(uint i = 0; i < seqres_.size(); ++i){
+ switch(seqres_[i]){
+ case 'S':{
+ dimensionality_.push_back(0);
+ break;
+ }
+ case 'V':{
+ dimensionality_.push_back(0);
+ break;
+ }
+ case 'T':{
+ dimensionality_.push_back(0);
+ break;
+ }
+ case 'C':{
+ dimensionality_.push_back(0);
+ break;
+ }
+ default:{
+ dimensionality_.push_back(1);
+ break;
+ }
+ }
+ }
+}
+
+void RotamerCruncher::Save(const String& filename){
+
+ std::ofstream out_stream(filename.c_str(), std::ios::binary);
+ uint32_t seqres_size = seqres_.size();
+ out_stream.write(reinterpret_cast<char*>(&seqres_size),4);
+ out_stream.write(reinterpret_cast<char*>(&seqres_[0]),seqres_size);
+ for(uint i = 0; i < merged_rot_densities_.size(); ++i){
+ uint32_t density_size = merged_rot_densities_[i]->Bins()*sizeof(double);
+ out_stream.write(reinterpret_cast<char*>(&density_size),4);
+ double* merged_density_data = merged_rot_densities_[i]->Data();
+ out_stream.write(reinterpret_cast<char*>(merged_density_data),density_size);
+ }
+
+
+}
+
+RotamerCruncherPtr RotamerCruncher::Load(const String& filename){
+
+ std::ifstream in_stream(filename.c_str(), std::ios::binary);
+ uint32_t seqres_size;
+ in_stream.read(reinterpret_cast<char*>(&seqres_size),4);
+ String seqres(seqres_size,'X');
+ in_stream.read(reinterpret_cast<char*>(&seqres[0]),seqres_size);
+
+ RotamerCruncherPtr loaded_cruncher(new RotamerCruncher);
+ std::vector<Real> empty_vec;
+ loaded_cruncher->Initialize(seqres,empty_vec);
+ for(uint i = 0; i < seqres_size; ++i){
+ uint32_t density_size;
+ in_stream.read(reinterpret_cast<char*>(&density_size),4);
+ if(loaded_cruncher->dimensionality_[i] == 0){
+ loaded_cruncher->merged_rot_densities_.push_back(RotamerDensityPtr(new RotamerDensity1D(density_size)));
+ }
+ else{
+ loaded_cruncher->merged_rot_densities_.push_back(RotamerDensityPtr(new RotamerDensity2D(sqrt(density_size))));
+ }
+ double* merged_density_data = loaded_cruncher->merged_rot_densities_.back()->Data();
+ in_stream.read(reinterpret_cast<char*>(merged_density_data),density_size);
+ }
+ loaded_cruncher->readonly_ = true;
+
+
+
+ return loaded_cruncher;
+}
+
+
+void RotamerCruncher::AddStructure(const ost::mol::ChainHandle& chain, const ost::seq::AlignmentHandle& aln, uint c){
+ if(readonly_){
+ throw promod3::Error("Cannot add structure to readonly RotamerCruncher!");
+ }
+
+ if(c >= c_weights_.size()){
+ throw promod3::Error("Invalid cluster index!");
+ }
+ int seqres_pos = 0;
+ ost::seq::ConstSequenceHandle seqres_seq = aln.GetSequence(0);
+ if(seqres_seq.GetGaplessString() != seqres_){
+ throw promod3::Error("Seqres and target sequence from the alignment are not the same!");
+ }
+ ost::seq::ConstSequenceHandle templ_seq = aln.GetSequence(1);
+ ost::mol::ResidueHandleList res_list = chain.GetResidueList();
+ RotamerLibEntryPtr rotamer;
+
+ // build naive cluster densities
+ for(int i = 0; i < aln.GetLength(); ++i){
+ if (seqres_seq[i] == '-'){
+ continue;
+ }
+ if (templ_seq[i] == '-'){
+ ++seqres_pos;
+ continue;
+ }
+ if (seqres_seq[i] == 'A' || seqres_seq[i] == 'G'){
+ ++seqres_pos;
+ continue;
+ }
+ if(seqres_seq[i] == templ_seq[i]){
+ int residue_index = templ_seq.GetResidueIndex(i);
+
+ ost::mol::ResidueHandle residue = res_list[residue_index];
+ try{
+ rotamer = RotamerLibEntry::FromResidue(residue);
+ } catch(promod3::Error& e){
+ std::cerr<<"just caught an error"<<std::endl;
+ ++seqres_pos;
+ continue;
+ }
+ extracted_rotamers_[seqres_pos][c].push_back(*rotamer);
+ }
+ // include similar sidechains
+ else if(seqres_seq[i] == 'N' && templ_seq[i] == 'D'){
+ int residue_index = templ_seq.GetResidueIndex(i);
+ ost::mol::ResidueHandle residue = res_list[residue_index];
+ try{
+ rotamer = RotamerLibEntry::FromResidue(residue);
+ } catch(promod3::Error& e){
+ ++seqres_pos;
+ continue;
+ }
+ extracted_rotamers_[seqres_pos][c].push_back(*rotamer);
+ // not sure if this works ...
+ rotamer->chi2 -= M_PI;
+ extracted_rotamers_[seqres_pos][c].push_back(*rotamer);
+ }
+ else if(seqres_seq[i] == 'D' && templ_seq[i] == 'N'){
+ int residue_index = templ_seq.GetResidueIndex(i);
+ ost::mol::ResidueHandle residue = res_list[residue_index];
+ try{
+ rotamer = RotamerLibEntry::FromResidue(residue);
+ } catch(promod3::Error& e){
+ ++seqres_pos;
+ continue;
+ }
+ extracted_rotamers_[seqres_pos][c].push_back(*rotamer);
+ }
+ else if(seqres_seq[i] == 'R' && templ_seq[i] == 'K'){
+ int residue_index = templ_seq.GetResidueIndex(i);
+ ost::mol::ResidueHandle residue = res_list[residue_index];
+ try{
+ rotamer = RotamerLibEntry::FromResidue(residue);
+ } catch(promod3::Error& e){
+ ++seqres_pos;
+ continue;
+ }
+ extracted_rotamers_[seqres_pos][c].push_back(*rotamer);
+ }
+ else if(seqres_seq[i] == 'K' && templ_seq[i] == 'R'){
+
+ int residue_index = templ_seq.GetResidueIndex(i);
+ ost::mol::ResidueHandle residue = res_list[residue_index];
+ try{
+ rotamer = RotamerLibEntry::FromResidue(residue);
+ } catch(promod3::Error& e){
+ ++seqres_pos;
+ continue;
+ }
+ extracted_rotamers_[seqres_pos][c].push_back(*rotamer);
+ }
+ else if(seqres_seq[i] == 'S' && templ_seq[i] == 'T'){
+ int residue_index = templ_seq.GetResidueIndex(i);
+ ost::mol::ResidueHandle residue = res_list[residue_index];
+ try{
+ rotamer = RotamerLibEntry::FromResidue(residue);
+ } catch(promod3::Error& e){
+ ++seqres_pos;
+ continue;
+ }
+ extracted_rotamers_[seqres_pos][c].push_back(*rotamer);
+ }
+ else if(seqres_seq[i] == 'T' && templ_seq[i] == 'S'){
+ int residue_index = templ_seq.GetResidueIndex(i);
+ ost::mol::ResidueHandle residue = res_list[residue_index];
+ try{
+ rotamer = RotamerLibEntry::FromResidue(residue);
+ } catch(promod3::Error& e){
+ ++seqres_pos;
+ continue;
+ }
+ extracted_rotamers_[seqres_pos][c].push_back(*rotamer);
+ }
+ else if(seqres_seq[i] == 'Y' && templ_seq[i] == 'F'){
+ int residue_index = templ_seq.GetResidueIndex(i);
+ ost::mol::ResidueHandle residue = res_list[residue_index];
+ try{
+ rotamer = RotamerLibEntry::FromResidue(residue);
+ } catch(promod3::Error& e){
+ ++seqres_pos;
+ continue;
+ }
+ extracted_rotamers_[seqres_pos][c].push_back(*rotamer);
+ }
+ else if(seqres_seq[i] == 'F' && templ_seq[i] == 'Y'){
+ int residue_index = templ_seq.GetResidueIndex(i);
+ ost::mol::ResidueHandle residue = res_list[residue_index];
+ try{
+ rotamer = RotamerLibEntry::FromResidue(residue);
+ } catch(promod3::Error& e){
+ ++seqres_pos;
+ continue;
+ }
+ extracted_rotamers_[seqres_pos][c].push_back(*rotamer);
+ }
+ else if(seqres_seq[i] == 'E' && templ_seq[i] == 'Q'){
+ int residue_index = templ_seq.GetResidueIndex(i);
+ ost::mol::ResidueHandle residue = res_list[residue_index];
+ try{
+ rotamer = RotamerLibEntry::FromResidue(residue);
+ } catch(promod3::Error& e){
+ ++seqres_pos;
+ continue;
+ }
+ extracted_rotamers_[seqres_pos][c].push_back(*rotamer);
+ }
+ else if(seqres_seq[i] == 'Q' && templ_seq[i] == 'E'){
+ int residue_index = templ_seq.GetResidueIndex(i);
+ ost::mol::ResidueHandle residue = res_list[residue_index];
+ try{
+ rotamer = RotamerLibEntry::FromResidue(residue);
+ } catch(promod3::Error& e){
+ ++seqres_pos;
+ continue;
+ }
+ extracted_rotamers_[seqres_pos][c].push_back(*rotamer);
+ }
+
+
+ ++seqres_pos;
+ }
+}
+void RotamerCruncher::MergeDensities(){
+ // normalize and apply weight to cluster densities
+ if(readonly_){
+ throw promod3::Error("Cannot merge densities of readonly RotamerCruncher!");
+ }
+ merged_rot_densities_.clear();
+
+ for(uint i = 0; i < extracted_rotamers_.size(); ++i){
+ if(dimensionality_[i] == 0){
+ merged_rot_densities_.push_back(RotamerDensityPtr(new RotamerDensity1D(72)));
+ }
+ else if(dimensionality_[i] == 1){
+ merged_rot_densities_.push_back(RotamerDensityPtr(new RotamerDensity2D(72)));
+ }
+ }
+
+ std::cerr<<"merging probs"<<std::endl;
+ RotamerLibEntry actual_rotamer;
+ for(uint i = 0; i < extracted_rotamers_.size(); ++i){
+ for(uint j = 0; j < extracted_rotamers_[i].size(); ++j){
+ Real p = c_weights_[j] / extracted_rotamers_[i][j].size();
+ for(uint k = 0; k < extracted_rotamers_[i][j].size(); ++k){
+ actual_rotamer = extracted_rotamers_[i][j][k];
+ actual_rotamer.probability = p;
+ merged_rot_densities_[i]->AddRotamer(actual_rotamer);
+ }
+ }
+ }
+
+ for(uint i = 0; i < extracted_rotamers_.size(); ++i){
+ merged_rot_densities_[i]->Normalize();
+ }
+}
+
+Real RotamerCruncher::GetExtractedProbability(uint seqres_pos, Real chi1, Real chi2){
+
+
+ if(seqres_pos >= seqres_.size()){
+ throw promod3::Error("Invalid index observed!");
+ }
+
+ Real chi_angles[2];
+ chi_angles[0] = chi1;
+
+ if(dimensionality_[seqres_pos] == 0){
+ return merged_rot_densities_[seqres_pos]->GetP(&chi_angles[0]);
+ }
+
+ if(chi2 != chi2){
+ throw promod3::Error("Chi2 angle must be given!");
+ }
+
+ chi_angles[1] = chi2;
+
+ return merged_rot_densities_[seqres_pos]->GetP(&chi_angles[0]);
+}
+
+void RotamerCruncher::UpdateRotamerProbabilities(std::vector<RotamerLibEntryPtr>& rotamers, uint seqres_pos, Real sigma){
+
+ std::cerr<<"start to update probs"<<std::endl;
+
+ String seqres = seqres_;
+ if(seqres_pos >= seqres.size()){
+ throw Error("Seqres_pos out of range!");
+ }
+
+ int N = merged_rot_densities_[seqres_pos]->GetNumRotamers();
+ Real sippl_prefac_one = (1.0/(N*0.05+1));
+ Real sippl_prefac_two = (N*0.05/(N*0.05+1));
+ Real chi_angles[2];
+ std::vector<Real> density_probs;
+ std::vector<Real> original_probs;
+
+ //in case of only chi1
+ if(dimensionality_[seqres_pos]==0){
+
+ for(std::vector<RotamerLibEntryPtr>::iterator i = rotamers.begin();
+ i != rotamers.end(); ++i){
+
+ Real chi1 = (*i)->chi1;
+ Real summed_prob = 0.0;
+
+ chi_angles[0] = chi1-_margin;
+
+ for(int j = 0; j < _num_sampling_points; ++j){
+ summed_prob += merged_rot_densities_[seqres_pos]->GetP(&chi_angles[0]);
+ chi_angles[0] += _bin_size;
+ }
+
+ density_probs.push_back(summed_prob);
+ original_probs.push_back((*i)->probability);
+ }
+ }
+
+ else{
+ for(std::vector<RotamerLibEntryPtr>::iterator i = rotamers.begin();
+ i != rotamers.end(); ++i){
+
+ Real chi1 = (*i)->chi1;
+ Real chi2 = (*i)->chi2;
+ Real summed_prob = 0.0;
+
+ chi_angles[0] = chi1-_margin;
+
+
+
+ for(int j = 0; j<_num_sampling_points; ++j){
+ chi_angles[1] = chi2-_margin;
+ for(int k = 0; k<_num_sampling_points; ++k){
+ summed_prob += merged_rot_densities_[seqres_pos]->GetP(&chi_angles[0]);
+ chi_angles[1] += _bin_size;
+ }
+ chi_angles[0] += _bin_size;
+ }
+
+ density_probs.push_back(summed_prob);
+ original_probs.push_back((*i)->probability);
+ }
+ }
+
+ Real sum = 0.0;
+
+ for(std::vector<Real>::iterator i = density_probs.begin(); i != density_probs.end(); ++i){
+ sum += (*i);
+ }
+
+ if(sum == 0.0) return;
+
+ for(uint i = 0; i < density_probs.size(); ++i){
+ density_probs[i] = density_probs[i]/sum;
+ }
+
+ std::vector<Real> sippl_probs;
+
+ for(uint i = 0; i < density_probs.size(); ++i){
+ sippl_probs.push_back(sippl_prefac_one*original_probs[i]+sippl_prefac_two*density_probs[i]);
+ }
+
+ sum = 0.0;
+ for(std::vector<Real>::iterator i = sippl_probs.begin(); i != sippl_probs.end(); ++i){
+ sum += (*i);
+ }
+
+ for(uint i = 0; i < sippl_probs.size(); ++i){
+ std::cerr<<"2D probability before updating "<<rotamers[i]->probability<<std::endl;
+ rotamers[i]->probability = sippl_probs[i]/sum;
+ std::cerr<<"2D probability after updating "<<rotamers[i]->probability<<std::endl;
+ }
+}
+
+}}//ns
diff --git a/sidechain/src/rotamer_cruncher.hh b/sidechain/src/rotamer_cruncher.hh
new file mode 100644
index 0000000000000000000000000000000000000000..dcf0453196463425d0db8d69c3f81bd735431be5
--- /dev/null
+++ b/sidechain/src/rotamer_cruncher.hh
@@ -0,0 +1,57 @@
+#ifndef PROMOD3_ROTAMER_CRUNCHER_HH
+#define PROMOD3_ROTAMER_CRUNCHER_HH
+
+#include <math.h>
+#include <vector>
+#include <limits.h>
+#include <fstream>
+
+#include <boost/shared_ptr.hpp>
+
+#include <ost/seq/alignment_handle.hh>
+#include <ost/mol/entity_handle.hh>
+#include <promod3/sidechain/tetrahedral_polytope.hh>
+#include <promod3/sidechain/pairwise_terms.hh>
+#include <promod3/sidechain/rotamer_lib_entry.hh>
+#include <promod3/sidechain/rotamer_density.hh>
+
+
+namespace promod3{ namespace sidechain{
+
+class RotamerCruncher;
+typedef boost::shared_ptr<RotamerCruncher> RotamerCruncherPtr;
+
+
+class RotamerCruncher{
+
+public:
+
+ RotamerCruncher();
+
+ void Save(const String& filename);
+
+ static RotamerCruncherPtr Load(const String& filename);
+
+ void Initialize(const String& seqres, const std::vector<Real>& c_weights);
+
+ void AddStructure(const ost::mol::ChainHandle& chain, const ost::seq::AlignmentHandle& aln, uint c);
+
+ void MergeDensities();
+
+ Real GetExtractedProbability(uint seqres_pos, Real chi1, Real chi2 = std::numeric_limits<Real>::quiet_NaN());
+
+ void UpdateRotamerProbabilities(std::vector<RotamerLibEntryPtr>& rotamers, uint seqres_pos, Real sigma);
+
+private:
+
+ String seqres_;
+ std::vector<Real> c_weights_;
+ std::vector<int> dimensionality_;
+ std::vector<std::vector<std::vector<RotamerLibEntry> > > extracted_rotamers_;
+ std::vector<RotamerDensityPtr> merged_rot_densities_;
+ bool readonly_;
+};
+
+}}//ns
+
+#endif
diff --git a/sidechain/src/rotamer_density.cc b/sidechain/src/rotamer_density.cc
index ce5752e309298cf7df93baa6da1833eb73df7ced..3806376abfba6898c04a77cc309c64478b1d778d 100644
--- a/sidechain/src/rotamer_density.cc
+++ b/sidechain/src/rotamer_density.cc
@@ -3,8 +3,7 @@
namespace promod3{ namespace sidechain{
RotamerDensity1D::RotamerDensity1D(uint num_bins):num_rotamers_(0),
- num_bins_(num_bins),
- normalized_(false){
+ num_bins_(num_bins){
if(num_bins_ > 1000){
std::stringstream ss;
@@ -19,8 +18,7 @@ RotamerDensity1D::RotamerDensity1D(uint num_bins):num_rotamers_(0),
}
RotamerDensity1D::RotamerDensity1D(const std::vector<RotamerDensity1DPtr>& densities,
- const std::vector<Real>& weights): num_rotamers_(0),
- normalized_(false){
+ const std::vector<Real>& weights): num_rotamers_(0){
if(densities.empty()){
std::stringstream ss;
@@ -65,17 +63,18 @@ RotamerDensity1D::RotamerDensity1D(const std::vector<RotamerDensity1DPtr>& densi
void RotamerDensity1D::Normalize(){
//sum everything up
- Real sum = 0.0;
+ float sum = 0.0;
for(uint i = 0; i < num_bins_; ++i){
sum += data_[i];
}
//and normalize
- Real factor = 1.0/sum;
- for(uint i = 0; i < num_bins_; ++i){
- data_[i] *= factor;
+ if(sum != 0.0){
+ float factor = 1.0/sum;
+ for(uint i = 0; i < num_bins_; ++i){
+ data_[i] *= factor;
+ }
}
- normalized_ = true;
}
void RotamerDensity1D::AddRotamer(const RotamerLibEntry& rot){
@@ -131,7 +130,15 @@ void RotamerDensity1D::AddRotamer(const RotamerLibEntry& rot){
++num_rotamers_;
}
-Real RotamerDensity1D::GetP(Real chi1){
+void RotamerDensity1D::ApplyFactor(Real factor){
+ for(uint i = 0; i < num_bins_; ++i){
+ data_[i] *= factor;
+ }
+}
+
+Real RotamerDensity1D::GetP(Real* chi_angles){
+
+ Real chi1 = chi_angles[0];
while(chi1 < -M_PI) chi1 += 2*M_PI;
while(chi1 >= M_PI) chi1 -= 2*M_PI;
@@ -140,9 +147,28 @@ Real RotamerDensity1D::GetP(Real chi1){
return data_[bin];
}
+void RotamerDensity1D::MergeProbabilities(RotamerDensityPtr other){
+
+ if(this->Bins() != other->Bins()){
+ throw promod3::Error("Inconsistent size when adding Rotamer Densities!");
+ }
+
+ double* other_data = other->Data();
+ uint other_num_rotamers = other->GetNumRotamers();
+ num_rotamers_ += other_num_rotamers;
+ for(uint i = 0; i < num_bins_; ++i){
+ data_[i] += other_data[i];
+ }
+
+}
+
+void RotamerDensity1D::Clear(){
+ num_rotamers_ = 0;
+ memset(data_,0,num_bins_*sizeof(double));
+}
+
RotamerDensity2D::RotamerDensity2D(uint num_bins):num_rotamers_(0),
- num_bins_(num_bins),
- normalized_(false){
+ num_bins_(num_bins){
if(num_bins > 1000){
std::stringstream ss;
@@ -158,8 +184,7 @@ RotamerDensity2D::RotamerDensity2D(uint num_bins):num_rotamers_(0),
}
RotamerDensity2D::RotamerDensity2D(const std::vector<RotamerDensity2DPtr>& densities,
- const std::vector<Real>& weights): num_rotamers_(0),
- normalized_(false){
+ const std::vector<Real>& weights): num_rotamers_(0){
if(densities.empty()){
std::stringstream ss;
@@ -206,9 +231,11 @@ void RotamerDensity2D::Normalize(){
for(uint i = 0; i < num_bins_*num_bins_; ++i){
sum += data_[i];
}
- double factor = 1.0/sum;
- for(uint i = 0; i < num_bins_*num_bins_; ++i){
- data_[i] *= factor;
+ if(sum != 0.0){
+ double factor = 1.0/sum;
+ for(uint i = 0; i < num_bins_*num_bins_; ++i){
+ data_[i] *= factor;
+ }
}
}
@@ -287,7 +314,15 @@ void RotamerDensity2D::AddRotamer(const RotamerLibEntry& rot){
++num_rotamers_;
}
-Real RotamerDensity2D::GetP(Real chi1, Real chi2){
+void RotamerDensity2D::ApplyFactor(Real factor){
+ for(uint i = 0; i < num_bins_*num_bins_; ++i){
+ data_[i] *= factor;
+ }
+}
+
+Real RotamerDensity2D::GetP(Real* chi_angles){
+ Real chi1 = chi_angles[0];
+ Real chi2 = chi_angles[1];
while(chi1 < -M_PI) chi1 += 2*M_PI;
while(chi1 >= M_PI) chi1 -= 2*M_PI;
while(chi2 < -M_PI) chi2 += 2*M_PI;
@@ -300,4 +335,24 @@ Real RotamerDensity2D::GetP(Real chi1, Real chi2){
return data_[bin];
}
+void RotamerDensity2D::MergeProbabilities(RotamerDensityPtr other){
+
+ if(this->Bins() != other->Bins()){
+ throw promod3::Error("Inconsistent size when adding Rotamer Densities!");
+ }
+
+ double* other_data = other->Data();
+ uint other_num_rotamers = other->GetNumRotamers();
+ num_rotamers_ += other_num_rotamers;
+ for(uint i = 0; i < num_bins_*num_bins_; ++i){
+ data_[i] += other_data[i];
+ }
+}
+
+void RotamerDensity2D::Clear(){
+ num_rotamers_ = 0;
+ memset(data_,0,num_bins_*num_bins_*sizeof(double));
+}
+
}} //ns
+
diff --git a/sidechain/src/rotamer_density.hh b/sidechain/src/rotamer_density.hh
index 67ce145539bb475064d8b5a80467e481d585cbf3..7e024cb471c447911a43be193023be5654295d32 100644
--- a/sidechain/src/rotamer_density.hh
+++ b/sidechain/src/rotamer_density.hh
@@ -10,11 +10,44 @@ namespace promod3{ namespace sidechain{
class RotamerDensity1D;
class RotamerDensity2D;
+class RotamerDensity;
typedef boost::shared_ptr<RotamerDensity1D> RotamerDensity1DPtr;
typedef boost::shared_ptr<RotamerDensity2D> RotamerDensity2DPtr;
+typedef boost::shared_ptr<RotamerDensity> RotamerDensityPtr;
-class RotamerDensity1D{
+
+class RotamerDensity{
+public:
+
+ RotamerDensity() { }
+
+ virtual ~RotamerDensity() { };
+
+ virtual void Normalize() = 0;
+
+ virtual void AddRotamer(const RotamerLibEntry& rot) = 0;
+
+ virtual uint GetNumRotamers() = 0;
+
+ virtual uint BinsPerDimension() = 0;
+
+ virtual uint Bins() = 0;
+
+ virtual void ApplyFactor(Real factor) = 0;
+
+ virtual Real GetP(Real* chi_angles) = 0;
+
+ virtual void MergeProbabilities(RotamerDensityPtr other) = 0;
+
+ virtual void Clear() = 0;
+
+ virtual double* Data() = 0;
+
+};
+
+
+class RotamerDensity1D : public RotamerDensity{
public:
RotamerDensity1D(uint num_bins);
@@ -22,17 +55,28 @@ public:
RotamerDensity1D(const std::vector<RotamerDensity1DPtr>& densities,
const std::vector<Real>& weights);
- ~RotamerDensity1D(){ delete [] data_; }
+ virtual ~RotamerDensity1D(){ delete [] data_; }
+
+ virtual void Normalize();
+
+ virtual void AddRotamer(const RotamerLibEntry& rot);
+
+ virtual uint GetNumRotamers() { return num_rotamers_;}
+
+ virtual uint BinsPerDimension() { return num_bins_; }
+
+ virtual uint Bins() { return num_bins_;}
+
+ virtual void ApplyFactor(Real factor);
- void Normalize();
+ virtual Real GetP(Real* chi_angles);
- void AddRotamer(const RotamerLibEntry& rot);
+ virtual void MergeProbabilities(RotamerDensityPtr other);
- uint GetNumRotamers() { return num_rotamers_;}
+ virtual void Clear();
- uint BinsPerDimension() { return num_bins_; }
+ virtual double* Data() { return data_; }
- Real GetP(Real chi1);
private:
@@ -40,11 +84,10 @@ private:
uint num_rotamers_;
Real bin_size_;
uint num_bins_;
- bool normalized_;
};
-class RotamerDensity2D{
+class RotamerDensity2D : public RotamerDensity{
public:
RotamerDensity2D(uint num_bins);
@@ -52,17 +95,27 @@ public:
RotamerDensity2D(const std::vector<RotamerDensity2DPtr>& densities,
const std::vector<Real>& weights);
- ~RotamerDensity2D(){ delete data_; }
+ virtual ~RotamerDensity2D(){ delete [] data_; }
- void Normalize();
+ virtual void Normalize();
- void AddRotamer(const RotamerLibEntry& rot);
+ virtual void AddRotamer(const RotamerLibEntry& rot);
- uint GetNumRotamers() { return num_rotamers_;}
+ virtual uint GetNumRotamers() { return num_rotamers_;}
- uint BinsPerDimension() { return num_bins_; }
+ virtual uint BinsPerDimension() { return num_bins_; }
- Real GetP(Real chi1, Real chi2);
+ virtual uint Bins() { return num_bins_*num_bins_;}
+
+ virtual void ApplyFactor(Real factor);
+
+ virtual Real GetP(Real* chi_angles);
+
+ virtual void MergeProbabilities(RotamerDensityPtr other);
+
+ virtual void Clear();
+
+ virtual double* Data() { return data_; }
private:
@@ -70,10 +123,9 @@ private:
uint num_rotamers_;
Real bin_size_;
uint num_bins_;
- bool normalized_;
};
-
}} //ns
#endif
+