port openstructure dependency from Eigen2 to Eigen3

changes:

-GetPrincipalAxes in Vec3List
 Results with Eigen3 have been verified with 10 randomly picked pdb structures.
 Differences only occur with changing signs of the vectors.

-CalcPrincipalAxis in BoundingBox
 No difference in algorithm, the SelfAdjointEigenSolver no lives in
 Eigenvalues. Does in principal the same thing as the function above,
 just using different functionality.

-SVD superposer
 Unit tests still run through with the changes introduced for Eigen3
 and also manual testing seems to work well.

-CartoonRenderer
 Principal axes calculation for rendering helices as cylinders had to be
 adapted. Rendering looks exactly the same...

-SuperposeFrames
 Not yet tested...

ToDo:
numeric.h and levenberg_marquardt.h in the img/alg module have to be
ported. This will be done in the next few days.

CMakeLists.txt

@@ -154,8 +154,6 @@ else()
   set(_UBUNTU_LAYOUT OFF)
 endif()
 
-add_definitions(-DEIGEN2_SUPPORT)
-
 if (COMPOUND_LIB)
   set(_COMP_LIB "${COMPOUND_LIB}")
   if (NOT IS_ABSOLUTE "${COMPOUND_LIB}")

modules/geom/src/vec3.cc

@@ -60,7 +60,7 @@ Mat3 Vec3List::GetPrincipalAxes() const
   Mat3 inertia=this->GetInertia();  
   EMat3 inertia_mat(inertia.Data());
 
-  Eigen::SVD<EMat3> svd(inertia_mat);
+  Eigen::JacobiSVD<EMat3> svd(inertia_mat,Eigen::ComputeThinU);
   EMat3 rot=svd.matrixU();
   Mat3 axes(rot.data());
   return axes;

modules/gfx/src/impl/cartoon_renderer.cc

@@ -23,10 +23,7 @@
 
 #include "cartoon_renderer.hh"
 
-#include <Eigen/Core>
-#include <Eigen/Array>
 #include <Eigen/SVD>
-#include <Eigen/LU>
 
 #include <ost/gfx/entity.hh>
 #include <ost/gfx/impl/tabulated_trig.hh>
@@ -138,7 +135,7 @@ namespace {
       A.row(i)=to_eigen(points[i]-cen);
     }
 
-    Eigen::SVD<EMatX> svd(A);
+    Eigen::JacobiSVD<EMatX> svd(A,Eigen::ComputeThinV);
     EMatX V=svd.matrixV();
     geom::Vec3 ax(V(0,0),V(1,0),V(2,0));
 

modules/mol/alg/src/superpose_frames.cc

@@ -17,10 +17,7 @@
 // 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 //------------------------------------------------------------------------------
 
-#include <Eigen/Core>
-#include <Eigen/Array>
-#include <Eigen/SVD>
-#include <Eigen/LU>
+#include <Eigen/Eigenvalues>
 #include <ost/message.hh>
 #include <ost/mol/mol.hh>
 #include <ost/mol/alg/superpose_frames.hh>
@@ -100,7 +97,7 @@ void AddSuperposedFrame(CoordGroupHandle& superposed, EMatX3& ref_mat,EMatX3& re
   frame_center=frame_mat.rowwise().sum()/frame_mat.cols();
   frame_centered=row_sub(frame_mat, frame_center);
   //single value decomposition
-  Eigen::SVD<EMat3> svd(frame_centered*ref_centered);
+  Eigen::JacobiSVD<EMat3> svd(frame_centered*ref_centered,Eigen::ComputeThinU | Eigen::ComputeThinV);
   EMat3 matrixVT=svd.matrixV().transpose();
   //determine rotation
   Real detv=matrixVT.determinant();

modules/mol/alg/src/svd_superpose.cc

@@ -21,10 +21,8 @@
 #include <iostream>
 #include <boost/bind.hpp>
 
-#include <Eigen/Core>
-#include <Eigen/Array>
 #include <Eigen/SVD>
-#include <Eigen/LU>
+#include <Eigen/Geometry>
 
 #include <ost/base.hh>
 #include <ost/geom/vec3.hh>
@@ -94,7 +92,7 @@ Real calc_rmsd_for_ematx(const EMatX& atoms1,
                          const EMatX& atoms2,
                          const EMat4& transformation)
 {
-  EMatX transformed_atoms1 = EMatX::Zero(atoms1.rows(), 3);
+  EMatX transformed = EMatX::Zero(atoms1.rows(), 3);
 
   EMatX vector = EMatX::Zero(4,1);
   EMatX transformed_vector = EMatX::Zero(4,1);
@@ -103,17 +101,12 @@ Real calc_rmsd_for_ematx(const EMatX& atoms1,
   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();
+    transformed.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());
-
+  transformed = transformed - atoms2;
+  transformed = transformed.array()*transformed.array();
+  return sqrt(transformed.rowwise().sum().sum()/atoms1.rows());
 }
 
 Real CalculateRMSD(const mol::EntityView& ev1,
@@ -204,7 +197,8 @@ public:
 
   EMatX TransformEMatX(const EMatX& mat, const EMat4& transformation) const;
 
-  std::vector<int> CreateMatchingSubsets(const EMatX& atoms, const EMatX& atoms_ref, Real distance_threshold) const;
+  void CreateMatchingSubsets(std::vector<int>& matching_subset, const EMatX& atoms, 
+                             const EMatX& atoms_ref, Real distance_threshold) const;
 
   SuperpositionResult IterativeMinimize(int ncycles, Real distance_threshold) const;
 
@@ -287,7 +281,7 @@ SuperpositionResult MeanSquareMinimizerImpl::IterativeMinimize(int max_cycles, R
   for(;cycles<max_cycles;++cycles){
 
     transformed_atoms = this->TransformEMatX(atoms1_, transformation_matrix);
-    matching_indices = this->CreateMatchingSubsets(transformed_atoms, atoms2_, distance_threshold);
+    this->CreateMatchingSubsets(matching_indices, transformed_atoms, atoms2_, distance_threshold);
 
     if(matching_indices.size()<3){
       std::stringstream ss;
@@ -312,7 +306,7 @@ SuperpositionResult MeanSquareMinimizerImpl::IterativeMinimize(int max_cycles, R
 
     diff = transformation_matrix_old-transformation_matrix;
 
-    if(diff.cwise().abs().sum()<0.0001){
+    if(diff.array().abs().sum()<0.0001){
       ++cycles;
       break;
     }
@@ -348,24 +342,25 @@ EMatX MeanSquareMinimizerImpl::TransformEMatX(const EMatX& mat, const EMat4& tra
   return transformed_mat;
 }
 
-std::vector<int> MeanSquareMinimizerImpl::CreateMatchingSubsets(const EMatX& atoms, const EMatX& atoms_ref, Real distance_threshold) const{
+void MeanSquareMinimizerImpl::CreateMatchingSubsets(std::vector<int>& matching_subset, 
+                                                    const EMatX& atoms, const EMatX& atoms_ref, 
+                                                    Real distance_threshold) const{
 
 
-  std::vector<int> return_vec;
+  matching_subset.clear();
   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();
+  dist = (diff.array()*diff.array()).rowwise().sum(); 
+
+  Real squared_dist_threshold = distance_threshold * distance_threshold;
   
   for(int i=0; i<dist.rows(); ++i){
-    if(dist(i,0) <= distance_threshold){
-      return_vec.push_back(i);
+    if(dist(i,0) <= squared_dist_threshold){
+      matching_subset.push_back(i);
     }
   }
-
-  return return_vec;
 }
 
 
@@ -381,7 +376,7 @@ SuperpositionResult MeanSquareMinimizerImpl::Minimize(const EMatX& atoms, const
   EMatX atoms_ref_shifted = MatrixShiftedBy(atoms_ref, avg_ref).transpose();
 
   // determine rotational component
-  Eigen::SVD<EMat3> svd(atoms_ref_shifted*atoms_shifted);
+  Eigen::JacobiSVD<EMat3> svd(atoms_ref_shifted*atoms_shifted,Eigen::ComputeThinU | Eigen::ComputeThinV);
   EMatX matrixVT=svd.matrixV().transpose();
 
   //determine rotation

modules/mol/base/src/bounding_box.cc

@@ -22,7 +22,7 @@
  */
 
 #include <limits>
-#include <Eigen/QR>
+#include <Eigen/Eigenvalues>
 #include <ost/mol/mol.hh>
 #include <ost/mol/bounding_box.hh>
 namespace ost { namespace mol {