Reapply 3+1D CDT support

docs/4d-standard-cdt-candidate.md

@@ -0,0 +1,78 @@
+# 4D Standard-CDT Candidate Validation
+
+This implementation treats `cdt -d4` as a standard-CDT candidate only when the
+state passes the validator in `Foliated_triangulation_4.hpp`. The validator
+requires periodic time, closed `S3` spatial-slice metadata, non-negative simplex
+counts, non-negative proposal multiplicities, valid reciprocal gluing for
+explicit combinatorial simplices, and exact agreement between `N4` and the sum
+of `N41`, `N32`, `N23`, and `N14`.
+
+The action convention follows the usual 4D CDT bare-coupling form
+
+```text
+S_E = -(kappa_0 + 6 Delta) N0
+      + kappa_4 N4
+      + Delta (2 N41_total + N32_total)
+      + epsilon (N4 - target_N4)^2
+```
+
+where `N41_total = N41 + N14` and `N32_total = N32 + N23`.
+The numerical setup follows the standard CDT requirements that Monte Carlo
+moves preserve fixed topology, causality and detailed balance; see the CDT
+transfer-matrix lecture notes for the seven-move 4D setup and for the
+`cos^3` de Sitter spatial-volume profile and `N4^(1/4)`, `N4^(3/4)` scaling.
+
+Sources:
+
+- https://indico.tpi.uni-jena.de/event/76/attachments/38/132/Goerlich_-_Causal_Dynamical_Triangulations.pdf
+- https://www.scholarpedia.org/article/Causal_Dynamical_Triangulation
+
+## Move Catalogue
+
+The source of truth is `Move_catalog_4.hpp`. Tests compare directly against this
+table so documentation and runtime behavior cannot silently diverge.
+
+| Move | Inverse | Proposal multiplicity | Delta `(N0,N1,N2,N3,N4,N41,N32,N23,N14)` |
+| --- | --- | --- | --- |
+| `TWO_FOUR` | `FOUR_TWO` | legal spacelike tetrahedra | `(0,+1,+4,+5,+2,+1,+1,0,0)` |
+| `FOUR_TWO` | `TWO_FOUR` | removable timelike edges | `(0,-1,-4,-5,-2,-1,-1,0,0)` |
+| `THREE_THREE` | `THREE_THREE` | legal mixed triangles | `(0,0,0,0,0,0,-1,+1,0)` or the time-reversed sign |
+| `FOUR_SIX` | `SIX_FOUR` | legal mixed triangles | `(0,+1,+3,+4,+2,0,+1,+1,0)` |
+| `SIX_FOUR` | `FOUR_SIX` | removable order-six local stars | `(0,-1,-3,-4,-2,0,-1,-1,0)` |
+| `TWO_EIGHT` | `EIGHT_TWO` | legal spacelike tetrahedra | `(+1,+6,+10,+10,+6,+2,+1,+1,+2)` |
+| `EIGHT_TWO` | `TWO_EIGHT` | removable spatial vertices | `(-1,-6,-10,-10,-6,-2,-1,-1,-2)` |
+
+## Detailed Balance
+
+For a proposed transition `A -> B`, the sampler uses
+
+```text
+acceptance = min(1, exp(-(S(B)-S(A))) * q(B -> A) / q(A -> B))
+```
+
+where `q` is computed from the current proposal inventory, not from historical
+move frequencies. The enumerable detailed-balance tests build a small reachable
+state graph, verify reverse transitions, and compare both sides of
+
+```text
+pi(A) q(A -> B) A(A -> B) = pi(B) q(B -> A) A(B -> A)
+```
+
+with `pi(T) = exp(-S(T))`.
+
+## Phase Diagnostics
+
+The single-run summary reports the conservative profile diagnostic for the
+measurements available in that run. A full `c_ds_supported` finite-size claim
+uses `diagnose_c_ds_finite_size()` and requires all of the following:
+
+- centered ensemble profile is better fit by `cos^3` than collapsed or
+  alternating-slice alternatives;
+- finite-size width scales as `N4^(1/4)`;
+- finite-size peak volume scales as `N4^(3/4)`;
+- covariance/effective-action extraction has enough decorrelated samples;
+- C_b diagnostics are not triggered.
+
+Synthetic tests exercise the analysis layer. Production claims still require
+real independent chains and decorrelated measurements written under
+`results/<run-id>/`.

include/Detailed_balance_4.hpp

@@ -0,0 +1,243 @@
+/*******************************************************************************
+ Causal Dynamical Triangulations in C++ using CGAL
+*******************************************************************************/
+
+/// @file Detailed_balance_4.hpp
+/// @brief Small-ensemble detailed-balance verifier for 4D CDT candidates.
+
+#ifndef CDT_PLUSPLUS_DETAILED_BALANCE_4_HPP
+#define CDT_PLUSPLUS_DETAILED_BALANCE_4_HPP
+
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <limits>
+#include <map>
+#include <queue>
+#include <string>
+#include <vector>
+
+#include "Ergodic_moves_4.hpp"
+
+namespace cdt::four_d
+{
+  struct DetailedBalanceEdge4D
+  {
+    std::string              from_hash;
+    std::string              to_hash;
+    move_tracker::MoveType4D move{move_tracker::MoveType4D::NO_MOVE};
+    long double              lhs{0.0L};
+    long double              rhs{0.0L};
+    long double              residual{0.0L};
+  };
+
+  struct DetailedBalanceReport4D
+  {
+    bool                               passed{true};
+    std::vector<DetailedBalanceEdge4D> edges;
+    std::vector<std::string>           errors;
+  };
+
+  [[nodiscard]] inline auto proposal_probability(
+      FoliatedTriangulation4 const&  triangulation,
+      move_tracker::MoveType4D const move) -> long double
+  {
+    auto const multiplicity = triangulation.candidate_multiplicity(move);
+    if (multiplicity <= 0) { return 0.0L; }
+    return 1.0L / (static_cast<long double>(move_tracker::NUMBER_OF_4D_MOVES) *
+                   static_cast<long double>(multiplicity));
+  }
+
+  [[nodiscard]] inline auto acceptance_probability(
+      FoliatedTriangulation4 const& before, FoliatedTriangulation4 const& after,
+      move_tracker::MoveType4D const move, S4Couplings const& couplings)
+      -> long double
+  {
+    auto const reverse   = move_tracker::reverse_move(move);
+    auto const forward_q = proposal_probability(before, move);
+    auto const reverse_q = proposal_probability(after, reverse);
+    if (forward_q == 0.0L || reverse_q == 0.0L) { return 0.0L; }
+    auto const delta =
+        S4_action_difference(before.counts(), after.counts(), couplings);
+    if (!std::isfinite(delta))
+    {
+      return std::numeric_limits<long double>::quiet_NaN();
+    }
+    auto const log_ratio = -delta + std::log(reverse_q) - std::log(forward_q);
+    if (!std::isfinite(log_ratio))
+    {
+      return std::numeric_limits<long double>::quiet_NaN();
+    }
+    if (log_ratio >= 0.0L) { return 1.0L; }
+    auto const probability = std::exp(log_ratio);
+    return std::isfinite(probability)
+             ? probability
+             : std::numeric_limits<long double>::quiet_NaN();
+  }
+
+  [[nodiscard]] inline auto boltzmann_weight(
+      FoliatedTriangulation4 const& triangulation, S4Couplings const& couplings)
+      -> long double
+  {
+    return std::exp(-S4_bulk_action(triangulation.counts(), couplings));
+  }
+
+  [[nodiscard]] inline auto log_boltzmann_weight(
+      FoliatedTriangulation4 const& triangulation, S4Couplings const& couplings)
+      -> long double
+  {
+    return -S4_bulk_action(triangulation.counts(), couplings);
+  }
+
+  [[nodiscard]] inline auto verify_detailed_balance(
+      FoliatedTriangulation4 const& seed, S4Couplings const& couplings,
+      int const max_depth, long double const tolerance = 1.0e-10L,
+      std::size_t const max_states = 1024) -> DetailedBalanceReport4D
+  {
+    DetailedBalanceReport4D report;
+    if (max_depth <= 0)
+    {
+      report.passed = false;
+      report.errors.emplace_back(
+          "Detailed-balance enumeration max_depth must be positive.");
+      return report;
+    }
+    if (max_states == 0)
+    {
+      report.passed = false;
+      report.errors.emplace_back(
+          "Detailed-balance enumeration max_states must be positive.");
+      return report;
+    }
+
+    std::map<std::string, FoliatedTriangulation4>      states;
+    std::map<std::string, int>                         depths;
+    std::queue<std::pair<FoliatedTriangulation4, int>> frontier;
+
+    states.emplace(seed.canonical_hash(), seed);
+    depths.emplace(seed.canonical_hash(), 0);
+    frontier.emplace(seed, 0);
+
+    auto cap_reached = false;
+    while (!frontier.empty() && !cap_reached)
+    {
+      auto [state, depth] = frontier.front();
+      frontier.pop();
+      if (depth >= max_depth) { continue; }
+      for (auto const descriptor : all_move_descriptors_4d())
+      {
+        auto moved = moves::apply(state, descriptor.move);
+        if (!moved) { continue; }
+        auto const hash = moved->triangulation.canonical_hash();
+        if (!states.contains(hash))
+        {
+          if (states.size() >= max_states)
+          {
+            report.passed = false;
+            report.errors.emplace_back(
+                "Detailed-balance enumeration reached max_states.");
+            cap_reached = true;
+            break;
+          }
+          states.emplace(hash, moved->triangulation);
+          depths.emplace(hash, depth + 1);
+          frontier.emplace(moved->triangulation, depth + 1);
+        }
+      }
+    }
+    if (cap_reached) { return report; }
+
+    for (auto const& [from_hash, from_state] : states)
+    {
+      auto const from_depth = depths.at(from_hash);
+      if (from_depth >= max_depth) { continue; }
+      for (auto const descriptor : all_move_descriptors_4d())
+      {
+        auto moved = moves::apply(from_state, descriptor.move);
+        if (!moved) { continue; }
+        auto const to_hash = moved->triangulation.canonical_hash();
+        auto const to_it   = states.find(to_hash);
+        if (to_it == states.end())
+        {
+          report.passed = false;
+          report.errors.emplace_back(
+              "Reachable transition escaped enumeration depth.");
+          continue;
+        }
+        auto const& to_state = to_it->second;
+        auto        reverse  = moves::apply(to_state, descriptor.inverse);
+        if (!reverse || reverse->triangulation.canonical_hash() != from_hash)
+        {
+          report.passed = false;
+          report.errors.emplace_back("Missing reverse transition.");
+          continue;
+        }
+
+        auto const forward_q =
+            proposal_probability(from_state, descriptor.move);
+        auto const reverse_q =
+            proposal_probability(to_state, descriptor.inverse);
+        auto const forward_acceptance = acceptance_probability(
+            from_state, to_state, descriptor.move, couplings);
+        auto const reverse_acceptance = acceptance_probability(
+            to_state, from_state, descriptor.inverse, couplings);
+        auto const from_log_weight =
+            log_boltzmann_weight(from_state, couplings);
+        auto const to_log_weight = log_boltzmann_weight(to_state, couplings);
+        auto const min_log = std::log(std::numeric_limits<long double>::min());
+        auto const max_log = std::log(std::numeric_limits<long double>::max());
+        if (!std::isfinite(from_log_weight) || !std::isfinite(to_log_weight) ||
+            from_log_weight <= min_log || to_log_weight <= min_log ||
+            from_log_weight >= max_log || to_log_weight >= max_log ||
+            forward_q <= 0.0L || reverse_q <= 0.0L ||
+            forward_acceptance <= 0.0L || reverse_acceptance <= 0.0L ||
+            !std::isfinite(forward_acceptance) ||
+            !std::isfinite(reverse_acceptance))
+        {
+          report.passed = false;
+          report.errors.emplace_back(
+              "Detailed-balance transition has non-finite or underflowed "
+              "weight.");
+          continue;
+        }
+
+        // Both sides use the same Metropolis-Hastings acceptance rule as the
+        // sampler, so this check is an algebraic detailed-balance identity
+        // evaluated in log space to avoid overflow and silent underflow.
+        auto const log_lhs = from_log_weight + std::log(forward_q) +
+                             std::log(forward_acceptance);
+        auto const log_rhs =
+            to_log_weight + std::log(reverse_q) + std::log(reverse_acceptance);
+        if (!std::isfinite(log_lhs) || !std::isfinite(log_rhs))
+        {
+          report.passed = false;
+          report.errors.emplace_back(
+              "Detailed-balance transition has non-finite log weight.");
+          continue;
+        }
+        if (log_lhs <= min_log || log_rhs <= min_log || log_lhs >= max_log ||
+            log_rhs >= max_log)
+        {
+          report.passed = false;
+          report.errors.emplace_back(
+              "Detailed-balance transition weight cannot be represented "
+              "without underflow or overflow.");
+          continue;
+        }
+        auto const residual = std::abs(log_lhs - log_rhs);
+        report.edges.push_back(DetailedBalanceEdge4D{
+            from_hash, to_hash, descriptor.move, std::exp(log_lhs),
+            std::exp(log_rhs), residual});
+        if (residual > tolerance)
+        {
+          report.passed = false;
+          report.errors.emplace_back(
+              "Detailed-balance log residual exceeds tolerance.");
+        }
+      }
+    }
+    return report;
+  }
+}  // namespace cdt::four_d
+
+#endif  // CDT_PLUSPLUS_DETAILED_BALANCE_4_HPP