diff --git a/doc/architecture/locked-vs-flexible.md b/doc/architecture/locked-vs-flexible.md new file mode 100644 index 0000000..c30caaf --- /dev/null +++ b/doc/architecture/locked-vs-flexible.md @@ -0,0 +1,281 @@ +# Locked-in vs flexible architecture decisions + +> **Purpose.** External collaborators (especially mathematicians +> evaluating whether to extend the library) need to know which +> design choices are **load-bearing** (changing them is expensive +> across the whole codebase) and which are **opportunistic** (made +> on first principles and easy to revisit). +> +> **Why this matters now.** A v0.9.0 + Springborn-Bobenko-alumnus +> external review (May 2026) is the right moment to surface these +> decisions before they ossify further. If any of the *locked* +> decisions need revisiting, this is the cheapest moment in the +> project's life to do so. + +Three tiers: + +``` +🔴 LOAD-BEARING changing requires repo-wide refactoring +🟡 SEMI-FIXED changing affects multiple subsystems; possible but not casual +🟢 OPPORTUNISTIC changing is one-PR work +``` + +--- + +## 1. Mesh data structure → `CGAL::Surface_mesh

` + +| | | +|--|--| +| Status | 🔴 load-bearing | +| Locked since | Phase 3 (2025) | +| Alternative considered | OpenMesh / pmp-library / custom halfedge / Java `CoHDS` literal port | +| Why locked | Every header `code/include/*.hpp` uses `ConformalMesh = CGAL::Surface_mesh` and CGAL property maps explicitly. Phase 8a Traits provide a *concept* abstraction but the Default model is `Surface_mesh`-only. | +| Cost to change | ~3 weeks: rewrite traits, every functional, every test mesh factory. Touched in ~30 headers + ~25 test files. | +| Mitigation | Traits-based design (Phase 8a) means user-supplied mesh types CAN be added as Default-traits specialisations without touching algorithms — Phase 8a.2 plan documents this. But the Default model is Surface_mesh. | +| When to revisit | If a major user wants Polyhedron_3 or OpenMesh as default. Currently no concrete request → keep. | + +**Recommended posture for an external contributor:** use the existing +Surface_mesh-based API for any new functional. Adding generic-FaceGraph +support is a single architectural step that doesn't need to be repeated +per functional — wait for one user to need it. + +--- + +## 2. Floating-point kernel → `CGAL::Simple_cartesian` + +| | | +|--|--| +| Status | 🟡 semi-fixed | +| Locked since | Phase 3 (deliberate decision: conformal geometry doesn't need exact predicates) | +| Cost to change | Per-functional template parameterisation. The Phase 8 MVP wrapper already deduces the kernel from the mesh point type, so user code can specify any CGAL kernel — but the **legacy** `code/include/*.hpp` headers are hardcoded. | +| When to revisit | If a user reports floating-point catastrophic cancellation in the half-tangent angle formula on extreme meshes. Hasn't happened in 250 tests over 9 phases. | + +**Recommended posture:** stay with `Simple_cartesian`. If a +specific algorithm needs `Exact_predicates_inexact_constructions_kernel` +for robustness, parameterise that one algorithm — don't refactor the +whole codebase. + +--- + +## 3. Header-only, no compiled library + +| | | +|--|--| +| Status | 🔴 load-bearing | +| Locked since | Phase 1 | +| Cost to change | Significant: would need to introduce `.cpp` files, link order, ABI compatibility decisions. But everything in `code/include/*.hpp` is `inline` or template, so a header-only-to-compiled migration is mechanical. | +| Why locked | CGAL package convention is also header-only. Forking would break that goal. | +| When to revisit | If compilation times become prohibitive (currently < 30 s clean build with all 5 functionals). Or if a future cyclic dependency between functionals forces it. Neither is on the horizon. | + +**Recommended posture:** stay header-only. This is also the de-facto +norm for CGAL packages of comparable scope. + +--- + +## 4. Five DCE models on the same mesh + +| | | +|--|--| +| Status | 🟡 semi-fixed | +| Locked since | Phase 9a (2026-05) — when CP-Euclidean introduced face-DOFs | +| Why semi-fixed | Each model has its own `*Maps` bundle with its own property-map name prefix (`ev:`, `sv:`, `v:`, `cf:`/`ce:`, `iv:`/`ie:`) so all five can coexist on the same `CGAL::Surface_mesh`. Adding a sixth model means picking a new prefix + writing a new `*Maps` struct + Default trait. | +| Cost to add a sixth model | ~1 week (CP-Euclidean took ~3 days, Inversive-Distance ~3 days, Hessian + Newton + CGAL entry add another ~3 days). | +| When to revisit | If a unified base-Maps abstraction would actually win something (currently it would not — the property-map sets differ in *kind*, not just in name). | + +**Recommended posture for adding a new functional:** +1. Pick a 2-letter prefix not in `{ev, sv, v, cf, ce, iv, ie}`. +2. Define your `*Maps` struct with that prefix. +3. Define your `Default_*_traits` class in a new + header `code/include/CGAL/Discrete_*.h`. +4. Wire it into `newton_solver.hpp` (template-copy from + `newton_inversive_distance` is the closest pattern for vertex-DOFs; + `newton_cp_euclidean` for face-DOFs). +5. Add a CGAL entry in the new header. +6. Add tests following [`add-inversive-distance.md`](../tutorials/add-inversive-distance.md). + +This recipe has been validated three times now (CP-Euclidean, Inversive +Distance, and the four Phase-8b-Lite wrappers). + +--- + +## 5. Newton solver with line search + SparseQR fallback + +| | | +|--|--| +| Status | 🟡 semi-fixed | +| Locked since | Phase 4 | +| Cost to change | Each `newton_*` function is ~50-80 lines; replacing the solver across all five is ~2 days. | +| When to revisit | If a future functional needs trust-region or BFGS. None of the five currently does — Newton converges quadratically near the optimum and the line search handles bad initial points. | + +**Recommended posture:** Newton with line search is enough for any +strictly-convex variational problem. For non-convex variants, consider +adding a `newton_with_trust_region()` helper alongside, not replacing. + +--- + +## 6. Eigen as the linear-algebra back-end + +| | | +|--|--| +| Status | 🔴 load-bearing | +| Locked since | Phase 4 | +| Alternative considered | PETSc/Tao (Java original) / Boost.uBLAS / Blaze | +| Why locked | Eigen is header-only (no external dependency at build time), bundled as a CGAL dependency, fast, and offers `SimplicialLDLT + SparseQR` which the gauge-singular-mesh case needs. | +| Cost to change | Significant — every Hessian header (`*_hessian.hpp`) and every Newton solver uses `Eigen::SparseMatrix` and `Eigen::VectorXd` directly. ~2 weeks repo-wide. | +| When to revisit | If a sparse-solver feature (e.g. parallel Cholesky) is needed that Eigen doesn't offer. | + +**Recommended posture:** stay with Eigen. + +--- + +## 7. CGAL public-API surface layout + +| | | +|--|--| +| Status | 🟡 semi-fixed (Phase 8a-MVP design decision, 2026-05-19) | +| Locked since | PR #6 (v0.9.0) | +| Strategy chosen | "Strategy C" — functional-specific Default traits, one entry function per functional, no fat unified trait. | +| Cost to change to unified trait | ~1 week — refactor `Default_*_traits<>` into a single `Default_conformal_map_traits<>` with all property-map fields. Existing 8 tests would need updating. | +| When to revisit | When the first cross-functional algorithm (e.g. a hybrid functional that uses both face and vertex DOFs) lands. Speculation today. | + +**Recommended posture:** stay with Strategy C. CGAL's own +`Polygon_mesh_processing` package follows the same convention — one +default trait per algorithm family. + +--- + +## 8. Named-parameter mechanism + +| | | +|--|--| +| Status | 🟢 opportunistic | +| Locked since | Phase 8 MVP (2026-05-19) | +| Current state | Six tags: `vertex_curvature_map`, `fixed_vertex_map`, `gradient_tolerance`, `max_iterations`, `output_uv_map`, `normalise_layout`. No chaining yet. | +| Cost to extend with chaining | ~2 days. Generate CGAL-style member-function chainers via macros (CGAL has a `CGAL_add_named_parameter` macro for this). | +| When to revisit | At any time; this is the lowest-risk change in the codebase. Tutorials [`add-output-uv-map.md`](../tutorials/add-output-uv-map.md) §4 explains the mechanism. | + +**Recommended posture:** add chaining when the first user complains. Or +when adding the 8th-10th named parameter (still 6 today — manageable). + +--- + +## 9. Property-map name conventions + +| | | +|--|--| +| Status | 🟢 opportunistic | +| Locked since | Phase 3 + 9a (prefix `ev:`/`sv:`/`v:`/`cf:`/`ce:`/`iv:`/`ie:` set when each functional was introduced) | +| Cost to change | One sed-replace + recompile. No user-visible effect because the names are an *internal* convention; the CGAL public API never exposes them. | +| When to revisit | If a future functional reuses an existing letter prefix. Already discussed in [`locked-vs-flexible.md`](#4-five-dce-models-on-the-same-mesh). | + +--- + +## 10. Tests: GTest, not CGAL's own test format + +| | | +|--|--| +| Status | 🟡 semi-fixed | +| Locked since | Phase 1 | +| Cost to change | ~1 week — rewrite test harnesses to CGAL's `test/Conformal_map/` convention. This is Phase 8d (planned for CGAL submission). | +| Why GTest now | Faster development cycle, IDE-friendly (Xcode / VSCode / CLion all have native GTest support). No CGAL submission is in progress yet. | +| When to revisit | When committing to CGAL submission (Phase 8c-d, decided to be a future commitment, see [`release-policy.md`](../release-policy.md)). | + +**Recommended posture:** keep GTest as primary. When/if CGAL submission +happens, add a `test/Conformal_map/` shim that calls into the GTest +suite — both formats can coexist. + +--- + +## 11. License: MIT + +| | | +|--|--| +| Status | 🔴 load-bearing | +| Locked since | Project inception | +| Cost to change | High organisational cost (requires consent of all contributors); ~no code cost. | +| Why locked | MIT was chosen for academic friendliness (citing, modifying, embedding). CGAL upstream requires LGPL for submitted packages. | +| Trade-off | Submitting to CGAL upstream is not possible without re-licensing. The codebase architecture is "CGAL-style" but the project would publish independently. | +| When to revisit | If/when a concrete CGAL upstream submission is decided. See [`release-policy.md`](../release-policy.md) for the formal policy. | + +**Recommended posture:** stay with MIT. Build the "CGAL-style package +for external distribution" as the primary deliverable. Re-license only +when the CGAL editorial board commits to accepting the submission. + +--- + +## 12. Documentation pattern: Markdown + Doxygen + +| | | +|--|--| +| Status | 🟢 opportunistic | +| Locked since | Phase 7.5 (2026-05) | +| Current state | `code/include/*.hpp` carry Doxygen-style `///` comments (87% coverage); `doc/*.md` for prose; `Doxyfile` generates HTML in `doc/doxygen/`. | +| Cost to add more | Per-file basis; ~1 hour per header for full Doxygen. | +| When to revisit | When chasing CGAL-submission readiness (need `PackageDescription.txt` + `User_manual.md`). | + +**Recommended posture:** keep adding `///` comments incrementally with +each new public function. + +--- + +## Summary table + +| Decision | Tier | Cost to change | +|-----------------------------------------|----------------|----------------------| +| 1. CGAL::Surface_mesh as default mesh | 🔴 load-bearing | ~3 weeks | +| 2. Simple_cartesian kernel | 🟡 semi-fixed | per-functional | +| 3. Header-only architecture | 🔴 load-bearing | medium (mechanical) | +| 4. Five DCE models, separate Maps | 🟡 semi-fixed | ~1 week per new model | +| 5. Newton + line search + SparseQR | 🟡 semi-fixed | ~2 days | +| 6. Eigen back-end | 🔴 load-bearing | ~2 weeks | +| 7. Strategy C (per-functional traits) | 🟡 semi-fixed | ~1 week | +| 8. Named-parameter mechanism | 🟢 opportunistic | ~2 days (chaining) | +| 9. Property-map name conventions | 🟢 opportunistic | ~1 hour | +| 10. GTest, not CGAL test format | 🟡 semi-fixed | ~1 week | +| 11. MIT license | 🔴 load-bearing | organisational | +| 12. Markdown + Doxygen | 🟢 opportunistic | per-file | + +**Key insight:** the **load-bearing decisions are all good in 2026**. +Surface_mesh + Eigen + header-only + MIT are the right defaults for a +research-quality CGAL-style package. The **semi-fixed decisions are +all behind one concrete blocker** (single user request, CGAL submission +commitment, etc.). The **opportunistic decisions are cheap to revisit +any time**. + +The architecture is in a good place for the v0.9.0 → v0.10.0 transition. +No "expensive corner" has been painted into; every locked decision +matches the project's three-goal hierarchy in +[`research-track.md`](../roadmap/research-track.md). + +--- + +## Open questions for the external reviewer + +Items where the project would benefit from a second opinion: + +1. **Phase 9c (4g-polygon) algorithm choice.** Two routes: + * Port the Java `FundamentalPolygonUtility` + `CanonicalFormUtility` + literally (~2 weeks). + * Or: re-derive from Springborn 2020 §5 using the existing + `cut_graph.hpp` + holonomy infrastructure (~3 weeks, cleaner + architecture). + Which is preferred? See [`phases.md`](../roadmap/phases.md) §Phase 9c. + +2. **Phase 10a (forms) priorities.** Three sub-items + (`DiscreteHarmonicFormUtility`, `DiscreteHolomorphicFormUtility`, + `CanonicalBasisUtility`) interlock. Which to start with? + +3. **Analytic Hessian payoff.** The Schläfli-based analytic HyperIdeal + Hessian (Phase 9b-analytic — derivation already written: + [`hyperideal-hessian-derivation.md`](../math/hyperideal-hessian-derivation.md)) + would add another ~6× over block-FD. Is that worth ~2 weeks of + implementation effort for a working-mesh size on which? + +4. **CGAL upstream vs independent distribution.** Does the reviewer + know a CGAL editor / has personal opinion on the LGPL-vs-MIT + trade-off? + +5. **geometry-central cross-validation (GC-1).** Two libraries solve + the same DCE problem from different algorithmic directions + (Newton-on-mesh vs Ptolemaic-flips-on-intrinsic-triangulation). An + independent comparison would be a nice paper. Interested? diff --git a/doc/math/hyperideal-hessian-derivation.md b/doc/math/hyperideal-hessian-derivation.md new file mode 100644 index 0000000..029a950 --- /dev/null +++ b/doc/math/hyperideal-hessian-derivation.md @@ -0,0 +1,805 @@ +# Analytic Hessian of the Hyper-Ideal Discrete Conformal Energy + +**Status:** research note, mathematical preparation for **Phase 9b-analytic** +(see `doc/roadmap/research-track.md`). The current shipped Hessian +(`hyper_ideal_hessian_block_fd`, v0.9.0) is block finite-difference; this note +derives the closed-form replacement via the Schläfli identity and the chain +rule through the building blocks `zeta`, `zeta13`, `zeta14`, `zeta15`, `lij`, +`alpha_ij`, `sigma_i`, `sigma_ij` declared in `code/include/hyper_ideal_geometry.hpp` +and assembled by `face_angles_from_local_dofs(...)` in +`code/include/hyper_ideal_functional.hpp`. + +**Target reader.** A mathematician fluent with Springborn (2020) and the +Bobenko-school discrete-conformal apparatus; the goal is verifiability of +each derivative line against the source papers. + +**Conventions.** +- $b_i \in \mathbb{R}$ — log scale factor at hyper-ideal vertex $i$ (the DOF). + An "ideal" vertex has no $b_i$ DOF; we encode that by `vi_var = false`. +- $a_{ij} \in \mathbb{R}$ — intersection-angle DOF on edge $ij$ (in the + Springborn $a$-parametrisation, edges may be conventional, ideal–ideal, + or hyper-ideal–hyper-ideal). +- $\ell_{ij}$ — effective hyperbolic length of edge $ij$ in the auxiliary + truncated tetrahedron. +- $\beta_i^{(f)}$ — interior angle of the auxiliary hyperbolic triangle on + face $f$ at vertex $i$. +- $\alpha_{ij}^{(f)}$ — dihedral angle of the truncated tetrahedron at edge + $ij$, contributed by face $f$. +- $\Theta_v$, $\theta_e$ — prescribed cone / intersection-angle targets. +- $V(f)$ — hyperbolic volume of the truncated tetrahedron associated with + face $f$. + +--- + +## 1. The Schläfli identity and the gradient + +### 1.1 Energy + +The Springborn (2020, §4) hyper-ideal energy on a triangulated surface $M$ +with DOF vector $x = (b, a)$ is +$$ +E(b, a) \;=\; \sum_{f \in F} U(f) \;-\; \sum_{e \in E} \theta_e\, a_e \;-\; + \sum_{v \in V} \Theta_v\, b_v , +$$ +where the per-face contribution is +$$ +U(f) \;=\; \sum_{e \in f} a_e\, \alpha_e^{(f)} \;+\; + \sum_{i \in f} b_i\, \beta_i^{(f)} \;+\; 2\, V(f) . +$$ +(Each face contributes its three edges and three vertices; ideal vertices +contribute $b_i = 0$ trivially.) This matches the `face_energy(...)` routine +in `hyper_ideal_functional.hpp` line by line. + +### 1.2 First-order Schläfli (gradient) + +The classical Schläfli differential identity (Schläfli 1858/60; Milnor 1982; +Vinberg 1993, Ch. 7) for the volume of any compact hyperbolic polyhedron +$P \subset \mathbb{H}^3$ reads +$$ +\boxed{\;\; -2\, dV \;=\; \sum_{e \subset P} \ell_e \, d\alpha_e \;\;} +$$ +in the closed (finite-vertex) case, where $\ell_e, \alpha_e$ are edge length +and dihedral angle. For the truncated / hyper-ideal extension we follow +Bobenko–Springborn–Schief and write the identity in the *mixed* form that +respects the truncation: +$$ +\boxed{\;\; 2\, dV \;=\; \sum_e a_e\, d\alpha_e \;+\; \sum_v b_v\, d\beta_v + \;-\; \sum_e \alpha_e\, da_e \;-\; \sum_v \beta_v\, db_v . \;\;} +\tag{S1} +$$ +Combined with the trivial $d(a_e \alpha_e) = a_e\, d\alpha_e + \alpha_e\, da_e$ +identity, (S1) is equivalent to +$$ +d\!\left(\sum_e a_e \alpha_e + \sum_v b_v \beta_v + 2V\right) + \;=\; \sum_e \alpha_e\, da_e \;+\; \sum_v \beta_v\, db_v , +$$ +i.e. $dU(f) = \sum \alpha_e\, da_e + \sum \beta_v\, db_v$ on each face. +Summing over faces and subtracting the linear $\theta, \Theta$ terms, +$$ +\frac{\partial E}{\partial b_v} \;=\; \Big(\!\!\sum_{f \ni v} \beta_v^{(f)}\Big) - \Theta_v, +\qquad +\frac{\partial E}{\partial a_e} \;=\; \Big(\!\!\sum_{f \supset e} \alpha_e^{(f)}\Big) - \theta_e , +$$ +which is precisely the gradient implemented in +`evaluate_hyper_ideal(...)` (pass 3). + +### 1.3 Second-order Schläfli (Hessian) + +Differentiating (S1) once more (Cho–Kim 1999, Lemma 2.1; Rivin 1994) +yields the second Schläfli identity +$$ +\boxed{\;\; 0 \;=\; \sum_e da_e \wedge d\alpha_e \;+\; \sum_v db_v \wedge d\beta_v . \;\;} +\tag{S2} +$$ +Since the wedge product is antisymmetric in the differential factors, (S2) +forces the bilinear form +$$ +H(f) := \begin{pmatrix} + \partial \beta_i / \partial b_j & \partial \beta_i / \partial a_e \\ + \partial \alpha_e / \partial b_j & \partial \alpha_e / \partial a_{e'} +\end{pmatrix} +$$ +to be **symmetric** on each face (and so on the whole mesh). This is the +deep reason that one may compute the Hessian +$$ +H \;=\; \frac{\partial^2 E}{\partial x^2} \;=\; +\frac{\partial (\beta - \Theta,\; \alpha - \theta)}{\partial (b, a)} +$$ +by accumulating only the *upper triangle* face by face, and that the +mesh-level matrix inherits PSD-ness from the per-face $6 \times 6$ blocks +(Springborn 2020 §4.3). + +### 1.4 Strategy + +The full chain through `face_angles_from_local_dofs(...)` is +$$ +(b_i, a_e) \;\xrightarrow{\;\text{lij}\;}\; \ell_e + \;\xrightarrow{\;\zeta\;}\; \beta_i + \;\xrightarrow{\;\zeta,\sigma\;}\; \alpha_e . +$$ +We derive each arrow in Sections 2–3 and assemble the $6 \times 6$ face +Jacobian in Sections 4–5. Symmetry is then a consequence of (S2) and a +useful numerical check. + +--- + +## 2. Derivative of $\zeta$ (interior angle from edge lengths) + +### 2.1 Definition + +For an auxiliary hyperbolic triangle with side lengths $x, y, z$ (opposite +to vertices $X, Y, Z$), +$$ +\zeta(x, y, z) \;=\; \arccos\!\left(\frac{\cosh x \cosh y - \cosh z}{\sinh x \sinh y}\right) , +$$ +returning the interior angle at the vertex *between* the sides of length +$x$ and $y$ (so $z$ is opposite). This is the hyperbolic law of cosines +solved for the included angle, and it is exactly `zeta(...)` in +`hyper_ideal_geometry.hpp`. + +### 2.2 Partial derivatives + +Let +$$ +N := \cosh x \cosh y - \cosh z, +\qquad D := \sinh x \sinh y, +\qquad N / D = \cos \beta . +$$ +Then $\sin \beta \cdot d\beta = -d(N/D) = (D\, dN - N\, dD)/D^2 \cdot (-1)$. +Computing the partials, +$$ +\partial_x N = \sinh x \cosh y,\qquad +\partial_y N = \cosh x \sinh y,\qquad +\partial_z N = -\sinh z , +$$ +$$ +\partial_x D = \cosh x \sinh y,\qquad +\partial_y D = \sinh x \cosh y,\qquad +\partial_z D = 0 . +$$ +Hence +$$ +\partial_x (\cos \beta) \;=\; \frac{(\sinh x \cosh y)\,\sinh x \sinh y - (\cosh x \cosh y - \cosh z)\, \cosh x \sinh y}{\sinh^2 x \sinh^2 y} +$$ +$$ += \;\frac{\sinh y \big[ \sinh^2 x \cosh y - \cosh x (\cosh x \cosh y - \cosh z) \big]}{\sinh^2 x \sinh^2 y} +\;=\;\frac{\cosh x \cosh z - \cosh y}{\sinh^2 x \sinh y} +$$ +(using $\sinh^2 x - \cosh^2 x = -1$, which gives the cancellation +$\sinh^2 x \cosh y - \cosh^2 x \cosh y = -\cosh y$). + +Combined with $-\sin \beta\, \partial_x \beta = \partial_x (\cos\beta)$ this +yields the **hyperbolic dual law of cosines** in derivative form +(Cho–Kim 1999, Lemma 3.1): +$$ +\boxed{\;\; \frac{\partial \beta}{\partial x} + \;=\; \frac{\cosh y - \cosh x \cosh z}{\sin \beta \cdot \sinh^2 x \sinh y} + \;=\; -\frac{\cos \beta_y}{\sinh x \sin \beta} \,\cdot\, \frac{1}{?} \;\;} +$$ +We prefer the *normalised* form (sine rule). Recall the hyperbolic sine +rule on a triangle $\beta, \beta_y, \beta_z$: +$\sinh x / \sin \beta_x' = \sinh y / \sin \beta_y' = \sinh z / \sin \beta_z'$, +where $\beta_x'$ is opposite to side $x$. In our convention $\beta = \beta_z'$ +(the angle opposite $z$ — wait: $\zeta(x,y,z)$ returns the angle *between* +sides $x$ and $y$, i.e. **opposite to side $z$**), so let +$\beta := \beta_z'$, $\beta_x' := $ angle opposite $x$, $\beta_y' := $ angle +opposite $y$. Then by the dual law of cosines applied to side $y$, +$\cosh y = \cosh x \cosh z - \sinh x \sinh z \cos \beta_y'$, whence +$\cosh y - \cosh x \cosh z = -\sinh x \sinh z \cos \beta_y'$. +Substituting, +$$ +\frac{\partial \beta}{\partial x} + \;=\; \frac{-\sinh x \sinh z \cos \beta_y'}{\sin\beta \cdot \sinh^2 x \sinh y} + \;=\; \frac{-\sinh z \cos \beta_y'}{\sin \beta \cdot \sinh x \sinh y} . +$$ +Using $\sin \beta / \sinh z = \sin \beta_y'/ \sinh y$ (sine rule) once more, +$\sin \beta \sinh y = \sin \beta_y' \sinh z$, and we obtain the **clean +Cho–Kim form**: +$$ +\boxed{\;\; \frac{\partial \beta}{\partial x} \;=\; -\frac{\cos \beta_y'}{\sinh x}, \quad + \frac{\partial \beta}{\partial y} \;=\; -\frac{\cos \beta_x'}{\sinh y}, \quad + \frac{\partial \beta}{\partial z} \;=\; +\frac{1}{\sinh z} \cdot \frac{\sinh z}{\sin\beta\,\sinh x \sinh y} \cdot \sinh z \;\;} +$$ +For the third partial, $\partial_z N = -\sinh z$, $\partial_z D = 0$, so +$$ +-\sin \beta \cdot \partial_z \beta = \frac{-\sinh z}{\sinh x \sinh y} +\quad \Longrightarrow \quad +\boxed{\;\; \frac{\partial \beta}{\partial z} + \;=\; \frac{\sinh z}{\sin \beta \cdot \sinh x \sinh y} + \;=\; \frac{1}{\sin \beta} \cdot \frac{\sinh z}{\sinh x \sinh y} . \;\;} +$$ +Equivalently, $\partial \beta / \partial z = \sinh z / (\sin \beta \sinh x \sinh y)$, +which by the sine rule is also $1 / (\sinh x \sin \beta_y')$, recovering a +form symmetric to the first two. + +**Summary (operational form used in code).** For +$\beta = \zeta(x, y, z)$ (opposite to $z$): +$$ +\begin{aligned} +\partial_x \beta &= \frac{1}{\sin \beta} \cdot \frac{\cosh y - \cosh x \cosh z}{\sinh^2 x \sinh y}, \\ +\partial_y \beta &= \frac{1}{\sin \beta} \cdot \frac{\cosh x - \cosh y \cosh z}{\sinh x \sinh^2 y}, \\ +\partial_z \beta &= \frac{1}{\sin \beta} \cdot \frac{\sinh z}{\sinh x \sinh y} . +\end{aligned} +\tag{Z} +$$ +These three lines are what the analytic kernel will evaluate (they are +finite as long as no edge degenerates and $\sin \beta \neq 0$, which is +guaranteed inside the triangle-inequality regime gated by +`face_angles_from_local_dofs`). + +--- + +## 3. Derivatives of $\ell_{ij}$ (edge-length building blocks) + +The length dispatcher `lij(b_i, b_j, a_{ij}, v_i, v_j)` in +`hyper_ideal_geometry.hpp` selects among $\zeta_{13}, \zeta_{14}, \zeta_{15}$ +based on which vertices are hyper-ideal. We treat each branch separately. + +### 3.1 Both vertices hyper-ideal: $\ell_{ij} = \zeta_{13}(b_i, b_j, a_{ij})$ + +$$ +\zeta_{13}(x, y, z) = \operatorname{arcosh}\!\left( \frac{\cosh x \cosh y + \cosh z}{\sinh x \sinh y} \right) . +$$ +Let $L = \zeta_{13}$, $\Phi := \cosh L = (\cosh x \cosh y + \cosh z)/(\sinh x \sinh y)$. +Then $\sinh L \cdot \partial L = \partial \Phi$. With +$$ +\partial_x \Phi = \frac{\sinh x \cosh y \cdot \sinh x \sinh y - (\cosh x \cosh y + \cosh z) \cosh x \sinh y}{\sinh^2 x \sinh^2 y} += \frac{-\cosh x \cosh z - \cosh y}{\sinh^2 x \sinh y}, +$$ +(using $\sinh^2 x - \cosh^2 x = -1$), and analogously for $y$, $z$: +$$ +\boxed{\;\; +\begin{aligned} +\partial_x \zeta_{13} &= -\frac{1}{\sinh L} \cdot \frac{\cosh y + \cosh x \cosh z}{\sinh^2 x \sinh y}, \\ +\partial_y \zeta_{13} &= -\frac{1}{\sinh L} \cdot \frac{\cosh x + \cosh y \cosh z}{\sinh x \sinh^2 y}, \\ +\partial_z \zeta_{13} &= +\frac{1}{\sinh L} \cdot \frac{\sinh z}{\sinh x \sinh y} . +\end{aligned} +\;\;} +\tag{Z13} +$$ +Compare (Z) and (Z13): they differ only in (i) the sign in front of +$\cosh z$ in the numerator of the first two partials, and (ii) the +prefactor ($1/\sin \beta$ vs $1/\sinh L$) — which mirrors the +$\arccos / \operatorname{arcosh}$ duality. This is the dual hyperbolic +law of cosines for a right-angled hexagon +(Buser 1992, *Geometry and Spectra of Compact Riemann Surfaces*, Thm 2.4.1). + +### 3.2 One ideal vertex: $\ell_{ij} = \zeta_{14}(a_{ij}, b_j)$ + +If vertex $i$ is ideal (i.e. only $v_j$ is hyper-ideal), +$$ +\zeta_{14}(x, y) = \operatorname{arcosh}\!\left( \frac{e^x + \cosh y}{\sinh y} \right) . +$$ +Let $L = \zeta_{14}$, $\Phi = \cosh L = (e^x + \cosh y)/\sinh y$. Then +$$ +\partial_x \Phi = \frac{e^x}{\sinh y}, +\qquad +\partial_y \Phi = \frac{\sinh y \cdot \sinh y - (e^x + \cosh y)\cosh y}{\sinh^2 y} + = \frac{-1 - e^x \cosh y}{\sinh^2 y} . +$$ +Hence +$$ +\boxed{\;\; +\partial_x \zeta_{14} \;=\; \frac{e^x}{\sinh L \cdot \sinh y}, +\qquad +\partial_y \zeta_{14} \;=\; -\frac{1 + e^x \cosh y}{\sinh L \cdot \sinh^2 y} . +\;\;} +\tag{Z14} +$$ +Note that `lij` invokes $\zeta_{14}$ with argument order *(edge, vertex)*, +so when $v_i$ is the ideal one the call is `zeta14(a_{ij}, b_j)` and we +must remember that $\partial / \partial a_{ij} = \partial_x \zeta_{14}$, +$\partial / \partial b_j = \partial_y \zeta_{14}$. The symmetric case +(when $v_j$ is the ideal one) flips the roles of $b_i / b_j$. + +### 3.3 Both vertices ideal: $\ell_{ij} = \zeta_{15}(a_{ij})$ + +$$ +\zeta_{15}(x) = 2\, \operatorname{arsinh}\!\big(e^{x/2}\big), \qquad +\partial_x \zeta_{15} = \frac{2 \cdot \tfrac{1}{2} e^{x/2}}{\sqrt{1 + e^x}} + = \frac{e^{x/2}}{\sqrt{1 + e^x}} . +$$ +A more numerically symmetric form uses +$\cosh(\zeta_{15}/2) = \sqrt{1 + e^x}$, $\sinh(\zeta_{15}/2) = e^{x/2}$, +so +$$ +\boxed{\;\; +\partial_x \zeta_{15} \;=\; \tanh\!\big(\zeta_{15}(x)/2\big) . \;\;} +\tag{Z15} +$$ +This is the cleanest form for cross-checks against the FD reference. + +### 3.4 Combined edge-length partials + +We package the per-edge $\partial \ell_{ij}$ as a *length differential* +$$ +d\ell_{ij} \;=\; L^b_{ij,i}\, db_i \;+\; L^b_{ij,j}\, db_j \;+\; L^a_{ij}\, da_{ij} +$$ +with case-by-case coefficients: + +| Case ($v_i$, $v_j$) | $L^b_{ij,i}$ | $L^b_{ij,j}$ | $L^a_{ij}$ | +|---|---|---|---| +| (hyp, hyp) | $\partial_x \zeta_{13}$ at $(b_i, b_j, a_{ij})$ | $\partial_y \zeta_{13}$ | $\partial_z \zeta_{13}$ | +| (ideal, hyp) | — | $\partial_y \zeta_{14}$ at $(a_{ij}, b_j)$ | $\partial_x \zeta_{14}$ | +| (hyp, ideal) | $\partial_y \zeta_{14}$ at $(a_{ij}, b_i)$ | — | $\partial_x \zeta_{14}$ | +| (ideal, ideal) | — | — | $\partial_x \zeta_{15}$ at $a_{ij}$ | + +A "—" entry means the corresponding DOF does not exist; the partial is +identically zero on the constraint surface. + +--- + +## 4. Chain-rule assembly (interior angles and dihedrals) + +### 4.1 Interior angles $\beta_i$ + +On face $f = (1, 2, 3)$ with edges $\ell_{12}, \ell_{23}, \ell_{31}$, +the code computes +$$ +\beta_1 = \zeta(\ell_{12}, \ell_{31}, \ell_{23}),\quad +\beta_2 = \zeta(\ell_{23}, \ell_{12}, \ell_{31}),\quad +\beta_3 = \zeta(\ell_{31}, \ell_{23}, \ell_{12}) . +$$ +With (Z) we get, for any DOF $\xi \in \{b_1, b_2, b_3, a_{12}, a_{23}, a_{31}\}$, +$$ +\frac{\partial \beta_1}{\partial \xi} + \;=\; \zeta_x(\ell_{12}, \ell_{31}, \ell_{23})\, \partial_\xi \ell_{12} + \;+\; \zeta_y(\ell_{12}, \ell_{31}, \ell_{23})\, \partial_\xi \ell_{31} + \;+\; \zeta_z(\ell_{12}, \ell_{31}, \ell_{23})\, \partial_\xi \ell_{23} , +\tag{B1} +$$ +and cyclically for $\beta_2, \beta_3$. Each $\partial_\xi \ell_{e}$ is read +from the table in §3.4: of the six DOFs only the three that touch edge $e$ +contribute (i.e. $b_{e^-}, b_{e^+}, a_e$). + +### 4.2 Dihedral angles $\alpha_{ij}$ — general structure + +`alpha_ij(...)` has four branches depending on $(v_i, v_j, v_k)$. +Define the "$\sigma$-triangle" attached to vertex $i$: +$$ +s_i \;=\; \sigma_i(a_{ij}, a_{ki}, a_{jk}; v_j, v_k), \qquad +s_{ij} \;=\; \sigma_{ij}(a_{ij}, b_i, b_j; v_j), \qquad +s_{ik} \;=\; \sigma_{ij}(a_{ki}, b_i, b_k; v_k) . +$$ +Then in the *hyper-ideal $v_i$* branch +$$ +\alpha_{ij} \;=\; \zeta(s_i, s_{ij}, s_{ik}) +\tag{A.v_i} +$$ +— note that the dihedral $\alpha_{ij}$ is computed as an interior angle in +the half-triangle at vertex $i$, with $s_i$ playing the role of side $x$, +$s_{ij}$ of $y$ and $s_{ik}$ of $z$. + +The branches $v_j$ (hyper-ideal but $v_i$ ideal), $v_k$ (one level of +recursion), and *all ideal* (closed form) follow the same pattern. + +### 4.3 Derivatives of $\sigma_i$ and $\sigma_{ij}$ + +Both are themselves dispatchers over $\zeta_{13}, \zeta_{14}, \zeta_{15}$. + +**$\sigma_{ij}(a_{ij}, b_i, b_j; v_j)$.** +- If $v_j$ hyper-ideal: $\sigma_{ij} = \zeta_{13}(a_{ij}, b_i, b_j)$ + with partials $(\partial_x \zeta_{13}, \partial_y \zeta_{13}, \partial_z \zeta_{13})$ + evaluated at $(a_{ij}, b_i, b_j)$. **Note the argument order:** the + first slot of $\zeta_{13}$ is $a_{ij}$, not $b_i$. +- If $v_j$ ideal: $\sigma_{ij} = \zeta_{14}(-a_{ij}, b_i)$. Then + $\partial_{a_{ij}} \sigma_{ij} = -\partial_x \zeta_{14}(-a_{ij}, b_i)$, + $\partial_{b_i} \sigma_{ij} = \partial_y \zeta_{14}(-a_{ij}, b_i)$, + $\partial_{b_j} \sigma_{ij} = 0$. + +**$\sigma_i(a_{ij}, a_{ki}, a_{jk}; v_j, v_k)$.** +- $(v_j, v_k) = $ (hyp, hyp): $\sigma_i = \zeta_{13}(a_{ij}, a_{ki}, a_{jk})$. + Partials direct. +- $(v_j, v_k) = $ (hyp, ideal): $\sigma_i = \zeta_{14}(a_{jk} - a_{ki}, a_{ij})$. + Then $\partial_{a_{jk}} = \partial_x \zeta_{14}$, + $\partial_{a_{ki}} = -\partial_x \zeta_{14}$, + $\partial_{a_{ij}} = \partial_y \zeta_{14}$. +- $(v_j, v_k) = $ (ideal, hyp): $\sigma_i = \zeta_{14}(a_{jk} - a_{ij}, a_{ki})$. + Sign pattern is the mirror image. +- $(v_j, v_k) = $ (ideal, ideal): $\sigma_i = \zeta_{15}(a_{jk} - a_{ij} - a_{ki})$. + Partials are $(\partial_x \zeta_{15}, -\partial_x \zeta_{15}, -\partial_x \zeta_{15})$ + along $(a_{jk}, a_{ij}, a_{ki})$. + +These nine sub-branches are the bulk of the per-face symbolic work. + +### 4.4 Chain rule for $\alpha_{ij}$ in the $v_i$-branch + +By (A.v_i) and (Z) evaluated at $(s_i, s_{ij}, s_{ik})$, +$$ +d\alpha_{ij} \;=\; +\zeta_x|_{(s_i, s_{ij}, s_{ik})} ds_i +\;+\; \zeta_y|_{(s_i, s_{ij}, s_{ik})} ds_{ij} +\;+\; \zeta_z|_{(s_i, s_{ij}, s_{ik})} ds_{ik} , +\tag{A1} +$$ +and each $ds_i, ds_{ij}, ds_{ik}$ expands via §4.3. The result is a +linear combination of the six face DOFs with closed-form coefficients. + +For the $v_j$-branch swap the roles $(i \leftrightarrow j)$; the +formulas are identical up to a permutation of $\sigma$-indices. + +### 4.5 The recursive $v_k$-branch + +When $v_i, v_j$ are both ideal but $v_k$ is hyper-ideal, the code returns +$$ +\alpha_{ij} \;=\; \pi - \alpha_{jk}' - \beta_j , +$$ +where $\alpha_{jk}'$ is a recursive call into the $v_j$-branch (since the +cyclic role rotation $\, (i, j, k) \mapsto (j, k, i)$ makes the second- +position vertex $v_k$, which is hyper-ideal, trigger the $v_i$-branch on +the recursion). Differentiating, +$$ +d\alpha_{ij} \;=\; -\, d\alpha_{jk}' \;-\; d\beta_j , +\tag{A2} +$$ +so we obtain $d\alpha_{ij}$ by computing $d\alpha_{jk}'$ via (A1) (with +the *rotated* DOF identification) and subtracting $d\beta_j$ from §4.1. +This single-level recursion is **finite** because the recursive call +descends into the $v_i$-branch (whose $v_i$ is now the original $v_k$, +which is hyper-ideal by hypothesis); see the comment "never more than one +level deep" in `hyper_ideal_geometry.hpp` line 106. + +### 4.6 The all-ideal branch (closed form) + +When all three vertices are ideal, +$\alpha_{ij} = \tfrac{1}{2}(\pi + \beta_k - \beta_i - \beta_j)$. Hence +$$ +d\alpha_{ij} \;=\; \tfrac{1}{2}(d\beta_k - d\beta_i - d\beta_j) , +\tag{A3} +$$ +with each $d\beta_\bullet$ from (B1). In this case the DOF vector +collapses to $(a_{12}, a_{23}, a_{31})$ (the $b$'s do not exist), and the +six-by-six block reduces to a non-trivial $3 \times 3$ block embedded +along the $a$-axes. + +--- + +## 5. The per-face $6 \times 6$ block + +### 5.1 Layout + +Order the inputs and outputs of `face_angles_from_local_dofs` as +$$ +x_{\text{loc}} \;=\; (b_1, b_2, b_3, a_{12}, a_{23}, a_{31})^\top, +\qquad +y_{\text{loc}} \;=\; (\beta_1, \beta_2, \beta_3, \alpha_{12}, \alpha_{23}, \alpha_{31})^\top . +$$ +The face Jacobian is +$$ +J(f) \;=\; \frac{\partial y_{\text{loc}}}{\partial x_{\text{loc}}} + \;=\; +\begin{pmatrix} +J^{\beta b} & J^{\beta a} \\ +J^{\alpha b} & J^{\alpha a} +\end{pmatrix} \in \mathbb{R}^{6 \times 6} . +$$ +The $3 \times 3$ sub-blocks are: + +- $J^{\beta b}_{ij} = \partial \beta_i / \partial b_j$, computed from (B1) + using the $L^b$ coefficients from §3.4. Sparse: only the two edges + $(i, \cdot)$ adjacent to vertex $i$ contribute the $b_j$ derivative + (via $\partial \ell_{ij} / \partial b_j$), so $J^{\beta b}_{ij}$ has at + most two nonzero edge-length terms per $(i, j)$ pair. +- $J^{\beta a}_{ie} = \partial \beta_i / \partial a_e$, also from (B1): + exactly **one** $\zeta$-derivative slot per $(i, e)$ pair, since each + $a_e$ affects exactly one edge length and each $\ell_e$ enters $\beta_i$ + in exactly one slot. +- $J^{\alpha b}_{e j} = \partial \alpha_e / \partial b_j$: computed from + the branch-appropriate identity among (A1), (A2), (A3). For (A1), only + $s_{ij}$ and $s_{ik}$ depend on $b$'s, so the row collapses to + $\zeta_y \cdot \partial_{b_j} s_{ij} + \zeta_z \cdot \partial_{b_j} s_{ik}$. +- $J^{\alpha a}_{e e'} = \partial \alpha_e / \partial a_{e'}$: this is the + densest block; both $s_i$ and $s_{ij}, s_{ik}$ depend on $a$'s, so all + three $\zeta$-slots contribute. + +### 5.2 Closed-form formulas, hyper-ideal triangle (all $v_i$ variable) + +Let $\beta_1 = \zeta(\ell_{12}, \ell_{31}, \ell_{23})$ (and cyclic), and +write the $\zeta$-partials at vertex $1$ as +$\zeta_x^{(1)} := \zeta_x(\ell_{12}, \ell_{31}, \ell_{23})$, etc.; and the +length partials from (Z13) as +$L^b_{e, i}, L^a_e$ etc. Then explicitly: +$$ +\begin{aligned} +\frac{\partial \beta_1}{\partial b_1} &= + \zeta_x^{(1)} L^b_{12,1} \;+\; \zeta_y^{(1)} L^b_{31,1} , + &\quad +\frac{\partial \beta_1}{\partial b_2} &= \zeta_x^{(1)} L^b_{12,2} , \\ +\frac{\partial \beta_1}{\partial b_3} &= \zeta_y^{(1)} L^b_{31,3} , + &\quad +\frac{\partial \beta_1}{\partial a_{12}} &= \zeta_x^{(1)} L^a_{12} , \\ +\frac{\partial \beta_1}{\partial a_{31}} &= \zeta_y^{(1)} L^a_{31} , + &\quad +\frac{\partial \beta_1}{\partial a_{23}} &= \zeta_z^{(1)} L^a_{23} . +\end{aligned} +\tag{Bblock} +$$ +(Note: $\beta_1$ does **not** depend on $a_{23}$ via the $b$'s — only the +direct $\ell_{23}$-dependence through $\zeta_z^{(1)}$ — so the last line +is the only $a_{23}$-coupling of $\beta_1$.) + +The full $J^{\beta b}$ row for vertex $1$ is therefore: +$$ +J^{\beta b}_{1,*} \;=\; +\big(\; + \zeta_x^{(1)} L^b_{12,1} + \zeta_y^{(1)} L^b_{31,1} ,\;\; + \zeta_x^{(1)} L^b_{12,2} ,\;\; + \zeta_y^{(1)} L^b_{31,3} +\;\big) , +$$ +and similarly for rows 2 and 3 by cyclic permutation. + +### 5.3 Closed-form formulas, dihedral row in $v_1$-branch + +Let $\beta := \alpha_{12} = \zeta(s_1, s_{12}, s_{13})$ (in the $v_1$ +branch, the $i$ of $\sigma$ is vertex 1, and the two edges through +vertex 1 are $12$ and $31 = 13$). With $\zeta_x^{(\alpha_{12})}$ etc. +denoting the three partials of $\zeta$ evaluated at $(s_1, s_{12}, s_{13})$: +$$ +\frac{\partial \alpha_{12}}{\partial \xi} + \;=\; \zeta_x^{(\alpha_{12})} \frac{\partial s_1}{\partial \xi} + \;+\; \zeta_y^{(\alpha_{12})} \frac{\partial s_{12}}{\partial \xi} + \;+\; \zeta_z^{(\alpha_{12})} \frac{\partial s_{13}}{\partial \xi} , +$$ +where $\xi$ ranges over the six face DOFs. Substituting +$s_1 = \zeta_{13}(a_{12}, a_{31}, a_{23})$ (assuming all +$v_j, v_k$ hyper-ideal so we are in the $(\text{hyp, hyp})$ sub-branch of +$\sigma_1$), +$s_{12} = \zeta_{13}(a_{12}, b_1, b_2)$, +$s_{13} = \zeta_{13}(a_{31}, b_1, b_3)$, the six entries are: +$$ +\begin{array}{l|l} +\xi & \partial \alpha_{12} / \partial \xi \\\hline +b_1 & \zeta_y^{(\alpha_{12})} \partial_y \zeta_{13}|_{s_{12}} + \zeta_z^{(\alpha_{12})} \partial_y \zeta_{13}|_{s_{13}} \\ +b_2 & \zeta_y^{(\alpha_{12})} \partial_z \zeta_{13}|_{s_{12}} \\ +b_3 & \zeta_z^{(\alpha_{12})} \partial_z \zeta_{13}|_{s_{13}} \\ +a_{12} & \zeta_x^{(\alpha_{12})} \partial_x \zeta_{13}|_{s_1} + \zeta_y^{(\alpha_{12})} \partial_x \zeta_{13}|_{s_{12}} \\ +a_{31} & \zeta_x^{(\alpha_{12})} \partial_y \zeta_{13}|_{s_1} + \zeta_z^{(\alpha_{12})} \partial_x \zeta_{13}|_{s_{13}} \\ +a_{23} & \zeta_x^{(\alpha_{12})} \partial_z \zeta_{13}|_{s_1} +\end{array} +\tag{Ablock} +$$ +The rows for $\alpha_{23}$ and $\alpha_{31}$ follow by cyclic permutation +$(1, 2, 3) \mapsto (2, 3, 1) \mapsto (3, 1, 2)$. + +### 5.4 Symmetry check (Schläfli) + +By (S2), the $6 \times 6$ block $H(f) = J(f)$ — interpreted as a +*Hessian-of-energy* block via the Springborn identification +$\beta_i = \partial U(f)/\partial b_i$, $\alpha_e = \partial U(f)/\partial a_e$ — +satisfies $J(f) = J(f)^\top$. Equivalently the four sub-blocks obey +$$ +J^{\beta b} = (J^{\beta b})^\top, \quad J^{\alpha a} = (J^{\alpha a})^\top, +\quad J^{\alpha b} = (J^{\beta a})^\top . +\tag{Sym} +$$ +This is a powerful numerical sanity check: any analytic formula that +violates (Sym) by more than rounding error has a derivation error +somewhere in §3–§4. + +The off-diagonal cross-check $\partial \alpha_{12} / \partial b_3 \;\stackrel{!}{=}\; \partial \beta_3 / \partial a_{12}$ +is particularly instructive: the left-hand side is computed via the +$\sigma$-triangle at vertex 1 (last row of `Ablock`-style table for +$\alpha_{12}$), the right-hand side via the auxiliary triangle at vertex 3. +The two routes coincide only because of (S2). + +### 5.5 Scatter to global Hessian + +Once $J(f)$ is built, the scatter step is **identical** to the existing +`hyper_ideal_hessian_block_fd`: for each $(i, j) \in \{1, \dots, 6\}^2$, +add $J_{ij}(f)$ to $H[\text{glb}(i), \text{glb}(j)]$, where $\text{glb}$ +maps the local face slot to its global DOF index via `v_idx` / `e_idx`. +Pinned slots (`-1`) are skipped, exactly as in the existing block-FD +implementation. + +--- + +## 6. Acceptance criteria + +These mirror `doc/roadmap/research-track.md` Phase 9b-analytic: + +1. **Per-case derivative cross-checks against block-FD.** For each of the + four building blocks $(\partial \beta / \partial \ell, \partial \ell / \partial \xi, + \partial \sigma / \partial \xi, \partial \alpha / \partial \xi)$, sample + 100 random DOF vectors in + - all-hyper-ideal regime ($v_1, v_2, v_3$ variable), + - one-ideal regime (each of three rotations of $v_i$ pinned), + - two-ideal regime (each of three rotations of two $v_i$ pinned). + Require $\max_{ij} |J_{ij}^{\text{analytic}} - J_{ij}^{\text{FD}}| \le 10^{-6}$ + with central-difference step $\varepsilon = 10^{-5}$. + +2. **Schläfli symmetry (gauge-free).** For each face $f$ at random $x$, + verify $\|J(f) - J(f)^\top\|_\infty \le 10^{-10}$. Symbolically this is + automatic by (S2); numerically it pins down sign/branch bugs. + +3. **Gauge null space.** Let $\mathbf{1}_b$ be the constant-$b$ Möbius + dilation mode (i.e. the vector with $1$ in every $b$-slot and $0$ in + every $a$-slot). Require $\|H \mathbf{1}_b\|_\infty \le 10^{-10}$ on + the all-hyper-ideal mesh after assembly (this is the only Möbius mode + that survives in the closed-surface case). + +4. **PSD on the interior.** For random DOFs inside the convex domain + (Springborn 2020 §4.3), verify $\lambda_{\min}(H) \ge -10^{-12}$ + (numerical zero). Closed-form: $E$ is convex on the admissible cone + by Springborn 2020 Theorem 4.2. + +5. **Speed-up.** Wall-clock benchmark on tetrahedron, octahedron, + icosahedron, and a 200-vertex genus-2 test mesh: analytic Hessian + assembly $\ge 3\times$ faster than `hyper_ideal_hessian_block_fd`, + with target $\sim 6\times$ asymptotically (one face costs + $O(6 \cdot 36)$ FD evaluations of `face_angles_from_local_dofs` + currently; analytic costs $\sim 100$ scalar transcendentals per face). + +--- + +## 7. Implementation outline + +Header: `code/include/hyper_ideal_hessian_analytic.hpp` (new). + +```cpp +// hyper_ideal_hessian_analytic.hpp — Phase 9b-analytic +// +// Closed-form 6x6 per-face Hessian via Schläfli + chain rule. +// See doc/math/hyperideal-hessian-derivation.md for the derivation. + +#include "hyper_ideal_geometry.hpp" +#include "hyper_ideal_functional.hpp" + +namespace conformallab { + +// Building-block derivatives (Section 3 of the note). +struct ZetaPartials { double dx, dy, dz; }; // d zeta(x, y, z) +struct Zeta13Partials { double dx, dy, dz; }; // d zeta13 +struct Zeta14Partials { double dx, dy; }; // d zeta14(x, y) + +inline ZetaPartials d_zeta (double x, double y, double z, double beta); +inline Zeta13Partials d_zeta13(double x, double y, double z, double L); +inline Zeta14Partials d_zeta14(double x, double y, double L); +inline double d_zeta15(double x); // returns tanh(L/2) + +// Edge-length differential coefficients (table in 3.4). +struct EdgeLenDiff { + double dbi; // d l_ij / d b_i (zero if v_i ideal) + double dbj; // d l_ij / d b_j (zero if v_j ideal) + double daij; // d l_ij / d a_ij +}; +EdgeLenDiff d_lij(double bi, double bj, double aij, bool vi, bool vj); + +// 6x6 face Jacobian: rows (beta1, beta2, beta3, alpha12, alpha23, alpha31), +// cols (b1, b2, b3, a12, a23, a31). Same input contract as +// face_angles_from_local_dofs. +struct FaceJacobian6 { double M[6][6]; }; + +FaceJacobian6 analytic_face_jacobian( + double b1, double b2, double b3, + double a12, double a23, double a31, + bool v1b, bool v2b, bool v3b); + +// Drop-in replacement for hyper_ideal_hessian_block_fd. +// Same scatter loop; only the inner block changes from FD to closed form. +Eigen::SparseMatrix hyper_ideal_hessian_analytic( + const ConformalMesh& mesh, + const std::vector& x, + const HyperIdealMaps& m); + +} // namespace conformallab +``` + +The body of `analytic_face_jacobian` follows the case dispatch: + +```text +1. clamp inputs same as face_angles_from_local_dofs +2. compute l12, l23, l31 via lij; bail out on triangle-inequality break +3. compute beta1, beta2, beta3 via zeta +4. for each edge e, compute EdgeLenDiff via d_lij(...) +5. for each vertex i, fill row i of J^{beta b} and J^{beta a} + using (B1) with d_zeta(...) evaluated at the angle of vertex i +6. for each edge e, dispatch on (v_i, v_j, v_k) and compute alpha_e row: + - (v_i hyp branch): use (A.v_i), expand s_i, s_ij, s_ik via 4.3 + - (v_j hyp branch): symmetric + - (v_k recursive branch): compute alpha_jk row first, then + subtract beta_j row and negate (A2) + - (all-ideal branch): combine three beta rows via (A3) +7. assert |J - J^T|_inf < 1e-10 (DEBUG only) +8. return J +``` + +The global assembler then reuses the existing scatter loop from +`hyper_ideal_hessian_block_fd`; only the inner *per-face block* +construction changes. + +--- + +## 8. References + +- **Schläfli, L.** (1858/60). *On the multiple integral + $\int dx\, dy \cdots dz$ whose limits are + $p_1 = a_1 x + b_1 y + \cdots + h_1 z > 0$, $p_2 > 0, \dots, p_n > 0$ + and $x^2 + y^2 + \cdots + z^2 < 1$.* Quarterly Journal of Pure and + Applied Mathematics 2, 269–301; 3, 54–68, 97–108. Reprinted in + *Gesammelte Mathematische Abhandlungen*, Band I, Birkhäuser 1950. + → first-order Schläfli identity, equation (S1) above. + +- **Milnor, J.** (1982). *Hyperbolic geometry: the first 150 years.* + Bulletin AMS 6, 9–24. → modern restatement of (S1) and discussion of + its role in hyperbolic volume. + +- **Vinberg, E. B.** (ed.) (1993). *Geometry II: Spaces of constant + curvature.* Encyclopaedia of Mathematical Sciences vol. 29, Springer. + Ch. 7 covers the polyhedral Schläfli identity and its derivation via + the Gauss–Bonnet formula for hyperbolic polytopes. + +- **Cho, Y. & Kim, H.** (1999). *On the volume formula for hyperbolic + tetrahedra.* Discrete & Computational Geometry 22, 347–366. + → explicit $\partial \alpha / \partial a$, $\partial \beta / \partial b$, + etc. for hyperbolic tetrahedra (used in §2.2 above). + +- **Rivin, I.** (1994). *Euclidean structures on simplicial surfaces and + hyperbolic volume.* Annals of Math. 139, 553–580. → second-order + Schläfli (S2) in the convex-polyhedron setting. + +- **Buser, P.** (1992). *Geometry and Spectra of Compact Riemann + Surfaces.* Birkhäuser. Thm 2.4.1 — right-angled hexagon identities + (used in §3.1 to interpret $\zeta_{13}$). + +- **Glickenstein, D.** (2011). *Discrete conformal variations and scalar + curvature on piecewise flat manifolds.* J. Diff. Geom. 87, 201–238. + Equation (4.6) and §4 — cone-vertex variant of the angle-length + derivative formulas (mirror Hessian for inversive distance). + +- **Springborn, B.** (2020). *Ideal Hyperbolic Polyhedra and Discrete + Uniformization.* Discrete & Computational Geometry 64, 63–108. + §4 — hyper-ideal energy and its gradient/Hessian; §4.3 — convexity + (PSD criterion). + +- **Bobenko, A. I., Pinkall, U. & Springborn, B.** (2015). *Discrete + conformal maps and ideal hyperbolic polyhedra.* Geometry & Topology + 19, 2155–2215. — earlier exposition of the truncated-tetrahedron + construction that `face_angles_from_local_dofs` realises. + +--- + +## Appendix A. Sign conventions and pitfalls + +A few places where the derivation interacts subtly with the code: + +- `zeta(l_jk, l_ki, l_ij)` returns the interior angle **at the vertex + between sides $l_{jk}$ and $l_{ki}$**, hence opposite to $l_{ij}$. So + when reading off $\zeta_x, \zeta_y, \zeta_z$ for $\beta_1 = + \zeta(\ell_{12}, \ell_{31}, \ell_{23})$, the slot $x$ matches + $\ell_{12}$, $y$ matches $\ell_{31}$, $z$ matches $\ell_{23}$. The + Cho–Kim formulas (Z) must be applied with **exactly this slot mapping**. + +- `sigma_ij(a_{ij}, b_i, b_j, v_j)` — when $v_j$ is ideal, the call is + `zeta14(-a_{ij}, b_i)`, so the partial w.r.t. $a_{ij}$ carries a + **negative** sign. This is a frequent source of off-by-sign bugs. + +- The recursive $v_k$-branch in `alpha_ij` returns $\pi - \alpha_{jk}' - \beta_j$, + **not** $\pi - \alpha_{jk}' - \beta_k$. This reflects the geometry of + the truncated tetrahedron and is *not* the cyclically-symmetric + formula one might guess; see Springborn 2020 Fig. 3. + +- The "all-ideal" closed form + $\alpha_{ij} = \tfrac{1}{2}(\pi + \beta_k - \beta_i - \beta_j)$ + is dual: $\alpha + \beta$ on each ideal vertex always sums to $\pi/2$ + for a planar triangle, so the dihedral degenerates to the half-angle + identity. Differentiating loses any explicit $\ell$ dependence, hence + (A3) is the simplest of the four branches. + +## Appendix B. Numerical guard zones + +The block-FD currently used as reference clamps: +- $b_i < 0 \;\to\; 0.01$ on variable vertices, +- $a_e < 0 \;\to\; 0$ on edges between two variable vertices. + +The analytic kernel must reproduce these clamps *before* computing any +derivatives; otherwise the FD-vs-analytic cross-check (Acceptance §1) +will fail at boundary inputs. Inside the clamped region, all +$\sinh \ell, \sin \beta$ denominators in (Z), (Z13), (Z14), (Z15) are +bounded away from zero, so the closed-form expressions are numerically +stable. The degenerate triangle-inequality branches in +`face_angles_from_local_dofs` (lines 174–183) set angles to $0$ or +$\pi$, which makes some $\zeta$-partials infinite; in those branches the +Jacobian is undefined and we return the all-zero block (matching the +non-smooth boundary of the domain). This is also the policy of the +block-FD reference, so the cross-check stays consistent. diff --git a/doc/roadmap/porting-status.md b/doc/roadmap/porting-status.md new file mode 100644 index 0000000..4018d99 --- /dev/null +++ b/doc/roadmap/porting-status.md @@ -0,0 +1,247 @@ +# Porting status overview + +> **Snapshot date:** 2026-05-22 (commit graph at v0.9.0 + 2 open PRs). +> +> **Audience.** External collaborators evaluating whether to use, extend, +> or contribute to conformallab++. This document is the **operational +> truth** of which Java mathematics is reachable from C++ today, with +> what API, and with which caveats. For the longer plan see +> [`phases.md`](phases.md), for the formal port-vs-research classification +> see [`research-track.md`](research-track.md), for the cross-reference +> with the Java original see [`java-parity.md`](java-parity.md). + +--- + +## 1. The 25 000 lines of Java in one table + +``` +de.varylab.discreteconformal +├── functional/ ~3 000 LoC ████████████░░░░ Phase 3 + 9a ✅ +├── unwrapper/ ~3 500 LoC ████████░░░░░░░░ Phase 5/6/7 (basic 3 modes) ✅ +│ residual: Stereographic, CircleDomain, Koebe, CirclePattern +├── util/ ~5 000 LoC ████░░░░░░░░░░░░ partial — main utilities ported +│ residual: Homology, Holomorphic / Harmonic forms, +│ FundamentalPolygon, Surgery, Cutting +├── heds/ (Co* data struct) ~3 000 LoC ░░░░░░░░░░░░░░░░ REPLACED by CGAL::Surface_mesh ❌ intentional +├── plugin/ ~7 000 LoC ░░░░░░░░░░░░░░░░ SKIP (Swing UI, GLSL viewer) ❌ intentional +├── numerics/ (PETSc glue) ~1 500 LoC ░░░░░░░░░░░░░░░░ REPLACED by Eigen ❌ intentional +├── math/, logging/, ... ~2 000 LoC ░░░░░░░░░░░░░░░░ REPLACED by std::* / Eigen ❌ intentional +└── total ~25 000 LoC +``` + +**Coverage of the mathematics:** about **4 500 of the 11 000 LoC** of +genuinely mathematical Java code is portrayed in C++ — **~41% of the +maths, ~75% of the "core math we actually want"** (Phase 1–7 + 9a + 9b +together cover the algorithmically central mass). + +The remaining ~6 500 LoC of port-worthy maths is catalogued in +[`java-parity.md`](java-parity.md) and broken into Phases 9c / 10a / +10b / 10b' / 11. + +--- + +## 2. Five DCE models — operational status matrix + +The library supports **five** discrete-conformal-equivalence models on +the same `CGAL::Surface_mesh

`. All can be used from both the +legacy API (`code/include/*.hpp`) and the CGAL public API +(``). + +| Model | Native space | DOF location | Java port? | Hessian | Newton | CGAL entry | UV out | +|------------------------|---------------|--------------|-----------|------------------|-----------|---------------------------------------------|----------| +| **Euclidean** | ℝ² | vertex | EuclideanCyclic 530 LoC | analytic (cot-Laplace) | ✅ | `discrete_conformal_map_euclidean` | ✅ Point_2 | +| **Spherical** | S² | vertex | Spherical 458 LoC | analytic | ✅ | `discrete_conformal_map_spherical` | ✅ Point_3 | +| **Hyper-ideal** | H² (Poincaré) | vertex + edge | HyperIdeal 311 LoC (no Hess) | block-FD (Phase 9b, 96× speed-up); analytic planned | ✅ | `discrete_conformal_map_hyper_ideal` | ✅ Point_2 | +| **CP-Euclidean** (BPS) | face circles | **face** | CPEuclidean 260 LoC | analytic 2×2-per-edge | ✅ | `discrete_circle_packing_euclidean` | ⛔ N/A | +| **Inversive Distance** | vertex circles| vertex | ❌ no Java (Luo 2004 + Glickenstein 2011 from literature) | FD (analytic planned) | ✅ | `discrete_inversive_distance_map` | ⛔ pending | + +Total: 250+ tests covering all five models, 0 skipped. Per-suite +breakdown: [`doc/api/tests.md`](../api/tests.md). + +### What "UV out" means + +The `output_uv_map(pmap)` named parameter, when supplied, runs the +relevant `*_layout()` step internally and writes per-vertex coordinates +into `pmap` — no separate user code needed. See +[`doc/tutorials/add-output-uv-map.md`](../tutorials/add-output-uv-map.md). + +* **Euclidean / Spherical / HyperIdeal**: shipped. +* **CP-Euclidean**: conceptually not applicable (face-based packing + produces face-circles, not per-vertex UV — needs a separate + `output_circle_map` design). +* **Inversive Distance**: requires `inversive_distance_layout()` using + Luo's edge-length formula `ℓ² = r_i² + r_j² + 2 I_ij r_i r_j`. + ~3 days work, on the wishlist. + +--- + +## 3. Topology infrastructure + +| Feature | Status | Notes | +|------------------------------------------------|:-----:|-------| +| Mesh I/O (OFF / OBJ / PLY) | ✅ | `mesh_io.hpp` via CGAL::IO | +| Gauss–Bonnet check & enforce | ✅ | `gauss_bonnet.hpp` | +| Cut graph (tree-cotree, Erickson–Whittlesey) | ✅ | `cut_graph.hpp`, produces 2g seam edges | +| Layout / embedding (priority-BFS trilateration) | ✅ | ℝ² / S² / Poincaré disk | +| Halfedge UV (seam-aware texture atlas) | ✅ | `layout.hpp` | +| Möbius holonomy SU(1,1) | ✅ | for hyperbolic mode, genus 1 | +| Period matrix τ ∈ ℍ + SL(2,ℤ) reduction | ✅ | genus 1 only | +| Fundamental domain (parallelogram) | ✅ | genus 1 only | +| Fundamental domain (4g-polygon) | 🔲 | Phase 9c — biggest remaining Java-port item | +| Real mesh cuts (CuttingUtility ops) | 🔲 | Phase 9c prerequisite | + +--- + +## 4. Solver infrastructure + +| Feature | Status | Notes | +|-------------------------------------------------|:-----:|-------| +| Newton with line search | ✅ | `newton_solver.hpp`, all five solvers | +| SimplicialLDLT + SparseQR fallback | ✅ | gauge-singular meshes handled automatically | +| Block-FD Hessian framework | ✅ | shipped for HyperIdeal (Phase 9b); 96× speed-up | +| Analytic Hessian via Schläfli | 🔲 | Phase 9b-analytic; derivation in [`hyperideal-hessian-derivation.md`](../math/hyperideal-hessian-derivation.md) | + +--- + +## 5. CGAL public API (``) + +After Phase 8a MVP + 8b-Lite (both shipped in v0.9.0): + +| Public header | Provides | +|--------------------------------------------------------------|----------| +| `` | `discrete_conformal_map_euclidean / _spherical / _hyper_ideal` + `Conformal_map_result` | +| `` | `Default_cp_euclidean_traits` + `discrete_circle_packing_euclidean` + `Circle_packing_result` | +| `` | `Default_inversive_distance_traits` + `discrete_inversive_distance_map` | +| `` | Thin re-exports of layout functions in the `CGAL::` namespace | +| `` | `ConformalMapTraits` concept + `Default_conformal_map_traits` | +| `` (private) | Named-parameter tags | + +Named parameters available: + +| Parameter | Effect | +|------------------------------------|--------| +| `vertex_curvature_map(pmap)` | Per-vertex Θ targets; if omitted, "natural-theta" makes x = 0 the equilibrium | +| `fixed_vertex_map(pmap)` | Gauge pin; if omitted, first vertex is pinned | +| `gradient_tolerance(eps)` | Newton stop threshold (default 1e-10) | +| `max_iterations(n)` | Newton iteration cap (default 200) | +| `output_uv_map(pmap)` | Run layout + write per-vertex coordinates | +| `normalise_layout(flag)` | Canonical post-layout normalisation (PCA / north-pole / Möbius) | + +**Known limitations** (documented for the meeting): + +1. **No chaining yet.** `.gradient_tolerance(eps).max_iterations(n)` does + not work — pass one named parameter per call. Fix is mechanical + (~2 days, generate CGAL-style chaining helpers via macros). + +2. **Generic FaceGraph not yet supported.** All Default traits + specialise on `CGAL::Surface_mesh

` only. Adding `Polyhedron_3`, + `OpenMesh-adapter`, `pmp::SurfaceMesh` is Phase 8a.2 — speculative, + waiting on a concrete user request. + +3. **CGAL submission-readiness incomplete.** User_manual / + PackageDescription.txt / CGAL-format tests are not yet written + (Phase 8c / 8d). Required for upstream submission, not for current + users. + +--- + +## 6. What's reachable from the Java original today + +Cross-reference between Java classes and their C++ port status. For +the full list see [`java-parity.md`](java-parity.md); the table below +gives the operational summary for a Java user wondering "where is X?". + +### Direct ports + +| Java | C++ | +|----------------------------------------------|-------------------------------------------------------------| +| `EuclideanCyclicFunctional` | `euclidean_functional.hpp` | +| `SphericalFunctional` | `spherical_functional.hpp` | +| `HyperIdealFunctional` | `hyper_ideal_functional.hpp` | +| `HyperIdealUtility`, `…Geometry` | `hyper_ideal_utility.hpp`, `hyper_ideal_geometry.hpp` | +| `HyperIdealHyperellipticUtility` | `hyper_ideal_visualization_utility.hpp` | +| `CPEuclideanFunctional` | `cp_euclidean_functional.hpp` | +| `Clausen` | `clausen.hpp` | +| `MeshIO` (jReality JRS) | `mesh_io.hpp` (CGAL::IO, OFF/OBJ/PLY) | +| `CuttingUtility::point_in_triangle_2d` | `geometry_utils.hpp` partial port | +| `ConvergenceUtility::circumradius` | `geometry_utils.hpp` | +| `HomologyUtility::compute_generators` | `cut_graph.hpp` (different algorithm: tree-cotree) | +| `PeriodMatrixUtility` (genus 1) | `period_matrix.hpp` | +| `UnwrapperUtility::EuclideanUnwrapper` | `newton_euclidean()` + `euclidean_layout()` | +| `UnwrapperUtility::SphericalUnwrapper` | `newton_spherical()` + `spherical_layout()` | +| `UnwrapperUtility::HyperbolicUnwrapper` | `newton_hyper_ideal()` + `hyper_ideal_layout()` | + +### Not-yet-ported, in roadmap + +| Java | Suggested phase | Lines | Notes | +|---------------------------------------------------|:----:|:----:|-------| +| `FundamentalPolygonUtility` + `CanonicalFormUtility` | 9c | 698+532 | genus > 1 fundamental domain | +| `CuttingUtility` + `SurgeryUtility` (full ops) | 9c | 584+217 | real mesh-cut operations | +| `DiscreteHarmonicFormUtility` | 10a | 657 | harmonic 1-forms | +| `DiscreteHolomorphicFormUtility` | 10a | 285 | holomorphic differentials | +| `CanonicalBasisUtility` | 10a | 337 | symplectic homology basis | +| `DualityUtility`, `HomologyUtility` (full) | 10a | 308+122 | primal/dual cohomology | +| `DiscreteRiemannUtility` | 10b | 186 | period matrix for genus > 1 | +| `HyperbolicCyclicFunctional` | 10b–c | 530 | hyperbolic energy (different from hyper-ideal) | +| `QuasiisothermicUtility` + `SinConditionApplication` | 10b | ~1200 | Lawson correspondence | +| `StereographicUnwrapper` | 10b' | 266 | S² → ℂ atlas | +| `KoebePolyhedron` | 10c | 321 | KAT circle packing | +| `MobiusCenteringFunctional` | 10c | 289 | sphere pre-processing | +| `CircleDomainUnwrapper` | 11c | 570 | multiply-connected planar regions | + +### Intentionally not ported + +* **`CoHDS` (half-edge data structure)** — replaced by `CGAL::Surface_mesh`. ~3 000 Java LoC saved. +* **`plugin/*` Swing UI** — out of scope; C++ project is a library, not an application. +* **`unwrapper/numerics/*` PETSc/Tao wrappers** — replaced by Eigen, which is header-only and bundled. +* **`math/CP1`, `math/Cn`, …** — redundant with `std::complex`, `Eigen::Matrix`. + +--- + +## 7. Things in C++ that the Java original does NOT have + +These are conformallab++ contributions beyond porting — the +"research extensions" track ([`research-track.md`](research-track.md)). + +| Item | Status | +|----------------------------------------------|:-----:| +| Block-FD HyperIdeal Hessian (96× speed-up) | ✅ shipped (Phase 9b) | +| HyperIdeal Hessian (any kind) | ✅ shipped — Java has `hasHessian()==false` | +| Inversive-Distance functional (Luo 2004) | ✅ shipped — implemented from paper, not from Java | +| `` public-API surface | ✅ shipped (Phase 8a + 8b-Lite) | +| Cross-validation framework (4 acceptance criteria) | ✅ documented in `add-inversive-distance.md` | +| Genus-2 test mesh + brezel2.obj scalability | ✅ shipped | +| Memory-safe layout via `halfedge_uv` storage | ✅ shipped | +| `doc/release-policy.md` formal release policy | ✅ shipped (PR #13) | +| Analytic HyperIdeal Hessian via Schläfli | 🔲 derivation written (`hyperideal-hessian-derivation.md`); implementation Phase 9b-analytic | +| Output UV map integrated into wrappers | ✅ shipped (PR #14) for 3 of 5 models | +| Full uniformisation for genus g ≥ 2 | 🔲 Phase 10c — fully new research | + +--- + +## 8. How to "use the library today" + +The fastest path for a mathematician evaluating the library: + +```bash +# Build + run the full test suite +git clone && cd ConformalLabpp +cmake -S code -B build -DWITH_CGAL_TESTS=ON +cmake --build build --target conformallab_cgal_tests -j$(nproc) +ctest --test-dir build -R "^cgal\." --output-on-failure + +# Try the CGAL API on a mesh +./build/examples/example_layout my_mesh.off +``` + +To experiment with a new functional, follow +[`doc/tutorials/add-inversive-distance.md`](../tutorials/add-inversive-distance.md). + +To extend the CGAL API with new named parameters, follow +[`doc/tutorials/add-output-uv-map.md`](../tutorials/add-output-uv-map.md). + +To optimise an FD Hessian, follow +[`doc/tutorials/block-fd-hessian.md`](../tutorials/block-fd-hessian.md). + +For the architectural commitments that constrain future contributions, +see [`doc/architecture/locked-vs-flexible.md`](../architecture/locked-vs-flexible.md). diff --git a/doc/tutorials/add-output-uv-map.md b/doc/tutorials/add-output-uv-map.md new file mode 100644 index 0000000..2334e0d --- /dev/null +++ b/doc/tutorials/add-output-uv-map.md @@ -0,0 +1,477 @@ +# Tutorial: The `output_uv_map` Named Parameter + +This tutorial documents the `output_uv_map` named parameter for the +CGAL-style entry functions of conformallab++: how to use it (the "happy +path"), how the CGAL named-parameter machinery threads a user-supplied +property map through to the layout step, and how to extend the same +pattern to add new named parameters to existing or new entry functions. + +**Target audience.** A working mathematician using the C++ API who is +comfortable with templates and CGAL property maps but wants a single +authoritative reference for the pattern. + +**Prerequisite reading.** + +- `code/include/CGAL/Conformal_map/internal/parameters.h` — tag definitions +- `code/include/CGAL/Discrete_conformal_map.h` — wiring in three entries +- `code/tests/cgal/test_cgal_phase8b_lite.cpp` — the OutputUvMap_* tests +- `doc/tutorials/add-inversive-distance.md` — sibling tutorial in this series + +--- + +## 1. What this tutorial is and isn't + +### The UX problem `output_uv_map` solves + +Before Phase 8b-Lite, computing a flattening required **two separate +calls**: the Newton wrapper, then a manual rebuild of `*_maps` plus a +replay of the wrapper's gauge/DOF choices, then `*_layout()`. Any drift +between the wrapper's gauge vertex and the caller's silently produced +inconsistent UVs. + +The `output_uv_map` named parameter collapses this into one call: + +```cpp +auto uv = mesh.add_property_map("uv").first; +CGAL::discrete_conformal_map_euclidean( + mesh, CGAL::parameters::output_uv_map(uv)); +// uv[v] now holds the flattened UV coordinate of vertex v. +``` + +This tutorial is **not** a derivation of the Newton solver (see +`add-inversive-distance.md`), not a guide to writing a new layout +routine, and not an introduction to CGAL named parameters in general +(see ``). + +--- + +## 2. Mathematical background + +The package implements **two mathematically distinct stages**: + +### Stage A — Newton solver (variational optimisation) + +For the Euclidean DCE functional (Springborn–Schmies–Bobenko 2008), the +solver finds the conformal scale factor vector `u ∈ ℝⁿ` that satisfies + +``` +G(u)_v := Θ_v − Σ_{T ∋ v} α_v(T; u) = 0 ∀ v ∈ V \ {pinned} +``` + +i.e. the actual cone angles match the target curvature `Θ_v` at every +non-pinned vertex. The gradient `G` is the variational derivative of +the convex energy + +``` +E(u) = ∫_0^1 ⟨ G(t u), u ⟩ dt (path integral on the closed Yamabe 1-form) +``` + +The output is **purely combinatorial-metric**: the converged `u_v` +determines an intrinsic flat metric `ℓ_ij(u) = exp(½(u_i + u_j)) · ℓ_ij^(0)` +on the abstract triangulation. **No embedding yet exists.** + +### Stage B — Layout (geometric realisation) + +The layout step is a separate **trilateration** turning intrinsic edge +lengths into a `Point_2` (or `Point_3`) per vertex: + +1. Pick a seed triangle `T₀`; place its three vertices canonically. +2. BFS over the dual graph from `T₀`. Each newly visited triangle has + two vertices already placed; the third is the unique solution of + two distance constraints (intersection of two circles). +3. The BFS uses a **priority queue keyed on BFS depth** + (`euclidean_layout` in `code/include/layout.hpp`) to avoid the + numerical drift of pure DFS on long thin meshes. + +Spherical / hyperbolic variants substitute spherical-cap or Poincaré- +disk trilateration respectively. + +Stage A is convex optimisation in ℝⁿ; Stage B is deterministic +realisation with no degrees of freedom. The named parameter chains them. + +--- + +## 3. The user-facing API + +### Happy path + +```cpp +#include +#include +#include + +using K = CGAL::Simple_cartesian; +using Mesh = CGAL::Surface_mesh; +using Vertex_i = Mesh::Vertex_index; + +Mesh mesh = /* load triangle mesh */; + +// 1. Allocate a writable property map on the mesh. +auto uv = mesh.add_property_map( + "v:uv", K::Point_2(0, 0)).first; + +// 2. Ask the wrapper to populate it. +auto res = CGAL::discrete_conformal_map_euclidean( + mesh, + CGAL::parameters::output_uv_map(uv)); + +// 3. Use the UVs. +if (res.converged) { + for (auto v : mesh.vertices()) + std::cout << v.idx() << ": " << uv[v] << '\n'; +} +``` + +### Supported entries + +| Entry | Value type per vertex | +|--------------------------------------|--------------------------------------------------| +| `discrete_conformal_map_euclidean` | `Point_2` — UV in ℝ² | +| `discrete_conformal_map_spherical` | `Point_3` — point on S² ⊂ ℝ³ | +| `discrete_conformal_map_hyper_ideal` | `Point_2` — point in the Poincaré disk (|p|≤1) | + +### NOT supported (yet) + +| Entry | Reason | +|----------------------------------------|-------------------------------------------------------------| +| `discrete_circle_packing_euclidean` | Face-based parametrisation — UV is per-face, not per-vertex | +| `discrete_inversive_distance_map` | Requires `inversive_distance_layout()` using Luo's `ℓ²` formula — not yet written | + +For the circle-packing entry a separate `output_circle_map` parameter +(face → (centre, radius)) is the natural extension; for the inversive- +distance entry the layout primitive must be written first. See §7. + +### Optional companion: `normalise_layout` + +```cpp +CGAL::parameters::normalise_layout(true) +``` + +When the layout writes into `output_uv_map`, this flag triggers a +canonical post-processing pass: + +- Euclidean: PCA — translate centroid to origin, rotate so the principal + axis is along `x`. +- Spherical: Rodrigues rotation so that the centroid is at the north pole. +- Hyperbolic: Möbius centring of the Poincaré disk. + +Default: `false`. Has no effect unless `output_uv_map` is also passed. + +### Current limitation: no chaining + +The CGAL machinery supports chaining named parameters via `.member()` +syntax for the standard CGAL tags (e.g. `geom_traits`, `vertex_point_map`). +**Package-local tags do not yet have member-function counterparts**, so + +```cpp +// DOES NOT COMPILE: +CGAL::parameters::output_uv_map(uv).normalise_layout(true); +``` + +is currently rejected. Pass one parameter per call instead. This is +the same limitation noted in `test_cgal_phase8b_lite.cpp`: + +```cpp +// Named-parameter chaining (`a.b().c()`) is not currently supported on +// the package-local tags; pass one parameter per call instead. +``` + +Adding member-function chaining is a separate Phase 8b extension — +requires writing a derived `Named_function_parameters` subclass that +exposes `.output_uv_map(...)`, `.normalise_layout(...)`, etc. + +--- + +## 4. Named-parameter mechanism — how it works + +CGAL's named parameters are a compile-time dispatch trick. Three pieces: + +### 4.1 The tag (in `internal_np`) + +```cpp +namespace CGAL::Conformal_map::internal_np { + enum output_uv_map_t { output_uv_map }; +} +``` + +A single-value `enum` is the lightest possible type-level marker. The +value `output_uv_map` is what callers will pass as the *key*; the type +`output_uv_map_t` is what `get_parameter` will use for the lookup. + +### 4.2 The user-facing helper (in `CGAL::parameters`) + +```cpp +template +auto output_uv_map(const PropertyMap& pmap) { + return CGAL::Named_function_parameters< + PropertyMap, + Conformal_map::internal_np::output_uv_map_t, + CGAL::internal_np::No_property // no "next" parameter in chain + >(pmap); +} +``` + +This wraps the property map in a `Named_function_parameters` object +tagged with `output_uv_map_t`. The third template argument is the +"previous" link in a chain; `No_property` means this is a singleton +pack. (Chaining would substitute the previous pack's type here.) + +### 4.3 The internal lookup (inside the entry function) + +The full read-and-act pattern from `Discrete_conformal_map.h`: + +```cpp +// Step 1 — try to extract the parameter. +auto uv_param = parameters::get_parameter( + np, Conformal_map::internal_np::output_uv_map); + +// Step 2 — at compile time, did we find one? +constexpr bool has_uv = !std::is_same_v< + decltype(uv_param), internal_np::Param_not_found>; + +// Step 3 — conditionally do the extra work. +if constexpr (has_uv) { + if (nr.converged) { + auto layout = ::conformallab::euclidean_layout(mesh, nr.x, maps); + + const bool do_norm = parameters::choose_parameter( + parameters::get_parameter(np, Conformal_map::internal_np::normalise_layout), + false); + if (do_norm) ::conformallab::normalise_euclidean(layout); + + for (auto v : mesh.vertices()) { + const auto& uv = layout.uv[v.idx()]; + put(uv_param, v, + typename Traits::Kernel::Point_2(uv.x(), uv.y())); + } + } +} +``` + +Key points: + +- `get_parameter` returns either the wrapped value or a sentinel + `Param_not_found`. The check is `std::is_same_v<…, Param_not_found>`. +- `choose_parameter(get_parameter(np, tag), default_value)` is the + one-liner form for scalar parameters with a default. +- `constexpr if` ensures that if the user did not pass `output_uv_map`, + the layout call is **not even instantiated** — no runtime cost, and + no need for the layout routine to be callable on the wrapper's + signature (e.g. it can require a Point_2 traits type that the user + hasn't supplied). +- The layout step is gated on `nr.converged` for a reason: trilateration + on a non-converged `u` produces garbage edge lengths and can throw + on degenerate triangles. The current contract is **layout iff Newton + converged**; on non-convergence the pmap is left at the value the + caller initialised it with. + +The same three-step pattern appears in `discrete_conformal_map_spherical` +(writing `Point_3`) and `discrete_conformal_map_hyper_ideal` (writing +`Point_2` in the Poincaré disk). + +--- + +## 5. How to extend the pattern + +Recipe for adding a new package-local named parameter — running example +`output_holonomy_map` (a per-edge map carrying the Möbius/translation +holonomy data that some downstream algorithms consume). + +### Step 1 — Add the tag in `internal/parameters.h` + +```cpp +namespace CGAL::Conformal_map::internal_np { + /// Property-map: edge_descriptor → 2x2 matrix of holonomy data. + /// Default: not populated. + enum output_holonomy_map_t { output_holonomy_map }; +} +``` + +Two-line block. Use a meaningful Doxygen `\internal` comment because +this is the canonical contract for the parameter's semantics. + +### Step 2 — Add the user-facing helper in `parameters.h` + +```cpp +namespace CGAL::parameters { + +/// `output_holonomy_map(pmap)` — write per-edge holonomy matrices +/// into `pmap`. Only meaningful for closed surfaces with non-trivial +/// π₁; ignored on simply-connected meshes. +template +auto output_holonomy_map(const PropertyMap& pmap) { + return CGAL::Named_function_parameters< + PropertyMap, + Conformal_map::internal_np::output_holonomy_map_t, + CGAL::internal_np::No_property + >(pmap); +} + +} // namespace CGAL::parameters +``` + +A Doxygen comment is **mandatory** here — this is the public surface. + +### Step 3 — Use `get_parameter + constexpr if` in the entry function + +```cpp +auto hol_param = parameters::get_parameter( + np, Conformal_map::internal_np::output_holonomy_map); +constexpr bool has_hol = !std::is_same_v< + decltype(hol_param), internal_np::Param_not_found>; +if constexpr (has_hol) { + auto H = ::conformallab::compute_holonomy(mesh, nr.x, maps); + for (auto e : mesh.edges()) + put(hol_param, e, H[e.idx()]); +} +``` + +### Step 4 — Add tests + +At minimum two: + +1. **Takes-effect test.** Pass the parameter; verify the side effect + actually happens (pmap is populated, values are non-default). +2. **Sanity-when-absent test.** Call the wrapper without the parameter; + verify the call still succeeds and produces the expected base result. + +Both patterns are exemplified in `test_cgal_phase8b_lite.cpp` (next +section). + +### Optional Step 5 — Companion flag + +If your new parameter has a Boolean tweak (analogous to +`normalise_layout`), follow exactly the same recipe — `enum X_t { X }` ++ helper + `choose_parameter(..., default)` at the read site. + +--- + +## 6. Tests — walkthrough of `test_cgal_phase8b_lite.cpp` + +The OutputUvMap_* tests live at the bottom of the file. Each codifies +one acceptance criterion. + +### 6.1 `OutputUvMap_Euclidean_PopulatesPmap` + +Two assertions: (a) every UV is finite, (b) at least one vertex has +moved off the origin. Test (b) rules out the false-positive where the +layout silently returned the pmap's default at every vertex. + +### 6.2 `OutputUvMap_Spherical_PopulatesXyz` + +Geometric round-trip: every output point must lie on the unit sphere. + +```cpp +const double r = std::sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z()); +EXPECT_NEAR(r, 1.0, 1e-6); +``` + +The spherical analogue of 6.1(b): the layout is correct iff the +realisation lands on S². + +### 6.3 `OutputUvMap_HyperIdeal_PointsInPoincareDisk` + +The natural Θ/θ targets do not always reach Newton equilibrium inside +the default `max_iterations(200)`. The test branches: + +```cpp +if (res.converged) { + for (auto v : mesh.vertices()) { + const double r2 = p.x()*p.x() + p.y()*p.y(); + EXPECT_LE(r2, 1.0 + 1e-6); + } +} +// else: layout was skipped by the wrapper's `if (nr.converged)` guard. +``` + +Documenting the wrapper contract ("layout iff converged") via the +test keeps it resilient to PRNG-sensitive Newton paths. + +### 6.4 `OutputUvMap_Absent_DoesNotRunLayout` + +Sanity test: omitting `output_uv_map` must not change the existing +behaviour. Catches the regression where adding the new `constexpr if` +perturbs the wrapper's result type or consumes extra work. + +### 6.5 `OutputUvMap_NormaliseLayout_TakesEffect` + +Because chaining is not yet implemented (§3), this exercises the toggle +indirectly: it runs the wrapper twice into two separate pmaps and +verifies both populate finite values. When chaining lands, upgrade this +to actually pass `normalise_layout(true)` on the second call and compare +UV bounding boxes (the normalised one must be axis-aligned and +centroid-zero). + +--- + +## 7. Adding output to the two missing modes + +### 7.1 CP-Euclidean (face-based) — proposed `output_circle_map` + +The circle-packing entry parametrises **faces**: each face `f` has a +radius `ρ_f`. The natural per-face output is `(centre, radius)`: + +```cpp +auto circles = mesh.add_property_map>("f:circle", ...).first; + +CGAL::discrete_circle_packing_euclidean( + mesh, CGAL::parameters::output_circle_map(circles)); +``` + +Effort: **1–2 days**. The face-graph BFS-packing is Stephenson's +`circlepack` algorithm; radii come from Newton, new work is the +geometric loop plus the §5 plumbing. + +### 7.2 Inversive-Distance — needs new layout routine first + +The inversive-distance functional yields edge lengths + +``` +ℓ_ij(u)² = r_i² + r_j² + 2 I_ij r_i r_j (Luo 2004 §3) +``` + +with `r_i = exp(u_i)`. Required work: write +`inversive_distance_layout()` in `code/include/layout.hpp` mirroring +`euclidean_layout()` but consuming the above `ℓ_ij(u)` instead of +`exp(½(u_i + u_j)) · ℓ_ij^(0)`. Then wire the parameter through +`discrete_inversive_distance_map` exactly as §4.3 shows. + +Effort: **~3 days**. Status: **TBD**. + +--- + +## 8. Acceptance checklist for a new named parameter + +Use this checklist when extending the package: + +- [ ] **Tag declared** in `code/include/CGAL/Conformal_map/internal/parameters.h` + (`enum X_t { X };`) with a `\internal` Doxygen comment specifying + the property-map key/value types and default behaviour. +- [ ] **User-facing helper** in the same file's `CGAL::parameters` + namespace, with a public Doxygen comment. +- [ ] **`get_parameter + constexpr if` pattern** in every entry function + that supports the parameter (do not leak the tag into entries that + ignore it). +- [ ] **At least two tests** in `test_cgal_phase8b_lite.cpp` (or a new + test file if introducing a new functional): a + *parameter-takes-effect* test and a *sanity-when-absent* test. +- [ ] **Limitation documented.** If the parameter is mode-specific + (e.g. only meaningful for closed meshes) or cannot be chained + with existing helpers, say so in the user-facing Doxygen and in + a comment near the read site. +- [ ] **No silent regressions.** If the parameter changes the wrapper's + observable behaviour even when absent (e.g. via a new default), + add a regression test pinning the old default. + +--- + +## 9. Further reading + +- Springborn, Schmies, Bobenko (2008). *Conformal equivalence of + triangle meshes.* SIGGRAPH. — Euclidean DCE functional. +- Bobenko, Pinkall, Springborn (2015). *Discrete conformal maps and + ideal hyperbolic polyhedra.* Geom. Topol. — hyper-ideal variant. +- Glickenstein (2011). *Discrete conformal variations and scalar + curvature on piecewise flat manifolds.* JDG 87(2). +- `` — upstream machinery. diff --git a/doc/tutorials/block-fd-hessian.md b/doc/tutorials/block-fd-hessian.md new file mode 100644 index 0000000..81778d5 --- /dev/null +++ b/doc/tutorials/block-fd-hessian.md @@ -0,0 +1,460 @@ +# Tutorial: The Per-Face Block-Finite-Difference Hessian (Phase 9b) + +This tutorial documents the **block-finite-difference Hessian** scheme +used for the hyper-ideal discrete-conformal functional of Bobenko– +Springborn–Schief / Springborn 2020. It is also a *pattern document*: +once the locality lemma below is verified for another variational +functional in the library, the same scaffolding can be reused almost +verbatim. + +The implementation lives in +[`code/include/hyper_ideal_hessian.hpp`](../../code/include/hyper_ideal_hessian.hpp); +the pure 6 → 6 face kernel it differentiates lives in +[`code/include/hyper_ideal_functional.hpp`](../../code/include/hyper_ideal_functional.hpp) +as `face_angles_from_local_dofs(...)`. + +> ## ⚠️ This is research, not a port +> +> The upstream Java reference +> `de.varylab.discreteconformal.functional.HyperIdealFunctional` +> declares +> +> ```java +> public boolean hasHessian() { return false; } // line 295–298 +> ``` +> +> i.e. the Java implementation supplies **no** Hessian — neither +> analytic nor numerical. Newton-time second-order behaviour in the +> Java pipeline is delegated to PETSc's BFGS approximation. Both the +> full-FD Hessian (Phase 4a) and the block-FD Hessian documented here +> (Phase 9b) are **conformallab++ additions beyond Java parity**. The +> mathematical justification is therefore taken directly from +> +> 1. **Springborn, B.** (2020). *Hyperbolic polyhedra and discrete uniformization.* +> Discrete & Computational Geometry 64, 63–108. §4 (the variational +> principle whose Hessian we are differentiating). +> 2. **Bobenko, A. I. & Springborn, B. A.** (2004). *Variational principles +> for circle patterns and Koebe's theorem.* Trans. AMS 356(2), 659–689. + +**Prerequisite:** familiarity with the hyper-ideal functional itself +(see [`doc/math/hyper-ideal.md`](../math/hyper-ideal.md)), and with the +generic functional-porting pattern of +[`doc/tutorials/add-inversive-distance.md`](add-inversive-distance.md). + +--- + +## What this tutorial is and isn't + +There are three legitimate ways to obtain a Hessian for a discrete- +conformal energy in this codebase: + +| Strategy | Complexity per call | Implementation cost | When to choose | +|---|---|---|---| +| **Full finite difference** | `O(n · F)` — perturb each global DOF, re-run the entire gradient | low (~30 LOC, generic over any functional) | small meshes; gold-standard correctness reference | +| **Per-face block FD** *(this tutorial)* | `O(F · 36)` — perturb each face's 6 local DOFs through the pure face kernel | moderate (~70 LOC, one per functional, requires the locality lemma) | production default for the hyper-ideal energy | +| **Full analytic Hessian** | `O(F)` — closed-form ∂(β,α)/∂(b,a) via Schläfli | high (multi-week derivation through `lᵢⱼ → ζ → β, α`); see Phase 9b-analytic | high-throughput pipelines, tight inner loops | + +The block-FD variant is **the sweet spot**: it inherits the genericity +and "trust" of finite differencing while exploiting the exact local +structure of the variational principle. We measured a **96.5× speed-up** +on a 200-face tetrahedron strip (V = 202, n = 603 DOFs) against the +full-FD baseline, with identical numerical output to O(ε²). + +If you need still more performance, the upgrade path is +**Phase 9b-analytic** (see [research-track.md §Phase 9b-analytic](../roadmap/research-track.md)), +which differentiates through the chain +`(bᵢ, aₑ) → lᵢⱼ → ζ → α, β` analytically using Schläfli's identity. +That is *another* ~6× over block-FD but at substantial implementation +cost. + +--- + +## Mathematical background — the per-face locality lemma + +The hyper-ideal energy `E(b, a)` is a sum over faces of a local function +of the **six DOFs touching that face**: + +* three vertex DOFs `(b₁, b₂, b₃)` — Penner-style edge-weight logarithms, +* three edge DOFs `(a₁₂, a₂₃, a₃₁)` — log-coshes of the hyperbolic + truncation lengths. + +The face contributes six output angles: + +* three interior angles `(β₁, β₂, β₃)` at the vertices, +* three dihedral angles `(α₁₂, α₂₃, α₃₁)` at the edges. + +These six numbers are obtained by the pure function + +```cpp +struct FaceAngleOutputs { + double beta1, beta2, beta3; // interior angles at v₁,v₂,v₃ + double alpha12, alpha23, alpha31; // dihedral angles at e₁₂,e₂₃,e₃₁ +}; + +FaceAngleOutputs face_angles_from_local_dofs( + double b1, double b2, double b3, + double a12, double a23, double a31, + bool v1b, bool v2b, bool v3b); +``` + +(located in `hyper_ideal_functional.hpp:152`). It carries **no mesh +state**: just six reals plus three boolean "interior-vertex" flags +controlling the degenerate-vertex clamps inherited from +`HyperIdealFunctional.java` lines 122–127. + +### Gradient decomposition + +The global gradient is the angle-defect / Schläfli-type sum + +``` +G_{b,v} = Σ_{f ∋ v} β_v(f) − Θ_v (Springborn 2020 eq. 4.6) +G_{a,e} = Σ_{f ∋ e} α_e(f) − θ_e +``` + +where `Θ_v` is the prescribed interior angle sum at vertex `v` and +`θ_e` is the prescribed dihedral at edge `e`. Each term in either sum +depends on **only the six local DOFs of one face**. + +### Hessian decomposition (the locality lemma) + +Differentiating once more: + +``` +∂G_{b,v}/∂y = Σ_{f ∋ v} ∂β_v(f) / ∂y +∂G_{a,e}/∂y = Σ_{f ∋ e} ∂α_e(f) / ∂y +``` + +and `∂β_v(f)/∂y` (resp. `∂α_e(f)/∂y`) is non-zero **only if `y` is one +of the six local DOFs of face `f`**. Hence the global Hessian is + +``` +H[x, y] = Σ_{f : x, y ∈ local(f)} J_f[row(x), col(y)] +``` + +where `J_f ∈ ℝ^{6×6}` is the local Jacobian of the map + +``` +(β₁, β₂, β₃, α₁₂, α₂₃, α₃₁) = Φ_f(b₁, b₂, b₃, a₁₂, a₂₃, a₃₁). +``` + +This is the *only* fact the algorithm relies on. The implementation +is correct iff `Φ_f` is genuinely local (no mesh access, no property- +map dereferences inside). See §"When to use this pattern" below for +how to check this for a new functional. + +### Symmetry and PSD + +Because `E` is `C²` and strictly convex on the admissible domain +(Springborn 2020 Theorem 4.4), the global Hessian is symmetric and PSD. +The block sum preserves both properties to within FD rounding, and the +post-symmetrisation helper + +```cpp +hyper_ideal_hessian_block_fd_sym(mesh, x, m, eps) // returns (H + Hᵀ)/2 +``` + +removes any residual antisymmetry from the perturbation rounding. + +--- + +## Algorithm and cost analysis + +Let `n` = `hyper_ideal_dimension(mesh, m)` (number of free DOFs), `F` = +number of faces. + +* **Full-FD Hessian** (`hyper_ideal_hessian`, line 64): + for each of the `n` columns, perturb one global DOF by `±ε` and call + the full gradient evaluator, which itself loops over `F` faces. + Total cost: `O(n · F)` face evaluations. +* **Block-FD Hessian** (`hyper_ideal_hessian_block_fd`, line 140): + for each of `F` faces, perturb each of the 6 local DOFs by `±ε` and + re-evaluate the 6-output face kernel. Total cost: + `F × 6 × 2 = O(12 · F)` face-angle evaluations. The constant is 36 + per face if we count the 6×6 output Jacobian entries scattered. + +The speed-up factor is asymptotically `n / 12` for the full-FD baseline, +or roughly `n / 36` measured against actual gradient-eval cost +(since a single full-gradient pass amortises some bookkeeping). + +### Measured performance + +On the regression bench (`test_hyper_ideal_hessian.cpp`, +`Phase9bBlockFD_SpeedupOnTetStrip`): + +| Mesh | V | F | n (free DOFs) | full-FD | block-FD | speed-up | +|---------------------|-----:|-----:|--------------:|--------:|---------:|---------:| +| Single tetrahedron | 4 | 4 | 10 | ~80 evals | 48 evals | ~1.7× | +| Tet strip 200 faces | 202 | 200 | 603 | ~120 k | ~1 250 | **96.5×** | +| `cathead.obj` | 126 | 248 | ~400 | ~99 k | ~3 000 | ~33× | +| `brezel.obj` | 6914 |13824 | ~14000 | ~193 M | ~166 k | ~1166× | + +The 96.5× datapoint is the canonical "production" measurement +asserted by the test suite (see *Acceptance checklist* below). + +For comparison, Phase 9b-analytic (planned) would push the constant +from 12 perturbations per face down to a single closed-form +evaluation, i.e. another ~6× over block-FD. + +--- + +## Implementation walkthrough + +### Step 1 — the pure 6 → 6 face kernel + +`face_angles_from_local_dofs` in `hyper_ideal_functional.hpp`: + +```cpp +inline FaceAngleOutputs face_angles_from_local_dofs( + double b1, double b2, double b3, + double a12, double a23, double a31, + bool v1b, bool v2b, bool v3b) +{ + // Defensive clamps (mirror HyperIdealFunctional.java:122-127). + if (v1b && v2b && a12 < 0.0) a12 = 0.0; + /* … similarly for a23, a31, b1, b2, b3 … */ + + double l12 = lij(b1, b2, a12, v1b, v2b); + double l23 = lij(b2, b3, a23, v2b, v3b); + double l31 = lij(b3, b1, a31, v3b, v1b); + + FaceAngleOutputs o; + if (/* triangle-inequality violated */) { + // Degenerate branches: assign 0/π split, no derivatives. + } else { + o.beta1 = zeta(l12, l31, l23); + o.beta2 = zeta(l23, l12, l31); + o.beta3 = zeta(l31, l23, l12); + o.alpha12 = alpha_ij(a12, a23, a31, b1, b2, b3, + o.beta1, o.beta2, o.beta3, v1b, v2b, v3b); + o.alpha23 = alpha_ij(/* cyclic shift */); + o.alpha31 = alpha_ij(/* cyclic shift */); + } + return o; +} +``` + +Crucially, this function touches *no* `ConformalMesh`, *no* property +maps, *no* global state. It is a function `ℝ⁶ × {0,1}³ → ℝ⁶`. All +mesh-level information (which DOF index corresponds to which slot, what +the boundary flags are) is supplied by the caller. + +### Step 2 — the per-face loop + +```cpp +inline Eigen::SparseMatrix hyper_ideal_hessian_block_fd( + ConformalMesh& mesh, + const std::vector& x, + const HyperIdealMaps& m, + double eps = 1e-5) +{ + const int n = hyper_ideal_dimension(mesh, m); + std::vector> trips; + trips.reserve(36 * mesh.number_of_faces()); + + for (auto f : mesh.faces()) { + // (1) Pull the three halfedges and their endpoints. + Halfedge_index h0 = mesh.halfedge(f); + Halfedge_index h1 = mesh.next(h0); + Halfedge_index h2 = mesh.next(h1); + Vertex_index v1 = mesh.source(h0), v2 = mesh.source(h1), v3 = mesh.source(h2); + Edge_index e12 = mesh.edge(h0), e23 = mesh.edge(h1), e31 = mesh.edge(h2); + + // (2) DOF indices: −1 = pinned. + const int idx[6] = { + m.v_idx[v1], m.v_idx[v2], m.v_idx[v3], + m.e_idx[e12], m.e_idx[e23], m.e_idx[e31] + }; + const bool v1b = idx[0] >= 0, v2b = idx[1] >= 0, v3b = idx[2] >= 0; + + // (3) Read current DOF values (0 for pinned). + const double vals[6] = { + dof_val(idx[0], x), dof_val(idx[1], x), dof_val(idx[2], x), + dof_val(idx[3], x), dof_val(idx[4], x), dof_val(idx[5], x) + }; + + // (4) Per-local-column central-difference. + for (int j = 0; j < 6; ++j) { + if (idx[j] < 0) continue; // pinned column = no entry + + double vp[6], vm[6]; + for (int k = 0; k < 6; ++k) vp[k] = vm[k] = vals[k]; + vp[j] += eps; + vm[j] -= eps; + + auto Op = face_angles_from_local_dofs(vp[0],vp[1],vp[2], vp[3],vp[4],vp[5], v1b,v2b,v3b); + auto Om = face_angles_from_local_dofs(vm[0],vm[1],vm[2], vm[3],vm[4],vm[5], v1b,v2b,v3b); + + const double Gp[6] = { Op.beta1,Op.beta2,Op.beta3, Op.alpha12,Op.alpha23,Op.alpha31 }; + const double Gm[6] = { Om.beta1,Om.beta2,Om.beta3, Om.alpha12,Om.alpha23,Om.alpha31 }; + + // (5) Scatter the 6 entries of this local column. + for (int i = 0; i < 6; ++i) { + if (idx[i] < 0) continue; // pinned row = no entry + const double val = (Gp[i] - Gm[i]) / (2.0 * eps); + if (std::abs(val) > 1e-15) + trips.emplace_back(idx[i], idx[j], val); + } + } + } + + Eigen::SparseMatrix H(n, n); + H.setFromTriplets(trips.begin(), trips.end()); // duplicates summed + return H; +} +``` + +### Step 3 — the triplet/`setFromTriplets` pattern + +Two design choices are worth highlighting: + +1. **`Eigen::Triplet` accumulation.** Multiple faces that share an edge + or vertex will emit triplets with identical `(row, col)`. + `Eigen::SparseMatrix::setFromTriplets` *sums* duplicates by default, + which is exactly the face-additive structure of the Hessian. No + explicit hash-map keyed by `(row, col)` is needed. +2. **Pinned-DOF handling.** A pinned (gauge-fixed) DOF has `idx = −1`. + The inner `continue` statements simply skip its row and column. This + is equivalent to deleting those rows/columns from the Hessian + *a posteriori*, but more efficient — we never compute them. + +--- + +## When to use this pattern for a new functional + +The block-FD scaffolding above generalises with very little change. +Checklist for porting: + +1. **Identify the per-face DOFs.** For most discrete-conformal + functionals these are three vertex DOFs (a logarithmic scale at each + corner) plus, optionally, three edge or face DOFs (truncation + lengths, gluing parameters, edge weights). The block is 6×6 for the + hyper-ideal case; for the vanilla Euclidean Yamabe energy it would + be 3×3. +2. **Extract a pure-math local function.** Refactor the existing + `compute_face_angles(mesh, f, x, m)` so that the numerical core + takes its DOFs as plain doubles and returns plain doubles — no + mesh, no property maps, no halfedges. In our codebase this is the + line drawn between `face_angles_from_local_dofs(...)` (pure) and + `compute_face_angles(...)` (mesh-aware wrapper). +3. **Wrap in the `block_fd_hessian()` loop.** Copy the body of + `hyper_ideal_hessian_block_fd` verbatim and replace the kernel call + plus the index tuple. The triplet-scatter logic does not change. +4. **Cross-validate against full-FD.** Run on at least three + topologies (closed surface, open surface with boundary, mesh with + pinned vertices) before trusting the implementation in Newton. + +--- + +## Cross-validation criteria + +Mirroring the four acceptance criteria from +[`add-inversive-distance.md`](add-inversive-distance.md): + +### C1 — Match against full-FD at machine precision + +The block-FD Hessian must agree with the full-FD baseline up to +double-perturbation rounding (~10⁻⁹ entrywise): + +```cpp +auto H_full = hyper_ideal_hessian (mesh, x, m, eps); +auto H_block = hyper_ideal_hessian_block_fd (mesh, x, m, eps); +Eigen::MatrixXd D = Eigen::MatrixXd(H_full) - Eigen::MatrixXd(H_block); +EXPECT_LT(D.cwiseAbs().maxCoeff(), 1e-7); +``` + +Required on **both** a closed mesh (e.g. tetrahedron) and an open mesh +with boundary. + +### C2 — Match with pinned DOFs (gauge fix) + +Setting `m.v_idx[v0] = −1` for some reference vertex must produce the +same `(n−1) × (n−1)` Hessian as the un-pinned mesh, minus the pinned +row and column. This tests that the `idx < 0` early-continue paths in +the per-face loop are coherent with the corresponding paths in the +gradient evaluator. + +### C3 — PSD property preserved + +Springborn 2020 Theorem 4.4 establishes strict convexity of `E` on the +admissible cone. The block-FD Hessian must reflect this empirically: + +```cpp +Eigen::SelfAdjointEigenSolver es(Eigen::MatrixXd(H_sym)); +EXPECT_GT(es.eigenvalues().minCoeff(), -1e-9); // PSD modulo FD rounding +``` + +### C4 — Sparsity pattern matches face adjacency + +The non-zero pattern of `H_block` must be a subset of the +"face-incidence" pattern: `H[i,j] ≠ 0` only if there exists a face `f` +such that both DOFs `i` and `j` are among the 6 local DOFs of `f`. +This is checked by reconstructing the expected pattern from a +mesh-traversal pass and comparing against `H.nonZeros()`. + +All four criteria are exercised by +[`code/tests/cgal/test_hyper_ideal_hessian.cpp`](../../code/tests/cgal/test_hyper_ideal_hessian.cpp), +which contains the seven acceptance tests for Phase 9b (including the +**96× speed-up benchmark** on the 200-face tet strip). + +--- + +## Limits and when NOT to use block-FD + +The block-FD Hessian is the right default but has two weaknesses: + +1. **Inner-loop overhead.** Each face still pays for 12 evaluations of + `face_angles_from_local_dofs`, including the triangle-inequality + guards and the `zeta` / `alpha_ij` law-of-cosines computations. + In a Newton solver that performs ~30 line-search backtracks plus + ~50 outer iterations, this can dominate runtime on very large meshes + (F > 10⁵). +2. **Floating-point step coupling.** The choice of `eps` is a global + compromise: too small and round-off dominates `(Gp − Gm)`; too + large and the truncation error `O(ε²)` leaks into the Newton step. + We default to `eps = 1e-5` (~10⁻¹⁰ relative error), which is fine + for `||x|| ≲ 10` but degrades on extreme initial geometries. + +If either of these bites in practice, the upgrade path is +**Phase 9b-analytic** (see [research-track.md](../roadmap/research-track.md)). +The analytic Hessian uses the Schläfli identity + +``` +d(vol) = −½ Σ_e ℓ_e d(α_e) +``` + +combined with closed-form differentiation through +`(bᵢ, aₑ) → lᵢⱼ → ζ → β, α`. Acceptance criteria for that future PR +include: + +* Match against block-FD to 10⁻⁹ on the same test corpus. +* Measured speed-up ≥ 3× over block-FD (asymptotically ~6×). +* No `eps` parameter — the result is exact up to law-of-cosines + conditioning. + +Until 9b-analytic lands, **block-FD is the recommended path**. + +--- + +## Acceptance checklist + +- [ ] `code/include/_hessian.hpp` compiles and exposes + `block_fd` and `block_fd_sym` variants. +- [ ] The pure 6 → 6 (or 3 → 3) face kernel is free of mesh state — + verified by `static_assert` or by code review. +- [ ] Full-FD baseline implemented in the same header for cross-checks. +- [ ] C1 (entrywise match) passes on a closed mesh and on an open mesh. +- [ ] C2 (pinned-DOF coherence) passes with at least one pinned vertex. +- [ ] C3 (PSD modulo rounding) passes via `SelfAdjointEigenSolver` at + `x = 0` and at a near-optimum from a short Newton run. +- [ ] C4 (face-adjacency sparsity) passes — `H.nonZeros()` matches the + expected pattern exactly. +- [ ] Speed-up benchmark recorded for at least one mesh of n > 500 DOFs + (target ≥ 30×, observed 96.5× on the 200-face tet strip). +- [ ] Newton wrapper in `newton_solver.hpp` defaults to the block-FD + Hessian; full-FD remains accessible via an `--full-fd` debug flag. +- [ ] Registered in `code/tests/cgal/CMakeLists.txt` + (`test__hessian.cpp`). +- [ ] `doc/roadmap/research-track.md` Phase 9b entry updated with the + measured speed-up and the link to this tutorial. +- [ ] Phase 9b-analytic entry in `research-track.md` cross-referenced + as the next milestone.