From bb5398870c591e4c8f37ee7c1d918cde6dd05def Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Tue, 21 Jul 2026 07:52:54 -0700 Subject: [PATCH 1/3] fix(moves): enforce causal ergodic move contracts - reject malformed cavities, invalid time assignments, and stale metadata - preserve source manifolds and rebuild handle-bearing state after mutation - validate complete CDT invariants and literature-derived move deltas - document admissibility, inverse relationships, and failure atomicity - update the repository UV pin to 0.11.30 Closes #106 --- Justfile | 2 +- README.md | 3 + docs/ergodic-moves.md | 103 ++++ include/Ergodic_moves_3.hpp | 191 ++++-- include/Metropolis.hpp | 3 +- include/Move_command.hpp | 3 +- tests/CMakeLists.txt | 1 + tests/Ergodic_moves_3_audit_test.cpp | 840 +++++++++++++++++++++++++++ 8 files changed, 1098 insertions(+), 48 deletions(-) create mode 100644 docs/ergodic-moves.md create mode 100644 tests/Ergodic_moves_3_audit_test.cpp diff --git a/Justfile b/Justfile index aefe49247..27eb2fbf7 100644 --- a/Justfile +++ b/Justfile @@ -5,7 +5,7 @@ set minimum-version := "1.57.0" set shell := ["bash", "-euo", "pipefail", "-c"] just_version := "1.57.0" -uv_version := "0.11.29" +uv_version := "0.11.30" pinact_version := "4.1.0" pinact_module := "github.com/suzuki-shunsuke/pinact/v4/cmd/pinact@v" + pinact_version llvm_version := "22" diff --git a/README.md b/README.md index 6e55efc67..91264b014 100644 --- a/README.md +++ b/README.md @@ -370,6 +370,9 @@ Online documentation is at . The scientific transition, proposal-ratio, geometry-delta, counter, and precision contracts are recorded in [`docs/metropolis-hastings.md`](docs/metropolis-hastings.md). +The literature-backed contracts, exact deltas, inverse relationships, and +failure-atomicity rules for the complete 2+1D move set are recorded in +[`docs/ergodic-moves.md`](docs/ergodic-moves.md). Seed replay, PCG stream ownership, checkpoint metadata, and the parallel stream policy are recorded in [`docs/reproducibility.md`](docs/reproducibility.md). The repository-wide scientific bibliography is maintained in diff --git a/docs/ergodic-moves.md b/docs/ergodic-moves.md new file mode 100644 index 000000000..2d6917059 --- /dev/null +++ b/docs/ergodic-moves.md @@ -0,0 +1,103 @@ +# 2+1D CDT ergodic move audit + +This note records the scientific and implementation contract for the five +supported 2+1-dimensional causal dynamical triangulation (CDT) moves. It is the +audit record for issue #106 and distinguishes three separate claims: + +1. the local replacement is one of the Lorentzian Monte Carlo moves; +2. CGAL preserves the validity of the combinatorial triangulation data + structure; and +3. CDT++ independently verifies the stronger causal, foliation, metadata, and + failure-atomicity invariants. + +The primary CDT source is Ambjørn, Jurkiewicz, and Loll, *Dynamically +triangulating Lorentzian quantum gravity*, section 7.1, equations (56)-(58), +with the corresponding f-vector changes in equations (71)-(73). See +@ref cdt-framework-2001 "CDT framework (2001)". The topological background is +@ref pachner-moves "Pachner's theorem", and the implementation primitive +contracts come from the @ref cgal-triangulations "CGAL 6.2 manual". + +## Conventions and common contract + +A `(p,q)` tetrahedron has `p` vertices on slice `t` and `q` vertices on slice +`t+1`. Every accepted finite tetrahedron must therefore have exactly two +adjacent time values and type `(3,1)`, `(2,2)`, or `(1,3)`. A spacelike +subsimplex has equal endpoint times; a timelike one spans adjacent slices. + +Every high-level `do_*_move()` and `propose_*_move()` operation mutates an +owning snapshot. An inapplicable site, a checked CGAL rejection, or a failed +postcondition leaves the source manifold's canonical topology, vertex and cell +metadata, geometry counts, time bounds, foliation parameters, and handle caches +unchanged. A successful move constructs a new `FoliatedTriangulation`, which +rebuilds all handle-bearing caches and scalar geometry from the moved canonical +triangulation. + +CGAL 6.2 documents that checked three-dimensional flips preserve vertex handles +and invalidate only affected cell handles. CDT++ never carries affected cell +handles past a successful flip. A copied triangulation has different handles; +the `(4,4)` and `(6,2)` implementations therefore capture stable point values +and re-resolve vertices and edges in their private copies before mutation. + +None of these moves is required to preserve the Euclidean empty-sphere +(Delaunay) property of the representative coordinates. The scientific state is +a valid causal combinatorial triangulation. `tds().is_valid()` establishes the +CGAL incidence and adjacency contract; `Manifold::is_correct()` adds foliation +and causal cell-metadata checks. + +## Move-by-move record + +| Move | Local cavity and time assignment | Independent delta `(N0, N1_SL, N1_TL, N2, N3_31, N3_22, N3_13, N3)` | Implementation and admissibility | +| --- | --- | --- | --- | +| `(2,3)` | Equation (58): `(3,1) 1345 + (2,2) 2345 -> (3,1) 1234 + (2,2) 1235 + (2,2) 1245`, or its time reflection. The shared triangle `345` is timelike and the new edge `12` is timelike. | `(0, 0, +1, +2, 0, +1, 0, +1)` | A checked `Triangulation_3::flip(facet)` is attempted only from correctly labelled `(2,2)` cells whose neighbor is a correctly labelled `(3,1)` or `(1,3)` cell and whose opposite vertices lie on adjacent slices. CGAL additionally rejects infinite, nonflippable, or geometrically inverted cavities. | +| `(3,2)` | The inverse of equation (58): two `(2,2)` cells and exactly one `(3,1)` or `(1,3)` cell meet at the timelike edge removed by the move. | `(0, 0, -1, -2, 0, -1, 0, -1)` | Before `Triangulation_3::flip(edge)`, CDT++ independently requires three finite incident cells spanning adjacent slices with the exact `2 x (2,2) + 1 x ((3,1) or (1,3))` causal composition. CGAL then rejects hull edges and enforces the total degree-three and geometric flippability contracts. Other degree-three timelike cavities can be combinatorially flippable but causally invalid. | +| `(2,6)` | Equation (56): `(1,3) 1345 + (3,1) 2345` share the spacelike triangle `345`. A new vertex `6` is inserted on the same slice and joined to all five old vertices, producing three tetrahedra above and three below. | `(+1, +3, +2, +8, +2, 0, +2, +4)` | `tds().insert_in_facet()` subdivides the common spacelike facet. The input cells and their metadata must be `(1,3)` and `(3,1)`, the shared vertices must have one time value, and the new star must contain six valid cells. The new point is the facet centroid and receives the facet time. | +| `(6,2)` | The inverse of equation (56): a degree-five vertex has six incident cells, exactly three `(3,1)` and three `(1,3)`, with no `(2,2)` cell. | `(-1, -3, -2, -8, -2, 0, -2, -4)` | On a private copy, a checked timelike edge flip reduces the candidate to degree four; `tds().remove_from_maximal_dimension_simplex()` then applies its documented degree-`dimension+1` removal. Exact finite incidence, causal types, metadata, output counts, and output cell types are checked before publication. | +| `(4,4)` | Equation (57): two `(1,3)` and two `(3,1)` tetrahedra form a diamond. The diagonal of the spatial quadrilateral is exchanged; the move is its own inverse. | `(0, 0, 0, 0, 0, 0, 0, 0)` | The pivot must be a spacelike edge with exactly four finite incident cells, two `(3,1)` and two `(1,3)`, all correctly labelled. A checked TDS facet flip creates the new diagonal, then a checked TDS edge flip removes the old one. The composition is performed on a private copy because the transient geometry can fail a `Triangulation_3` geometric flip even when the final abstract diamond is valid. | + +## Independent delta derivations + +For `(2,3)`, two tetrahedra sharing one face initially have five vertices, nine +edges, seven faces, and two cells. Replacing the shared face by the dual edge +adds that one timelike edge. The three resulting tetrahedra have nine distinct +faces, so `(N1, N2, N3)` changes by `(+1, +2, +1)`. Equation (58) changes one +`(3,1)` plus one `(2,2)` into one `(3,1)` plus two `(2,2)`. + +For `(2,6)`, the inserted vertex is joined to all five old vertices. Its three +edges to the shared spatial triangle are spacelike and its two edges to the +opposite vertices are timelike. Thus `N0` changes by `+1` and `N1` by `+5`. +The local complex changes from two to six tetrahedra. With the boundary fixed, +the Euler relation for the three-ball gives +`delta N2 = delta N1 + delta N3 - delta N0 = 5 + 4 - 1 = 8`. Equation (56) +changes one tetrahedron of each time orientation into three of each, giving +`delta N3_31 = delta N3_13 = +2`. The `(6,2)` changes are the negatives. + +For `(4,4)`, the six vertices and the boundary of the four-cell diamond are +unchanged. One spacelike diagonal replaces another, four tetrahedra replace +four tetrahedra, and the two tetrahedra of each time orientation remain two. +The fixed boundary and Euler relation then force every tracked count delta to +zero. + +These derivations agree with the source's f-vectors (71)-(73) but do not call +`ergodic_moves::check_move()` or any production count-delta helper. + +## Verification record + +`tests/Ergodic_moves_3_audit_test.cpp` supplies: + +- deterministic minimal successful fixtures for all five moves; +- direct enumeration of finite vertices, edges, facets, cells, edge types, and + tetrahedron types from vertex time labels; +- reciprocal adjacency, Euler, TDS, foliation, time-bound, metadata, and cache + checks after every successful move; +- canonical point/time/simplex representations independent of CGAL handle + identity; +- exact inverse round trips for `(2,3)/(3,2)`, `(2,6)/(6,2)`, and `(4,4)`; +- malformed-handle, stale-metadata, wrong-cavity, non-applicable, and empty-state + rejection checks with canonical failure-atomicity comparisons; and +- replay of every inverse pair with seeds `0`, `1`, `2`, `92`, `106`, and + `20260721`. + +The wrong-cavity `(3,2)` regression deliberately constructs a geometrically +flippable CGAL edge whose incident CDT types are invalid. It verifies that raw +CGAL accepts the flip while CDT++ rejects before mutation, directly testing the +boundary between combinatorial and scientific validity. diff --git a/include/Ergodic_moves_3.hpp b/include/Ergodic_moves_3.hpp index 771057c99..f28c49bfe 100644 --- a/include/Ergodic_moves_3.hpp +++ b/include/Ergodic_moves_3.hpp @@ -14,12 +14,11 @@ /// @see [Pachner moves](../REFERENCES.md#pachner-moves) /// @see [Three-dimensional CDT move /// set](../REFERENCES.md#three-dimensional-cdt-2001) +/// @see [2+1D CDT ergodic move audit](../docs/ergodic-moves.md) #ifndef CDT_PLUSPLUS_ERGODIC_MOVES_3_HPP #define CDT_PLUSPLUS_ERGODIC_MOVES_3_HPP -#include - #include #include @@ -65,6 +64,32 @@ namespace ergodic_moves valid_index(edge.third) && edge.second != edge.third; } + /// @brief Collect the finite cells incident to a checked edge. + /// @details The returned handles borrow from @p triangulation and remain + /// valid only until an affected-cell mutation is performed. + [[nodiscard]] inline auto finite_incident_cells( + Delaunay const& triangulation, Edge_handle const& edge) + -> std::optional + { + if (!is_well_formed_edge(edge) || triangulation.dimension() != 3 || + !triangulation.tds().is_edge(edge.first, edge.second, edge.third)) + { + return std::nullopt; + } + + auto circulator = triangulation.incident_cells(edge, edge.first); + Cell_container incident_cells; + do + { // NOLINT(cppcoreguidelines-avoid-do-while) + if (!triangulation.is_infinite(circulator)) + { + incident_cells.emplace_back(circulator); + } + } + while (++circulator != edge.first); + return incident_cells; + } + /// Select exactly one raw proposal site uniformly. template [[nodiscard]] inline auto random_element(Container const& candidates, @@ -85,7 +110,9 @@ namespace ergodic_moves [[nodiscard]] inline auto null_move(Manifold const& t_manifold) -> Expected { return t_manifold; } // null_move - /// @brief Perform a TriangulationDataStructure_3::flip on a facet + /// @brief Perform a checked Triangulation_3 flip on a causal CDT facet + /// @details On success, affected cell handles are invalidated; vertex handles + /// remain valid. Rejection occurs before mutation. /// @param triangulation The triangulation containing the cell to flip /// @param to_be_moved The cell on which to try the move /// @returns True if move succeeded @@ -94,7 +121,14 @@ namespace ergodic_moves [[nodiscard]] inline auto try_23_move(Delaunay& triangulation, Cell_handle const& to_be_moved) -> bool { - if (to_be_moved->info() != 22) { return false; } // NOLINT + if (to_be_moved == nullptr || triangulation.dimension() != 3 || + !triangulation.tds().is_cell(to_be_moved) || + !foliated_triangulations::is_cell_type_correct<3>(to_be_moved) || + foliated_triangulations::expected_cell_type<3>(to_be_moved) != + Cell_type::TWO_TWO) + { + return false; + } auto flipped = false; // Try every facet of the (2,2) cell for (auto i = 0; i < 4; ++i) @@ -102,6 +136,15 @@ namespace ergodic_moves auto const neighbor = to_be_moved->neighbor(i); if (triangulation.is_infinite(neighbor)) { continue; } + auto const neighbor_type = + foliated_triangulations::expected_cell_type<3>(neighbor); + if (!foliated_triangulations::is_cell_type_correct<3>(neighbor) || + (neighbor_type != Cell_type::THREE_ONE && + neighbor_type != Cell_type::ONE_THREE)) + { + continue; + } + // A causal (2,3) move must replace the facet with a timelike edge. // CGAL also permits topological flips that create a spacelike edge. auto const mirror_index = neighbor->index(to_be_moved); @@ -131,13 +174,15 @@ namespace ergodic_moves /// This function calls try_23_move on (2,2) simplices drawn from a /// randomly shuffled container until it succeeds or runs out of simplices. /// - /// If successful, the triangulation is no longer Delaunay. + /// The move guarantees a valid causal combinatorial triangulation. It does + /// not preserve or require the empty-sphere property of the coordinates. /// /// @tparam Generator A uniform random bit generator type /// @param t_manifold The simplicial manifold /// @param generator Caller-owned generator whose state advances during the /// move /// @returns The Expected (2,3) moved manifold or an Unexpected + /// @note The source manifold is unchanged on success and failure. template [[nodiscard]] inline auto do_23_move(Manifold const& t_manifold, Generator& generator) -> Expected @@ -180,15 +225,58 @@ namespace ergodic_moves return std::unexpected("Selected (2,3) proposal site is not movable.\n"); } - /// @brief Perform a TriangulationDataStructure_3::flip on an edge + namespace detail + { + /// @brief Check for the causal cavity inverse to a (2,3) move. + [[nodiscard]] inline auto is_32_movable(Delaunay const& triangulation, + Edge_handle const& candidate) + -> bool + { + auto const incident_cells = + finite_incident_cells(triangulation, candidate); + if (!incident_cells || incident_cells->size() != 3) { return false; } + + auto const first_time = static_cast( + candidate.first->vertex(candidate.second)->info()); + auto const second_time = static_cast( + candidate.first->vertex(candidate.third)->info()); + auto const time_difference = first_time > second_time + ? first_time - second_time + : second_time - first_time; + if (time_difference != 1) { return false; } + + if (!std::ranges::all_of(*incident_cells, [](auto const& cell) { + return foliated_triangulations::is_cell_type_correct<3>(cell); + })) + { + return false; + } + + auto const cell_type_count = [&](Cell_type const type) { + return std::ranges::count_if(*incident_cells, [&](auto const& cell) { + return foliated_triangulations::expected_cell_type<3>(cell) == type; + }); + }; + auto const incident_31 = cell_type_count(Cell_type::THREE_ONE); + auto const incident_22 = cell_type_count(Cell_type::TWO_TWO); + auto const incident_13 = cell_type_count(Cell_type::ONE_THREE); + return incident_22 == 2 && ((incident_31 == 1 && incident_13 == 0) || + (incident_31 == 0 && incident_13 == 1)); + } + } // namespace detail + + /// @brief Perform a checked (3,2) flip on a causal CDT edge cavity. + /// @details On success, affected cell handles are invalidated; vertex handles + /// remain valid. Rejection occurs before mutation. /// @param triangulation The triangulation containing the edge to flip /// @param to_be_moved The edge on which to try the move - /// @returns True if move succeeded + /// @returns True if the CDT cavity is admissible and CGAL performs the flip /// @see /// https://doc.cgal.org/latest/TDS_3/classTriangulationDataStructure__3.html#a5837d666e4198f707f862003c1ffa033 [[nodiscard]] inline auto try_32_move(Delaunay& triangulation, Edge_handle const& to_be_moved) -> bool { + if (!detail::is_32_movable(triangulation, to_be_moved)) { return false; } return triangulation.flip(to_be_moved.first, to_be_moved.second, to_be_moved.third); } // try_32_move @@ -198,12 +286,14 @@ namespace ergodic_moves /// This removes a (2,2) simplex and the timelike edge. /// This function calls try_32_move on timelike edges drawn from a /// randomly shuffled container until it succeeds or runs out of edges. - /// If successful, the triangulation is no longer Delaunay. + /// The move guarantees a valid causal combinatorial triangulation. It does + /// not preserve or require the empty-sphere property of the coordinates. /// @tparam Generator A uniform random bit generator type /// @param t_manifold The simplicial manifold /// @param generator Caller-owned generator whose state advances during the /// move /// @returns The Expected (3,2) moved manifold or an Unexpected + /// @note The source manifold is unchanged on success and failure. template [[nodiscard]] inline auto do_32_move(Manifold const& t_manifold, Generator& generator) -> Expected @@ -251,11 +341,19 @@ namespace ergodic_moves [[nodiscard]] inline auto find_adjacent_31_cell(Cell_handle const& t_cell) -> std::optional { - if (t_cell->info() != 13) { return std::nullopt; } // NOLINT + if (t_cell == nullptr || + !foliated_triangulations::is_cell_type_correct<3>(t_cell) || + foliated_triangulations::expected_cell_type<3>(t_cell) != + Cell_type::ONE_THREE) + { + return std::nullopt; + } for (auto i = 0; i < 4; ++i) { - if (foliated_triangulations::expected_cell_type<3>(t_cell->neighbor(i)) == - Cell_type::THREE_ONE) + auto const neighbor = t_cell->neighbor(i); + if (foliated_triangulations::is_cell_type_correct<3>(neighbor) && + foliated_triangulations::expected_cell_type<3>(neighbor) == + Cell_type::THREE_ONE) { return std::make_optional(i); } @@ -272,7 +370,8 @@ namespace ergodic_moves /// vertex. /// This function calls find_adjacent_31_cell on (1,3) simplices drawn from a /// randomly shuffled container until it succeeds or runs out of simplices. - /// If successful, the triangulation is no longer Delaunay. + /// The move guarantees a valid causal combinatorial triangulation. It does + /// not preserve or require the empty-sphere property of the coordinates. /// @image html 26.png /// @image latex 26.eps width=7cm /// @tparam Generator A uniform random bit generator type @@ -281,6 +380,7 @@ namespace ergodic_moves /// move /// @param only_first_site Whether to examine only one uniformly selected site /// @returns The Expected (2,6) moved manifold or an Unexpected + /// @note The source manifold is unchanged on success and failure. template [[nodiscard]] inline auto do_26_move_impl(Manifold const& t_manifold, Generator& generator, @@ -460,10 +560,13 @@ namespace ergodic_moves auto const incident_22 = cell_type_count(Cell_type::TWO_TWO); auto const incident_13 = cell_type_count(Cell_type::ONE_THREE); - // All cells should be classified - if (incident_13 + incident_22 + incident_31 != 6) // NOLINT + // All cells should be causally classified and carry matching metadata. + if (incident_13 + incident_22 + incident_31 != 6 || // NOLINT + !std::ranges::all_of(incident_cells, [](auto const& cell) { + return foliated_triangulations::is_cell_type_correct<3>(cell); + })) { - spdlog::warn("Some incident cells on this vertex need to be fixed.\n"); + return false; } return incident_31 == 3 && incident_22 == 0 && incident_13 == 3; @@ -555,14 +658,15 @@ namespace ergodic_moves /// This function calls is_62_movable() on a randomly shuffled container /// of vertices until it succeeds or runs out of vertices. /// - /// If successful, the triangulation remains Delaunay. (Other moves may - /// change this, however.) + /// The move guarantees a valid causal combinatorial triangulation. It does + /// not preserve or require the empty-sphere property of the coordinates. /// /// @tparam Generator A uniform random bit generator type /// @param t_manifold The simplicial manifold /// @param generator Caller-owned generator whose state advances during the /// move /// @returns The Expected (6,2) moved manifold or Unexpected + /// @note The source manifold is unchanged on success and failure. template [[nodiscard]] inline auto do_62_move(Manifold const& t_manifold, Generator& generator) -> Expected @@ -613,25 +717,7 @@ namespace ergodic_moves Delaunay_t<3> const& triangulation, Edge_handle const& edge) -> std::optional { - if (!detail::is_well_formed_edge(edge) || - !triangulation.tds().is_edge(edge.first, edge.second, edge.third)) - { - return std::nullopt; - } - // Create the circulator of cells around the edge, starting with the cell - // containing the edge - auto circulator = triangulation.incident_cells(edge, edge.first); - Cell_container incident_cells; - // Add cells to the container until we get back to the first one in the - // circulator - do - { // NOLINT(cppcoreguidelines-avoid-do-while) - // Ignore cells containing the infinite vertex - if (triangulation.is_infinite(circulator)) { continue; } - incident_cells.emplace_back(circulator); - } - while (++circulator != edge.first); - return incident_cells; + return detail::finite_incident_cells(triangulation, edge); } // incident_cells_from_edge() /// @brief Find a bistellar flip location @@ -646,13 +732,21 @@ namespace ergodic_moves Delaunay const& triangulation, Edge_handle const& t_edge_candidate) -> std::optional { + if (!detail::is_well_formed_edge(t_edge_candidate)) { return std::nullopt; } + auto const first_time = + t_edge_candidate.first->vertex(t_edge_candidate.second)->info(); + auto const second_time = + t_edge_candidate.first->vertex(t_edge_candidate.third)->info(); + if (first_time != second_time) { return std::nullopt; } + if (auto incident_cells = incident_cells_from_edge(triangulation, t_edge_candidate); incident_cells.has_value() && incident_cells->size() == 4) { auto const cell_type_count = [&](Cell_type const type) { return std::ranges::count_if(*incident_cells, [&](auto const cell) { - return cell->info() == static_cast(type); + return foliated_triangulations::is_cell_type_correct<3>(cell) && + foliated_triangulations::expected_cell_type<3>(cell) == type; }); }; if (cell_type_count(Cell_type::THREE_ONE) == 2 && @@ -866,6 +960,7 @@ namespace ergodic_moves /// @param generator Caller-owned generator whose state advances during the /// move /// @return The Expected (4,4) moved manifold or Unexpected + /// @note The source manifold is unchanged on success and failure. template [[nodiscard]] inline auto do_44_move(Manifold const& t_manifold, Generator& generator) -> Expected @@ -979,7 +1074,10 @@ namespace ergodic_moves "Selected (4,4) proposal site could not be flipped.\n"); } - /// @brief Check move correctness + /// @brief Check tracked move deltas and the complete CDT manifold invariant + /// @details This verifies causal foliation, cell metadata, TDS validity, + /// geometry counts, time bounds, and preserved foliation parameters. It does + /// not require Euclidean Delaunayhood. /// @param t_before The manifold before the move /// @param t_after The manifold after the move /// @param t_move The type of move @@ -989,10 +1087,17 @@ namespace ergodic_moves move_tracker::move_type const& t_move) -> bool { + if (!t_after.is_correct() || + t_after.initial_radius() != t_before.initial_radius() || + t_after.foliation_spacing() != t_before.foliation_spacing()) + { + return false; + } + switch (t_move) { case move_tracker::move_type::FOUR_FOUR: - return t_after.is_valid() && t_after.N3() == t_before.N3() && + return t_after.N3() == t_before.N3() && t_after.N3_31() == t_before.N3_31() && t_after.N3_22() == t_before.N3_22() && t_after.N3_13() == t_before.N3_13() && @@ -1003,7 +1108,7 @@ namespace ergodic_moves t_after.max_time() == t_before.max_time() && t_after.min_time() == t_before.min_time(); case move_tracker::move_type::TWO_THREE: - return t_after.is_valid() && t_after.N3() == t_before.N3() + 1 && + return t_after.N3() == t_before.N3() + 1 && t_after.N3_31() == t_before.N3_31() && t_after.N3_22() == t_before.N3_22() + 1 && t_after.N3_13() == t_before.N3_13() && @@ -1015,7 +1120,7 @@ namespace ergodic_moves t_after.max_time() == t_before.max_time() && t_after.min_time() == t_before.min_time(); case move_tracker::move_type::THREE_TWO: - return t_after.is_valid() && t_after.N3() == t_before.N3() - 1 && + return t_after.N3() == t_before.N3() - 1 && t_after.N3_31() == t_before.N3_31() && t_after.N3_22() == t_before.N3_22() - 1 && t_after.N3_13() == t_before.N3_13() && @@ -1027,7 +1132,7 @@ namespace ergodic_moves t_after.max_time() == t_before.max_time() && t_after.min_time() == t_before.min_time(); case move_tracker::move_type::TWO_SIX: - return t_after.is_valid() && t_after.N3() == t_before.N3() + 4 && + return t_after.N3() == t_before.N3() + 4 && t_after.N3_31() == t_before.N3_31() + 2 && t_after.N3_22() == t_before.N3_22() && t_after.N3_13() == t_before.N3_13() + 2 && @@ -1039,7 +1144,7 @@ namespace ergodic_moves t_after.max_time() == t_before.max_time() && t_after.min_time() == t_before.min_time(); case move_tracker::move_type::SIX_TWO: - return t_after.is_valid() && t_after.N3() == t_before.N3() - 4 && + return t_after.N3() == t_before.N3() - 4 && t_after.N3_31() == t_before.N3_31() - 2 && t_after.N3_22() == t_before.N3_22() && t_after.N3_13() == t_before.N3_13() - 2 && diff --git a/include/Metropolis.hpp b/include/Metropolis.hpp index cfa9b5ea3..cfa2d6894 100644 --- a/include/Metropolis.hpp +++ b/include/Metropolis.hpp @@ -344,8 +344,7 @@ class MoveStrategy ++m_attempted_moves[move]; auto candidate = propose_candidate(current, move); - if (!candidate || !candidate->is_correct() || - !ergodic_moves::check_move(current, *candidate, move)) + if (!candidate || !ergodic_moves::check_move(current, *candidate, move)) { ++m_failed_moves[move]; ++m_rejected_moves[move]; diff --git a/include/Move_command.hpp b/include/Move_command.hpp index 69dd92b1e..6c211e02f 100644 --- a/include/Move_command.hpp +++ b/include/Move_command.hpp @@ -132,8 +132,7 @@ class MoveCommand generator); result) { - if (result->is_correct() && - ergodic_moves::check_move(m_manifold, *result, move_type)) + if (ergodic_moves::check_move(m_manifold, *result, move_type)) { swap(result.value(), m_manifold); ++m_succeeded[as_integer(move_type)]; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index d9f64b528..b50eb472c 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -3,6 +3,7 @@ add_executable( ${PROJECT_SOURCE_DIR}/tests/main.cpp Apply_move_test.cpp Bistellar_flip_test.cpp + Ergodic_moves_3_audit_test.cpp Ergodic_moves_3_test.cpp Foliated_triangulation_test.cpp Function_ref_test.cpp diff --git a/tests/Ergodic_moves_3_audit_test.cpp b/tests/Ergodic_moves_3_audit_test.cpp new file mode 100644 index 000000000..6cbbc9652 --- /dev/null +++ b/tests/Ergodic_moves_3_audit_test.cpp @@ -0,0 +1,840 @@ +/******************************************************************************* + Causal Dynamical Triangulations in C++ using CGAL + + Copyright © 2026 Adam Getchell + ******************************************************************************/ + +/// @file Ergodic_moves_3_audit_test.cpp +/// @brief Independent scientific audit tests for the 2+1D CDT move set + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Ergodic_moves_3.hpp" + +namespace +{ + using Manifold = manifolds::Manifold_3; + using Delaunay = ergodic_moves::Delaunay; + using Cell_handle = ergodic_moves::Cell_handle; + using Edge_handle = ergodic_moves::Edge_handle; + using Vertex_handle = ergodic_moves::Vertex_handle; + + using Vertex_key = std::tuple; + + template + using Simplex_key = std::array; + + using Cell_record = std::pair, Int_precision>; + + struct Direct_counts + { + Int_precision n0{}; + Int_precision n1_sl{}; + Int_precision n1_tl{}; + Int_precision n2{}; + Int_precision n3_31{}; + Int_precision n3_22{}; + Int_precision n3_13{}; + Int_precision n3{}; + + auto operator==(Direct_counts const&) const -> bool = default; + }; + + struct Count_delta + { + Int_precision n0{}; + Int_precision n1_sl{}; + Int_precision n1_tl{}; + Int_precision n2{}; + Int_precision n3_31{}; + Int_precision n3_22{}; + Int_precision n3_13{}; + Int_precision n3{}; + }; + + struct Canonical_state + { + std::vector vertices; + std::vector> edges; + std::vector> facets; + std::vector cells; + Direct_counts counts; + std::vector> spacelike_facets; + Int_precision min_time{}; + Int_precision max_time{}; + double initial_radius{}; + double foliation_spacing{}; + + auto operator==(Canonical_state const&) const -> bool = default; + }; + + [[nodiscard]] auto vertex_key(Vertex_handle const vertex) -> Vertex_key + { + return {CGAL::to_double(vertex->point().x()), + CGAL::to_double(vertex->point().y()), + CGAL::to_double(vertex->point().z()), vertex->info()}; + } + + template + [[nodiscard]] auto sorted_simplex( + std::array const& vertices) -> Simplex_key + { + Simplex_key result{}; + std::ranges::transform(vertices, result.begin(), vertex_key); + std::ranges::sort(result); + return result; + } + + [[nodiscard]] auto independent_cell_type(Cell_handle const cell) -> Cell_type + { + std::array times{}; + for (int index = 0; index < 4; ++index) + { + times.at(static_cast(index)) = cell->vertex(index)->info(); + } + auto const [minimum, maximum] = std::ranges::minmax(times); + if (maximum - minimum != 1) { return Cell_type::ACAUSAL; } + + auto const lower_count = std::ranges::count(times, minimum); + if (lower_count == 3) { return Cell_type::THREE_ONE; } + if (lower_count == 2) { return Cell_type::TWO_TWO; } + if (lower_count == 1) { return Cell_type::ONE_THREE; } + return Cell_type::UNCLASSIFIED; + } + + void assign_independent_cell_metadata(Delaunay& triangulation) + { + for (auto cell = triangulation.finite_cells_begin(); + cell != triangulation.finite_cells_end(); ++cell) + { + cell->info() = static_cast(independent_cell_type(cell)); + } + } + + [[nodiscard]] auto direct_counts(Delaunay const& triangulation) + -> Direct_counts + { + Direct_counts result{}; + result.n0 = static_cast(triangulation.number_of_vertices()); + result.n2 = + static_cast(triangulation.number_of_finite_facets()); + result.n3 = + static_cast(triangulation.number_of_finite_cells()); + + for (auto edge = triangulation.finite_edges_begin(); + edge != triangulation.finite_edges_end(); ++edge) + { + auto const first_time = edge->first->vertex(edge->second)->info(); + auto const second_time = edge->first->vertex(edge->third)->info(); + if (first_time == second_time) { ++result.n1_sl; } + else + { + ++result.n1_tl; + } + } + + for (auto cell = triangulation.finite_cells_begin(); + cell != triangulation.finite_cells_end(); ++cell) + { + switch (independent_cell_type(cell)) + { + case Cell_type::THREE_ONE: ++result.n3_31; break; + case Cell_type::TWO_TWO: ++result.n3_22; break; + case Cell_type::ONE_THREE: ++result.n3_13; break; + case Cell_type::ACAUSAL: + case Cell_type::UNCLASSIFIED: break; + } + } + return result; + } + + [[nodiscard]] auto canonical_state(Manifold const& manifold) + -> Canonical_state + { + auto triangulation = manifold.delaunay_snapshot(); + Canonical_state result{}; + result.counts = direct_counts(triangulation); + result.min_time = manifold.min_time(); + result.max_time = manifold.max_time(); + result.initial_radius = manifold.initial_radius(); + result.foliation_spacing = manifold.foliation_spacing(); + + for (auto vertex = triangulation.finite_vertices_begin(); + vertex != triangulation.finite_vertices_end(); ++vertex) + { + result.vertices.emplace_back(vertex_key(vertex)); + } + + for (auto edge = triangulation.finite_edges_begin(); + edge != triangulation.finite_edges_end(); ++edge) + { + result.edges.emplace_back(sorted_simplex<2>({ + edge->first->vertex(edge->second), + edge->first->vertex(edge->third), + })); + } + + std::map spacelike_facets; + for (auto facet = triangulation.finite_facets_begin(); + facet != triangulation.finite_facets_end(); ++facet) + { + std::array vertices{}; + auto output = vertices.begin(); + for (int index = 0; index < 4; ++index) + { + if (index != facet->second) { *output++ = facet->first->vertex(index); } + } + auto key = sorted_simplex<3>(vertices); + if (std::get<3>(key[0]) == std::get<3>(key[1]) && + std::get<3>(key[1]) == std::get<3>(key[2])) + { + ++spacelike_facets[std::get<3>(key[0])]; + } + result.facets.emplace_back(std::move(key)); + } + + for (auto cell = triangulation.finite_cells_begin(); + cell != triangulation.finite_cells_end(); ++cell) + { + result.cells.emplace_back( + sorted_simplex<4>({cell->vertex(0), cell->vertex(1), cell->vertex(2), + cell->vertex(3)}), + cell->info()); + } + + std::ranges::sort(result.vertices); + std::ranges::sort(result.edges); + std::ranges::sort(result.facets); + std::ranges::sort(result.cells); + result.spacelike_facets.assign(spacelike_facets.begin(), + spacelike_facets.end()); + return result; + } + + [[nodiscard]] auto canonical_triangulation(Delaunay const& triangulation) + -> Canonical_state + { + Canonical_state result{}; + result.counts = direct_counts(triangulation); + std::map spacelike_facets; + + for (auto vertex = triangulation.finite_vertices_begin(); + vertex != triangulation.finite_vertices_end(); ++vertex) + { + result.vertices.emplace_back(vertex_key(vertex)); + } + if (!result.vertices.empty()) + { + auto const times = + result.vertices | + std::views::transform([](auto const& v) { return std::get<3>(v); }); + auto const [minimum, maximum] = std::ranges::minmax(times); + result.min_time = minimum; + result.max_time = maximum; + } + + for (auto edge = triangulation.finite_edges_begin(); + edge != triangulation.finite_edges_end(); ++edge) + { + result.edges.emplace_back(sorted_simplex<2>({ + edge->first->vertex(edge->second), + edge->first->vertex(edge->third), + })); + } + + for (auto facet = triangulation.finite_facets_begin(); + facet != triangulation.finite_facets_end(); ++facet) + { + std::array vertices{}; + auto output = vertices.begin(); + for (int index = 0; index < 4; ++index) + { + if (index != facet->second) { *output++ = facet->first->vertex(index); } + } + auto key = sorted_simplex<3>(vertices); + if (std::get<3>(key[0]) == std::get<3>(key[1]) && + std::get<3>(key[1]) == std::get<3>(key[2])) + { + ++spacelike_facets[std::get<3>(key[0])]; + } + result.facets.emplace_back(std::move(key)); + } + + for (auto cell = triangulation.finite_cells_begin(); + cell != triangulation.finite_cells_end(); ++cell) + { + result.cells.emplace_back( + sorted_simplex<4>({cell->vertex(0), cell->vertex(1), cell->vertex(2), + cell->vertex(3)}), + cell->info()); + } + + std::ranges::sort(result.vertices); + std::ranges::sort(result.edges); + std::ranges::sort(result.facets); + std::ranges::sort(result.cells); + result.spacelike_facets.assign(spacelike_facets.begin(), + spacelike_facets.end()); + return result; + } + + void check_independent_invariants(Manifold const& manifold) + { + auto const triangulation = manifold.delaunay_snapshot(); + auto const counts = direct_counts(triangulation); + + CHECK(triangulation.tds().is_valid()); + CHECK(manifold.is_foliated()); + CHECK(manifold.check_vertices()); + CHECK(manifold.check_simplices()); + CHECK(manifold.is_correct()); + + CHECK_EQ(counts.n0, manifold.N0()); + CHECK_EQ(counts.n0, manifold.vertices()); + CHECK_EQ(counts.n1_sl + counts.n1_tl, manifold.N1()); + CHECK_EQ(counts.n1_sl + counts.n1_tl, manifold.edges()); + CHECK_EQ(counts.n1_sl, manifold.N1_SL()); + CHECK_EQ(counts.n1_tl, manifold.N1_TL()); + CHECK_EQ(counts.n2, manifold.N2()); + CHECK_EQ(counts.n2, manifold.faces()); + CHECK_EQ(counts.n3, manifold.N3()); + CHECK_EQ(counts.n3, manifold.simplices()); + CHECK_EQ(counts.n3_31, manifold.N3_31()); + CHECK_EQ(counts.n3_22, manifold.N3_22()); + CHECK_EQ(counts.n3_13, manifold.N3_13()); + CHECK_EQ(counts.n3_31 + counts.n3_13, manifold.N3_31_13()); + CHECK_EQ(counts.n0 - counts.n1_sl - counts.n1_tl + counts.n2 - counts.n3, + 1); + + std::map spacelike_facets; + for (auto facet = triangulation.finite_facets_begin(); + facet != triangulation.finite_facets_end(); ++facet) + { + std::optional time; + auto spacelike = true; + for (int index = 0; index < 4; ++index) + { + if (index == facet->second) { continue; } + auto const vertex_time = facet->first->vertex(index)->info(); + if (!time) { time = vertex_time; } + else if (*time != vertex_time) { spacelike = false; } + } + if (spacelike) { ++spacelike_facets[*time]; } + } + for (auto time = manifold.min_time(); time <= manifold.max_time(); ++time) + { + CHECK_EQ(manifold.spacelike_face_count(time), spacelike_facets[time]); + } + + for (auto cell = triangulation.finite_cells_begin(); + cell != triangulation.finite_cells_end(); ++cell) + { + Cell_handle const cell_handle = cell; + auto const type = independent_cell_type(cell); + CHECK((type == Cell_type::THREE_ONE || type == Cell_type::TWO_TWO || + type == Cell_type::ONE_THREE)); + CHECK_EQ(cell->info(), static_cast(type)); + for (int index = 0; index < 4; ++index) + { + auto const neighbor = cell->neighbor(index); + REQUIRE(neighbor != nullptr); + REQUIRE(neighbor->has_neighbor(cell_handle)); + CHECK(neighbor->neighbor(neighbor->index(cell_handle)) == cell_handle); + } + } + } + + void check_delta(Manifold const& before, Manifold const& after, + Count_delta const expected) + { + auto const lhs = direct_counts(before.delaunay_snapshot()); + auto const rhs = direct_counts(after.delaunay_snapshot()); + CHECK_EQ(rhs.n0 - lhs.n0, expected.n0); + CHECK_EQ(rhs.n1_sl - lhs.n1_sl, expected.n1_sl); + CHECK_EQ(rhs.n1_tl - lhs.n1_tl, expected.n1_tl); + CHECK_EQ(rhs.n2 - lhs.n2, expected.n2); + CHECK_EQ(rhs.n3_31 - lhs.n3_31, expected.n3_31); + CHECK_EQ(rhs.n3_22 - lhs.n3_22, expected.n3_22); + CHECK_EQ(rhs.n3_13 - lhs.n3_13, expected.n3_13); + CHECK_EQ(rhs.n3 - lhs.n3, expected.n3); + } + + [[nodiscard]] auto make_23_fixture() -> Manifold + { + static auto constexpr radius_2 = 2.0 * std::numbers::inv_sqrt3_v; + static auto constexpr sqrt_2 = std::numbers::sqrt2_v; + std::vector> vertices{ + { 1, 0, 0}, + { 0, 1, 0}, + { 0, 0, 1}, + {radius_2, radius_2, radius_2}, + { sqrt_2, sqrt_2, 0}, + }; + std::vector const times{1, 1, 1, 2, 2}; + return Manifold{manifolds::make_causal_vertices<3>(vertices, times)}; + } + + [[nodiscard]] auto make_26_fixture() -> Manifold + { + static auto constexpr radius_2 = 2.0 * std::numbers::inv_sqrt3_v; + std::vector> vertices{ + { 0, 0, 0}, + { 1, 0, 0}, + { 0, 1, 0}, + { 0, 0, 1}, + {radius_2, radius_2, radius_2}, + }; + std::vector const times{0, 1, 1, 1, 2}; + return Manifold{manifolds::make_causal_vertices<3>(vertices, times)}; + } + + [[nodiscard]] auto make_44_fixture() -> Manifold + { + static auto constexpr inverse_sqrt_2 = 1.0 / std::numbers::sqrt2_v; + std::vector> vertices{ + { 0, 0, 0}, + { inverse_sqrt_2, 0, inverse_sqrt_2}, + { 0, inverse_sqrt_2, inverse_sqrt_2}, + {-inverse_sqrt_2, 0, inverse_sqrt_2}, + { 0, -inverse_sqrt_2, inverse_sqrt_2}, + { 0, 0, 2}, + }; + std::vector const times{0, 1, 1, 1, 1, 2}; + return Manifold{manifolds::make_causal_vertices<3>(vertices, times), 0, 1}; + } + + [[nodiscard]] auto finite_incident_cells(Delaunay const& triangulation, + Edge_handle const edge) + -> std::optional + { + auto circulator = triangulation.incident_cells(edge, edge.first); + auto const done = circulator; + ergodic_moves::Cell_container cells; + do + { // NOLINT(cppcoreguidelines-avoid-do-while) + if (triangulation.is_infinite(circulator)) { return std::nullopt; } + cells.emplace_back(circulator); + } + while (++circulator != done); + return cells; + } + + [[nodiscard]] auto find_inverse_32_edge(Delaunay const& triangulation) + -> std::optional + { + for (auto edge = triangulation.finite_edges_begin(); + edge != triangulation.finite_edges_end(); ++edge) + { + auto const first_time = edge->first->vertex(edge->second)->info(); + auto const second_time = edge->first->vertex(edge->third)->info(); + if (first_time == second_time) { continue; } + auto const cells = finite_incident_cells(triangulation, *edge); + if (!cells || cells->size() != 3) { continue; } + auto const count = [&](Cell_type const type) { + return std::ranges::count_if(*cells, [&](auto const cell) { + return independent_cell_type(cell) == type; + }); + }; + if (count(Cell_type::TWO_TWO) == 2 && + count(Cell_type::THREE_ONE) + count(Cell_type::ONE_THREE) == 1) + { + return *edge; + } + } + return std::nullopt; + } + + [[nodiscard]] auto find_44_pivot(Delaunay const& triangulation) + -> std::optional + { + for (auto edge = triangulation.finite_edges_begin(); + edge != triangulation.finite_edges_end(); ++edge) + { + auto const first_time = edge->first->vertex(edge->second)->info(); + auto const second_time = edge->first->vertex(edge->third)->info(); + if (first_time != second_time) { continue; } + auto const cells = finite_incident_cells(triangulation, *edge); + if (!cells || cells->size() != 4) { continue; } + auto const count = [&](Cell_type const type) { + return std::ranges::count_if(*cells, [&](auto const cell) { + return independent_cell_type(cell) == type; + }); + }; + if (count(Cell_type::THREE_ONE) == 2 && count(Cell_type::ONE_THREE) == 2) + { + return *edge; + } + } + return std::nullopt; + } + + [[nodiscard]] auto find_vertex(Delaunay const& triangulation, + Point_t<3> const& point) + -> std::optional + { + for (auto vertex = triangulation.finite_vertices_begin(); + vertex != triangulation.finite_vertices_end(); ++vertex) + { + if (vertex->point() == point) { return vertex; } + } + return std::nullopt; + } + + void check_23_32_round_trip(std::uint64_t const seed) + { + auto const original = make_23_fixture(); + auto const before = canonical_state(original); + cdt::Random random{seed}; + CAPTURE(seed); + auto forward = ergodic_moves::do_23_move(original, random); + REQUIRE(forward.has_value()); + check_independent_invariants(*forward); + auto inverse = ergodic_moves::do_32_move(*forward, random); + REQUIRE(inverse.has_value()); + check_independent_invariants(*inverse); + CHECK_EQ(canonical_state(*inverse), before); + } + + void check_26_62_round_trip(std::uint64_t const seed) + { + auto const original = make_26_fixture(); + auto const before = canonical_state(original); + cdt::Random random{seed}; + CAPTURE(seed); + auto forward = ergodic_moves::do_26_move(original, random); + REQUIRE(forward.has_value()); + check_independent_invariants(*forward); + auto inverse = ergodic_moves::do_62_move(*forward, random); + REQUIRE(inverse.has_value()); + check_independent_invariants(*inverse); + CHECK_EQ(canonical_state(*inverse), before); + } + + void check_44_round_trip(std::uint64_t const seed) + { + auto const original = make_44_fixture(); + auto const before = canonical_state(original); + cdt::Random random{seed}; + CAPTURE(seed); + auto forward = ergodic_moves::do_44_move(original, random); + REQUIRE(forward.has_value()); + check_independent_invariants(*forward); + auto inverse = ergodic_moves::do_44_move(*forward, random); + REQUIRE(inverse.has_value()); + check_independent_invariants(*inverse); + CHECK_EQ(canonical_state(*inverse), before); + } +} // namespace + +SCENARIO("Every 2+1D CDT move has the literature-derived local delta" * + doctest::test_suite("ergodic-audit")) +{ + GIVEN("the minimal causal cavity for a (2,3) move") + { + auto const before = make_23_fixture(); + cdt::Random random{10623}; + WHEN("the cavity is replaced by three tetrahedra") + { + auto const after = ergodic_moves::do_23_move(before, random); + REQUIRE(after.has_value()); + THEN("the independent f-vector delta is (0,0,1,2,0,1,0,1)") + { + check_delta(before, *after, {0, 0, 1, 2, 0, 1, 0, 1}); + check_independent_invariants(*after); + } + AND_WHEN("the inverse (3,2) move is applied") + { + auto const restored = ergodic_moves::do_32_move(*after, random); + REQUIRE(restored.has_value()); + THEN("the canonical complex is restored without comparing handles") + { + check_delta(*after, *restored, {0, 0, -1, -2, 0, -1, 0, -1}); + CHECK_EQ(canonical_state(*restored), canonical_state(before)); + } + } + } + } + + GIVEN("the minimal causal cavity for a (2,6) move") + { + auto const before = make_26_fixture(); + cdt::Random random{10626}; + WHEN("a vertex is inserted into the shared spacelike triangle") + { + auto const after = ergodic_moves::do_26_move(before, random); + REQUIRE(after.has_value()); + THEN("the independent f-vector delta is (1,3,2,8,2,0,2,4)") + { + check_delta(before, *after, {1, 3, 2, 8, 2, 0, 2, 4}); + check_independent_invariants(*after); + } + AND_WHEN("the inverse (6,2) move removes the inserted vertex") + { + auto const restored = ergodic_moves::do_62_move(*after, random); + REQUIRE(restored.has_value()); + THEN("the canonical complex is restored without comparing handles") + { + check_delta(*after, *restored, {-1, -3, -2, -8, -2, 0, -2, -4}); + CHECK_EQ(canonical_state(*restored), canonical_state(before)); + } + } + } + } + + GIVEN("the minimal four-tetrahedron diamond for a (4,4) move") + { + auto const before = make_44_fixture(); + cdt::Random random{10644}; + WHEN("the spacelike diagonal is exchanged") + { + auto const after = ergodic_moves::do_44_move(before, random); + REQUIRE(after.has_value()); + THEN("every tracked count is unchanged") + { + check_delta(before, *after, {}); + check_independent_invariants(*after); + CHECK_NE(canonical_state(*after), canonical_state(before)); + } + AND_WHEN("the inverse diagonal exchange is applied") + { + auto const restored = ergodic_moves::do_44_move(*after, random); + REQUIRE(restored.has_value()); + THEN("the original canonical diamond is restored") + { CHECK_EQ(canonical_state(*restored), canonical_state(before)); } + } + } + } +} + +SCENARIO("CDT move rejection is causal and failure-atomic" * + doctest::test_suite("ergodic-audit")) +{ + GIVEN("an empty source manifold") + { + Manifold const source; + auto const before = canonical_state(source); + cdt::Random random{106}; + WHEN("each unsupported move is attempted") + { + CHECK_FALSE(ergodic_moves::do_23_move(source, random)); + CHECK_FALSE(ergodic_moves::do_32_move(source, random)); + CHECK_FALSE(ergodic_moves::do_26_move(source, random)); + CHECK_FALSE(ergodic_moves::do_62_move(source, random)); + CHECK_FALSE(ergodic_moves::do_44_move(source, random)); + THEN("the complete canonical state and caches are unchanged") + { CHECK_EQ(canonical_state(source), before); } + } + } + + GIVEN("a geometrically flippable timelike edge with the wrong CDT cavity") + { + auto triangulation = make_23_fixture().delaunay_snapshot(); + auto two_two = foliated_triangulations::filter_cells<3>( + foliated_triangulations::collect_cells<3>(triangulation), + Cell_type::TWO_TWO); + REQUIRE_EQ(two_two.size(), 1); + REQUIRE(ergodic_moves::try_23_move(triangulation, two_two.front())); + assign_independent_cell_metadata(triangulation); + auto const inverse_edge = find_inverse_32_edge(triangulation); + REQUIRE(inverse_edge.has_value()); + + auto const first = inverse_edge->first->vertex(inverse_edge->second); + auto const second = inverse_edge->first->vertex(inverse_edge->third); + auto const lower_time = std::min(first->info(), second->info()); + auto const incident = finite_incident_cells(triangulation, *inverse_edge); + REQUIRE(incident.has_value()); + for (auto const cell : *incident) + { + for (int index = 0; index < 4; ++index) + { + auto const vertex = cell->vertex(index); + if (vertex != first && vertex != second) + { + vertex->info() = lower_time; + } + } + } + assign_independent_cell_metadata(triangulation); + auto const before = canonical_triangulation(triangulation); + + auto control = triangulation; + auto const control_first = find_vertex(control, first->point()); + auto const control_second = find_vertex(control, second->point()); + REQUIRE(control_first.has_value()); + REQUIRE(control_second.has_value()); + Cell_handle control_cell = nullptr; + int control_i{}; + int control_j{}; + REQUIRE(control.is_edge(*control_first, *control_second, control_cell, + control_i, control_j)); + REQUIRE(control.flip(control_cell, control_i, control_j)); + + WHEN("the CDT (3,2) helper checks the site") + { + CHECK_FALSE(ergodic_moves::try_32_move(triangulation, *inverse_edge)); + THEN("it rejects before CGAL can create an acausal tetrahedron") + { CHECK_EQ(canonical_triangulation(triangulation), before); } + } + } + + GIVEN("valid topology with stale or incompatible causal metadata") + { + auto const corrupt_metadata = [](Cell_handle const& cell) { + auto const expected = + foliated_triangulations::expected_cell_type<3>(cell); + cell->info() = static_cast(expected == Cell_type::TWO_TWO + ? Cell_type::THREE_ONE + : Cell_type::TWO_TWO); + }; + + auto move_23 = make_23_fixture().delaunay_snapshot(); + auto two_two = foliated_triangulations::filter_cells<3>( + foliated_triangulations::collect_cells<3>(move_23), Cell_type::TWO_TWO); + REQUIRE_EQ(two_two.size(), 1); + corrupt_metadata(two_two.front()); + auto const before_23 = canonical_triangulation(move_23); + + auto move_32 = make_23_fixture().delaunay_snapshot(); + auto move_32_cells = foliated_triangulations::filter_cells<3>( + foliated_triangulations::collect_cells<3>(move_32), Cell_type::TWO_TWO); + REQUIRE_EQ(move_32_cells.size(), 1); + REQUIRE(ergodic_moves::try_23_move(move_32, move_32_cells.front())); + assign_independent_cell_metadata(move_32); + auto const inverse_edge = find_inverse_32_edge(move_32); + REQUIRE(inverse_edge.has_value()); + auto const inverse_cells = finite_incident_cells(move_32, *inverse_edge); + REQUIRE(inverse_cells.has_value()); + corrupt_metadata(inverse_cells->front()); + auto const before_32 = canonical_triangulation(move_32); + + auto move_26 = make_26_fixture().delaunay_snapshot(); + auto one_three = foliated_triangulations::filter_cells<3>( + foliated_triangulations::collect_cells<3>(move_26), + Cell_type::ONE_THREE); + REQUIRE_EQ(one_three.size(), 1); + corrupt_metadata(one_three.front()); + auto const before_26 = canonical_triangulation(move_26); + + auto const source_62 = make_26_fixture(); + cdt::Random random_62{10662}; + auto const expanded_62 = ergodic_moves::do_26_move(source_62, random_62); + REQUIRE(expanded_62.has_value()); + auto move_62 = expanded_62->delaunay_snapshot(); + auto vertices_62 = foliated_triangulations::collect_vertices<3>(move_62); + auto const candidate_62 = + std::ranges::find_if(vertices_62, [&](auto const& vertex) { + return ergodic_moves::is_62_movable(move_62, vertex); + }); + REQUIRE(candidate_62 != vertices_62.end()); + ergodic_moves::Cell_container incident_62; + move_62.tds().incident_cells(*candidate_62, + std::back_inserter(incident_62)); + REQUIRE_EQ(incident_62.size(), 6); + corrupt_metadata(incident_62.front()); + auto const before_62 = canonical_triangulation(move_62); + + auto move_44 = make_44_fixture().delaunay_snapshot(); + auto const pivot = find_44_pivot(move_44); + REQUIRE(pivot.has_value()); + auto const pivot_cells = finite_incident_cells(move_44, *pivot); + REQUIRE(pivot_cells.has_value()); + corrupt_metadata(pivot_cells->front()); + auto const before_44 = canonical_triangulation(move_44); + + WHEN("the corresponding cavity selectors inspect the sites") + { + CHECK_FALSE(ergodic_moves::try_23_move(move_23, two_two.front())); + CHECK_FALSE(ergodic_moves::try_32_move(move_32, *inverse_edge)); + CHECK_FALSE(ergodic_moves::find_adjacent_31_cell(one_three.front())); + CHECK_FALSE(ergodic_moves::is_62_movable(move_62, *candidate_62)); + CHECK_FALSE(ergodic_moves::find_bistellar_flip_location(move_44, *pivot)); + THEN("no topology is mutated") + { + CHECK_EQ(canonical_triangulation(move_23), before_23); + CHECK_EQ(canonical_triangulation(move_32), before_32); + CHECK_EQ(canonical_triangulation(move_26), before_26); + CHECK_EQ(canonical_triangulation(move_62), before_62); + CHECK_EQ(canonical_triangulation(move_44), before_44); + } + } + } + + GIVEN("a topological cavity with an invalid time assignment") + { + auto triangulation = make_23_fixture().delaunay_snapshot(); + auto two_two = foliated_triangulations::filter_cells<3>( + foliated_triangulations::collect_cells<3>(triangulation), + Cell_type::TWO_TWO); + REQUIRE_EQ(two_two.size(), 1); + two_two.front()->vertex(0)->info() += 3; + assign_independent_cell_metadata(triangulation); + auto const before = canonical_triangulation(triangulation); + + WHEN("the (2,3) selector checks the now-acausal tetrahedron") + { + CHECK_FALSE(ergodic_moves::try_23_move(triangulation, two_two.front())); + THEN("the invalid site is rejected without mutation") + { CHECK_EQ(canonical_triangulation(triangulation), before); } + } + } + + GIVEN("causal manifolds without the required move incidences") + { + auto const source_23 = make_23_fixture(); + auto const source_26 = make_26_fixture(); + auto const before_23 = canonical_state(source_23); + auto const before_26 = canonical_state(source_26); + cdt::Random random{1060}; + + WHEN("degree, composition, or flippability preconditions are absent") + { + CHECK_FALSE(ergodic_moves::do_32_move(source_23, random)); + CHECK_FALSE(ergodic_moves::do_26_move(source_23, random)); + CHECK_FALSE(ergodic_moves::do_44_move(source_23, random)); + CHECK_FALSE(ergodic_moves::do_62_move(source_26, random)); + THEN("each source remains canonically unchanged") + { + CHECK_EQ(canonical_state(source_23), before_23); + CHECK_EQ(canonical_state(source_26), before_26); + } + } + } + + GIVEN("malformed and nonflippable edge representations") + { + auto triangulation = make_44_fixture().delaunay_snapshot(); + Edge_handle const malformed{nullptr, 0, 1}; + auto const before = canonical_triangulation(triangulation); + WHEN("the checked edge helpers receive them") + { + CHECK_FALSE(ergodic_moves::get_incident_cells(triangulation, malformed)); + CHECK_FALSE(ergodic_moves::try_32_move(triangulation, malformed)); + THEN("the triangulation remains unchanged") + { CHECK_EQ(canonical_triangulation(triangulation), before); } + } + } +} + +TEST_CASE("Seeded CDT move replay preserves invariants and inverse structure" * + doctest::test_suite("ergodic-audit")) +{ + for (auto const seed : + std::array{0, 1, 2, 92, 106, 20'260'721}) + { + check_23_32_round_trip(seed); + check_26_62_round_trip(seed); + check_44_round_trip(seed); + } +} From f4c7df9637592cafaf41014d20fba89288000249 Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Tue, 21 Jul 2026 08:35:43 -0700 Subject: [PATCH 2/3] fix(ci): restore sanitizer and coverage reliability - compare preserved move configuration values by representation to satisfy strict floating-point warnings without weakening drift detection - allow coverage-instrumented unit tests enough time to complete --- .github/workflows/codecov-upload.yml | 4 +++- include/Ergodic_moves_3.hpp | 20 ++++++++++++++-- tests/Ergodic_moves_3_audit_test.cpp | 35 ++++++++++++++++++++++------ 3 files changed, 49 insertions(+), 10 deletions(-) diff --git a/.github/workflows/codecov-upload.yml b/.github/workflows/codecov-upload.yml index 45d64d811..91cd69cbb 100644 --- a/.github/workflows/codecov-upload.yml +++ b/.github/workflows/codecov-upload.yml @@ -56,7 +56,9 @@ jobs: - name: Test working-directory: build - run: ctest --no-tests=error --output-on-failure -j 2 + # Coverage instrumentation makes the monolithic unit-test binary + # substantially slower than it is in normal CI. + run: ctest --no-tests=error --output-on-failure --timeout 600 -j 2 - name: Submit to codecov.io uses: codecov/codecov-action@fb8b3582c8e4def4969c97caa2f19720cb33a72f # v7.0.0 diff --git a/include/Ergodic_moves_3.hpp b/include/Ergodic_moves_3.hpp index f28c49bfe..9a9b94095 100644 --- a/include/Ergodic_moves_3.hpp +++ b/include/Ergodic_moves_3.hpp @@ -19,6 +19,9 @@ #ifndef CDT_PLUSPLUS_ERGODIC_MOVES_3_HPP #define CDT_PLUSPLUS_ERGODIC_MOVES_3_HPP +#include +#include +#include #include #include @@ -39,6 +42,17 @@ namespace ergodic_moves namespace detail { + /// @brief Compare preserved floating-point configuration state exactly. + /// @details Move construction copies these values; arithmetic tolerance is + /// inappropriate because any representation change indicates state drift. + [[nodiscard]] inline auto same_configuration_value( + double const first, double const second) noexcept -> bool + { + using Representation = std::array; + return std::bit_cast(first) == + std::bit_cast(second); + } + /// @brief Rebuild all derived topology and geometry state around a value. [[nodiscard]] inline auto make_manifold(Delaunay triangulation, Manifold const& source) -> Manifold @@ -1088,8 +1102,10 @@ namespace ergodic_moves -> bool { if (!t_after.is_correct() || - t_after.initial_radius() != t_before.initial_radius() || - t_after.foliation_spacing() != t_before.foliation_spacing()) + !detail::same_configuration_value(t_after.initial_radius(), + t_before.initial_radius()) || + !detail::same_configuration_value(t_after.foliation_spacing(), + t_before.foliation_spacing())) { return false; } diff --git a/tests/Ergodic_moves_3_audit_test.cpp b/tests/Ergodic_moves_3_audit_test.cpp index 6cbbc9652..17ec2c797 100644 --- a/tests/Ergodic_moves_3_audit_test.cpp +++ b/tests/Ergodic_moves_3_audit_test.cpp @@ -11,6 +11,8 @@ #include #include +#include +#include #include #include #include @@ -32,6 +34,11 @@ namespace using Vertex_handle = ergodic_moves::Vertex_handle; using Vertex_key = std::tuple; + using Double_representation = std::array; + + [[nodiscard]] auto double_representation(double const value) noexcept + -> Double_representation + { return std::bit_cast(value); } template using Simplex_key = std::array; @@ -74,8 +81,8 @@ namespace std::vector> spacelike_facets; Int_precision min_time{}; Int_precision max_time{}; - double initial_radius{}; - double foliation_spacing{}; + Double_representation initial_radius{}; + Double_representation foliation_spacing{}; auto operator==(Canonical_state const&) const -> bool = default; }; @@ -165,11 +172,12 @@ namespace { auto triangulation = manifold.delaunay_snapshot(); Canonical_state result{}; - result.counts = direct_counts(triangulation); - result.min_time = manifold.min_time(); - result.max_time = manifold.max_time(); - result.initial_radius = manifold.initial_radius(); - result.foliation_spacing = manifold.foliation_spacing(); + result.counts = direct_counts(triangulation); + result.min_time = manifold.min_time(); + result.max_time = manifold.max_time(); + result.initial_radius = double_representation(manifold.initial_radius()); + result.foliation_spacing = + double_representation(manifold.foliation_spacing()); for (auto vertex = triangulation.finite_vertices_begin(); vertex != triangulation.finite_vertices_end(); ++vertex) @@ -838,3 +846,16 @@ TEST_CASE("Seeded CDT move replay preserves invariants and inverse structure" * check_44_round_trip(seed); } } + +TEST_CASE("Move validation rejects exact configuration-state drift" * + doctest::test_suite("ergodic-audit")) +{ + auto const before = make_44_fixture(); + auto const after = Manifold{ + foliated_triangulations::FoliatedTriangulation_3{ + before.delaunay_snapshot(), -0.0, before.foliation_spacing()} + }; + REQUIRE(after.is_correct()); + CHECK_FALSE(ergodic_moves::check_move(before, after, + move_tracker::move_type::FOUR_FOUR)); +} From cec15b57fd6b49338fe1bdfc5ffdf1f227a1aa4e Mon Sep 17 00:00:00 2001 From: Adam Getchell Date: Tue, 21 Jul 2026 10:02:37 -0700 Subject: [PATCH 3/3] fix(moves): harden ergodic move rejection Reject boundary and foreign edges before accessing triangulation metadata, and treat any infinite incident cell as an invalid move site. Preserve source manifolds when validation rejects post-mutation states for (2,6), (6,2), and (4,4) transitions. --- include/Ergodic_moves_3.hpp | 244 ++++++++++++++++++--------- tests/Ergodic_moves_3_audit_test.cpp | 182 +++++++++++++++++++- 2 files changed, 340 insertions(+), 86 deletions(-) diff --git a/include/Ergodic_moves_3.hpp b/include/Ergodic_moves_3.hpp index 9a9b94095..94cadb469 100644 --- a/include/Ergodic_moves_3.hpp +++ b/include/Ergodic_moves_3.hpp @@ -21,6 +21,7 @@ #include #include +#include #include #include #include @@ -53,6 +54,11 @@ namespace ergodic_moves std::bit_cast(second); } + /// @brief Default validator for internal post-mutation test seams. + [[nodiscard]] inline auto accept_post_mutation( + [[maybe_unused]] Delaunay const& triangulation) noexcept -> bool + { return true; } + /// @brief Rebuild all derived topology and geometry state around a value. [[nodiscard]] inline auto make_manifold(Delaunay triangulation, Manifold const& source) -> Manifold @@ -95,10 +101,8 @@ namespace ergodic_moves Cell_container incident_cells; do { // NOLINT(cppcoreguidelines-avoid-do-while) - if (!triangulation.is_infinite(circulator)) - { - incident_cells.emplace_back(circulator); - } + if (triangulation.is_infinite(circulator)) { return std::nullopt; } + incident_cells.emplace_back(circulator); } while (++circulator != edge.first); return incident_cells; @@ -375,31 +379,25 @@ namespace ergodic_moves return std::nullopt; } // find_26_move() - /// @brief Perform a (2,6) move - /// @details A (2,6) move inserts a vertex into the spacelike face between a - /// (1,3) simplex on the bottom connected to a (3,1) simplex on top. - /// This adds 2 (1,3) simplices and 2 (3,1) simplices. - /// It adds 2 spacelike faces and 6 timelike faces. - /// It also adds 2 timelike edges and 3 spacelike edges, as well as the - /// vertex. - /// This function calls find_adjacent_31_cell on (1,3) simplices drawn from a - /// randomly shuffled container until it succeeds or runs out of simplices. - /// The move guarantees a valid causal combinatorial triangulation. It does - /// not preserve or require the empty-sphere property of the coordinates. - /// @image html 26.png - /// @image latex 26.eps width=7cm - /// @tparam Generator A uniform random bit generator type - /// @param t_manifold The simplicial manifold - /// @param generator Caller-owned generator whose state advances during the - /// move - /// @param only_first_site Whether to examine only one uniformly selected site - /// @returns The Expected (2,6) moved manifold or an Unexpected - /// @note The source manifold is unchanged on success and failure. - template - [[nodiscard]] inline auto do_26_move_impl(Manifold const& t_manifold, - Generator& generator, - bool const only_first_site) - -> Expected + namespace detail + { + template + requires std::predicate + [[nodiscard]] inline auto do_26_move_impl( + Manifold const& t_manifold, Generator& generator, + bool const only_first_site, + Post_mutation_validator post_mutation_validator) -> Expected; + } // namespace detail + + // Internal validation seam used to test rejection after mutation. + template + requires std::predicate + [[nodiscard]] inline auto detail::do_26_move_impl( + Manifold const& t_manifold, Generator& generator, + bool const only_first_site, + Post_mutation_validator post_mutation_validator) -> Expected { static auto constexpr INCIDENT_CELLS_FOR_6_2_MOVE = 6; Delaunay triangulation{t_manifold.delaunay_snapshot()}; @@ -417,7 +415,8 @@ namespace ergodic_moves } else { - // Shuffle the container to pick a random sequence of (1,3) cells to try. + // Shuffle the container to pick a random sequence of (1,3) cells to + // try. std::ranges::shuffle(one_three, generator); } for (auto const& bottom : one_three) @@ -473,7 +472,8 @@ namespace ergodic_moves return std::unexpected(msg); } - // Each incident cell should be combinatorially and geometrically valid + // Each incident cell should be combinatorially and geometrically + // valid if (auto check_cells = std::ranges::all_of(incident_cells, [&triangulation](auto const& cell) { @@ -496,19 +496,52 @@ namespace ergodic_moves // Final checks // is_valid() checks for combinatorial and geometric validity - if (!triangulation.tds().is_valid(v_center, true, 1)) + if (!post_mutation_validator( + static_cast(triangulation)) || + !triangulation.tds().is_valid(v_center, true, 1)) { std::string const msg = "v_center is invalid.\n"; return std::unexpected(msg); } - return detail::make_manifold(std::move(triangulation), t_manifold); + return make_manifold(std::move(triangulation), t_manifold); } } // We've run out of (1,3) simplices to try std::string const msg = "No (2,6) move possible.\n"; return std::unexpected(msg); - } // do_26_move() + } // do_26_move_impl() + + /// @brief Perform a (2,6) move + /// @details A (2,6) move inserts a vertex into the spacelike face between a + /// (1,3) simplex on the bottom connected to a (3,1) simplex on top. + /// This adds 2 (1,3) simplices and 2 (3,1) simplices. + /// It adds 2 spacelike faces and 6 timelike faces. + /// It also adds 2 timelike edges and 3 spacelike edges, as well as the + /// vertex. + /// This function calls find_adjacent_31_cell on (1,3) simplices drawn from a + /// randomly shuffled container until it succeeds or runs out of simplices. + /// The move guarantees a valid causal combinatorial triangulation. It does + /// not preserve or require the empty-sphere property of the coordinates. + /// @image html 26.png + /// @image latex 26.eps width=7cm + /// @tparam Generator A uniform random bit generator type + /// @param t_manifold The simplicial manifold + /// @param generator Caller-owned generator whose state advances during the + /// move + /// @param only_first_site Whether to examine only one uniformly selected site + /// @returns The Expected (2,6) moved manifold or an Unexpected + /// @note The source manifold is unchanged on success and failure. + + template + [[nodiscard]] inline auto do_26_move_impl(Manifold const& t_manifold, + Generator& generator, + bool const only_first_site) + -> Expected + { + return detail::do_26_move_impl(t_manifold, generator, only_first_site, + detail::accept_post_mutation); + } /// @brief Perform a (2,6) move using a caller-owned random stream. template @@ -587,20 +620,26 @@ namespace ergodic_moves } // find_62_moves() - /// @brief Apply the combinatorial (6,2) retriangulation on a private copy. - /// @details A (3,2) flip of either timelike edge incident to the removable - /// vertex leaves that vertex with degree four. Removing it from its maximal - /// simplex then replaces the original six cells with the required two cells. - /// Keeping both operations in a private copy makes rejection failure-atomic. - /// @param source_triangulation The triangulation containing the candidate - /// @param source_candidate The degree-five vertex to remove - /// @param generator The caller-owned random engine used to order flip paths - /// @returns The moved triangulation, or nullopt when the topology is not - /// flippable or the result violates a triangulation or causal-cell invariant - template - [[nodiscard]] inline auto try_62_move(Delaunay const& source_triangulation, - Vertex_handle const source_candidate, - Generator& generator) + namespace detail + { + template + requires std::predicate + [[nodiscard]] inline auto try_62_move_impl( + Delaunay const& source_triangulation, + Vertex_handle const source_candidate, Generator& generator, + Post_mutation_validator post_mutation_validator) + -> std::optional; + } // namespace detail + + // Internal validation seam used to test rejection after mutation. + template + requires std::predicate + [[nodiscard]] inline auto detail::try_62_move_impl( + Delaunay const& source_triangulation, + Vertex_handle const source_candidate, Generator& generator, + Post_mutation_validator post_mutation_validator) -> std::optional { if (!source_triangulation.tds().is_vertex(source_candidate) || @@ -639,7 +678,8 @@ namespace ergodic_moves } tds.remove_from_maximal_dimension_simplex(candidate); - if (!tds.is_valid() || + if (!post_mutation_validator(static_cast(triangulation)) || + !tds.is_valid() || triangulation.number_of_finite_cells() + 4 != old_cells || triangulation.number_of_vertices() + 1 != old_vertices) { @@ -657,6 +697,26 @@ namespace ergodic_moves } return triangulation; + } // try_62_move_impl() + + /// @brief Apply the combinatorial (6,2) retriangulation on a private copy. + /// @details A (3,2) flip of either timelike edge incident to the removable + /// vertex leaves that vertex with degree four. Removing it from its maximal + /// simplex then replaces the original six cells with the required two cells. + /// Keeping both operations in a private copy makes rejection failure-atomic. + /// @param source_triangulation The triangulation containing the candidate + /// @param source_candidate The degree-five vertex to remove + /// @param generator The caller-owned random engine used to order flip paths + /// @returns The moved triangulation, or nullopt when the topology is not + /// flippable or the result violates a triangulation or causal-cell invariant + template + [[nodiscard]] inline auto try_62_move(Delaunay const& source_triangulation, + Vertex_handle const source_candidate, + Generator& generator) + -> std::optional + { + return detail::try_62_move_impl(source_triangulation, source_candidate, + generator, detail::accept_post_mutation); } // try_62_move() /// @brief Perform a (6,2) move @@ -747,27 +807,26 @@ namespace ergodic_moves -> std::optional { if (!detail::is_well_formed_edge(t_edge_candidate)) { return std::nullopt; } + auto incident_cells = + incident_cells_from_edge(triangulation, t_edge_candidate); + if (!incident_cells || incident_cells->size() != 4) { return std::nullopt; } + auto const first_time = t_edge_candidate.first->vertex(t_edge_candidate.second)->info(); auto const second_time = t_edge_candidate.first->vertex(t_edge_candidate.third)->info(); if (first_time != second_time) { return std::nullopt; } - if (auto incident_cells = - incident_cells_from_edge(triangulation, t_edge_candidate); - incident_cells.has_value() && incident_cells->size() == 4) + auto const cell_type_count = [&](Cell_type const type) { + return std::ranges::count_if(*incident_cells, [&](auto const cell) { + return foliated_triangulations::is_cell_type_correct<3>(cell) && + foliated_triangulations::expected_cell_type<3>(cell) == type; + }); + }; + if (cell_type_count(Cell_type::THREE_ONE) == 2 && + cell_type_count(Cell_type::ONE_THREE) == 2) { - auto const cell_type_count = [&](Cell_type const type) { - return std::ranges::count_if(*incident_cells, [&](auto const cell) { - return foliated_triangulations::is_cell_type_correct<3>(cell) && - foliated_triangulations::expected_cell_type<3>(cell) == type; - }); - }; - if (cell_type_count(Cell_type::THREE_ONE) == 2 && - cell_type_count(Cell_type::ONE_THREE) == 2) - { - return incident_cells.value(); - } + return incident_cells; } return std::nullopt; } // find_bistellar_flip_location() @@ -783,22 +842,27 @@ namespace ergodic_moves return incident_cells_from_edge(triangulation, edge); } // get_incident_cells() - /// @brief Perform a bistellar flip on triangulation via the given edge - /// @details Pass by value to avoid modifying the original triangulation - /// in the event that the flip is unsuccessful. - /// @param source_triangulation The triangulation to flip - /// @param source_edge The edge to pivot on - /// @param source_top Top vertex of the cells being flipped - /// @param source_bottom Bottom vertex of the cells being flipped - /// @returns A flipped triangulation or nullopt - /// @see [Pachner moves](../REFERENCES.md#pachner-moves) - [[nodiscard]] inline auto bistellar_flip(Delaunay const& source_triangulation, - Edge_handle const source_edge, - Vertex_handle const source_top, - Vertex_handle const source_bottom) + namespace detail + { + template + requires std::predicate + [[nodiscard]] inline auto bistellar_flip_impl( + Delaunay const& source_triangulation, Edge_handle const source_edge, + Vertex_handle const source_top, Vertex_handle const source_bottom, + Post_mutation_validator post_mutation_validator) + -> std::optional; + } // namespace detail + + // Internal validation seam used to test rejection after mutation. + template + requires std::predicate + [[nodiscard]] inline auto detail::bistellar_flip_impl( + Delaunay const& source_triangulation, Edge_handle const source_edge, + Vertex_handle const source_top, Vertex_handle const source_bottom, + Post_mutation_validator post_mutation_validator) -> std::optional { - if (!detail::is_well_formed_edge(source_edge) || source_top == nullptr || + if (!is_well_formed_edge(source_edge) || source_top == nullptr || source_bottom == nullptr || !source_triangulation.tds().is_edge( source_edge.first, source_edge.second, source_edge.third) || @@ -911,7 +975,11 @@ namespace ergodic_moves return std::nullopt; } - if (!triangulation.tds().is_valid()) { return std::nullopt; } + if (!post_mutation_validator(static_cast(triangulation)) || + !triangulation.tds().is_valid()) + { + return std::nullopt; + } for (auto const cell : triangulation.finite_cell_handles()) { @@ -920,6 +988,26 @@ namespace ergodic_moves } return triangulation; + } // bistellar_flip_impl() + + /// @brief Perform a bistellar flip on triangulation via the given edge + /// @details Pass by value to avoid modifying the original triangulation + /// in the event that the flip is unsuccessful. + /// @param source_triangulation The triangulation to flip + /// @param source_edge The edge to pivot on + /// @param source_top Top vertex of the cells being flipped + /// @param source_bottom Bottom vertex of the cells being flipped + /// @returns A flipped triangulation or nullopt + /// @see [Pachner moves](../REFERENCES.md#pachner-moves) + [[nodiscard]] inline auto bistellar_flip(Delaunay const& source_triangulation, + Edge_handle const source_edge, + Vertex_handle const source_top, + Vertex_handle const source_bottom) + -> std::optional + { + return detail::bistellar_flip_impl(source_triangulation, source_edge, + source_top, source_bottom, + detail::accept_post_mutation); } // bistellar_flip() /// @return The center edge of a 4-cell complex diff --git a/tests/Ergodic_moves_3_audit_test.cpp b/tests/Ergodic_moves_3_audit_test.cpp index 17ec2c797..26e7a3d81 100644 --- a/tests/Ergodic_moves_3_audit_test.cpp +++ b/tests/Ergodic_moves_3_audit_test.cpp @@ -300,14 +300,39 @@ namespace void check_independent_invariants(Manifold const& manifold) { - auto const triangulation = manifold.delaunay_snapshot(); - auto const counts = direct_counts(triangulation); - + auto const triangulation = manifold.delaunay_snapshot(); + auto const counts = direct_counts(triangulation); + + auto const production_is_foliated = manifold.is_foliated(); + auto const production_checks_vertices = manifold.check_vertices(); + auto const production_checks_simplices = manifold.check_simplices(); + auto const production_is_correct = manifold.is_correct(); + CAPTURE(production_is_foliated); + CAPTURE(production_checks_vertices); + CAPTURE(production_checks_simplices); + CAPTURE(production_is_correct); + + CHECK_EQ(triangulation.dimension(), 3); CHECK(triangulation.tds().is_valid()); - CHECK(manifold.is_foliated()); - CHECK(manifold.check_vertices()); - CHECK(manifold.check_simplices()); - CHECK(manifold.is_correct()); + + std::optional raw_min_time; + std::optional raw_max_time; + std::size_t raw_vertex_count{}; + for (auto vertex = triangulation.finite_vertices_begin(); + vertex != triangulation.finite_vertices_end(); ++vertex) + { + CHECK(triangulation.tds().is_vertex(vertex)); + raw_min_time = raw_min_time ? std::min(*raw_min_time, vertex->info()) + : vertex->info(); + raw_max_time = raw_max_time ? std::max(*raw_max_time, vertex->info()) + : vertex->info(); + ++raw_vertex_count; + } + REQUIRE(raw_min_time.has_value()); + REQUIRE(raw_max_time.has_value()); + CHECK_EQ(raw_vertex_count, static_cast(counts.n0)); + CHECK_EQ(*raw_min_time, manifold.min_time()); + CHECK_EQ(*raw_max_time, manifold.max_time()); CHECK_EQ(counts.n0, manifold.N0()); CHECK_EQ(counts.n0, manifold.vertices()); @@ -346,22 +371,27 @@ namespace CHECK_EQ(manifold.spacelike_face_count(time), spacelike_facets[time]); } + std::size_t classified_cell_count{}; for (auto cell = triangulation.finite_cells_begin(); cell != triangulation.finite_cells_end(); ++cell) { Cell_handle const cell_handle = cell; auto const type = independent_cell_type(cell); + CHECK(triangulation.tds().is_cell(cell_handle)); CHECK((type == Cell_type::THREE_ONE || type == Cell_type::TWO_TWO || type == Cell_type::ONE_THREE)); CHECK_EQ(cell->info(), static_cast(type)); + ++classified_cell_count; for (int index = 0; index < 4; ++index) { auto const neighbor = cell->neighbor(index); REQUIRE(neighbor != nullptr); + CHECK(triangulation.tds().is_cell(neighbor)); REQUIRE(neighbor->has_neighbor(cell_handle)); CHECK(neighbor->neighbor(neighbor->index(cell_handle)) == cell_handle); } } + CHECK_EQ(classified_cell_count, static_cast(counts.n3)); } void check_delta(Manifold const& before, Manifold const& after, @@ -798,6 +828,122 @@ SCENARIO("CDT move rejection is causal and failure-atomic" * } } + GIVEN("a (2,6) move rejected by validation after vertex insertion") + { + auto const source = make_26_fixture(); + auto const before = canonical_state(source); + auto const source_counts = direct_counts(source.delaunay_snapshot()); + cdt::Random random{106'260}; + bool observed_mutation = false; + + WHEN("the injected count validator rejects the private copy") + { + auto const rejected = ergodic_moves::detail::do_26_move_impl( + source, random, false, [&](Delaunay const& mutated) { + auto const mutated_counts = direct_counts(mutated); + observed_mutation = true; + CHECK_EQ(mutated_counts.n0, source_counts.n0 + 1); + CHECK_EQ(mutated_counts.n3, source_counts.n3 + 4); + return mutated_counts.n0 == source_counts.n0; + }); + + THEN("the post-mutation rejection leaves the source unchanged") + { + REQUIRE(observed_mutation); + CHECK_FALSE(rejected.has_value()); + CHECK_EQ(canonical_state(source), before); + } + } + } + + GIVEN("a (6,2) move rejected by validation after vertex removal") + { + auto const seed = make_26_fixture(); + cdt::Random setup_random{106'620}; + auto const expanded = ergodic_moves::do_26_move(seed, setup_random); + REQUIRE(expanded.has_value()); + + auto source = expanded->delaunay_snapshot(); + auto vertices = foliated_triangulations::collect_vertices<3>(source); + auto const candidate = + std::ranges::find_if(vertices, [&](auto const& vertex) { + return ergodic_moves::is_62_movable(source, vertex); + }); + REQUIRE(candidate != vertices.end()); + auto const source_counts = direct_counts(source); + auto const before = canonical_triangulation(source); + cdt::Random random{106'621}; + bool observed_mutation = false; + + WHEN("the injected count validator rejects the private copy") + { + auto const rejected = ergodic_moves::detail::try_62_move_impl( + source, *candidate, random, [&](Delaunay const& mutated) { + auto const mutated_counts = direct_counts(mutated); + observed_mutation = true; + CHECK_EQ(mutated_counts.n0, source_counts.n0 - 1); + CHECK_EQ(mutated_counts.n3, source_counts.n3 - 4); + return mutated_counts.n0 == source_counts.n0; + }); + + THEN("the post-mutation rejection leaves the source unchanged") + { + REQUIRE(observed_mutation); + CHECK_FALSE(rejected.has_value()); + CHECK_EQ(canonical_triangulation(source), before); + } + } + } + + GIVEN("a bistellar flip rejected by validation after both TDS flips") + { + auto source = make_44_fixture().delaunay_snapshot(); + auto const pivot = find_44_pivot(source); + REQUIRE(pivot.has_value()); + auto const incident = finite_incident_cells(source, *pivot); + REQUIRE(incident.has_value()); + + auto const pivot_first = pivot->first->vertex(pivot->second); + auto const pivot_second = pivot->first->vertex(pivot->third); + Vertex_handle top = nullptr; + Vertex_handle bottom = nullptr; + for (auto const cell : *incident) + { + for (int index = 0; index < 4; ++index) + { + auto const vertex = cell->vertex(index); + if (vertex == pivot_first || vertex == pivot_second) { continue; } + if (top == nullptr || vertex->info() > top->info()) { top = vertex; } + if (bottom == nullptr || vertex->info() < bottom->info()) + { + bottom = vertex; + } + } + } + REQUIRE(top != nullptr); + REQUIRE(bottom != nullptr); + auto const before = canonical_triangulation(source); + bool observed_mutation = false; + + WHEN("the injected TDS validator rejects the private copy") + { + auto const rejected = ergodic_moves::detail::bistellar_flip_impl( + source, *pivot, top, bottom, [&](Delaunay const& mutated) { + observed_mutation = true; + CHECK(mutated.tds().is_valid()); + CHECK_NE(canonical_triangulation(mutated), before); + return !mutated.tds().is_valid(); + }); + + THEN("the post-mutation rejection leaves the source unchanged") + { + REQUIRE(observed_mutation); + CHECK_FALSE(rejected.has_value()); + CHECK_EQ(canonical_triangulation(source), before); + } + } + } + GIVEN("causal manifolds without the required move incidences") { auto const source_23 = make_23_fixture(); @@ -822,12 +968,32 @@ SCENARIO("CDT move rejection is causal and failure-atomic" * GIVEN("malformed and nonflippable edge representations") { - auto triangulation = make_44_fixture().delaunay_snapshot(); + auto triangulation = make_44_fixture().delaunay_snapshot(); + auto foreign = make_44_fixture().delaunay_snapshot(); + auto const foreign_pivot = find_44_pivot(foreign); + REQUIRE(foreign_pivot.has_value()); + + std::optional boundary_edge; + for (auto edge = triangulation.finite_edges_begin(); + edge != triangulation.finite_edges_end(); ++edge) + { + if (!finite_incident_cells(triangulation, *edge)) + { + boundary_edge = *edge; + break; + } + } + REQUIRE(boundary_edge.has_value()); + Edge_handle const malformed{nullptr, 0, 1}; auto const before = canonical_triangulation(triangulation); WHEN("the checked edge helpers receive them") { CHECK_FALSE(ergodic_moves::get_incident_cells(triangulation, malformed)); + CHECK_FALSE( + ergodic_moves::get_incident_cells(triangulation, *boundary_edge)); + CHECK_FALSE(ergodic_moves::find_bistellar_flip_location(triangulation, + *foreign_pivot)); CHECK_FALSE(ergodic_moves::try_32_move(triangulation, malformed)); THEN("the triangulation remains unchanged") { CHECK_EQ(canonical_triangulation(triangulation), before); }