Unit 1: Micro-optimizations (stack alloc + So3Y const ref)

include/robin/core.hpp

@@ -48,6 +48,11 @@ template <typename P> class SetBase {
    * @brief Return the ith point in the set
    * @param i
    * @return
+   *
+   * Returns by value intentionally: matrix-backed types (VectorY, Points3d,
+   * Points3dWithNormals) return Eigen column expressions which are temporaries
+   * and cannot be returned by reference. Types backed by std::vector (e.g. So3Y)
+   * can provide a separate const-ref accessor (ref()) for zero-copy access.
    */
   virtual P operator[](size_t i) const = 0;
 };
@@ -95,8 +100,7 @@ class BaseCompGraphConstructor {
     size_t nr_subsets = robin::Choose(N, CompCheckFuncArity);
 
     // for each possible subset, check compatibility
-    auto* subset_indices = new size_t[CompCheckFuncArity];
-    auto* edge_pair_indices = new size_t[2];
+    size_t subset_indices[CompCheckFuncArity];
     for (size_t i = 0; i < nr_subsets; ++i) {
       // compute current subset indices
       robin::CombinationDecode(N, CompCheckFuncArity, i, subset_indices);
@@ -114,9 +118,6 @@ class BaseCompGraphConstructor {
         }
       }
     }
-
-    delete[] edge_pair_indices;
-    delete[] subset_indices;
   }
 
   /**
@@ -139,8 +140,7 @@ class BaseCompGraphConstructor {
     // for each possible subset, check compatibility
 #pragma omp parallel default(none) shared(N, nr_subsets, edge_queue)
     {
-      auto* subset_indices = new size_t[CompCheckFuncArity];
-      auto* edge_pair_indices = new size_t[2];
+      size_t subset_indices[CompCheckFuncArity];
 #pragma omp for
       for (size_t i = 0; i < nr_subsets; ++i) {
         // compute current subset indices
@@ -165,9 +165,6 @@ class BaseCompGraphConstructor {
           }
         }
       }
-
-      delete[] edge_pair_indices;
-      delete[] subset_indices;
     };
 
     // add edges to graph
@@ -217,7 +214,7 @@ class CompGraphConstructor<Measurements, CompCheckFunc, 2>
     // 2. a critical section for dequeuing edges into graph
 #pragma omp parallel default(none) shared(g, N)
     {
-      auto* subset_indices = new size_t[2];
+      size_t subset_indices[2];
       std::queue<EdgeType> edge_buffer;
 
 #pragma omp for
@@ -243,7 +240,6 @@ class CompGraphConstructor<Measurements, CompCheckFunc, 2>
           edge_buffer.pop();
         }
       };
-      delete[] subset_indices;
     }
   }
 
@@ -255,7 +251,7 @@ class CompGraphConstructor<Measurements, CompCheckFunc, 2>
     // parallel compatibility graph building
     // 1. each thread initialize a queue for storing edges
     // 2. a critical section for dequeuing edges into graph
-    auto* subset_indices = new size_t[2];
+    size_t subset_indices[2];
 
     for (size_t k = 0; k < N * (N - 1) / 2; ++k) {
       size_t i = k / N;
@@ -271,8 +267,6 @@ class CompGraphConstructor<Measurements, CompCheckFunc, 2>
         g->AddEdgeUnsafe(i, j);
       }
     }
-
-    delete[] subset_indices;
   }
 
   /**
@@ -289,7 +283,7 @@ class CompGraphConstructor<Measurements, CompCheckFunc, 2>
     // 2. a critical section for dequeuing edges into graph
 #pragma omp parallel default(none) shared(g, N, adj_list, edge_count)
     {
-      auto* subset_indices = new size_t[2];
+      size_t subset_indices[2];
 
 #pragma omp for
       // visit each edge twice
@@ -304,8 +298,6 @@ class CompGraphConstructor<Measurements, CompCheckFunc, 2>
           }
         }
       }
-
-      delete[] subset_indices;
     }
 
     g->Clear();
@@ -328,7 +320,7 @@ class CompGraphConstructor<Measurements, CompCheckFunc, 2>
 #pragma omp parallel default(none) shared(g, N, offsets)
     {
       // subset_indices: a buffer for each thread
-      auto* subset_indices = new size_t[2];
+      size_t subset_indices[2];
 
 #pragma omp for
       // for loop that goes through each (i,j) edge twice
@@ -343,8 +335,6 @@ class CompGraphConstructor<Measurements, CompCheckFunc, 2>
           }
         }
       }
-
-      delete[] subset_indices;
     }
 
     // prefix sum to convert vertex degrees to offsets array for CSR graph
@@ -362,7 +352,7 @@ class CompGraphConstructor<Measurements, CompCheckFunc, 2>
 #pragma omp parallel default(none) shared(g, N, edges, offsets, vertex_local_offsets)
     {
       // subset_indices: a buffer for each thread
-      auto* subset_indices = new size_t[2];
+      size_t subset_indices[2];
 
 #pragma omp for
       // for loop that goes through each (i,j) edge twice
@@ -379,8 +369,6 @@ class CompGraphConstructor<Measurements, CompCheckFunc, 2>
           }
         }
       }
-
-      delete[] subset_indices;
     }
 
     // now the edge & offsets array are ready to be moved into a CSR graph
@@ -408,7 +396,7 @@ class CompGraphConstructor<Measurements, CompCheckFunc, 2>
 #pragma omp parallel default(none) shared(g, N, offsets)
     {
       // subset_indices: a buffer for each thread
-      auto* subset_indices = new size_t[2];
+      size_t subset_indices[2];
 
 #pragma omp for
       for (size_t k = 0; k < N * (N - 1) / 2; ++k) {
@@ -426,8 +414,6 @@ class CompGraphConstructor<Measurements, CompCheckFunc, 2>
           offsets[j + 1]++;
         }
       }
-
-      delete[] subset_indices;
     }
 
     // prefix sum to convert vertex degrees to offsets array for CSR graph
@@ -448,7 +434,7 @@ class CompGraphConstructor<Measurements, CompCheckFunc, 2>
 #pragma omp parallel default(none) shared(g, N, edges, offsets, vertex_local_offsets)
     {
       // subset_indices: a buffer for each thread
-      auto* subset_indices = new size_t[2];
+      size_t subset_indices[2];
 
 #pragma omp for
       for (size_t k = 0; k < N * (N - 1) / 2; ++k) {
@@ -470,10 +456,9 @@ class CompGraphConstructor<Measurements, CompCheckFunc, 2>
           // update edge from side of vertex j
           // put the edge in the edge array & increment local edge index offset
           auto idx_j = vertex_local_offsets[j]++;
-          edges[offsets[j].load() + idx_j].store(i);        }
+          edges[offsets[j].load() + idx_j].store(i);
+        }
       }
-
-      delete[] subset_indices;
     }
 
     // now the edge & offsets array are ready to be moved into a CSR graph

include/robin/problems.hpp

@@ -80,6 +80,8 @@ class So3Y : public robin::SetBase<Eigen::Matrix3d> {
   So3Y() = delete;
   explicit So3Y(std::vector<Eigen::Matrix3d> y_set) : y_set_(std::move(y_set)) {}
   Eigen::Matrix3d operator[](size_t i) const override { return y_set_[i]; }
+  /// Zero-copy const ref access. Safe because So3Y is backed by std::vector (no temporary).
+  const Eigen::Matrix3d& ref(size_t i) const { return y_set_[i]; }
   size_t size() const override { return y_set_.size(); }
 
 private:
@@ -132,7 +134,7 @@ template <typename DistFn> struct SraCompCheck {
    */
   bool operator()(So3Y* measurements, const size_t* subset_indices) {
     double dist =
-        (*dist_fn)((*measurements)[subset_indices[0]], (*measurements)[subset_indices[1]]);
+        (*dist_fn)(measurements->ref(subset_indices[0]), measurements->ref(subset_indices[1]));
     return dist <= noise_threshold;
   }
 

tests/core_test.cpp

@@ -283,6 +283,65 @@ TEST_CASE("large compatibility graph adj list") {
   std::cout << "Core time   (ms): " << core_finding_time / trials << std::endl;
 }
 
+TEST_CASE("graph construction methods produce same edge sets") {
+  // Validates that stack-allocation refactor doesn't break any construction path.
+  // Uses a small deterministic problem so results are reproducible.
+  robin::test::VecAveragingProblemFixture fixture;
+  Eigen::Vector2d exp_vector;
+  exp_vector << 2, 2;
+  double noise_bound = 0.1;
+  size_t N = 5;
+  auto measurements = fixture.GenerateMeasurements(exp_vector, N, 0, noise_bound);
+  robin::VectorY measurements_set(measurements);
+  robin::SvaCompCheck sva_comp_check(2 * noise_bound);
+
+  robin::CompGraphConstructor<robin::VectorY, robin::SvaCompCheck, 2> graph_constructor;
+  graph_constructor.SetCompCheckFunction(&sva_comp_check);
+  graph_constructor.SetMeasurements(&measurements_set);
+
+  // Helper: extract sorted edge set from an AdjListGraph
+  auto extract_edges = [](robin::AdjListGraph& g) {
+    std::set<std::pair<size_t, size_t>> edges;
+    for (size_t i = 0; i < g.VertexCount(); ++i) {
+      for (size_t k = 0; k < g.GetVertexDegree(i); ++k) {
+        size_t j = g.GetVertexEdge(i, k);
+        // store each edge with smaller index first
+        edges.insert(i < j ? std::make_pair(i, j) : std::make_pair(j, i));
+      }
+    }
+    return edges;
+  };
+
+  SECTION("serial vs parallel_edge_buffer vs parallel_vertex") {
+    robin::AdjListGraph g_serial;
+    graph_constructor.BuildCompGraph_serial(&g_serial);
+    auto edges_serial = extract_edges(g_serial);
+
+    robin::AdjListGraph g_edge_buf;
+    graph_constructor.BuildCompGraph_parallel_edge_buffer(&g_edge_buf);
+    auto edges_edge_buf = extract_edges(g_edge_buf);
+
+    robin::AdjListGraph g_vertex;
+    graph_constructor.BuildCompGraph_parallel_vertex(&g_vertex);
+    auto edges_vertex = extract_edges(g_vertex);
+
+    // All methods must produce the same vertex and edge counts
+    REQUIRE(g_serial.VertexCount() == N);
+    REQUIRE(g_edge_buf.VertexCount() == N);
+    REQUIRE(g_vertex.VertexCount() == N);
+
+    REQUIRE(edges_serial.size() == edges_edge_buf.size());
+    REQUIRE(edges_serial.size() == edges_vertex.size());
+
+    // All methods must produce identical edge sets
+    REQUIRE(edges_serial == edges_edge_buf);
+    REQUIRE(edges_serial == edges_vertex);
+
+    // For N=5 with no outliers and tight noise bound, expect complete graph
+    REQUIRE(edges_serial.size() == N * (N - 1) / 2);
+  }
+}
+
 TEST_CASE("sample comp graph construction") {
   SECTION("2d vector averaging") {
     //

tests/problems_test.cpp

@@ -16,5 +16,22 @@
 
 #include <robin/core.hpp>
 #include <robin/graph.hpp>
+#include <robin/problems.hpp>
 #include <robin/utils.hpp>
 
+TEST_CASE("So3Y ref() returns const reference") {
+  std::vector<Eigen::Matrix3d> rotations = {Eigen::Matrix3d::Identity(),
+                                             Eigen::Matrix3d::Identity() * 2};
+  robin::So3Y y_set(rotations);
+  const Eigen::Matrix3d& ref1 = y_set.ref(0);
+  const Eigen::Matrix3d& ref2 = y_set.ref(0);
+  REQUIRE(&ref1 == &ref2); // Same address -- true reference, not temporary
+}
+
+TEST_CASE("So3Y operator[] still works by value") {
+  std::vector<Eigen::Matrix3d> rotations = {Eigen::Matrix3d::Identity()};
+  robin::So3Y y_set(rotations);
+  Eigen::Matrix3d val = y_set[0];
+  REQUIRE(val.isApprox(Eigen::Matrix3d::Identity()));
+}
+