ConformalLabpp/doc/architecture/phase-9a-validation.md

Phase 9a Validation Report

This document records the mathematical and code-level validation of the two Phase 9a functionals against their reference sources. It is produced as part of the Phase 9a acceptance procedure (decided 2026-05-19).

Functional	Reference	Java original	Lines
9a.1 cp_euclidean_functional.hpp	Bobenko-Pinkall-Springborn 2010	CPEuclideanFunctional.java (260)	320
9a.2 inversive_distance_functional.hpp	Luo 2004 + Glickenstein 2011	none	290

---

1. CP-Euclidean (9a.1) — line-by-line Java mapping

Reference file: /Users/tarikmoussa/Desktop/conformallab/src/de/varylab/discreteconformal/functional/CPEuclideanFunctional.java

Java line(s)	Java construct	C++ counterpart	Status
55–58	public class CPEuclideanFunctional<V,E,F>	namespace conformallab (header-only, monomorphic on ConformalMesh)	✓
61–69	thetaMap, phiMap ctor args	CPEuclideanMaps { f_idx, theta_e, phi_f } struct	✓
95–97	getDimension = hds.numFaces() (incl. face 0)	cp_euclidean_dimension(mesh, m) counts only f_idx ≥ 0	refined
103–123	getNonZeroPattern returns face-to-face adjacency	implicit in cp_euclidean_hessian triplet construction	✓
135–165	evaluateHessian — interior edge h_jk = sin θ / (cosh Δρ − cos θ)	cp_euclidean_hessian() lines 224-247	✓
146–151	"skip if k == 0 / j == 0" gauge fix	if (j >= 0) / if (k >= 0) guard in triplet emission	equivalent
178–193	per-face term +φ_f · ρ_f to energy and gradient	cp_euclidean_energy/gradient loop over faces	✓
196–232	per-(directed)-edge term: boundary or interior branch	cp_euclidean_energy/gradient halfedge loop with mesh.is_border(opposite) check	✓
224–227	E += ½ p Δρ + Λ(θ*+p) − θ* ρ_left	identical formula via clausen2()	✓
230	grad[left] -= p + θ*	identical	✓
243–247	p(θ*, Δρ) = 2 atan(tan(θ*/2) tanh(Δρ/2))	cp_detail::p_function	✓

1.1 Gauge convention divergence

The Java code uses "face index 0 is pinned" as an implicit convention: hard-coded if (face.getIndex() == 0) short-circuits scatter throughout. The C++ port keeps the gauge user-selectable via assign_cp_euclidean_face_dof_indices(mesh, m, pinned_face). The default convenience overload picks the first face in iteration order, matching the Java behaviour for any mesh whose face 0 is the first iterated face.

1.2 Test parity

Java test	C++ test	Reproduces what
init() lines 53-86 dodecahedron-minus-face setup	make_open_tetrahedron() (3 faces, 3 bdy edges)	open-mesh boundary code path
setXGradient(rho) (base class FunctionalTest)	FDGradientCheck_ClosedTetrahedron_RandomRho	FD-vs-analytic gradient
setXHessian(rho) (base class FunctionalTest)	FDHessianCheck_ClosedTetrahedron_RandomRho	FD-vs-analytic Hessian
(none — implicit from FunctionalTest convexity check)	HessianIsPSD	BPS-2010 §6 convexity
(none)	TangentialLimitGradientEqualsPhi	θ=0 analytic check
(none)	PFunctionKnownValues, NaturalPhiMakesZeroTheEquilibrium	helper-function unit tests

All Java tests reproduced; C++ adds 4 mathematical-property tests (PSD, tangential limit, p-function unit, natural-φ equilibrium).

1.3 Numerical equivalence — seed-1 random ρ, FD gradient

Java FunctionalTest uses Random(1) + uniform [-0.5, +0.5]. C++ uses std::mt19937(1) + uniform_real_distribution<double>(-0.5, 0.5). The two RNG streams differ; numerical equivalence is therefore tested via the property "FD matches analytic to 1e-6 absolute" rather than via identical numerical traces. Both implementations satisfy this property.

---

2. Inversive Distance (9a.2) — line-by-line literature mapping

No Java original exists. The implementation derives directly from three published papers; each formula in the header is annotated with its source line in the paper.

2.1 Edge-length formula (Luo 2004 §3 / Glickenstein 2011 eq. 2.1)

ℓ_ij(u)²  =  exp(2 u_i) + exp(2 u_j) + 2 I_ij exp(u_i + u_j)
          =  r_i² + r_j² + 2 I_ij r_i r_j

Implementation: id_detail::edge_length_squared (lines 130-138 of inversive_distance_functional.hpp).

Tested at three special-case limits in test_inversive_distance_functional.cpp §1:

Inversive distance	Geometric meaning	Closed-form ℓ	Test
I = +1	tangential circles	r_i + r_j	EdgeLengthFormula_TangentialLimit
I = 0	orthogonal circles	√(r_i² + r_j²)	EdgeLengthFormula_OrthogonalLimit
I = −1	inside-tangent / inverted	`	r_i − r_j
I < −1	impossible packing (no real ℓ)	nan/signal	EdgeLengthFormula_DegenerateReturnsMinusOne

2.2 Bowers-Stephenson identity (Bowers-Stephenson 2004)

I_ij  =  ( ℓ_ij² − r_i² − r_j² )  /  ( 2 r_i r_j )

Implementation: compute_inversive_distance_init_from_mesh lines 152-172. The round-trip property — given (ℓ_3d, r0_i, r0_j) recover ℓ via the Luo formula — is tested in BowersStephensonRoundTrip.

2.3 Gradient (Luo 2004 Lemma 3.1)

∂E/∂u_v  =  Θ_v  −  Σ_{T ∋ v}  α_v(T)

Implementation: inversive_distance_gradient lines 191-241. This is structurally identical to the Euclidean functional (same halfedge convention, same euclidean_angles() law-of-cosines reused). The only difference is the edge-length input: Luo's ℓ²(I) instead of Springborn's exp(λ° + u_i + u_j).

2.4 Energy — Luo's closed 1-form

Luo (2004) proves the gradient 1-form ω = Σ_v (Θ_v − Σ α_v) du_v is closed on the open domain where every triangle is valid. Hence the energy exists as a path integral. No general closed-form expression is known (Glickenstein 2011 surveys partial results). The implementation uses the same 10-point Gauss-Legendre quadrature as euclidean_functional.hpp (lines 243-274).

This shared-quadrature choice keeps both functionals identical in the energy-evaluation cost and the FD-gradient validation pattern.

2.5 Hessian — finite difference for now

Glickenstein 2011 eq. (4.6) gives an analytic Hessian for the inversive- distance variational principle, but is more involved than the Springborn cot-Laplacian. For MVP we rely on FD; the analytic form is a future optimisation (joining the Phase 9b roadmap for HyperIdeal as a sibling optimisation task).

---

3. Cross-validation between 9a.1 and 9a.2

The two functionals describe geometrically distinct circle packings (face-dual vs vertex-based) but agree in the angle-defect structure they expose to the Newton solver: both report Θ − Σ α_actual at every free DOF. This is tested directly:

TEST(InversiveDistanceFunctional, AngleDefectAtU0_AgreesWithEuclideanAtU0)
{
    auto mesh = make_quad_strip();

    auto G_id = inversive_distance_gradient(mesh, x_id, m_id);
    auto G_eu = euclidean_gradient            (mesh, x_eu, m_eu);

    for (i in DOFs) EXPECT_NEAR(G_id[i], G_eu[i], 1e-10);
}

Why this works: at u = 0, both initialisation procedures reconstruct the input 3-D edge length exactly (compute_*_lambda0_from_mesh for Euclidean, compute_inversive_distance_init_from_mesh for inversive distance). Therefore the actual angle sums per vertex are identical, and both gradients reduce to the same Θ_default − Σ α(ℓ_3d).

This is the operational consequence of Glickenstein 2011 §5: "different parametrisations of the same initial discrete metric produce the same Newton-time-zero gradient".

---

4. Comparison with euclidean_functional.hpp

Aspect	Euclidean	CP-Euclidean (9a.1)	Inversive Distance (9a.2)
DOF basis	per vertex (u_v = log scale)	per face (ρ_f = log radius)	per vertex (u_v = log radius)
Per-edge constant	λ°_ij = 2 log ℓ°_ij	θ_e (intersection angle)	I_ij (inversive distance)
Edge-length formula	ℓ = exp((λ° + u_i + u_j)/2)	(no per-edge ℓ — face-circle radii)	ℓ² = r_i² + r_j² + 2 I r_i r_j
Angles	half-tangent law of cosines	(face-internal angles via p(θ*, Δρ))	half-tangent law of cosines (reuses euclidean_angles)
Gradient	Θ_v − Σ α_v(f)	φ_f − Σ_{h:face(h)=f} (p+θ*)	Θ_v − Σ α_v(f)
Energy form	path integral (10-pt GL)	closed form via ½ p Δρ + Λ(θ*+p) − θ* ρ	path integral (10-pt GL)
Hessian	analytic (cotangent Laplacian)	analytic (sin θ / (cosh Δρ − cos θ))	finite difference (analytic deferred — Glickenstein 2011 eq. 4.6)
Convexity	strictly convex after gauge fix	strictly convex after gauge fix	locally convex on triangle-inequality domain (Luo 2004 Thm 1.2)
Gauge fix	pin one vertex (v_idx = −1)	pin one face (f_idx = −1)	pin one vertex (v_idx = −1)

4.1 Structural reuse

inversive_distance_functional.hpp reuses euclidean_angles() from euclidean_geometry.hpp unchanged. The only difference from euclidean_functional.hpp is the four lines that map (u_i, u_j, I_ij) to ℓ_ij² before invoking euclidean_angles().

cp_euclidean_functional.hpp shares no algorithmic code with the Euclidean functional — it operates on a different DOF basis and a different energy form. It does share the Clausen function via clausen2() from clausen.hpp.

---

5. Test counts and acceptance criteria

Suite	Tests	Status
cgal.CPEuclideanFunctional.*	10	✅ all pass
cgal.InversiveDistanceFunctional.*	11	✅ all pass
Full CGAL regression	205	✅ all pass
Non-CGAL fast suite	36	✅ all pass

Acceptance criteria — met

• FD-vs-analytic gradient check passes on all test meshes (both functionals)
• FD-vs-analytic Hessian check passes on closed and open tetrahedron (9a.1)
• Hessian PSD-property test passes (9a.1)
• Bowers-Stephenson round-trip identity verified (9a.2)
• Three special-case limits of Luo's edge formula verified at machine ε (9a.2)
• Cross-validation against euclidean_functional at u = 0 (9a.2)
• All formulas cite their source paper or Java line number

---

6. References

• Bobenko, A. I., Pinkall, U. & Springborn, B. (2010). Discrete conformal maps and ideal hyperbolic polyhedra. Geometry & Topology 14, 379–426.
• Bowers, P. L. & Stephenson, K. (2004). Uniformizing dessins and Belyĭ maps via circle packing. Memoirs of the AMS 170(805).
• Glickenstein, D. (2011). Discrete conformal variations and scalar curvature on piecewise flat manifolds. J. Differential Geometry 87(2), 201–238.
• Luo, F. (2004). Combinatorial Yamabe Flow on Surfaces. Communications in Contemporary Mathematics 6(5), 765–780.
• Java original (CPEuclidean only): /Users/tarikmoussa/Desktop/conformallab/src/de/varylab/discreteconformal/functional/CPEuclideanFunctional.java
• Java test: /Users/tarikmoussa/Desktop/conformallab/src-test/de/varylab/discreteconformal/functional/CPEuclideanFunctionalTest.java