ConformalLabpp/doc/math/validation.md

# Mathematical Validation

This document lists analytically known results and explains how to verify
them against conformallab++ output.  It is the primary tool for an independent
mathematician to check the correctness of the implementation.

---

## How to run the examples

```bash
cmake -S code -B build -DWITH_CGAL=ON -DCMAKE_BUILD_TYPE=Release
cmake --build build --target conformallab_cgal_tests
ctest --test-dir build -R cgal --output-on-failure
```

All tests pass, 0 skipped (per-suite breakdown: `doc/api/tests.md`).

---

## 1 — Gauss–Bonnet (topology)

**Theorem.** For any closed triangulated surface M,

```
Σᵥ (2π − Θᵥ) = 2π · χ(M)
```

where χ(M) = 2 − 2g is the Euler characteristic.

| Surface | g | χ | Σ(2π − Θᵥ) |
|---|---|---|---|
| Sphere (tetrahedron, cube, …) | 0 | 2 | 4π |
| Torus | 1 | 0 | 0 |
| Double torus | 2 | −2 | −4π |

**How to check:**

```cpp
#include "gauss_bonnet.hpp"
auto defect = gauss_bonnet_sum(mesh, maps);   // Σ(2π − Θᵥ)
auto chi    = euler_characteristic(mesh);     // free function, not a member
EXPECT_NEAR(defect, 2.0 * M_PI * chi, 1e-10);
```

Covered by: `cgal.GaussBonnet.*` tests in `test_phase6.cpp`.

---

## 2 — Period matrix: fundamental domain invariants

**Theorem (SL(2,ℤ)-reduction).** Every lattice τ ∈ ℍ has a unique representative
in the standard fundamental domain

```
F = { τ ∈ ℍ :  |τ| ≥ 1,  |Re(τ)| ≤ 1/2,  Im(τ) > 0 }
```

**After calling `compute_period_matrix(hol)`, the returned τ must satisfy:**

| Condition | Invariant |
|---|---|
| `pd.tau_reduced.imag() > 0` | τ lies in the upper half-plane |
| `std::abs(pd.tau_reduced) >= 1.0 - 1e-10` | τ outside unit disk |
| `std::abs(pd.tau_reduced.real()) <= 0.5 + 1e-10` | τ in vertical strip |

These three conditions hold for **any** closed genus-1 triangulated surface
processed through Euclidean uniformization — they are topology, not geometry.

Covered by: `cgal.PeriodMatrix.TauInFundamentalDomain_*` tests in `test_phase7.cpp`.

---

## 3 — Torus of revolution: conformal modulus

**Setup.** The bundled torus meshes are surfaces of **revolution** (a tube of
minor radius `r` swept around a major circle of radius `R > r`), *not* abstract
square/hexagonal flat tori:

| mesh | major R | minor r | cross section |
|---|---|---|---|
| `torus_4x4.off`     | 2 | 1 | square (4-gon)   |
| `torus_hex_6x6.off` | 3 | 1 | hexagon (6-gon)  |
| `torus_8x8.off`     | 3 | 1 | octagon (8-gon)  |

**Expected.** The induced metric `ds² = (R + r cos φ)² dθ² + r² dφ²` is made
flat by the conformal change of variable `dψ = r/(R + r cos φ) dφ`.  The
ψ-period is `∮ r/(R + r cos φ) dφ = 2πr/√(R²−r²)`, so the flat torus is the
rectangular lattice `2π·ℤ × (2πr/√(R²−r²))·ℤ`.  Its period ratio, reduced so
that `|τ| ≥ 1`, is therefore **purely imaginary**:

```
Re(τ) = 0                  (meridian ⟂ longitude reflection symmetry)
Im(τ) = √(R² − r²) / r     (reduced conformal modulus)
```

giving the analytic targets

| mesh | analytic τ = i·√(R²−r²)/r |
|---|---|
| `torus_4x4.off`     | i·√3 ≈ **1.732 i** |
| `torus_hex_6x6.off` | i·√8 ≈ **2.828 i** |
| `torus_8x8.off`     | i·√8 ≈ **2.828 i** |

> A common pitfall is to expect τ = i (or the order-6 fixed point e^{iπ/3})
> from the 4-fold (6-fold) symmetry.  That reasoning is **wrong** here: a torus
> of revolution's rotational symmetry is a rotation about the axis, which acts
> on the surface as a *fixed-point-free* translation along the longitude — it is
> not an order-4 conformal automorphism with a fixed point, so it does not pin τ
> to a modular fixed point.  The correct invariant is the rectangular modulus
> above.

**Reproduced end-to-end** (solve → `compute_cut_graph` → `euclidean_layout(…,
&cg, &hol)` → `compute_period_matrix`).  The coarse polygonal cross sections
approximate the circular modulus from above; the gap shrinks as the cross
section gains sides:

| mesh | analytic | computed τ | rel. error |
|---|---|---|---|
| `torus_4x4.off`     | 1.732 i | ≈ 1.79 i | ~3 % (4-gon) |
| `torus_hex_6x6.off` | 2.828 i | ≈ 2.85 i | ~0.8 % (6-gon) |
| `torus_8x8.off`     | 2.828 i | ≈ 2.84 i | ~0.4 % (8-gon) |

Covered by: `cgal.HolonomyEndToEnd.Torus*_TauMatchesRevolutionModulus` in
`test_phase7.cpp`.

**Conformal flattening** of the torus is wired end-to-end and converges in a
handful of Newton steps (3–4 on the bundled meshes):

```bash
./bin/conformallab_core -i code/data/off/torus_4x4.off -g euclidean -v -o lay.off
# → topology: closed, free DOFs=15, genus=1
# → Euclidean: converged=yes  iter=3  |grad|_inf≈1e-12
```

Equivalent C++ (current API — note `newton_euclidean` takes an `x0` vector and
the DOF indices must be assigned first):

```cpp
ConformalMesh mesh = load_mesh("code/data/off/torus_4x4.off");
EuclideanMaps  maps = setup_euclidean_maps(mesh);   // Θ_v = 2π (flat target)
compute_euclidean_lambda0_from_mesh(mesh, maps);

// Pin one vertex (scale gauge), free the rest; make the flat target
// Gauss-Bonnet-consistent.
int idx = 0; bool pinned = false;
for (auto v : mesh.vertices())
    maps.v_idx[v] = (!pinned ? (pinned = true, -1) : idx++);
enforce_gauss_bonnet(mesh, maps);

std::vector<double> x0(idx, 0.0);
auto res = newton_euclidean(mesh, x0, maps);        // converged after ~3 iters
```

The CLI reports the period ratio for genus-1 inputs, e.g.

```bash
./bin/conformallab_core -i code/data/off/torus_4x4.off -g euclidean -v
# → period ratio τ = 0.000000 + 1.793... i  (genus 1, reduced to fundamental domain)
```

---

## 4 — Higher-resolution cross sections converge to the circular modulus

`torus_hex_6x6.off` (R=3, r=1, hexagon) and `torus_8x8.off` (R=3, r=1, octagon)
share the same analytic modulus `i·√8 ≈ 2.828 i` (see §3).  Because both
approximate the *circular* tube cross section, the computed τ approaches the
analytic value as the polygon gains sides: the 8-gon (≈ 2.84 i, ~0.4 %) is
closer than the 6-gon (≈ 2.85 i, ~0.8 %), which is closer than the 4-gon of
`torus_4x4.off` (~3 %).  All three are asserted in
`cgal.HolonomyEndToEnd.Torus*_TauMatchesRevolutionModulus`.

---

## 5 — Newton convergence rate

**Theorem.** Because the Euclidean and hyper-ideal energies are strictly convex
(after gauge-fixing), Newton's method converges quadratically near the optimum.

**Expected:** for any mesh with up to a few hundred faces, Newton converges in
**fewer than 30 iterations** starting from u = 0.

```cpp
std::vector<double> x0(n_dofs, 0.0);
auto res = newton_euclidean(mesh, x0, maps);
EXPECT_LT(res.iterations, 30);
EXPECT_LT(res.grad_inf_norm, 1e-8);   // field is grad_inf_norm, default tol 1e-8
```

Covered by: `cgal.NewtonSolver.*` (Euclidean ×3, Spherical ×4, HyperIdeal ×4)
and `cgal.NewtonPhase9a.*` (the two circle-packing solvers).

---

## 6 — Gradient check (finite differences)

For each functional F(u), the gradient G = ∂F/∂u is verified by:

```
|G(u)ᵢ  −  (F(u + εeᵢ) − F(u − εeᵢ)) / (2ε)| < 1e-6
```

with ε = 1e-5.  This check is run **inside the test suite** for all three
geometries (Euclidean, Spherical, HyperIdeal) at u = 0 and at random u.

Relevant test suites:
```
cgal.EuclideanFunctional.GradientCheck_*
cgal.SphericalFunctional.GradientCheck_*
cgal.HyperIdealFunctional.GradientCheck_*
```

A failing gradient check means the energy and its derivative are inconsistent —
the Newton solver will converge to the wrong point.

---

## 7 — Holonomy composition (Möbius maps)

For a closed surface, the composition of holonomies around any contractible cycle
must be the identity.  In genus 1 with a single handle:

```
T₁ · T₂ · T₁⁻¹ · T₂⁻¹ = Id    (commutator = Id for a torus)
```

because π₁(T²) = ℤ × ℤ is abelian.

For genus g ≥ 2, the fundamental group is non-abelian and this check does not hold,
but the representation ρ: π₁(Σ_g) → SU(1,1) must still satisfy the relation

```
[T₁, T₂] · [T₃, T₄] · … = Id    (product of g commutators = Id)
```

The Möbius arithmetic these checks rely on (identity, inverse, composition,
`from_three`) is verified in `test_phase7.cpp` under `cgal.MobiusMap.*`.  The
Euclidean end-to-end holonomy extraction is now validated against the analytic
torus-of-revolution modulus (§3, `cgal.HolonomyEndToEnd.*`).  A standalone
`cgal.HolonomyData.*` commutator-closes-up suite for the *hyperbolic* (Möbius)
holonomy is still future work.

---

## 9 — Cross-validation with geometry-central *(optional / hypothetical)*

> **Note:** This section describes a possible external cross-validation that is not
> a prerequisite for the correctness of the implementation.
> It is of interest because geometry-central implements the same mathematical core
> (Gillespie, Springborn, Crane — SIGGRAPH 2021, building on
> Springborn 2020), but with a different algorithmic strategy
> (Ptolemaic flips + intrinsic triangulations instead of Newton on the
> original triangulation).

### Which outputs are comparable?

| Output | conformallab++ | geometry-central | Comparable? |
|---|---|---|---|
| u-vector (scale parameters) | `res.x` | `u` after Yamabe flow | ✓ after normalisation |
| UV coordinates | `layout.uv[v]` | conformal parameterisation | ✓ up to Möbius transformation |
| Gauss-Bonnet deficit | `gauss_bonnet_sum()` | implicit via curvature flow | ✓ (analytically identical) |
| Number of Newton iterations | `res.iterations` | Yamabe steps | ~ (different algorithm) |
| Period matrix τ | `pd.tau_reduced` | **not available** | ✗ |
| Möbius holonomy | `hol.T_a, T_b` | **not available** | ✗ |

### Normalisation alignment

The u-vector in conformallab++ has one degree of freedom (global additive constant —
gauge freedom after pin-fixing). geometry-central may use a different convention.
Normalise before comparing:

```cpp
// conformallab++: centre u
double mean_u = std::accumulate(x.begin(), x.end(), 0.0) / x.size();
std::vector<double> x_norm(x.size());
for (int i = 0; i < x.size(); ++i) x_norm[i] = x[i] - mean_u;

// Then compare with the geometry-central u-vector (also centred):
// max|x_norm[i] - gc_u[i]| < 1e-8  →  identical convergence point
```

### When is the comparison useful?

| Point in time | What is possible |
|---|---|
| **Now (Phase 7)** | Manual comparison using the same `.off`/`.obj` test meshes |
| **After Phase 8** | Automated comparison script (Python or separate C++ binary) |
| **Phase 10 (research)** | Algorithm comparison: Newton vs. Ptolemaic flips on difficult meshes |

### Connection to the literature

The Springborn 2020 paper ("Ideal Hyperbolic Polyhedra and Discrete Uniformization")
is **already implemented in conformallab++** — it is the mathematical foundation
for the HyperIdeal geometry mode (Phase 2/3). The geometry-central implementation
is based on the extension by Gillespie, Springborn & Crane (2021), which uses the
same variational principle of Bobenko–Springborn 2004 but additionally applies
Ptolemaic flips to improve the triangulation during optimisation — an idea not yet
implemented in conformallab++ (→ GC-2 in the phase roadmap).

---

## 8 — Checklist for an independent reviewer

Run these in order to validate the implementation:

- [ ] `ctest --test-dir build -R cgal --output-on-failure` → all pass, 0 skipped (count: `doc/api/tests.md`)
- [ ] `cgal.GaussBonnet.*` all pass → topology is correctly read from mesh
- [ ] `cgal.EuclideanFunctional.GradientCheck_*` pass → energy = integral of gradient
- [ ] `cgal.PeriodMatrix.*` pass → SL(2,ℤ) reduction correct (on prescribed holonomy)
- [ ] `cgal.MobiusMap.*` pass → Möbius arithmetic (identity, inverse, compose) correct

All of the above are **deterministic, analytic tests** — no mesh loading, no
file I/O, no floating-point non-determinism beyond standard IEEE-754.