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 176 tests pass, 0 skipped (see `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    = mesh.euler_characteristic();
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 — Square-symmetric torus

**Setup.** Take a torus mesh with 4-fold rotational symmetry around the z-axis
(e.g. `code/data/off/torus_4x4.off`, which has M=4 columns of vertices).

**Expected.** The symmetry group Z₄ acts conformally.  Conformal automorphisms
of the torus correspond to SL(2,ℤ) symmetries of τ.  The unique fixed point of
a rotation of order 4 in the modular group is τ = i.  Therefore:

```
For a mesh with exact 4-fold symmetry and uniform edge lengths:
    Re(τ) = 0  (to machine precision, by symmetry)
    Im(τ) ≈ 1  (approaches 1 as mesh is refined)
```

The coarse 4×4 mesh (`torus_4x4.off`) gives Im(τ) in (0.7, 1.3) depending on
the 3D embedding (R=2, r=1 torus of revolution has unequal inner/outer edge lengths).
The uniformization algorithm finds the conformal class of the *abstract* metric
encoded in the edge lengths.

**Manual verification** (run from the build directory after adding a small
program or reading from the test output):

```cpp
ConformalMesh mesh; load_mesh(mesh, "code/data/off/torus_4x4.off");
EuclideanMaps maps = setup_euclidean_maps(mesh);
compute_euclidean_lambda0_from_mesh(mesh, maps);
enforce_gauss_bonnet(mesh, maps);
auto res = newton_euclidean(mesh, maps);
CutGraph cg = compute_cut_graph(mesh);
HolonomyData hol;
euclidean_layout(mesh, res.x, maps, &cg, &hol, true);
PeriodData pd = compute_period_matrix(hol);
// pd.tau_reduced satisfies the fundamental domain invariants above
```

---

## 4 — Hexagonal-symmetric torus

**Setup.** Take a torus mesh with 6-fold rotational symmetry
(`code/data/off/torus_hex_6x6.off`, M=6).

**Expected.** The unique τ fixed under a rotation of order 6 in SL(2,ℤ) is
τ = e^{iπ/3} = ½ + i√3/2.  So:

```
Re(τ) = 0.5  (to machine precision, by symmetry)
Im(τ) = √3/2 ≈ 0.8660
```

The coarse 6×6 torus of revolution approximates this: Re(τ) ≈ 0.5 by symmetry,
Im(τ) approaches √3/2 as the mesh is refined toward a flat hexagonal lattice.

---

## 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
auto res = newton_euclidean(mesh, maps);
EXPECT_LT(res.iterations, 30);
EXPECT_LT(res.gradient_norm, 1e-10);
```

Covered by: `cgal.EuclideanPipeline.ConvRates_*` and similar tests.

---

## 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)
```

These are the **holonomy consistency** checks implemented in `test_phase7.cpp`
(`cgal.HolonomyData.*`).

---

## 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` → 176 tests pass, 0 skipped
- [ ] `cgal.GaussBonnet.*` all pass → topology is correctly read from mesh
- [ ] `cgal.EuclideanFunctional.GradientCheck_*` pass → energy = integral of gradient
- [ ] `cgal.PeriodMatrix.TauInFundamentalDomain_*` pass → SL(2,ℤ) reduction correct
- [ ] `cgal.MobiusMap.Compose_*` and `Inverse_*` pass → Möbius arithmetic correct
- [ ] `cgal.HolonomyData.*` pass → holonomy loops close up

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