Merge pull request 'feat(geometry): pn_geometry.hpp — Pn projective-metric substrate (jReality port)' (#30) from feat/lawson-hyperideal-golden into main

code/tests/cgal/CMakeLists.txt

@@ -86,6 +86,10 @@ add_executable(conformallab_cgal_tests
     #    newton_inversive_distance (FD Hessian).
     test_newton_phase9a.cpp
 
+    # ── Tier-3 Java cross-validation: Lawson square-tiled HyperIdeal ─────────
+    #    Low-level half-edge genus-2 generator + golden-vector convergence.
+    test_lawson_hyperideal.cpp
+
     # ── Phase 8b-Lite: CGAL entry wrappers for the 4 non-Euclidean modes ─────
     #    Spherical, HyperIdeal, CircleP-Euclidean, Inversive-Distance via
     #    <CGAL/Discrete_*.h> public API + Conformal_layout.h wrapper.

code/tests/cgal/test_lawson_hyperideal.cpp

@@ -0,0 +1,257 @@
+// Copyright (c) 2024-2026 Tarik Moussa.
+// SPDX-License-Identifier: MIT
+
+// test_lawson_hyperideal.cpp
+//
+// Tier-3 Java cross-validation: HyperIdealConvergenceTest (Lawson square-tiled).
+//
+// Ports de.varylab.discreteconformal.functional.HyperIdealGenerator. The base is
+// a genus-2 surface with 4 vertices and 12 edges → MULTI-EDGES (≥2 edges per
+// vertex pair), so CGAL add_face/OFF/polygon-soup cannot build it; we use the
+// low-level Surface_mesh half-edge API directly.
+//
+// Variants ported (vs HyperIdealConvergenceTest golden vectors):
+//   1. createLawsonSquareTiled()              → diagonal triangulation
+//   2. createLawsonSquareTiledWithBranchPoints() → stellar subdivision + ideal
+//      branch-point centres
+//
+// NOT ported — createLawsonHyperelliptic(): loads `lawson_curve_source.xml`
+// (a Java conformal-data HalfedgeEmbedding, not OBJ/OFF) and derives θ via
+// HyperIdealHyperellipticUtility.calculateCircleIntersections. It needs a reader
+// for that XML format + a port of that utility, and its golden vector is not
+// class-symmetric (mixed ±values), so the symmetry shortcut does not apply.
+// Deferred as its own task; see doc/reviewer/java-ignore-crossvalidation.md.
+
+#include "conformal_mesh.hpp"
+#include "hyper_ideal_functional.hpp"
+#include "newton_solver.hpp"
+#include <CGAL/Polygon_mesh_processing/triangulate_faces.h>
+#include <CGAL/boost/graph/Euler_operations.h>
+#include <gtest/gtest.h>
+#include <array>
+#include <vector>
+#include <set>
+#include <cmath>
+
+using namespace conformallab;
+
+namespace {
+
+// Build the Lawson square-tiled BASE (4 vertices, 12 edges, 6 quad faces,
+// genus 2) via the low-level half-edge API — multi-edges (≥2 edges per vertex
+// pair) make add_face/OFF/polygon-soup impossible.  Returns the 4 original
+// vertices in `V`.  Not yet triangulated.
+void build_lawson_base(ConformalMesh& m, std::array<Vertex_index, 4>& V)
+{
+    V = {
+        m.add_vertex(Point3(0, 0, 0)),  // A
+        m.add_vertex(Point3(1, 0, 0)),  // B
+        m.add_vertex(Point3(0, 1, 0)),  // C
+        m.add_vertex(Point3(0, 0, 1))   // D
+    };
+
+    // Java half-edge target vertices (A=0,B=1,C=2,D=3), indices 0..23.
+    const int tgt[24] = {
+        0,1,3,2,  1,0,2,3,  3,2,0,1,  2,3,1,0,  0,1,3,2,  1,0,2,3
+    };
+    const int pairs[12][2] = {            // Java linkOppositeEdge
+        {0,6},{1,21},{2,4},{3,23},{5,11},{7,9},
+        {8,14},{10,12},{13,19},{15,17},{16,22},{18,20}
+    };
+    const int cyc[6][4] = {               // Java next-edge cycles → 6 quads
+        {0,1,2,3},{4,5,6,7},{8,9,10,11},{12,13,14,15},{16,17,18,19},{20,21,22,23}
+    };
+
+    std::array<Halfedge_index, 24> he;
+    for (const auto& p : pairs) {
+        Halfedge_index h = m.add_edge();
+        he[static_cast<std::size_t>(p[0])] = h;
+        he[static_cast<std::size_t>(p[1])] = m.opposite(h);
+    }
+    for (int j = 0; j < 24; ++j)
+        m.set_target(he[static_cast<std::size_t>(j)], V[static_cast<std::size_t>(tgt[j])]);
+    for (const auto& c : cyc) {
+        Face_index f = m.add_face();
+        m.set_halfedge(f, he[static_cast<std::size_t>(c[0])]);
+        for (int k = 0; k < 4; ++k) {
+            m.set_face(he[static_cast<std::size_t>(c[k])], f);
+            m.set_next(he[static_cast<std::size_t>(c[k])],
+                       he[static_cast<std::size_t>(c[(k + 1) % 4])]);
+        }
+    }
+    for (int j = 0; j < 24; ++j)
+        m.set_halfedge(V[static_cast<std::size_t>(tgt[j])], he[static_cast<std::size_t>(j)]);
+}
+
+// Plain Lawson: base + diagonal triangulation (6 quads → 12 triangles).
+// `original_edges` receives the 12 base edges (the 6 diagonals are "aux").
+ConformalMesh make_lawson_square_tiled(std::set<Edge_index>* original_edges = nullptr)
+{
+    ConformalMesh m;
+    std::array<Vertex_index, 4> V;
+    build_lawson_base(m, V);
+
+    if (original_edges) {
+        original_edges->clear();
+        for (auto e : m.edges()) original_edges->insert(e);
+    }
+    CGAL::Polygon_mesh_processing::triangulate_faces(m);
+    return m;
+}
+
+// Branch-points Lawson: base + STELLAR subdivision (Java StellarLinear) — a
+// center vertex per quad fan-connected to its 4 corners (6 ideal "branch"
+// centers, 24 triangles, 36 edges).  `orig_verts` = the 4 real vertices;
+// `base_edges` = the 12 base edges (θ=π); the 24 spokes are θ=π/2.
+ConformalMesh make_lawson_branch_points(std::array<Vertex_index, 4>& orig_verts,
+                                        std::set<Edge_index>&        base_edges)
+{
+    ConformalMesh m;
+    build_lawson_base(m, orig_verts);
+
+    base_edges.clear();
+    for (auto e : m.edges()) base_edges.insert(e);
+
+    // Stellar-subdivide each of the 6 quads (collect halfedges first; the loop
+    // mutates the mesh).
+    std::vector<Halfedge_index> face_he;
+    for (auto f : m.faces()) face_he.push_back(m.halfedge(f));
+    for (auto h : face_he)
+        CGAL::Euler::add_center_vertex(h, m);
+
+    return m;
+}
+
+} // namespace
+
+TEST(LawsonHyperIdeal, BuildsValidGenus2Mesh)
+{
+    std::set<Edge_index> original;
+    ConformalMesh m = make_lawson_square_tiled(&original);
+
+    EXPECT_TRUE(m.is_valid(false)) << "Lawson square-tiled mesh is not a valid halfedge structure";
+    EXPECT_EQ(m.number_of_vertices(), 4u);
+    EXPECT_EQ(m.number_of_faces(),    12u);   // 6 quads → 12 triangles
+    EXPECT_EQ(m.number_of_edges(),    18u);   // 12 original + 6 diagonals
+    EXPECT_EQ(original.size(),         12u);
+
+    // Euler characteristic χ = V − E + F = 4 − 18 + 12 = −2 ⇒ genus 2.
+    const int chi = static_cast<int>(m.number_of_vertices())
+                  - static_cast<int>(m.number_of_edges())
+                  + static_cast<int>(m.number_of_faces());
+    EXPECT_EQ(chi, -2) << "expected genus 2 (χ = −2), got χ = " << chi;
+}
+
+// ════════════════════════════════════════════════════════════════════════════
+// Golden-vector cross-validation against Java HyperIdealConvergenceTest
+//
+// Java sets Θ_v = 2π (vertices), θ_e = π/2 (the 12 original edges) and θ_e = π
+// (the 6 triangulation/aux edges), then solves with TAO/BLMVM and asserts the
+// converged solution vector.  The solution is perfectly symmetric:
+//     vertices (×4)            → 1.1462158341786262
+//     original edges (×12)     → 1.7627471737467797
+//     aux/diagonal edges (×6)  → 2.633915794495759
+// We assert membership in these three classes (robust to DOF ordering, which
+// differs between Java and CGAL).
+// ════════════════════════════════════════════════════════════════════════════
+TEST(LawsonHyperIdeal, ConvergenceGoldenVector_JavaXVal)
+{
+    std::set<Edge_index> original;
+    ConformalMesh m = make_lawson_square_tiled(&original);
+
+    HyperIdealMaps maps = setup_hyper_ideal_maps(m);     // Θ_v=2π, θ_e=π
+    const int n = assign_all_dof_indices(m, maps);       // all vertices + edges
+    ASSERT_EQ(n, 4 + 18);                                // 4 b + 18 a = 22 DOFs
+
+    // θ_e = π/2 for the 12 original edges; the 6 diagonals keep the default π.
+    for (auto e : m.edges())
+        if (original.count(e)) maps.theta_e[e] = PI / 2.0;
+
+    // Solve from a positive interior point (HyperIdeal variables b,a > 0).
+    std::vector<double> x0(static_cast<std::size_t>(n), 1.0);
+    auto res = newton_hyper_ideal(m, x0, maps, /*tol=*/1e-10, /*max_iter=*/200);
+    ASSERT_TRUE(res.converged)
+        << "HyperIdeal Newton did not converge; ||G||=" << res.grad_inf_norm;
+
+    constexpr double b_gold      = 1.1462158341786262;
+    constexpr double a_orig_gold = 1.7627471737467797;
+    constexpr double a_aux_gold  = 2.633915794495759;
+    const double tol = 1e-5;
+
+    for (auto v : m.vertices()) {
+        const int iv = maps.v_idx[v];
+        ASSERT_GE(iv, 0);
+        EXPECT_NEAR(res.x[static_cast<std::size_t>(iv)], b_gold, tol)
+            << "vertex DOF " << iv << " off golden b";
+    }
+    for (auto e : m.edges()) {
+        const int ie = maps.e_idx[e];
+        ASSERT_GE(ie, 0);
+        const double gold = original.count(e) ? a_orig_gold : a_aux_gold;
+        EXPECT_NEAR(res.x[static_cast<std::size_t>(ie)], gold, tol)
+            << "edge DOF " << ie << " off golden a ("
+            << (original.count(e) ? "original π/2" : "aux π") << ")";
+    }
+}
+
+// ════════════════════════════════════════════════════════════════════════════
+// Branch-points variant — Java HyperIdealConvergenceTest...WithBranchPoints
+//
+// Java applies StellarLinear (stellar subdivision: a center vertex per quad) to
+// the base, makes the 4 original vertices variable and the 6 centers IDEAL
+// (b = 0, solver index −1), with Θ_v = 2π, θ_e = π on the 12 base edges and
+// θ_e = π/2 on the 24 spokes.  The converged solution (symmetric per class):
+//     original vertices (×4) → 1.3169579…
+//     base edges (×12)       → 2.2924317…
+//     spoke edges (×24)      → 0          (the π/2 spokes collapse to ideal)
+// ════════════════════════════════════════════════════════════════════════════
+TEST(LawsonHyperIdeal, BranchPointsGoldenVector_JavaXVal)
+{
+    std::array<Vertex_index, 4> orig;
+    std::set<Edge_index>        base;
+    ConformalMesh m = make_lawson_branch_points(orig, base);
+
+    EXPECT_TRUE(m.is_valid(false));
+    EXPECT_EQ(m.number_of_vertices(), 10u);   // 4 original + 6 stellar centers
+    EXPECT_EQ(m.number_of_faces(),    24u);   // 6 quads × 4
+    EXPECT_EQ(m.number_of_edges(),    36u);   // 12 base + 24 spokes
+    EXPECT_EQ(base.size(),            12u);
+
+    HyperIdealMaps maps = setup_hyper_ideal_maps(m);   // Θ_v=2π, θ_e=π
+    const std::set<Vertex_index> origset(orig.begin(), orig.end());
+
+    // 4 original vertices variable; 6 centers ideal (−1); all 36 edges variable.
+    int idx = 0;
+    for (auto v : m.vertices()) maps.v_idx[v] = origset.count(v) ? idx++ : -1;
+    for (auto e : m.edges())    maps.e_idx[e] = idx++;
+    const int n = idx;
+    ASSERT_EQ(n, 4 + 36);
+
+    // θ_e = π/2 on the 24 spokes (base edges keep the default π).
+    for (auto e : m.edges())
+        if (!base.count(e)) maps.theta_e[e] = PI / 2.0;
+
+    std::vector<double> x0(static_cast<std::size_t>(n), 1.0);
+    auto res = newton_hyper_ideal(m, x0, maps, /*tol=*/1e-10, /*max_iter=*/300);
+    ASSERT_TRUE(res.converged)
+        << "branch-points Newton did not converge; ||G||=" << res.grad_inf_norm;
+
+    constexpr double b_gold       = 1.3169579;
+    constexpr double a_base_gold  = 2.2924317;
+    constexpr double a_spoke_gold = 0.0;
+    const double tol = 1e-4;   // golden per-class spread is ~1e-7
+
+    for (auto v : m.vertices()) {
+        const int iv = maps.v_idx[v];
+        if (iv < 0) continue;                     // ideal centers
+        EXPECT_NEAR(res.x[static_cast<std::size_t>(iv)], b_gold, tol)
+            << "branch-points vertex DOF " << iv << " off golden b";
+    }
+    for (auto e : m.edges()) {
+        const int ie = maps.e_idx[e];
+        const double gold = base.count(e) ? a_base_gold : a_spoke_gold;
+        EXPECT_NEAR(res.x[static_cast<std::size_t>(ie)], gold, tol)
+            << "branch-points edge DOF " << ie << " off golden a ("
+            << (base.count(e) ? "base π" : "spoke π/2") << ")";
+    }
+}