optimize computation and memory (#43)

* optimize computation and memory * access tiles by const reference * add const version of operator[] in tile * missing const * solve comment * replace map with array * initialize to 0 array elements * avoid expensive floating point operations * adding declaration in header file --------- Co-authored-by: Erica Brondolin <[email protected]>
key4hep · Aug 15, 2023 · 4dafad2 · 4dafad2
1 parent 6f75053
commit 4dafad2
Show file tree

Hide file tree

Showing 3 changed files with 42 additions and 20 deletions.
diff --git a/include/CLUEAlgo.h b/include/CLUEAlgo.h
@@ -150,6 +150,8 @@ class CLUEAlgo_T {
   void calculateDistanceToHigher(TILES & );
   void findAndAssignClusters();
   inline float distance(int i, int j, bool isPhi = false, float r = 0.0) const ;
+  inline float distance2(int i, int j, bool isPhi = false, float r = 0.0) const ;
+
 };
 
 using CLUEAlgo = CLUEAlgo_T<LayerTiles>;

diff --git a/include/LayerTiles.h b/include/LayerTiles.h
@@ -80,7 +80,7 @@ class LayerTiles_T {
       return phiBin + yBin*T::nColumnsPhi;
     }
 
-    std::array<int,4> searchBox(float xMin, float xMax, float yMin, float yMax){
+    std::array<int,4> searchBox(float xMin, float xMax, float yMin, float yMax) const {
       int xBinMin = getXBin(xMin);
       int xBinMax = getXBin(xMax);
       int yBinMin = getYBin(yMin);
@@ -118,6 +118,10 @@ class LayerTiles_T {
     std::vector<int>& operator[](int globalBinId) {
       return layerTiles_[globalBinId];
     }
+
+    const std::vector<int>& operator[](int globalBinId) const {
+      return layerTiles_[globalBinId];
+    }
 
   private:
     std::vector< std::vector<int>> layerTiles_;

diff --git a/src/CLUEAlgo.cc b/src/CLUEAlgo.cc
@@ -1,4 +1,5 @@
 #include "CLUEAlgo.h"
+#include <array>
 
 template <typename TILES>
 void CLUEAlgo_T<TILES>::makeClusters(){
@@ -54,24 +55,26 @@ template <typename TILES>
 void CLUEAlgo_T<TILES>::calculateLocalDensity( TILES & allLayerTiles ){
 
   std::array<int,4> search_box = {0, 0, 0, 0};
+  auto dc2 = dc_*dc_;
 
   // loop over all points
   for(unsigned i = 0; i < points_.n; i++) {
-    auto lt = allLayerTiles[points_.layer[i]];
+    const auto& lt = allLayerTiles[points_.layer[i]];
     float ri = points_.r[i];
-    float phi_i = points_.x[i]/(1.*ri);
+    float inv_ri = 1.f/ri;
+    float phi_i = points_.x[i]*inv_ri;
 
     // get search box
     search_box = lt.searchBox(points_.x[i]-dc_, points_.x[i]+dc_, points_.y[i]-dc_, points_.y[i]+dc_);
 
     if(!TILES::constants_type_t::endcap){
-      float dc_phi = dc_/ri;
+      float dc_phi = dc_*inv_ri;
       search_box = lt.searchBoxPhiZ(phi_i-dc_phi, phi_i+dc_phi, points_.y[i]-dc_, points_.y[i]+dc_);
     }
 
     // loop over bins in the search box
-    for(int xBin = search_box[0]; xBin < search_box[1]+1; ++xBin) {
-      for(int yBin = search_box[2]; yBin < search_box[3]+1; ++yBin) {
+    for(int xBin = search_box[0]; xBin <= search_box[1]; ++xBin) {
+      for(int yBin = search_box[2]; yBin <= search_box[3]; ++yBin) {
 
         // get the id of this bin
         int binId = lt.getGlobalBinByBin(xBin,yBin);
@@ -86,9 +89,9 @@ void CLUEAlgo_T<TILES>::calculateLocalDensity( TILES & allLayerTiles ){
         for (unsigned int binIter = 0; binIter < binSize; binIter++) {
           int j = lt[binId][binIter];
           // query N_{dc_}(i)
-          float dist_ij = TILES::constants_type_t::endcap ?
-           distance(i, j) : distance(i, j, true, ri);
-          if(dist_ij <= dc_) {
+          float dist2_ij = TILES::constants_type_t::endcap ?
+           distance2(i, j) : distance2(i, j, true, ri);
+          if(dist2_ij <= dc2) {
             // sum weights within N_{dc_}(i)
             points_.rho[i] += (i == j ? 1.f : 0.5f) * points_.weight[j];
           }
@@ -111,19 +114,20 @@ void CLUEAlgo_T<TILES>::calculateDistanceToHigher( TILES & allLayerTiles ){
     float xi = points_.x[i];
     float yi = points_.y[i];
     float ri = points_.r[i];
-    float phi_i = points_.x[i]/(1.*ri);
+    float inv_ri = 1.f/ri;
+    float phi_i = points_.x[i]*inv_ri;
     float rho_i = points_.rho[i];
 
     //get search box
-    auto lt = allLayerTiles[points_.layer[i]];
-    float dm_phi = dm/ri;
+    const auto& lt = allLayerTiles[points_.layer[i]];
+    float dm_phi = dm*inv_ri;
     std::array<int,4> search_box = TILES::constants_type_t::endcap ? 
      lt.searchBox(xi-dm, xi+dm, yi-dm, yi+dm):
      lt.searchBoxPhiZ(phi_i-dm_phi, phi_i+dm_phi, points_.y[i]-dm, points_.y[i]+dm);
 
     // loop over all bins in the search box
-    for(int xBin = search_box[0]; xBin < search_box[1]+1; ++xBin) {
-      for(int yBin = search_box[2]; yBin < search_box[3]+1; ++yBin) {
+    for(int xBin = search_box[0]; xBin <= search_box[1]; ++xBin) {
+      for(int yBin = search_box[2]; yBin <= search_box[3]; ++yBin) {
 
         // get the id of this bin
         int phi = (xBin % TILES::constants_type_t::nColumnsPhi);
@@ -166,7 +170,7 @@ void CLUEAlgo_T<TILES>::findAndAssignClusters(){
 
   auto start = std::chrono::high_resolution_clock::now();
 
-  std::map<int,int> nClustersPerLayer;
+  std::array<int,TILES::constants_type_t::nLayers> nClustersPerLayer{};
 
   // find cluster seeds and outlier
   std::vector<int> localStack;
@@ -227,27 +231,39 @@ void CLUEAlgo_T<TILES>::findAndAssignClusters(){
 }
 
 template <typename TILES>
-inline float CLUEAlgo_T<TILES>::distance(int i, int j, bool isPhi, float r ) const {
+inline float CLUEAlgo_T<TILES>::distance2(int i, int j, bool isPhi, float r ) const {
 
   // 2-d distance on the layer
   if(points_.layer[i] == points_.layer[j] ) {
     if (isPhi) {
-      const float phi_i = points_.x[i]/(1.*points_.r[i]);
-      const float phi_j = points_.x[j]/(1.*points_.r[j]);
+      const float phi_i = points_.x[i]/(points_.r[i]);
+      const float phi_j = points_.x[j]/(points_.r[j]);
       const float drphi = r * reco::deltaPhi(phi_i, phi_j);
       const float dy = points_.y[i] - points_.y[j];
-      return std::sqrt(dy * dy + drphi * drphi);
+      return dy * dy + drphi * drphi;
     } else {
       const float dx = points_.x[i] - points_.x[j];
       const float dy = points_.y[i] - points_.y[j];
-      return std::sqrt(dx * dx + dy * dy);
+      return dx * dx + dy * dy;
     }
   } else {
     return std::numeric_limits<float>::max();
   }
 
 }
 
+template <typename TILES>
+inline float CLUEAlgo_T<TILES>::distance(int i, int j, bool isPhi, float r ) const {
+
+  // 2-d distance on the layer
+  if(points_.layer[i] == points_.layer[j] ) {
+    return std::sqrt(distance2(i, j, isPhi, r));
+  } else {
+    return std::numeric_limits<float>::max();
+  }
+
+}
+
 // explicit template instantiation
 template class CLUEAlgo_T<LayerTiles>;
 template class CLUEAlgo_T<CLICdetEndcapLayerTiles>;