diff --git a/include/robin/core.hpp b/include/robin/core.hpp index 5477d1d..674c83a 100644 --- a/include/robin/core.hpp +++ b/include/robin/core.hpp @@ -48,6 +48,11 @@ template 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 // 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 edge_buffer; #pragma omp for @@ -243,7 +240,6 @@ class CompGraphConstructor edge_buffer.pop(); } }; - delete[] subset_indices; } } @@ -255,7 +251,7 @@ class CompGraphConstructor // 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 g->AddEdgeUnsafe(i, j); } } - - delete[] subset_indices; } /** @@ -289,7 +283,7 @@ class CompGraphConstructor // 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 } } } - - delete[] subset_indices; } g->Clear(); @@ -328,7 +320,7 @@ class CompGraphConstructor #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 } } } - - delete[] subset_indices; } // prefix sum to convert vertex degrees to offsets array for CSR graph @@ -362,7 +352,7 @@ class CompGraphConstructor #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 } } } - - delete[] subset_indices; } // now the edge & offsets array are ready to be moved into a CSR graph @@ -408,7 +396,7 @@ class CompGraphConstructor #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 offsets[j + 1]++; } } - - delete[] subset_indices; } // prefix sum to convert vertex degrees to offsets array for CSR graph @@ -448,7 +434,7 @@ class CompGraphConstructor #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 // 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 diff --git a/include/robin/problems.hpp b/include/robin/problems.hpp index 926ac1c..a18ec00 100644 --- a/include/robin/problems.hpp +++ b/include/robin/problems.hpp @@ -80,6 +80,8 @@ class So3Y : public robin::SetBase { So3Y() = delete; explicit So3Y(std::vector 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 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; } diff --git a/tests/core_test.cpp b/tests/core_test.cpp index 08614e1..f63e22a 100644 --- a/tests/core_test.cpp +++ b/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 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> 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") { // diff --git a/tests/problems_test.cpp b/tests/problems_test.cpp index 2afb322..3a9b76d 100644 --- a/tests/problems_test.cpp +++ b/tests/problems_test.cpp @@ -16,5 +16,22 @@ #include #include +#include #include +TEST_CASE("So3Y ref() returns const reference") { + std::vector 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 rotations = {Eigen::Matrix3d::Identity()}; + robin::So3Y y_set(rotations); + Eigen::Matrix3d val = y_set[0]; + REQUIRE(val.isApprox(Eigen::Matrix3d::Identity())); +} +