docs: restructure documentation into focused files

README.md: reduced from 703 to ~75 lines — what/why, status, quick start, minimal usage example, navigation table to doc/ files. doc/architecture/overall_pipeline.md: trimmed — roadmap, extension points, declarative pipeline YAML, and references sections removed (each now has its own dedicated file). Replaced with a link table. New files: doc/getting-started.md — build modes, single-test invocation, CLI doc/api/pipeline.md — full pipeline API with code for all 3 geometries doc/api/extending.md — new functionals, geometry modes, Java porting guide doc/api/contracts.md — processing unit preconditions/provides table doc/api/cgal-package.md — Phase 8 CGAL package design + YAML pipeline (TODO) doc/math/geometry-modes.md — Euclidean/Spherical/HyperIdeal comparison doc/math/references.md — all papers by module doc/roadmap/phases.md — Phases 1–10 with porting/research boundary doc/roadmap/java-parity.md — Java vs C++ feature parity table doc/contributing.md — language policy, test standards, release flow Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
2026-05-17 21:17:15 +02:00
parent 279f964b96
commit 14134b99ce
12 changed files with 1229 additions and 835 deletions
--- a/doc/api/cgal-package.md
+++ b/doc/api/cgal-package.md
@@ -0,0 +1,163 @@
+# Phase 8 — CGAL Package Design
+
+> **Status: planned.** This document describes the target architecture for Phase 8.
+> No code has been written yet. The design is informed by the CGAL package submission
+> guidelines at https://www.cgal.org/developers.html
+
+---
+
+## Goal
+
+Integrate conformallab++ into the CGAL library as a proper CGAL package:
+`Discrete_conformal_map`. The package must satisfy all CGAL submission requirements:
+traits-class design, Doxygen documentation, CGAL-format test suite, and coverage of
+the CGAL coding conventions.
+
+---
+
+## 8a — Traits class & concepts
+
+The current code is tightly coupled to `CGAL::Surface_mesh<Point3>`. Phase 8a introduces
+a traits class that separates the mesh type from the algorithm:
+
+```cpp
+// TODO(Phase 8a): implement this header
+// include/CGAL/Conformal_map_traits.h
+
+template<
+    typename MeshType,      // any CGAL halfedge mesh
+    typename KernelType,    // CGAL kernel
+    typename ScalarType = double
+>
+struct Conformal_map_traits {
+    using Mesh   = MeshType;
+    using Kernel = KernelType;
+    using FT     = ScalarType;
+    // ... vertex/edge/face descriptor types
+    // ... property map access
+};
+```
+
+Concept checks will ensure any user-provided mesh satisfies the halfedge mesh concept.
+
+---
+
+## 8b — Public header hierarchy
+
+A clean public API separate from the internal implementation:
+
+```
+include/CGAL/
+  Discrete_conformal_map.h        ← single user-facing include
+  Conformal_map_traits.h
+  Conformal_newton_solver.h
+  Conformal_layout.h
+  Conformal_cut_graph.h
+  conformal_map_package.h         ← PackageDescription
+```
+
+All existing `include/*.hpp` headers remain as internal implementation details,
+not part of the public CGAL API.
+
+---
+
+## 8c — CGAL-style documentation
+
+```
+doc/Conformal_map/
+  PackageDescription.txt
+  User_manual.md
+  Reference_manual.md
+  fig/                            ← pipeline diagrams, mathematical figures
+```
+
+All public functions and concepts require Doxygen comments following the CGAL style.
+
+---
+
+## 8d — CGAL test format
+
+```
+test/Conformal_map/
+  CMakeLists.txt                  ← CGAL-format, uses find_package(CGAL)
+  test_euclidean_functional.cpp
+  test_newton_solver.cpp
+  ...
+```
+
+The existing GTest suite remains. CGAL-format tests are added alongside as a separate
+target, following the CGAL test infrastructure conventions.
+
+---
+
+## 8e — Declarative YAML pipeline
+
+A lightweight YAML format for reproducible experiments. The CLI accepts
+`--pipeline experiment.yml`; the validator checks `require`/`provide` tokens
+before execution.
+
+Example:
+
+```yaml
+pipeline:
+  name: flat_torus_period
+  geometry: euclidean
+
+  input:
+    source: data/torus.off
+
+  steps:
+    - id: setup
+      unit: setup_euclidean_maps
+      provide: [maps_initialised]
+
+    - id: gauss_bonnet
+      unit: enforce_gauss_bonnet
+      require: [maps_initialised]
+      provide: [gauss_bonnet_satisfied]
+
+    - id: solve
+      unit: newton_euclidean
+      require: [gauss_bonnet_satisfied]
+      params:
+        tol: 1.0e-10
+        max_iter: 200
+      provide: [x_converged]
+
+    - id: cut
+      unit: compute_cut_graph
+      require: [mesh_closed]
+      provide: [cut_graph]
+
+    - id: layout
+      unit: euclidean_layout
+      require: [x_converged, cut_graph]
+      params:
+        normalise: true
+      provide: [layout_uv, holonomy]
+
+    - id: period
+      unit: compute_period_matrix
+      require: [holonomy]
+      provide: [tau]
+
+  output:
+    layout: out/torus_layout.off
+    json:   out/torus_result.json
+    tau:    out/torus_tau.txt
+```
+
+The contract table in [contracts.md](contracts.md) defines the valid `require`/`provide`
+token vocabulary.
+
+---
+
+## TODO
+
+- [ ] Design `Conformal_map_traits.h` interface (8a)
+- [ ] Define concept requirements for `MeshType` (8a)
+- [ ] Create `include/CGAL/` header skeleton (8b)
+- [ ] Write `PackageDescription.txt` (8c)
+- [ ] Port GTest tests to CGAL format (8d)
+- [ ] Implement YAML validator (8e)
+- [ ] CLI: `--pipeline` flag (8e)
--- a/doc/api/contracts.md
+++ b/doc/api/contracts.md
@@ -0,0 +1,69 @@
+# Processing Unit Contracts
+
+Each pipeline unit has explicit preconditions and guarantees.
+A pipeline is valid if every unit's preconditions are satisfied by the outputs
+of all preceding units.
+
+---
+
+## Contract table
+
+| Unit | Preconditions | Provides |
+|---|---|---|
+| `load_mesh()` | Valid file path, supported format (OFF/OBJ/PLY) | Manifold, oriented, triangulated `ConformalMesh` |
+| `setup_*_maps()` | Triangulated mesh | Initialised property maps; `lambda0` zeroed |
+| `compute_*_lambda0_from_mesh()` | `setup_*_maps()` called | `lambda0[e]` set from 3-D edge lengths |
+| `check_gauss_bonnet()` | `theta_v[v]` set for all vertices | Throws `std::runtime_error` if Σ(2π−Θᵥ) ≠ 2π·χ(M) |
+| `enforce_gauss_bonnet()` | `theta_v[v]` set | Σ(2π−Θᵥ) = 2π·χ(M) guaranteed; `theta_v` modified |
+| `newton_euclidean()` | GB satisfied · DOFs assigned · `lambda0` initialised | `NewtonResult.x` — converged scale factors; `.converged`, `.iterations`, `.grad_inf_norm` |
+| `newton_spherical()` | GB satisfied · DOFs assigned · gauge vertex pinned | Same as above |
+| `newton_hyper_ideal()` | `assign_all_dof_indices()` called · `lambda0` initialised | Same as above |
+| `compute_cut_graph()` | Closed, orientable, triangulated mesh | `CutGraph.cut_edge_flags` — 2g seam edges; `.genus` |
+| `euclidean_layout()` | `newton_euclidean()` converged | `Layout2D.uv[v]` · `.halfedge_uv[h]` · `HolonomyData.translations` |
+| `spherical_layout()` | `newton_spherical()` converged | `Layout3D.xyz[v]` |
+| `hyper_ideal_layout()` | `newton_hyper_ideal()` converged | `Layout2D.uv[v]` · `.halfedge_uv[h]` · `HolonomyData.mobius_maps` |
+| `normalise_euclidean()` | `layout.success == true` | Centroid at origin, major axis = x-axis |
+| `normalise_hyperbolic()` | `layout.success == true` | Weighted centroid at Poincaré disk origin |
+| `normalise_spherical()` | `layout.success == true` | Area centroid at north pole |
+| `compute_period_matrix()` | `hol.translations.size() >= 2` | `PeriodData.tau ∈ ℍ`; optionally SL(2,ℤ)-reduced |
+| `compute_fundamental_domain()` | `HolonomyData` (genus 1) | CCW parallelogram `{0, ω₁, ω₁+ω₂, ω₂}` + edge identifications |
+| `tiling_neighbourhood()` | `Layout2D` + `HolonomyData` (genus 1) | `(2m+1)·(2n+1)` translated layout copies |
+| `save_layout_off()` | `Layout2D` + mesh | OFF file with UV coordinates |
+| `save_result_json/xml()` | `NewtonResult` + optional `Layout2D` | JSON/XML serialised result |
+
+---
+
+## DOF assignment conventions
+
+```cpp
+// Euclidean / Spherical: pin one vertex, assign sequential indices
+maps.v_idx[first_vertex] = -1;          // pinned: u_v = 0
+int idx = 0;
+for (auto v : remaining_vertices)
+    maps.v_idx[v] = idx++;
+
+// HyperIdeal: all vertices and edges are free DOFs
+int n_dofs = assign_all_dof_indices(mesh, maps);
+// maps.v_idx[v] ∈ [0, n_v)
+// maps.e_idx[e] ∈ [n_v, n_v + n_e)
+```
+
+`-1` in `v_idx` or `e_idx` means the DOF is pinned (fixed at its `lambda0` value).
+The solver never writes to pinned DOFs. All indexing is 0-based and contiguous.
+
+---
+
+## Gauss–Bonnet — the most common source of failure
+
+Prescribing angles that violate Gauss–Bonnet means no conformal factor can realise the
+target metric — the Newton solver will iterate indefinitely without converging.
+
+```cpp
+// Option A: verify before solving (throws on violation)
+check_gauss_bonnet(mesh, maps);
+
+// Option B: auto-correct (redistributes defect uniformly across all vertices)
+enforce_gauss_bonnet(mesh, maps);
+```
+
+The target violation is `|Σ(2π−Θᵥ) − 2π·χ(M)| > tol`. Default tolerance: `1e-10`.
--- a/doc/api/extending.md
+++ b/doc/api/extending.md
@@ -0,0 +1,141 @@
+# Extending conformallab++
+
+## Adding a new functional
+
+All three existing functionals (`euclidean_functional.hpp`, `spherical_functional.hpp`,
+`hyper_ideal_functional.hpp`) follow the same pattern. Copy the structure of the
+simplest one (`euclidean_functional.hpp`) as a template.
+
+### Step 1 — Maps struct
+
+```cpp
+// my_functional.hpp
+struct MyMaps {
+    ConformalMesh::Property_map<Vertex_index, double> lambda0;  // initial log-lengths
+    ConformalMesh::Property_map<Vertex_index, double> theta_v;  // target angles
+    ConformalMesh::Property_map<Vertex_index, int>    v_idx;    // DOF index, -1 = pinned
+    // add edge DOFs if needed:
+    ConformalMesh::Property_map<Edge_index,   int>    e_idx;
+};
+
+MyMaps setup_my_maps(ConformalMesh& mesh);
+void compute_my_lambda0_from_mesh(ConformalMesh& mesh, MyMaps& maps);
+```
+
+### Step 2 — Gradient
+
+The gradient must satisfy: `G_v = Σ(angle contributions at v) − theta_v[v]`.
+
+```cpp
+std::vector<double> my_gradient(
+    const ConformalMesh& mesh,
+    const std::vector<double>& x,
+    const MyMaps& maps)
+{
+    std::vector<double> G(n_dofs, 0.0);
+    for (auto f : mesh.faces()) {
+        // compute angles in face f given current x
+        // accumulate into G[maps.v_idx[v]] for each vertex v of f
+    }
+    // subtract target angles
+    for (auto v : mesh.vertices())
+        if (maps.v_idx[v] >= 0)
+            G[maps.v_idx[v]] -= maps.theta_v[v];
+    return G;
+}
+```
+
+### Step 3 — Gradient check test
+
+Before claiming correctness, verify with finite differences. Copy any `GradientCheck_*`
+test suite from `tests/cgal/test_*_functional.cpp`:
+
+```cpp
+TEST(MyFunctional, GradientCheck_Triangle) {
+    ConformalMesh mesh = make_single_triangle();
+    MyMaps maps = setup_my_maps(mesh);
+    // ... assign DOFs, set natural theta ...
+
+    const double eps = 1e-5;
+    auto G = my_gradient(mesh, x0, maps);
+    for (int i = 0; i < n; ++i) {
+        std::vector<double> xp = x0, xm = x0;
+        xp[i] += eps;  xm[i] -= eps;
+        double fd = (my_energy(mesh, xp, maps) - my_energy(mesh, xm, maps)) / (2*eps);
+        EXPECT_NEAR(G[i], fd, 1e-7) << "DOF " << i;
+    }
+}
+```
+
+### Step 4 — Hessian and Newton wrapper
+
+```cpp
+// Reuse solve_linear_system from newton_solver.hpp
+Eigen::SparseMatrix<double> H = my_hessian(mesh, x, maps);
+bool used_fallback;
+auto dx = conformallab::solve_linear_system(H, -G, &used_fallback);
+```
+
+Or write a thin `newton_my()` wrapper following the structure of `newton_euclidean()`.
+
+---
+
+## Adding a new geometry mode
+
+To add a new space (e.g. de Sitter, flat 3-torus):
+
+1. **Trilateration function** — implement `trilaterate_my(p1, p2, l12, l13, l23)` that
+   places a third point given two placed points and three edge lengths.
+   This is the only geometry-specific part of the layout.
+
+2. **BFS structure** — reuse the priority-BFS loop from `euclidean_layout()` in `layout.hpp`.
+   The loop itself is geometry-agnostic; swap in your trilateration function.
+
+3. **Holonomy** — if the space has a holonomy group, record the transition maps at seam edges
+   the same way `HolonomyData` does for translations (Euclidean) or Möbius maps (hyperbolic).
+
+---
+
+## Adding a new processing unit
+
+1. Declare preconditions and outputs explicitly (see [contracts.md](contracts.md)).
+
+2. Write a header-only implementation in `code/include/my_unit.hpp`.
+
+3. Add tests in `code/tests/cgal/test_my_unit.cpp` and register in
+   `code/tests/cgal/CMakeLists.txt`.
+
+4. Pipeline validation is currently manual (documented contracts). A compile-time or
+   runtime check via the declarative YAML pipeline (Phase 8e) is the planned upgrade —
+   see [cgal-package.md](cgal-package.md).
+
+---
+
+## Porting from Java
+
+When porting a class from the Java library:
+
+1. Locate the original at [github.com/varylab/conformallab](https://github.com/varylab/conformallab)
+   under `src/main/java/de/varylab/discreteconformal/`.
+
+2. The Java library uses `CoHDS` (half-edge data structure) with intrusive vertex/edge data.
+   Map these to CGAL property maps:
+   ```
+   Java: vertex.getLambda()     → C++: maps.lambda0[v]
+   Java: edge.getAlpha()        → C++: maps.e_alpha[e]
+   Java: vertex.getTheta()      → C++: maps.theta_v[v]
+   Java: vertex.getSolverIndex() → C++: maps.v_idx[v]
+   ```
+
+3. Java uses `HalfedgeInterface` adapters with named accessors like `getOppositeVertex()`.
+   C++ equivalent:
+   ```cpp
+   // Java: h.getOppositeVertex()
+   Vertex_index v_opp = mesh.target(mesh.next(h));
+
+   // Java: h.getNextHalfedge().getOppositeVertex()
+   Vertex_index v_opp2 = mesh.target(mesh.next(mesh.next(h)));
+   ```
+
+4. Port the test cases from the Java `@Test` methods. The "natural theta" trick
+   (`theta_v[v] = actual_angle_sum`) works the same way in both languages.
--- a/doc/api/pipeline.md
+++ b/doc/api/pipeline.md
@@ -0,0 +1,214 @@
+# Pipeline API Reference
+
+The full pipeline from mesh loading to layout and serialisation.
+See [doc/architecture/overall_pipeline.md](../architecture/overall_pipeline.md) for the
+Mermaid diagram and stage descriptions.
+
+---
+
+## Minimal Euclidean pipeline
+
+```cpp
+#include "conformal_mesh.hpp"
+#include "mesh_io.hpp"
+#include "euclidean_functional.hpp"
+#include "gauss_bonnet.hpp"
+#include "newton_solver.hpp"
+#include "layout.hpp"
+
+using namespace conformallab;
+
+ConformalMesh mesh = load_mesh("input.off");
+
+// 1. Set up property maps and initialise log-edge-lengths from 3-D positions
+EuclideanMaps maps = setup_euclidean_maps(mesh);
+compute_euclidean_lambda0_from_mesh(mesh, maps);
+
+// 2. Assign DOF indices — pin first vertex (gauge fix for open meshes)
+auto vit = mesh.vertices().begin();
+maps.v_idx[*vit++] = -1;           // pinned: u_v = 0
+int idx = 0;
+for (; vit != mesh.vertices().end(); ++vit)
+    maps.v_idx[*vit] = idx++;
+
+// 3. Set target angles — "natural equilibrium" makes x* = 0
+std::vector<double> x0(idx, 0.0);
+auto G0 = euclidean_gradient(mesh, x0, maps);
+for (auto v : mesh.vertices())
+    if (maps.v_idx[v] >= 0)
+        maps.theta_v[v] -= G0[maps.v_idx[v]];
+
+// 4. Check Gauss–Bonnet (mandatory before solving)
+check_gauss_bonnet(mesh, maps);     // throws if Σ(2π−Θᵥ) ≠ 2π·χ(M)
+
+// 5. Solve
+NewtonResult res = newton_euclidean(mesh, x0, maps);
+// res.converged, res.x, res.iterations, res.grad_inf_norm
+
+// 6. Layout
+Layout2D layout = euclidean_layout(mesh, res.x, maps);
+// layout.uv[v.idx()], layout.halfedge_uv[h.idx()]
+```
+
+---
+
+## Layout with holonomy (closed surfaces)
+
+```cpp
+#include "cut_graph.hpp"
+#include "period_matrix.hpp"
+#include "fundamental_domain.hpp"
+
+// Tree-cotree cut graph — 2g seam edges
+CutGraph cg = compute_cut_graph(mesh);
+// cg.cut_edge_flags[e.idx()] — true if seam edge
+// cg.genus                   — topological genus g
+
+// Layout with holonomy tracking
+HolonomyData hol;
+Layout2D layout = euclidean_layout(mesh, res.x, maps, &cg, &hol, /*normalise=*/true);
+// hol.translations[i] — lattice generator ωᵢ ∈ ℂ (Euclidean)
+
+// Period matrix (genus 1 flat torus)
+PeriodData pd = compute_period_matrix(hol);
+// pd.tau               — complex period ratio τ = ω₂/ω₁ ∈ ℍ
+// pd.omega             — {ω₁, ω₂} as complex<double>
+// pd.in_fundamental_domain — true if SL(2,ℤ)-reduced to |τ|≥1, |Re(τ)|<½
+
+// Fundamental domain parallelogram
+FundamentalDomain fd = compute_fundamental_domain(hol);
+// fd.vertices[0..3]    — CCW corners: {0, ω₁, ω₁+ω₂, ω₂}
+// fd.generators[0..1]  — ω₁, ω₂
+
+// Tiling copies for universal cover visualisation
+auto tiles = tiling_neighbourhood(layout, hol, /*m_max=*/2, /*n_max=*/2);
+// returns (2·m_max+1)·(2·n_max+1) translated layout copies
+```
+
+---
+
+## HyperIdeal pipeline
+
+```cpp
+#include "hyper_ideal_functional.hpp"
+#include "newton_solver.hpp"
+#include "layout.hpp"
+#include "cut_graph.hpp"
+
+HyperIdealMaps maps = setup_hyper_ideal_maps(mesh);
+compute_hyper_ideal_lambda0_from_mesh(mesh, maps);
+
+// HyperIdeal has both vertex and edge DOFs — assign all automatically
+int n_dofs = assign_all_dof_indices(mesh, maps);
+// No vertex needs to be pinned — strictly convex functional
+
+std::vector<double> x0(n_dofs, 0.0);
+NewtonResult res = newton_hyper_ideal(mesh, x0, maps);
+
+// Layout with Möbius holonomy
+CutGraph cg = compute_cut_graph(mesh);
+HolonomyData hol;
+Layout2D layout = hyper_ideal_layout(mesh, res.x, maps, &cg, &hol);
+// hol.mobius_maps[i] — Möbius isometry Tᵢ ∈ SU(1,1) per cut edge
+```
+
+---
+
+## Spherical pipeline
+
+```cpp
+#include "spherical_functional.hpp"
+#include "newton_solver.hpp"
+#include "layout.hpp"
+
+SphericalMaps maps = setup_spherical_maps(mesh);
+compute_spherical_lambda0_from_mesh(mesh, maps);
+
+// Pin one vertex (gauge fix) + assign DOFs
+auto vit = mesh.vertices().begin();
+maps.v_idx[*vit++] = -1;
+int idx = 0;
+for (; vit != mesh.vertices().end(); ++vit)
+    maps.v_idx[*vit] = idx++;
+
+std::vector<double> x0(idx, 0.0);
+NewtonResult res = newton_spherical(mesh, x0, maps);
+
+// 3-D layout on the unit sphere
+Layout3D slayout = spherical_layout(mesh, res.x, maps);
+// slayout.xyz[v.idx()] — 3-D position on S²
+```
+
+---
+
+## SparseQR fallback — direct use
+
+```cpp
+#include "newton_solver.hpp"
+
+bool used_fallback = false;
+Eigen::VectorXd dx = conformallab::solve_linear_system(H, rhs, &used_fallback);
+if (used_fallback)
+    std::cerr << "Warning: H is rank-deficient, used SparseQR\n";
+```
+
+The fallback is automatically invoked by all three `newton_*` functions when
+`SimplicialLDLT` fails. It finds the minimum-norm Newton step orthogonal to the null space.
+
+---
+
+## Serialisation
+
+```cpp
+#include "serialization.hpp"
+#include "mesh_io.hpp"
+
+// Save layout as OFF mesh file
+save_layout_off("layout.off", mesh, layout);
+
+// Full result: DOF vector + metadata + layout UVs
+save_result_json("result.json", res, "euclidean", mesh, &layout);
+save_result_xml ("result.xml",  res, "euclidean", mesh, &layout);
+
+// Round-trip load
+NewtonResult res2;
+std::string geom;
+Layout2D uv2;
+load_result_json("result.json", &res2, &geom, &uv2);
+```
+
+---
+
+## Attaching custom property maps
+
+```cpp
+// Add a custom per-vertex map
+auto [curv, created] = mesh.add_property_map<Vertex_index, double>("v:my_curv", 0.0);
+curv[v] = 3.14;
+
+// Retrieve an existing map
+auto [curv2, found] = mesh.property_map<Vertex_index, double>("v:my_curv");
+```
+
+---
+
+## Halfedge traversal conventions
+
+```cpp
+for (auto f : mesh.faces()) {
+    auto h0 = mesh.halfedge(f);   // canonical halfedge of face f
+    auto h1 = mesh.next(h0);
+    auto h2 = mesh.next(h1);
+
+    Vertex_index v1 = mesh.source(h0);  // = mesh.target(h2)
+    Vertex_index v2 = mesh.source(h1);
+    Vertex_index v3 = mesh.source(h2);
+
+    // Angle at v3 is opposite to halfedge h0 (edge v1–v2)
+    // h_alpha[h0] = α₃,  h_alpha[h1] = α₁,  h_alpha[h2] = α₂
+
+    Edge_index e0 = mesh.edge(h0);
+    bool is_seam = cg.cut_edge_flags[e0.idx()];
+    bool is_boundary = mesh.is_border(mesh.opposite(h0));
+}
+```