diff --git a/modelling/pymod/_closegaps.py b/modelling/pymod/_closegaps.py
index 9523221847481551a8fd706476a7cb708f6eb00d..a94e1c180498e129f7cf1bfb025121c169000d37 100644
--- a/modelling/pymod/_closegaps.py
+++ b/modelling/pymod/_closegaps.py
@@ -33,7 +33,7 @@ def _GetBestLC(mhandle, loop_candidates, start_resnum, chain_idx,
all_scores = db_scores.Copy() # copy to avoid changing db_scores
# add backbone scores
- loop_candidates.CalculateBackboneScores(all_scores, mhandle.backbone_scorer,
+ loop_candidates.CalculateBackboneScores(all_scores, mhandle,
start_resnum, chain_idx)
# add all atom scores?
@@ -385,7 +385,12 @@ def CloseSmallDeletions(mhandle, max_extension=9, clash_thresh=1.0,
# check/setup scoring
if not IsBackboneScoringSetUp(mhandle):
SetupDefaultBackboneScoring(mhandle)
- clash_scorer = mhandle.backbone_scorer.Get("clash")
+
+ # we're only calculating clash scores here, so we copy the default
+ # env and only attach a clash scorer.
+ clash_scorer_env = mhandle.backbone_scorer_env.Copy()
+ clash_scorer = scoring.ClashScorer()
+ clash_scorer.AttachEnvironment(clash_scorer_env)
if ff_lookup == None:
ff_lookup = loop.ForcefieldLookup.GetDefault()
@@ -444,9 +449,18 @@ def CloseSmallDeletions(mhandle, max_extension=9, clash_thresh=1.0,
bb_relaxer = BackboneRelaxer(bb_list, ff_lookup)
potential_e = bb_relaxer.Run(bb_list, 300, 0.1)
# check for clashes
- score = clash_scorer.CalculateScore(bb_list,
- n_stem_resnum.GetNum(),
+
+ # its a bit problematic since we score loops that are shifting
+ # around... so we perform a stash operation for each scoring step
+
+ clash_scorer_env.Stash(n_stem_resnum.GetNum(),
+ len(bb_list), current_chain_index)
+ clash_scorer_env.SetEnvironment(bb_list, n_stem_resnum.GetNum(),
+ current_chain_index)
+ score = clash_scorer.CalculateScore(n_stem_resnum.GetNum(),
+ len(bb_list),
current_chain_index)
+ clash_scorer_env.Pop()
# if there is no clash and potential energy is low enough we
# just solved a gap, delete it and update the scorer for the
@@ -457,9 +471,13 @@ def CloseSmallDeletions(mhandle, max_extension=9, clash_thresh=1.0,
ost.LogInfo("Closed: %s by relaxing %s" % \
(mhandle.gaps[current_gap_index], current_gap))
InsertLoopClearGaps(mhandle, bb_list, current_gap)
+ clash_scorer_env.SetEnvironment(bb_list,
+ n_stem_resnum.GetNum(),
+ current_chain_index)
success = True
break
+
# On closed gap, it is removed so the no. of gaps goes down by itself.
# In case of no success, counter needs to be increased to jump to the
# next gap.
@@ -775,13 +793,15 @@ def FillLoopsByDatabase(mhandle, fragment_db, structure_db,
# remove clashing ones if desired
if clash_thresh > 0:
- scorer = mhandle.backbone_scorer.Get("clash")
+ clash_score_container = ScoreContainer()
+ candidates.CalculateBackboneScores(clash_score_container,
+ mhandle, ["clash"],
+ n_stem.GetNumber().GetNum(),
+ actual_gap.GetChainIndex())
+ clash_scores = clash_score_container.Get("clash")
to_keep = []
for i, bb_list in enumerate(candidates):
- score = scorer.CalculateScore(bb_list,
- n_stem.GetNumber().GetNum(),
- actual_gap.GetChainIndex())
- if score < clash_thresh:
+ if clash_scores[i] < clash_thresh:
to_keep.append(i)
candidates = candidates.Extract(to_keep)
orig_indices = [orig_indices[i] for i in to_keep]
@@ -1025,7 +1045,9 @@ def FillLoopsByMonteCarlo(mhandle, torsion_sampler, max_loops_to_search=6,
weights = _GetMCWeights()
mc_scorer = LinearScorer(mhandle.backbone_scorer,
+ mhandle.backbone_scorer_env,
actual_gap.before.GetNumber().GetNum(),
+ len(actual_gap.full_seq),
actual_chain_idx, weights)
start_temperature = 100
@@ -1047,7 +1069,8 @@ def FillLoopsByMonteCarlo(mhandle, torsion_sampler, max_loops_to_search=6,
candidates = LoopCandidates.FillFromMonteCarloSampler(
actual_gap.full_seq, mc_num_loops, mc_steps, mc_sampler,
mc_closer, mc_scorer, mc_cooler)
- except RuntimeError:
+ #except RuntimeError:
+ except Exception, e:
# monte carlo cannot be initialized when the stems are too far
# apart => we need further loop extension
ost.LogVerbose("Failed in setting up " + str(actual_gap) +
@@ -1274,7 +1297,9 @@ def ModelTermini(mhandle, torsion_sampler, fragger_handles=None,
# setup scorer
weights = _GetMCWeights()
- mc_scorer = LinearScorer(mhandle.backbone_scorer, start_resnum,
+ mc_scorer = LinearScorer(mhandle.backbone_scorer,
+ mhandle.backbone_scorer_env,
+ start_resnum, len(actual_gap.full_seq),
actual_chain_idx, weights)
# setup cooler
@@ -1295,8 +1320,7 @@ def ModelTermini(mhandle, torsion_sampler, fragger_handles=None,
if len(candidates) > 0:
# score candidates
bb_scores = ScoreContainer()
- candidates.CalculateBackboneScores(bb_scores,
- mhandle.backbone_scorer,
+ candidates.CalculateBackboneScores(bb_scores, mhandle,
start_resnum, actual_chain_idx)
scores = bb_scores.LinearCombine(ScoringWeights.GetWeights())
min_score = min(scores)
@@ -1480,8 +1504,10 @@ def CloseLargeDeletions(mhandle, structure_db, linker_length=8,
closer.Close(bb_list, bb_list)
# a simple score gets calculated and best chosen
+ mhandle.backbone_scorer_env.SetEnvironment(bb_list, first_num,
+ actual_chain_idx)
score = mhandle.backbone_scorer.CalculateLinearCombination(\
- scorer_weights, bb_list, first_num, actual_chain_idx)
+ scorer_weights, first_num, len(bb_list), actual_chain_idx)
if score < best_score:
best_score = score
diff --git a/modelling/pymod/_denovo.py b/modelling/pymod/_denovo.py
index 94a06c8ca623df214a99c9f88e3cf8ab474a2292..51dd26d0ca0b540c9b6523c66337f5aa0a96fa52 100644
--- a/modelling/pymod/_denovo.py
+++ b/modelling/pymod/_denovo.py
@@ -8,6 +8,7 @@ def GenerateDeNovoTrajectories(sequence,
psipred_prediction = None,
fragment_handler = None,
scorer = None,
+ scorer_env = None,
scoring_weights = None):
'''Example de novo modelling pipeline based on Fragment sampling and
backbone scoring. Take this as a starting point for more advanced
@@ -47,6 +48,10 @@ def GenerateDeNovoTrajectories(sequence,
torsion and pairwise
:type scorer: :class:`promod3.scoring.BackboneOverallScorer`
+ :param scorer_env: The scoring env that relates to **scorer**
+ This environment will be changed!
+ :type scorer_env: :class:`promod3.scoring.BackboneScoreEnv`
+
:param scoring_weights: Linear weights for different scores. If not provided,
the output of ScoringWeights.GetWeights() is used.
Please note, that the weights must be consistent
@@ -71,6 +76,9 @@ def GenerateDeNovoTrajectories(sequence,
err = "psipred_prediction must be consistent with the input sequence!"
raise RuntimeError(err)
+ if scorer != None and scorer_env == None:
+ raise RuntimeError("Must provide an environment when you provide a scorer!")
+
if fragment_handler == None:
# we first have to build our own handler
fragment_handler = modelling.FraggerHandle(sequence,
@@ -104,7 +112,8 @@ def GenerateDeNovoTrajectories(sequence,
mc_sampler = FragmentSampler(sequence, fragger_list,
init_fragments = 5)
mc_closer = DeNovoCloser()
- mc_scorer = LinearScorer(scorer, 1, 0, scoring_weights)
+ mc_scorer = LinearScorer(scorer, scorer_env, 1, len(sequence), 0,
+ scoring_weights)
mc_cooler = ExponentialCooler(start_temperature,
cooling_interval,
cooling_factor)
diff --git a/modelling/pymod/export_loop_candidate.cc b/modelling/pymod/export_loop_candidate.cc
index 3ef94b937e1d04c803c0b87c1edac864dd69e1a5..939c3a170e7fb62695cada011bce8f504b0d0a73 100644
--- a/modelling/pymod/export_loop_candidate.cc
+++ b/modelling/pymod/export_loop_candidate.cc
@@ -135,18 +135,18 @@ boost::python::list WrapGetClusteredCandidates(LoopCandidatesPtr p,
void WrapCalculateBackboneScores(LoopCandidatesPtr p,
ScoreContainer& score_container,
- scoring::BackboneOverallScorerPtr scorer,
+ const ModellingHandle& mhandle,
uint start_resnum, uint chain_idx) {
- p->CalculateBackboneScores(score_container, scorer, start_resnum, chain_idx);
+ p->CalculateBackboneScores(score_container, mhandle, start_resnum, chain_idx);
}
void WrapCalculateBackboneScoresK(LoopCandidatesPtr p,
ScoreContainer& score_container,
- scoring::BackboneOverallScorerPtr scorer,
+ const ModellingHandle& mhandle,
const boost::python::list& keys,
uint start_resnum, uint chain_idx) {
std::vector<String> v_keys;
core::ConvertListToVector(keys, v_keys);
- p->CalculateBackboneScores(score_container, scorer, v_keys, start_resnum,
+ p->CalculateBackboneScores(score_container, mhandle, v_keys, start_resnum,
chain_idx);
}
void WrapCalculateAllAtomScores(LoopCandidatesPtr p,
@@ -268,10 +268,10 @@ void export_loop_candidate() {
(arg("n_stem"), arg("c_stem"), arg("pivot_one"), arg("pivot_two"),
arg("pivot_three")))
.def("CalculateBackboneScores", &WrapCalculateBackboneScores,
- (arg("score_container"), arg("scorer"),
+ (arg("score_container"), arg("mhandle"),
arg("start_resnum"), arg("chain_idx")=0))
.def("CalculateBackboneScores", &WrapCalculateBackboneScoresK,
- (arg("score_container"), arg("scorer"), arg("keys"),
+ (arg("score_container"), arg("mhandle"), arg("keys"),
arg("start_resnum"), arg("chain_idx")=0))
.def("CalculateAllAtomScores", &WrapCalculateAllAtomScores,
(arg("score_container"), arg("mhandle"),
diff --git a/modelling/pymod/export_monte_carlo.cc b/modelling/pymod/export_monte_carlo.cc
index 84d614fd599d31afc8f58e7cfda948e06168aae6..7d296ccf9550ec7db9ebf53c06aa4e4792370e0f 100644
--- a/modelling/pymod/export_monte_carlo.cc
+++ b/modelling/pymod/export_monte_carlo.cc
@@ -225,12 +225,14 @@ CCDCloserPtr CCDCloserWrapperList(const ost::mol::ResidueHandle& n_stem,
}
LinearScorerPtr LinearScorerInitWrapper(scoring::BackboneOverallScorerPtr scorer,
- uint start_resnum, uint chain_index,
+ scoring::BackboneScoreEnvPtr scorer_env,
+ uint start_resnum, uint num_residues,
+ uint chain_index,
const boost::python::dict& weights) {
std::map<String,Real> m_weights;
core::ConvertDictToMap(weights, m_weights);
- LinearScorerPtr p(new LinearScorer(scorer, start_resnum, chain_index,
- m_weights));
+ LinearScorerPtr p(new LinearScorer(scorer, scorer_env, start_resnum,
+ num_residues, chain_index, m_weights));
return p;
}
diff --git a/modelling/src/loop_candidate.cc b/modelling/src/loop_candidate.cc
index 82fda8923c1cab9c909f45903756889f6de2f577..e622837b2d5aa435f5b1c5b5aa1b15dc898318bb 100644
--- a/modelling/src/loop_candidate.cc
+++ b/modelling/src/loop_candidate.cc
@@ -414,23 +414,38 @@ LoopCandidatesPtr LoopCandidates::GetLargestCluster(Real max_dist,
void LoopCandidates::CalculateBackboneScores(
ScoreContainer& score_container,
- scoring::BackboneOverallScorerPtr scorer,
+ const ModellingHandle& mhandle,
uint start_resnum, uint chain_idx) const {
+
const std::vector<String>&
keys = ScoringWeights::GetInstance().GetBackboneScoringKeys();
- CalculateBackboneScores(score_container, scorer, keys, start_resnum,
+ CalculateBackboneScores(score_container, mhandle, keys, start_resnum,
chain_idx);
}
void LoopCandidates::CalculateBackboneScores(
ScoreContainer& score_container,
- scoring::BackboneOverallScorerPtr scorer,
+ const ModellingHandle& mhandle,
const std::vector<String>& keys,
uint start_resnum, uint chain_idx) const {
+
+
+ scoring::BackboneOverallScorerPtr scorer = mhandle.backbone_scorer;
+ scoring::BackboneScoreEnvPtr env = mhandle.backbone_scorer_env;
+
+ uint N = candidates_.size();
+ uint L = sequence_.size();
+
+ env->Stash(start_resnum, L, chain_idx);
+
for (uint i = 0; i < keys.size(); ++i) {
- std::vector<Real> scores;
- scores = scorer->Calculate(keys[i], candidates_, start_resnum, chain_idx);
+ std::vector<Real> scores(N);
+ for(uint j = 0; j < N; ++j) {
+ env->SetEnvironment(candidates_[j], start_resnum, chain_idx);
+ scores[j] = scorer->Calculate(keys[i], start_resnum, L, chain_idx);
+ }
score_container.Set(keys[i], scores);
}
+ env->Pop();
}
void LoopCandidates::CalculateAllAtomScores(
diff --git a/modelling/src/loop_candidate.hh b/modelling/src/loop_candidate.hh
index 39b47b8a53663ab18dd47b177e813849ca075737..06cfea97b9ab7afed41c9b6471efcdeefac77486 100644
--- a/modelling/src/loop_candidate.hh
+++ b/modelling/src/loop_candidate.hh
@@ -88,10 +88,10 @@ public:
bool superposed_rmsd = false) const;
void CalculateBackboneScores(ScoreContainer& score_container,
- scoring::BackboneOverallScorerPtr scorer,
+ const ModellingHandle& mhandle,
uint start_resnum, uint chain_idx = 0) const;
void CalculateBackboneScores(ScoreContainer& score_container,
- scoring::BackboneOverallScorerPtr scorer,
+ const ModellingHandle& mhandle,
const std::vector<String>& keys,
uint start_resnum, uint chain_idx = 0) const;
void CalculateAllAtomScores(ScoreContainer& score_container,
diff --git a/modelling/src/monte_carlo_scorer.cc b/modelling/src/monte_carlo_scorer.cc
index f8dd0ac6cfe01f7e72a1aedbcdbb1c3f84132dd4..ba536e7d772fc8b8d506c8375269e4d9fb528da1 100644
--- a/modelling/src/monte_carlo_scorer.cc
+++ b/modelling/src/monte_carlo_scorer.cc
@@ -5,21 +5,34 @@ namespace promod3 { namespace modelling {
MonteCarloScorer::~MonteCarloScorer() { }
LinearScorer::LinearScorer(scoring::BackboneOverallScorerPtr scorer,
- uint start_resnum, uint chain_index,
- const std::map<String,Real>& weights)
- : scorer_(scorer)
- , chain_index_(chain_index)
- , start_resnum_(start_resnum)
+ scoring::BackboneScoreEnvPtr scorer_env,
+ uint start_resnum, uint num_residues,
+ uint chain_index,
+ const std::map<String,Real>& weights):
+ scorer_(scorer),
+ scorer_env_(scorer_env),
+ start_resnum_(start_resnum),
+ num_residues_(num_residues),
+ chain_index_(chain_index)
{
w_scorers_ = scorer->GetWeightedScorers(weights);
}
Real LinearScorer::GetScore(const loop::BackboneList& positions) {
- Real score = 0;
+
+ if(positions.size() != num_residues_) {
+ throw promod3::Error("Input BackboneList is inconsistent with loop length "
+ "that you defined at initialization of the "
+ "LinearScorer");
+ }
+
+ scorer_env_->SetEnvironment(positions, start_resnum_, chain_index_);
+ Real score = 0.0;
for (uint i = 0; i < w_scorers_.size(); ++i) {
const Real weight = w_scorers_[i].first;
scoring::BackboneScorerPtr one_scorer = w_scorers_[i].second;
- score += weight * one_scorer->CalculateScore(positions, start_resnum_,
+ score += weight * one_scorer->CalculateScore(start_resnum_,
+ num_residues_,
chain_index_);
}
return score;
diff --git a/modelling/src/monte_carlo_scorer.hh b/modelling/src/monte_carlo_scorer.hh
index 6bd3e47b34675df0578affed4ccc289ea294a47b..a6707b359e0a138f296120815e524bfd9f4488b8 100644
--- a/modelling/src/monte_carlo_scorer.hh
+++ b/modelling/src/monte_carlo_scorer.hh
@@ -25,15 +25,24 @@ public:
class LinearScorer : public MonteCarloScorer {
public:
- LinearScorer(scoring::BackboneOverallScorerPtr scorer, uint start_resnum,
+ LinearScorer(scoring::BackboneOverallScorerPtr scorer,
+ scoring::BackboneScoreEnvPtr scorer_env,
+ uint start_resnum, uint num_residues,
uint chain_index, const std::map<String,Real>& weights);
Real GetScore(const loop::BackboneList& positions);
+ void StashScoringEnv() { scorer_env_->Stash(start_resnum_, num_residues_,
+ chain_index_); }
+
+ void PopScoringEnv() { scorer_env_->Pop(); }
+
private:
scoring::BackboneOverallScorerPtr scorer_;
- uint chain_index_;
+ scoring::BackboneScoreEnvPtr scorer_env_;
uint start_resnum_;
+ uint num_residues_;
+ uint chain_index_;
std::vector<std::pair<Real, scoring::BackboneScorerPtr> > w_scorers_;
};
diff --git a/scoring/pymod/CMakeLists.txt b/scoring/pymod/CMakeLists.txt
index 278cc7c320fa13c02baeeb1021241e60e9a50d36..204b23c99faca4e12550be667881f6bf5e3cb867 100644
--- a/scoring/pymod/CMakeLists.txt
+++ b/scoring/pymod/CMakeLists.txt
@@ -5,6 +5,7 @@ export_object_loader.cc
export_pairwise_functions.cc
export_constraint_constructor.cc
export_scwrl3_energy_functions.cc
+export_density_scorer.cc
wrap_scoring.cc
)
diff --git a/scoring/pymod/export_backbone_score.cc b/scoring/pymod/export_backbone_score.cc
index 9bb90032a28d07d791913eb156448baaab0cd3bf..ae4415881b08dc8a6009ef09837c6f16eefd96fe 100644
--- a/scoring/pymod/export_backbone_score.cc
+++ b/scoring/pymod/export_backbone_score.cc
@@ -8,7 +8,6 @@
#include <promod3/scoring/clash_score.hh>
#include <promod3/scoring/hbond_score.hh>
#include <promod3/scoring/ss_agreement_score.hh>
-#include <promod3/scoring/density_score.hh>
#include <promod3/scoring/torsion_score.hh>
#include <promod3/scoring/pairwise_score.hh>
@@ -68,6 +67,27 @@ namespace {
env->SetPsipredPrediction(pp_v);
}
+ void WrapStashSingle(BackboneScoreEnvPtr env, uint rnum, uint num_residues,
+ uint chain_idx) {
+ env->Stash(rnum, num_residues, chain_idx);
+ }
+
+
+ void WrapStashList(BackboneScoreEnvPtr env, const boost::python::list& rnum,
+ const boost::python::list& num_residues,
+ const boost::python::list& chain_idx) {
+
+ std::vector<uint> v_rnum;
+ std::vector<uint> v_num_residues;
+ std::vector<uint> v_chain_idx;
+
+ core::ConvertListToVector(rnum, v_rnum);
+ core::ConvertListToVector(num_residues, v_num_residues);
+ core::ConvertListToVector(chain_idx, v_chain_idx);
+
+ env->Stash(v_rnum, v_num_residues, v_chain_idx);
+ }
+
BackboneScorerPtr bos_getitem(BackboneOverallScorer& s, const String& key) {
return s[key];
}
@@ -76,52 +96,18 @@ namespace {
s[key] = item;
}
Real WrapCalculate(BackboneOverallScorer& s, const String& key,
- const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) {
- return s.Calculate(key, bb_list, start_resnum, chain_idx);
- }
- boost::python::list WrapCalculateList(
- BackboneOverallScorer& s, const String& key,
- const boost::python::list& py_list,
- uint start_resnum, uint chain_idx) {
- // get bb_lists
- std::vector<loop::BackboneList> bb_lists;
- core::ConvertListToVector(py_list, bb_lists);
- // get scores
- std::vector<Real> scores = s.Calculate(key, bb_lists, start_resnum,
- chain_idx);
- // convert and return
- boost::python::list return_list;
- core::AppendVectorToList(scores, return_list);
- return return_list;
+ uint start_resnum, uint num_residues, uint chain_idx) {
+ return s.Calculate(key, start_resnum, num_residues, chain_idx);
}
Real WrapCalculateLinearCombination(BackboneOverallScorer& s,
const boost::python::dict& weights,
- const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) {
+ uint start_resnum, uint num_residues,
+ uint chain_idx) {
std::map<String, Real> w;
core::ConvertDictToMap(weights, w);
- return s.CalculateLinearCombination(w, bb_list, start_resnum, chain_idx);
- }
- boost::python::list WrapCalculateLinearCombinationList(
- BackboneOverallScorer& s,
- const boost::python::dict& weights,
- const boost::python::list& py_list,
- uint start_resnum, uint chain_idx) {
- // get bb_lists
- std::vector<loop::BackboneList> bb_lists;
- core::ConvertListToVector(py_list, bb_lists);
- // get scores
- std::map<String, Real> w;
- core::ConvertDictToMap(weights, w);
- std::vector<Real> scores =
- s.CalculateLinearCombination(w, bb_lists, start_resnum, chain_idx);
- // convert and return
- boost::python::list return_list;
- core::AppendVectorToList(scores, return_list);
- return return_list;
+ return s.CalculateLinearCombination(w, start_resnum, num_residues,
+ chain_idx);
}
-
TorsionScorerPtr WrapTorsionScorerInit(
const boost::python::list& group_identifier,
uint bins_per_dimension) {
@@ -132,10 +118,10 @@ namespace {
}
boost::python::list WrapScoreProfile(BackboneScorerPtr scorer,
- const promod3::loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx){
+ uint start_resnum, uint num_residues,
+ uint chain_idx){
std::vector<Real> profile;
- scorer->CalculateScoreProfile(bb_list, start_resnum, chain_idx, profile);
+ scorer->CalculateScoreProfile(start_resnum, num_residues, chain_idx, profile);
boost::python::list return_list;
core::AppendVectorToList(profile, return_list);
return return_list;
@@ -166,13 +152,16 @@ void export_BackboneScore() {
(arg("start_resnum"), arg("num_residues"), arg("chain_idx") = 0))
.def("SetPsipredPrediction", &WrapSetPsipredPrediction, (arg("pred")))
.def("SetPsipredPrediction", &WrapSetPsipredPredictionList, (arg("pred")))
- .def("SetMap", &BackboneScoreEnv::SetMap,
- (arg("map"), arg("resolution"), arg("all_atom")=false))
.def("AddPairwiseFunction", &BackboneScoreEnv::AddPairwiseFunction,
(arg("function"), arg("function_type")))
.def("ApplyPairwiseFunction", &BackboneScoreEnv::ApplyPairwiseFunction,
(arg("chain_idx_one"), arg("resnum_one"), arg("chain_idx_two"),
arg("resnum_two"), arg("f_idx")))
+ .def("Stash", &WrapStashSingle, (arg("start_rnum"), arg("num_residues"),
+ arg("chain_idx")=0))
+ .def("Stash", &WrapStashList, (arg("start_rnums"), arg("num_residues"),
+ arg("chain_indices")=0))
+ .def("Pop", &BackboneScoreEnv::Pop)
.def("GetSeqres", &BackboneScoreEnv::GetSeqres)
;
@@ -190,14 +179,9 @@ void export_BackboneScore() {
.def("AttachEnvironment", &BackboneOverallScorer::AttachEnvironment,
(arg("env")))
.def("Calculate", &WrapCalculate,
- (arg("key"), arg("bb_list"), arg("start_resnum"), arg("chain_idx")=0))
- .def("Calculate", &WrapCalculateList,
- (arg("key"), arg("bb_lists"), arg("start_resnum"), arg("chain_idx")=0))
+ (arg("key"), arg("start_resnum"), arg("num_residues"), arg("chain_idx")=0))
.def("CalculateLinearCombination", &WrapCalculateLinearCombination,
- (arg("linear_weights"), arg("bb_list"), arg("start_resnum"),
- arg("chain_idx")=0))
- .def("CalculateLinearCombination", &WrapCalculateLinearCombinationList,
- (arg("linear_weights"), arg("bb_lists"), arg("start_resnum"),
+ (arg("linear_weights"), arg("start_resnum"), arg("num_residues"),
arg("chain_idx")=0))
;
@@ -209,10 +193,8 @@ void export_BackboneScore() {
.def("LoadPortable", &CBPackingScorer::LoadPortable,(arg("filename"))).staticmethod("LoadPortable")
.def("SavePortable", &CBPackingScorer::SavePortable,(arg("filename")))
.def("SetEnergy", &CBPackingScorer::SetEnergy,(arg("aa"),arg("count"),arg("energy")))
- .def("CalculateScore", &CBPackingScorer::CalculateScore,(arg("bb_list"), arg("start_resnum"), arg("chain_idx") = 0))
- .def("CalculateScoreProfile", &WrapScoreProfile,(arg("bb_list"), arg("start_resnum"), arg("chain_idx") = 0))
- .def("UseClassicMode", &CBPackingScorer::UseClassicMode)
- .def("UseIncludeEnvMode",&CBPackingScorer::UseIncludeEnvMode)
+ .def("CalculateScore", &CBPackingScorer::CalculateScore,(arg("start_resnum"), arg("num_residues"), arg("chain_idx") = 0))
+ .def("CalculateScoreProfile", &WrapScoreProfile,(arg("start_resnum"), arg("num_residues"), arg("chain_idx") = 0))
.def("DoNormalize", &CBPackingScorer::DoNormalize)
.def("AttachEnvironment", &CBPackingScorer::AttachEnvironment,(arg("env")))
;
@@ -224,8 +206,8 @@ void export_BackboneScore() {
.def("LoadPortable", &CBetaScorer::LoadPortable,(arg("filename"))).staticmethod("LoadPortable")
.def("SavePortable", &CBetaScorer::SavePortable,(arg("filename")))
.def("SetEnergy", &CBetaScorer::SetEnergy,(arg("aa1"),arg("aa2"),arg("bin"),arg("energy")))
- .def("CalculateScore", &CBetaScorer::CalculateScore,(arg("bb_list"), arg("start_resnum"), arg("chain_idx") = 0))
- .def("CalculateScoreProfile", &WrapScoreProfile,(arg("bb_list"), arg("start_resnum"), arg("chain_idx") = 0))
+ .def("CalculateScore", &CBetaScorer::CalculateScore,(arg("start_resnum"), arg("num_residues"), arg("chain_idx") = 0))
+ .def("CalculateScoreProfile", &WrapScoreProfile,(arg("start_resnum"), arg("num_residues"), arg("chain_idx") = 0))
.def("DoInternalScores", &CBetaScorer::DoInternalScores)
.def("DoExternalScores", &CBetaScorer::DoExternalScores)
.def("DoNormalize", &CBetaScorer::DoNormalize)
@@ -239,8 +221,8 @@ void export_BackboneScore() {
.def("LoadPortable", &ReducedScorer::LoadPortable,(arg("filename"))).staticmethod("LoadPortable")
.def("SavePortable", &ReducedScorer::SavePortable,(arg("filename")))
.def("SetEnergy", &ReducedScorer::SetEnergy,(arg("aa1"),arg("aa2"),arg("dist_bin"),arg("alpha_bin"),arg("beta_bin"),arg("gamma_bin"),arg("energy")))
- .def("CalculateScore", &ReducedScorer::CalculateScore,(arg("bb_list"), arg("start_resnum"), arg("chain_idx") = 0))
- .def("CalculateScoreProfile", &WrapScoreProfile,(arg("bb_list"), arg("start_resnum"), arg("chain_idx") = 0))
+ .def("CalculateScore", &ReducedScorer::CalculateScore,(arg("start_resnum"), arg("num_residues"), arg("chain_idx") = 0))
+ .def("CalculateScoreProfile", &WrapScoreProfile,(arg("start_resnum"), arg("num_residues"), arg("chain_idx") = 0))
.def("DoInternalScores", &ReducedScorer::DoInternalScores)
.def("DoExternalScores", &ReducedScorer::DoExternalScores)
.def("DoNormalize", &ReducedScorer::DoNormalize)
@@ -249,21 +231,14 @@ void export_BackboneScore() {
class_<ClashScorer, ClashScorerPtr, bases<BackboneScorer> >
("ClashScorer",init<>())
- .def("CalculateScore", &ClashScorer::CalculateScore,(arg("bb_list"), arg("start_resnum"), arg("chain_idx") = 0))
- .def("CalculateScoreProfile", &WrapScoreProfile,(arg("bb_list"), arg("start_resnum"), arg("chain_idx") = 0))
+ .def("CalculateScore", &ClashScorer::CalculateScore,(arg("start_resnum"), arg("num_residues"), arg("chain_idx") = 0))
+ .def("CalculateScoreProfile", &WrapScoreProfile,(arg("start_resnum"), arg("num_residues"), arg("chain_idx") = 0))
.def("AttachEnvironment", &ClashScorer::AttachEnvironment,(arg("env")))
.def("DoInternalScores", &ClashScorer::DoInternalScores)
.def("DoExternalScores", &ClashScorer::DoExternalScores)
.def("DoNormalize", &ClashScorer::DoNormalize)
;
- class_<DensityScorer, DensityScorerPtr, bases<BackboneScorer> >
- ("DensityScorer", init<>())
- .def("CalculateScore", &DensityScorer::CalculateScore,(arg("bb_list"), arg("start_resnum"), arg("chain_idx") = 0))
- .def("CalculateScoreProfile", &WrapScoreProfile,(arg("bb_list"), arg("start_resnum"), arg("chain_idx") = 0))
- .def("AttachEnvironment", &DensityScorer::AttachEnvironment,(arg("env")))
- ;
-
class_<HBondScorer, HBondScorerPtr, bases<BackboneScorer> >
("HBondScorer", init<Real,Real,Real,Real,Real,Real,Real,Real,uint,uint,uint,uint>())
.def("Load", &HBondScorer::Load,(arg("filename"))).staticmethod("Load")
@@ -271,8 +246,8 @@ void export_BackboneScore() {
.def("LoadPortable", &HBondScorer::LoadPortable,(arg("filename"))).staticmethod("LoadPortable")
.def("SavePortable", &HBondScorer::SavePortable,(arg("filename")))
.def("SetEnergy", &HBondScorer::SetEnergy,(arg("state"),arg("d_bin"),arg("alpha_bin"),arg("beta_bin"),arg("gamma_bin"),arg("energy")))
- .def("CalculateScore", &HBondScorer::CalculateScore,(arg("bb_list"), arg("start_resnum"), arg("chain_idx") = 0))
- .def("CalculateScoreProfile", &WrapScoreProfile,(arg("bb_list"), arg("start_resnum"), arg("chain_idx") = 0))
+ .def("CalculateScore", &HBondScorer::CalculateScore,(arg("start_resnum"), arg("num_residues"), arg("chain_idx") = 0))
+ .def("CalculateScoreProfile", &WrapScoreProfile,(arg("start_resnum"), arg("num_residues"), arg("chain_idx") = 0))
.def("DoInternalScores", &HBondScorer::DoInternalScores)
.def("DoExternalScores", &HBondScorer::DoExternalScores)
.def("DoNormalize", &HBondScorer::DoNormalize)
@@ -286,8 +261,8 @@ void export_BackboneScore() {
.def("LoadPortable", &SSAgreementScorer::LoadPortable,(arg("filename"))).staticmethod("LoadPortable")
.def("SavePortable", &SSAgreementScorer::SavePortable,(arg("filename")))
.def("SetScore", &SSAgreementScorer::SetScore,(arg("psipred_state"),arg("psipred_cfi"),arg("dssp_state"),arg("score")))
- .def("CalculateScore", &SSAgreementScorer::CalculateScore,(arg("bb_list"), arg("start_resnum"), arg("chain_idx") = 0))
- .def("CalculateScoreProfile", &WrapScoreProfile,(arg("bb_list"), arg("start_resnum"), arg("chain_idx") = 0))
+ .def("CalculateScore", &SSAgreementScorer::CalculateScore,(arg("start_resnum"), arg("num_residues"), arg("chain_idx") = 0))
+ .def("CalculateScoreProfile", &WrapScoreProfile,(arg("start_resnum"), arg("num_residues"), arg("chain_idx") = 0))
.def("DoNormalize", &SSAgreementScorer::DoNormalize)
.def("AttachEnvironment", &SSAgreementScorer::AttachEnvironment,(arg("env")))
;
@@ -300,16 +275,16 @@ void export_BackboneScore() {
.def("LoadPortable", &TorsionScorer::LoadPortable,(arg("filename"))).staticmethod("LoadPortable")
.def("SavePortable", &TorsionScorer::SavePortable,(arg("filename")))
.def("SetEnergy", &TorsionScorer::SetEnergy,(arg("group_idx"),arg("phi_bin"),arg("psi_bin"),arg("energy")))
- .def("CalculateScore", &TorsionScorer::CalculateScore,(arg("bb_list"), arg("start_resnum"), arg("chain_idx") = 0))
- .def("CalculateScoreProfile", &WrapScoreProfile,(arg("bb_list"), arg("start_resnum"), arg("chain_idx") = 0))
+ .def("CalculateScore", &TorsionScorer::CalculateScore,(arg("start_resnum"), arg("num_residues"), arg("chain_idx") = 0))
+ .def("CalculateScoreProfile", &WrapScoreProfile,(arg("start_resnum"), arg("num_residues"), arg("chain_idx") = 0))
.def("DoNormalize", &TorsionScorer::DoNormalize)
.def("AttachEnvironment", &TorsionScorer::AttachEnvironment,(arg("env")))
;
class_<PairwiseScorer, PairwiseScorerPtr, bases<BackboneScorer> >
("PairwiseScorer", init<>())
- .def("CalculateScore", &PairwiseScorer::CalculateScore,(arg("bb_list"), arg("start_resnum"), arg("chain_idx") = 0))
- .def("CalculateScoreProfile", &WrapScoreProfile,(arg("bb_list"), arg("start_resnum"), arg("chain_idx") = 0))
+ .def("CalculateScore", &PairwiseScorer::CalculateScore,(arg("start_resnum"), arg("num_residues"), arg("chain_idx") = 0))
+ .def("CalculateScoreProfile", &WrapScoreProfile,(arg("start_resnum"), arg("num_residues"), arg("chain_idx") = 0))
.def("AttachEnvironment", &PairwiseScorer::AttachEnvironment,(arg("env")))
.def("DoInternalScores", &PairwiseScorer::DoInternalScores)
.def("DoExternalScores", &PairwiseScorer::DoExternalScores)
diff --git a/scoring/pymod/export_density_scorer.cc b/scoring/pymod/export_density_scorer.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f022f963aee63b41120ca950a1db90d89b14ffec
--- /dev/null
+++ b/scoring/pymod/export_density_scorer.cc
@@ -0,0 +1,11 @@
+#include <boost/python.hpp>
+#include <promod3/scoring/density_score.hh>
+
+using namespace boost::python;
+
+void export_DensityScoring() {
+ def("DensityScore", &promod3::scoring::DensityScore, (arg("bb_list"),
+ arg("map"),
+ arg("resolution"),
+ arg("all_atom")));
+}
diff --git a/scoring/pymod/wrap_scoring.cc b/scoring/pymod/wrap_scoring.cc
index 779239c08317d9973b2c44f82145a1e1e9ebedd3..4b717f9cbc40743c1e84fcef141dc8e95bdb874f 100644
--- a/scoring/pymod/wrap_scoring.cc
+++ b/scoring/pymod/wrap_scoring.cc
@@ -6,6 +6,7 @@ void export_PairwiseFunction();
void export_constraint_constructor();
void export_AllAtomScore();
void export_SCWRL3EnergyFunctions();
+void export_DensityScoring();
using namespace boost::python;
@@ -17,4 +18,5 @@ BOOST_PYTHON_MODULE(_scoring)
export_PairwiseFunction();
export_constraint_constructor();
export_SCWRL3EnergyFunctions();
+ export_DensityScoring();
}
diff --git a/scoring/src/backbone_overall_scorer.cc b/scoring/src/backbone_overall_scorer.cc
index 1933a0d73f0d12bc9865d7e9d5c80d96bee71f67..0991f64c0e602088a9c1394cbee6e0fff45c1c29 100644
--- a/scoring/src/backbone_overall_scorer.cc
+++ b/scoring/src/backbone_overall_scorer.cc
@@ -30,30 +30,15 @@ BackboneOverallScorer::GetWeightedScorers(
}
Real BackboneOverallScorer::Calculate(const String& key,
- const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const {
- return this->Get(key)->CalculateScore(bb_list, start_resnum, chain_idx);
-}
-
-std::vector<Real> BackboneOverallScorer::Calculate(
- const String& key,
- const std::vector<loop::BackboneList>& bb_lists,
- uint start_resnum, uint chain_idx) const {
- // get scorer once (avoid check / lookup)
- BackboneScorerPtr scorer = this->Get(key);
- // score all
- const uint N = bb_lists.size();
- std::vector<Real> scores(N);
- for (uint i = 0; i < N; ++i) {
- scores[i] = scorer->CalculateScore(bb_lists[i], start_resnum, chain_idx);
- }
- return scores;
+ uint start_resnum, uint num_residues,
+ uint chain_idx) const {
+ return this->Get(key)->CalculateScore(start_resnum, num_residues, chain_idx);
}
Real BackboneOverallScorer::CalculateLinearCombination(
const std::map<String, Real>& linear_weights,
- const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const {
+ uint start_resnum, uint num_residues,
+ uint chain_idx) const {
// loop and sum up
Real score = 0;
for (std::map<String, Real>::const_iterator i = linear_weights.begin();
@@ -62,35 +47,11 @@ Real BackboneOverallScorer::CalculateLinearCombination(
score += i->second;
} else {
// note: error checking done by Calculate
- score += i->second * this->Calculate(i->first, bb_list, start_resnum,
+ score += i->second * this->Calculate(i->first, start_resnum, num_residues,
chain_idx);
}
}
return score;
}
-std::vector<Real> BackboneOverallScorer::CalculateLinearCombination(
- const std::map<String, Real>& linear_weights,
- const std::vector<loop::BackboneList>& bb_lists,
- uint start_resnum, uint chain_idx) const {
- // get scorers once (avoid check / lookup)
- std::vector<std::pair<Real, BackboneScorerPtr> > scorers;
- scorers = this->GetWeightedScorers(linear_weights);
- const uint Nscores = scorers.size();
- // score all
- const uint N = bb_lists.size();
- std::vector<Real> scores(N);
- for (uint i = 0; i < N; ++i) {
- Real score = 0;
- for (uint j = 0; j < Nscores; ++j) {
- const Real weight = scorers[j].first;
- BackboneScorerPtr scorer = scorers[j].second;
- score += weight * scorer->CalculateScore(bb_lists[i], start_resnum,
- chain_idx);
- }
- scores[i] = score;
- }
- return scores;
-}
-
}} // ns
\ No newline at end of file
diff --git a/scoring/src/backbone_overall_scorer.hh b/scoring/src/backbone_overall_scorer.hh
index 17016a75de74d4eaef1703e3598c0e5f3011a90f..c9d9aac48961b686e099979ac50719a7846b0c87 100644
--- a/scoring/src/backbone_overall_scorer.hh
+++ b/scoring/src/backbone_overall_scorer.hh
@@ -13,15 +13,14 @@ typedef boost::shared_ptr<BackboneOverallScorer> BackboneOverallScorerPtr;
// dummy score for intercept -> always returns 1
class InterceptScorer : public BackboneScorer {
public:
- virtual Real CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum,
+ virtual Real CalculateScore(uint start_resnum, uint num_residues,
uint chain_idx) const { return 1.0; }
- virtual void CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum,
+ virtual void CalculateScoreProfile(uint start_resnum,
+ uint num_residues,
uint chain_idx,
std::vector<Real>& profile) const{
- profile.assign(bb_list.size(), 1.0);
+ profile.assign(num_residues, 1.0);
}
};
@@ -56,21 +55,14 @@ public:
GetWeightedScorers(const std::map<String, Real>& linear_weights) const;
// calculate one single score
- Real Calculate(const String& key, const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx = 0) const;
- std::vector<Real>
- Calculate(const String& key, const std::vector<loop::BackboneList>& bb_lists,
- uint start_resnum, uint chain_idx = 0) const;
+ Real Calculate(const String& key, uint start_resnum, uint num_residues,
+ uint chain_idx = 0) const;
// linear combination of weighted scores
// -> use special key "intercept" to add a constant value to the score
Real CalculateLinearCombination(const std::map<String, Real>& linear_weights,
- const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx = 0) const;
- std::vector<Real>
- CalculateLinearCombination(const std::map<String, Real>& linear_weights,
- const std::vector<loop::BackboneList>& bb_lists,
- uint start_resnum, uint chain_idx = 0) const;
+ uint start_resnum, uint num_residues,
+ uint chain_idx = 0) const;
private:
std::map<String, BackboneScorerPtr> map_;
diff --git a/scoring/src/backbone_score_base.cc b/scoring/src/backbone_score_base.cc
index c1bcaa59398baf2d51bcf6bb5c559d3f19618a91..b01b9dcf00049663edf1564519d9cd2b78679b4e 100644
--- a/scoring/src/backbone_score_base.cc
+++ b/scoring/src/backbone_score_base.cc
@@ -63,9 +63,6 @@ BackboneScoreEnvPtr BackboneScoreEnv::Copy() const{
return_ptr->psipred_pred_ = psipred_pred_;
return_ptr->psipred_cfi_ = psipred_cfi_;
return_ptr->env_set_ = env_set_;
- if(density_map_data_.map.IsValid()) {
- return_ptr->density_map_data_ = density_map_data_;
- }
return_ptr->n_pos_ = n_pos_;
return_ptr->ca_pos_ = ca_pos_;
return_ptr->cb_pos_ = cb_pos_;
@@ -324,13 +321,6 @@ void BackboneScoreEnv::SetPsipredPrediction(std::vector<loop::PsipredPredictionP
}
}
-void BackboneScoreEnv::SetMap(const ost::img::MapHandle& map, Real resolution,
- bool all_atom) {
- density_map_data_.map = map;
- density_map_data_.resolution = resolution;
- density_map_data_.all_atom = all_atom;
-}
-
uint BackboneScoreEnv::AddPairwiseFunction(PairwiseFunctionPtr f,
PairwiseFunctionType ft) {
pairwise_functions_.push_back(f);
@@ -371,6 +361,127 @@ void BackboneScoreEnv::ApplyPairwiseFunction(uint chain_idx_one, uint resnum_one
applied_pairwise_functions_[idx_two].push_back(data_two);
}
+void BackboneScoreEnv::Stash(uint start_resnum, uint num_residues,
+ uint chain_idx) {
+ std::vector<uint> indices = idx_handler_.ToIdxVector(start_resnum,
+ num_residues,
+ chain_idx);
+ this->Stash(indices);
+}
+
+void BackboneScoreEnv::Stash(const std::vector<uint>& start_resnum,
+ const std::vector<uint>& num_residues,
+ const std::vector<uint>& chain_idx) {
+
+ std::vector<uint> start_indices;
+ std::vector<uint> lengths;
+
+ std::vector<uint> indices = idx_handler_.ToIdxVector(start_resnum,
+ num_residues,
+ chain_idx,
+ start_indices,
+ lengths);
+ this->Stash(indices);
+}
+
+void BackboneScoreEnv::Stash(const std::vector<uint>& indices) {
+
+ if(stashed_indices_.size() > 100) {
+ throw promod3::Error("You already performed a 100 stash events... "
+ "Cant do more....");
+ }
+
+ std::vector<uint> stash_indices;
+ std::vector<int> stash_env_set;
+ std::vector<geom::Vec3> stash_n_pos;
+ std::vector<geom::Vec3> stash_ca_pos;
+ std::vector<geom::Vec3> stash_cb_pos;
+ std::vector<geom::Vec3> stash_c_pos;
+ std::vector<geom::Vec3> stash_o_pos;
+
+ uint n_indices = indices.size();
+ stash_indices.reserve(n_indices);
+ stash_env_set.reserve(n_indices);
+ stash_n_pos.reserve(n_indices);
+ stash_ca_pos.reserve(n_indices);
+ stash_cb_pos.reserve(n_indices);
+ stash_c_pos.reserve(n_indices);
+ stash_o_pos.reserve(n_indices);
+
+ for(std::vector<uint>::const_iterator i = indices.begin();
+ i != indices.end(); ++i) {
+ stash_indices.push_back(*i);
+ stash_env_set.push_back(env_set_[*i]);
+ stash_n_pos.push_back(n_pos_[*i]);
+ stash_ca_pos.push_back(ca_pos_[*i]);
+ stash_cb_pos.push_back(cb_pos_[*i]);
+ stash_c_pos.push_back(c_pos_[*i]);
+ stash_o_pos.push_back(o_pos_[*i]);
+ }
+
+ stashed_indices_.push(stash_indices);
+ stashed_env_set_.push(stash_env_set);
+ stashed_n_pos_.push(stash_n_pos);
+ stashed_ca_pos_.push(stash_ca_pos);
+ stashed_cb_pos_.push(stash_cb_pos);
+ stashed_c_pos_.push(stash_c_pos);
+ stashed_o_pos_.push(stash_o_pos);
+}
+
+void BackboneScoreEnv::Pop() {
+
+ if(stashed_indices_.empty()) {
+ throw promod3::Error("Nothing stashed! Cannot pop!");
+ }
+
+ std::vector<uint> indices_to_set;
+ std::vector<uint> indices_to_reset;
+ std::vector<uint> indices_to_clear;
+
+ const std::vector<uint>& stash_indices = stashed_indices_.top();
+ const std::vector<int>& stash_env_set = stashed_env_set_.top();
+ const std::vector<geom::Vec3>& stash_n_pos = stashed_n_pos_.top();
+ const std::vector<geom::Vec3>& stash_ca_pos = stashed_ca_pos_.top();
+ const std::vector<geom::Vec3>& stash_cb_pos = stashed_cb_pos_.top();
+ const std::vector<geom::Vec3>& stash_c_pos = stashed_c_pos_.top();
+ const std::vector<geom::Vec3>& stash_o_pos = stashed_o_pos_.top();
+
+ for(uint i = 0; i < stash_indices.size(); ++i) {
+ uint env_idx = stash_indices[i];
+ n_pos_[env_idx] = stash_n_pos[i];
+ ca_pos_[env_idx] = stash_ca_pos[i];
+ cb_pos_[env_idx] = stash_cb_pos[i];
+ c_pos_[env_idx] = stash_c_pos[i];
+ o_pos_[env_idx] = stash_o_pos[i];
+
+ if(env_set_[env_idx] == 1 && stash_env_set[i] == 1) {
+ indices_to_reset.push_back(env_idx);
+ } else if(env_set_[env_idx] == 1 && stash_env_set[i] == 0) {
+ indices_to_clear.push_back(env_idx);
+ } else if(env_set_[env_idx] == 0 && stash_env_set[i] == 1) {
+ indices_to_set.push_back(env_idx);
+ }
+
+ env_set_[env_idx] = stash_env_set[i];
+ }
+
+ //and call all listeners
+ for(std::vector<BackboneScoreEnvListener*>::iterator i = listener_.begin();
+ i != listener_.end(); ++i){
+ (*i)->SetEnvironment(*this, indices_to_set);
+ (*i)->ResetEnvironment(*this, indices_to_reset);
+ (*i)->ClearEnvironment(*this, indices_to_clear);
+ }
+
+ stashed_indices_.pop();
+ stashed_env_set_.pop();
+ stashed_n_pos_.pop();
+ stashed_ca_pos_.pop();
+ stashed_cb_pos_.pop();
+ stashed_c_pos_.pop();
+ stashed_o_pos_.pop();
+}
+
BackboneScoreEnvListener* BackboneScoreEnv::GetListener(const String& id) const {
for (uint i = 0; i < listener_.size(); ++i) {
if (listener_[i]->WhoAmI() == id) return listener_[i];
diff --git a/scoring/src/backbone_score_base.hh b/scoring/src/backbone_score_base.hh
index feade33ead9a28e021a79d4f72f852f4bdce878b..0676f6d2da8928c67fe0ab5585b3f0872e183010 100644
--- a/scoring/src/backbone_score_base.hh
+++ b/scoring/src/backbone_score_base.hh
@@ -11,6 +11,7 @@
#include <promod3/scoring/pairwise_scoring_function.hh>
#include <algorithm>
+#include <stack>
namespace promod3 { namespace scoring {
@@ -37,13 +38,6 @@ struct PairwiseFunctionData {
};
-struct DensityMapData {
- ost::img::MapHandle map;
- Real resolution;
- bool all_atom;
-};
-
-
class BackboneScoreEnv {
public:
@@ -68,9 +62,6 @@ public:
void SetPsipredPrediction(std::vector<loop::PsipredPredictionPtr>& pp);
- void SetMap(const ost::img::MapHandle& map, Real resolution,
- bool all_atom = false);
-
uint AddPairwiseFunction(PairwiseFunctionPtr f,
PairwiseFunctionType ft);
@@ -78,6 +69,14 @@ public:
uint chain_idx_two, uint resnum_two,
uint f_idx);
+ void Stash(uint start_resnum, uint num_residues, uint chain_idx = 0);
+
+ void Stash(const std::vector<uint>& start_resnum,
+ const std::vector<uint>& num_residues,
+ const std::vector<uint>& chain_idx);
+
+ void Pop();
+
ost::seq::SequenceList GetSeqres() const { return seqres_; }
// methods to access internal data intended to be used by the scorer objects
@@ -105,9 +104,6 @@ public:
const std::vector<PairwiseFunctionData>*
GetPairwiseFunctionData() const { return &applied_pairwise_functions_[0]; }
- const DensityMapData* GetDensityMapData() const { return &density_map_data_; }
-
-
// listener related methods
// return listener with that ID or NULL if not existing
@@ -133,6 +129,8 @@ public:
private:
+ void Stash(const std::vector<uint>& indices);
+
ost::seq::SequenceList seqres_;
loop::IdxHandler idx_handler_;
std::vector<ost::conop::AminoAcid> aa_seqres_;
@@ -141,12 +139,19 @@ private:
// the reason to use a vector with int here is that the STL has an optimized
// vector<bool> (bitset, min. memory usage) which can't be used as C arrays
std::vector<int> env_set_;
- DensityMapData density_map_data_;
std::vector<geom::Vec3> n_pos_;
std::vector<geom::Vec3> ca_pos_;
std::vector<geom::Vec3> cb_pos_;
std::vector<geom::Vec3> c_pos_;
std::vector<geom::Vec3> o_pos_;
+ // same as above, but to stash things
+ std::stack<std::vector<uint> > stashed_indices_;
+ std::stack<std::vector<int> > stashed_env_set_;
+ std::stack<std::vector<geom::Vec3> > stashed_n_pos_;
+ std::stack<std::vector<geom::Vec3> > stashed_ca_pos_;
+ std::stack<std::vector<geom::Vec3> > stashed_cb_pos_;
+ std::stack<std::vector<geom::Vec3> > stashed_c_pos_;
+ std::stack<std::vector<geom::Vec3> > stashed_o_pos_;
// every residue gets a vector of pairs
// every pair contains:
// - a function ptr as a first element, i.e. the pairwise
@@ -188,13 +193,11 @@ public:
virtual ~BackboneScorer() { }
- virtual Real CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum,
+ virtual Real CalculateScore(uint start_resnum, uint num_residues,
uint chain_idx) const = 0;
- virtual void CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum,
- uint chain_idx,
+ virtual void CalculateScoreProfile(uint start_resnum, uint num_residues,
+ uint chain_idx,
std::vector<Real>& profile) const = 0;
virtual void AttachEnvironment(BackboneScoreEnv& env) { }
diff --git a/scoring/src/cb_packing_score.cc b/scoring/src/cb_packing_score.cc
index 8f971b73b83a39c7460add0a0589a00ba5e14c40..794777d32801be15d9d91196047179573cc7cbd0 100644
--- a/scoring/src/cb_packing_score.cc
+++ b/scoring/src/cb_packing_score.cc
@@ -7,7 +7,6 @@ CBPackingScorer::CBPackingScorer(Real cutoff, uint max_count):
cutoff_(cutoff),
max_count_(max_count),
attached_environment_(false),
- include_env_mode_(false),
normalize_(true),
env_(NULL) {
@@ -182,53 +181,44 @@ void CBPackingScorer::SetEnergy(ost::conop::AminoAcid aa, uint count, Real e) {
energies_[this->EIdx(aa, count)] = e;
}
-Real CBPackingScorer::CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const {
+Real CBPackingScorer::CalculateScore(uint start_resnum, uint num_residues,
+ uint chain_idx) const {
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
"CBPackingScorer::CalculateScore", 2);
- if(include_env_mode_) {
- return this->CalculateScoreIncludeEnv(bb_list, start_resnum,
- chain_idx);
- }
- else {
- std::vector<Real> scores;
- this->Score(bb_list, start_resnum, chain_idx, scores);
- Real score = 0.0;
- for(uint i = 0; i < scores.size(); ++i) {
- score += scores[i];
- }
+ std::vector<Real> scores;
+ std::vector<uint> idx_vector =
+ idx_handler_->ToIdxVector(start_resnum, num_residues, chain_idx);
- if(normalize_ && !scores.empty()) {
- score /= scores.size();
- }
+ this->Score(idx_vector, scores);
+ Real score = 0.0;
+ for(uint i = 0; i < scores.size(); ++i) {
+ score += scores[i];
+ }
- return score;
+ if(normalize_ && !scores.empty()) {
+ score /= scores.size();
}
+
+ return score;
}
-void CBPackingScorer::CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
+void CBPackingScorer::CalculateScoreProfile(uint start_resnum, uint num_residues,
+ uint chain_idx,
std::vector<Real>& profile) const {
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
"CBPackingScorer::CalculateScoreProfile", 2);
- if(include_env_mode_) {
- throw promod3::Error("CBPackingScorer cannot generate score profile when "
- "include_env_mode is active! Call the UseClassicMode()"
- "function to deactivate and try again...");
- }
-
- this->Score(bb_list, start_resnum, chain_idx, profile);
+ std::vector<uint> idx_vector =
+ idx_handler_->ToIdxVector(start_resnum, num_residues, chain_idx);
+ this->Score(idx_vector, profile);
}
-void CBPackingScorer::Score(const loop::BackboneList& bb_list,
- uint start_resnum,
- uint chain_idx,
+void CBPackingScorer::Score(const std::vector<uint>& indices,
std::vector<Real>& scores) const {
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
@@ -239,157 +229,46 @@ void CBPackingScorer::Score(const loop::BackboneList& bb_list,
String err = "Can only calculate score when environment is attached!";
throw promod3::Error(err);
}
- uint bb_list_size, start_idx, end_idx;
- idx_handler_->SetupScoreCalculation(bb_list, start_resnum, chain_idx,
- bb_list_size, start_idx, end_idx);
-
- scores.assign(bb_list_size, 0.0);
- if(bb_list_size == 0) return;
-
- // get counts
- std::vector<uint> count_buffer(bb_list_size, 0);
-
- // first internal
- const Real squared_cutoff = cutoff_ * cutoff_;
- for (uint i = 0; i < bb_list_size; ++i) {
- for (uint j = i+1; j < bb_list_size; ++j) {
- if (geom::Length2(bb_list.GetCB(i) - bb_list.GetCB(j)) < squared_cutoff) {
- count_buffer[i] += 1;
- count_buffer[j] += 1;
- }
- }
- }
- // then external and sum up score as we go
- CBetaSpatialOrganizer::WithinResult within_result;
- std::pair<CBetaSpatialOrganizerItem*,Real>* a;
- for (uint i = 0; i < bb_list_size; ++i) {
- // add neighbors to internal count
- within_result = env_->FindWithin(bb_list.GetCB(i), cutoff_);
- a = within_result.first;
- uint counter = count_buffer[i];
- for (uint j = 0; j < within_result.second; ++j) {
- if (a[j].first->idx < start_idx || a[j].first->idx > end_idx) {
- ++counter;
- }
- }
- // score it
- counter = std::min(counter, max_count_);
- scores[i] = energies_[this->EIdx(bb_list.GetAA(i), counter)];
- }
-}
-Real CBPackingScorer::CalculateScoreIncludeEnv(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const {
+ uint n_indices = indices.size();
- promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
- "CBPackingScorer::CalculateScoreIncludeEnv", 2);
+ scores.assign(n_indices, 0.0);
- // setup
- if(!attached_environment_) {
- String err = "Can only calculate score when environment is attached!";
- throw promod3::Error(err);
- }
- uint bb_list_size, start_idx, end_idx;
- std::map<uint,uint> close_env_residues;
- std::map<uint,uint>::iterator close_env_it;
- uint env_idx;
-
- idx_handler_->SetupScoreCalculation(bb_list, start_resnum, chain_idx,
- bb_list_size, start_idx, end_idx);
-
- // get counts
- std::vector<uint> count_buffer(bb_list_size,0);
-
- // first internal
- Real squared_cutoff = cutoff_ * cutoff_;
- for (uint i = 0; i < bb_list_size; ++i) {
- for (uint j = i+1; j < bb_list_size; ++j) {
- if (geom::Length2(bb_list.GetCB(i) - bb_list.GetCB(j)) < squared_cutoff) {
- count_buffer[i] += 1;
- count_buffer[j] += 1;
- }
- }
+ if(n_indices == 0) {
+ return;
}
- // score the bb_list items and keep track of the surroundings
- Real score = 0.0;
+ // then external and sum up score as we go
CBetaSpatialOrganizer::WithinResult within_result;
std::pair<CBetaSpatialOrganizerItem*,Real>* a;
- for (uint i = 0; i < bb_list_size; ++i) {
- // add neighbors to internal count
- within_result = env_->FindWithin(bb_list.GetCB(i), cutoff_);
- a = within_result.first;
- uint counter = count_buffer[i];
- for (uint j = 0; j < within_result.second; ++j) {
- env_idx = a[j].first->idx;
- if (env_idx < start_idx || env_idx > end_idx) {
- ++counter;
- close_env_it = close_env_residues.find(env_idx);
- if(close_env_it == close_env_residues.end()) {
- close_env_residues[env_idx] = 1;
- }
- else{
- close_env_it->second += 1;
- }
- }
- }
- // score it
- counter = std::min(counter, max_count_);
- score += energies_[this->EIdx(bb_list.GetAA(i), counter)];
- }
-
- // We have now all energies of the residues in the BackboneList...
- // The problem is, that the residues closeby also get affected
- // by the positions of bb_list. The effect gets honoured by calculating
- // a delta in energy for those residue between the energy they would
- // have in the environment and the new energy, given the positions of the
- // bb_list
- Real diff_score = 0.0;
- uint old_env_counter, new_env_counter;
+ for (uint i = 0; i < indices.size(); ++i) {
+ uint my_idx = indices[i];
- for(close_env_it = close_env_residues.begin();
- close_env_it != close_env_residues.end(); ++close_env_it) {
+ if(!env_set_[my_idx]) {
+ throw promod3::Error("Must set environment before scoring!");
+ }
- const CBetaSpatialOrganizerItem& env_data =
- env_->GetEnvData(close_env_it->first);
+ ost::conop::AminoAcid my_aa = env_->GetEnvData(my_idx).aa;
- within_result = env_->FindWithin(env_data.pos, cutoff_);
+ // count neighbours
+ uint count = 0;
+ within_result = env_->FindWithin(env_->GetEnvData(my_idx).pos, cutoff_);
a = within_result.first;
-
- // the number of neighbours before looking at the input bb_list
- old_env_counter = within_result.second - 1; // -1, because the position
- // itself also gets found...
-
- // the number of neighbours when bb_list would be added to env
- new_env_counter = close_env_it->second; // the number of close residues
- // it has from the bb_list
-
- // we still have to add all the residues from the environment not
- // part of the bb_list
- for(uint j = 0; j < within_result.second; ++j) {
- env_idx = a[j].first->idx;
- if(env_idx < start_idx || env_idx > end_idx) ++new_env_counter;
+ for (uint j = 0; j < within_result.second; ++j) {
+ if (my_idx != a[j].first->idx) {
+ ++count;
+ }
}
- new_env_counter -= 1;// The position itself will also be found
-
- new_env_counter = std::min(new_env_counter, max_count_);
- old_env_counter = std::min(old_env_counter, max_count_);
- diff_score += energies_[this->EIdx(env_data.aa,new_env_counter)];
- diff_score -= energies_[this->EIdx(env_data.aa,old_env_counter)];
- }
-
- score += diff_score;
-
- if(normalize_ && bb_list_size) {
- score /= bb_list_size;
+ // score it
+ count = std::min(count, max_count_);
+ scores[i] = energies_[this->EIdx(my_aa, count)];
}
-
- return score;
}
void CBPackingScorer::AttachEnvironment(BackboneScoreEnv& env) {
env_ = env.UniqueListener<CBetaEnvListener>("CBetaEnvListener");
idx_handler_ = env.GetIdxHandlerData();
+ env_set_ = env.GetEnvSetData();
attached_environment_ = true;
}
diff --git a/scoring/src/cb_packing_score.hh b/scoring/src/cb_packing_score.hh
index 1aaf510b9bc04a9f729d6b74f78340393d075b5d..51041352820fca2032ad87b63c6c6b7f282d51e7 100644
--- a/scoring/src/cb_packing_score.hh
+++ b/scoring/src/cb_packing_score.hh
@@ -30,17 +30,13 @@ public:
void SetEnergy(ost::conop::AminoAcid aa, uint count, Real e);
- virtual Real CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const;
+ virtual Real CalculateScore(uint start_resnum, uint num_residues,
+ uint chain_idx) const;
- void CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
+ void CalculateScoreProfile(uint start_resnum, uint num_residues,
+ uint chain_idx,
std::vector<Real>& profile) const;
- void UseClassicMode() { include_env_mode_ = false; }
-
- void UseIncludeEnvMode() { include_env_mode_ = true; }
-
void DoNormalize(bool do_it) {normalize_ = do_it; }
void AttachEnvironment(BackboneScoreEnv& env);
@@ -48,7 +44,6 @@ public:
private:
CBPackingScorer(): attached_environment_(false),
- include_env_mode_(false),
normalize_(true),
env_(NULL) { }
@@ -56,14 +51,9 @@ private:
return aa * (max_count_ + 1) + count;
}
- void Score(const loop::BackboneList& bb_list,
- uint start_resnum,
- uint chain_idx,
+ void Score(const std::vector<uint>& indices,
std::vector<Real>& scores) const;
- Real CalculateScoreIncludeEnv(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const;
-
// portable serialization (only stores data relevant for Load/Save!)
// (cleanly element by element with fixed-width base-types)
template <typename DS>
@@ -88,9 +78,9 @@ private:
//environment specific information, that will be set upon calling
//AttachEnvironment
bool attached_environment_;
- bool include_env_mode_;
bool normalize_;
CBetaEnvListener* env_;
+ const int* env_set_;
const loop::IdxHandler* idx_handler_;
};
diff --git a/scoring/src/cbeta_score.cc b/scoring/src/cbeta_score.cc
index b1d0917619592360779aef788baf1c82fc5357f8..6c0481eb0e66e0d5915b27639cdcda01189284a7 100644
--- a/scoring/src/cbeta_score.cc
+++ b/scoring/src/cbeta_score.cc
@@ -185,8 +185,8 @@ void CBetaScorer::SetEnergy(ost::conop::AminoAcid a, ost::conop::AminoAcid b,
energies_[this->EIdx(b, a, bin)] = e;
}
-Real CBetaScorer::CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const {
+Real CBetaScorer::CalculateScore(uint start_resnum, uint num_residues,
+ uint chain_idx) const {
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
"CBetaScorer::CalculateScore", 2);
@@ -194,7 +194,10 @@ Real CBetaScorer::CalculateScore(const loop::BackboneList& bb_list,
// get scores
std::vector<Real> scores;
std::vector<Real> internal_scores;
- this->Score(bb_list, start_resnum, chain_idx, scores, internal_scores);
+ std::vector<uint> idx_vector =
+ idx_handler_->ToIdxVector(start_resnum, num_residues, chain_idx);
+
+ this->Score(idx_vector, scores, internal_scores);
// sum them
Real score = 0.0;
Real internal_score = 0.0;
@@ -213,8 +216,8 @@ Real CBetaScorer::CalculateScore(const loop::BackboneList& bb_list,
return score;
}
-void CBetaScorer::CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
+void CBetaScorer::CalculateScoreProfile(uint start_resnum, uint num_residues,
+ uint chain_idx,
std::vector<Real>& profile) const {
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
@@ -222,7 +225,10 @@ void CBetaScorer::CalculateScoreProfile(const loop::BackboneList& bb_list,
std::vector<Real> scores;
std::vector<Real> internal_scores;
- this->Score(bb_list, start_resnum, chain_idx, scores, internal_scores);
+ std::vector<uint> idx_vector =
+ idx_handler_->ToIdxVector(start_resnum, num_residues, chain_idx);
+
+ this->Score(idx_vector, scores, internal_scores);
profile.resize(scores.size());
for(uint i = 0; i < scores.size(); ++i) {
@@ -231,10 +237,8 @@ void CBetaScorer::CalculateScoreProfile(const loop::BackboneList& bb_list,
}
-void CBetaScorer::Score(const loop::BackboneList& bb_list,
- uint start_resnum,
- uint chain_idx,
- std::vector<Real>& scores,
+void CBetaScorer::Score(const std::vector<uint>& indices,
+ std::vector<Real>& external_scores,
std::vector<Real>& internal_scores) const {
@@ -246,69 +250,56 @@ void CBetaScorer::Score(const loop::BackboneList& bb_list,
String err = "Can only calculate score when environment is attached!";
throw promod3::Error(err);
}
- uint bb_list_size, start_idx, end_idx;
- idx_handler_->SetupScoreCalculation(bb_list, start_resnum, chain_idx,
- bb_list_size, start_idx, end_idx);
- scores.assign(bb_list_size, 0.0);
- internal_scores.assign(bb_list_size, 0.0);
+ uint n_indices = indices.size();
+
+ external_scores.assign(n_indices, 0.0);
+ internal_scores.assign(n_indices, 0.0);
+
+ if(n_indices == 0) return;
- if(bb_list_size == 0) return;
+ // we need to store which interactions are internal and which are external
+ std::vector<bool> occupied(idx_handler_->GetNumResidues(), false);
+ for(std::vector<uint>::const_iterator i = indices.begin();
+ i != indices.end(); ++i) {
+ occupied[*i] = true;
+ }
const Real inv_bin_size = bins_ / cutoff_;
// let's reduce the cutoff a bit -> hack to ensure bin always within range
- const Real squared_cutoff = cutoff_*cutoff_ - 0.0001;
-
- if(do_external_scores_) {
-
- const Real cutoff = cutoff_ - 0.0001;
-
- // let's precalculate the idx ranges we have to neglect
- int min_idx_chain = idx_handler_->GetChainStartIdx(chain_idx);
- int max_idx_chain = min_idx_chain +
- idx_handler_->GetChainSize(chain_idx) - 1;
- std::vector<int> start_neglect(bb_list_size, start_idx);
- std::vector<int> end_neglect(bb_list_size, end_idx);
- int neglect_range = int(seq_sep_) - 1;
-
- int max_i = std::min(neglect_range, static_cast<int>(bb_list_size));
- for (int i = 0; i < max_i; ++i) {
- start_neglect[i] =
- std::max(min_idx_chain, start_neglect[i]-(neglect_range-i));
- end_neglect[bb_list_size-1-i] =
- std::min(max_idx_chain, end_neglect[i]+(neglect_range-i));
+ const Real cutoff = cutoff_ - 0.0001;
+
+ // let's go over every loop residue
+ CBetaSpatialOrganizer::WithinResult within_result;
+ std::pair<CBetaSpatialOrganizerItem*,Real>* a;
+ for (uint i = 0; i < n_indices; ++i) {
+ int my_idx = indices[i];
+
+ if(!env_set_[my_idx]) {
+ throw promod3::Error("Must set environment before scoring!");
}
- // let's go over every loop residue
- CBetaSpatialOrganizer::WithinResult within_result;
- std::pair<CBetaSpatialOrganizerItem*,Real>* a;
- for (uint i = 0; i < bb_list_size; ++i) {
- within_result = env_->FindWithin(bb_list.GetCB(i), cutoff);
- a = within_result.first;
- const ost::conop::AminoAcid myaa = bb_list.GetAA(i);
- for (uint j = 0; j < within_result.second; ++j) {
- const int my_idx = int(a[j].first->idx);
- if (my_idx < start_neglect[i] || my_idx > end_neglect[i]) {
- const uint bin = a[j].second * inv_bin_size;
- scores[i] += energies_[this->EIdx(myaa, a[j].first->aa, bin)];
+ int my_chain_idx = idx_handler_->GetChainIdx(my_idx);
+ ost::conop::AminoAcid my_aa = env_->GetEnvData(my_idx).aa;
+
+ within_result = env_->FindWithin(env_->GetEnvData(my_idx).pos, cutoff);
+ a = within_result.first;
+
+ for (uint j = 0; j < within_result.second; ++j) {
+ int other_idx = int(a[j].first->idx);
+ int other_chain_idx = idx_handler_->GetChainIdx(other_idx);
+
+ if(my_chain_idx != other_chain_idx ||
+ std::abs(my_idx - other_idx) >= seq_sep_) {
+ uint bin = a[j].second * inv_bin_size;
+ Real e = energies_[this->EIdx(my_aa, a[j].first->aa, bin)];
+
+ if(occupied[other_idx] && do_internal_scores_) {
+ internal_scores[i] += e;
}
- }
- }
- }
-
- if(do_internal_scores_) {
- if (bb_list_size >= seq_sep_) {
- for (uint i = 0; i < bb_list_size-seq_sep_; ++i) {
- const ost::conop::AminoAcid myaa = bb_list.GetAA(i);
- for (uint j = i+seq_sep_; j < bb_list_size; ++j) {
- const Real dist2 = geom::Length2(bb_list.GetCB(i) -
- bb_list.GetCB(j));
- if (dist2 < squared_cutoff) {
- const uint bin = std::sqrt(dist2) * inv_bin_size;
- const ost::conop::AminoAcid aaj = bb_list.GetAA(j);
- internal_scores[i] += energies_[this->EIdx(myaa, aaj, bin)];
- internal_scores[j] += energies_[this->EIdx(aaj, myaa, bin)];
- }
+
+ if(!occupied[other_idx] && do_external_scores_) {
+ external_scores[i] += e;
}
}
}
@@ -319,6 +310,7 @@ void CBetaScorer::Score(const loop::BackboneList& bb_list,
void CBetaScorer::AttachEnvironment(BackboneScoreEnv& env) {
env_ = env.UniqueListener<CBetaEnvListener>("CBetaEnvListener");
idx_handler_ = env.GetIdxHandlerData();
+ env_set_ = env.GetEnvSetData();
attached_environment_ = true;
}
diff --git a/scoring/src/cbeta_score.hh b/scoring/src/cbeta_score.hh
index ce7643900c9686e9fe44d2bc262761c6e21fe150..589b2078a02af56cd7052d2ac29ce38f2c6fa526 100644
--- a/scoring/src/cbeta_score.hh
+++ b/scoring/src/cbeta_score.hh
@@ -37,11 +37,11 @@ public:
void DoNormalize(bool do_it) { normalize_ = do_it; }
- virtual Real CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const;
+ virtual Real CalculateScore(uint start_resnum, uint num_residues,
+ uint chain_idx) const;
- virtual void CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum,
+ virtual void CalculateScoreProfile(uint start_resnum,
+ uint num_residues,
uint chain_idx,
std::vector<Real>& profile) const;
@@ -57,10 +57,8 @@ private:
return (a * ost::conop::XXX + b) * bins_ + bin;
}
- void Score(const loop::BackboneList& bb_list,
- uint start_resnum,
- uint chain_idx,
- std::vector<Real>& scores,
+ void Score(const std::vector<uint>& indices,
+ std::vector<Real>& external_scores,
std::vector<Real>& internal_scores) const;
// portable serialization (only stores data relevant for Load/Save!)
@@ -89,6 +87,7 @@ private:
//AttachEnvironment
bool attached_environment_;
CBetaEnvListener* env_;
+ const int* env_set_;
const loop::IdxHandler* idx_handler_;
};
diff --git a/scoring/src/clash_env_listener.cc b/scoring/src/clash_env_listener.cc
index 12d56562d8f5f09b5ea0852e55c698436b4046d0..37fa788768092c5b52a46ac9206600872fca198a 100644
--- a/scoring/src/clash_env_listener.cc
+++ b/scoring/src/clash_env_listener.cc
@@ -7,7 +7,10 @@ ClashEnvListener::ClashEnvListener(): env_(8.0),
env_data_(NULL) { }
ClashEnvListener::~ClashEnvListener() {
- if (env_data_ != NULL) delete [] env_data_;
+ if (env_data_ != NULL) {
+ delete [] is_glycine_;
+ delete [] env_data_;
+ }
}
void ClashEnvListener::Init(const BackboneScoreEnv& base_env) {
@@ -18,14 +21,17 @@ void ClashEnvListener::Init(const BackboneScoreEnv& base_env) {
if (env_data_ != NULL) {
// the env listener has already been initialized at some point...
// let's clean up first...
+ delete [] is_glycine_;
delete [] env_data_;
env_data_ = NULL;
env_.Clear();
}
+ const ost::conop::AminoAcid* aa_data = base_env.GetAAData();
const loop::IdxHandler* idx_handler = base_env.GetIdxHandlerData();
uint num_residues = idx_handler->GetNumResidues();
env_data_ = new ClashSpatialOrganizerItem[5*num_residues];
+ is_glycine_ = new bool[num_residues];
for (uint i = 0; i < num_residues; ++i) {
// N
@@ -43,6 +49,13 @@ void ClashEnvListener::Init(const BackboneScoreEnv& base_env) {
// CB
env_data_[5*i+4].clash_type = BB_CLASH_C;
env_data_[5*i+4].idx = i;
+
+ if(aa_data[i] == ost::conop::GLY) {
+ is_glycine_[i] = true;
+ }
+ else {
+ is_glycine_[i] = false;
+ }
}
}
diff --git a/scoring/src/clash_env_listener.hh b/scoring/src/clash_env_listener.hh
index 6dce4b0112a1df4bbe86c10354aae8270a9d4f93..85761433ace7cf0695d29ec6ada1047c1bc3c1b2 100644
--- a/scoring/src/clash_env_listener.hh
+++ b/scoring/src/clash_env_listener.hh
@@ -52,10 +52,16 @@ public:
return env_.FindWithin(pos, cutoff);
}
+ // for internal use in the actual scorer to access required per residue data
+ // NO RANGE CHECK!
+ ClashSpatialOrganizerItem* GetParticleData(uint idx) { return &env_data_[5*idx]; };
+ bool IsGlycine(uint idx) { return is_glycine_[idx]; }
+
private:
ClashSpatialOrganizer env_;
// here: 5 atoms stored for each index: use env_data_[5*idx + sub_idx]
+ bool* is_glycine_;
ClashSpatialOrganizerItem* env_data_;
};
diff --git a/scoring/src/clash_score.cc b/scoring/src/clash_score.cc
index c90a29a0b6228096ba336b49ac9de6153f90dfa5..b5909713fa354f915e7d0b1ce4ae2390478830dd 100644
--- a/scoring/src/clash_score.cc
+++ b/scoring/src/clash_score.cc
@@ -43,16 +43,19 @@ ClashScorer::~ClashScorer() {
}
-Real ClashScorer::CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const {
+Real ClashScorer::CalculateScore(uint start_resnum, uint num_residues,
+ uint chain_idx) const {
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
"ClashScorer::CalculateScore", 2);
+ // get scores
std::vector<Real> scores;
std::vector<Real> internal_scores;
+ std::vector<uint> idx_vector =
+ idx_handler_->ToIdxVector(start_resnum, num_residues, chain_idx);
- this->Score(bb_list, start_resnum, chain_idx, scores, internal_scores);
+ this->Score(idx_vector, scores, internal_scores);
Real score = 0.0;
Real internal_score = 0.0;
@@ -72,8 +75,8 @@ Real ClashScorer::CalculateScore(const loop::BackboneList& bb_list,
}
-void ClashScorer::CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
+void ClashScorer::CalculateScoreProfile(uint start_resnum, uint num_residues,
+ uint chain_idx,
std::vector<Real>& profile) const {
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
@@ -81,8 +84,10 @@ void ClashScorer::CalculateScoreProfile(const loop::BackboneList& bb_list,
std::vector<Real> scores;
std::vector<Real> internal_scores;
+ std::vector<uint> idx_vector =
+ idx_handler_->ToIdxVector(start_resnum, num_residues, chain_idx);
- this->Score(bb_list, start_resnum, chain_idx, scores, internal_scores);
+ this->Score(idx_vector, scores, internal_scores);
profile.resize(scores.size());
for(uint i = 0; i < scores.size(); ++i) {
@@ -91,9 +96,8 @@ void ClashScorer::CalculateScoreProfile(const loop::BackboneList& bb_list,
}
-void ClashScorer::Score(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
- std::vector<Real>& scores,
+void ClashScorer::Score(const std::vector<uint>& indices,
+ std::vector<Real>& external_scores,
std::vector<Real>& internal_scores) const {
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
@@ -104,266 +108,57 @@ void ClashScorer::Score(const loop::BackboneList& bb_list,
String err = "Can only calculate score when environment is attached!";
throw promod3::Error(err);
}
- uint bb_list_size, start_idx, end_idx;
- idx_handler_->SetupScoreCalculation(bb_list, start_resnum, chain_idx,
- bb_list_size, start_idx, end_idx);
-
- // setup result
- scores.assign(bb_list_size, 0.0);
- internal_scores.assign(bb_list_size, 0.0);
- if (bb_list_size == 0) return;
+
+ uint n_indices = indices.size();
+
+ external_scores.assign(n_indices, 0.0);
+ internal_scores.assign(n_indices, 0.0);
+
+ if(n_indices == 0) return;
+
+ // we need to store which interactions are internal and which are external
+ std::vector<bool> occupied(idx_handler_->GetNumResidues(), false);
+ for(std::vector<uint>::const_iterator i = indices.begin();
+ i != indices.end(); ++i) {
+ occupied[*i] = true;
+ }
+
uint seq_sep = 2;
uint bin;
+ Real score;
- if(do_external_scores_) {
- // let's precalculate the idx ranges we have to neglect
- int min_idx_chain = idx_handler_->GetChainStartIdx(chain_idx);
- int max_idx_chain = min_idx_chain +
- idx_handler_->GetChainSize(chain_idx) - 1;
- std::vector<int> start_neglect(bb_list_size, start_idx);
- std::vector<int> end_neglect(bb_list_size, end_idx);
- // we use a fixed seq. separation of 2 residues
- int neglect_range = int(seq_sep) - 1;
- int max_i = std::min(neglect_range, static_cast<int>(bb_list_size));
- for (int i = 0; i < max_i; ++i) {
- start_neglect[i] =
- std::max(min_idx_chain, start_neglect[i]-(neglect_range-i));
- end_neglect[bb_list_size-1-i] =
- std::min(max_idx_chain, end_neglect[i]+(neglect_range-i));
- }
+ // let's go over every loop residue
+ ClashSpatialOrganizer::WithinResult within_result;
+ std::pair<ClashSpatialOrganizerItem*,Real>* a;
- // let's go over every loop residue
- ClashSpatialOrganizer::WithinResult within_result;
- std::pair<ClashSpatialOrganizerItem*,Real>* a;
- for (uint i = 0; i < bb_list_size; ++i) {
- // N
- within_result = env_->FindWithin(bb_list.GetN(i), Real(2.9));
- a = within_result.first;
- for (uint j = 0; j < within_result.second; ++j) {
- const int my_idx = int(a[j].first->idx);
- if (my_idx < start_neglect[i] || my_idx > end_neglect[i]) {
- bin = a[j].second * bins_per_a_;
- scores[i] += clash_scores_[BB_CLASH_N][a[j].first->clash_type][bin];
- }
- }
- // CA
- within_result = env_->FindWithin(bb_list.GetCA(i), Real(3.2));
- a = within_result.first;
- for (uint j = 0; j < within_result.second; ++j) {
- const int my_idx = int(a[j].first->idx);
- if (my_idx < start_neglect[i] || my_idx > end_neglect[i]) {
- bin = a[j].second * bins_per_a_;
- scores[i] += clash_scores_[BB_CLASH_C][a[j].first->clash_type][bin];
- }
- }
- // C
- within_result = env_->FindWithin(bb_list.GetC(i), Real(3.2));
- a = within_result.first;
- for (uint j = 0; j < within_result.second; ++j) {
- const int my_idx = int(a[j].first->idx);
- if (my_idx < start_neglect[i] || my_idx > end_neglect[i]) {
- bin = a[j].second * bins_per_a_;
- scores[i] += clash_scores_[BB_CLASH_C][a[j].first->clash_type][bin];
- }
- }
- // O
- within_result = env_->FindWithin(bb_list.GetO(i), Real(2.9));
- a = within_result.first;
- for (uint j = 0; j < within_result.second; ++j) {
- const int my_idx = int(a[j].first->idx);
- if (my_idx < start_neglect[i] || my_idx > end_neglect[i]) {
- bin = a[j].second * bins_per_a_;
- scores[i] += clash_scores_[BB_CLASH_O][a[j].first->clash_type][bin];
- }
- }
- // CB (if needed)
- if (bb_list.GetAA(i) != ost::conop::GLY) {
- within_result = env_->FindWithin(bb_list.GetCB(i), Real(3.2));
- a = within_result.first;
- for (uint j = 0; j < within_result.second; ++j) {
- const int my_idx = int(a[j].first->idx);
- if (my_idx < start_neglect[i] || my_idx > end_neglect[i]) {
- bin = a[j].second * bins_per_a_;
- scores[i] += clash_scores_[BB_CLASH_C][a[j].first->clash_type][bin];
- }
- }
- }
+ for (uint i = 0; i < n_indices; ++i) {
+ int my_idx = indices[i];
+
+ if(!env_set_[my_idx]) {
+ throw promod3::Error("Must set environment before scoring!");
}
- }
- // skip stuff below if we dont do internal scores or size is too small
- if (!do_internal_scores_ || bb_list_size <= seq_sep) {
- return;
- }
+ uint n_particles = env_->IsGlycine(my_idx) ? 4 : 5;
+ ClashSpatialOrganizerItem* particle_data = env_->GetParticleData(my_idx);
- // get internal scores -> N^2 loop
- Real clash_score;
- Real d;
- for (uint i = 0; i < bb_list_size-seq_sep; ++i) {
- // do not consider neighbours
- for (uint j = i+seq_sep; j < bb_list_size; ++j) {
- // skip this residue-residue interaction if ca-distance is above 7A
- d = geom::Length2(bb_list.GetCA(i) - bb_list.GetCA(j));
- if (d > Real(49)) continue;
-
- bin = std::sqrt(d) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_C][BB_CLASH_C][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- bin = geom::Distance(bb_list.GetCA(i), bb_list.GetN(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_C][BB_CLASH_N][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- bin = geom::Distance(bb_list.GetCA(i), bb_list.GetC(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_C][BB_CLASH_C][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- bin = geom::Distance(bb_list.GetCA(i), bb_list.GetO(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_C][BB_CLASH_O][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
-
- bin = geom::Distance(bb_list.GetN(i), bb_list.GetN(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_N][BB_CLASH_N][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- bin = geom::Distance(bb_list.GetN(i), bb_list.GetCA(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_N][BB_CLASH_C][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- bin = geom::Distance(bb_list.GetN(i), bb_list.GetC(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_N][BB_CLASH_C][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- bin = geom::Distance(bb_list.GetN(i), bb_list.GetO(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_N][BB_CLASH_O][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
+ for(uint j = 0; j < n_particles; ++j) {
+ int clash_type = particle_data[j].clash_type;
+ within_result = env_->FindWithin(particle_data[j].pos, 3.2);
+ a = within_result.first;
- bin = geom::Distance(bb_list.GetC(i), bb_list.GetN(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_C][BB_CLASH_N][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- bin = geom::Distance(bb_list.GetC(i), bb_list.GetCA(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_C][BB_CLASH_C][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- bin = geom::Distance(bb_list.GetC(i), bb_list.GetC(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_C][BB_CLASH_C][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- bin = geom::Distance(bb_list.GetC(i), bb_list.GetO(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_C][BB_CLASH_O][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
+ for(uint k = 0; k < within_result.second; ++k) {
+ int other_idx = int(a[k].first->idx);
- bin = geom::Distance(bb_list.GetO(i), bb_list.GetN(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_O][BB_CLASH_N][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- bin = geom::Distance(bb_list.GetO(i), bb_list.GetCA(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_O][BB_CLASH_C][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- bin = geom::Distance(bb_list.GetO(i), bb_list.GetC(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_O][BB_CLASH_C][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- bin = geom::Distance(bb_list.GetO(i), bb_list.GetO(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_O][BB_CLASH_O][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
+ if(std::abs(my_idx - other_idx) >= seq_sep) {
+ bin = a[k].second * bins_per_a_;
+ score = clash_scores_[clash_type][a[k].first->clash_type][bin];
- if (bb_list.GetAA(i) != ost::conop::GLY) {
- bin = geom::Distance(bb_list.GetCB(i), bb_list.GetN(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_C][BB_CLASH_N][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- bin = geom::Distance(bb_list.GetCB(i), bb_list.GetCA(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_C][BB_CLASH_C][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- bin = geom::Distance(bb_list.GetCB(i), bb_list.GetC(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_C][BB_CLASH_C][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- bin = geom::Distance(bb_list.GetCB(i), bb_list.GetO(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_C][BB_CLASH_O][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- }
+ if(occupied[other_idx] && do_internal_scores_) {
+ internal_scores[i] += score;
+ }
- if (bb_list.GetAA(j) != ost::conop::GLY) {
- bin = geom::Distance(bb_list.GetN(i), bb_list.GetCB(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_N][BB_CLASH_C][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- bin = geom::Distance(bb_list.GetCA(i), bb_list.GetCB(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_C][BB_CLASH_C][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- bin = geom::Distance(bb_list.GetC(i), bb_list.GetCB(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_C][BB_CLASH_C][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- bin = geom::Distance(bb_list.GetO(i), bb_list.GetCB(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_O][BB_CLASH_C][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
- }
- if (bb_list.GetAA(i) != ost::conop::GLY) {
- bin = geom::Distance(bb_list.GetCB(i), bb_list.GetCB(j)) * bins_per_a_;
- if(bin < bins_) {
- clash_score = clash_scores_[BB_CLASH_C][BB_CLASH_C][bin];
- internal_scores[i] += clash_score;
- internal_scores[j] += clash_score;
+ if(!occupied[other_idx] && do_external_scores_) {
+ external_scores[i] += score;
}
}
}
@@ -374,6 +169,7 @@ void ClashScorer::Score(const loop::BackboneList& bb_list,
void ClashScorer::AttachEnvironment(BackboneScoreEnv& env) {
env_ = env.UniqueListener<ClashEnvListener>("ClashEnvListener");
idx_handler_ = env.GetIdxHandlerData();
+ env_set_ = env.GetEnvSetData();
attached_environment_ = true;
}
diff --git a/scoring/src/clash_score.hh b/scoring/src/clash_score.hh
index 1f1ad1c7ac109623cf746660fd1ec3ed8190545b..71b5d7405f97c3f4b19f4e4009d24066095b145d 100644
--- a/scoring/src/clash_score.hh
+++ b/scoring/src/clash_score.hh
@@ -25,12 +25,12 @@ public:
void DoNormalize(bool do_it) { normalize_ = do_it; }
- virtual Real CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum,
+ virtual Real CalculateScore(uint start_resnum,
+ uint num_residues,
uint chain_idx) const;
- virtual void CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum,
+ virtual void CalculateScoreProfile(uint start_resnum,
+ uint num_residues,
uint chain_idx,
std::vector<Real>& profile) const;
@@ -38,10 +38,8 @@ public:
private:
- void Score(const loop::BackboneList& bb_list,
- uint start_resnum,
- uint chain_idx,
- std::vector<Real>& scores,
+ void Score(const std::vector<uint>& indices,
+ std::vector<Real>& external_scores,
std::vector<Real>& internal_scores) const;
// hard coded bounds (safe for dist. < 3.21A)
@@ -51,10 +49,12 @@ private:
bool do_internal_scores_;
bool do_external_scores_;
bool normalize_;
+
//environment specific information, that will be set upon calling
//AttachEnvironment
bool attached_environment_;
ClashEnvListener* env_;
+ const int* env_set_;
const loop::IdxHandler* idx_handler_;
};
diff --git a/scoring/src/density_score.cc b/scoring/src/density_score.cc
index 8a8c77aa7816db5afcce9aaa785b713f1e0d2bcb..803a843b68e1dd8396d27b3cd67ef916c5232649 100644
--- a/scoring/src/density_score.cc
+++ b/scoring/src/density_score.cc
@@ -2,28 +2,14 @@
namespace promod3{ namespace scoring{
-Real DensityScorer::CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const {
- promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
- "DensityScorer::CalculateScore", 2);
-
-
- if (!attached_environment_) {
- String err = "Can only calculate score when environment is attached!";
- throw promod3::Error(err);
- }
-
- if(!map_data_->map.IsValid()){
- String err = "Environment attached to density scorer must have a valid ";
- err += "density map for score calculation!";
- throw promod3::Error(err);
- }
+Real DensityScore(const loop::BackboneList& bb_list, ost::img::MapHandle& map,
+ Real resolution, bool all_atom) {
geom::AlignedCuboid bounding_box = bb_list.GetBounds(true);
geom::Vec3 min = bounding_box.GetMin();
geom::Vec3 max = bounding_box.GetMax();
- geom::Vec3 sampling = map_data_->map.GetSpatialSampling();
+ geom::Vec3 sampling = map.GetSpatialSampling();
//thats the extend of the image we're actually interested in
uint x_extent = std::ceil((max[0] - min[0]) / sampling[0]);
@@ -31,9 +17,9 @@ Real DensityScorer::CalculateScore(const loop::BackboneList& bb_list,
uint z_extent = std::ceil((max[2] - min[2]) / sampling[2]);
//let's enlarge a bit...
- Real x_overlap = Real(2 * map_data_->resolution) / sampling[0];
- Real y_overlap = Real(2 * map_data_->resolution) / sampling[1];
- Real z_overlap = Real(2 * map_data_->resolution) / sampling[2];
+ Real x_overlap = Real(2 * resolution) / sampling[0];
+ Real y_overlap = Real(2 * resolution) / sampling[1];
+ Real z_overlap = Real(2 * resolution) / sampling[2];
x_extent += static_cast<uint>(2 * x_overlap);
y_extent += static_cast<uint>(2 * y_overlap);
z_extent += static_cast<uint>(2 * z_overlap);
@@ -46,20 +32,19 @@ Real DensityScorer::CalculateScore(const loop::BackboneList& bb_list,
ost::img::MapHandle bb_list_map = ost::img::CreateImage(s);
bb_list_map.SetSpatialSampling(sampling);
- geom::Vec3 pixel_coord = map_data_->map.CoordToIndex(min -
- geom::Vec3(x_overlap,
- y_overlap,
- z_overlap));
+ geom::Vec3 pixel_coord = map.CoordToIndex(min - geom::Vec3(x_overlap,
+ y_overlap,
+ z_overlap));
ost::img::Point rounded_pixel_coord(round(pixel_coord[0]),
round(pixel_coord[1]),
round(pixel_coord[2]));
- min = map_data_->map.IndexToCoord(rounded_pixel_coord);
+ min = map.IndexToCoord(rounded_pixel_coord);
bb_list_map.SetAbsoluteOrigin(min);
//let's fill the density of the BackboneList
- bb_list.InsertInto(bb_list_map, map_data_->resolution, map_data_->all_atom);
+ bb_list.InsertInto(bb_list_map, resolution, all_atom);
uint num_voxels = (s[0] * s[1] * s[2]);
@@ -74,164 +59,13 @@ Real DensityScorer::CalculateScore(const loop::BackboneList& bb_list,
p = ost::img::Point(i,j,k);
p_orig = p + rounded_pixel_coord;
//will return 0.0 if outside of map
- map_values[vec_position] = map_data_->map.GetReal(p_orig);
+ map_values[vec_position] = map.GetReal(p_orig);
bb_list_values[vec_position] = bb_list_map.GetReal(p);
++vec_position;
}
}
}
- return this->Correlate(map_values, bb_list_values, num_voxels);
-}
-
-
-void DensityScorer::CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
- std::vector<Real>& profile) const{
-
- promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
- "DensityScorer::CalculateScore", 2);
-
-
- if (!attached_environment_) {
- String err = "Can only calculate score when environment is attached!";
- throw promod3::Error(err);
- }
-
- if(!map_data_->map.IsValid()){
- String err = "Environment attached to density scorer must have a valid ";
- err += "density map for score calculation!";
- throw promod3::Error(err);
- }
-
- geom::Vec3 sampling = map_data_->map.GetSpatialSampling();
- Real x_overlap = Real(2 * map_data_->resolution) / sampling[0];
- Real y_overlap = Real(2 * map_data_->resolution) / sampling[1];
- Real z_overlap = Real(2 * map_data_->resolution) / sampling[2];
-
- profile.clear();
- for(uint res_idx = 0; res_idx < bb_list.size(); ++res_idx){
-
- //we have to calculate the bounding box for a single residue manually
-
- Real min_x = std::numeric_limits<Real>::max();
- Real max_x = -std::numeric_limits<Real>::max();
- Real min_y = std::numeric_limits<Real>::max();
- Real max_y = -std::numeric_limits<Real>::max();
- Real min_z = std::numeric_limits<Real>::max();
- Real max_z = -std::numeric_limits<Real>::max();
-
- geom::Vec3 pos;
- pos = bb_list.GetN(res_idx);
- min_x = std::min(pos[0], min_x);
- max_x = std::max(pos[0], max_x);
- min_y = std::min(pos[1], min_y);
- max_y = std::max(pos[1], max_y);
- min_z = std::min(pos[2], min_z);
- max_z = std::max(pos[2], max_z);
-
- pos = bb_list.GetCA(res_idx);
- min_x = std::min(pos[0], min_x);
- max_x = std::max(pos[0], max_x);
- min_y = std::min(pos[1], min_y);
- max_y = std::max(pos[1], max_y);
- min_z = std::min(pos[2], min_z);
- max_z = std::max(pos[2], max_z);
-
- if (bb_list.GetAA(res_idx) != ost::conop::GLY) {
- pos = bb_list.GetCB(res_idx);
- min_x = std::min(pos[0], min_x);
- max_x = std::max(pos[0], max_x);
- min_y = std::min(pos[1], min_y);
- max_y = std::max(pos[1], max_y);
- min_z = std::min(pos[2], min_z);
- max_z = std::max(pos[2], max_z);
- }
-
- pos = bb_list.GetC(res_idx);
- min_x = std::min(pos[0], min_x);
- max_x = std::max(pos[0], max_x);
- min_y = std::min(pos[1], min_y);
- max_y = std::max(pos[1], max_y);
- min_z = std::min(pos[2], min_z);
- max_z = std::max(pos[2], max_z);
-
- pos = bb_list.GetO(res_idx);
- min_x = std::min(pos[0], min_x);
- max_x = std::max(pos[0], max_x);
- min_y = std::min(pos[1], min_y);
- max_y = std::max(pos[1], max_y);
- min_z = std::min(pos[2], min_z);
- max_z = std::max(pos[2], max_z);
-
- //thats the extend of the image we're actually interested in
- uint x_extent = std::ceil((max_x - min_x) / sampling[0]);
- uint y_extent = std::ceil((max_y - min_y) / sampling[1]);
- uint z_extent = std::ceil((max_z - min_z) / sampling[2]);
-
- //let's enlarge a bit...
- x_extent += static_cast<uint>(2 * x_overlap);
- y_extent += static_cast<uint>(2 * y_overlap);
- z_extent += static_cast<uint>(2 * z_overlap);
-
- //let's create a new image, where we will fill the density of the bb_list
- //we make sure, that the grid cubes exactly match the ones in the
- //set map
- ost::img::Size s(x_extent,y_extent,z_extent);
-
- ost::img::MapHandle bb_list_map = ost::img::CreateImage(s);
- bb_list_map.SetSpatialSampling(sampling);
-
- geom::Vec3 pixel_coord = map_data_->map.CoordToIndex(geom::Vec3(min_x,
- min_y,
- min_z) -
- geom::Vec3(x_overlap,
- y_overlap,
- z_overlap));
-
- ost::img::Point rounded_pixel_coord(round(pixel_coord[0]),
- round(pixel_coord[1]),
- round(pixel_coord[2]));
-
- geom::Vec3 min = map_data_->map.IndexToCoord(rounded_pixel_coord);
- bb_list_map.SetAbsoluteOrigin(min);
-
- //let's fill the density of the BackboneList
- bb_list.InsertInto(bb_list_map, map_data_->resolution, map_data_->all_atom);
-
- uint num_voxels = (s[0] * s[1] * s[2]);
-
- std::vector<Real> map_values(num_voxels);
- std::vector<Real> bb_list_values(num_voxels);
-
- uint vec_position = 0;
- ost::img::Point p, p_orig;
- for(uint i = 0; i < s[0]; ++i){
- for(uint j = 0; j < s[1]; ++j){
- for(uint k = 0; k < s[2]; ++k){
- p = ost::img::Point(i,j,k);
- p_orig = p + rounded_pixel_coord;
- //will return 0.0 if outside of map
- map_values[vec_position] = map_data_->map.GetReal(p_orig);
- bb_list_values[vec_position] = bb_list_map.GetReal(p);
- ++vec_position;
- }
- }
- }
-
- profile.push_back(this->Correlate(map_values, bb_list_values, num_voxels));
- }
-}
-
-void DensityScorer::AttachEnvironment(BackboneScoreEnv& env) {
- map_data_ = env.GetDensityMapData();
- attached_environment_ = true;
-}
-
-Real DensityScorer::Correlate(const std::vector<Real>& map_values,
- const std::vector<Real>& bb_list_values,
- uint num_voxels) const{
-
//we have to get the means, the standart deviations and the correlation itself...
Real m_one = 0.0;
Real m_two = 0.0;
diff --git a/scoring/src/density_score.hh b/scoring/src/density_score.hh
index 5a8d7bcb961e3ec7c3c0ab68e23f9dc8fc60817a..d5e5123c3f733f2976e47c8928d32b6d11dc8f57 100644
--- a/scoring/src/density_score.hh
+++ b/scoring/src/density_score.hh
@@ -3,42 +3,14 @@
#include <ost/img/map.hh>
#include <promod3/loop/density_creator.hh>
-#include <promod3/scoring/backbone_score_base.hh>
+#include <promod3/loop/backbone.hh>
#include <promod3/core/runtime_profiling.hh>
namespace promod3 { namespace scoring {
-class DensityScorer;
-typedef boost::shared_ptr<DensityScorer> DensityScorerPtr;
-
-
-class DensityScorer : public BackboneScorer {
-
-public:
-
- DensityScorer(): attached_environment_(false) { }
-
- virtual ~DensityScorer() { }
-
- virtual Real CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const;
-
- void CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
- std::vector<Real>& profile) const;
-
-void AttachEnvironment(BackboneScoreEnv& env);
-
-private:
-
- Real Correlate(const std::vector<Real>& map_values,
- const std::vector<Real>& bb_list_values,
- uint num_voxels) const;
-
- const DensityMapData* map_data_;
- bool attached_environment_;
-};
-
+Real DensityScore(const promod3::loop::BackboneList& bb_list,
+ ost::img::MapHandle& map,
+ Real resolution, bool all_atom);
}}//ns
diff --git a/scoring/src/hbond_env_listener.cc b/scoring/src/hbond_env_listener.cc
index ca4fbd672d9899f6da2332d119d399bc454391bb..e6143ac5eb6a5360dba33d82ab3eeef93ac0fdeb 100644
--- a/scoring/src/hbond_env_listener.cc
+++ b/scoring/src/hbond_env_listener.cc
@@ -1,4 +1,5 @@
#include <promod3/scoring/hbond_env_listener.hh>
+#include <set>
namespace promod3{ namespace scoring{
@@ -95,113 +96,63 @@ void HBondEnvListener::ClearEnvironment(const BackboneScoreEnv& base_env,
this->UpdateStatesAndHPos(idx, base_env);
}
-void HBondEnvListener::UpdateStatesAndHPos(uint chain_idx,
+void HBondEnvListener::UpdateStatesAndHPos(const std::vector<uint>& idx,
const BackboneScoreEnv& env) {
const loop::IdxHandler* idx_handler = env.GetIdxHandlerData();
const int* env_set = env.GetEnvSetData();
- //the h positions and states of the hbond term get calculated as a
- //postprocessing step
- //note, that we reconstruct the stuff for the full chain, since the state of
- //one residue incluences the neighbouring residues...
- //we dont really know how far this influence goes
- Real phi,psi;
- HBondSpatialOrganizerItem* current_item;
- std::vector<int> states(idx_handler->GetChainSize(chain_idx),0);
- int actual_idx = idx_handler->GetChainStartIdx(chain_idx);
-
- for(uint i = 0; i < idx_handler->GetChainSize(chain_idx); ++i){
- if(env_set[actual_idx]){
- current_item = &env_data_[actual_idx];
- phi = -1.0472;
- psi = -0.78540;
- if(i > 0 && i < idx_handler->GetChainSize(chain_idx)-1){
- if(env_set[actual_idx - 1] &&
- env_set[actual_idx + 1]){
-
- phi = geom::DihedralAngle(env_data_[actual_idx-1].c_pos,
- current_item->n_pos,
- current_item->pos,
- current_item->c_pos);
- psi = geom::DihedralAngle(current_item->n_pos,
- current_item->pos,
- current_item->c_pos,
- env_data_[actual_idx+1].n_pos);
-
- states[i] = GetHBondState(phi,psi);
- }
- }
-
- core::ConstructAtomPos(current_item->c_pos, current_item->pos,
- current_item->n_pos, 0.9970, 2.0420,
- phi + Real(M_PI), current_item->h_pos);
- }
- ++actual_idx;
- }
-
- //we still need to check consistency of consecutive residues
- //if state is 1 or 2, at least one neighbouring residues must
- //have the same state, the state is set to zero otherwise
-
- //first do the first and last residue
- if(idx_handler->GetChainSize(chain_idx) > 2){
- if(states[0] == 1){
- if(states[1] != 1) states[0] = 0;
- }
- if(states[0] == 2){
- if(states[1] != 2) states[0] = 0;
+ // we have no guarentee that the indices are consecutive...
+ // we have to do the inefficient way and always add a full triplett
+ // to a set to get a unique set for further processing.
+ std::set<uint> stuff_to_check;
+
+ for(std::vector<uint>::const_iterator i = idx.begin();
+ i != idx.end(); ++i){
+
+ uint chain_idx = idx_handler->GetChainIdx(*i);
+ uint chain_start_idx = idx_handler->GetChainIdx(chain_idx);
+ uint chain_last_idx = idx_handler->GetChainLastIdx(chain_idx);
+
+ if(env_set[*i]) {
+ stuff_to_check.insert(*i);
}
- if(states.back() == 1){
- if(states[states.size()-2] != 1) states.back() = 0;
+
+ if(*i > chain_start_idx && env_set[*i - 1]) {
+ stuff_to_check.insert(*i - 1);
}
- if(states.back() == 2){
- if(states[states.size()-2] != 2) states.back() = 0;
+
+ if(*i < chain_last_idx && env_set[*i + 1]) {
+ stuff_to_check.insert(*i + 1);
}
}
- //let's do the remaining residues in between
- for(uint i = 1; i < states.size()-1; ++i){
- if(states[i] == 0) continue;
- if(states[i] == 1){
- if(states[i-1] == 1 || states[i+1] == 1) continue;
- }
- if(states[i] == 2){
- if(states[i-1] == 2 || states[i+1] == 2) continue;
- }
- states[i] = 0;
- }
-
- //let's finally map the states to the hbond_env
- actual_idx = idx_handler->GetChainStartIdx(chain_idx);
- for(uint i = 0; i < idx_handler->GetChainSize(chain_idx); ++i){
- env_data_[actual_idx].state = states[i];
- ++actual_idx;
- }
-}
-void HBondEnvListener::UpdateStatesAndHPos(const std::vector<uint>& idx,
- const BackboneScoreEnv& env) {
+ for(std::set<uint>::iterator i = stuff_to_check.begin();
+ i != stuff_to_check.end(); ++i) {
+
+ Real phi = -1.0472;
+ Real psi = -0.78540;
- const loop::IdxHandler* idx_handler = env.GetIdxHandlerData();
+ uint chain_idx = idx_handler->GetChainIdx(*i);
+ uint chain_start_idx = idx_handler->GetChainIdx(chain_idx);
+ uint chain_last_idx = idx_handler->GetChainLastIdx(chain_idx);
+
+ if(*i > chain_start_idx && *i < chain_last_idx &&
+ env_set[*i - 1] && env_set[*i + 1]) {
- //we have to figure out to which chains the residue belong to
- std::set<uint> affected_chains;
- uint chain_start, chain_end;
- for(std::vector<uint>::const_iterator i = idx.begin(); i != idx.end(); ++i){
- for(uint j = 0; j < idx_handler->GetNumChains(); ++j){
- chain_start = idx_handler->GetChainStartIdx(j);
- chain_end = chain_start + idx_handler->GetChainSize(j);
- if( (*i >= chain_start) && (*i < chain_end) ){
- affected_chains.insert(j);
- break;
- }
+ phi = geom::DihedralAngle(env_data_[*i-1].c_pos, env_data_[*i].n_pos,
+ env_data_[*i].pos, env_data_[*i].c_pos);
+ psi = geom::DihedralAngle(env_data_[*i].n_pos, env_data_[*i].pos,
+ env_data_[*i].c_pos, env_data_[*i+1].n_pos);
+
+ // we only assign a state if both, phi and psi, are set
+ env_data_[*i].state = GetHBondState(phi,psi);
+
}
- }
- //call according function for every affected chain
- for(std::set<uint>::iterator i = affected_chains.begin();
- i != affected_chains.end(); ++i){
- this->UpdateStatesAndHPos(*i, env);
+ core::ConstructAtomPos(env_data_[*i].c_pos, env_data_[*i].pos,
+ env_data_[*i].n_pos, 0.9970, 2.0420,
+ phi + Real(M_PI), env_data_[*i].h_pos);
}
}
diff --git a/scoring/src/hbond_env_listener.hh b/scoring/src/hbond_env_listener.hh
index efcc9c5ff118332dd2599b18d699504158c00a78..95ff88717153f8a442350dab749b9a0af346cb3c 100644
--- a/scoring/src/hbond_env_listener.hh
+++ b/scoring/src/hbond_env_listener.hh
@@ -42,7 +42,7 @@ struct HBondSpatialOrganizerItem {
geom::Vec3 h_pos; //donor H
geom::Vec3 c_pos; //acceptor C
geom::Vec3 o_pos; //acceptor O
- uint state;
+ uint8_t state;
uint idx;
bool is_proline;
};
@@ -82,11 +82,8 @@ public:
private:
- //update the hbond states and h positions of one particular chain
- void UpdateStatesAndHPos(uint chain_idx, const BackboneScoreEnv& env);
-
- //figures out which chains are affected and calls the function above
- void UpdateStatesAndHPos(const std::vector<uint>& idx, const BackboneScoreEnv& env);
+ void UpdateStatesAndHPos(const std::vector<uint>& idx,
+ const BackboneScoreEnv& base_env);
HBondSpatialOrganizer env_;
HBondSpatialOrganizerItem* env_data_;
diff --git a/scoring/src/hbond_score.cc b/scoring/src/hbond_score.cc
index d686098764dd32be6d5bd853a20e2805559c9e12..dec216111f2c7d79037e1c75f8801c50c46518db 100644
--- a/scoring/src/hbond_score.cc
+++ b/scoring/src/hbond_score.cc
@@ -275,16 +275,18 @@ void HBondScorer::SetEnergy(uint state, uint d_bin, uint alpha_bin,
energies_[this->EIdx(state, d_bin, alpha_bin, beta_bin, gamma_bin)] = e;
}
-Real HBondScorer::CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const {
+Real HBondScorer::CalculateScore(uint start_resnum, uint num_residues,
+ uint chain_idx) const {
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
"HBondScorer::CalculateScore", 2);
std::vector<Real> scores;
std::vector<Real> internal_scores;
+ std::vector<uint> idx_vector =
+ idx_handler_->ToIdxVector(start_resnum, num_residues, chain_idx);
- this->Score(bb_list, start_resnum, chain_idx, scores, internal_scores);
+ this->Score(idx_vector, scores, internal_scores);
Real score = 0.0;
Real internal_score = 0.0;
@@ -304,8 +306,8 @@ Real HBondScorer::CalculateScore(const loop::BackboneList& bb_list,
}
-void HBondScorer::CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
+void HBondScorer::CalculateScoreProfile(uint start_resnum, uint num_residues,
+ uint chain_idx,
std::vector<Real>& profile) const {
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
@@ -313,8 +315,10 @@ void HBondScorer::CalculateScoreProfile(const loop::BackboneList& bb_list,
std::vector<Real> scores;
std::vector<Real> internal_scores;
+ std::vector<uint> idx_vector =
+ idx_handler_->ToIdxVector(start_resnum, num_residues, chain_idx);
- this->Score(bb_list, start_resnum, chain_idx, scores, internal_scores);
+ this->Score(idx_vector, scores, internal_scores);
profile.resize(scores.size());
for(uint i = 0; i < scores.size(); ++i) {
@@ -322,9 +326,8 @@ void HBondScorer::CalculateScoreProfile(const loop::BackboneList& bb_list,
}
}
-void HBondScorer::Score(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
- std::vector<Real>& scores,
+void HBondScorer::Score(const std::vector<uint>& indices,
+ std::vector<Real>& external_scores,
std::vector<Real>& internal_scores) const {
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
@@ -336,185 +339,109 @@ void HBondScorer::Score(const loop::BackboneList& bb_list,
}
// setup
- uint bb_list_size, start_idx, end_idx;
- idx_handler_->SetupScoreCalculation(bb_list, start_resnum, chain_idx,
- bb_list_size, start_idx, end_idx);
+ uint n_indices = indices.size();
- scores.assign(bb_list_size, 0.0);
- internal_scores.assign(bb_list_size, 0.0);
- if(bb_list_size == 0) return;
+ external_scores.assign(n_indices, 0.0);
+ internal_scores.assign(n_indices, 0.0);
+ if(n_indices == 0) return;
+
+ // we need to store which interactions are internal and which are external
+ std::vector<bool> occupied(idx_handler_->GetNumResidues(), false);
+ for(std::vector<uint>::const_iterator i = indices.begin();
+ i != indices.end(); ++i) {
+ occupied[*i] = true;
+ }
HBondSpatialOrganizer::WithinResult within_result;
std::pair<HBondSpatialOrganizerItem*,Real>* a;
- //note, that the n-terminal backbone hydrogen will be constructed below
- std::vector<geom::Vec3> h_positions(bb_list_size);
+ for(uint i = 0; i < n_indices; ++i) {
+ int my_idx = indices[i];
- for(uint i = 1; i < bb_list_size; ++i) {
- Real phi = bb_list.GetPhiTorsion(i);
- core::ConstructAtomPos(bb_list.GetC(i), bb_list.GetCA(i),
- bb_list.GetN(i), 0.9970, 2.0420,
- phi + Real(M_PI), h_positions[i]);
- }
-
- std::vector<uint> bb_list_states(bb_list_size);
- Real n_phi = -1.0472;
- Real c_psi = -0.78540;
- if(start_resnum > 1) {
- if(env_set_[start_idx - 1]) {
- n_phi = geom::DihedralAngle(env_->GetEnvData(start_idx - 1).c_pos,
- bb_list.GetN(0),
- bb_list.GetCA(0),
- bb_list.GetC(0));
+ if(!env_set_[my_idx]) {
+ throw promod3::Error("Must set environment before scoring!");
}
- }
- if((start_resnum + bb_list_size - 1) < idx_handler_->GetChainSize(chain_idx)) {
- if(env_set_[end_idx + 1]) {
- c_psi = geom::DihedralAngle(bb_list.GetN(bb_list_size-1),
- bb_list.GetCA(bb_list_size-1),
- bb_list.GetC(bb_list_size-1),
- env_->GetEnvData(end_idx + 1).n_pos);
+ int my_chain_idx = idx_handler_->GetChainIdx(my_idx);
+ const HBondSpatialOrganizerItem& my_env_data = env_->GetEnvData(my_idx);
+ int my_state = my_env_data.state;
+
+ if(my_state != 0) {
+
+ int my_chain_start_idx = idx_handler_->GetChainIdx(my_chain_idx);
+ int my_chain_last_idx = idx_handler_->GetChainLastIdx(my_chain_idx);
+ bool consistent_state = false;
+
+ if(my_idx > my_chain_start_idx && env_set_[my_idx - 1] &&
+ env_->GetEnvData(my_idx - 1).state == my_state) {
+ consistent_state = true;
+ } else if(my_idx < my_chain_last_idx && env_set_[my_idx + 1] &&
+ env_->GetEnvData(my_idx + 1).state == my_state) {
+ consistent_state = true;
+ }
+
+ if(!consistent_state) {
+ my_state = 0;
+ }
}
- }
- //we still have to calculate the position of the n-terminal hydrogen
- core::ConstructAtomPos(bb_list.GetC(0), bb_list.GetCA(0),
- bb_list.GetN(0), 0.9970, 2.0420,
- n_phi + Real(M_PI), h_positions[0]);
-
- if(bb_list_size == 1) {
- bb_list_states[0] = GetHBondState(n_phi, c_psi);
- }
- else{
- bb_list_states[0] = GetHBondState(n_phi,bb_list.GetPsiTorsion(0));
- bb_list_states[bb_list_size-1] =
- GetHBondState(bb_list.GetPhiTorsion(bb_list_size-1),c_psi);
- }
-
- for(uint i = 1; i < bb_list_size-1; ++i) {
- bb_list_states[i] = GetHBondState(bb_list.GetPhiTorsion(i),
- bb_list.GetPsiTorsion(i));
- }
+ within_result = env_->FindWithin(my_env_data.pos, 9.0);
+ a = within_result.first;
- //handle the terminal states
- int prev_state = 0;
- int next_state = 0;
- if(start_resnum > 1) {
- if(env_set_[start_idx - 1]) {
- prev_state = env_->GetEnvData(start_idx - 1).state;
- }
- }
- if((start_resnum + bb_list_size - 1) < idx_handler_->GetChainSize(chain_idx)) {
- if(env_set_[end_idx + 1]) {
- next_state = env_->GetEnvData(end_idx + 1).state;
- }
- }
+ for(uint j = 0; j < within_result.second; ++j) {
- if(bb_list_states[0] == 1) {
- if(prev_state != 1 && bb_list_states[1] != 1) {
- bb_list_states[0] = 0;
- }
- }
- if(bb_list_states[0] == 2) {
- if(prev_state != 2 && bb_list_states[1] != 2) {
- bb_list_states[0] = 0;
- }
- }
- if(bb_list_states[bb_list_size-1] == 1) {
- if(next_state != 1 && bb_list_states[bb_list_size-2] != 1) {
- bb_list_states[bb_list_size-1] = 0;
- }
- }
- if(bb_list_states[bb_list_size-1] == 2) {
- if(next_state != 2 && bb_list_states[bb_list_size-2] != 2) {
- bb_list_states[bb_list_size-1] = 0;
- }
- }
+ int state = 0;
+ int other_idx = a[j].first->idx;
+ int other_state = a[j].first->state;
- for(uint i = 1; i < bb_list_size-1; ++i) {
- if(bb_list_states[i] == 1) {
- if(bb_list_states[i-1] == 1 || bb_list_states[i+1] == 1) continue;
- }
- if(bb_list_states[i] == 2) {
- if(bb_list_states[i-1] == 2 || bb_list_states[i+1] == 2) continue;
- }
- bb_list_states[i] = 0;
- }
+ if(other_state != 0 && other_state == my_state) {
+
+ int other_chain_idx = idx_handler_->GetChainIdx(other_idx);
+ int other_chain_start_idx = idx_handler_->GetChainIdx(other_chain_idx);
+ int other_chain_last_idx = idx_handler_->GetChainLastIdx(other_chain_idx);
+
+ bool consistent_state = false;
- //let's go over every bb_list residue
- int state;
-
- if(do_external_scores_) {
- for(uint i = 0; i < bb_list_size; ++i) {
- within_result = env_->FindWithin(bb_list.GetCA(i), 9.0);
- a = within_result.first;
- for(uint j = 0; j < within_result.second; ++j) {
-
- if(a[j].first->idx < start_idx || a[j].first->idx > end_idx) {
-
- if(bb_list_states[i] == a[j].first->state) state = bb_list_states[i];
- else state = 0;
-
- //case of being donor for current bb item
- if(bb_list.GetOLC(i) != 'P') {
- scores[i] += this->EvalHBondPotential(state,bb_list.GetN(i),
- h_positions[i],
- a[j].first->pos,
- a[j].first->c_pos,
- a[j].first->o_pos);
- }
-
- //case of being acceptor for current bb item
- if(!a[j].first->is_proline) {
- scores[i] += this->EvalHBondPotential(state, a[j].first->n_pos,
- a[j].first->h_pos,
- bb_list.GetCA(i),
- bb_list.GetC(i),
- bb_list.GetO(i));
- }
+ if(other_idx > other_chain_start_idx && env_set_[other_idx - 1] &&
+ env_->GetEnvData(other_idx - 1).state == other_state) {
+ consistent_state = true;
+ } else if(other_idx < other_chain_last_idx && env_set_[other_idx + 1] &&
+ env_->GetEnvData(other_idx + 1).state == other_state) {
+ consistent_state = true;
+ }
+
+ if(consistent_state) {
+ state = my_state;
}
}
- }
- }
- if(!do_internal_scores_) {
- return;
- }
+ Real e = 0.0;
+
+ //case of being donor for current bb item
+ if(!my_env_data.is_proline) {
+ e += this->EvalHBondPotential(state, my_env_data.n_pos,
+ my_env_data.h_pos,
+ a[j].first->pos,
+ a[j].first->c_pos,
+ a[j].first->o_pos);
+ }
+
+ //case of being acceptor for current bb item
+ if(!a[j].first->is_proline) {
+ e += this->EvalHBondPotential(state, a[j].first->n_pos,
+ a[j].first->h_pos,
+ my_env_data.pos,
+ my_env_data.c_pos,
+ my_env_data.o_pos);
+ }
- // we do not have the loop internal interactions yet...
- for (uint i = 0; i < bb_list_size; ++i) {
- for (uint j = i+3; j < bb_list_size; ++j) {
- // hbonds closer than 3 residues do not really make sense...
- if (geom::Length2(bb_list.GetCA(i) - bb_list.GetCA(j)) > 81) continue;
-
- if (bb_list_states[i] == bb_list_states[j]) state = bb_list_states[i];
- else state = 0;
-
- // case of i being donor
- if(bb_list.GetOLC(i) != 'P') {
- Real score = this->EvalHBondPotential(state,
- bb_list.GetN(i),
- h_positions[i],
- bb_list.GetCA(j),
- bb_list.GetC(j),
- bb_list.GetO(j));
- internal_scores[i] += score;
- internal_scores[j] += score;
+ if(occupied[other_idx] && do_internal_scores_) {
+ internal_scores[i] += e;
}
- // case of i being acceptor
- if(bb_list.GetOLC(j) != 'P') {
- Real score = this->EvalHBondPotential(state,
- bb_list.GetN(j),
- h_positions[j],
- bb_list.GetCA(i),
- bb_list.GetC(i),
- bb_list.GetO(i));
-
- internal_scores[i] += score;
- internal_scores[j] += score;
+ if(!occupied[other_idx] && do_external_scores_) {
+ external_scores[i] += e;
}
}
}
diff --git a/scoring/src/hbond_score.hh b/scoring/src/hbond_score.hh
index 1e4ea8d5b4ff2d1435a06e0e55c7aa42f52c0b37..276d5b9fe100cb871384d1b348e5d81e5eb36ab6 100644
--- a/scoring/src/hbond_score.hh
+++ b/scoring/src/hbond_score.hh
@@ -45,11 +45,11 @@ public:
void DoNormalize(bool do_it) { normalize_ = do_it; }
- virtual Real CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const;
+ virtual Real CalculateScore(uint start_resnum, uint num_residues,
+ uint chain_idx) const;
- virtual void CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
+ virtual void CalculateScoreProfile(uint start_resnum, uint num_residues,
+ uint chain_idx,
std::vector<Real>& profile) const;
void AttachEnvironment(BackboneScoreEnv& env);
@@ -87,11 +87,8 @@ private:
return energies_[this->EIdx(state, d_bin, alpha_bin, beta_bin, gamma_bin)];
}
-
- void Score(const loop::BackboneList& bb_list,
- uint start_resnum,
- uint chain_idx,
- std::vector<Real>& scores,
+ void Score(const std::vector<uint>& indices,
+ std::vector<Real>& external_scores,
std::vector<Real>& internal_scores) const;
// portable serialization (only stores data relevant for Load/Save!)
diff --git a/scoring/src/pairwise_score.cc b/scoring/src/pairwise_score.cc
index 8362efb65d94a8fdbcdda78ee2a31857a3a39d51..372741965c6dd10cd8b978736b4cb44df7cf3c6f 100644
--- a/scoring/src/pairwise_score.cc
+++ b/scoring/src/pairwise_score.cc
@@ -9,57 +9,60 @@ PairwiseScorer::PairwiseScorer(): do_internal_scores_(true),
PairwiseScorer::~PairwiseScorer() { }
-Real PairwiseScorer::CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const {
+Real PairwiseScorer::CalculateScore(uint start_resnum, uint num_residues,
+ uint chain_idx) const {
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
"PairwiseScorer::CalculateScore", 2);
- std::vector<Real> scores;
+ std::vector<Real> external_scores;
std::vector<Real> internal_scores;
+ std::vector<uint> idx_vector =
+ idx_handler_->ToIdxVector(start_resnum, num_residues, chain_idx);
- this->Score(bb_list, start_resnum, chain_idx, scores, internal_scores);
+ this->Score(idx_vector, external_scores, internal_scores);
- Real score = 0.0;
+ Real external_score = 0.0;
Real internal_score = 0.0;
- for(uint i = 0; i < scores.size(); ++i) {
- score += scores[i];
+ for(uint i = 0; i < external_scores.size(); ++i) {
+ external_score += external_scores[i];
internal_score += internal_scores[i];
}
// all interactions in internal_score have been observed twice!
- score += Real(0.5)*internal_score;
+ external_score += Real(0.5)*internal_score;
- if(normalize_ && !scores.empty()) {
- score /= scores.size();
+ if(normalize_ && !external_scores.empty()) {
+ external_score /= external_scores.size();
}
- return score;
+ return external_score;
}
-void PairwiseScorer::CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
+void PairwiseScorer::CalculateScoreProfile(uint start_resnum, uint num_residues,
+ uint chain_idx,
std::vector<Real>& profile) const {
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
"PairwiseScorer::CalculateScoreProfile", 2);
- std::vector<Real> scores;
+ std::vector<Real> external_scores;
std::vector<Real> internal_scores;
+ std::vector<uint> idx_vector =
+ idx_handler_->ToIdxVector(start_resnum, num_residues, chain_idx);
- this->Score(bb_list, start_resnum, chain_idx, scores, internal_scores);
+ this->Score(idx_vector, external_scores, internal_scores);
- profile.resize(scores.size());
- for(uint i = 0; i < scores.size(); ++i) {
- profile[i] = scores[i] + internal_scores[i];
+ profile.resize(external_scores.size());
+ for(uint i = 0; i < external_scores.size(); ++i) {
+ profile[i] = external_scores[i] + internal_scores[i];
}
}
-void PairwiseScorer::Score(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
- std::vector<Real>& scores,
+void PairwiseScorer::Score(const std::vector<uint>& indices,
+ std::vector<Real>& external_scores,
std::vector<Real>& internal_scores) const {
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
@@ -71,52 +74,54 @@ void PairwiseScorer::Score(const loop::BackboneList& bb_list,
throw promod3::Error(err);
}
- uint bb_list_size, start_idx, end_idx;
- idx_handler_->SetupScoreCalculation(bb_list, start_resnum, chain_idx,
- bb_list_size, start_idx, end_idx);
+ uint n_indices = indices.size();
+
+ external_scores.assign(n_indices, 0.0);
+ internal_scores.assign(n_indices, 0.0);
+
+ if(n_indices == 0) return;
+
+ // we need to store which interactions are internal and which are external
+ std::vector<bool> occupied(idx_handler_->GetNumResidues(), false);
+ for(std::vector<uint>::const_iterator i = indices.begin();
+ i != indices.end(); ++i) {
+ occupied[*i] = true;
+ }
- scores.assign(bb_list_size, 0.0);
- internal_scores.assign(bb_list_size, 0.0);
- if(bb_list_size == 0) return;
- Real d;
- uint other_idx;
- for(uint i = 0; i < bb_list_size; ++i) {
+ for(uint i = 0; i < n_indices; ++i) {
+ uint my_idx = indices[i];
+
+ if(!env_set_[my_idx]) {
+ throw promod3::Error("Must set environment before scoring!");
+ }
+
const std::vector<PairwiseFunctionData>& functions =
- pairwise_functions_[start_idx + i];
+ pairwise_functions_[my_idx];
if(!functions.empty()) {
for(uint j = 0; j < functions.size(); ++j) {
- other_idx = functions[j].other_idx;
- if(do_external_scores_ &&
- (other_idx < start_idx || other_idx > end_idx)) {
- //it's a constraint to the environment!
- //let's see wether the env residue is set...
- if(env_set_[other_idx]) {
- if(functions[j].ft == CB_PAIRWISE_FUNCTION) {
- d = geom::Distance(bb_list.GetCB(i), cb_pos_[other_idx]);
- }
- else{
- d = geom::Distance(bb_list.GetCA(i), ca_pos_[other_idx]);
- }
- scores[i] += functions[j].f->Score(d);
+
+ uint other_idx = functions[j].other_idx;
+
+ if(env_set_[other_idx]) {
+ Real val;
+
+ if(functions[j].ft == CB_PAIRWISE_FUNCTION) {
+ val = geom::Distance(cb_pos_[my_idx], cb_pos_[other_idx]);
}
- }
- else if(do_internal_scores_) {
- //it's an internal contact!
- if(other_idx < start_idx + i) {
- //make sure, every constraint gets only evaluated once...
- if(functions[j].ft == CB_PAIRWISE_FUNCTION) {
- d = geom::Distance(bb_list.GetCB(i),
- bb_list.GetCB(other_idx - start_idx));
- }
- else{
- d = geom::Distance(bb_list.GetCA(i),
- bb_list.GetCA(other_idx - start_idx));
- }
- Real score = functions[j].f->Score(d);
- internal_scores[i] += score;
- internal_scores[other_idx - start_idx] += score;
+ else{
+ val = geom::Distance(ca_pos_[my_idx], ca_pos_[other_idx]);
+ }
+
+ val = functions[j].f->Score(val);
+
+ if(occupied[other_idx] && do_internal_scores_) {
+ internal_scores[i] += val;
+ }
+
+ if(!occupied[other_idx] && do_external_scores_) {
+ external_scores[i] += val;
}
}
}
diff --git a/scoring/src/pairwise_score.hh b/scoring/src/pairwise_score.hh
index 28e7a52f9d05bda3a543e5498e1755ddd804004d..6931fba504d33e343990712e722e6b2a71fd9933 100644
--- a/scoring/src/pairwise_score.hh
+++ b/scoring/src/pairwise_score.hh
@@ -25,21 +25,19 @@ public:
void DoNormalize(bool do_it) { normalize_ = do_it; }
- virtual Real CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const;
+ virtual Real CalculateScore(uint start_resnum, uint num_residues,
+ uint chain_idx) const;
- virtual void CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
+ virtual void CalculateScoreProfile(uint start_resnum, uint num_residues,
+ uint chain_idx,
std::vector<Real>& profile) const;
void AttachEnvironment(BackboneScoreEnv& env);
private:
- void Score(const loop::BackboneList& bb_list,
- uint start_resnum,
- uint chain_idx,
- std::vector<Real>& scores,
+ void Score(const std::vector<uint>& indices,
+ std::vector<Real>& external_scores,
std::vector<Real>& internal_scores) const;
bool do_internal_scores_;
diff --git a/scoring/src/reduced_score.cc b/scoring/src/reduced_score.cc
index 72e0d4e4d8d5cfcf159a5a6efb4c70faae1a8881..4c07f1344bd788cf4b708284e4cd0db1d8e1e84d 100644
--- a/scoring/src/reduced_score.cc
+++ b/scoring/src/reduced_score.cc
@@ -248,8 +248,8 @@ void ReducedScorer::SetEnergy(ost::conop::AminoAcid a, ost::conop::AminoAcid b,
}
-Real ReducedScorer::CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const {
+Real ReducedScorer::CalculateScore(uint start_resnum, uint num_residues,
+ uint chain_idx) const {
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
"ReducedScorer::CalculateScore", 2);
@@ -257,7 +257,10 @@ Real ReducedScorer::CalculateScore(const loop::BackboneList& bb_list,
// get scores
std::vector<Real> scores;
std::vector<Real> internal_scores;
- this->Score(bb_list, start_resnum, chain_idx, scores, internal_scores);
+ std::vector<uint> idx_vector =
+ idx_handler_->ToIdxVector(start_resnum, num_residues, chain_idx);
+
+ this->Score(idx_vector, scores, internal_scores);
// sum them
Real score = 0.0;
Real internal_score = 0.0;
@@ -277,8 +280,8 @@ Real ReducedScorer::CalculateScore(const loop::BackboneList& bb_list,
}
-void ReducedScorer::CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
+void ReducedScorer::CalculateScoreProfile(uint start_resnum, uint num_residues,
+ uint chain_idx,
std::vector<Real>& profile) const {
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
@@ -286,8 +289,10 @@ void ReducedScorer::CalculateScoreProfile(const loop::BackboneList& bb_list,
std::vector<Real> scores;
std::vector<Real> internal_scores;
+ std::vector<uint> idx_vector =
+ idx_handler_->ToIdxVector(start_resnum, num_residues, chain_idx);
- this->Score(bb_list, start_resnum, chain_idx, scores, internal_scores);
+ this->Score(idx_vector, scores, internal_scores);
profile.resize(scores.size());
for(uint i = 0; i < scores.size(); ++i) {
@@ -296,10 +301,10 @@ void ReducedScorer::CalculateScoreProfile(const loop::BackboneList& bb_list,
}
-void ReducedScorer::Score(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
- std::vector<Real>& scores,
+void ReducedScorer::Score(const std::vector<uint>& indices,
+ std::vector<Real>& external_scores,
std::vector<Real>& internal_scores) const {
+
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
"ReducedScorer::Score", 2);
@@ -309,35 +314,29 @@ void ReducedScorer::Score(const loop::BackboneList& bb_list,
String err = "Can only calculate score when environment is attached!";
throw promod3::Error(err);
}
- uint bb_list_size, start_idx, end_idx;
- idx_handler_->SetupScoreCalculation(bb_list, start_resnum, chain_idx,
- bb_list_size, start_idx, end_idx);
- scores.assign(bb_list_size, 0.0);
- internal_scores.assign(bb_list_size, 0.0);
- if(bb_list_size == 0) return;
+ uint n_indices = indices.size();
+
+ external_scores.assign(n_indices, 0.0);
+ internal_scores.assign(n_indices, 0.0);
+
+ if(n_indices == 0) return;
+
+ // we need to store which interactions are internal and which are external
+ std::vector<bool> occupied(idx_handler_->GetNumResidues(), false);
+ for(std::vector<uint>::const_iterator i = indices.begin();
+ i != indices.end(); ++i) {
+ occupied[*i] = true;
+ }
// let's reduce the cutoff a bit -> hack to ensure bin always within range
- const Real squared_cutoff = cutoff_*cutoff_ - 0.001;
const Real cutoff = cutoff_ - 0.001;
// get inv sizes for binning
const Real inv_dist_bin_size = dist_bins_ / cutoff_;
const Real inv_angle_bin_size = angle_bins_ / Real(M_PI);
const Real inv_dihedral_bin_size = dihedral_bins_ / (2 * Real(M_PI));
- // extract bb_list data
- std::vector<geom::Vec3> bb_list_axis(bb_list_size);
- std::vector<geom::Vec3> bb_list_ca_pos(bb_list_size);
- for (uint i = 0; i < bb_list_size; ++i) {
- bb_list_ca_pos[i] = bb_list.GetCA(i);
- bb_list_axis[i] = geom::Normalize(geom::Normalize(bb_list_ca_pos[i]
- -bb_list.GetN(i)) +
- geom::Normalize(bb_list_ca_pos[i]
- -bb_list.GetC(i)));
- }
-
// some param required in score calculation
- Real dist;
Real alpha;
Real beta;
Real gamma;
@@ -345,84 +344,62 @@ void ReducedScorer::Score(const loop::BackboneList& bb_list,
uint alpha_bin;
uint beta_bin;
uint gamma_bin;
+ Real e;
geom::Vec3 ca_vec;
+ // let's go over every loop residue
+ ReducedSpatialOrganizer::WithinResult within_result;
+ std::pair<ReducedSpatialOrganizerItem*,Real>* a;
+ for (uint i = 0; i < n_indices; ++i) {
+ int my_idx = indices[i];
- if(do_external_scores_) {
- // let's precalculate the idx ranges we have to neglect
- int min_idx_chain = idx_handler_->GetChainStartIdx(chain_idx);
- int max_idx_chain = min_idx_chain +
- idx_handler_->GetChainSize(chain_idx) - 1;
- std::vector<int> start_neglect(bb_list_size, start_idx);
- std::vector<int> end_neglect(bb_list_size, end_idx);
- int neglect_range = int(seq_sep_) - 1;
-
- int max_i = std::min(neglect_range, static_cast<int>(bb_list_size));
- for (int i = 0; i < max_i; ++i) {
- start_neglect[i] =
- std::max(min_idx_chain, start_neglect[i]-(neglect_range-i));
- end_neglect[bb_list_size-1-i] =
- std::min(max_idx_chain, end_neglect[i]+(neglect_range-i));
+ if(!env_set_[my_idx]) {
+ throw promod3::Error("Must set environment before scoring!");
}
- // let's go over every loop residue
- ReducedSpatialOrganizer::WithinResult within_result;
- std::pair<ReducedSpatialOrganizerItem*,Real>* a;
- for (uint i = 0; i < bb_list_size; ++i) {
- within_result = env_->FindWithin(bb_list_ca_pos[i], cutoff);
- a = within_result.first;
- const ost::conop::AminoAcid myaa = bb_list.GetAA(i);
- for (uint j = 0; j < within_result.second; ++j) {
- const int my_idx = int(a[j].first->idx);
- if (my_idx < start_neglect[i] || my_idx > end_neglect[i]) {
- ca_vec = a[j].first->pos - bb_list_ca_pos[i];
- ca_vec *= (Real(1.0) / a[j].second);
- GetReducedAngles(bb_list_axis[i], ca_vec, a[j].first->axis,
- alpha, beta, gamma);
- dist_bin = a[j].second * inv_dist_bin_size;
- alpha_bin = alpha * inv_angle_bin_size;
- beta_bin = beta * inv_angle_bin_size;
- gamma_bin = gamma * inv_dihedral_bin_size;
- scores[i] += energies_[this->EIdx(myaa, a[j].first->aa, dist_bin,
- alpha_bin, beta_bin, gamma_bin)];
+ int my_chain_idx = idx_handler_->GetChainIdx(my_idx);
+ ost::conop::AminoAcid my_aa = env_->GetEnvData(my_idx).aa;
+ geom::Vec3 my_pos = env_->GetEnvData(my_idx).pos;
+ geom::Vec3 my_axis = env_->GetEnvData(my_idx).axis;
+
+ within_result = env_->FindWithin(my_pos, cutoff);
+ a = within_result.first;
+
+ for (uint j = 0; j < within_result.second; ++j) {
+ int other_idx = int(a[j].first->idx);
+ int other_chain_idx = idx_handler_->GetChainIdx(other_idx);
+
+ if (my_chain_idx != other_chain_idx ||
+ std::abs(my_idx - other_idx) >= seq_sep_) {
+ ca_vec = a[j].first->pos - my_pos;
+ ca_vec *= (Real(1.0) / a[j].second);
+ GetReducedAngles(my_axis, ca_vec, a[j].first->axis,
+ alpha, beta, gamma);
+ dist_bin = a[j].second * inv_dist_bin_size;
+ alpha_bin = alpha * inv_angle_bin_size;
+ beta_bin = beta * inv_angle_bin_size;
+ gamma_bin = gamma * inv_dihedral_bin_size;
+
+ e = energies_[this->EIdx(my_aa, a[j].first->aa, dist_bin,
+ alpha_bin, beta_bin, gamma_bin)];
+
+ if(occupied[other_idx] && do_internal_scores_) {
+ internal_scores[i] += e;
}
- }
- }
- }
- if(do_internal_scores_) {
- if (bb_list_size >= seq_sep_) {
- for (uint i = 0; i < bb_list_size-seq_sep_; ++i) {
- const ost::conop::AminoAcid myaa = bb_list.GetAA(i);
- for (uint j = i+seq_sep_; j < bb_list_size; ++j) {
- ca_vec = bb_list_ca_pos[j] - bb_list_ca_pos[i];
- dist = geom::Length2(ca_vec);
- if (dist < squared_cutoff) {
- dist = std::sqrt(dist);
- ca_vec *= (Real(1.0) / dist);
- GetReducedAngles(bb_list_axis[i], ca_vec, bb_list_axis[j],
- alpha, beta, gamma);
- dist_bin = dist * inv_dist_bin_size;
- alpha_bin = alpha * inv_angle_bin_size;
- beta_bin = beta * inv_angle_bin_size;
- gamma_bin = gamma * inv_dihedral_bin_size;
-
- internal_scores[i] += energies_[this->EIdx(myaa, bb_list.GetAA(j),
- dist_bin, alpha_bin,
- beta_bin, gamma_bin)];
- internal_scores[j] += energies_[this->EIdx(bb_list.GetAA(j), myaa,
- dist_bin, beta_bin,
- alpha_bin, gamma_bin)];
- }
+ if(!occupied[other_idx] && do_external_scores_) {
+ external_scores[i] += e;
}
}
}
}
+
}
void ReducedScorer::AttachEnvironment(BackboneScoreEnv& env) {
env_ = env.UniqueListener<ReducedEnvListener>("ReducedEnvListener");
idx_handler_ = env.GetIdxHandlerData();
+ env_set_ = env.GetEnvSetData();
attached_environment_ = true;
}
diff --git a/scoring/src/reduced_score.hh b/scoring/src/reduced_score.hh
index f6586f8ac965df2737684b9b59e7c21bbcb9f729..20b478df27ccf0bf6a85191e37d767c69f2998e1 100644
--- a/scoring/src/reduced_score.hh
+++ b/scoring/src/reduced_score.hh
@@ -39,11 +39,11 @@ public:
void DoNormalize(bool do_it) { normalize_ = do_it; }
- virtual Real CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const;
+ virtual Real CalculateScore(uint start_resnum, uint num_residues,
+ uint chain_idx) const;
- virtual void CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
+ virtual void CalculateScoreProfile(uint start_resnum, uint num_residues,
+ uint chain_idx,
std::vector<Real>& profile) const;
void AttachEnvironment(BackboneScoreEnv& env);
@@ -65,10 +65,8 @@ private:
* angle_bins_ * dihedral_bins_;
}
- void Score(const loop::BackboneList& bb_list,
- uint start_resnum,
- uint chain_idx,
- std::vector<Real>& scores,
+ void Score(const std::vector<uint>& indices,
+ std::vector<Real>& external_scores,
std::vector<Real>& internal_scores) const;
// portable serialization (only stores data relevant for Load/Save!)
@@ -101,6 +99,7 @@ private:
//AttachEnvironment
bool attached_environment_;
ReducedEnvListener* env_;
+ const int* env_set_;
const loop::IdxHandler* idx_handler_;
};
diff --git a/scoring/src/scoring_object_loader.cc b/scoring/src/scoring_object_loader.cc
index 043225ba2c95a271f79f632fd8ce1a24ee403543..6bbc5a96c4770ea424f0085c45da515324b8028f 100644
--- a/scoring/src/scoring_object_loader.cc
+++ b/scoring/src/scoring_object_loader.cc
@@ -3,7 +3,6 @@
#include <promod3/core/runtime_profiling.hh>
#include <promod3/scoring/clash_score.hh>
#include <promod3/scoring/pairwise_score.hh>
-#include <promod3/scoring/density_score.hh>
#include <promod3/scoring/all_atom_clash_score.hh>
namespace promod3{ namespace scoring {
@@ -84,7 +83,6 @@ BackboneOverallScorerPtr LoadDefaultBackboneOverallScorer() {
(*scorer)["ss_agreement"] = LoadSSAgreementScorer();
(*scorer)["torsion"] = LoadTorsionScorer();
(*scorer)["pairwise"] = PairwiseScorerPtr(new PairwiseScorer);
- (*scorer)["density"] = DensityScorerPtr(new DensityScorer);
return scorer;
}
diff --git a/scoring/src/ss_agreement_env_listener.cc b/scoring/src/ss_agreement_env_listener.cc
index 7456098a432ae9de23a7fe66030dd39daef0f12c..3cf2e41a5a516dc593b6aa220ca490b1d23a3cc8 100644
--- a/scoring/src/ss_agreement_env_listener.cc
+++ b/scoring/src/ss_agreement_env_listener.cc
@@ -48,11 +48,6 @@ void SSAgreementEnvListener::SetEnvironment(const BackboneScoreEnv& base_env,
core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
"SSAgreementEnvListener::SetEnvironment", 2);
- if(!temp_indices_.empty()) {
- throw promod3::Error("Cannot set SSAgreement environment when temporary"
- " data is set!");
- }
-
const geom::Vec3* n_pos_data = base_env.GetNPosData();
const geom::Vec3* ca_pos_data = base_env.GetCAPosData();
const geom::Vec3* c_pos_data = base_env.GetCPosData();
@@ -76,11 +71,6 @@ void SSAgreementEnvListener::ResetEnvironment(const BackboneScoreEnv& base_env,
core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
"SSAgreementEnvListener::ResetEnvironment", 2);
- if(!temp_indices_.empty()) {
- throw promod3::Error("Cannot reset SSAgreement environment when temporary"
- " data is set!");
- }
-
const geom::Vec3* n_pos_data = base_env.GetNPosData();
const geom::Vec3* ca_pos_data = base_env.GetCAPosData();
const geom::Vec3* c_pos_data = base_env.GetCPosData();
@@ -106,11 +96,6 @@ void SSAgreementEnvListener::ClearEnvironment(const BackboneScoreEnv& base_env,
core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
"SSAgreementEnvListener::ClearEnvironment", 2);
- if(!temp_indices_.empty()) {
- throw promod3::Error("Cannot clear SSAgreement environment when temporary"
- " data is set!");
- }
-
for(std::vector<uint>::const_iterator i = idx.begin();
i != idx.end(); ++i){
env_.Remove(&env_data_[*i], env_data_[*i].pos);
@@ -147,161 +132,6 @@ void SSAgreementEnvListener::ClearEnvironment(const BackboneScoreEnv& base_env,
}
}
-void SSAgreementEnvListener::SetTemporaryInteractions(
- const promod3::loop::BackboneList& bb_list,
- uint start_idx) {
-
- temp_start_idx_ = start_idx;
- temp_size_ = bb_list.size();
-
- // lets first copy over all the data from the bb_list, as well as one before
- // and one after
- int idx = std::max(0, temp_start_idx_ - 1);
- int end = std::min(static_cast<int>(idx_handler_->GetNumResidues()),
- temp_start_idx_ + temp_size_ + 1);
-
- for(; idx < end; ++idx) {
- temp_indices_.push_back(idx);
- original_env_data_.push_back(env_data_[idx]);
- original_ca_pos_.push_back(ca_positions_[idx]);
- original_donor_for_one_.push_back(donor_for_one_[idx]);
- original_donor_for_two_.push_back(donor_for_two_[idx]);
- original_connected_to_next_.push_back(connected_to_next_[idx]);
- }
-
- // update the env_data_ vector and the spatial organizer
- for(int i = 0; i < temp_size_; ++i) {
-
- // if the env is set, we have to reset the spatial organizer
- // we just add it otherwise
- if(env_set_[temp_start_idx_+i]) {
- env_.Reset(&env_data_[temp_start_idx_+i],
- env_data_[temp_start_idx_+i].pos, bb_list.GetCA(i));
- } else {
- env_.Add(&env_data_[temp_start_idx_+i], bb_list.GetCA(i));
- }
-
- SSAgreementSpatialOrganizerItem& organizer_item = env_data_[temp_start_idx_+i];
- organizer_item.pos = bb_list.GetCA(i);
- organizer_item.n_pos = bb_list.GetN(i);
- organizer_item.c_pos = bb_list.GetC(i);
- organizer_item.o_pos = bb_list.GetO(i);
-
- // and directly copy over the ca position
- ca_positions_[temp_start_idx_+i] = bb_list.GetCA(i);
- }
-
- // lets update connectivity
-
- // first the residue before
- if(temp_start_idx_ > 0 && env_set_[temp_start_idx_-1] &&
- geom::Length2(env_data_[temp_start_idx_-1].c_pos -
- env_data_[temp_start_idx_].n_pos) <= 6.25) {
- connected_to_next_[temp_start_idx_-1] = 1;
- } else {
- connected_to_next_[temp_start_idx_-1] = 0;
- }
-
- // connection of bb_list itself
- for(int i = 0; i < temp_size_ - 1; ++i) {
- if(geom::Length2(bb_list.GetC(i) - bb_list.GetN(i+1)) <= 6.25) {
- connected_to_next_[temp_start_idx_+i] = 1;
- } else {
- connected_to_next_[temp_start_idx_+i] = 0;
- }
- }
-
- // lets see whether the last residue is connected
- if(temp_start_idx_ + temp_size_ <
- static_cast<int>(idx_handler_->GetNumResidues()) &&
- env_set_[temp_start_idx_ + temp_size_] &&
- geom::Length2(env_data_[temp_start_idx_+temp_size_-1].c_pos -
- env_data_[temp_start_idx_+temp_size_].n_pos) <= 6.25) {
- connected_to_next_[temp_start_idx_+temp_size_-1] = 1;
- } else {
- connected_to_next_[temp_start_idx_+temp_size_-1] = 0;
- }
-
- // set h positions of all loop residues, as well as the one
- // residue after (if present)
- idx = temp_start_idx_;
- end = temp_start_idx_ + temp_size_;
- if(static_cast<uint>(end) < idx_handler_->GetNumResidues()) ++end;
- for(; idx < end; ++idx) {
- // h_pos of this residue
- if(!env_data_[idx].is_proline) {
- if(idx > 0 && connected_to_next_[idx-1]) {
- geom::Vec3 dir_vec = geom::Normalize(env_data_[idx-1].c_pos -
- env_data_[idx-1].o_pos);
- env_data_[idx].h_pos = env_data_[idx].n_pos + dir_vec;
- env_data_[idx].valid_h_pos = true;
- }
- else {
- env_data_[idx].valid_h_pos = false;
- }
- }
- }
-
- // find the optimal donor for all bb_list residues
- for(int i = 0; i < temp_size_; ++i) {
- this->FindOptimalDonors(temp_start_idx_ + i, false);
- }
-
- // we have now everything for the bb_list itself! but we also need
- // to update the donor information of close residue and store
- // their original values, so they can be restored!
-
- std::set<int> close_residues;
- for(int i = 0; i < temp_size_; ++i) {
- this->AddCloseResidues(temp_start_idx_ + i, close_residues);
- }
-
- for(std::set<int>::iterator i = close_residues.begin();
- i != close_residues.end(); ++i) {
- // the residues of the bb_list have already been handled
- int idx = *i;
- if(idx >= temp_start_idx_ && idx < temp_start_idx_ + temp_size_) continue;
-
- temp_indices_.push_back(idx);
- original_env_data_.push_back(env_data_[idx]);
- original_ca_pos_.push_back(ca_positions_[idx]);
- original_donor_for_one_.push_back(donor_for_one_[idx]);
- original_donor_for_two_.push_back(donor_for_two_[idx]);
- original_connected_to_next_.push_back(connected_to_next_[idx]);
-
- this->FindOptimalDonors(idx);
- }
-}
-
-void SSAgreementEnvListener::Restore() {
-
- for(uint i = 0; i < temp_indices_.size(); ++i) {
- int temp_idx = temp_indices_[i];
- if(temp_idx >= temp_start_idx_ && temp_idx < temp_start_idx_ + temp_size_) {
- // its in the temporarily set bb_list, so we have to handle the
- // spatial organizer
- if(env_set_[temp_idx]) {
- env_.Reset(&env_data_[temp_idx],
- env_data_[temp_idx].pos, original_ca_pos_[i]);
- } else {
- env_.Remove(&env_data_[temp_idx], env_data_[temp_idx].pos);
- }
- }
- env_data_[temp_idx] = original_env_data_[i];
- ca_positions_[temp_idx] = original_ca_pos_[i];
- donor_for_one_[temp_idx] = original_donor_for_one_[i];
- donor_for_two_[temp_idx] = original_donor_for_two_[i];
- connected_to_next_[temp_idx] = original_connected_to_next_[i];
- }
-
- temp_indices_.clear();
- original_env_data_.clear();
- original_ca_pos_.clear();
- original_donor_for_one_.clear();
- original_donor_for_two_.clear();
- original_connected_to_next_.clear();
-}
-
void SSAgreementEnvListener::AddCloseResidues(uint idx,
std::set<int>& close_residues) const {
SSAgreementSpatialOrganizer::WithinResult within_result =
@@ -312,13 +142,12 @@ void SSAgreementEnvListener::AddCloseResidues(uint idx,
}
}
-void SSAgreementEnvListener::FindOptimalDonors(uint idx, bool check_env_set) {
+void SSAgreementEnvListener::FindOptimalDonors(uint idx) {
donor_for_one_[idx] = -1;
donor_for_two_[idx] = -1;
if(!env_data_[idx].valid_h_pos) return;
- if(check_env_set && !env_set_[idx]) return;
std::pair<Real, Real> donor_energies =
std::make_pair(std::numeric_limits<Real>::max(),
@@ -330,7 +159,7 @@ void SSAgreementEnvListener::FindOptimalDonors(uint idx, bool check_env_set) {
std::pair<SSAgreementSpatialOrganizerItem*,Real>* a = within_result.first;
for (uint j = 0; j < within_result.second; ++j) {
const uint other_idx = a[j].first->idx;
- if((other_idx < idx) || (other_idx > idx + 1)) {
+ if((idx > 0 && (other_idx < idx - 1)) || (other_idx > (idx + 1))) {
Real e = ost::mol::alg::DSSPHBondEnergy(env_data_[idx].h_pos,
env_data_[idx].n_pos,
env_data_[other_idx].c_pos,
@@ -363,9 +192,10 @@ void SSAgreementEnvListener::Update(const std::vector<uint>& idx) {
for(std::vector<uint>::const_iterator i = idx.begin(); i != idx.end(); ++i){
ca_positions_[*i] = env_data_[*i].pos;
+ bool last_residue = (*i >= idx_handler_->GetNumResidues() - 1);
// set connectivity of the residue itself
- if(*i < idx_handler_->GetNumResidues() - 1 && env_set_[*i + 1] &&
+ if(!last_residue && env_set_[*i + 1] &&
geom::Length2(env_data_[*i].c_pos - env_data_[*i+1].n_pos) <= 6.25) {
connected_to_next_[*i] = 1;
} else {
@@ -396,8 +226,6 @@ void SSAgreementEnvListener::Update(const std::vector<uint>& idx) {
}
// the h_pos of the next residue might have changed
- // no check for index validity of idx + 1... this informatin is contained
- // in connected_to_next_
if(connected_to_next_[*i]) {
if(!env_data_[*i+1].is_proline) {
geom::Vec3 dir_vec = geom::Normalize(env_data_[*i].c_pos -
@@ -407,6 +235,8 @@ void SSAgreementEnvListener::Update(const std::vector<uint>& idx) {
} else {
env_data_[*i+1].valid_h_pos = false;
}
+ } else if(!last_residue) {
+ env_data_[*i+1].valid_h_pos = false;
}
}
diff --git a/scoring/src/ss_agreement_env_listener.hh b/scoring/src/ss_agreement_env_listener.hh
index 3d924e23419591bf88567186eaa6f1949cee9687..ca537387da4e38f99fef186188ef975384636936 100644
--- a/scoring/src/ss_agreement_env_listener.hh
+++ b/scoring/src/ss_agreement_env_listener.hh
@@ -51,11 +51,6 @@ public:
virtual String WhoAmI() const { return "SSAgreementEnvListener"; }
- void SetTemporaryInteractions(const promod3::loop::BackboneList& bb_list,
- uint start_idx);
-
- void Restore();
-
const std::vector<geom::Vec3>& GetCAPositions() const { return ca_positions_; }
const std::vector<int>& GetDonorForOne() const { return donor_for_one_; }
const std::vector<int>& GetDonorForTwo() const { return donor_for_two_; }
@@ -70,31 +65,18 @@ private:
void AddCloseResidues(uint idx, std::set<int>& close_residues) const;
- void FindOptimalDonors(uint idx, bool check_env_set = true);
+ void FindOptimalDonors(uint idx);
SSAgreementSpatialOrganizer env_;
SSAgreementSpatialOrganizerItem* env_data_;
const int* env_set_;
const loop::IdxHandler* idx_handler_;
+ // thats the input for the ost dssp implementation
std::vector<geom::Vec3> ca_positions_;
std::vector<int> donor_for_one_;
std::vector<int> donor_for_two_;
std::vector<int> connected_to_next_;
-
- // when temporary stuff gets set, the original values get
- // stored in the vectors below. When calling Restore(),
- // the stuff goes back into the original vectors
-
- // the derived data
- int temp_start_idx_;
- int temp_size_;
- std::vector<int> temp_indices_;
- std::vector<SSAgreementSpatialOrganizerItem> original_env_data_;
- std::vector<geom::Vec3> original_ca_pos_;
- std::vector<int> original_donor_for_one_;
- std::vector<int> original_donor_for_two_;
- std::vector<int> original_connected_to_next_;
};
}} // ns
diff --git a/scoring/src/ss_agreement_score.cc b/scoring/src/ss_agreement_score.cc
index a90d155978412b5f8402ea302767c8ad2b1251ed..568eec17677accf9669e3dc1b69145275d639a16 100644
--- a/scoring/src/ss_agreement_score.cc
+++ b/scoring/src/ss_agreement_score.cc
@@ -160,43 +160,69 @@ void SSAgreementScorer::SetScore(char psipred_state, int psipred_cfi,
}
-Real SSAgreementScorer::CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const{
+Real SSAgreementScorer::CalculateScore(uint start_resnum, uint num_residues,
+ uint chain_idx) const{
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
"SSAgreementScorer::CalculateScore", 2);
- std::vector<Real> scores;
+ std::vector<uint> v_start_indices(1, idx_handler_->ToIdx(start_resnum,
+ chain_idx));
+ std::vector<uint> v_num_residues(1, num_residues);
- this->Score(bb_list, start_resnum, chain_idx, scores);
+ if(v_start_indices[0] + v_num_residues[0] - 1 >
+ idx_handler_->GetChainLastIdx(chain_idx)) {
+ throw promod3::Error("Invalid start_resnum and loop length!");
+ }
+
+ std::vector<std::vector<Real> > scores;
+
+ this->Score(v_start_indices, v_num_residues, scores);
Real score = 0.0;
+ uint n = 0;
for(uint i = 0; i < scores.size(); ++i){
- score += scores[i];
+ for(uint j = 0; j < scores[i].size(); ++j) {
+ score += scores[i][j];
+ }
+ n += scores[i].size();
}
- if(normalize_ && !scores.empty()) {
- score /= scores.size();
+ if(normalize_ && n > 0) {
+ score /= n;
}
return score;
}
-void SSAgreementScorer::CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
+void SSAgreementScorer::CalculateScoreProfile(uint start_resnum, uint num_residues,
+ uint chain_idx,
std::vector<Real>& profile) const{
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
"SSAgreementScorer::CalculateScoreProfile", 2);
- this->Score(bb_list, start_resnum, chain_idx, profile);
+ std::vector<uint> v_start_indices(1, idx_handler_->ToIdx(start_resnum,
+ chain_idx));
+ std::vector<uint> v_num_residues(1, num_residues);
+
+ if(v_start_indices[0] + v_num_residues[0] - 1 >
+ idx_handler_->GetChainLastIdx(chain_idx)) {
+ throw promod3::Error("Invalid start_resnum and loop length!");
+ }
+
+ std::vector<std::vector<Real> > scores;
+
+ this->Score(v_start_indices, v_num_residues, scores);
+
+ profile = scores[0];
}
-void SSAgreementScorer::Score(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
- std::vector<Real>& scores) const{
+void SSAgreementScorer::Score(const std::vector<uint>& start_indices,
+ const std::vector<uint>& num_residues,
+ std::vector<std::vector<Real> >& scores) const{
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
"SSAgreementScorer::Score", 2);
@@ -206,27 +232,39 @@ void SSAgreementScorer::Score(const loop::BackboneList& bb_list,
throw promod3::Error(err);
}
- // setup
- uint bb_list_size, start_idx, end_idx;
- idx_handler_->SetupScoreCalculation(bb_list, start_resnum, chain_idx,
- bb_list_size, start_idx, end_idx);
+ uint n_profiles = start_indices.size();
- scores.assign(bb_list_size, 0.0);
- if(bb_list_size == 0) return;
+ if(n_profiles == 0) return;
- env_->SetTemporaryInteractions(bb_list, start_idx);
+ scores.assign(n_profiles, std::vector<Real>());
- String ss = ost::mol::alg::RawEstimateSS(env_->GetCAPositions(),
- start_idx, bb_list_size,
- env_->GetDonorForOne(),
- env_->GetDonorForTwo(),
- env_->GetConnection());
- env_->Restore();
+ // check whether the full env is set
+ for(uint i = 0; i < start_indices.size(); ++i) {
+ for(uint j = 0; j < num_residues[i]; ++j) {
+ if(!env_set_[start_indices[i] + j]) {
+ throw promod3::Error("Must set environment before scoring!");
+ }
+ }
+ scores[i].assign(num_residues[i], 0.0);
+ }
- for(uint i = 0; i < bb_list_size; ++i){
- scores[i] = scores_[this->SIdx(GetPsipredIdx(psipred_pred_[start_idx + i]),
- psipred_cfi_[start_idx + i],
- GetDSSPIdx(ss[i]))];
+
+ for(uint i = 0; i < n_profiles; ++i) {
+
+ if(num_residues[i] == 0) continue;
+
+ String ss = ost::mol::alg::RawEstimateSS(env_->GetCAPositions(),
+ start_indices[i],
+ num_residues[i],
+ env_->GetDonorForOne(),
+ env_->GetDonorForTwo(),
+ env_->GetConnection());
+
+ for(uint j = 0; j < num_residues[i]; ++j){
+ scores[i][j] = scores_[this->SIdx(GetPsipredIdx(psipred_pred_[start_indices[i] + j]),
+ psipred_cfi_[start_indices[i] + j],
+ GetDSSPIdx(ss[j]))];
+ }
}
}
diff --git a/scoring/src/ss_agreement_score.hh b/scoring/src/ss_agreement_score.hh
index f696e1d2e77f99c04750bc8328c36a62de34614e..d77e80e2fd33aceecd9d274089a6635ddab04795 100644
--- a/scoring/src/ss_agreement_score.hh
+++ b/scoring/src/ss_agreement_score.hh
@@ -59,11 +59,10 @@ public:
void SetScore(char psipred_state, int psipred_cfi,
char dssp_state, Real score);
- virtual Real CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const;
+ virtual Real CalculateScore(uint start_resnum, uint num_residues,
+ uint chain_idx) const;
- virtual void CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum,
+ virtual void CalculateScoreProfile(uint start_resnum, uint num_residues,
uint chain_idx,
std::vector<Real>& profile) const;
@@ -80,10 +79,9 @@ private:
dssp_idx;
}
- void Score(const loop::BackboneList& bb_list,
- uint start_resnum,
- uint chain_idx,
- std::vector<Real>& scores) const;
+ void Score(const std::vector<uint>& start_indices,
+ const std::vector<uint>& num_residues,
+ std::vector<std::vector<Real> >& scores) const;
// portable serialization (only stores data relevant for Load/Save!)
// (cleanly element by element with fixed-width base-types)
diff --git a/scoring/src/torsion_score.cc b/scoring/src/torsion_score.cc
index 2dfb83819261499d3d6cbeaf056141a4e3e70cf6..0b6cfec7c5352a866d26a6c8d7484d886b1e29bc 100644
--- a/scoring/src/torsion_score.cc
+++ b/scoring/src/torsion_score.cc
@@ -238,14 +238,16 @@ void TorsionScorer::SetEnergy(uint group_idx, uint phi_bin, uint psi_bin, Real e
}
-Real TorsionScorer::CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const {
+Real TorsionScorer::CalculateScore(uint start_resnum, uint num_residues,
+ uint chain_idx) const {
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
"TorsionScorer::CalculateScore", 2);
std::vector<Real> scores;
- this->Score(bb_list, start_resnum, chain_idx, scores);
+ std::vector<uint> idx_vector =
+ idx_handler_->ToIdxVector(start_resnum, num_residues, chain_idx);
+ this->Score(idx_vector, scores);
Real score = 0.0;
for(uint i = 0; i < scores.size(); ++i) {
@@ -260,19 +262,20 @@ Real TorsionScorer::CalculateScore(const loop::BackboneList& bb_list,
}
-void TorsionScorer::CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
+void TorsionScorer::CalculateScoreProfile(uint start_resnum, uint num_residues,
+ uint chain_idx,
std::vector<Real>& profile) const {
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
"TorsionScorer::CalculateScoreProfile", 2);
- this->Score(bb_list, start_resnum, chain_idx, profile);
+ std::vector<uint> idx_vector =
+ idx_handler_->ToIdxVector(start_resnum, num_residues, chain_idx);
+ this->Score(idx_vector, profile);
}
-void TorsionScorer::Score(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
+void TorsionScorer::Score(const std::vector<uint>& indices,
std::vector<Real>& scores) const {
promod3::core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
@@ -284,107 +287,46 @@ void TorsionScorer::Score(const loop::BackboneList& bb_list,
}
// setup
- uint bb_list_size, start_idx, end_idx;
- idx_handler_->SetupScoreCalculation(bb_list, start_resnum, chain_idx,
- bb_list_size, start_idx, end_idx);
+ uint n_indices = indices.size();
- scores.assign(bb_list_size, 0.0);
- if(bb_list_size == 0) return;
+ scores.assign(n_indices, 0.0);
+
+ if(n_indices == 0) return;
Real angle1, angle2;
uint bin1, bin2;
- uint torsion_group_idx;
-
- for(uint i = 1; i < bb_list_size - 1; ++i) {
- torsion_group_idx = torsion_group_indices_[start_idx + i];
- angle1 = bb_list.GetPhiTorsion(i);
- angle2 = bb_list.GetPsiTorsion(i);
- bin1 = std::min(static_cast<uint>((angle1+Real(M_PI))/bin_size_),
- bins_per_dimension_-1);
- bin2 = std::min(static_cast<uint>((angle2+Real(M_PI))/bin_size_),
- bins_per_dimension_-1);
-
- scores[i] =
- energies_[this->EIdx(torsion_group_idx, bin1 * bins_per_dimension_ + bin2)];
- }
- // first handle the special case of bb_list being of size one
- if(bb_list_size == 1) {
+ for(uint i = 0; i < n_indices; ++i) {
+ int my_idx = indices[i];
- if(start_resnum <= 1 || start_resnum >= idx_handler_->GetChainSize(chain_idx)) {
- //there's nothing we can do...
- return;
+ if(!env_set_[my_idx]) {
+ throw promod3::Error("Must set environment before scoring!");
}
- torsion_group_idx = torsion_group_indices_[start_idx];
-
- angle1 = geom::DihedralAngle(c_pos_[start_idx - 1],
- bb_list.GetN(0),
- bb_list.GetCA(0),
- bb_list.GetC(0));
-
- angle2 = geom::DihedralAngle(bb_list.GetN(0),
- bb_list.GetCA(0),
- bb_list.GetC(0),
- n_pos_[end_idx + 1]);
-
- bin1 = std::min(static_cast<uint>((angle1+Real(M_PI))/bin_size_),
- bins_per_dimension_-1);
- bin2 = std::min(static_cast<uint>((angle2+Real(M_PI))/bin_size_),
- bins_per_dimension_-1);
-
- scores[0] =
- energies_[this->EIdx(torsion_group_idx, bin1 * bins_per_dimension_ + bin2)];
-
- // that's it... since we just dealt with the only residue, the stuff below
- // is not necessary anymore
- return;
- }
+ int my_chain_idx = idx_handler_->GetChainIdx(my_idx);
+ int my_chain_first_idx = idx_handler_->GetChainStartIdx(my_chain_idx);
+ int my_chain_last_idx = idx_handler_->GetChainLastIdx(my_chain_idx);
+ if(my_idx > my_chain_first_idx && env_set_[my_idx - 1] &&
+ my_idx < my_chain_last_idx && env_set_[my_idx + 1]) {
- if(start_resnum > 1) {
- //check, whether the environment for the residue before the nstem is set
- if(env_set_[start_idx - 1]) {
- torsion_group_idx = torsion_group_indices_[start_idx];
-
- //the data for the residues before and after get extracted from the env
- angle1 = geom::DihedralAngle(c_pos_[start_idx - 1],
- bb_list.GetN(0),
- bb_list.GetCA(0),
- bb_list.GetC(0));
-
- angle2 = bb_list.GetPsiTorsion(0);
-
- bin1 = std::min(static_cast<uint>((angle1+Real(M_PI))/bin_size_),
- bins_per_dimension_-1);
- bin2 = std::min(static_cast<uint>((angle2+Real(M_PI))/bin_size_),
- bins_per_dimension_-1);
+ angle1 = geom::DihedralAngle(c_pos_[my_idx - 1],
+ n_pos_[my_idx],
+ ca_pos_[my_idx],
+ c_pos_[my_idx]);
- scores[0] =
- energies_[this->EIdx(torsion_group_idx, bin1 * bins_per_dimension_ +
- bin2)];
- }
- }
+ angle2 = geom::DihedralAngle(n_pos_[my_idx],
+ ca_pos_[my_idx],
+ c_pos_[my_idx],
+ n_pos_[my_idx + 1]);
- if((start_resnum + bb_list_size - 1) < idx_handler_->GetChainSize(chain_idx)) {
- //check, whether the environment for the residue after the cstem is set
- if(env_set_[end_idx + 1]) {
- torsion_group_idx = torsion_group_indices_[end_idx];
-
- //the data for the residues before and after get extracted from the env
- angle1 = bb_list.GetPhiTorsion(bb_list_size-1);
- angle2 = geom::DihedralAngle(bb_list.GetN(bb_list_size-1),
- bb_list.GetCA(bb_list_size-1),
- bb_list.GetC(bb_list_size-1),
- n_pos_[end_idx + 1]);
bin1 = std::min(static_cast<uint>((angle1+Real(M_PI))/bin_size_),
bins_per_dimension_-1);
bin2 = std::min(static_cast<uint>((angle2+Real(M_PI))/bin_size_),
bins_per_dimension_-1);
- scores[bb_list_size - 1] =
- energies_[this->EIdx(torsion_group_idx, bin1 * bins_per_dimension_ +
- bin2)];
+ scores[i] = energies_[this->EIdx(torsion_group_indices_[my_idx],
+ bin1 * bins_per_dimension_ + bin2)];
}
}
}
@@ -394,6 +336,7 @@ void TorsionScorer::AttachEnvironment(BackboneScoreEnv& env) {
idx_handler_ = env.GetIdxHandlerData();
env_set_ = env.GetEnvSetData();
n_pos_ = env.GetNPosData();
+ ca_pos_ = env.GetCAPosData();
c_pos_ = env.GetCPosData();
torsion_group_indices_.assign(idx_handler_->GetNumResidues(),0);
diff --git a/scoring/src/torsion_score.hh b/scoring/src/torsion_score.hh
index f7ea8efce672d7792077db9d11b2fd08995acf3f..631dc2637f90f13ad3a03ed5e1c76f1488f01d7b 100644
--- a/scoring/src/torsion_score.hh
+++ b/scoring/src/torsion_score.hh
@@ -33,11 +33,11 @@ public:
void SetEnergy(uint group_idx, uint phi_bin, uint psi_bin, Real e);
- virtual Real CalculateScore(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx) const;
+ virtual Real CalculateScore(uint start_resnum, uint num_residues,
+ uint chain_idx) const;
- virtual void CalculateScoreProfile(const loop::BackboneList& bb_list,
- uint start_resnum, uint chain_idx,
+ virtual void CalculateScoreProfile(uint start_resnum, uint num_residues,
+ uint chain_idx,
std::vector<Real>& profile) const;
void DoNormalize(bool do_it) { normalize_ = do_it; }
@@ -52,9 +52,7 @@ private:
return group_idx * bins_squared_ + torsion_bin;
}
- void Score(const loop::BackboneList& bb_list,
- uint start_resnum,
- uint chain_idx,
+ void Score(const std::vector<uint>& indices,
std::vector<Real>& scores) const;
// portable serialization (only stores data relevant for Load/Save!)
@@ -92,6 +90,7 @@ private:
const loop::IdxHandler* idx_handler_;
const int* env_set_;
const geom::Vec3* n_pos_;
+ const geom::Vec3* ca_pos_;
const geom::Vec3* c_pos_;
std::vector<uint> torsion_group_indices_;
};
diff --git a/scoring/tests/test_backbone_scorer.cc b/scoring/tests/test_backbone_scorer.cc
index d1080b45c831ec12ffe2c25f431628ad8213af50..db23d2a95364cd1f91ac6d2ebe7b203751543854 100644
--- a/scoring/tests/test_backbone_scorer.cc
+++ b/scoring/tests/test_backbone_scorer.cc
@@ -120,25 +120,25 @@ void SetPsipredPrediction(BackboneScoreEnv& env) {
void CheckScorerThrows(const promod3::loop::BackboneList& bb_list,
BackboneScorerPtr scorer) {
// test check of chain_idx
- BOOST_CHECK_THROW(scorer->CalculateScore(bb_list, 214, -1), promod3::Error);
- BOOST_CHECK_THROW(scorer->CalculateScore(bb_list, 214, 1), promod3::Error);
+ BOOST_CHECK_THROW(scorer->CalculateScore(214, bb_list.size(), -1), promod3::Error);
+ BOOST_CHECK_THROW(scorer->CalculateScore(214, bb_list.size(), 1), promod3::Error);
// test check of start resnum
- BOOST_CHECK_THROW(scorer->CalculateScore(bb_list, 0, 0), promod3::Error);
- BOOST_CHECK_THROW(scorer->CalculateScore(bb_list, -1, 0), promod3::Error);
- BOOST_CHECK_THROW(scorer->CalculateScore(bb_list, 1000, 0), promod3::Error);
+ BOOST_CHECK_THROW(scorer->CalculateScore(0, bb_list.size(), 0), promod3::Error);
+ BOOST_CHECK_THROW(scorer->CalculateScore(-1, bb_list.size(), 0), promod3::Error);
+ BOOST_CHECK_THROW(scorer->CalculateScore(1000, bb_list.size(), 0), promod3::Error);
// the first resnum, where the bb_list doesn't fit anymore...
- BOOST_CHECK_THROW(scorer->CalculateScore(bb_list, 255, 0), promod3::Error);
+ BOOST_CHECK_THROW(scorer->CalculateScore(255, bb_list.size(), 0), promod3::Error);
// check that empty BB list gives score of 0 (always)
promod3::loop::BackboneList sub_bb_list;
- BOOST_CHECK_EQUAL(scorer->CalculateScore(sub_bb_list, 214, 0), Real(0));
+ BOOST_CHECK_EQUAL(scorer->CalculateScore(214, sub_bb_list.size(), 0), Real(0));
// check that border cases don't break
if (bb_list.size() > 1) {
sub_bb_list = bb_list.Extract(0, 1);
- BOOST_CHECK_NO_THROW(scorer->CalculateScore(sub_bb_list, 214, 0));
+ BOOST_CHECK_NO_THROW(scorer->CalculateScore(214, sub_bb_list.size(), 0));
sub_bb_list = bb_list.Extract(0, 2);
- BOOST_CHECK_NO_THROW(scorer->CalculateScore(sub_bb_list, 214, 0));
+ BOOST_CHECK_NO_THROW(scorer->CalculateScore(214, sub_bb_list.size(), 0));
}
}
@@ -148,12 +148,12 @@ void SanityCheckOverallScorer(BackboneOverallScorer& scorer,
const promod3::loop::BackboneList& bb_list,
const String& key) {
// ref = single score
- Real ref_score = single_scorer->CalculateScore(bb_list, 214, 0);
- Real score = scorer.Calculate(key, bb_list, 214, 0);
+ Real ref_score = single_scorer->CalculateScore(214, bb_list.size(), 0);
+ Real score = scorer.Calculate(key, 214, bb_list.size(), 0);
BOOST_CHECK_CLOSE(score, ref_score, Real(1e-3));
std::map<String, Real> weights;
weights[key] = 1.0;
- score = scorer.CalculateLinearCombination(weights, bb_list, 214, 0);
+ score = scorer.CalculateLinearCombination(weights, 214, bb_list.size(), 0);
BOOST_CHECK_CLOSE(score, ref_score, Real(1e-3));
}
@@ -509,17 +509,17 @@ BOOST_AUTO_TEST_CASE(qmean_comparison) {
std::vector<Real> ss_agreement_profile;
// compare the orig values
- cb_packing_scorer->CalculateScoreProfile(bb_list, 26, 0,
+ cb_packing_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
cb_packing_profile);
- cb_scorer->CalculateScoreProfile(bb_list, 26, 0,
+ cb_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
cb_profile);
- reduced_scorer->CalculateScoreProfile(bb_list, 26, 0,
+ reduced_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
reduced_profile);
- torsion_scorer->CalculateScoreProfile(bb_list, 26, 0,
+ torsion_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
torsion_profile);
- hbond_scorer->CalculateScoreProfile(bb_list, 26, 0,
+ hbond_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
hbond_profile);
- ss_agreement_scorer->CalculateScoreProfile(bb_list, 26, 0,
+ ss_agreement_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
ss_agreement_profile);
Real eps = 0.0001;
@@ -548,17 +548,18 @@ BOOST_AUTO_TEST_CASE(qmean_comparison) {
bb_list.ReplaceFragment(bb_list_helix, 0, false);
// compare the reset values
- cb_packing_scorer->CalculateScoreProfile(bb_list, 26, 0,
+ cb_packing_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
cb_packing_profile);
- cb_scorer->CalculateScoreProfile(bb_list, 26, 0,
+ cb_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
cb_profile);
- reduced_scorer->CalculateScoreProfile(bb_list, 26, 0,
+ reduced_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
reduced_profile);
- torsion_scorer->CalculateScoreProfile(bb_list, 26, 0,
+ torsion_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
torsion_profile);
- hbond_scorer->CalculateScoreProfile(bb_list, 26, 0,
+ hbond_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
hbond_profile);
- ss_agreement_scorer->CalculateScoreProfile(bb_list, 26, 0,
+
+ ss_agreement_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
ss_agreement_profile);
@@ -581,17 +582,17 @@ BOOST_AUTO_TEST_CASE(qmean_comparison) {
env.ClearEnvironment(83, 6);
// compare the delete values
- cb_packing_scorer->CalculateScoreProfile(bb_list, 26, 0,
+ cb_packing_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
cb_packing_profile);
- cb_scorer->CalculateScoreProfile(bb_list, 26, 0,
+ cb_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
cb_profile);
- reduced_scorer->CalculateScoreProfile(bb_list, 26, 0,
+ reduced_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
reduced_profile);
- torsion_scorer->CalculateScoreProfile(bb_list, 26, 0,
+ torsion_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
torsion_profile);
- hbond_scorer->CalculateScoreProfile(bb_list, 26, 0,
+ hbond_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
hbond_profile);
- ss_agreement_scorer->CalculateScoreProfile(bb_list, 26, 0,
+ ss_agreement_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
ss_agreement_profile);
@@ -611,6 +612,377 @@ BOOST_AUTO_TEST_CASE(qmean_comparison) {
}
}
+BOOST_AUTO_TEST_CASE(test_env_stashing) {
+
+ // Same test as above, but this time we use some stash / pop
+ // magic from the BackboneScoreEnv
+ BackboneScoreEnv env = GenerateBackboneEnv();
+ promod3::loop::BackboneList bb_list = GenerateTestBBList("data/3DEFA.pdb",
+ 26, 40);
+
+ CBPackingScorerPtr cb_packing_scorer = LoadCBPackingScorer();
+ CBetaScorerPtr cb_scorer = LoadCBetaScorer();
+ ReducedScorerPtr reduced_scorer = LoadReducedScorer();
+ HBondScorerPtr hbond_scorer = LoadHBondScorer();
+ SSAgreementScorerPtr ss_agreement_scorer = LoadSSAgreementScorer();
+ SetPsipredPrediction(env);
+ TorsionScorerPtr torsion_scorer = LoadTorsionScorer();
+
+ cb_packing_scorer->AttachEnvironment(env);
+ cb_scorer->AttachEnvironment(env);
+ reduced_scorer->AttachEnvironment(env);
+ hbond_scorer->AttachEnvironment(env);
+ ss_agreement_scorer->AttachEnvironment(env);
+ torsion_scorer->AttachEnvironment(env);
+
+ // orig means: the current environment
+ // reset means: 8 residues starting from resnum 26 are replaced by a helix as
+ // it is constructed in the default constructor of the
+ // BackboneList. Positioning is performed by superposing first
+ // helix residue on current residue 26
+
+ std::vector<Real> orig_cb_packing;
+ orig_cb_packing.push_back(-1.0796);
+ orig_cb_packing.push_back(-0.5889);
+ orig_cb_packing.push_back(-0.8458);
+ orig_cb_packing.push_back(-0.9514);
+ orig_cb_packing.push_back(-0.6117);
+ orig_cb_packing.push_back(-0.7100);
+ orig_cb_packing.push_back(-0.9639);
+ orig_cb_packing.push_back(-0.5944);
+ orig_cb_packing.push_back(-0.5149);
+ orig_cb_packing.push_back(-1.0095);
+ orig_cb_packing.push_back(-0.8998);
+ orig_cb_packing.push_back(-0.9830);
+ orig_cb_packing.push_back(-0.4551);
+ orig_cb_packing.push_back(-0.6673);
+
+ std::vector<Real> orig_cb;
+ orig_cb.push_back(-0.342035);
+ orig_cb.push_back(0.238565);
+ orig_cb.push_back(-0.410386);
+ orig_cb.push_back(0.340058);
+ orig_cb.push_back(-0.341305);
+ orig_cb.push_back(-0.0710439);
+ orig_cb.push_back(-0.0911381);
+ orig_cb.push_back(-0.0181578);
+ orig_cb.push_back(0.808517);
+ orig_cb.push_back(0.461206);
+ orig_cb.push_back(0.139356);
+ orig_cb.push_back(-0.482617);
+ orig_cb.push_back(-0.497085);
+ orig_cb.push_back(0.0492029);
+
+
+ std::vector<Real> orig_reduced;
+ orig_reduced.push_back(-3.55399);
+ orig_reduced.push_back(0.52569);
+ orig_reduced.push_back(-5.0972);
+ orig_reduced.push_back(-5.45686);
+ orig_reduced.push_back(-1.86842);
+ orig_reduced.push_back(-0.992303);
+ orig_reduced.push_back(0.543362);
+ orig_reduced.push_back(-3.62267);
+ orig_reduced.push_back(-0.268999);
+ orig_reduced.push_back(-4.16786);
+ orig_reduced.push_back(-2.89717);
+ orig_reduced.push_back(-4.37486);
+ orig_reduced.push_back(-7.23112);
+ orig_reduced.push_back(-10.5035);
+
+ std::vector<Real> orig_torsion;
+ orig_torsion.push_back(0.4922);
+ orig_torsion.push_back(0.4896);
+ orig_torsion.push_back(0.9574);
+ orig_torsion.push_back(-0.4831);
+ orig_torsion.push_back(-0.5004);
+ orig_torsion.push_back(-0.4875);
+ orig_torsion.push_back(-0.4739);
+ orig_torsion.push_back(-0.7882);
+ orig_torsion.push_back(0.3369);
+ orig_torsion.push_back(-0.8233);
+ orig_torsion.push_back(-0.4457);
+ orig_torsion.push_back(-0.4084);
+ orig_torsion.push_back(-1.1914);
+ orig_torsion.push_back(-1.1170);
+
+ std::vector<Real> orig_hbond;
+ orig_hbond.push_back(-1.5817);
+ orig_hbond.push_back(-1.4568);
+ orig_hbond.push_back(-2.7138);
+ orig_hbond.push_back(-2.9179);
+ orig_hbond.push_back(-1.4514);
+ orig_hbond.push_back(-1.1016);
+ orig_hbond.push_back(-2.6790);
+ orig_hbond.push_back(-0.4062);
+ orig_hbond.push_back(-0.9520);
+ orig_hbond.push_back(-0.9263);
+ orig_hbond.push_back(0.0000);
+ orig_hbond.push_back(-1.6571);
+ orig_hbond.push_back(-1.0162);
+ orig_hbond.push_back(-1.7538);
+
+ std::vector<Real> orig_ss_agreement;
+ orig_ss_agreement.push_back(1.0960);
+ orig_ss_agreement.push_back(1.0960);
+ orig_ss_agreement.push_back(1.0960);
+ orig_ss_agreement.push_back(0.9231);
+ orig_ss_agreement.push_back(0.9231);
+ orig_ss_agreement.push_back(0.6663);
+ orig_ss_agreement.push_back(0.2756);
+ orig_ss_agreement.push_back(0.4928);
+ orig_ss_agreement.push_back(0.7852);
+ orig_ss_agreement.push_back(0.7852);
+ orig_ss_agreement.push_back(-0.9896);
+ orig_ss_agreement.push_back(-0.1106);
+ orig_ss_agreement.push_back(1.3440);
+ orig_ss_agreement.push_back(1.4422);
+
+ std::vector<Real> reset_cb_packing;
+ reset_cb_packing.push_back(-1.0796);
+ reset_cb_packing.push_back(-0.4539);
+ reset_cb_packing.push_back(-1.0485);
+ reset_cb_packing.push_back(-0.9514);
+ reset_cb_packing.push_back(-0.6117);
+ reset_cb_packing.push_back(-0.7100);
+ reset_cb_packing.push_back(-0.9639);
+ reset_cb_packing.push_back(-0.9962);
+ reset_cb_packing.push_back(-0.5149);
+ reset_cb_packing.push_back(-0.7859);
+ reset_cb_packing.push_back(-0.8998);
+ reset_cb_packing.push_back(-0.9830);
+ reset_cb_packing.push_back(-0.4551);
+ reset_cb_packing.push_back(-0.6673);
+
+ std::vector<Real> reset_cb;
+ reset_cb.push_back(-0.32877);
+ reset_cb.push_back(0.371232);
+ reset_cb.push_back(-0.479475);
+ reset_cb.push_back(-0.369955);
+ reset_cb.push_back(-0.428097);
+ reset_cb.push_back(-0.0879088);
+ reset_cb.push_back(-0.43112);
+ reset_cb.push_back(-0.324168);
+ reset_cb.push_back(-0.0822164);
+ reset_cb.push_back(0.228412);
+ reset_cb.push_back(0.0871677);
+ reset_cb.push_back(-0.510713);
+ reset_cb.push_back(-0.497085);
+ reset_cb.push_back(0.0492029);
+
+ std::vector<Real> reset_reduced;
+ reset_reduced.push_back(-3.32411);
+ reset_reduced.push_back(1.60945);
+ reset_reduced.push_back(-1.03455);
+ reset_reduced.push_back(-4.02522);
+ reset_reduced.push_back(-0.961522);
+ reset_reduced.push_back(-0.435582);
+ reset_reduced.push_back(0.653847);
+ reset_reduced.push_back(1.73842);
+ reset_reduced.push_back(-0.985268);
+ reset_reduced.push_back(-2.20225);
+ reset_reduced.push_back(-2.83298);
+ reset_reduced.push_back(-4.28751);
+ reset_reduced.push_back(-7.23112);
+ reset_reduced.push_back(-10.5035);
+
+ std::vector<Real> reset_torsion;
+ reset_torsion.push_back(1.0555);
+ reset_torsion.push_back(0.5257);
+ reset_torsion.push_back(0.6707);
+ reset_torsion.push_back(-0.4831);
+ reset_torsion.push_back(-0.4842);
+ reset_torsion.push_back(-0.4924);
+ reset_torsion.push_back(0.1419);
+ reset_torsion.push_back(4.7277);
+ reset_torsion.push_back(4.5157);
+ reset_torsion.push_back(-0.8233);
+ reset_torsion.push_back(-0.4457);
+ reset_torsion.push_back(-0.4084);
+ reset_torsion.push_back(-1.1914);
+ reset_torsion.push_back(-1.1170);
+
+ std::vector<Real> reset_hbond;
+ reset_hbond.push_back(-1.3101);
+ reset_hbond.push_back(-0.9559);
+ reset_hbond.push_back(-1.2444);
+ reset_hbond.push_back(-1.8388);
+ reset_hbond.push_back(-0.9559);
+ reset_hbond.push_back(-0.9559);
+ reset_hbond.push_back(-0.9559);
+ reset_hbond.push_back(-0.4811);
+ reset_hbond.push_back(0.0000);
+ reset_hbond.push_back(-0.9263);
+ reset_hbond.push_back(0.0000);
+ reset_hbond.push_back(-1.6571);
+ reset_hbond.push_back(-1.0162);
+ reset_hbond.push_back(-1.7538);
+
+ std::vector<Real> reset_ss_agreement;
+ reset_ss_agreement.push_back(1.0960);
+ reset_ss_agreement.push_back(1.0960);
+ reset_ss_agreement.push_back(1.0960);
+ reset_ss_agreement.push_back(0.9231);
+ reset_ss_agreement.push_back(0.9231);
+ reset_ss_agreement.push_back(0.6663);
+ reset_ss_agreement.push_back(-0.7839);
+ reset_ss_agreement.push_back(0.5699);
+ reset_ss_agreement.push_back(0.7852);
+ reset_ss_agreement.push_back(0.7852);
+ reset_ss_agreement.push_back(-0.9896);
+ reset_ss_agreement.push_back(-0.1106);
+ reset_ss_agreement.push_back(1.3440);
+ reset_ss_agreement.push_back(1.4422);
+
+ std::vector<Real> cb_packing_profile;
+ std::vector<Real> cb_profile;
+ std::vector<Real> reduced_profile;
+ std::vector<Real> torsion_profile;
+ std::vector<Real> hbond_profile;
+ std::vector<Real> ss_agreement_profile;
+
+ // compare the orig values
+ cb_packing_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ cb_packing_profile);
+ cb_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ cb_profile);
+ reduced_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ reduced_profile);
+ torsion_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ torsion_profile);
+ hbond_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ hbond_profile);
+ ss_agreement_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ ss_agreement_profile);
+
+ Real eps = 0.0001;
+ for(uint i = 0; i < bb_list.size(); ++i) {
+ BOOST_CHECK_SMALL(std::abs(orig_cb_packing[i] -
+ cb_packing_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(orig_cb[i] -
+ cb_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(orig_reduced[i] -
+ reduced_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(orig_torsion[i] -
+ torsion_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(orig_hbond[i] -
+ hbond_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(orig_ss_agreement[i] -
+ ss_agreement_profile[i]), eps);
+ }
+
+
+ //////////////////////////////////////////////////////////////////////////
+ // PERFORM A STASH OPERATION AND CHECK, WHETHER THE SCORING RESULTS ARE //
+ // STILL THE SAME //
+ //////////////////////////////////////////////////////////////////////////
+
+ env.Stash(26, bb_list.size(), 0);
+
+ cb_packing_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ cb_packing_profile);
+ cb_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ cb_profile);
+ reduced_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ reduced_profile);
+ torsion_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ torsion_profile);
+ hbond_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ hbond_profile);
+ ss_agreement_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ ss_agreement_profile);
+
+ for(uint i = 0; i < bb_list.size(); ++i) {
+ BOOST_CHECK_SMALL(std::abs(orig_cb_packing[i] -
+ cb_packing_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(orig_cb[i] -
+ cb_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(orig_reduced[i] -
+ reduced_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(orig_torsion[i] -
+ torsion_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(orig_hbond[i] -
+ hbond_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(orig_ss_agreement[i] -
+ ss_agreement_profile[i]), eps);
+ }
+
+ // lets perform the described reset action
+ String helix_sequence = "FGKLKQKD";
+ promod3::loop::BackboneList bb_list_helix =
+ promod3::loop::BackboneList(helix_sequence);
+ geom::Mat4 transform = bb_list_helix.GetTransform(0, bb_list, 0);
+ bb_list_helix.ApplyTransform(transform);
+ env.SetEnvironment(bb_list_helix, 26);
+ bb_list.ReplaceFragment(bb_list_helix, 0, false);
+
+ // compare the reset values
+ cb_packing_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ cb_packing_profile);
+ cb_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ cb_profile);
+ reduced_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ reduced_profile);
+ torsion_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ torsion_profile);
+ hbond_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ hbond_profile);
+
+ ss_agreement_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ ss_agreement_profile);
+
+
+ for(uint i = 0; i < bb_list.size(); ++i) {
+ BOOST_CHECK_SMALL(std::abs(reset_cb_packing[i] -
+ cb_packing_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(reset_cb[i] -
+ cb_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(reset_reduced[i] -
+ reduced_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(reset_torsion[i] -
+ torsion_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(reset_hbond[i] -
+ hbond_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(reset_ss_agreement[i] -
+ ss_agreement_profile[i]), eps);
+ }
+
+ ////////////////////////////////////////////////////////////////////////////
+ // WE PERFORM A POP ACTION AND CHECK WHETHER THE SCORES CORRESPOND TO THE //
+ // ORIGINAL SCORES AGAIN //
+ ////////////////////////////////////////////////////////////////////////////
+
+ env.Pop();
+
+ cb_packing_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ cb_packing_profile);
+ cb_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ cb_profile);
+ reduced_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ reduced_profile);
+ torsion_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ torsion_profile);
+ hbond_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ hbond_profile);
+ ss_agreement_scorer->CalculateScoreProfile(26, bb_list.size(), 0,
+ ss_agreement_profile);
+
+ for(uint i = 0; i < bb_list.size(); ++i) {
+ BOOST_CHECK_SMALL(std::abs(orig_cb_packing[i] -
+ cb_packing_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(orig_cb[i] -
+ cb_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(orig_reduced[i] -
+ reduced_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(orig_torsion[i] -
+ torsion_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(orig_hbond[i] -
+ hbond_profile[i]), eps);
+ BOOST_CHECK_SMALL(std::abs(orig_ss_agreement[i] -
+ ss_agreement_profile[i]), eps);
+ }
+}
+
BOOST_AUTO_TEST_CASE(test_clash_score) {
BackboneScoreEnv env = GenerateBackboneEnv();
promod3::loop::BackboneList bb_list = GenerateTestBBList();
@@ -627,7 +999,7 @@ BOOST_AUTO_TEST_CASE(test_clash_score2) {
promod3::loop::BackboneList bb_list = GenerateTestBBList(pdb_name, 35, 40);
ClashScorerPtr scorer(new ClashScorer);
scorer->AttachEnvironment(env);
- Real clash_score = scorer->CalculateScore(bb_list, 35, 0);
+ Real clash_score = scorer->CalculateScore(35, bb_list.size(), 0);
BOOST_CHECK_CLOSE(clash_score, 0.107386872172, Real(1e-3));
}
@@ -640,7 +1012,7 @@ BOOST_AUTO_TEST_CASE(test_pairwise_score) {
//CheckScorer(bb_list, scorer, 0.0);
CheckScorerThrows(bb_list, scorer);
- BOOST_CHECK_CLOSE(scorer->CalculateScore(bb_list, 214, 0), 0.0,
+ BOOST_CHECK_CLOSE(scorer->CalculateScore(214, bb_list.size(), 0), 0.0,
Real(1e-3));
@@ -665,7 +1037,7 @@ BOOST_AUTO_TEST_CASE(test_pairwise_score) {
//scorer should give 0.0
//CheckScorer(bb_list, scorer, 0.0);
- BOOST_CHECK_CLOSE(scorer->CalculateScore(bb_list, 214, 0), 0.0,
+ BOOST_CHECK_CLOSE(scorer->CalculateScore(214, bb_list.size(), 0), 0.0,
Real(1e-3));
@@ -676,7 +1048,7 @@ BOOST_AUTO_TEST_CASE(test_pairwise_score) {
//CheckScorer(bb_list, scorer, 0.36586231);
- BOOST_CHECK_CLOSE(scorer->CalculateScore(bb_list, 214, 0), 0.36586231,
+ BOOST_CHECK_CLOSE(scorer->CalculateScore(214, bb_list.size(), 0), 0.36586231,
Real(1e-3));
@@ -684,7 +1056,7 @@ BOOST_AUTO_TEST_CASE(test_pairwise_score) {
//that should also have an influence
//CheckScorer(bb_list, scorer, 1.12155807);
- BOOST_CHECK_CLOSE(scorer->CalculateScore(bb_list, 214, 0), 1.12155807,
+ BOOST_CHECK_CLOSE(scorer->CalculateScore(214, bb_list.size(), 0), 1.12155807,
Real(1e-3));
//let's finally try to add a contact function
@@ -697,7 +1069,7 @@ BOOST_AUTO_TEST_CASE(test_pairwise_score) {
//compared to previous check
//CheckScorer(bb_list, scorer, 1.12155807);
- BOOST_CHECK_CLOSE(scorer->CalculateScore(bb_list, 214, 0), 1.12155807,
+ BOOST_CHECK_CLOSE(scorer->CalculateScore(214, bb_list.size(), 0), 1.12155807,
Real(1e-3));
@@ -706,7 +1078,7 @@ BOOST_AUTO_TEST_CASE(test_pairwise_score) {
//therefore doesn't influence the score
//CheckScorer(bb_list, scorer, 1.12155807);
- BOOST_CHECK_CLOSE(scorer->CalculateScore(bb_list, 214, 0), 1.12155807,
+ BOOST_CHECK_CLOSE(scorer->CalculateScore(214, bb_list.size(), 0), 1.12155807,
Real(1e-3));
env.ApplyPairwiseFunction(0,216,0,228,3); //the cb positions of residues 216
@@ -714,7 +1086,7 @@ BOOST_AUTO_TEST_CASE(test_pairwise_score) {
//therefore influences the score
//CheckScorer(bb_list, scorer, 6.67711401);
- BOOST_CHECK_CLOSE(scorer->CalculateScore(bb_list, 214, 0), 6.67711401,
+ BOOST_CHECK_CLOSE(scorer->CalculateScore(214, bb_list.size(), 0), 6.67711401,
Real(1e-3));
}
@@ -734,21 +1106,8 @@ BOOST_AUTO_TEST_CASE(test_density_score){
ost::img::MapHandle map = full_bb_list.ToDensity(10.0,geom::Vec3(1.0,1.0,1.0),
5.0, true);
- DensityScorerPtr scorer(new DensityScorer);
-
- //check if error gets thrown when no environment is attached
- BOOST_CHECK_THROW(scorer->CalculateScore(bb_list,5,0), promod3::Error);
-
- scorer->AttachEnvironment(env);
- //check if error gets thrown when no density map is set in the environment
- BOOST_CHECK_THROW(scorer->CalculateScore(bb_list,5,0), promod3::Error);
-
- env.SetMap(map, 5.0, true);
-
- //The CheckScorer function cannot be used, since the resnum and chain idx
- //checks are not permormed in the density scorer
- BOOST_CHECK_CLOSE(scorer->CalculateScore(bb_list, 1, 0), 0.352861583,
- Real(1e-3));
+ BOOST_CHECK_CLOSE(promod3::scoring::DensityScore(bb_list, map, 5.0, true),
+ 0.352861583, Real(1e-3));
}
BOOST_AUTO_TEST_CASE(test_overall_scorer) {
@@ -770,8 +1129,8 @@ BOOST_AUTO_TEST_CASE(test_overall_scorer) {
weights["torsion"] = 0.7;
weights["intercept"] = 1.0;
// check for errors
- BOOST_CHECK_THROW(scorer.Calculate("hbond", bb_list, 214, 0), promod3::Error);
- BOOST_CHECK_THROW(scorer.CalculateLinearCombination(weights, bb_list, 214, 0),
+ BOOST_CHECK_THROW(scorer.Calculate("hbond", 214, bb_list.size(), 0), promod3::Error);
+ BOOST_CHECK_THROW(scorer.CalculateLinearCombination(weights, 214, bb_list.size(), 0),
promod3::Error);
// attach and sanity check single scores
@@ -787,10 +1146,10 @@ BOOST_AUTO_TEST_CASE(test_overall_scorer) {
Real exp_score = 0.1 * -0.657944619656 + 0.2 * 0.0189835
+ 0.3 * -0.64918 + 0.4 * 0.0 + 0.5 * -0.482397
+ 0.6 * 0.498321 + 0.7 * -0.143922114538 + 1.0;
- BOOST_CHECK_CLOSE(scorer.CalculateLinearCombination(weights, bb_list, 214, 0),
+ BOOST_CHECK_CLOSE(scorer.CalculateLinearCombination(weights, 214, bb_list.size(), 0),
exp_score, Real(1e-3));
// check for more errors
- BOOST_CHECK_THROW(scorer.Calculate("meh", bb_list, 214, 0), promod3::Error);
+ BOOST_CHECK_THROW(scorer.Calculate("meh", 214, bb_list.size(), 0), promod3::Error);
}
BOOST_AUTO_TEST_SUITE_END();