diff --git a/modules/gui/pymod/dng/menu.py b/modules/gui/pymod/dng/menu.py index 59e3e15f13d0ded4577fb42a24f09dab1011ac8f..875d43cf2a4fd213839a3f292f91f30a39d43670 100644 --- a/modules/gui/pymod/dng/menu.py +++ b/modules/gui/pymod/dng/menu.py @@ -190,12 +190,22 @@ class SceneMenu(QMenu): i += 1 gfx_ent_2 = sel.GetActiveNode(i) sd = superpositiondialog.SuperpositionDialog(gfx_ent_1, gfx_ent_2) - if sd.rmsd != None: + if sd.rmsd_superposed_atoms != None: if sd.reference == 0: - gfx_ent_1.UpdatePositions() + gfx_ent_2.UpdatePositions() gfx.Scene().CenterOn(gfx_ent_1) else: + gfx_ent_1.UpdatePositions() + gfx.Scene().CenterOn(gfx_ent_2) + LogScript('RMSD: %.3f'%sd.rmsd) + LogScript('RMSD superposed atoms: %.3f'%sd.rmsd_superposed_atoms) + LogScript('fraction superposed: %.3f'%sd.fraction_superposed) + elif sd.rmsd != None: + if sd.reference == 0: gfx_ent_2.UpdatePositions() + gfx.Scene().CenterOn(gfx_ent_1) + else: + gfx_ent_1.UpdatePositions() gfx.Scene().CenterOn(gfx_ent_2) LogScript('RMSD: %.3f'%sd.rmsd) diff --git a/modules/gui/pymod/dng/superpositiondialog.py b/modules/gui/pymod/dng/superpositiondialog.py index 9ad8bc22fd0cf0d30c290f864367e4e015e2627f..7eb3918038f74ab47302c27011770e5c04b69e5d 100644 --- a/modules/gui/pymod/dng/superpositiondialog.py +++ b/modules/gui/pymod/dng/superpositiondialog.py @@ -121,7 +121,9 @@ class SuperpositionDialog(QDialog): def __init__(self, ent_one, ent_two, parent=None): # class variables + self.rmsd_superposed_atoms = None self.rmsd = None + self.fraction_superposed = None self._mmethod_dict = {'number': 'number', 'index': 'index', 'local alignment': 'local-aln', @@ -189,6 +191,11 @@ class SuperpositionDialog(QDialog): layout.addWidget(QLabel('match residues by'), grow, 0) grow += 1 layout.addWidget(self._methods) + + # iterative + self._iterative=None + self._it_box, self._it_in, self._dist_in = self._ItBox() + layout.addWidget(self._it_box, grow, 0) # atoms self._atoms = self._FetchAtoms(self._methods.size(), self.ent_one, @@ -222,8 +229,12 @@ class SuperpositionDialog(QDialog): atoms = self._GetAtomSelection() sp = Superpose(view_two, view_one, self._mmethod_dict[str(self._methods.currentText())], - atoms) + atoms, iterative=self._iterative, + max_iterations=self._it_in.value(), distance_threshold=self._dist_in.value()) self.rmsd = sp.rmsd + if self._iterative: + self.rmsd_superposed_atoms = sp.rmsd_superposed_atoms + self.fraction_superposed = sp.fraction_superposed def _toggle_atoms(self, checked): if checked: @@ -296,6 +307,44 @@ class SuperpositionDialog(QDialog): cbox.setCurrentIndex(0) return cbox + def _toggle_iterative(self, checked): + if checked: + self._it_in.setEnabled(True) + self._dist_in.setEnabled(True) + self._iterative=True + else: + self._it_in.setEnabled(False) + self._dist_in.setEnabled(False) + self._iterative=False + + def _ItBox(self): + bt1 = QRadioButton("On") + iteration_label=QLabel("Max Iterations: ") + distance_label=QLabel("Dist Thresh: ") + iteration_in=QSpinBox() + iteration_in.setRange(1,30) + iteration_in.setValue(8) + distance_in=QDoubleSpinBox() + distance_in.setRange(1.0,10.0) + distance_in.setValue(3.0) + distance_in.setDecimals(1) + distance_in.setSingleStep(0.5) + iteration_in.setEnabled(False) + distance_in.setEnabled(False) + bt1.setChecked(False) + self._iterative=False + vbox_layout = QVBoxLayout() + vbox_layout.addWidget(bt1) + vbox_layout.addWidget(iteration_label) + vbox_layout.addWidget(iteration_in) + vbox_layout.addWidget(distance_label) + vbox_layout.addWidget(distance_in) + vbox_layout.addSpacing(50) + QObject.connect(bt1, SIGNAL('toggled(bool)'), self._toggle_iterative) + box = QGroupBox("Iterative") + box.setLayout(vbox_layout) + return box,iteration_in, distance_in + def _ChangeChainSelection(self, index): if index == 0: self._chain_one.SetItems(self.ent_one, self.gfx_one) diff --git a/modules/mol/alg/pymod/export_svd_superpose.cc b/modules/mol/alg/pymod/export_svd_superpose.cc index 443ca71bccad8d13cda6fae3223d611bc6efab44..9c4087021f959a7906384d77dda7e6ca9a104541 100644 --- a/modules/mol/alg/pymod/export_svd_superpose.cc +++ b/modules/mol/alg/pymod/export_svd_superpose.cc @@ -39,6 +39,8 @@ void export_svdSuperPose() class_<SuperpositionResult>("SuperpositionResult", init<>()) .def_readwrite("ncycles", &SuperpositionResult::ncycles) .def_readwrite("rmsd", &SuperpositionResult::rmsd) + .def_readwrite("fraction_superposed", &SuperpositionResult::fraction_superposed) + .def_readwrite("rmsd_superposed_atoms", &SuperpositionResult::rmsd_superposed_atoms) .def_readwrite("transformation", &SuperpositionResult::transformation) .def_readwrite("view1", @@ -49,6 +51,7 @@ void export_svdSuperPose() def("SuperposeAtoms", &SuperposeAtoms,(arg("apply_transform")=true)); def("SuperposeSVD", sup1, (arg("apply_transform")=true)); def("SuperposeSVD", sup2, (arg("apply_transform")=true)); + def("IterativeSuperposeSVD", &IterativeSuperposeSVD,(arg("max_iterations")=5, arg("distance_threshold")=3.0,arg("apply_transform")=true)); def("CalculateRMSD", &CalculateRMSD, (arg("transformation")=geom::Mat4())); } diff --git a/modules/mol/alg/pymod/superpose.py b/modules/mol/alg/pymod/superpose.py index 76d4afe98d53178bba1cf357e2e4f827ccef7f01..4f1ac1e0ac315ae8d430292634082dd568c2a503 100644 --- a/modules/mol/alg/pymod/superpose.py +++ b/modules/mol/alg/pymod/superpose.py @@ -259,7 +259,7 @@ def MatchResidueByGlobalAln(ent_a, ent_b, atoms='all'): return _MatchResidueByAln(ent_a, ent_b, atoms, ost.seq.alg.GlobalAlign) -def Superpose(ent_a, ent_b, match='number', atoms='all'): +def Superpose(ent_a, ent_b, match='number', atoms='all', iterative=False, max_iterations=5, distance_threshold=3.0): """ Superposes the model entity onto the reference. To do so, two views are created, returned with the result. **atoms** describes what goes in to these @@ -310,5 +310,8 @@ def Superpose(ent_a, ent_b, match='number', atoms='all'): else: raise ValueError(not_supported) ## action - res=ost.mol.alg.SuperposeSVD(view_a, view_b) + if iterative: + res=ost.mol.alg.IterativeSuperposeSVD(view_a, view_b, max_iterations, distance_threshold) + else: + res=ost.mol.alg.SuperposeSVD(view_a, view_b) return res diff --git a/modules/mol/alg/src/svd_superpose.cc b/modules/mol/alg/src/svd_superpose.cc index 56c0796c6f2462837143307b1bb33bd17ed905bf..cb6ee69304edb62132e5f7a7c5d37d2243f1de41 100644 --- a/modules/mol/alg/src/svd_superpose.cc +++ b/modules/mol/alg/src/svd_superpose.cc @@ -25,10 +25,12 @@ #include <Eigen/Array> #include <Eigen/SVD> #include <Eigen/LU> +#include <Eigen/Dense> #include <ost/base.hh> #include <ost/geom/vec3.hh> +#include <ost/geom/mat4.hh> #include <ost/mol/alg/svd_superpose.hh> #include <ost/mol/xcs_editor.hh> #include <ost/mol/residue_handle.hh> @@ -44,6 +46,7 @@ namespace ost { namespace mol { namespace alg { using boost::bind; typedef Eigen::Matrix<Real,3,1> EVec3; typedef Eigen::Matrix<Real,3,3> EMat3; +typedef Eigen::Matrix<Real,4,4> EMat4; typedef Eigen::Matrix<Real,1,3> ERVec3; typedef Eigen::Matrix<Real,Eigen::Dynamic,Eigen::Dynamic> EMatX; typedef Eigen::Matrix<Real,1,Eigen::Dynamic> ERVecX; @@ -89,6 +92,32 @@ Real calc_rmsd_for_atom_lists(const mol::AtomViewList& atoms1, return rmsd; } +Real calc_rmsd_for_ematx(const EMatX& atoms1, + const EMatX& atoms2, + const EMat4& transformation) +{ + EMatX transformed_atoms1 = EMatX::Zero(atoms1.rows(), 3); + + EMatX vector = EMatX::Zero(4,1); + EMatX transformed_vector = EMatX::Zero(4,1); + vector(3,0)=1; + + for(int i=0;i<atoms1.rows();++i){ + vector.block<3,1>(0,0)=atoms1.block<1,3>(i,0).transpose(); + transformed_vector = transformation*vector; + transformed_atoms1.block<1,3>(i,0)=transformed_vector.block<3,1>(0,0).transpose(); + } + + EMatX diff = EMatX::Zero(atoms1.rows(),atoms1.cols()); + EMatX squared_dist = EMatX::Zero(atoms1.rows(),1); + + diff = transformed_atoms1-atoms2; + squared_dist = (diff.cwise()*diff).rowwise().sum(); + + return sqrt(squared_dist.sum()/squared_dist.rows()); + +} + Real CalculateRMSD(const mol::EntityView& ev1, const mol::EntityView& ev2, const geom::Mat4& transformation) { @@ -104,7 +133,36 @@ geom::Vec3 EigenVec3ToVec3(const EVec3 &vec) geom::Mat3 EigenMat3ToMat3(const EMat3 &mat) { - return geom::Mat3(mat.data()); + geom::Mat3 return_mat; + for(int i=0;i<3;++i){ + for(int j=0;j<3;++j){ + return_mat(i,j) = mat(i,j); + } + } + return return_mat; + //return geom::Mat3(mat.data()); +} + +geom::Mat4 EigenMat4ToMat4(const EMat4 &mat) +{ + geom::Mat4 return_mat; + for(int i=0;i<4;++i){ + for(int j=0;j<4;++j){ + return_mat(i,j) = mat(i,j); + } + } + return return_mat; + //return geom::Mat4(mat.data()); +} + +EMatX Mat4ToEigenMat4(const geom::Mat4 &mat){ + EMat4 res = EMat4::Zero(); + for(int i=0;i<4;++i){ + for(int j=0;j<4;++j){ + res(i,j)=mat.At(i,j); + } + } + return res; } EVec3 Vec3ToEigenRVec(const geom::Vec3 &vec) @@ -126,6 +184,7 @@ EMatX MatrixShiftedBy(EMatX mat, ERVecX vec) return result; } + class MeanSquareMinimizerImpl { public: MeanSquareMinimizerImpl(int n_atoms, bool alloc_atoms): @@ -145,7 +204,16 @@ public: atoms1_.row(index) = ERVec3(Vec3ToEigenVec(pos)); } + SuperpositionResult Minimize(const EMatX& atoms, const EMatX& atoms_ref) const; + + EMatX TransformEMatX(const EMatX& mat, const EMat4& transformation) const; + + std::pair<EMatX,EMatX> CreateMatchingSubsets(const EMatX& atoms, const EMatX& atoms_ref, Real distance_threshold) const; + + SuperpositionResult IterativeMinimize(int ncycles, Real distance_threshold) const; + SuperpositionResult MinimizeOnce() const; + private: int n_atoms_; bool alloc_atoms_; @@ -155,9 +223,11 @@ private: MeanSquareMinimizer MeanSquareMinimizer::FromAtomLists(const AtomViewList& atoms, - const AtomViewList& atoms_ref) { + const AtomViewList& atoms_ref) +{ int n_atoms = atoms.size(); - if (n_atoms != atoms_ref.size()) { + int n_atoms_ref = atoms_ref.size(); + if (n_atoms != n_atoms_ref) { throw Error("atom counts do not match"); } if (n_atoms<3) { @@ -174,9 +244,11 @@ MeanSquareMinimizer MeanSquareMinimizer::FromAtomLists(const AtomViewList& atoms } MeanSquareMinimizer MeanSquareMinimizer::FromPointLists(const std::vector<geom::Vec3>& points, - const std::vector<geom::Vec3>& points_ref) { + const std::vector<geom::Vec3>& points_ref) +{ int n_points = points.size(); - if (n_points != points.size()) { + int n_points_ref = points_ref.size(); + if (n_points != n_points_ref) { throw Error("point counts do not match"); } if (n_points<3) { @@ -192,19 +264,126 @@ MeanSquareMinimizer MeanSquareMinimizer::FromPointLists(const std::vector<geom:: return msm; } -SuperpositionResult MeanSquareMinimizerImpl::MinimizeOnce() const { - ERVec3 avg1 = atoms1_.colwise().sum()/atoms1_.rows(); - ERVec3 avg2 = atoms2_.colwise().sum()/atoms2_.rows(); +SuperpositionResult MeanSquareMinimizerImpl::MinimizeOnce() const{ + return this->Minimize(atoms1_,atoms2_); +} + +SuperpositionResult MeanSquareMinimizerImpl::IterativeMinimize(int max_cycles, Real distance_threshold) const{ + + std::vector<EMat4> transformation_matrices; + EMat4 transformation_matrix; + EMatX atoms = atoms1_; + SuperpositionResult res; + EMat4 diff; + std::pair<EMatX,EMatX> subsets; + EMat4 identity_matrix = EMat4::Identity(); + + //do initial superposition + + res = this->Minimize(atoms, atoms2_); + transformation_matrices.push_back(Mat4ToEigenMat4(res.transformation)); + + //note, that the initial superposition is the first cycle... + int cycles=1; + + for(;cycles<max_cycles;++cycles){ + atoms = this->TransformEMatX(atoms, transformation_matrices.back()); + subsets = this->CreateMatchingSubsets(atoms, atoms2_, distance_threshold); + res = this->Minimize(subsets.first,subsets.second); + transformation_matrix = Mat4ToEigenMat4(res.transformation); + transformation_matrices.push_back(transformation_matrix); + + diff = transformation_matrix-identity_matrix; + + if(diff.cwise().abs().sum()<0.001){ + break; + } + + } + + res.rmsd_superposed_atoms = calc_rmsd_for_ematx(subsets.first, subsets.second, transformation_matrices.back()); + res.fraction_superposed = float(subsets.first.rows())/atoms1_.rows(); + + //combine the transformations into one transformation + transformation_matrix = transformation_matrices.back(); + transformation_matrices.pop_back(); + while(!transformation_matrices.empty()){ + transformation_matrix*=transformation_matrices.back(); + transformation_matrices.pop_back(); + } + + res.transformation = EigenMat4ToMat4(transformation_matrix); + res.ncycles=cycles; + + return res; + +} + +EMatX MeanSquareMinimizerImpl::TransformEMatX(const EMatX& mat, const EMat4& transformation) const { + + EMatX transformed_mat = EMatX::Zero(mat.rows(), 3); + + EMatX vector = EMatX::Zero(4,1); + EMatX transformed_vector = EMatX::Zero(4,1); + vector(3,0)=1; + + for(int i=0;i<mat.rows();++i){ + vector.block<3,1>(0,0)=mat.block<1,3>(i,0).transpose(); + transformed_vector = transformation*vector; + transformed_mat.block<1,3>(i,0)=transformed_vector.block<3,1>(0,0).transpose(); + } + + return transformed_mat; +} + +std::pair<EMatX, EMatX> MeanSquareMinimizerImpl::CreateMatchingSubsets(const EMatX& atoms, const EMatX& atoms_ref, Real distance_threshold) const{ + + + EMatX diff = EMatX::Zero(atoms.rows(),atoms.cols()); + EMatX dist = EMatX::Zero(atoms.rows(),1); + + diff = atoms-atoms_ref; + dist = (diff.cwise()*diff).rowwise().sum(); + dist = dist.cwise().sqrt(); + + for(int i = 0; i < dist.rows(); ++i){ + if(dist(i,0) <= distance_threshold){ + dist(i,0) = 1; + } + else{ + dist(i,0) = 0; + } + } + + EMatX atoms_subset = EMatX::Zero(int(dist.sum()),3); + EMatX atoms_ref_subset = EMatX::Zero(int(dist.sum()),3); + + int actual_pos=0; + + for(int i = 0; i < dist.rows() ; ++i){ + if(dist(i,0)==1){ + atoms_subset.row(actual_pos) = atoms.row(i); + atoms_ref_subset.row(actual_pos) = atoms_ref.row(i); + ++actual_pos; + } + } + return std::make_pair(atoms_subset, atoms_ref_subset); +} + + +SuperpositionResult MeanSquareMinimizerImpl::Minimize(const EMatX& atoms, const EMatX& atoms_ref) const { + ERVec3 avg = atoms.colwise().sum()/atoms.rows(); + ERVec3 avg_ref = atoms_ref.colwise().sum()/atoms_ref.rows(); // SVD only determines the rotational component of the superposition // we need to manually shift the centers of the two point sets on onto // origin - EMatX atoms1 = MatrixShiftedBy(atoms1_, avg1); - EMatX atoms2 = MatrixShiftedBy(atoms2_, avg2).transpose(); + EMatX atoms_shifted = MatrixShiftedBy(atoms, avg); + EMatX atoms_ref_shifted = MatrixShiftedBy(atoms_ref, avg_ref).transpose(); // determine rotational component - Eigen::SVD<EMat3> svd(atoms2*atoms1); + Eigen::SVD<EMat3> svd(atoms_ref_shifted*atoms_shifted); EMatX matrixVT=svd.matrixV().transpose(); //determine rotation @@ -222,9 +401,10 @@ SuperpositionResult MeanSquareMinimizerImpl::MinimizeOnce() const { SuperpositionResult res; - geom::Vec3 shift = EigenVec3ToVec3(avg2); - geom::Vec3 com_vec = -EigenVec3ToVec3(avg1); - geom::Mat3 rot = EigenMat3ToMat3(rotation.transpose()); + geom::Vec3 shift = EigenVec3ToVec3(avg_ref); + geom::Vec3 com_vec = -EigenVec3ToVec3(avg); + //geom::Mat3 rot = EigenMat3ToMat3(rotation.transpose()); + geom::Mat3 rot = EigenMat3ToMat3(rotation); geom::Mat4 mat4_com, mat4_rot, mat4_shift; mat4_rot.PasteRotation(rot); mat4_shift.PasteTranslation(shift); @@ -256,6 +436,10 @@ SuperpositionResult MeanSquareMinimizer::MinimizeOnce() const { return impl_->MinimizeOnce(); } +SuperpositionResult MeanSquareMinimizer::IterativeMinimize(int ncycles, Real distance_threshold) const { + return impl_->IterativeMinimize(ncycles, distance_threshold); +} + SuperpositionResult SuperposeAtoms(const mol::AtomViewList& atoms1, const mol::AtomViewList& atoms2, bool apply_transform=true) @@ -265,6 +449,8 @@ SuperpositionResult SuperposeAtoms(const mol::AtomViewList& atoms1, result.ncycles=1; result.rmsd = calc_rmsd_for_atom_lists(atoms1, atoms2, result.transformation); + result.rmsd_superposed_atoms = result.rmsd; + result.fraction_superposed = 1.0; if (apply_transform) { mol::AtomView jv=atoms1.front(); mol::XCSEditor ed=jv.GetEntity().GetHandle().EditXCS(); @@ -297,5 +483,28 @@ SuperpositionResult SuperposeSVD(const std::vector<geom::Vec3>& pl1, return result; } +SuperpositionResult IterativeSuperposeSVD(const mol::EntityView& ev, + const mol::EntityView& ev_ref, + int max_cycles, + Real distance_threshold, + bool apply_transform){ + + AtomViewList atoms = ev.GetAtomList(); + AtomViewList atoms_ref = ev_ref.GetAtomList(); + + MeanSquareMinimizer msm = MeanSquareMinimizer::FromAtomLists(atoms, atoms_ref); + SuperpositionResult result = msm.IterativeMinimize(max_cycles, distance_threshold); + result.rmsd = calc_rmsd_for_atom_lists(atoms, atoms_ref, result.transformation); + if (apply_transform) { + mol::AtomView jv=atoms.front(); + mol::XCSEditor ed=jv.GetEntity().GetHandle().EditXCS(); + ed.ApplyTransform(result.transformation); + } + + result.entity_view1 = ev; + result.entity_view2 = ev_ref; + return result; +} + }}} //ns diff --git a/modules/mol/alg/src/svd_superpose.hh b/modules/mol/alg/src/svd_superpose.hh index 8839a8ef5fb15f1c9655cd2f665fb9f595b96766..bdbaf505de02ef8ab6928f268ba9fd63de0d74ee 100644 --- a/modules/mol/alg/src/svd_superpose.hh +++ b/modules/mol/alg/src/svd_superpose.hh @@ -45,6 +45,8 @@ struct DLLEXPORT_OST_MOL_ALG SuperpositionResult { } int ncycles; Real rmsd; + Real rmsd_superposed_atoms; + Real fraction_superposed; geom::Mat4 transformation; mol::EntityView entity_view1; mol::EntityView entity_view2; @@ -64,6 +66,8 @@ public: const std::vector<geom::Vec3>& points_ref); SuperpositionResult MinimizeOnce() const; + SuperpositionResult IterativeMinimize(int ncycles=5, Real distance_threshold=2.0) const; + ~MeanSquareMinimizer(); @@ -92,16 +96,16 @@ SuperpositionResult DLLEXPORT_OST_MOL_ALG SuperposeSVD(const std::vector<geom::V const std::vector<geom::Vec3>& pl2); /// \brief iterative superposition -SuperpositionResult DLLEXPORT_OST_MOL_ALG IterativeSuperposition(mol::EntityView& ev1, - mol::EntityView& ev2, - int ncycles, - Real distance_threshold, - bool apply_transform); +SuperpositionResult DLLEXPORT_OST_MOL_ALG IterativeSuperposeSVD(const mol::EntityView& ev1, + const mol::EntityView& ev2, + int max_cycles, + Real distance_threshold, + bool apply_transform); /// \brief calculates RMSD for two entity view Real DLLEXPORT_OST_MOL_ALG CalculateRMSD(const mol::EntityView& ev1, - const mol::EntityView& ev2, - const geom::Mat4& transformation); + const mol::EntityView& ev2, + const geom::Mat4& transformation);