Extract Eigen positions only once and then directly fill a pairwise RMSD matric for clustering

core/src/superpose.cc

@@ -48,6 +48,38 @@ inline Real GetRMSD(const promod3::core::EMatX3& pos_one,
   return std::sqrt(rmsd/pos_one.rows());
 }
 
+// calculate rmsd of the entries of two position matrices with applying a
+// rotation before
+inline Real GetRMSD(const promod3::core::EMatX3& pos_one, 
+                    const promod3::core::EMatX3& pos_two,
+                    const promod3::core::EMat3& rot){
+  Real rmsd = 0.0;
+
+  for(uint i = 0; i < pos_one.rows(); ++i){
+    rmsd += (pos_two.row(i) - 
+             (rot * pos_one.row(i).transpose()).transpose()).squaredNorm();
+  }
+  return std::sqrt(rmsd/pos_one.rows());
+}
+
+// calculate average row of the point set
+inline promod3::core::ERVec3 AverageRow(promod3::core::EMatX3& pos) {
+  // return pos.colwise().sum()/pos.rows();
+  // NOTE: the colwise-operator is opt-flag-dependent and we hence do it manually
+  promod3::core::ERVec3 avg_row = promod3::core::ERVec3::Zero();
+  for (uint i = 0; i < pos.rows(); ++i) {
+    avg_row += pos.row(i);
+  }
+  return avg_row / pos.rows();
+}
+
+inline void CenterPositions(promod3::core::EMatX3& pos){
+  promod3::core::ERVec3 avg = AverageRow(pos);
+  for(uint i = 0; i < pos.rows(); ++i){
+    pos.row(i) -= avg;
+  }
+}
+
 inline void ApplySuperposition(promod3::core::EMatX3& pos,
                                promod3::core::ERVec3& avg_one,
                                promod3::core::ERVec3& avg_two,
@@ -66,17 +98,6 @@ inline void ApplySuperposition(promod3::core::EMatX3& pos,
   } 
 }
 
-// calculate average row of the point set
-inline promod3::core::ERVec3 AverageRow(promod3::core::EMatX3& pos) {
-  // return pos.colwise().sum()/pos.rows();
-  // NOTE: the colwise-operator is opt-flag-dependent and we hence do it manually
-  promod3::core::ERVec3 avg_row = promod3::core::ERVec3::Zero();
-  for (uint i = 0; i < pos.rows(); ++i) {
-    avg_row += pos.row(i);
-  }
-  return avg_row / pos.rows();
-}
-
 // The Superpose function performs all tasks required for a superposition.
 // After calling the function, you can construct the transformation matrix 
 // using avg_one, avg_two and rotation.
@@ -299,6 +320,40 @@ std::pair<geom::Mat4,Real> Superposition(EMatX3& pos_one, EMatX3& pos_two,
   return std::make_pair(transform, rmsd); 
 }
 
+void FillRMSDMatrix(EMatX3List& position_list, Real** data){
+
+  if(position_list.empty()){
+    throw promod3::Error("Empty position list in FillRMSDMatrix!");
+  }
+
+  // check whether all positions have the same size
+  int num_pos = position_list.size();
+  int pos_size = position_list[0].rows(); 
+  for(int i = 1; i < num_pos; ++i){
+    if(position_list[i].rows() != pos_size){
+      throw promod3::Error("Size inconsistency in FillRMSDMatrix!");
+    }
+  }
+
+  // we first center all the positions
+  for(int i = 0; i < num_pos; ++i){
+    CenterPositions(position_list[i]);
+  }
+
+  // fill the data
+  Real rmsd;
+  EMat3 rot;
+  for(int i = 0; i < num_pos; ++i){
+    data[i][i] = 0.0;
+    for(int j = i + 1; j < num_pos; ++j){
+      GetRotationalComponent(position_list[i], position_list[j], rot);
+      rmsd = GetRMSD(position_list[i], position_list[j], rot);
+      data[i][j] = rmsd;
+      data[j][i] = rmsd;
+    }
+  }
+}
+
 Real SuperposedRMSD(EMatX3& pos_one, EMatX3& pos_two, 
                     uint max_iterations,
                     Real distance_thresh, 

core/src/superpose.hh

@@ -7,6 +7,9 @@
 
 #include <Eigen/SVD>
 #include <Eigen/Geometry>
+#include <Eigen/StdVector>
+
+#include <vector>
 
 namespace promod3 { namespace core {
 
@@ -15,6 +18,11 @@ typedef Eigen::Matrix<Real,3,1> EVec3;
 typedef Eigen::Matrix<Real,1,3> ERVec3;
 typedef Eigen::Matrix<Real,Eigen::Dynamic,3> EMatX3;
 
+// If you want to use stl containers of fixed sized Eigentype (e.g. EMatX3)
+// you must use a custom allocator provided by Eigen to ensure proper memory
+// alignment (more on that in the Eigen documentation)
+typedef std::vector<EMatX3,Eigen::aligned_allocator<EMatX3> > EMatX3List;
+
 // Following functions are intended to be used by other code. Fill the eigen
 // matrices by themselves to avoid shouffling around too much data.
 // NOTE: the EMatX3 objects get altered inside some of these functions!
@@ -33,6 +41,9 @@ geom::Mat4 MinRMSDSuperposition(EMatX3& pos_one, EMatX3& pos_two,
 std::pair<geom::Mat4,Real> Superposition(EMatX3& pos_one, EMatX3& pos_two,
 	                                     bool apply_superposition = false);
 
+// Fill matrix with pairwise RMSD values
+// The positions you pass will be centered, their mean will be the origin...
+void FillRMSDMatrix(EMatX3List& position_list, Real** data);
 
 // ITERATIVE MIN RMSD SUPERPOSITION ALGORITHMS
 

modelling/src/loop_candidate.cc

@@ -268,24 +268,30 @@ LoopCandidates::GetClusters(Real max_dist,
     clusters.resize(1);
     clusters[0].push_back(0);
   } else if (num_loops > 0) {
-    // get RMSD matrices
-    float** rmsd_matrix = new float*[num_loops];
+
+    // prepare the Eigen positions
+    promod3::core::EMatX3List position_list;
+    position_list.reserve(num_loops);
+    uint loop_size = sequence_.size();
     for(uint i = 0; i < num_loops; ++i){
-      rmsd_matrix[i] = new float[num_loops];
+      promod3::core::EMatX3 bb_list_positions = 
+      promod3::core::EMatX3::Zero(loop_size, 3);
+
+      for(uint j = 0; j < loop_size; ++j){
+        promod3::core::EMatFillRow(bb_list_positions, j, 
+                                   candidates_[i].GetCA(j));
       }
-    // set diagonal elements to zero
-    for (uint i = 0; i < num_loops; ++i) {
-      rmsd_matrix[i][i] = 0.0;
+
+      position_list.push_back(bb_list_positions);
     }
-    // set rest
-    Real rmsd;
+
+    Real** rmsd_matrix = new Real*[num_loops];
     for (uint i= 0; i < num_loops; ++i) {
-      for (uint j = i+1; j < num_loops; ++j) {
-        rmsd = candidates_[i].CARMSD(candidates_[j]);
-        rmsd_matrix[i][j] = rmsd;
-        rmsd_matrix[j][i] = rmsd;
-      }
+      rmsd_matrix[i] = new Real[num_loops];
     }
+
+    promod3::core::FillRMSDMatrix(position_list, rmsd_matrix);
+
     // get clusters
     promod3::core::DoClustering(rmsd_matrix, num_loops,
                                 promod3::core::AVG_DIST_METRIC,