From b41397e839251d3b5b925f7162248414ddc5eb12 Mon Sep 17 00:00:00 2001
From: Marco Biasini <mvbiasini@gmail.com>
Date: Tue, 29 Jan 2013 10:25:03 +0100
Subject: [PATCH] removed AtomsBegin/End, ResiduesBegin/End from entity
they added very little benefits and were in general very slow.
The call sites of AtomsBegin have been replaced by a triple
for loop, ResiduesBegin by a double for loop. While that adds
some typing overhead, it gives a speed up compared to the old
version. It also happens to fix a problem noticed by Niklaus.
the superposition code has also been refactored to reduce
redundancy and make it easier to read.
---
modules/conop/src/model_check.cc | 21 +-
modules/conop/src/model_check.hh | 14 +-
modules/conop/tests/test_heuristic_builder.cc | 9 +-
.../conop/tests/test_rule_based_builder.cc | 6 +-
modules/mol/alg/pymod/export_svd_superpose.cc | 4 -
modules/mol/alg/src/clash_score.cc | 45 +-
modules/mol/alg/src/lddt.cc | 98 ++--
modules/mol/alg/src/svd_superpose.cc | 449 ++++++------------
modules/mol/alg/src/svd_superpose.hh | 41 +-
modules/mol/alg/tests/test_superposition.cc | 10 +-
modules/mol/base/src/CMakeLists.txt | 3 -
modules/mol/base/src/chain_handle.cc | 70 +--
modules/mol/base/src/chain_handle.hh | 23 -
modules/mol/base/src/coord_frame.cc | 11 +-
modules/mol/base/src/entity_handle.cc | 87 +---
modules/mol/base/src/entity_handle.hh | 33 --
modules/mol/base/src/entity_view.cc | 157 +++---
modules/mol/base/src/entity_view.hh | 15 -
modules/mol/base/src/iterator.cc | 198 --------
modules/mol/base/src/iterator.hh | 237 ---------
modules/mol/base/src/iterator_fw.hh | 35 --
modules/mol/base/src/mol.hh | 1 -
modules/mol/base/src/query_state.cc | 20 +-
modules/mol/base/src/residue_handle.cc | 32 +-
modules/mol/base/src/residue_handle.hh | 12 -
modules/mol/base/src/view_op.cc | 53 ++-
modules/mol/base/tests/CMakeLists.txt | 1 -
27 files changed, 463 insertions(+), 1222 deletions(-)
delete mode 100644 modules/mol/base/src/iterator.cc
delete mode 100644 modules/mol/base/src/iterator.hh
delete mode 100644 modules/mol/base/src/iterator_fw.hh
diff --git a/modules/conop/src/model_check.cc b/modules/conop/src/model_check.cc
index cb5d710b6..184112c17 100644
--- a/modules/conop/src/model_check.cc
+++ b/modules/conop/src/model_check.cc
@@ -1,4 +1,3 @@
-#include <ost/mol/iterator.hh>
#include "model_check.hh"
#include "amino_acids.hh"
@@ -8,8 +7,8 @@ namespace ost { namespace conop {
void Checker::CheckForCompleteness(bool require_hydrogens)
{
- for (ResidueHandleIter i=ent_.ResiduesBegin(), e=ent_.ResiduesEnd();
- i!=e; ++i) {
+ for (ResidueHandleList::const_iterator i=residues_.begin(),
+ e = residues_.end(); i!=e; ++i) {
ResidueHandle res=*i;
String anames="";
CompoundPtr compound=lib_->FindCompound(res.GetName(),Compound::PDB);
@@ -52,8 +51,8 @@ void Checker::CheckForCompleteness(bool require_hydrogens)
mol::AtomHandleList Checker::GetHydrogens()
{
AtomHandleList hydlist;
- for (ResidueHandleIter i=ent_.ResiduesBegin(), e=ent_.ResiduesEnd();
- i!=e; ++i) {
+ for (ResidueHandleList::const_iterator i=residues_.begin(),
+ e = residues_.end(); i!=e; ++i) {
ResidueHandle res=*i;
CompoundPtr compound=lib_->FindCompound(res.GetName(),Compound::PDB);
if (!compound) {
@@ -76,8 +75,8 @@ mol::AtomHandleList Checker::GetHydrogens()
mol::AtomHandleList Checker::GetZeroOccupancy()
{
AtomHandleList zerolist;
- for (ResidueHandleIter i=ent_.ResiduesBegin(), e=ent_.ResiduesEnd();
- i!=e; ++i) {
+ for (ResidueHandleList::const_iterator i=residues_.begin(),
+ e = residues_.end(); i!=e; ++i) {
ResidueHandle res=*i;
AtomHandleList atoms=res.GetAtomList();
for (AtomHandleList::const_iterator j=atoms.begin(), e2=atoms.end(); j!=e2; ++j) {
@@ -91,8 +90,8 @@ mol::AtomHandleList Checker::GetZeroOccupancy()
void Checker::CheckForNonStandard()
{
- for (ResidueHandleIter i=ent_.ResiduesBegin(), e=ent_.ResiduesEnd();
- i!=e; ++i) {
+ for (ResidueHandleList::const_iterator i=residues_.begin(),
+ e = residues_.end(); i!=e; ++i) {
ResidueHandle res=*i;
CompoundPtr compound=lib_->FindCompound(res.GetName(),Compound::PDB);
if (!compound) {
@@ -113,8 +112,8 @@ void Checker::CheckForNonStandard()
void Checker::CheckForUnknownAtoms()
{
- for (ResidueHandleIter i=ent_.ResiduesBegin(), e=ent_.ResiduesEnd();
- i!=e; ++i) {
+ for (ResidueHandleList::const_iterator i=residues_.begin(),
+ e = residues_.end(); i!=e; ++i) {
ResidueHandle res=*i;
String anames="";
CompoundPtr compound=lib_->FindCompound(res.GetName(),Compound::PDB);
diff --git a/modules/conop/src/model_check.hh b/modules/conop/src/model_check.hh
index 9053c2391..3e56f45e9 100644
--- a/modules/conop/src/model_check.hh
+++ b/modules/conop/src/model_check.hh
@@ -10,8 +10,9 @@ class DLLEXPORT_OST_CONOP Checker {
public:
Checker(CompoundLibPtr lib, const mol::EntityHandle& ent,
DiagEngine& diags): lib_(lib), ent_(ent), diags_(diags),
- checked_unk_res_(false)
- {}
+ checked_unk_res_(false),
+ residues_(ent_.GetResidueList())
+ { }
void CheckForUnknownAtoms();
void CheckForCompleteness(bool require_hydrogens=false);
void CheckForNonStandard();
@@ -20,10 +21,11 @@ public:
const std::vector<Diag*>& GetDiags() const { return diags_.GetDiags(); }
private:
- CompoundLibPtr lib_;
- mol::EntityHandle ent_;
- DiagEngine& diags_;
- bool checked_unk_res_;
+ CompoundLibPtr lib_;
+ mol::EntityHandle ent_;
+ DiagEngine& diags_;
+ bool checked_unk_res_;
+ mol::ResidueHandleList residues_;
};
}} /* ost::conop */
diff --git a/modules/conop/tests/test_heuristic_builder.cc b/modules/conop/tests/test_heuristic_builder.cc
index bca1b67ee..cd13802ca 100644
--- a/modules/conop/tests/test_heuristic_builder.cc
+++ b/modules/conop/tests/test_heuristic_builder.cc
@@ -166,7 +166,8 @@ BOOST_AUTO_TEST_CASE(name_based_connect)
ResidueHandle ile=make_leu(c);
ResidueHandle arg=make_arg(c);
HeuristicBuilder heuristic_builder;
- for (AtomHandleIter i=e.AtomsBegin(),x=e.AtomsEnd(); i!=x; ++i) {
+ AtomHandleList atoms =e .GetAtomList();
+ for (AtomHandleList::const_iterator i = atoms.begin(), e = atoms.end(); i!=e; ++i ){
heuristic_builder.FillAtomProps(*i);
}
@@ -175,7 +176,8 @@ BOOST_AUTO_TEST_CASE(name_based_connect)
ResidueHandle dile=make_defective_leu(dc);
HeuristicBuilder dheuristic_builder;
dheuristic_builder.SetBondFeasibilityCheck(false);
- for (AtomHandleIter i=de.AtomsBegin(),x=de.AtomsEnd(); i!=x; ++i) {
+ atoms = de.GetAtomList();
+ for (AtomHandleList::const_iterator i = atoms.begin(), e = atoms.end(); i!=e; ++i ){
dheuristic_builder.FillAtomProps(*i);
}
@@ -201,7 +203,8 @@ BOOST_AUTO_TEST_CASE(test_assign_torsions){
a2.SetChemClass(ChemClass(ChemClass::L_PEPTIDE_LINKING));
l3.SetChemClass(ChemClass(ChemClass::L_PEPTIDE_LINKING));
HeuristicBuilder heuristic_builder;
- for (AtomHandleIter i=e.AtomsBegin(),x=e.AtomsEnd(); i!=x; ++i) {
+ AtomHandleList atoms = e.GetAtomList();
+ for (AtomHandleList::const_iterator i = atoms.begin(), e = atoms.end(); i!=e; ++i ){
heuristic_builder.FillAtomProps(*i);
}
heuristic_builder.ConnectAtomsOfResidue(l1);
diff --git a/modules/conop/tests/test_rule_based_builder.cc b/modules/conop/tests/test_rule_based_builder.cc
index bc89d94a5..a287bc124 100644
--- a/modules/conop/tests/test_rule_based_builder.cc
+++ b/modules/conop/tests/test_rule_based_builder.cc
@@ -394,11 +394,13 @@ BOOST_AUTO_TEST_CASE(nucleotide_based_connect)
ResidueHandle du2=make_defective_uracil2(dc);
- for (AtomHandleIter i=e.AtomsBegin(),x=e.AtomsEnd(); i!=x; ++i) {
+ AtomHandleList atoms = e.GetAtomList();
+ for (AtomHandleList::const_iterator i = atoms.begin(), e = atoms.end(); i!=e; ++i ){
rb_builder.FillAtomProps(*i);
}
- for (AtomHandleIter i=de.AtomsBegin(),x=de.AtomsEnd(); i!=x; ++i) {
+ atoms = de.GetAtomList();
+ for (AtomHandleList::const_iterator i = atoms.begin(), e = atoms.end(); i!=e; ++i ){
drb_builder.FillAtomProps(*i);
}
diff --git a/modules/mol/alg/pymod/export_svd_superpose.cc b/modules/mol/alg/pymod/export_svd_superpose.cc
index a02b7ff63..443ca71bc 100644
--- a/modules/mol/alg/pymod/export_svd_superpose.cc
+++ b/modules/mol/alg/pymod/export_svd_superpose.cc
@@ -25,7 +25,6 @@
#include <ost/geom/mat4.hh>
#include <ost/mol/alg/svd_superpose.hh>
#include <ost/mol/entity_handle.hh>
-#include <ost/mol/iterator.hh>
#include <boost/python/suite/indexing/vector_indexing_suite.hpp>
@@ -51,8 +50,5 @@ void export_svdSuperPose()
def("SuperposeSVD", sup1, (arg("apply_transform")=true));
def("SuperposeSVD", sup2, (arg("apply_transform")=true));
def("CalculateRMSD", &CalculateRMSD, (arg("transformation")=geom::Mat4()));
- def("IterativeSuperposition", &IterativeSuperposition, (arg("ncycles")=200,
- arg("dist_thres")=4.0,
- arg("apply_transform")=true));
}
diff --git a/modules/mol/alg/src/clash_score.cc b/modules/mol/alg/src/clash_score.cc
index 81e97c632..3bac0cd84 100644
--- a/modules/mol/alg/src/clash_score.cc
+++ b/modules/mol/alg/src/clash_score.cc
@@ -18,7 +18,8 @@
//------------------------------------------------------------------------------
#include <ost/mol/entity_view.hh>
#include <ost/mol/atom_view.hh>
-#include <ost/mol/iterator.hh>
+#include <ost/mol/residue_view.hh>
+#include <ost/mol/chain_view.hh>
#include "clash_score.hh"
namespace ost { namespace mol { namespace alg {
@@ -64,19 +65,55 @@ Real StericEnergy(const geom::Vec3& pos1, Real r1,
Real ClashScore(const EntityView& ent_a, const EntityView& ent_b)
{
- return do_clash_score<EntityView, AtomViewIter>(ent_a, ent_b);
+ Real energy=0.0;
+ for (ChainViewList::const_iterator ci = ent_a.GetChainList().begin(),
+ ce = ent_a.GetChainList().end(); ci != ce; ++ci) {
+ for (ResidueViewList::const_iterator ri = ci->GetResidueList().begin(),
+ re = ci->GetResidueList().end(); ri != re; ++re) {
+ for (AtomViewList::const_iterator ai = ri->GetAtomList().begin(),
+ ae = ri->GetAtomList().end(); ai!= ae; ++ai) {
+ AtomViewList clashees=ent_b.FindWithin(ai->GetPos(),
+ ai->GetRadius()+1.7);
+
+ for (AtomViewList::iterator j=clashees.begin(),
+ e2=clashees.end(); j!=e2; ++j) {
+ energy+=StericEnergy((*j).GetPos(), (*j).GetRadius()-0.25,
+ ai->GetPos(), ai->GetRadius()-0.25);
+ }
+ }
+ }
+ }
+ return energy;
}
Real ClashScore(const EntityHandle& ent_a, const EntityView& ent_b)
{
- return do_clash_score<EntityHandle, AtomHandleIter>(ent_a, ent_b);
+ Real energy=0.0;
+ for (ChainViewList::const_iterator ci = ent_b.GetChainList().begin(),
+ ce = ent_b.GetChainList().end(); ci != ce; ++ci) {
+ for (ResidueViewList::const_iterator ri = ci->GetResidueList().begin(),
+ re = ci->GetResidueList().end(); ri != re; ++re) {
+ for (AtomViewList::const_iterator ai = ri->GetAtomList().begin(),
+ ae = ri->GetAtomList().end(); ai!= ae; ++ai) {
+ AtomHandleList clashees=ent_a.FindWithin(ai->GetPos(),
+ ai->GetRadius()+1.7);
+
+ for (AtomHandleList::iterator j=clashees.begin(),
+ e2=clashees.end(); j!=e2; ++j) {
+ energy+=StericEnergy((*j).GetPos(), (*j).GetRadius()-0.25,
+ ai->GetPos(), ai->GetRadius()-0.25);
+ }
+ }
+ }
+ }
+ return energy;
}
Real ClashScore(const AtomHandle& atom, const EntityView& ent_b)
{
Real energy=0.0;
AtomViewList clashees=ent_b.FindWithin(atom.GetPos(),
- atom.GetRadius()+2.0);
+ atom.GetRadius()+2.0);
for (AtomViewList::iterator j=clashees.begin(),
e2=clashees.end(); j!=e2; ++j) {
energy+=StericEnergy((*j).GetPos(), (*j).GetRadius(),
diff --git a/modules/mol/alg/src/lddt.cc b/modules/mol/alg/src/lddt.cc
index ebf1a6e89..1787018b2 100644
--- a/modules/mol/alg/src/lddt.cc
+++ b/modules/mol/alg/src/lddt.cc
@@ -32,7 +32,6 @@
#include <ost/conop/conop.hh>
#include <ost/string_ref.hh>
#include <ost/conop/amino_acids.hh>
-#include <ost/mol/iterator.hh>
#include <ost/platform.hh>
#include <ost/log.hh>
#include <ost/mol/alg/consistency_checks.hh>
@@ -417,55 +416,62 @@ int main (int argc, char **argv)
} else {
std::cout << "Chain\tResName\tResNum\tAsses.\tScore\t(Conserved/Total, over " << cutoffs.size() << " thresholds)" << std::endl;
}
- for (ResidueViewIter rit=outv.ResiduesBegin();rit!=outv.ResiduesEnd();++rit){
- ResidueView ritv=*rit;
- ResNum rnum = ritv.GetNumber();
- bool assessed = false;
- String assessed_string="No";
- bool quality_problems = false;
- String quality_problems_string="No";
- Real lddt_local = -1;
- String lddt_local_string="-";
- int conserved_dist = -1;
- int total_dist = -1;
- String dist_string = "-";
- if (is_resnum_in_globalrdmap(rnum,glob_dist_list)) {
- assessed = true;
- assessed_string="Yes";
- }
- if (ritv.HasProp("stereo_chemical_violation_sidechain") || ritv.HasProp("steric_clash_sidechain")) {
- quality_problems = true;
- quality_problems_string="Yes";
- }
- if (ritv.HasProp("stereo_chemical_violation_backbone") || ritv.HasProp("steric_clash_backbone")) {
- quality_problems = true;
- quality_problems_string="Yes+";
- }
+ for (ChainViewList::const_iterator ci = outv.GetChainList().begin(),
+ ce = outv.GetChainList().end(); ci != ce; ++ci) {
+ for (ResidueViewList::const_iterator rit = ci->GetResidueList().begin(),
+ re = ci->GetResidueList().end(); rit != re; ++rit) {
+
+ ResidueView ritv=*rit;
+ ResNum rnum = ritv.GetNumber();
+ bool assessed = false;
+ String assessed_string="No";
+ bool quality_problems = false;
+ String quality_problems_string="No";
+ Real lddt_local = -1;
+ String lddt_local_string="-";
+ int conserved_dist = -1;
+ int total_dist = -1;
+ String dist_string = "-";
+ if (is_resnum_in_globalrdmap(rnum,glob_dist_list)) {
+ assessed = true;
+ assessed_string="Yes";
+ }
+ if (ritv.HasProp("stereo_chemical_violation_sidechain") ||
+ ritv.HasProp("steric_clash_sidechain")) {
+ quality_problems = true;
+ quality_problems_string="Yes";
+ }
+ if (ritv.HasProp("stereo_chemical_violation_backbone") ||
+ ritv.HasProp("steric_clash_backbone")) {
+ quality_problems = true;
+ quality_problems_string="Yes+";
+ }
- if (assessed==true) {
- if (ritv.HasProp(label)) {
- lddt_local=ritv.GetFloatProp(label);
- std::stringstream stkeylddt;
- stkeylddt << std::fixed << std::setprecision(4) << lddt_local;
- lddt_local_string=stkeylddt.str();
- conserved_dist=ritv.GetIntProp(label+"_conserved");
- total_dist=ritv.GetIntProp(label+"_total");
- std::stringstream stkeydist;
- stkeydist << "("<< conserved_dist << "/" << total_dist << ")";
- dist_string=stkeydist.str();
+ if (assessed==true) {
+ if (ritv.HasProp(label)) {
+ lddt_local=ritv.GetFloatProp(label);
+ std::stringstream stkeylddt;
+ stkeylddt << std::fixed << std::setprecision(4) << lddt_local;
+ lddt_local_string=stkeylddt.str();
+ conserved_dist=ritv.GetIntProp(label+"_conserved");
+ total_dist=ritv.GetIntProp(label+"_total");
+ std::stringstream stkeydist;
+ stkeydist << "("<< conserved_dist << "/" << total_dist << ")";
+ dist_string=stkeydist.str();
+ } else {
+ lddt_local = 0;
+ lddt_local_string="0.0000";
+ conserved_dist = 0;
+ total_dist = 0;
+ dist_string="(0/0)";
+ }
+ }
+ if (structural_checks) {
+ std::cout << ritv.GetChain() << "\t" << ritv.GetName() << "\t" << ritv.GetNumber() << '\t' << assessed_string << '\t' << quality_problems_string << '\t' << lddt_local_string << "\t" << dist_string << std::endl;
} else {
- lddt_local = 0;
- lddt_local_string="0.0000";
- conserved_dist = 0;
- total_dist = 0;
- dist_string="(0/0)";
+ std::cout << ritv.GetChain() << "\t" << ritv.GetName() << "\t" << ritv.GetNumber() << '\t' << assessed_string << '\t' << lddt_local_string << "\t" << dist_string << std::endl;
}
}
- if (structural_checks) {
- std::cout << ritv.GetChain() << "\t" << ritv.GetName() << "\t" << ritv.GetNumber() << '\t' << assessed_string << '\t' << quality_problems_string << '\t' << lddt_local_string << "\t" << dist_string << std::endl;
- } else {
- std::cout << ritv.GetChain() << "\t" << ritv.GetName() << "\t" << ritv.GetNumber() << '\t' << assessed_string << '\t' << lddt_local_string << "\t" << dist_string << std::endl;
- }
}
std::cout << std::endl;
}
diff --git a/modules/mol/alg/src/svd_superpose.cc b/modules/mol/alg/src/svd_superpose.cc
index f6b0a353d..56c0796c6 100644
--- a/modules/mol/alg/src/svd_superpose.cc
+++ b/modules/mol/alg/src/svd_superpose.cc
@@ -17,9 +17,6 @@
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
//------------------------------------------------------------------------------
-/*
- * Author Juergen Haas
- */
#include <stdexcept>
#include <iostream>
#include <boost/bind.hpp>
@@ -34,9 +31,13 @@
#include <ost/geom/vec3.hh>
#include <ost/mol/alg/svd_superpose.hh>
#include <ost/mol/xcs_editor.hh>
+#include <ost/mol/residue_handle.hh>
+#include <ost/mol/residue_view.hh>
+#include <ost/mol/chain_view.hh>
+#include <ost/mol/chain_handle.hh>
#include <ost/mol/view_op.hh>
#include <ost/mol/atom_view.hh>
-#include <ost/mol/iterator.hh>
+
namespace ost { namespace mol { namespace alg {
@@ -47,50 +48,36 @@ typedef Eigen::Matrix<Real,1,3> ERVec3;
typedef Eigen::Matrix<Real,Eigen::Dynamic,Eigen::Dynamic> EMatX;
typedef Eigen::Matrix<Real,1,Eigen::Dynamic> ERVecX;
-Real CalculateRMSD(const mol::EntityView& ev1,
- const mol::EntityView& ev2,
- const geom::Mat4& transformation) {
- mol::EntityView evt;
- mol::EntityHandle eht, eh1;
- int counta=ev1.GetAtomCount();
- int countb=ev2.GetAtomCount();
- if (counta!=countb){
- throw Error("atom counts of the two views are not equal");
- }
- eht=ev1.GetHandle();
- int natoms=ev1.GetAtomCount();
+
+
+
+
+Real calc_rmsd_for_point_lists(const std::vector<geom::Vec3>& points1,
+ const std::vector<geom::Vec3>& points2,
+ const geom::Mat4& transformation)
+{
+
Real rmsd=0.0;
- mol::AtomViewIter a_ev1=ev1.AtomsBegin();
- mol::AtomViewIter a_ev1_end=ev1.AtomsEnd();
- mol::AtomViewIter a_ev2=ev2.AtomsBegin();
- for (int counter=0; a_ev1_end!=a_ev1; ++a_ev1, ++a_ev2, ++counter) {
- mol::AtomView av1=*a_ev1;
- geom::Vec3 apos=av1.GetPos();
- mol::AtomView av2=*a_ev2;
- geom::Vec3 bpos=av2.GetPos();
- Real val=geom::Length2(geom::Vec3(transformation*geom::Vec4(apos))-bpos);
+ std::vector<geom::Vec3>::const_iterator a_ev1=points1.begin();
+ std::vector<geom::Vec3>::const_iterator a_ev1_end=points1.end();
+ std::vector<geom::Vec3>::const_iterator a_ev2=points2.begin();
+ for (; a_ev1_end!=a_ev1; ++a_ev1, ++a_ev2) {
+ Real val=geom::Length2(geom::Vec3(transformation*geom::Vec4(*a_ev1))-*a_ev2);
rmsd+=val;
}
- rmsd=sqrt(rmsd/natoms);
+ rmsd=sqrt(rmsd/points1.size());
return rmsd;
}
-
-
-
Real calc_rmsd_for_atom_lists(const mol::AtomViewList& atoms1,
- const mol::AtomViewList& atoms2,
- const geom::Mat4& transformation)
+ const mol::AtomViewList& atoms2,
+ const geom::Mat4& transformation)
{
- mol::EntityView evt;
- mol::EntityHandle eht, eh1;
-// convert view ev1 to handle
- int natoms=static_cast<int>(atoms1.size());
Real rmsd=0.0;
mol::AtomViewList::const_iterator a_ev1=atoms1.begin();
mol::AtomViewList::const_iterator a_ev1_end=atoms1.end();
mol::AtomViewList::const_iterator a_ev2=atoms2.begin();
- for (int counter=0; a_ev1_end!=a_ev1; ++a_ev1, ++a_ev2, ++counter) {
+ for (; a_ev1_end!=a_ev1; ++a_ev1, ++a_ev2) {
mol::AtomView av1=*a_ev1;
geom::Vec3 apos=av1.GetPos();
mol::AtomView av2=*a_ev2;
@@ -98,335 +85,215 @@ Real calc_rmsd_for_atom_lists(const mol::AtomViewList& atoms1,
Real val=geom::Length2(geom::Vec3(transformation*geom::Vec4(apos))-bpos);
rmsd+=val;
}
- rmsd=sqrt(rmsd/natoms);
+ rmsd=sqrt(rmsd/atoms1.size());
return rmsd;
}
-class SuperposerSVDImpl {
-public:
- SuperposerSVDImpl(int natoms, bool alloc_atoms):
- natoms_(natoms), alloc_atoms_(alloc_atoms)
- {
-
- ERVec3 avg1_=EVec3::Zero();
- ERVec3 avg2_=EVec3::Zero();
-
- if (alloc_atoms_) {
- atoms1_=EMatX::Zero(natoms,3);
- atoms2_=EMatX::Zero(natoms,3);
- }
-
- }
-private:
- int natoms_;
- bool alloc_atoms_;
- EMatX atoms1_;
- EMatX atoms2_;
-public:
- void EntToMatrices(const mol::EntityView& ev1, const mol::EntityView& ev2);
- geom::Vec3 EigenVec3ToVec3(const EVec3 &vec);
- geom::Mat3 EigenMat3ToMat3(const EMat3 &mat);
- EVec3 Vec3ToEigenRVec(const geom::Vec3 &vec);
- EVec3 Vec3ToEigenVec(const geom::Vec3 &vec);
- EMatX SubtractVecFromMatrixRows(EMatX Mat,
- ERVecX Vec);
- SuperpositionResult Run(const mol::AtomViewList& atoms1,
- const mol::AtomViewList& atoms2);
- SuperpositionResult Run(const std::vector<geom::Vec3>& pl1,
- const std::vector<geom::Vec3>& pl2);
- SuperpositionResult DoSuperposition();
- EIGEN_MAKE_ALIGNED_OPERATOR_NEW //needed only if *fixed* sized matrices are
- //members of classes
-};
-
-SuperposerSVD::~SuperposerSVD()
-{
- assert(impl_);
- delete impl_;
-}
-
+Real CalculateRMSD(const mol::EntityView& ev1,
+ const mol::EntityView& ev2,
+ const geom::Mat4& transformation) {
-SuperposerSVD::SuperposerSVD(int natoms, bool alloc_atoms):
- impl_(new SuperposerSVDImpl(natoms, alloc_atoms))
-{
-
+ return calc_rmsd_for_atom_lists(ev1.GetAtomList(), ev2.GetAtomList(),
+ transformation);
}
-geom::Vec3 SuperposerSVDImpl::EigenVec3ToVec3(const EVec3 &vec)
+geom::Vec3 EigenVec3ToVec3(const EVec3 &vec)
{
return geom::Vec3(vec.data());
}
-geom::Mat3 SuperposerSVDImpl::EigenMat3ToMat3(const EMat3 &mat)
+geom::Mat3 EigenMat3ToMat3(const EMat3 &mat)
{
return geom::Mat3(mat.data());
}
-EVec3 SuperposerSVDImpl::Vec3ToEigenRVec(const geom::Vec3 &vec)
+EVec3 Vec3ToEigenRVec(const geom::Vec3 &vec)
{
return EVec3(&vec[0]);
}
-EVec3 SuperposerSVDImpl::Vec3ToEigenVec(const geom::Vec3 &vec)
+EVec3 Vec3ToEigenVec(const geom::Vec3 &vec)
{
return EVec3(&vec[0]);
}
-EMatX SuperposerSVDImpl::SubtractVecFromMatrixRows(EMatX Mat,
- ERVecX Vec)
+EMatX MatrixShiftedBy(EMatX mat, ERVecX vec)
{
+ EMatX result = mat;
+ for (int row=0; row<mat.rows();++row) {
+ result.row(row) -= vec;
+ }
+ return result;
+}
- for (int row=0; row<Mat.rows();++row) {
- Mat.row(row)-=Vec;
+class MeanSquareMinimizerImpl {
+public:
+ MeanSquareMinimizerImpl(int n_atoms, bool alloc_atoms):
+ n_atoms_(n_atoms), alloc_atoms_(alloc_atoms)
+ {
+ if (alloc_atoms_) {
+ atoms1_=EMatX::Zero(n_atoms_, 3);
+ atoms2_=EMatX::Zero(n_atoms_, 3);
+ }
+ }
+
+ void SetRefPos(size_t index, const geom::Vec3& pos) {
+ atoms2_.row(index) = ERVec3(Vec3ToEigenVec(pos));
+ }
+
+ void SetPos(size_t index, const geom::Vec3& pos) {
+ atoms1_.row(index) = ERVec3(Vec3ToEigenVec(pos));
}
- return Mat;
+
+ SuperpositionResult MinimizeOnce() const;
+private:
+ int n_atoms_;
+ bool alloc_atoms_;
+ EMatX atoms1_;
+ EMatX atoms2_;
+};
+
+
+MeanSquareMinimizer MeanSquareMinimizer::FromAtomLists(const AtomViewList& atoms,
+ const AtomViewList& atoms_ref) {
+ int n_atoms = atoms.size();
+ if (n_atoms != atoms_ref.size()) {
+ throw Error("atom counts do not match");
+ }
+ if (n_atoms<3) {
+ throw Error("at least 3 atoms are required for superposition");
+ }
+ MeanSquareMinimizer msm;
+ msm.impl_ = new MeanSquareMinimizerImpl(n_atoms, true);
+
+ for (size_t i = 0; i < atoms.size(); ++i ) {
+ msm.impl_->SetRefPos(i, atoms_ref[i].GetPos());
+ msm.impl_->SetPos(i, atoms[i].GetPos());
+ }
+ return msm;
}
-void SuperposerSVDImpl::EntToMatrices(const mol::EntityView& ev1,
- const mol::EntityView& ev2)
-{
- //Now iterate over atoms in entities and extract coords into Nx3 matrices
- //for superposition
- mol::AtomViewIter a_ev1=ev1.AtomsBegin();
- mol::AtomViewIter a_ev1_end=ev1.AtomsEnd();
- mol::AtomViewIter a_ev2=ev2.AtomsBegin();
- ERVec3 e;
- geom::Vec3 pos;
- mol::AtomView av1, av2;
-
- for (int counter=0; a_ev1_end!=a_ev1; ++a_ev1, ++a_ev2, ++counter) {
- av1=*a_ev1;
- pos=av1.GetPos();
- e=SuperposerSVDImpl::Vec3ToEigenVec(pos);
- atoms1_.row(counter)=e ;
- av2=*a_ev2;
- pos=av2.GetPos();
- e=SuperposerSVDImpl::Vec3ToEigenVec(pos);
- atoms2_.row(counter)=e;
+MeanSquareMinimizer MeanSquareMinimizer::FromPointLists(const std::vector<geom::Vec3>& points,
+ const std::vector<geom::Vec3>& points_ref) {
+ int n_points = points.size();
+ if (n_points != points.size()) {
+ throw Error("point counts do not match");
}
+ if (n_points<3) {
+ throw Error("at least 3 points are required for superposition");
+ }
+ MeanSquareMinimizer msm;
+ msm.impl_ = new MeanSquareMinimizerImpl(n_points, true);
+
+ for (size_t i = 0; i < points.size(); ++i ) {
+ msm.impl_->SetRefPos(i, points_ref[i]);
+ msm.impl_->SetPos(i, points[i]);
+ }
+ return msm;
}
+SuperpositionResult MeanSquareMinimizerImpl::MinimizeOnce() const {
+ ERVec3 avg1 = atoms1_.colwise().sum()/atoms1_.rows();
+ ERVec3 avg2 = atoms2_.colwise().sum()/atoms2_.rows();
-SuperpositionResult SuperposerSVDImpl::DoSuperposition()
-{
- ERVec3 avg1_=atoms1_.colwise().sum()/atoms1_.rows();
- ERVec3 avg2_=atoms2_.colwise().sum()/atoms2_.rows();
- //center the structures
- atoms1_=this->SubtractVecFromMatrixRows(atoms1_, avg1_);
- atoms2_=this->SubtractVecFromMatrixRows(atoms2_, avg2_);
- EMatX atoms2=atoms2_.transpose();
+ // SVD only determines the rotational component of the superposition
+ // we need to manually shift the centers of the two point sets on onto
+ // origin
- //single value decomposition
- Eigen::SVD<EMat3> svd(atoms2*atoms1_);
+ EMatX atoms1 = MatrixShiftedBy(atoms1_, avg1);
+ EMatX atoms2 = MatrixShiftedBy(atoms2_, avg2).transpose();
+ // determine rotational component
+ Eigen::SVD<EMat3> svd(atoms2*atoms1);
EMatX matrixVT=svd.matrixV().transpose();
+
//determine rotation
Real detv=matrixVT.determinant();
Real dett=svd.matrixU().determinant();
Real det=detv*dett;
- EMat3 ERot;
+ EMat3 rotation;
if (det<0) {
EMat3 tmat=EMat3::Identity();
tmat(2,2)=-1;
- ERot=(svd.matrixU()*tmat)*matrixVT;
- }else{
- ERot=svd.matrixU()*matrixVT;
+ rotation = (svd.matrixU()*tmat)*matrixVT;
+ } else {
+ rotation = svd.matrixU()*matrixVT;
}
-
+
SuperpositionResult res;
- // prepare rmsd calculation
- geom::Vec3 shift=SuperposerSVDImpl::EigenVec3ToVec3(avg2_);
- geom::Vec3 com_vec=-SuperposerSVDImpl::EigenVec3ToVec3(avg1_);
- geom::Mat3 rot=SuperposerSVDImpl::EigenMat3ToMat3(ERot.transpose());
+
+ geom::Vec3 shift = EigenVec3ToVec3(avg2);
+ geom::Vec3 com_vec = -EigenVec3ToVec3(avg1);
+ geom::Mat3 rot = EigenMat3ToMat3(rotation.transpose());
geom::Mat4 mat4_com, mat4_rot, mat4_shift;
mat4_rot.PasteRotation(rot);
mat4_shift.PasteTranslation(shift);
mat4_com.PasteTranslation(com_vec);
//save whole transformation
- res.transformation=geom::Mat4(mat4_shift*mat4_rot*mat4_com);
+ res.transformation = geom::Mat4(mat4_shift*mat4_rot*mat4_com);
return res;
-
}
-
-SuperpositionResult SuperposerSVD::Run(const mol::EntityView& ev1,
- const mol::EntityView& ev2)
-{
- impl_->EntToMatrices(ev1, ev2);
- SuperpositionResult r=impl_->DoSuperposition();
- r.rmsd=CalculateRMSD(ev1, ev2, r.transformation);
- return r;
+MeanSquareMinimizer& MeanSquareMinimizer::operator=(const MeanSquareMinimizer& rhs) {
+ MeanSquareMinimizer tmp(rhs);
+ this->swap(tmp);
+ return *this;
}
-SuperpositionResult SuperposerSVD::Run(const mol::AtomViewList& atoms1,
- const mol::AtomViewList& atoms2)
-{
- return impl_->Run(atoms1, atoms2);
+MeanSquareMinimizer::MeanSquareMinimizer(const MeanSquareMinimizer& rhs) {
+ if (rhs.impl_)
+ impl_ = new MeanSquareMinimizerImpl(*rhs.impl_);
+ else
+ impl_ = NULL;
}
-SuperpositionResult SuperposerSVD::Run(const std::vector<geom::Vec3>& pl1,
- const std::vector<geom::Vec3>& pl2)
-{
- return impl_->Run(pl1,pl2);
-}
-
-SuperpositionResult SuperposerSVDImpl::Run(const mol::AtomViewList& atoms1,
- const mol::AtomViewList& atoms2)
-{
- mol::AtomViewList::const_iterator a_ev1=atoms1.begin();
- mol::AtomViewList::const_iterator a_ev1_end=atoms1.end();
- mol::AtomViewList::const_iterator a_ev2=atoms2.begin();
-
- for (int counter=0; a_ev1_end!=a_ev1; ++a_ev1, ++a_ev2, ++counter) {
- mol::AtomView av=*a_ev1;
- geom::Vec3 pos=av.GetPos();
- ERVec3 e=SuperposerSVDImpl::Vec3ToEigenVec(pos);
- atoms1_.row(counter)=e ;
- av=*a_ev2;
- pos=av.GetPos();
- e=SuperposerSVDImpl::Vec3ToEigenVec(pos);
- atoms2_.row(counter)=e;
- }
- return this->DoSuperposition();
-}
-SuperpositionResult SuperposerSVDImpl::Run(const std::vector<geom::Vec3>& pl1,
- const std::vector<geom::Vec3>& pl2)
-{
- for (unsigned int counter=0; counter<pl1.size();++counter) {
- atoms1_.row(counter)=ERVec3(SuperposerSVDImpl::Vec3ToEigenVec(pl1[counter]));
- atoms2_.row(counter)=ERVec3(SuperposerSVDImpl::Vec3ToEigenVec(pl2[counter]));
- }
- return this->DoSuperposition();
+MeanSquareMinimizer::~MeanSquareMinimizer() {
+ if (impl_) delete impl_;
}
+SuperpositionResult MeanSquareMinimizer::MinimizeOnce() const {
+ return impl_->MinimizeOnce();
+}
SuperpositionResult SuperposeAtoms(const mol::AtomViewList& atoms1,
const mol::AtomViewList& atoms2,
bool apply_transform=true)
{
- int counta=static_cast<int>(atoms1.size());
- int countb=static_cast<int>(atoms2.size());
- if (counta!=countb){
- throw Error("atom counts do not match");
- }
- if ((counta<3)){
- throw Error("at least 3 atoms are required");
- }
- SuperposerSVD sp(counta, true);
-
- SuperpositionResult res=sp.Run(atoms1, atoms2);
- //save rmsd info
- res.rmsd=calc_rmsd_for_atom_lists(atoms1, atoms2, res.transformation);
- res.ncycles=1;
- mol::AtomView jv=atoms1.front();
- if (apply_transform){
- mol::XCSEditor ed=jv.GetResidue().GetChain().GetEntity().GetHandle().EditXCS();
- ed.ApplyTransform(res.transformation);
+ MeanSquareMinimizer msm = MeanSquareMinimizer::FromAtomLists(atoms1, atoms2);
+ SuperpositionResult result = msm.MinimizeOnce();
+
+ result.ncycles=1;
+ result.rmsd = calc_rmsd_for_atom_lists(atoms1, atoms2, result.transformation);
+ if (apply_transform) {
+ mol::AtomView jv=atoms1.front();
+ mol::XCSEditor ed=jv.GetEntity().GetHandle().EditXCS();
+ ed.ApplyTransform(result.transformation);
}
- return res;
+ return result;
+
}
SuperpositionResult SuperposeSVD(const mol::EntityView& ev1,
const mol::EntityView& ev2,
bool apply_transform=true) {
- int counta=ev1.GetAtomCount();
- int countb=ev2.GetAtomCount();
- if (counta!=countb){
- throw Error("atom counts of the two views are not equal");
- }
- if ((counta<3)){
- throw Error("atom counts of any view must be larger or "
- "equal 3");
- }
- int nrows=counta;
- SuperposerSVD sp(nrows, true);
- SuperpositionResult res=sp.Run(ev1, ev2);
- res.entity_view1=ev1;
- res.entity_view2=ev2;
- res.ncycles=1;
- //save rmsd info
- res.rmsd=CalculateRMSD(ev1, ev2, res.transformation);
- if (apply_transform){
- mol::XCSEditor ed=ev1.GetHandle().EditXCS();
- ed.ApplyTransform(res.transformation);
- }
- return res;
+ AtomViewList atoms1 = ev1.GetAtomList();
+ AtomViewList atoms2 = ev2.GetAtomList();
+
+ SuperpositionResult result = SuperposeAtoms(atoms1, atoms2, apply_transform);
+ result.entity_view1 = ev1;
+ result.entity_view2 = ev2;
+ return result;
}
SuperpositionResult SuperposeSVD(const std::vector<geom::Vec3>& pl1,
const std::vector<geom::Vec3>& pl2)
{
- if (pl1.size()!=pl2.size()){
- throw Error("pointlist lengths not equal");
- }
- if ((pl1.size()<3)){
- throw Error("three or more points required");
- }
- SuperposerSVD sp(pl1.size(), true);
- SuperpositionResult res=sp.Run(pl1, pl2);
- res.entity_view1=EntityView();
- res.entity_view2=EntityView();
- res.ncycles=1;
- res.rmsd=0.0;
- return res;
-}
-
-SuperpositionResult IterativeSuperposition(mol::EntityView& ev1,
- mol::EntityView& ev2,
- int ncycles,
- Real dist_thres,
- bool apply_transform=true) {
- SuperpositionResult result;
- int counter=0;
- mol::AtomViewList atoms_a, atoms_b;
- std::copy(ev1.AtomsBegin(), ev1.AtomsEnd(), std::back_inserter(atoms_a));
- std::copy(ev2.AtomsBegin(), ev2.AtomsEnd(), std::back_inserter(atoms_b));
- bool converged=false;
-
- if (dist_thres<=0){
- dist_thres=5.0;
- }
- counter=0;
- while (!converged) {
- if (atoms_a.size()<3) {
- return result;
- }
- result=SuperposeAtoms(atoms_a, atoms_b, false);
- converged=ncycles<=counter;
- if (converged) {
- continue;
- }
- converged=true;
- for (mol::AtomViewList::iterator i=atoms_a.begin(), j=atoms_b.begin(),
- e=atoms_a.end(); i!=e; ++i, ++j) {
- mol::AtomView iv=*i;
- mol::AtomView jv=*j;
-
- if (geom::Length(geom::Vec3(result.transformation*geom::Vec4(iv.GetPos()))-jv.GetPos())>dist_thres) {
- *i=*j=mol::AtomView();
- converged=false;
- }
- }
- if (!converged) {
- mol::AtomViewList::iterator end;
- mol::AtomView invalid_atom;
- end=std::remove(atoms_a.begin(), atoms_a.end(), invalid_atom);
- atoms_a.erase(end, atoms_a.end());
- end=std::remove(atoms_b.begin(), atoms_b.end(), invalid_atom);
- atoms_b.erase(end, atoms_b.end());
- }
- counter++;
-
- }
- result.ncycles=counter;
- result.entity_view1=CreateViewFromAtomList(atoms_a);
- result.entity_view2=CreateViewFromAtomList(atoms_b);
-
- if (apply_transform){
- mol::XCSEditor ed=ev1.GetHandle().EditXCS();
- ed.ApplyTransform(result.transformation);
- }
+ MeanSquareMinimizer msm = MeanSquareMinimizer::FromPointLists(pl1, pl2);
+ SuperpositionResult result = msm.MinimizeOnce();
+
+ result.ncycles=1;
+ result.rmsd = calc_rmsd_for_point_lists(pl1, pl2, result.transformation);
return result;
}
diff --git a/modules/mol/alg/src/svd_superpose.hh b/modules/mol/alg/src/svd_superpose.hh
index dad42504c..8839a8ef5 100644
--- a/modules/mol/alg/src/svd_superpose.hh
+++ b/modules/mol/alg/src/svd_superpose.hh
@@ -27,6 +27,7 @@
#include <ost/geom/geom.hh>
#include <ost/mol/entity_view.hh>
#include <ost/mol/alg/module_config.hh>
+#include <ost/mol/atom_view.hh>
namespace ost { namespace mol {
@@ -49,25 +50,32 @@ struct DLLEXPORT_OST_MOL_ALG SuperpositionResult {
mol::EntityView entity_view2;
};
-class SuperposerSVDImpl;
-/// \brief efficiently superpose a bunch of models with the same number of atoms
-/// Choose either two EntityViews or two AtomViewLists.
-class DLLEXPORT_OST_MOL_ALG SuperposerSVD {
+// the mean square minimizer is split into a private and a public class to avoid
+// exposing the lib Eigen classes to the outside world
+class MeanSquareMinimizerImpl;
+
+
+class DLLEXPORT_OST_MOL_ALG MeanSquareMinimizer {
public:
- SuperposerSVD(int natoms, bool alloc_atoms);
- ~SuperposerSVD();
- SuperpositionResult Run(const mol::EntityView& ev1,
- const mol::EntityView& ev2);
+ static MeanSquareMinimizer FromAtomLists(const mol::AtomViewList& atoms,
+ const mol::AtomViewList& atoms_ref);
+ static MeanSquareMinimizer FromPointLists(const std::vector<geom::Vec3>& points,
+ const std::vector<geom::Vec3>& points_ref);
+ SuperpositionResult MinimizeOnce() const;
- SuperpositionResult Run(const mol::AtomViewList& atoms1,
- const mol::AtomViewList& atoms2);
+ ~MeanSquareMinimizer();
+
- SuperpositionResult Run(const std::vector<geom::Vec3>& pl1,
- const std::vector<geom::Vec3>& pl2);
+ MeanSquareMinimizer& operator=(const MeanSquareMinimizer& rhs);
+ MeanSquareMinimizer(const MeanSquareMinimizer& rhs);
+ void swap(MeanSquareMinimizer& rhs) {
+ std::swap(rhs.impl_, impl_);
+ }
-private:
- SuperposerSVDImpl* impl_;
+protected:
+ MeanSquareMinimizer(): impl_(NULL) {}
+ MeanSquareMinimizerImpl* impl_;
};
/// \brief takes the corresponding atoms and superposes them
@@ -96,9 +104,6 @@ Real DLLEXPORT_OST_MOL_ALG CalculateRMSD(const mol::EntityView& ev1,
const geom::Mat4& transformation);
-/// \example superpose.py
-///
-/// Superpose 1AKE and 4AKE by a least squares fitting algorithm
-/// \sa \ref superpose.py "Superposing entities"
+
}}}//ns
#endif
diff --git a/modules/mol/alg/tests/test_superposition.cc b/modules/mol/alg/tests/test_superposition.cc
index 83bc902b3..6ef76a0d2 100644
--- a/modules/mol/alg/tests/test_superposition.cc
+++ b/modules/mol/alg/tests/test_superposition.cc
@@ -87,15 +87,11 @@ BOOST_AUTO_TEST_CASE(superposition_svd)
XCSEditor ed=f1.e.EditXCS();
EntityView ev2 = f2.e.CreateFullView();
ChainHandle ch1=f1.e.GetChainList()[0];
- ResidueHandleIter rit=ch1.ResiduesBegin();
- for (;rit!=ch1.ResiduesEnd();++rit){
- ResidueHandle resh=*rit;
- AtomHandleIter at=resh.AtomsBegin();
- for (;at!=resh.AtomsEnd();++at){
- AtomHandle atom=*at;
+ AtomHandleList atoms = f1.e.GetAtomList();
+ for (AtomHandleList::const_iterator i = atoms.begin(), e = atoms.end(); i !=e ; ++i) {
+ AtomHandle atom=*i;
ed.SetAtomPos(atom, geom::EulerTransformation(0.1, 0, 0)*atom.GetPos()
+geom::Vec3(0, 5, 0));
- }
}
SuperpositionResult res;
res=SuperposeSVD(ev1, ev2, true);
diff --git a/modules/mol/base/src/CMakeLists.txt b/modules/mol/base/src/CMakeLists.txt
index 8d206d3c6..86d5ccca6 100644
--- a/modules/mol/base/src/CMakeLists.txt
+++ b/modules/mol/base/src/CMakeLists.txt
@@ -19,7 +19,6 @@ entity_property_mapper.cc
entity_view.cc
entity_visitor.cc
ics_editor.cc
-iterator.cc
not_connected_error.cc
property_id.cc
query.cc
@@ -68,8 +67,6 @@ entity_visitor.hh
entity_visitor_fw.hh
handle_type_fw.hh
ics_editor.hh
-iterator.hh
-iterator_fw.hh
not_connected_error.hh
property_id.hh
query.hh
diff --git a/modules/mol/base/src/chain_handle.cc b/modules/mol/base/src/chain_handle.cc
index c1e745f2c..b15f431e0 100644
--- a/modules/mol/base/src/chain_handle.cc
+++ b/modules/mol/base/src/chain_handle.cc
@@ -21,8 +21,8 @@
#include <ost/mol/residue_handle.hh>
#include <ost/mol/impl/chain_impl.hh>
#include <ost/mol/impl/residue_impl.hh>
-#include <ost/mol/iterator.hh>
#include <ost/mol/impl/entity_impl.hh>
+
namespace ost { namespace mol {
@@ -96,7 +96,12 @@ AtomHandleList ChainHandle::GetAtomList() const
{
this->CheckValidity();
AtomHandleList atoms;
- std::copy(AtomsBegin(), AtomsEnd(), std::back_inserter(atoms));
+ atoms.reserve(this->GetAtomCount());
+ for (impl::ResidueImplList::const_iterator i= Impl()->GetResidueList().begin(),
+ e = Impl()->GetResidueList().end(); i !=e; ++i) {
+ std::copy((*i)->GetAtomList().begin(), (*i)->GetAtomList().end(),
+ std::back_inserter(atoms));
+ }
return atoms;
}
@@ -104,6 +109,7 @@ ResidueHandleList ChainHandle::GetResidueList() const
{
ResidueHandleList reslist;
this->CheckValidity();
+
const impl::ResidueImplList& rl = Impl()->GetResidueList();
reslist.reserve(rl.size());
for(impl::ResidueImplList::const_iterator it=rl.begin();it!=rl.end();++it) {
@@ -122,66 +128,6 @@ bool ChainHandle::operator!=(const ChainHandle& ref) const
return Impl()!=ref.Impl();
}
-ResidueHandleIter ChainHandle::ResiduesBegin() const {
- this->CheckValidity();
- impl::ChainImplPtr c=Impl();
- return ResidueHandleIter(c->GetEntity()->GetChainIter(this->GetName()),
- impl::begin(c->GetResidueList()),
- c->GetEntity(), true);
-}
-
-ResidueHandleIter ChainHandle::ResiduesEnd() const {
- this->CheckValidity();
- impl::ChainImplPtr c=Impl();
- impl::pointer_it<impl::ChainImplPtr> cc=c->GetEntity()->GetChainIter(this->GetName());
- /*impl::pointer_it<impl::ChainImplPtr> nc=cc; ++nc;
-
- if (nc!=impl::end(c->GetEntity()->GetChainList())) {
- return ResidueHandleIter(nc, (*nc)->GetResidueList().begin(),
- c->GetEntity(), true);
- } else {
- return ResidueHandleIter(cc, c->GetResidueList().end(),
- c->GetEntity(), true);
- }*/
- return ResidueHandleIter(cc, impl::end(c->GetResidueList()),
- c->GetEntity(), true);
-}
-
-AtomHandleIter ChainHandle::AtomsBegin() const
-{
- this->CheckValidity();
- impl::ChainImplPtr c=Impl();
- if (c->GetResidueList().empty()) {
-
- return AtomHandleIter();
- }
- impl::pointer_it<impl::ChainImplPtr> cc=c->GetEntity()->GetChainIter(this->GetName());
- return AtomHandleIter(cc, impl::begin(c->GetResidueList()),
- impl::begin(c->GetResidueList().front()->GetAtomList()),
- c->GetEntity(), true);
-
-}
-
-AtomHandleIter ChainHandle::AtomsEnd() const {
- this->CheckValidity();
- impl::ChainImplPtr c=Impl();
- if (c->GetResidueList().empty()) {
- return AtomHandleIter();
- }
- impl::pointer_it<impl::ChainImplPtr> cc=c->GetEntity()->GetChainIter(this->GetName());
- impl::pointer_it<impl::ChainImplPtr> nc=cc; ++nc;
- impl::ResidueImplList& rc=(*nc)->GetResidueList();
- if (nc!=impl::end(c->GetEntity()->GetChainList())) {
- return AtomHandleIter(nc, impl::begin(rc),
- impl::begin(rc.front()->GetAtomList()),
- c->GetEntity(), false);
- } else {
- return AtomHandleIter(cc, impl::end(c->GetResidueList()),
- impl::end(c->GetResidueList().back()->GetAtomList()),
- c->GetEntity(), false);
- }
-}
-
ResidueHandle ChainHandle::GetResidueByIndex(int index) const {
this->CheckValidity();
if (index<0 || !(static_cast<int>(Impl()->GetResidueList().size())>index))
diff --git a/modules/mol/base/src/chain_handle.hh b/modules/mol/base/src/chain_handle.hh
index 8f88c6e37..5ab00787c 100644
--- a/modules/mol/base/src/chain_handle.hh
+++ b/modules/mol/base/src/chain_handle.hh
@@ -32,7 +32,6 @@
#include "handle_type_fw.hh"
#include "residue_prop.hh"
#include "entity_visitor_fw.hh"
-#include "iterator_fw.hh"
#include "sec_structure.hh"
namespace ost { namespace mol {
@@ -134,28 +133,6 @@ public:
/// \c ChainHandle::GetResidueCount
ResidueHandleList GetResidueList() const;
- /// \brief Get iterator pointing to begin of residues
- ///
- /// \note Iterators are not fail-safe, meaning that the behaviour is
- /// undefined when using a ResidueHandleIter and residues are inserted
- /// or removed from the chain
- /// \code
- /// ResidueHandleIter res_it=chain.ResiduesBegin(),
- /// res_end=chain.ResiduesEnd();
- /// for (; res_it!=res_end; ++res_it) {
- /// ResidueHandle res=*res_it;
- /// // do something with residue
- /// }
- /// \endcode
- ResidueHandleIter ResiduesBegin() const;
-
- /// \brief Get iterator pointing to the end of the residues
- /// \sa ResiduesBegin
- ResidueHandleIter ResiduesEnd() const;
-
- AtomHandleIter AtomsBegin() const;
- AtomHandleIter AtomsEnd() const;
-
/// \brief Get entity's mass
Real GetMass() const;
diff --git a/modules/mol/base/src/coord_frame.cc b/modules/mol/base/src/coord_frame.cc
index 52907d7f7..8ba796305 100644
--- a/modules/mol/base/src/coord_frame.cc
+++ b/modules/mol/base/src/coord_frame.cc
@@ -420,8 +420,15 @@ namespace ost { namespace mol {
void GetPositions(const EntityView& sele, std::vector<geom::Vec3>& ref_pos){
ref_pos.reserve(sele.GetAtomCount());
- for (mol::AtomViewIter a=sele.AtomsBegin(),e=sele.AtomsEnd(); a!=e; ++a) {
- ref_pos.push_back((*a).GetPos());
+ for (ChainViewList::const_iterator ci = sele.GetChainList().begin(),
+ ce = sele.GetChainList().end(); ci != ce; ++ci) {
+ for (ResidueViewList::const_iterator ri = ci->GetResidueList().begin(),
+ re = ci->GetResidueList().end(); ri != re; ++ri) {
+ for (AtomViewList::const_iterator ai = ri->GetAtomList().begin(),
+ ae = ri->GetAtomList().end(); ai != ae; ++ai) {
+ ref_pos.push_back(ai->GetPos());
+ }
+ }
}
}
diff --git a/modules/mol/base/src/entity_handle.cc b/modules/mol/base/src/entity_handle.cc
index 36afa61b3..a463aa851 100644
--- a/modules/mol/base/src/entity_handle.cc
+++ b/modules/mol/base/src/entity_handle.cc
@@ -20,13 +20,15 @@
#include <ost/log.hh>
#include "impl/entity_impl.hh"
+
+#include "atom_view.hh"
#include "bond_handle.hh"
#include "torsion_handle.hh"
#include "entity_visitor.hh"
#include "atom_handle.hh"
+#include "residue_handle.hh"
#include "chain_handle.hh"
#include "entity_handle.hh"
-#include "iterator.hh"
#include "impl/chain_impl.hh"
#include "impl/atom_impl.hh"
#include "impl/residue_impl.hh"
@@ -315,68 +317,6 @@ AtomHandle EntityHandle::FindAtom(const String& chain_name,
}
-ResidueHandleIter EntityHandle::ResiduesBegin() const {
- this->CheckValidity();
- if (Impl()->GetChainList().empty()) {
- return ResidueHandleIter();
- }
- impl::EntityImplPtr i=Impl();
- impl::ChainImplPtr chain=i->GetChainList().front();
- return ResidueHandleIter(impl::begin(i->GetChainList()),
- impl::begin(chain->GetResidueList()), i, true);
-}
-
-ResidueHandleIter EntityHandle::ResiduesEnd() const {
- this->CheckValidity();
- if (Impl()->GetChainList().empty()) {
- return ResidueHandleIter();
- }
- impl::EntityImplPtr i=Impl();
- impl::ChainImplPtr chain=i->GetChainList().back();
- return ResidueHandleIter(impl::end(i->GetChainList()),
- impl::end(chain->GetResidueList()), i, false);
-}
-
-ChainHandleIter EntityHandle::ChainsBegin() const {
- this->CheckValidity();
- return ChainHandleIter(Impl()->GetChainList().begin());
-}
-
-
-ChainHandleIter EntityHandle::ChainsEnd() const {
- return ChainHandleIter(Impl()->GetChainList().end());
-}
-
-
-AtomHandleIter EntityHandle::AtomsBegin() const {
- this->CheckValidity();
- impl::EntityImplPtr ent=Impl();
- if (ent->GetChainList().empty()) {
- return AtomHandleIter();
- }
- impl::ResidueImplList& r=ent->GetChainList().front()->GetResidueList();
- if (r.empty()) {
- return AtomHandleIter();
- }
- return AtomHandleIter(impl::begin(ent->GetChainList()), impl::begin(r),
- impl::begin(r.front()->GetAtomList()), ent, true);
-}
-
-AtomHandleIter EntityHandle::AtomsEnd() const
-{
- this->CheckValidity();
- impl::EntityImplPtr ent=Impl();
- if (ent->GetChainList().empty()) {
- return AtomHandleIter();
- }
- impl::ResidueImplList& r=ent->GetChainList().back()->GetResidueList();
- if (r.empty()) {
- return AtomHandleIter();
- }
- return AtomHandleIter(impl::end(ent->GetChainList()), impl::end(r),
- impl::end(r.back()->GetAtomList()), ent, false);
-}
-
XCSEditor EntityHandle::EditXCS(EditMode mode) const
{
this->CheckValidity();
@@ -395,7 +335,14 @@ ResidueHandleList EntityHandle::GetResidueList() const
{
this->CheckValidity();
ResidueHandleList residues;
- std::copy(ResiduesBegin(), ResiduesEnd(), std::back_inserter(residues));
+ residues.reserve(this->GetResidueCount());
+ for (impl::ChainImplList::const_iterator i = Impl()->GetChainList().begin(),
+ e = Impl()->GetChainList().end(); i != e; ++i) {
+ for (impl::ResidueImplList::const_iterator j = (*i)->GetResidueList().begin(),
+ e2 = (*i)->GetResidueList().end(); j != e2; ++j) {
+ residues.push_back(ResidueHandle(*j));
+ }
+ }
return residues;
}
@@ -404,7 +351,17 @@ AtomHandleList EntityHandle::GetAtomList() const
{
this->CheckValidity();
AtomHandleList atoms;
- std::copy(AtomsBegin(), AtomsEnd(), std::back_inserter(atoms));
+ atoms.reserve(this->GetAtomCount());
+ for (impl::ChainImplList::const_iterator i = Impl()->GetChainList().begin(),
+ e = Impl()->GetChainList().end(); i != e; ++i) {
+ for (impl::ResidueImplList::const_iterator j = (*i)->GetResidueList().begin(),
+ e2 = (*i)->GetResidueList().end(); j != e2; ++j) {
+ for (impl::AtomImplList::const_iterator k = (*j)->GetAtomList().begin(),
+ e3 = (*j)->GetAtomList().end(); k != e3; ++k) {
+ atoms.push_back(AtomHandle(*k));
+ }
+ }
+ }
return atoms;
}
diff --git a/modules/mol/base/src/entity_handle.hh b/modules/mol/base/src/entity_handle.hh
index d165a997d..5cb26134b 100644
--- a/modules/mol/base/src/entity_handle.hh
+++ b/modules/mol/base/src/entity_handle.hh
@@ -29,7 +29,6 @@
#include "entity_view.hh"
#include "chain_handle.hh"
#include "handle_type_fw.hh"
-#include "iterator_fw.hh"
#include "editor_type_fw.hh"
@@ -125,38 +124,6 @@ public:
/// \brief Get list of chains
ChainHandleList GetChainList() const;
- /// \brief Iterator pointing to the beginning of the residues.
- ///
- /// The residues are first ordered by chain and then ascending residue number.
- /// \note The residue range is not fail-safe. If residues are removed or
- /// the iterators become invalid and the behaviour undefined.
- ResidueHandleIter ResiduesBegin() const;
-
- /// \brief Iterator pointing to the end of th residues
- ///
- /// The two iterators returned by ResiduesBegin() and ResiduesEnd() form a
- /// half-closed range. It is cheaper to cache the iterator returned by
- /// ResiduesEnd() than to call ResiduesEnd() after every loop, i.e. like
- /// \code
- /// // e is an instance of EntityHandle
- /// for (ResidueHandleIter i=e.ResiduesBegin(), x=e.ResiduesEnd(); i!=x; ++i) {
- /// . . .
- /// }
- /// \endcode
- ResidueHandleIter ResiduesEnd() const;
-
- /// \brief Iterator pointing to the beginning of the chains
- ChainHandleIter ChainsBegin() const;
-
- /// \brief Iterator pointing to the end of the chains
- ChainHandleIter ChainsEnd() const;
-
- /// \brief Iterator pointing to the beginning of atoms
- AtomHandleIter AtomsBegin() const;
-
- /// \brief Iterator pointing to the end of atoms
- AtomHandleIter AtomsEnd() const;
-
///\name Single item addressing
//@{
/// \brief Find chain by name.
diff --git a/modules/mol/base/src/entity_view.cc b/modules/mol/base/src/entity_view.cc
index 0c2d6db41..bb585270e 100644
--- a/modules/mol/base/src/entity_view.cc
+++ b/modules/mol/base/src/entity_view.cc
@@ -34,7 +34,6 @@
#include "residue_view.hh"
#include "atom_view.hh"
#include "bond_table.hh"
-#include "iterator.hh"
#include "query_state.hh"
#include "entity_property_mapper.hh"
@@ -161,29 +160,44 @@ int EntityView::GetResidueCount() const
geom::Vec3 EntityView::GetCenterOfAtoms() const
{
geom::Vec3 center;
- if (this->HasAtoms()) {
- unsigned int counter=0;
- AtomViewIter it=this->AtomsBegin();
- for(; it!=this->AtomsEnd(); ++it, ++counter) {
- center+=(*it).GetPos();
+
+ if (!this->HasAtoms()) {
+ return center;
+ }
+ unsigned int counter = 0;
+ for (ChainViewList::const_iterator ci = this->GetChainList().begin(),
+ ce = this->GetChainList().end(); ci != ce; ++ci) {
+ for (ResidueViewList::const_iterator ri = ci->GetResidueList().begin(),
+ re = ci->GetResidueList().end(); ri != re; ++ri) {
+ for (AtomViewList::const_iterator ai = ri->GetAtomList().begin(),
+ ae = ri->GetAtomList().end(); ai != ae; ++ai) {
+ center += ai->GetPos();
+ counter += 1;
+ }
}
- center/=static_cast<Real>(counter);
}
- return center;
+ return center/static_cast<Real>(counter);
}
geom::Vec3 EntityView::GetCenterOfMass() const
{
geom::Vec3 center;
Real mass = this->GetMass();
- if (this->HasAtoms() && mass>0) {
- AtomViewIter it=this->AtomsBegin();
- for(; it!=this->AtomsEnd(); ++it) {
- center+=(*it).GetPos()*(*it).GetMass();
+ if (mass==0) {
+ return center;
+ }
+
+ for (ChainViewList::const_iterator ci = this->GetChainList().begin(),
+ ce = this->GetChainList().end(); ci != ce; ++ci) {
+ for (ResidueViewList::const_iterator ri = ci->GetResidueList().begin(),
+ re = ci->GetResidueList().end(); ri != re; ++ri) {
+ for (AtomViewList::const_iterator ai = ri->GetAtomList().begin(),
+ ae = ri->GetAtomList().end(); ai != ae; ++ai) {
+ center += ai->GetPos()*ai->GetMass();
+ }
}
- center/=mass;
}
- return center;
+ return center / mass;
}
Real EntityView::GetMass() const
@@ -457,57 +471,6 @@ const ChainViewList& EntityView::GetChainList() const
return data_->chains;
}
-AtomViewIter EntityView::AtomsBegin() const
-{
- this->CheckValidity();
- if (data_->chains.empty()) {
- return AtomViewIter();
- }
- const ResidueViewList& rvl=data_->chains.front().GetResidueList();
- if (rvl.empty()) {
- return AtomViewIter();
- }
- return AtomViewIter(impl::begin(data_->chains), impl::begin(rvl),
- impl::begin(rvl.front().GetAtomList()),
- *this, true);
-}
-
-AtomViewIter EntityView::AtomsEnd() const {
- this->CheckValidity();
- if (data_->chains.empty()) {
- return AtomViewIter();
- }
- const ResidueViewList& rvl=data_->chains.back().GetResidueList();
- if (rvl.empty()) {
- return AtomViewIter();
- }
- return AtomViewIter(impl::end(data_->chains), impl::end(rvl),
- impl::end(rvl.back().GetAtomList()), *this, false);
-}
-
-ResidueViewIter EntityView::ResiduesBegin() const
-{
- this->CheckValidity();
- if (data_->chains.empty()) {
- return ResidueViewIter();
- }
- const ResidueViewList& rvl=data_->chains.front().GetResidueList();
- return ResidueViewIter(impl::begin(data_->chains),
- impl::begin(rvl), *this, true);
-}
-
-ResidueViewIter EntityView::ResiduesEnd() const
-{
- this->CheckValidity();
- if (data_->chains.empty()) {
- return ResidueViewIter();
- }
- const ResidueViewList& rvl=data_->chains.back().GetResidueList();
- return ResidueViewIter(impl::end(data_->chains),
- impl::end(rvl), *this, false);
-}
-
-
void EntityView::RemoveAtom(AtomView view) {
this->CheckValidity();
if (!view.IsValid())
@@ -606,7 +569,12 @@ ResidueViewList EntityView::GetResidueList() const
{
this->CheckValidity();
ResidueViewList residues;
- std::copy(ResiduesBegin(), ResiduesEnd(), std::back_inserter(residues));
+ residues.reserve(this->GetResidueCount());
+ ChainViewList::const_iterator i;
+ for (i=data_->chains.begin(); i!=data_->chains.end(); ++i) {
+ std::copy(i->GetResidueList().begin(), i->GetResidueList().end(),
+ std::back_inserter(residues));
+ }
return residues;
}
@@ -615,7 +583,15 @@ AtomViewList EntityView::GetAtomList() const
{
this->CheckValidity();
AtomViewList atoms;
- std::copy(AtomsBegin(), AtomsEnd(), std::back_inserter(atoms));
+ atoms.reserve(this->GetAtomCount());
+ ChainViewList::const_iterator i;
+ for (i=data_->chains.begin(); i!=data_->chains.end(); ++i) {
+ for (ResidueViewList::const_iterator j = (*i).GetResidueList().begin(),
+ e = (*i).GetResidueList().end(); j!=e; ++j) {
+ std::copy(j->GetAtomList().begin(), j->GetAtomList().end(),
+ std::back_inserter(atoms));
+ }
+ }
return atoms;
}
@@ -687,16 +663,25 @@ std::pair<Real,Real> EntityView::GetMinMax(const String& prop,
Prop::Level hint) const
{
EntityPropertyMapper epm(prop, hint);
+
Real min_v=std::numeric_limits<Real>::max();
Real max_v=-std::numeric_limits<Real>::max();
- for(AtomViewIter it=AtomsBegin(); it!=this->AtomsEnd(); ++it) {
- try {
- Real v=epm.Get(*it);
- max_v=std::max(v, max_v);
- min_v=std::min(v, min_v);
- } catch(...) {
- // do nothing in case of missing property
- continue;
+
+ for (ChainViewList::const_iterator ci = this->GetChainList().begin(),
+ ce = this->GetChainList().end(); ci != ce; ++ci) {
+ for (ResidueViewList::const_iterator ri = ci->GetResidueList().begin(),
+ re = ci->GetResidueList().end(); ri != re; ++ri) {
+ for (AtomViewList::const_iterator ai = ri->GetAtomList().begin(),
+ ae = ri->GetAtomList().end(); ai != ae; ++ai) {
+ try {
+ Real v=epm.Get(*ai);
+ max_v=std::max(v, max_v);
+ min_v=std::min(v, min_v);
+ } catch(...) {
+ // do nothing in case of missing property
+ continue;
+ }
+ }
}
}
return std::make_pair(min_v,max_v);
@@ -755,15 +740,23 @@ geom::AlignedCuboid EntityView::GetBounds() const
this->CheckValidity();
geom::Vec3 mmin( std::numeric_limits<Real>::max());
geom::Vec3 mmax(-std::numeric_limits<Real>::max());
- if (this->HasAtoms()) {
- for(AtomViewIter it=AtomsBegin(); it!=this->AtomsEnd(); ++it) {
- mmax=geom::Max(mmax, (*it).GetPos());
- mmin=geom::Min(mmin, (*it).GetPos());
- }
+
+ if (!this->HasAtoms()) {
return geom::AlignedCuboid(mmin, mmax);
- } else {
- return geom::AlignedCuboid(geom::Vec3(), geom::Vec3());
}
+
+ for (ChainViewList::const_iterator ci = this->GetChainList().begin(),
+ ce = this->GetChainList().end(); ci != ce; ++ci) {
+ for (ResidueViewList::const_iterator ri = ci->GetResidueList().begin(),
+ re = ci->GetResidueList().end(); ri != re; ++ri) {
+ for (AtomViewList::const_iterator ai = ri->GetAtomList().begin(),
+ ae = ri->GetAtomList().end(); ai != ae; ++ai) {
+ mmax=geom::Max(mmax, ai->GetPos());
+ mmin=geom::Min(mmin, ai->GetPos());
+ }
+ }
+ }
+ return geom::AlignedCuboid(mmin, mmax);
}
diff --git a/modules/mol/base/src/entity_view.hh b/modules/mol/base/src/entity_view.hh
index 458e6d40f..e88845d8b 100644
--- a/modules/mol/base/src/entity_view.hh
+++ b/modules/mol/base/src/entity_view.hh
@@ -30,7 +30,6 @@
#include <ost/mol/entity_visitor_fw.hh>
#include <ost/mol/handle_type_fw.hh>
#include <ost/mol/query.hh>
-#include <ost/mol/iterator_fw.hh>
#include <ost/mol/property_id.hh>
#include <ost/geom/geom.hh>
@@ -281,20 +280,6 @@ public:
const ChainViewList& GetChainList() const;
- /// \brief Get iterator pointing to the beginning of atoms
- AtomViewIter AtomsBegin() const;
-
- /// \brief Get iterator pointing to the end of atoms
- AtomViewIter AtomsEnd() const;
-
- /// \brief Get iterator pointing to the beginnin of residues
- /// \sa ResiduesEnd
- ResidueViewIter ResiduesBegin() const;
-
- /// \brief get iterator pointing to the end of residues
- /// \sa ResiduesBegin
- ResidueViewIter ResiduesEnd() const;
-
/// \brief Get list of all residues included in the view
ResidueViewList GetResidueList() const;
diff --git a/modules/mol/base/src/iterator.cc b/modules/mol/base/src/iterator.cc
deleted file mode 100644
index f72aefe62..000000000
--- a/modules/mol/base/src/iterator.cc
+++ /dev/null
@@ -1,198 +0,0 @@
-//------------------------------------------------------------------------------
-// This file is part of the OpenStructure project <www.openstructure.org>
-//
-// Copyright (C) 2008-2011 by the OpenStructure authors
-//
-// This library is free software; you can redistribute it and/or modify it under
-// the terms of the GNU Lesser General Public License as published by the Free
-// Software Foundation; either version 3.0 of the License, or (at your option)
-// any later version.
-// This library is distributed in the hope that it will be useful, but WITHOUT
-// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
-// details.
-//
-// You should have received a copy of the GNU Lesser General Public License
-// along with this library; if not, write to the Free Software Foundation, Inc.,
-// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
-//------------------------------------------------------------------------------
-#include <ost/mol/impl/residue_impl.hh>
-#include <ost/mol/impl/chain_impl.hh>
-#include <ost/mol/impl/entity_impl.hh>
-#include <ost/mol/impl/atom_impl.hh>
-
-#include "iterator.hh"
-
-namespace ost { namespace mol {
-
-void ResidueHandleIter::SkipEmpty()
-{
- if (cur_res_==impl::end((*cur_chain_)->GetResidueList())) {
- // we have to skip over empty chains otherwise we end up pointing to an
- // invalid residue.
- do {
- ++cur_chain_;
- if (impl::end(ent_->GetChainList())==cur_chain_) {
- break;
- }
- cur_res_=impl::begin((*cur_chain_)->GetResidueList());
- } while ((*cur_chain_)->GetResidueList().empty());
- }
-}
-
-ResidueHandleIter& ResidueHandleIter::operator++()
-{
- ++cur_res_;
- this->SkipEmpty();
- return *this;
-}
-
-ResidueHandleIter::ResidueHandleIter(impl::pointer_it<impl::ChainImplPtr> chain_it,
- impl::pointer_it<impl::ResidueImplPtr> res_it,
- impl::EntityImplPtr ent, bool skip_empty)
- : cur_chain_(chain_it), cur_res_(res_it),
- ent_(ent) {
- if (skip_empty) {
- this->SkipEmpty();
- }
-}
-
-ResidueHandle ResidueHandleIter::operator*() {
- return ResidueHandle(*cur_res_);
-}
-
-
-
-ResidueView ResidueViewIter::operator*() {
- return ResidueView(*cur_res_);
-}
-
-void ResidueViewIter::SkipEmpty()
-{
- if (cur_res_==impl::end(cur_chain_->GetResidueList())) {
- // we have to skip over empty chains otherwise we end up pointing to an
- // invalid residue.
- do {
- ++cur_chain_;
- if (impl::end(ent_.GetChainList())==cur_chain_) {
- break;
- }
- cur_res_=impl::begin(cur_chain_->GetResidueList());
- } while (cur_chain_->GetResidueList().empty());
- }
-}
-ResidueViewIter& ResidueViewIter::operator++() {
- ++cur_res_;
- this->SkipEmpty();
- return *this;
-}
-
-ResidueViewIter::ResidueViewIter(impl::pointer_it<ChainView> chain_it,
- impl::pointer_it<ResidueView> res_it,
- EntityView ent, bool skip_empty)
- : cur_chain_(chain_it), cur_res_(res_it),
- ent_(ent) {
- if (skip_empty) {
- this->SkipEmpty();
- }
-}
-
-AtomHandleIter::AtomHandleIter(impl::pointer_it<impl::ChainImplPtr> chain_it,
- impl::pointer_it<impl::ResidueImplPtr> res_it,
- impl::pointer_it<impl::AtomImplPtr> atom_it,
- impl::EntityImplPtr ent, bool skip_empty)
- : cur_chain_(chain_it), cur_res_(res_it), cur_atom_(atom_it),
- ent_(ent)
-{
- if (skip_empty) {
- this->SkipEmpty();
- }
-}
-
-AtomHandle AtomHandleIter::operator*() {
- return AtomHandle(*cur_atom_);
-}
-
-void AtomHandleIter::SkipEmpty()
-{
- if (impl::end((*cur_res_)->GetAtomList())==cur_atom_) {
- // we have to skip over empty chains and residues otherwise we end up
- // pointing to an invalid atom.
- do {
- ++cur_res_;
- if (impl::end((*cur_chain_)->GetResidueList())==cur_res_) {
- do {
- ++cur_chain_;
- if (impl::end(ent_->GetChainList())==cur_chain_) {
- return;
- }
-
- cur_res_=impl::begin((*cur_chain_)->GetResidueList());
- if (!(*cur_chain_)->GetResidueList().empty())
- cur_atom_=impl::begin((*cur_res_)->GetAtomList());
- else
- cur_atom_=impl::pointer_it<impl::AtomImplPtr>(NULL);
- } while ((*cur_chain_)->GetResidueList().empty());
- } else {
- cur_atom_=impl::begin((*cur_res_)->GetAtomList());
- }
- } while ((*cur_res_)->GetAtomList().empty());
- }
-}
-
-AtomHandleIter& AtomHandleIter::operator++()
-{
- ++cur_atom_;
- this->SkipEmpty();
- return *this;
-}
-
-AtomViewIter::AtomViewIter(impl::pointer_it<ChainView> chain_it,
- impl::pointer_it<ResidueView> res_it,
- impl::pointer_it<AtomView> atom_it,
- EntityView ent, bool skip_empty)
- : cur_chain_(chain_it), cur_res_(res_it), cur_atom_(atom_it),
- ent_(ent) {
- if (skip_empty) {
- this->SkipEmpty();
- }
-}
-
-void AtomViewIter::SkipEmpty()
-{
- if (impl::end(cur_res_->GetAtomList())==cur_atom_) {
- do {
- ++cur_res_;
- if (impl::end(cur_chain_->GetResidueList())==cur_res_) {
- do {
- ++cur_chain_;
- if (impl::end(ent_.GetChainList())==cur_chain_) {
- return;
- }
- cur_res_=impl::begin(cur_chain_->GetResidueList());
- if (!cur_chain_->GetResidueList().empty()) {
- cur_atom_=impl::begin(cur_res_->GetAtomList());
- } else {
- cur_atom_=impl::pointer_it<AtomView>(NULL);
- }
- } while (cur_chain_->GetResidueList().empty());
- } else {
- cur_atom_=impl::begin(cur_res_->GetAtomList());
- }
- } while (cur_res_->GetAtomList().empty());
- }
-}
-
-AtomView AtomViewIter::operator*()
-{
- return AtomView(*cur_atom_);
-}
-
-AtomViewIter& AtomViewIter::operator++()
-{
- ++cur_atom_;
- this->SkipEmpty();
- return *this;
-}
-
-}} // ns
diff --git a/modules/mol/base/src/iterator.hh b/modules/mol/base/src/iterator.hh
deleted file mode 100644
index 868f0b180..000000000
--- a/modules/mol/base/src/iterator.hh
+++ /dev/null
@@ -1,237 +0,0 @@
-//------------------------------------------------------------------------------
-// This file is part of the OpenStructure project <www.openstructure.org>
-//
-// Copyright (C) 2008-2011 by the OpenStructure authors
-//
-// This library is free software; you can redistribute it and/or modify it under
-// the terms of the GNU Lesser General Public License as published by the Free
-// Software Foundation; either version 3.0 of the License, or (at your option)
-// any later version.
-// This library is distributed in the hope that it will be useful, but WITHOUT
-// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
-// details.
-//
-// You should have received a copy of the GNU Lesser General Public License
-// along with this library; if not, write to the Free Software Foundation, Inc.,
-// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
-//------------------------------------------------------------------------------
-#ifndef OST_ITERATOR_HH
-#define OST_ITERATOR_HH
-
-#include <map>
-
-
-#include <ost/mol/chain_handle.hh>
-#include <ost/mol/chain_view.hh>
-#include <ost/mol/residue_handle.hh>
-#include <ost/mol/residue_view.hh>
-#include <ost/mol/atom_handle.hh>
-#include <ost/mol/atom_view.hh>
-#include <ost/mol/entity_handle.hh>
-#include <ost/mol/entity_view.hh>
-#include <ost/mol/impl/entity_impl_fw.hh>
-#include <ost/mol/impl/chain_impl_fw.hh>
-#include <ost/mol/impl/residue_impl_fw.hh>
-#include <ost/mol/impl/pointer_iterator.hh>
-
-namespace ost { namespace mol {
-
-/// \brief adaptor to turn iterator into an iterator range
-///
-/// To turn two iterators marking the half-closed range [beg, end) into a
-/// iterator range use:
-/// \code
-/// // ent is an instance of entity
-/// IterRange<AtomHandleIter> atom_range(ent.AtomsBegin(), ent.AtomsEnd());
-/// for (; !atom_range.AtEnd(); ++atom_range) {
-/// }
-/// \endcode
-template <typename I>
-class IterRange {
-public:
- typedef I iterator_type;
- typedef typename I::value_type value_type;
- IterRange() {}
- IterRange(I beg, I end)
- : end_(end), cur_(beg) {
- }
- /// \brief evaluates to true when end of range is reached
- bool AtEnd() const {
- return cur_==end_;
- }
- /// \brief step to next element in range
- IterRange<I>& Next() {
- ++cur_;
- return *this;
- }
- typename I::value_type& Get() {
- return *cur_;
- }
- IterRange<I>& operator++() {
- ++cur_;
- return *this;
- }
- typename I::value_type& operator*() {
- return this->Get();
- }
-private:
- I end_;
- I cur_;
-};
-
-
-namespace impl {
- // forward declearation of chain impl list.
- typedef std::vector<ChainImplPtr> ChainImplList;
-}
-
-class DLLEXPORT_OST_MOL ChainHandleIter : public std::iterator<std::forward_iterator_tag,
- ChainHandle> {
-public: // internally used constructors
- ChainHandleIter(impl::ChainImplList::iterator chain_it)
- : cur_(chain_it) {
- }
-public:
- ChainHandle operator*() {
- return ChainHandle(*cur_);
- }
-
- ChainHandleIter& operator++() {
- ++cur_;
- return *this;
- }
-
- bool operator==(const ChainHandleIter& rhs) const {
- return cur_==rhs.cur_;
- }
-
- bool operator!=(const ChainHandleIter& rhs) const {
- return !this->operator==(rhs);
- }
-private:
- impl::ChainImplList::iterator cur_;
-};
-
-class DLLEXPORT_OST_MOL ResidueHandleIter : public std::iterator<std::forward_iterator_tag,
- ResidueHandle> {
-public:
- ResidueHandleIter(): cur_chain_(NULL), cur_res_(NULL) { }
- ResidueHandleIter(impl::pointer_it<impl::ChainImplPtr> chain_it,
- impl::pointer_it<impl::ResidueImplPtr> res_it,
- impl::EntityImplPtr ent, bool skip_empty);
-
- bool operator==(const ResidueHandleIter& rhs) const
- {
- return cur_res_==rhs.cur_res_;
- }
-
- void SkipEmpty();
-
- bool operator!=(const ResidueHandleIter& rhs) const
- {
- return !this->operator==(rhs);
- }
- ResidueHandleIter& operator++();
-
- ResidueHandle operator*();
-
-private:
- impl::pointer_it<impl::ChainImplPtr> cur_chain_;
- impl::pointer_it<impl::ResidueImplPtr> cur_res_;
- impl::EntityImplPtr ent_;
-};
-
-class DLLEXPORT_OST_MOL ResidueViewIter : public std::iterator<std::forward_iterator_tag,
- ResidueView> {
-public:
- ResidueViewIter(): cur_chain_(NULL), cur_res_(NULL) { }
-
- ResidueViewIter(impl::pointer_it<ChainView> chain_it,
- impl::pointer_it<ResidueView> res_it,
- EntityView ent, bool skip_empty);
-
- bool operator==(const ResidueViewIter& rhs) const {
- return cur_res_==rhs.cur_res_;
- }
- void SkipEmpty();
-
- bool operator!=(const ResidueViewIter& rhs) const {
- return !this->operator==(rhs);
- }
- ResidueViewIter& operator++();
-
- ResidueView operator*();
-
-private:
- impl::pointer_it<ChainView> cur_chain_;
- impl::pointer_it<ResidueView> cur_res_;
- EntityView ent_;
-};
-
-class DLLEXPORT_OST_MOL AtomHandleIter : public std::iterator<std::forward_iterator_tag,
- AtomHandle> {
-public:
- AtomHandleIter(): cur_chain_(NULL), cur_res_(NULL), cur_atom_(NULL) { }
-
- AtomHandleIter(impl::pointer_it<impl::ChainImplPtr> chain_it,
- impl::pointer_it<impl::ResidueImplPtr> res_it,
- impl::pointer_it<impl::AtomImplPtr> atom_it,
- impl::EntityImplPtr ent, bool skip_empty);
-
- bool operator==(const AtomHandleIter& rhs) const
- {
- return cur_atom_==rhs.cur_atom_;
- }
-
- bool operator!=(const AtomHandleIter& rhs) const
- {
- return !this->operator==(rhs);
- }
-
- AtomHandleIter& operator++();
-
- AtomHandle operator*();
-private:
- void SkipEmpty();
- impl::pointer_it<impl::ChainImplPtr> cur_chain_;
- impl::pointer_it<impl::ResidueImplPtr> cur_res_;
- impl::pointer_it<impl::AtomImplPtr> cur_atom_;
- impl::EntityImplPtr ent_;
-};
-
-class DLLEXPORT_OST_MOL AtomViewIter : public std::iterator<std::forward_iterator_tag,
- AtomView> {
-public:
- AtomViewIter(): cur_chain_(NULL), cur_res_(NULL), cur_atom_(NULL) { }
-
- AtomViewIter(impl::pointer_it<ChainView> chain_it,
- impl::pointer_it<ResidueView> res_it,
- impl::pointer_it<AtomView> atom_it,
- EntityView ent, bool skip_empty);
-
- bool operator==(const AtomViewIter& rhs) const
- {
- return cur_atom_==rhs.cur_atom_;
- }
-
- bool operator!=(const AtomViewIter& rhs) const
- {
- return !this->operator==(rhs);
- }
-
- AtomViewIter& operator++();
-
- AtomView operator*();
-private:
- void SkipEmpty();
- impl::pointer_it<ChainView> cur_chain_;
- impl::pointer_it<ResidueView> cur_res_;
- impl::pointer_it<AtomView> cur_atom_;
- EntityView ent_;
-};
-
-}} //ns
-
-
-#endif
diff --git a/modules/mol/base/src/iterator_fw.hh b/modules/mol/base/src/iterator_fw.hh
deleted file mode 100644
index c5fd2ef1d..000000000
--- a/modules/mol/base/src/iterator_fw.hh
+++ /dev/null
@@ -1,35 +0,0 @@
-//------------------------------------------------------------------------------
-// This file is part of the OpenStructure project <www.openstructure.org>
-//
-// Copyright (C) 2008-2011 by the OpenStructure authors
-//
-// This library is free software; you can redistribute it and/or modify it under
-// the terms of the GNU Lesser General Public License as published by the Free
-// Software Foundation; either version 3.0 of the License, or (at your option)
-// any later version.
-// This library is distributed in the hope that it will be useful, but WITHOUT
-// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
-// details.
-//
-// You should have received a copy of the GNU Lesser General Public License
-// along with this library; if not, write to the Free Software Foundation, Inc.,
-// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
-//------------------------------------------------------------------------------
-#ifndef OST_ITERATOR_FW_HH
-#define OST_ITERATOR_FW_HH
-
-namespace ost { namespace mol {
-
-
-class ChainHandleIter;
-class ResidueHandleIter;
-class AtomHandleIter;
-
-class ResidueViewIter;
-class AtomViewIter;
-class ChainViewIter;
-
-}} //ns
-
-#endif
diff --git a/modules/mol/base/src/mol.hh b/modules/mol/base/src/mol.hh
index cea4586e3..5e7b7a333 100644
--- a/modules/mol/base/src/mol.hh
+++ b/modules/mol/base/src/mol.hh
@@ -31,7 +31,6 @@
#include <ost/mol/chain_view.hh>
#include <ost/mol/residue_view.hh>
#include <ost/mol/atom_view.hh>
-#include <ost/mol/iterator.hh>
#include <ost/mol/xcs_editor.hh>
#include <ost/mol/ics_editor.hh>
#endif
diff --git a/modules/mol/base/src/query_state.cc b/modules/mol/base/src/query_state.cc
index 5f4ef1295..74dac4306 100644
--- a/modules/mol/base/src/query_state.cc
+++ b/modules/mol/base/src/query_state.cc
@@ -44,6 +44,18 @@ struct LazilyBoundData {
};
+void add_points_to_organizer(SpatialOrganizer<bool>& org, const EntityView& view) {
+ for (ChainViewList::const_iterator ci = view.GetChainList().begin(),
+ ce = view.GetChainList().end(); ci != ce; ++ci) {
+ for (ResidueViewList::const_iterator ri = ci->GetResidueList().begin(),
+ re = ci->GetResidueList().end(); ri != re; ++ri) {
+ for (AtomViewList::const_iterator ai = ri->GetAtomList().begin(),
+ ae = ri->GetAtomList().end(); ai != ae; ++ai) {
+ org.Add(true, ai->GetPos());
+ }
+ }
+ }
+}
bool cmp_string(CompOP op,const String& lhs, const StringOrRegexParam& rhs) {
switch (op) {
@@ -106,9 +118,7 @@ QueryState::QueryState(const QueryImpl& query, const EntityHandle& ref)
r_->refs.resize(query.bracketed_expr_.size());
for (size_t i=0;i<query.bracketed_expr_.size(); ++i) {
EntityView view=ref.Select(Query(query.bracketed_expr_[i]));
- for (AtomViewIter j=view.AtomsBegin(), e=view.AtomsEnd(); j!=e; ++j) {
- r_->refs[i].points.Add(true, (*j).GetPos());
- }
+ add_points_to_organizer(r_->refs[i].points, view);
}
}
}
@@ -126,9 +136,7 @@ QueryState::QueryState(const QueryImpl& query, const EntityView& ref)
r_->refs.resize(query.bracketed_expr_.size());
for (size_t i=0;i<query.bracketed_expr_.size(); ++i) {
EntityView view=ref.Select(Query(query.bracketed_expr_[i]));
- for (AtomViewIter j=view.AtomsBegin(), e=view.AtomsEnd(); j!=e; ++j) {
- r_->refs[i].points.Add(true, (*j).GetPos());
- }
+ add_points_to_organizer(r_->refs[i].points, view);
}
}
}
diff --git a/modules/mol/base/src/residue_handle.cc b/modules/mol/base/src/residue_handle.cc
index d2ef00da1..117979652 100644
--- a/modules/mol/base/src/residue_handle.cc
+++ b/modules/mol/base/src/residue_handle.cc
@@ -24,8 +24,8 @@
#include <ost/mol/impl/residue_impl.hh>
#include <ost/mol/impl/chain_impl.hh>
#include <ost/mol/impl/entity_impl.hh>
-#include <ost/mol/iterator.hh>
#include <ost/mol/bond_handle.hh>
+
namespace ost { namespace mol {
ResidueHandle::ResidueHandle()
@@ -163,36 +163,6 @@ bool ResidueHandle::SwitchAtomPos(const String& group)
}
-AtomHandleIter ResidueHandle::AtomsBegin() const
-{
- this->CheckValidity();
- impl::ResidueImplPtr r=Impl();
- int index=this->GetIndex();
- impl::ChainImplPtr c=r->GetChain();
- impl::pointer_it<impl::ChainImplPtr> cc=r->GetEntity()->GetChainIter(c->GetName());
- return AtomHandleIter(cc, impl::begin(c->GetResidueList())+index,
- impl::begin(r->GetAtomList()), r->GetEntity(), true);
-}
-
-AtomHandleIter ResidueHandle::AtomsEnd() const
-{
- this->CheckValidity();
- impl::ResidueImplPtr r=Impl();
- int index=this->GetIndex();
- impl::ChainImplPtr c=r->GetChain();
- impl::pointer_it<impl::ChainImplPtr> cc=r->GetEntity()->GetChainIter(c->GetName());
- if (c->GetResidueList().begin()+index+1==c->GetResidueList().end()) {
- return AtomHandleIter(cc, impl::begin(c->GetResidueList())+index,
- impl::begin(r->GetAtomList()), r->GetEntity(), false);
- } else {
- impl::pointer_it<impl::ResidueImplPtr> x=(impl::begin(c->GetResidueList())+index+1);
- return AtomHandleIter(cc, impl::begin(c->GetResidueList())+index+1,
- impl::begin((*x)->GetAtomList()),
- r->GetEntity(), false);
- }
-
-}
-
bool InSequence(const ResidueHandle& r1, const ResidueHandle& r2)
{
if (!r1 || !r2) return false; // protect against invalid handles
diff --git a/modules/mol/base/src/residue_handle.hh b/modules/mol/base/src/residue_handle.hh
index 5048a0804..6462036c8 100644
--- a/modules/mol/base/src/residue_handle.hh
+++ b/modules/mol/base/src/residue_handle.hh
@@ -26,7 +26,6 @@
#include <ost/mol/module_config.hh>
#include <ost/mol/residue_base.hh>
#include <ost/mol/entity_visitor_fw.hh>
-#include <ost/mol/iterator_fw.hh>
namespace ost { namespace mol {
@@ -167,22 +166,11 @@ public:
bool SwitchAtomPos(const String& group);
//@}
- /// \brief Get iterator pointing to the beginning of the atoms
- ///
- /// Atom iterators are not fail-safe, meaning that when new atoms are
- /// inserted the iterators become invalid and the behaviour is undefined. To
- /// iterate over the residue while inserting/deleting atoms use GetAtomList().
- /// \sa AtomsEnd(), GetAtomList()
- AtomHandleIter AtomsBegin() const;
-
/// \brief get this handle
///
/// This method exists as a convenience duck-typing in Python.
ResidueHandle GetHandle() const;
unsigned long GetHashCode() const;
- /// \brief Get iterator pointing to the beginning of the atoms
- /// \sa AtomsBegin(), GetAtomList()
- AtomHandleIter AtomsEnd() const;
/// \brief return view based on a query object
/// \sa Query
diff --git a/modules/mol/base/src/view_op.cc b/modules/mol/base/src/view_op.cc
index 6b8aebad8..83e02b41c 100644
--- a/modules/mol/base/src/view_op.cc
+++ b/modules/mol/base/src/view_op.cc
@@ -180,32 +180,37 @@ mol::EntityView Intersection(const mol::EntityView& ev1,
BondMap bm;
BondHandleList bl;
BondHandleList bonds_to_add;
- for (AtomViewIter a=ev1.AtomsBegin(),e=ev1.AtomsEnd(); a!=e; ++a) {
- AtomView av=*a;
- AtomView av2=ev2.ViewForHandle(av.GetHandle());
- if (av2.IsValid()) {
- intersection.AddAtom(av.GetHandle());
- bl=av.GetBondList();
- for (BondHandleList::iterator i=bl.begin(), e2=bl.end(); i!=e2; ++i) {
- BondMap::iterator x=bm.find(*i);
- if (x==bm.end()) {
- bm.insert(std::make_pair(*i, true));
- }
- }
- bl=av2.GetBondList();
- for (BondHandleList::iterator i=bl.begin(), e2=bl.end(); i!=e2; ++i) {
- BondMap::iterator x=bm.find(*i);
- if (x!=bm.end()) {
- if (x->second) {
- bonds_to_add.push_back(*i);
- x->second=false;
+ for (ChainViewList::const_iterator ci = ev1.GetChainList().begin(),
+ ce = ev1.GetChainList().end(); ci != ce; ++ci) {
+ for (ResidueViewList::const_iterator ri = ci->GetResidueList().begin(),
+ re = ci->GetResidueList().end(); ri != re; ++ri) {
+ for (AtomViewList::const_iterator ai = ri->GetAtomList().begin(),
+ ae = ri->GetAtomList().end(); ai != ae; ++ai) {
+ AtomView av=*ai;
+ AtomView av2=ev2.ViewForHandle(av.GetHandle());
+ if (av2.IsValid()) {
+ intersection.AddAtom(av.GetHandle());
+ bl=av.GetBondList();
+ for (BondHandleList::iterator i=bl.begin(), e2=bl.end(); i!=e2; ++i) {
+ BondMap::iterator x=bm.find(*i);
+ if (x==bm.end()) {
+ bm.insert(std::make_pair(*i, true));
+ }
+ }
+ bl=av2.GetBondList();
+ for (BondHandleList::iterator i=bl.begin(), e2=bl.end(); i!=e2; ++i) {
+ BondMap::iterator x=bm.find(*i);
+ if (x!=bm.end()) {
+ if (x->second) {
+ bonds_to_add.push_back(*i);
+ x->second=false;
+ }
+ } else {
+ bm.insert(std::make_pair(*i, true));
+ }
}
- } else {
- bm.insert(std::make_pair(*i, true));
}
- }
- }
- }
+ }}}
for (BondHandleList::iterator i=bonds_to_add.begin(),
e=bonds_to_add.end(); i!=e; ++i) {
intersection.AddBond(*i);
diff --git a/modules/mol/base/tests/CMakeLists.txt b/modules/mol/base/tests/CMakeLists.txt
index 176f6ad24..de287b52c 100644
--- a/modules/mol/base/tests/CMakeLists.txt
+++ b/modules/mol/base/tests/CMakeLists.txt
@@ -7,7 +7,6 @@ set(OST_MOL_BASE_UNIT_TESTS
test_delete.cc
test_entity.cc
test_ics.cc
- test_iterators.cc
test_query.cc
test_surface.cc
test_residue.cc
--
GitLab