#pragma once
// Copyright (c) 2024-2026 Tarik Moussa.
// SPDX-License-Identifier: MIT

// layout.hpp
//
// Phase 5/6/7 — Layout / embedding: DOF vector → vertex coordinates in the
// target geometry via BFS-trilateration.
//
// ┌──────────────────────────────────────────────────────────────────────────┐
// │  Algorithm (all three geometries)                                        │
// │                                                                          │
// │  1. For every edge e compute the updated length l_e from the DOF x.     │
// │  2. Choose root face = largest 3-D area face (quality heuristic).       │
// │  3. Priority BFS over the dual graph: faces are processed in order of   │
// │     shortest distance (BFS depth) from the root face, minimising error  │
// │     accumulation.  Equal depth: arbitrary tie-break.                    │
// │  4. For each face: trilaterate the one unplaced vertex from the two     │
// │     already-placed vertices on the shared edge.                          │
// │  5. If a CutGraph is supplied, cut edges are treated as boundary;       │
// │     holonomy is recorded for each cut edge.                              │
// │  6. After the first connected component, unvisited faces start new      │
// │     BFS sweeps (multi-component support).                               │
// │                                                                          │
// │  For open meshes: globally consistent placement.                        │
// │  For closed meshes without CutGraph: first (shallowest) visit wins,    │
// │     has_seam = true.                                                    │
// │  For closed meshes with CutGraph: each vertex placed exactly once;     │
// │     holonomy = translation (Euclidean) or Möbius map (hyperbolic).     │
// └──────────────────────────────────────────────────────────────────────────┘
//
// Trilateration
//   Euclidean   — exact (analytic formula in ℝ²)
//   Spherical   — exact (spherical law of cosines on S²)
//   Hyperbolic  — exact (hyperbolic law of cosines + Möbius maps,
//                 Poincaré disk model)
//
// Normalisation
//   normalise_euclidean(layout)   — centroid → origin, major axis → x-axis (PCA)
//   normalise_hyperbolic(layout)  — face-area-weighted iterative Möbius centering
//   normalise_spherical(layout)   — rotate centroid to north pole (Rodrigues)
//
// Holonomy (genus-g surfaces with CutGraph)
//   HolonomyData.translations[i]  — translation ω_i   (Euclidean / spherical)
//   HolonomyData.mobius_maps[i]   — Möbius map T_i    (hyperbolic)
//   For a flat torus: ω_1, ω_2 are the lattice generators; τ = ω_2/ω_1 ∈ ℍ.
//
// Texture atlas
//   Layout2D.halfedge_uv[h.idx()] — UV of source(h) as seen from face(h).
//   At seam edges the two opposite halfedges carry different UV values,
//   enabling proper per-halfedge UV for GPU texture atlasing.
//
// Möbius map
//   MobiusMap::from_three(z1→w1, z2→w2, z3→w3) — fit a Möbius transformation
//   to three point correspondences (via 3×3 complex linear system).

#include "conformal_mesh.hpp"
#include "euclidean_functional.hpp"
#include "spherical_functional.hpp"
#include "hyper_ideal_functional.hpp"
#include "cut_graph.hpp"
#include <Eigen/Dense>
#include <vector>
#include <queue>
#include <complex>
#include <cmath>
#include <algorithm>
#include <limits>
#include <fstream>
#include <string>

namespace conformallab {

// ── Möbius map ────────────────────────────────────────────────────────────────
//
// T(z) = (a·z + b) / (c·z + d)
//
// For hyperbolic holonomy the map is an orientation-preserving isometry of
// the Poincaré disk (SU(1,1) element).
// ─────────────────────────────────────────────────────────────────────────────
/// Möbius transformation `T(z) = (a·z + b) / (c·z + d)` of the Riemann
/// sphere; restricted to SU(1,1) for hyperbolic holonomy on the
/// Poincaré disk.
struct MobiusMap {
    /// Complex scalar used for all entries.
    using C = std::complex<double>;
    C a{1.0, 0.0}; ///< Top-left coefficient.
    C b{0.0, 0.0}; ///< Top-right coefficient.
    C c{0.0, 0.0}; ///< Bottom-left coefficient.
    C d{1.0, 0.0}; ///< Bottom-right coefficient.

    /// Apply the transformation to a complex point.
    C apply(C z) const { return (a * z + b) / (c * z + d); }

    /// Apply the transformation to a 2-D real point (interpreted as `x + iy`).
    Eigen::Vector2d apply(const Eigen::Vector2d& p) const {
        C w = apply(C(p.x(), p.y()));
        return Eigen::Vector2d(w.real(), w.imag());
    }

    /// Identity map.
    static MobiusMap identity() { return {C(1), C(0), C(0), C(1)}; }

    /// Inverse map.
    MobiusMap inverse()               const { return {d, -b, -c, a}; }
    /// Composition `(*this) ∘ T`, i.e. apply `T` first then `*this`.
    MobiusMap compose(const MobiusMap& T) const {
        return { a*T.a + b*T.c,  a*T.b + b*T.d,
                 c*T.a + d*T.c,  c*T.b + d*T.d };
    }

    /// `true` iff the map is the identity up to tolerance `tol`.
    bool is_identity(double tol = 1e-9) const {
        if (std::abs(d) < 1e-14) return false;
        C a_ = a/d, b_ = b/d, c_ = c/d;
        return std::abs(a_ - C(1)) < tol && std::abs(b_) < tol && std::abs(c_) < tol;
    }

    /// Fit T to three point correspondences  T(z_i) = w_i.
    /// Sets d=1 and solves the 3×3 complex linear system.
    /// Returns identity when the system is (near-)singular.
    static MobiusMap from_three(C z1, C w1, C z2, C w2, C z3, C w3) {
        Eigen::Matrix3cd M;
        M << z1, C(1), -w1*z1,
             z2, C(1), -w2*z2,
             z3, C(1), -w3*z3;
        Eigen::Vector3cd rhs; rhs << w1, w2, w3;
        Eigen::ColPivHouseholderQR<Eigen::Matrix3cd> qr(M);
        if (qr.rank() < 3) return identity();
        Eigen::Vector3cd sol = qr.solve(rhs);
        return { sol[0], sol[1], sol[2], C(1.0) };
    }
};

// ── Result types ──────────────────────────────────────────────────────────────

/// Result of a 2-D layout (`euclidean_layout`, `hyper_ideal_layout`):
/// per-vertex UV coordinates plus a per-half-edge UV atlas for seamed
/// textures.
struct Layout2D {
    /// Primary 2-D position per vertex (first / shallowest-BFS-depth visit).
    ///
    /// **Indexing:** `uv[v.idx()]` — indexed by the raw integer vertex index.
    /// **Length:** `mesh.number_of_vertices()`.
    ///
    /// \pre  No vertices have been removed from `mesh` after loading (i.e.
    ///       `mesh.is_valid()` and no compaction was performed).  On a fresh
    ///       `Surface_mesh` loaded from file, `v.idx()` is always in
    ///       `[0, number_of_vertices())` and contiguous.  If vertices were
    ///       deleted and `mesh.collect_garbage()` was called, re-run the
    ///       layout — indices will have shifted.
    ///
    /// Access pattern:
    /// ```cpp
    /// for (auto v : mesh.vertices())
    ///     Eigen::Vector2d uv_v = layout.uv[v.idx()];
    /// ```
    std::vector<Eigen::Vector2d> uv;

    /// UV of `source(h)` as seen from `face(h)`, indexed by `h.idx()`.
    ///
    /// **Indexing:** `halfedge_uv[h.idx()]` — raw integer halfedge index.
    /// **Length:** `mesh.number_of_halfedges()`.  Same no-compaction precondition
    /// as `uv` (see above).
    ///
    /// For interior (non-seam) halfedges: equals `uv[source(h).idx()]`.
    /// For seam halfedges: carries the UV from the virtual unfolding across the
    ///   cut — the trilaterated position that was NOT used as the primary `uv`.
    ///   This gives each face its own copy of a seam vertex, enabling a proper
    ///   per-halfedge GPU texture atlas without vertex duplication.
    ///
    /// Border halfedges (outer face) hold `(0, 0)`.
    std::vector<Eigen::Vector2d> halfedge_uv;

    bool success  = false; ///< `true` iff the BFS placed every vertex.
    bool has_seam = false; ///< `true` when a vertex was reached via two paths.
};

/// Result of a 3-D layout (`spherical_layout`): per-vertex positions on S².
struct Layout3D {
    std::vector<Eigen::Vector3d> pos; ///< Per-vertex spherical positions.
    bool success  = false; ///< `true` iff the BFS placed every vertex.
    bool has_seam = false; ///< `true` when a vertex was reached via two paths.
};

/// Per-cut-edge holonomy.
///
/// Euclidean: translations[i] = ω_i  (translation vector for cut_edge_indices[i]).
///   For a flat torus ω_1, ω_2 are the lattice generators; τ = ω_2/ω_1 ∈ ℍ.
///
/// Hyperbolic: mobius_maps[i] = T_i  (Möbius isometry of the Poincaré disk).
///   T_i(p_actual) = p_virtual  — maps the placed vertex position to the
///   trilaterated virtual position obtained by continuing the unfolding across
///   the cut.
struct HolonomyData {
    std::vector<Eigen::Vector2d> translations;     ///< Euclidean / spherical translation per cut edge.
    std::vector<MobiusMap>       mobius_maps;      ///< Hyperbolic Möbius isometry per cut edge (Phase 7).
    std::vector<std::size_t>     cut_edge_indices; ///< Index (in the cut-graph edge list) of each holonomy entry.

    /// Residual rotation |arg(a)| (radians) of the Euclidean deck isometry
    /// z ↦ a·z + b per cut edge.  Zero for a perfectly flat cone metric;
    /// a non-zero value signals under-convergence or a genuine cone-angle
    /// defect, in which case `translations[i]` (the isometry's translation
    /// part b) is only an approximate lattice generator.  Empty for the
    /// hyperbolic (Möbius) path.
    std::vector<double>          residual_rotation;
};

// ── Internal helpers ──────────────────────────────────────────────────────────

namespace detail {

// ─────────────────────────────────────────────────────────────────────────────
// Priority-BFS queue entry
// Faces closer to the root (lower depth) are processed first, reducing the
// accumulation of trilateration errors for later-placed vertices.
// ─────────────────────────────────────────────────────────────────────────────
struct BFSEntry {
    int            depth;  ///< BFS depth = max(depth[v_src], depth[v_tgt]) + 1
    Halfedge_index h;      ///< halfedge pointing INTO the face to be placed
    /// min-heap: smaller depth = higher priority
    bool operator>(const BFSEntry& o) const { return depth > o.depth; }
};

using BFSQueue = std::priority_queue<BFSEntry,
                                     std::vector<BFSEntry>,
                                     std::greater<BFSEntry>>;

// ─────────────────────────────────────────────────────────────────────────────
// Euclidean trilateration — LEFT (CCW) side of p_a → p_b
// ─────────────────────────────────────────────────────────────────────────────
inline Eigen::Vector2d trilaterate_2d(
    const Eigen::Vector2d& pa, const Eigen::Vector2d& pb,
    double da, double db)
{
    Eigen::Vector2d ab = pb - pa;
    double d = ab.norm();
    if (d < 1e-14) return pa;
    Eigen::Vector2d e1 = ab / d;
    Eigen::Vector2d e2(-e1.y(), e1.x());
    double t  = (da*da - db*db + d*d) / (2.0 * d);
    double s2 = da*da - t*t;
    double s  = (s2 > 0.0) ? std::sqrt(s2) : 0.0;
    return pa + t*e1 + s*e2;
}

// ─────────────────────────────────────────────────────────────────────────────
// Spherical trilateration (S²) — LEFT side of geodesic arc pa → pb
// ─────────────────────────────────────────────────────────────────────────────
inline Eigen::Vector3d trilaterate_sph(
    const Eigen::Vector3d& pa, const Eigen::Vector3d& pb,
    double da, double db)
{
    double c     = pa.dot(pb);
    double denom = 1.0 - c*c;
    if (denom < 1e-14) return pa;
    double cda = std::cos(da), cdb = std::cos(db);
    double alpha = (cda - c*cdb) / denom;
    double beta  = (cdb - c*cda) / denom;
    Eigen::Vector3d cross  = pa.cross(pb);
    double          cn2    = cross.squaredNorm();
    double gamma2 = 1.0 - alpha*alpha - beta*beta - 2.0*alpha*beta*c;
    double gamma  = (gamma2 > 0.0 && cn2 > 1e-28) ? std::sqrt(gamma2 / cn2) : 0.0;
    Eigen::Vector3d p = alpha*pa + beta*pb + gamma*cross;
    double n = p.norm();
    return (n > 1e-14) ? (p / n) : pa;
}

// ─────────────────────────────────────────────────────────────────────────────
// Hyperbolic trilateration — exact via Möbius + hyperbolic law of cosines.
// LEFT (CCW) side of pa → pb in the Poincaré disk.
// ─────────────────────────────────────────────────────────────────────────────
inline Eigen::Vector2d trilaterate_hyp(
    const Eigen::Vector2d& pa, const Eigen::Vector2d& pb,
    double D, double da, double db)
{
    using C = std::complex<double>;
    if (D < 1e-12 || da < 1e-12) return pa;
    if (std::sinh(da) * std::sinh(D) < 1e-14) return pa;
    C a(pa.x(), pa.y()), b(pb.x(), pb.y());
    auto fwd = [](C z, C c_) { return (z-c_)/(C(1)-std::conj(c_)*z); };
    auto inv = [](C w, C c_) { return (w+c_)/(C(1)+std::conj(c_)*w); };
    double theta_b  = std::arg(fwd(b, a));
    double cos_alpha = (std::cosh(da)*std::cosh(D) - std::cosh(db))
                     / (std::sinh(da)*std::sinh(D));
    cos_alpha = std::max(-1.0, std::min(1.0, cos_alpha));
    C pc_origin = std::tanh(da*0.5) * std::exp(C(0.0, theta_b + std::acos(cos_alpha)));
    C pc_c = inv(pc_origin, a);
    return Eigen::Vector2d(pc_c.real(), pc_c.imag());
}

// ─────────────────────────────────────────────────────────────────────────────
// 3-D face area (cross product / 2)
// ─────────────────────────────────────────────────────────────────────────────
inline double face_3d_area(const ConformalMesh& mesh, Face_index f)
{
    Halfedge_index h = mesh.halfedge(f);
    auto pA = mesh.point(mesh.source(h));
    auto pB = mesh.point(mesh.target(h));
    auto pC = mesh.point(mesh.target(mesh.next(h)));
    Eigen::Vector3d ab(pB.x()-pA.x(), pB.y()-pA.y(), pB.z()-pA.z());
    Eigen::Vector3d ac(pC.x()-pA.x(), pC.y()-pA.y(), pC.z()-pA.z());
    return 0.5 * ab.cross(ac).norm();
}

// ─────────────────────────────────────────────────────────────────────────────
// Best root face: largest 3-D area, with 1.5× bonus for interior faces.
// ─────────────────────────────────────────────────────────────────────────────
inline Face_index best_root_face(const ConformalMesh& mesh)
{
    Face_index best;
    double best_score = -1.0;
    for (auto f : mesh.faces()) {
        double area = face_3d_area(mesh, f);
        bool interior = true;
        for (auto h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh))
            if (mesh.is_border(mesh.opposite(h))) { interior = false; break; }
        double score = area * (interior ? 1.5 : 1.0);
        if (score > best_score) { best_score = score; best = f; }
    }
    return best;
}

// ─────────────────────────────────────────────────────────────────────────────
// Euclidean bounding box of placed vertices (for multi-component offset)
// ─────────────────────────────────────────────────────────────────────────────
inline void bbox_2d(const std::vector<Eigen::Vector2d>& uv,
                    const std::vector<bool>& placed,
                    double& xmax)
{
    xmax = 0.0;
    for (std::size_t i = 0; i < uv.size(); ++i)
        if (placed[i]) xmax = std::max(xmax, uv[i].x());
}

// ─────────────────────────────────────────────────────────────────────────────
// Möbius centering — unweighted (fallback / Phase 6 compat.)
// ─────────────────────────────────────────────────────────────────────────────
inline void center_poincare_disk(std::vector<Eigen::Vector2d>& uv)
{
    if (uv.empty()) return;
    using C = std::complex<double>;
    C c(0);
    for (auto& p : uv) c += C(p.x(), p.y());
    c /= static_cast<double>(uv.size());
    double r = std::abs(c);
    if (r > 0.999) c *= 0.999/r;
    if (r < 1e-12) return;
    for (auto& p : uv) {
        C z(p.x(), p.y());
        C w = (z-c)/(C(1)-std::conj(c)*z);
        p = {w.real(), w.imag()};
    }
}

// ─────────────────────────────────────────────────────────────────────────────
// Face-area-weighted iterative Möbius centering (Fréchet mean, Phase 7).
// weights[i] = Voronoi area of vertex i.
// ─────────────────────────────────────────────────────────────────────────────
inline void center_poincare_disk_weighted(
    std::vector<Eigen::Vector2d>& uv, const std::vector<double>& weights)
{
    if (uv.empty()) return;
    using C = std::complex<double>;
    double total_w = 0.0; for (double w : weights) total_w += w;
    if (total_w < 1e-14) { center_poincare_disk(uv); return; }
    C c(0);
    for (int iter = 0; iter < 30; ++iter) {
        C mt(0);
        for (std::size_t i = 0; i < uv.size(); ++i) {
            C z(uv[i].x(), uv[i].y());
            mt += weights[i] * (z-c)/(C(1)-std::conj(c)*z);
        }
        mt /= total_w;
        C c_new = (mt+c)/(C(1)+std::conj(c)*mt);
        if (std::abs(c_new-c) < 1e-12) { c = c_new; break; }
        c = c_new;
    }
    double r = std::abs(c);
    if (r > 0.999) c *= 0.999/r;
    if (r < 1e-12) return;
    for (auto& p : uv) {
        C z(p.x(), p.y());
        C w = (z-c)/(C(1)-std::conj(c)*z);
        p = {w.real(), w.imag()};
    }
}

} // namespace detail

/// Compute per-vertex area weights (sum of 1/3 of each adjacent triangle area).
/// Used by area-weighted layout normalisation routines.
inline std::vector<double> compute_vertex_area_weights(const ConformalMesh& mesh)
{
    std::vector<double> w(mesh.number_of_vertices(), 0.0);
    for (auto f : mesh.faces()) {
        double area = detail::face_3d_area(mesh, f);
        for (auto h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh))
            w[static_cast<std::size_t>(mesh.target(h).idx())] += area / 3.0;
    }
    return w;
}

// ── Layout normalisation ──────────────────────────────────────────────────────

/// Euclidean canonical normalisation: translate centroid to the origin
/// and rotate the principal axis of the UV cloud onto the x-axis.
inline void normalise_euclidean(Layout2D& layout)
{
    if (!layout.success || layout.uv.empty()) return;
    const std::size_t n = layout.uv.size();
    Eigen::Vector2d mean = Eigen::Vector2d::Zero();
    for (auto& p : layout.uv) mean += p;
    mean /= static_cast<double>(n);
    for (auto& p : layout.uv) p -= mean;
    for (auto& p : layout.halfedge_uv) p -= mean;
    Eigen::Matrix2d cov = Eigen::Matrix2d::Zero();
    for (auto& p : layout.uv) cov += p * p.transpose();
    cov /= static_cast<double>(n);
    Eigen::SelfAdjointEigenSolver<Eigen::Matrix2d> eig(cov);
    Eigen::Vector2d major = eig.eigenvectors().col(1);
    double angle = -std::atan2(major.y(), major.x());
    Eigen::Matrix2d R;
    R << std::cos(angle), -std::sin(angle),
         std::sin(angle),  std::cos(angle);
    for (auto& p : layout.uv)         p = R * p;
    for (auto& p : layout.halfedge_uv) p = R * p;
}

/// Hyperbolic — face-area-weighted iterative Möbius centering.
inline void normalise_hyperbolic(Layout2D& layout, const ConformalMesh& mesh)
{
    if (!layout.success || layout.uv.empty()) return;
    auto weights = compute_vertex_area_weights(mesh);
    detail::center_poincare_disk_weighted(layout.uv, weights);
    // halfedge_uv follows the same Möbius map
    detail::center_poincare_disk(layout.halfedge_uv);
}
/// Hyperbolic canonical normalisation, mesh-free fallback: uniform
/// (unweighted) iterative Möbius centring of the Poincaré disk.
inline void normalise_hyperbolic(Layout2D& layout) // fallback without mesh
{
    if (!layout.success || layout.uv.empty()) return;
    detail::center_poincare_disk(layout.uv);
    detail::center_poincare_disk(layout.halfedge_uv);
}

/// Spherical canonical normalisation: rotate the layout so that the
/// per-vertex centroid (projected back to S²) coincides with the north
/// pole (Rodrigues rotation).
inline void normalise_spherical(Layout3D& layout)
{
    if (!layout.success || layout.pos.empty()) return;
    Eigen::Vector3d mean = Eigen::Vector3d::Zero();
    for (auto& p : layout.pos) mean += p;
    double len = mean.norm(); if (len < 1e-12) return;
    mean /= len;
    Eigen::Vector3d north(0, 0, 1);
    Eigen::Vector3d axis = mean.cross(north);
    double sin_a = axis.norm(), cos_a = mean.dot(north);
    if (sin_a < 1e-12) return;
    axis /= sin_a;
    Eigen::Matrix3d K;
    K <<    0, -axis.z(),  axis.y(),
       axis.z(),       0, -axis.x(),
      -axis.y(),  axis.x(),       0;
    Eigen::Matrix3d Rot = Eigen::Matrix3d::Identity() + sin_a*K + (1-cos_a)*K*K;
    for (auto& p : layout.pos) p = Rot * p;
}

// ─────────────────────────────────────────────────────────────────────────────
// BFS placement helpers shared across all three layout functions.
//
// set_face_huv: populate halfedge_uv for the three halfedges of face f,
//               given that the face was entered via halfedge h_enter.
//   h_enter: source=v_src, target=v_tgt   → huv = uv[v_src]
//   next(h_enter): source=v_tgt, target=v_new → huv = uv[v_tgt]
//   prev(h_enter): source=v_new, target=v_src → huv = p_new  (trilaterated)
// ─────────────────────────────────────────────────────────────────────────────
namespace detail {

inline void set_face_huv_2d(
    std::vector<Eigen::Vector2d>& huv,
    const ConformalMesh& mesh,
    Halfedge_index h,
    const std::vector<Eigen::Vector2d>& uv,
    const Eigen::Vector2d& p_new)
{
    auto s = [](std::size_t x){ return x; };
    huv[s(static_cast<std::size_t>(h.idx()))]               = uv[static_cast<std::size_t>(mesh.source(h).idx())];
    huv[s(static_cast<std::size_t>(mesh.next(h).idx()))]    = uv[static_cast<std::size_t>(mesh.target(h).idx())];
    huv[s(static_cast<std::size_t>(mesh.prev(h).idx()))]    = p_new;
}

inline void set_root_huv_2d(
    std::vector<Eigen::Vector2d>& huv,
    const ConformalMesh& mesh,
    Face_index f,
    const std::vector<Eigen::Vector2d>& uv)
{
    for (auto hf : CGAL::halfedges_around_face(mesh.halfedge(f), mesh))
        huv[static_cast<std::size_t>(hf.idx())] = uv[static_cast<std::size_t>(mesh.source(hf).idx())];
}

// ─────────────────────────────────────────────────────────────────────────────
// Euclidean holonomy via the developing map (genus-g closed surfaces).
//
// The per-vertex layout produced by euclidean_layout() places every face in a
// SINGLE consistent global frame (each vertex is placed once).  In that frame
// the holonomy is identically trivial — it is exactly the obstruction to such a
// single frame existing on the uncut surface.  To recover it we develop every
// face INDEPENDENTLY along a spanning tree of the dual graph that does not cross
// any cut edge, storing each face's own copy of its three corner positions
// (`hpos[h]` = global position of source(h) as developed inside face(h)).
//
// Two faces adjacent across a cut edge are NOT tree-adjacent, so they are
// developed via different tree branches and place the shared edge at two
// different locations.  The rigid motion identifying those two copies is the
// deck transformation of the generator loop (cut edge + tree path) — i.e. the
// holonomy.  For a flat cone metric (Θ ≡ 2π) the linear part is trivial, so the
// holonomy is the pure translation that identifies the two developments of the
// shared edge.  We recover it as a full rigid motion z ↦ a·z + b fitted to the
// correspondence (Bs,Bt) ↦ (As,At): the shared edge endpoints S,T as developed
// in face B map to the same endpoints as developed in face A.  Because both
// faces use the same edge length, |a| = 1 and the fit is an exact
// orientation-preserving isometry.  Its translation part b is the lattice
// generator ω consumed by compute_period_matrix(); the rotation magnitude
// |arg(a)| is returned as a convergence diagnostic.  For a perfectly flat cone
// metric (Θ ≡ 2π) a = 1 and b reduces to the midpoint displacement midA − midB
// (the previous formula), so converged flat tori are bit-for-bit unchanged.
struct EuclideanHolonomyResult {
    std::vector<Eigen::Vector2d> translations;       ///< ω_i = translation part b
    std::vector<double>          residual_rotation;  ///< |arg(a_i)| per cut edge
};

template <typename EdgeLenFn>
inline EuclideanHolonomyResult euclidean_holonomy(
    const ConformalMesh& mesh,
    const CutGraph&      cut,
    EdgeLenFn&&          edge_len)
{
    using C = std::complex<double>;
    const std::size_t nh = mesh.number_of_halfedges();
    const std::size_t nf = mesh.number_of_faces();

    std::vector<C>    hpos(nh, C(0.0, 0.0));   // pos of source(h) inside face(h)
    std::vector<bool> face_done(nf, false);

    auto place_root = [&](Face_index f) {
        Halfedge_index h0 = mesh.halfedge(f);
        Halfedge_index h1 = mesh.next(h0), h2 = mesh.next(h1);
        double lAB = edge_len(h0), lBC = edge_len(h1), lCA = edge_len(h2);
        Eigen::Vector2d A(0.0, 0.0), B(lAB, 0.0);
        Eigen::Vector2d Cc = trilaterate_2d(A, B, lCA, lBC);   // apex = source(h2)
        hpos[static_cast<std::size_t>(h0.idx())] = C(A.x(),  A.y());
        hpos[static_cast<std::size_t>(h1.idx())] = C(B.x(),  B.y());
        hpos[static_cast<std::size_t>(h2.idx())] = C(Cc.x(), Cc.y());
        face_done[static_cast<std::size_t>(f.idx())] = true;
    };

    auto develop = [&](Face_index root) {
        place_root(root);
        std::queue<Halfedge_index> q;     // halfedges pointing INTO unplaced faces
        auto enqueue = [&](Face_index f) {
            for (auto hf : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) {
                Halfedge_index ho = mesh.opposite(hf);
                if (mesh.is_border(ho)) continue;
                // Develop across the dual spanning tree T* ONLY.  Crossing any
                // other edge (a cut/generator edge OR a primal-tree edge) would
                // over-connect the development: the surface minus the 2g cut
                // edges is still non-simply-connected, so the immersion would
                // wrap around and place the two copies of a generator edge on
                // top of each other (zero/garbage holonomy).  Crossing only T*
                // unfolds the surface onto a disk (the fundamental polygon).
                if (!cut.is_dual_tree(mesh.edge(hf))) continue;
                Face_index fa = mesh.face(ho);
                if (!face_done[static_cast<std::size_t>(fa.idx())]) q.push(ho);
            }
        };
        enqueue(root);
        while (!q.empty()) {
            Halfedge_index h = q.front(); q.pop();
            Face_index f = mesh.face(h);
            if (face_done[static_cast<std::size_t>(f.idx())]) continue;

            Halfedge_index ho = mesh.opposite(h);
            // Shared edge endpoints, taken from the PARENT face's development:
            //   source(h) = target(ho) → parent pos hpos[next(ho)]
            //   target(h) = source(ho) → parent pos hpos[ho]
            C ps = hpos[static_cast<std::size_t>(mesh.next(ho).idx())];
            C pt = hpos[static_cast<std::size_t>(ho.idx())];
            Eigen::Vector2d A(ps.real(), ps.imag()), B(pt.real(), pt.imag());
            Eigen::Vector2d apex = trilaterate_2d(
                A, B, edge_len(mesh.prev(h)), edge_len(mesh.next(h)));

            hpos[static_cast<std::size_t>(h.idx())]            = ps;
            hpos[static_cast<std::size_t>(mesh.next(h).idx())] = pt;
            hpos[static_cast<std::size_t>(mesh.prev(h).idx())] = C(apex.x(), apex.y());
            face_done[static_cast<std::size_t>(f.idx())] = true;
            enqueue(f);
        }
    };

    develop(best_root_face(mesh));
    for (auto f : mesh.faces())
        if (!face_done[static_cast<std::size_t>(f.idx())]) develop(f);

    // ── Holonomy isometry per cut edge ─────────────────────────────────────────
    EuclideanHolonomyResult out;
    out.translations.reserve(cut.cut_edge_indices.size());
    out.residual_rotation.reserve(cut.cut_edge_indices.size());
    for (std::size_t ce : cut.cut_edge_indices) {
        Edge_index e = *std::next(mesh.edges().begin(), static_cast<std::ptrdiff_t>(ce));
        Halfedge_index h  = mesh.halfedge(e);
        Halfedge_index ho = mesh.opposite(h);
        if (mesh.is_border(h) || mesh.is_border(ho)
            || !face_done[static_cast<std::size_t>(mesh.face(h).idx())]
            || !face_done[static_cast<std::size_t>(mesh.face(ho).idx())]) {
            out.translations.push_back(Eigen::Vector2d::Zero());
            out.residual_rotation.push_back(0.0);
            continue;
        }
        // Face A = face(h) places edge e as halfedge h:  S=source(h), T=target(h)
        C As = hpos[static_cast<std::size_t>(h.idx())];
        C At = hpos[static_cast<std::size_t>(mesh.next(h).idx())];
        // Face B = face(ho) places the same edge as halfedge ho: source(ho)=T,
        //   target(ho)=S  →  S=pos of source(next(ho)), T=pos of source(ho)
        C Bt = hpos[static_cast<std::size_t>(ho.idx())];
        C Bs = hpos[static_cast<std::size_t>(mesh.next(ho).idx())];
        // Fit the deck isometry g(z) = a·z + b with (Bs,Bt) ↦ (As,At).
        // |a| = 1 exactly (shared edge length); b is the translation part.
        C edgeB = Bt - Bs;
        C a = (std::abs(edgeB) > 1e-300) ? (At - As) / edgeB : C(1.0, 0.0);
        C b = As - a * Bs;
        out.translations.push_back(Eigen::Vector2d(b.real(), b.imag()));
        out.residual_rotation.push_back(std::abs(std::arg(a)));
    }
    return out;
}

} // namespace detail

// ── Euclidean layout ──────────────────────────────────────────────────────────

/// Embed the mesh in ℝ² using the Euclidean metric encoded in x.
///
/// Runs priority-BFS trilateration: places vertices in order of BFS depth
/// from the root face (largest 3D area), so errors accumulate last.
/// For closed genus-g surfaces a CutGraph must be supplied — otherwise
/// the layout will have a seam discontinuity (Layout2D::has_seam = true).
///
/// \param mesh      Input surface mesh.  Must have lambda0 and v_idx set in maps.
/// \param x         DOF vector returned by newton_euclidean().
/// \param maps      EuclideanMaps (lambda0, v_idx, e_idx).
/// \param cut       Optional cut graph (compute_cut_graph()).  Pass nullptr for
///                  open meshes or if seams are acceptable.
/// \param holonomy  If non-null and cut != nullptr, receives the lattice
///                  translations ω_i ∈ ℂ per cut edge.
///                  Pass to compute_period_matrix() for the conformal modulus τ.
/// \param normalise If true, calls normalise_euclidean() on the result:
///                  centroid → origin, major axis → x-axis (PCA).
/// \return          Layout2D with .uv[v] (per-vertex UV) and
///                  .halfedge_uv[h] (per-halfedge UV for texture atlasing).
///
/// \note halfedge_uv differs from uv at seam edges: the two sides of a cut
///       carry different UV coordinates for proper GPU texture atlasing.
inline Layout2D euclidean_layout(
    ConformalMesh&             mesh,
    const std::vector<double>& x,
    const EuclideanMaps&       maps,
    const CutGraph*            cut      = nullptr,
    HolonomyData*              holonomy  = nullptr,
    bool                       normalise = false)
{
    const std::size_t nv = mesh.number_of_vertices();
    const std::size_t nf = mesh.number_of_faces();
    const std::size_t nh = mesh.number_of_halfedges();
    Layout2D result;
    result.uv.assign(nv, Eigen::Vector2d::Zero());
    result.halfedge_uv.assign(nh, Eigen::Vector2d::Zero());
    if (nf == 0) return result;

    auto get_u  = [&](Vertex_index v) {
        int iv = maps.v_idx[v]; return (iv>=0) ? x[static_cast<std::size_t>(iv)] : 0.0; };
    auto get_ue = [&](Edge_index e) {
        int ie = maps.e_idx[e]; return (ie>=0) ? x[static_cast<std::size_t>(ie)] : 0.0; };
    auto edge_len = [&](Halfedge_index h) {
        Edge_index e = mesh.edge(h);
        double lam = maps.lambda0[e] + get_u(mesh.source(h)) + get_u(mesh.target(h)) + get_ue(e);
        return std::exp(lam * 0.5); };

    std::vector<bool> vertex_placed(nv, false);
    std::vector<bool> face_placed(nf, false);
    std::vector<int>  vertex_depth(nv, std::numeric_limits<int>::max());

    auto place_component = [&](Face_index f_root, Eigen::Vector2d offset) {
        Halfedge_index h0 = mesh.halfedge(f_root);
        Halfedge_index h1 = mesh.next(h0), h2 = mesh.next(h1);
        Vertex_index vA = mesh.source(h0), vB = mesh.source(h1), vC = mesh.source(h2);
        double lAB = edge_len(h0), lBC = edge_len(h1), lCA = edge_len(h2);
        result.uv[vA.idx()] = offset;
        result.uv[vB.idx()] = offset + Eigen::Vector2d(lAB, 0.0);
        result.uv[vC.idx()] = detail::trilaterate_2d(result.uv[vA.idx()], result.uv[vB.idx()], lCA, lBC);
        vertex_placed[vA.idx()] = vertex_placed[vB.idx()] = vertex_placed[vC.idx()] = true;
        vertex_depth[vA.idx()]  = vertex_depth[vB.idx()]  = vertex_depth[vC.idx()]  = 0;
        face_placed[f_root.idx()] = true;
        detail::set_root_huv_2d(result.halfedge_uv, mesh, f_root, result.uv);

        detail::BFSQueue pq;
        auto enqueue = [&](Face_index f) {
            for (auto hf : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) {
                Halfedge_index ho = mesh.opposite(hf);
                if (mesh.is_border(ho)) continue;
                if (cut && cut->is_cut(mesh.edge(hf))) continue;
                Face_index fadj = mesh.face(ho);
                if (!face_placed[fadj.idx()]) {
                    int d = std::max(vertex_depth[mesh.source(ho).idx()],
                                     vertex_depth[mesh.target(ho).idx()]) + 1;
                    pq.push({d, ho});
                }
            }
        };
        enqueue(f_root);

        while (!pq.empty()) {
            auto [depth, h] = pq.top(); pq.pop();
            Face_index f = mesh.face(h);
            if (face_placed[f.idx()]) continue;

            Vertex_index v_src = mesh.source(h), v_tgt = mesh.target(h);
            Vertex_index v_new = mesh.target(mesh.next(h));
            Eigen::Vector2d p = detail::trilaterate_2d(
                result.uv[v_src.idx()], result.uv[v_tgt.idx()],
                edge_len(mesh.prev(h)), edge_len(mesh.next(h)));

            if (!vertex_placed[v_new.idx()]) {
                result.uv[v_new.idx()]     = p;
                vertex_placed[v_new.idx()] = true;
                vertex_depth[v_new.idx()]  = depth;
            } else {
                result.has_seam = true;
            }
            // halfedge_uv: p is the UV of v_new in THIS face (may differ from uv[v_new] at seam)
            detail::set_face_huv_2d(result.halfedge_uv, mesh, h, result.uv, p);
            face_placed[f.idx()] = true;
            enqueue(f);
        }
    };

    place_component(detail::best_root_face(mesh), Eigen::Vector2d::Zero());
    for (auto f : mesh.faces()) {
        if (face_placed[f.idx()]) continue;
        double xmax; detail::bbox_2d(result.uv, vertex_placed, xmax);
        place_component(f, Eigen::Vector2d(xmax * 1.1 + 1.0, 0.0));
    }
    result.success = true;

    // ── Holonomy ──────────────────────────────────────────────────────────────
    // The single-frame BFS layout above puts every face in one consistent frame,
    // so the holonomy read off it is identically trivial.  Instead develop each
    // face independently along a dual spanning tree that never crosses a cut edge
    // (detail::euclidean_holonomy): the displacement of each cut edge between the
    // two faces that share it is the lattice generator ω_i.
    if (cut && holonomy) {
        holonomy->cut_edge_indices = cut->cut_edge_indices;
        holonomy->mobius_maps.clear();
        auto eh = detail::euclidean_holonomy(mesh, *cut, edge_len);
        holonomy->translations      = std::move(eh.translations);
        holonomy->residual_rotation = std::move(eh.residual_rotation);

        // Preserve the per-cut-edge seam UV in halfedge_uv (texture atlas), as
        // before, so HalfedgeUV-based tests still see the seam-crossing layout.
        for (std::size_t ce_idx : cut->cut_edge_indices) {
            Edge_index e = *std::next(mesh.edges().begin(), static_cast<std::ptrdiff_t>(ce_idx));
            Halfedge_index h  = mesh.halfedge(e);
            Halfedge_index ho = mesh.opposite(h);
            Halfedge_index hx = mesh.is_border(ho) ? h : ho;
            if (mesh.is_border(hx)) continue;
            Vertex_index vs = mesh.source(hx), vt = mesh.target(hx), vn = mesh.target(mesh.next(hx));
            if (!vertex_placed[vn.idx()]) continue;
            Eigen::Vector2d p_tri = detail::trilaterate_2d(
                result.uv[vs.idx()], result.uv[vt.idx()],
                edge_len(mesh.prev(hx)), edge_len(mesh.next(hx)));
            detail::set_face_huv_2d(result.halfedge_uv, mesh, hx, result.uv, p_tri);
        }
    }

    if (normalise) normalise_euclidean(result);
    return result;
}

// ── Spherical layout ──────────────────────────────────────────────────────────

/// Embed the mesh on the unit sphere S² using the spherical metric encoded in x.
///
/// Runs priority-BFS trilateration using the spherical law of cosines.
/// Typical use: genus-0 (sphere-like) surfaces after newton_spherical().
///
/// \param mesh      Input genus-0 surface mesh.
/// \param x         DOF vector returned by newton_spherical().
/// \param maps      SphericalMaps.
/// \param cut       Optional cut graph (rarely needed for genus-0).
/// \param holonomy  If non-null, receives rotational holonomies (spherical).
/// \param normalise If true, calls normalise_spherical(): rotates the centroid
///                  to the north pole (Rodrigues rotation formula).
/// \return          Layout3D with .xyz[v] ∈ S² ⊂ ℝ³ for each vertex.
inline Layout3D spherical_layout(
    ConformalMesh&             mesh,
    const std::vector<double>& x,
    const SphericalMaps&       maps,
    const CutGraph*            cut      = nullptr,
    HolonomyData*              holonomy  = nullptr,
    bool                       normalise = false)
{
    const std::size_t nv = mesh.number_of_vertices();
    const std::size_t nf = mesh.number_of_faces();
    Layout3D result;
    result.pos.assign(nv, Eigen::Vector3d::Zero());
    if (nf == 0) return result;

    auto get_u  = [&](Vertex_index v) {
        int iv = maps.v_idx[v]; return (iv>=0) ? x[static_cast<std::size_t>(iv)] : 0.0; };
    auto get_ue = [&](Edge_index e) {
        int ie = maps.e_idx[e]; return (ie>=0) ? x[static_cast<std::size_t>(ie)] : 0.0; };
    auto arc_len = [&](Halfedge_index h) {
        Edge_index e = mesh.edge(h);
        double lam = maps.lambda0[e] + get_u(mesh.source(h)) + get_u(mesh.target(h)) + get_ue(e);
        return 2.0 * std::asin(std::min(std::exp(lam*0.5), 1.0)); };

    std::vector<bool> vertex_placed(nv, false);
    std::vector<bool> face_placed(nf, false);
    std::vector<int>  vertex_depth(nv, std::numeric_limits<int>::max());

    auto place_component = [&](Face_index f_root) {
        Halfedge_index h0 = mesh.halfedge(f_root);
        Halfedge_index h1 = mesh.next(h0), h2 = mesh.next(h1);
        Vertex_index vA = mesh.source(h0), vB = mesh.source(h1), vC = mesh.source(h2);
        double lAB = arc_len(h0), lBC = arc_len(h1), lCA = arc_len(h2);
        result.pos[vA.idx()] = Eigen::Vector3d(0, 0, 1);
        result.pos[vB.idx()] = Eigen::Vector3d(std::sin(lAB), 0, std::cos(lAB));
        result.pos[vC.idx()] = detail::trilaterate_sph(result.pos[vA.idx()], result.pos[vB.idx()], lCA, lBC);
        vertex_placed[vA.idx()] = vertex_placed[vB.idx()] = vertex_placed[vC.idx()] = true;
        vertex_depth[vA.idx()]  = vertex_depth[vB.idx()]  = vertex_depth[vC.idx()]  = 0;
        face_placed[f_root.idx()] = true;

        detail::BFSQueue pq;
        auto enqueue = [&](Face_index f) {
            for (auto hf : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) {
                Halfedge_index ho = mesh.opposite(hf);
                if (mesh.is_border(ho)) continue;
                if (cut && cut->is_cut(mesh.edge(hf))) continue;
                Face_index fadj = mesh.face(ho);
                if (!face_placed[fadj.idx()]) {
                    int d = std::max(vertex_depth[mesh.source(ho).idx()],
                                     vertex_depth[mesh.target(ho).idx()]) + 1;
                    pq.push({d, ho});
                }
            }
        };
        enqueue(f_root);

        while (!pq.empty()) {
            auto [depth, h] = pq.top(); pq.pop();
            Face_index f = mesh.face(h);
            if (face_placed[f.idx()]) continue;
            Vertex_index vs = mesh.source(h), vt = mesh.target(h), vn = mesh.target(mesh.next(h));
            Eigen::Vector3d p = detail::trilaterate_sph(
                result.pos[vs.idx()], result.pos[vt.idx()],
                arc_len(mesh.prev(h)), arc_len(mesh.next(h)));
            if (!vertex_placed[vn.idx()]) {
                result.pos[vn.idx()] = p; vertex_placed[vn.idx()] = true; vertex_depth[vn.idx()] = depth;
            } else result.has_seam = true;
            face_placed[f.idx()] = true; enqueue(f);
        }
    };

    place_component(detail::best_root_face(mesh));
    for (auto f : mesh.faces()) if (!face_placed[f.idx()]) place_component(f);
    result.success = true;

    // KNOWN LIMITATION (latent — every current caller passes holonomy == nullptr).
    // This block extracts holonomy from a single full-surface development (the BFS
    // above crosses every non-cut edge), then reads off an apex trilateration on one
    // side of each seam.  That is the same flawed pattern that produced garbage τ for
    // the Euclidean path; the correct approach is detail::euclidean_holonomy, which
    // develops across only the dual spanning tree (is_dual_tree) and measures the
    // shared-edge displacement between two independent developments.  Until a
    // detail::spherical_holonomy mirror exists (Phase 9c/10, see research-track.md),
    // these spherical translations are not trustworthy for genus g ≥ 1.
    if (cut && holonomy) {
        holonomy->cut_edge_indices = cut->cut_edge_indices;
        holonomy->translations.clear();
        holonomy->translations.reserve(cut->cut_edge_indices.size());
        holonomy->mobius_maps.clear();
        for (std::size_t ce_idx : cut->cut_edge_indices) {
            Edge_index e = *std::next(mesh.edges().begin(), static_cast<std::ptrdiff_t>(ce_idx));
            Halfedge_index hx = mesh.is_border(mesh.opposite(mesh.halfedge(e)))
                               ? mesh.halfedge(e) : mesh.opposite(mesh.halfedge(e));
            if (mesh.is_border(hx)) { holonomy->translations.push_back(Eigen::Vector2d::Zero()); continue; }
            Vertex_index vs = mesh.source(hx), vt = mesh.target(hx), vn = mesh.target(mesh.next(hx));
            if (!vertex_placed[vn.idx()]) { holonomy->translations.push_back(Eigen::Vector2d::Zero()); continue; }
            Eigen::Vector3d p_tri = detail::trilaterate_sph(
                result.pos[vs.idx()], result.pos[vt.idx()],
                arc_len(mesh.prev(hx)), arc_len(mesh.next(hx)));
            Eigen::Vector3d diff = p_tri - result.pos[vn.idx()];
            // Note: spherical holonomy is geometrically a 3-D rotation, not a 2-D
            // translation.  The Vector2d here stores only the (x,y) component of the
            // S²-position difference across the cut, which is an approximation.
            // For accurate spherical holonomy (rotation axis + angle) use the full
            // 3-D positions in result.pos[] directly.  Phase 10+ will replace this
            // with a proper SO(3) representation.
            holonomy->translations.push_back(Eigen::Vector2d(diff.x(), diff.y()));
        }
    }

    if (normalise) normalise_spherical(result);
    return result;
}

// ── HyperIdeal layout (Poincaré disk) — exact trilateration ──────────────────

/// Embed the mesh in the Poincaré disk (H²) using the hyperbolic metric encoded in x.
///
/// Runs priority-BFS trilateration using exact Möbius-isometric placement:
/// each new vertex is located by solving the hyperbolic law of cosines and
/// applying a Möbius map to position it in the disk.
/// For closed genus-g surfaces (g ≥ 1) a CutGraph is required.
///
/// \param mesh      Input genus-g surface mesh (g ≥ 1).
/// \param x         DOF vector returned by newton_hyper_ideal()
///                  (vertex b_v and edge a_e variables).
/// \param maps      HyperIdealMaps (lambda0, v_idx, e_idx).
/// \param cut       CutGraph from compute_cut_graph().  Required for closed surfaces.
/// \param holonomy  If non-null, receives the Möbius maps T_i ∈ SU(1,1) per cut edge.
///                  Pass to compute_period_matrix() for holonomy analysis.
/// \param normalise If true, calls normalise_hyperbolic(): iterative face-area-weighted
///                  Möbius centring (Fréchet mean, 30 iterations) → disk origin.
/// \return          Layout2D with .uv[v] ∈ Poincaré disk (|uv| < 1) for each vertex,
///                  and .halfedge_uv[h] for seam-aware texture atlasing.
///
/// \note All vertex positions satisfy |uv[v]| < 1 (inside the Poincaré disk)
///       if the metric is hyperbolic.  Points on or outside the boundary indicate
///       a non-hyperbolic metric (Gauss–Bonnet violation).
inline Layout2D hyper_ideal_layout(
    ConformalMesh&             mesh,
    const std::vector<double>& x,
    const HyperIdealMaps&      maps,
    const CutGraph*            cut      = nullptr,
    HolonomyData*              holonomy  = nullptr,
    bool                       normalise = false)
{
    const std::size_t nv = mesh.number_of_vertices();
    const std::size_t nf = mesh.number_of_faces();
    const std::size_t nh = mesh.number_of_halfedges();
    Layout2D result;
    result.uv.assign(nv, Eigen::Vector2d::Zero());
    result.halfedge_uv.assign(nh, Eigen::Vector2d::Zero());
    if (nf == 0) return result;

    auto get_b = [&](Vertex_index v) {
        int iv = maps.v_idx[v]; return (iv>=0) ? x[static_cast<std::size_t>(iv)] : 0.0; };
    auto get_a = [&](Edge_index e) {
        int ie = maps.e_idx[e]; return (ie>=0) ? x[static_cast<std::size_t>(ie)] : 0.0; };
    auto hyp_dist = [&](Halfedge_index h) {
        double bi = get_b(mesh.source(h)), bj = get_b(mesh.target(h)), a = get_a(mesh.edge(h));
        return std::acosh(std::max(1.0, std::cosh(bi+a*0.5)*std::cosh(bj+a*0.5)-std::sinh(bi)*std::sinh(bj))); };

    std::vector<bool> vertex_placed(nv, false);
    std::vector<bool> face_placed(nf, false);
    std::vector<int>  vertex_depth(nv, std::numeric_limits<int>::max());

    auto place_component = [&](Face_index f_root) {
        Halfedge_index h0 = mesh.halfedge(f_root);
        Halfedge_index h1 = mesh.next(h0), h2 = mesh.next(h1);
        Vertex_index vA = mesh.source(h0), vB = mesh.source(h1), vC = mesh.source(h2);
        double dAB = hyp_dist(h0), dBC = hyp_dist(h1), dCA = hyp_dist(h2);
        result.uv[vA.idx()] = Eigen::Vector2d::Zero();
        result.uv[vB.idx()] = Eigen::Vector2d(std::tanh(dAB*0.5), 0.0);
        result.uv[vC.idx()] = detail::trilaterate_hyp(result.uv[vA.idx()], result.uv[vB.idx()], dAB, dCA, dBC);
        vertex_placed[vA.idx()] = vertex_placed[vB.idx()] = vertex_placed[vC.idx()] = true;
        vertex_depth[vA.idx()]  = vertex_depth[vB.idx()]  = vertex_depth[vC.idx()]  = 0;
        face_placed[f_root.idx()] = true;
        detail::set_root_huv_2d(result.halfedge_uv, mesh, f_root, result.uv);

        detail::BFSQueue pq;
        auto enqueue = [&](Face_index f) {
            for (auto hf : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) {
                Halfedge_index ho = mesh.opposite(hf);
                if (mesh.is_border(ho)) continue;
                if (cut && cut->is_cut(mesh.edge(hf))) continue;
                Face_index fadj = mesh.face(ho);
                if (!face_placed[fadj.idx()]) {
                    int d = std::max(vertex_depth[mesh.source(ho).idx()],
                                     vertex_depth[mesh.target(ho).idx()]) + 1;
                    pq.push({d, ho});
                }
            }
        };
        enqueue(f_root);

        while (!pq.empty()) {
            auto [depth, h] = pq.top(); pq.pop();
            Face_index f = mesh.face(h);
            if (face_placed[f.idx()]) continue;
            Vertex_index vs = mesh.source(h), vt = mesh.target(h), vn = mesh.target(mesh.next(h));
            double D = hyp_dist(h), da = hyp_dist(mesh.prev(h)), db = hyp_dist(mesh.next(h));
            Eigen::Vector2d p = detail::trilaterate_hyp(result.uv[vs.idx()], result.uv[vt.idx()], D, da, db);
            if (!vertex_placed[vn.idx()]) {
                result.uv[vn.idx()] = p; vertex_placed[vn.idx()] = true; vertex_depth[vn.idx()] = depth;
            } else result.has_seam = true;
            detail::set_face_huv_2d(result.halfedge_uv, mesh, h, result.uv, p);
            face_placed[f.idx()] = true; enqueue(f);
        }
    };

    place_component(detail::best_root_face(mesh));
    for (auto f : mesh.faces()) if (!face_placed[f.idx()]) place_component(f);
    result.success = true;

    // ── Möbius-map holonomy ───────────────────────────────────────────────────
    // KNOWN LIMITATION (latent — every current caller passes holonomy == nullptr).
    // Like the spherical block above, this reads the Möbius deck transformation from
    // a single full-surface development (BFS crosses all non-cut edges) and one-sided
    // apex trilateration — the same flawed pattern fixed for the Euclidean path by
    // detail::euclidean_holonomy (develop across the dual tree only, measure seam
    // displacement between two independent developments).  A faithful
    // detail::hyperbolic_holonomy is Phase 9c/10 work and additionally requires
    // cpp_dec_float_50: products of these generators grow exponentially, so verifying
    // the group relation ∏gᵢ = Id overflows double (see CLAUDE.md, research-track.md).
    // Until then these mobius_maps are NOT correct for genus g ≥ 2 uniformization.
    if (cut && holonomy) {
        holonomy->cut_edge_indices = cut->cut_edge_indices;
        holonomy->translations.clear();
        holonomy->mobius_maps.clear();
        holonomy->mobius_maps.reserve(cut->cut_edge_indices.size());
        for (std::size_t ce_idx : cut->cut_edge_indices) {
            Edge_index e = *std::next(mesh.edges().begin(), static_cast<std::ptrdiff_t>(ce_idx));
            Halfedge_index hcut = mesh.halfedge(e), ho = mesh.opposite(hcut);
            Halfedge_index hx   = mesh.is_border(ho) ? hcut : ho;
            if (mesh.is_border(hx)) { holonomy->mobius_maps.push_back(MobiusMap::identity()); continue; }
            Vertex_index vs = mesh.source(hx), vt = mesh.target(hx), vn = mesh.target(mesh.next(hx));
            if (!vertex_placed[vn.idx()]) { holonomy->mobius_maps.push_back(MobiusMap::identity()); continue; }
            double D = hyp_dist(hx), da = hyp_dist(mesh.prev(hx)), db = hyp_dist(mesh.next(hx));
            Eigen::Vector2d p_tri = detail::trilaterate_hyp(result.uv[vs.idx()], result.uv[vt.idx()], D, da, db);
            detail::set_face_huv_2d(result.halfedge_uv, mesh, hx, result.uv, p_tri);
            using C = std::complex<double>;
            // Möbius deck transformation T across cut edge (vs,vt):
            // T is the unique Möbius isometry of the Poincaré disk that:
            //   - fixes vs and vt  (z1=w1, z2=w2: the cut-edge endpoints are
            //     identified across the seam, so T maps each to itself)
            //   - maps vn (placed side) → p_tri (virtual side)
            // This uniquely determines the hyperbolic translation/rotation
            // along the geodesic through vs and vt.
            holonomy->mobius_maps.push_back(MobiusMap::from_three(
                C(result.uv[vs.idx()].x(), result.uv[vs.idx()].y()),
                C(result.uv[vs.idx()].x(), result.uv[vs.idx()].y()),
                C(result.uv[vt.idx()].x(), result.uv[vt.idx()].y()),
                C(result.uv[vt.idx()].x(), result.uv[vt.idx()].y()),
                C(result.uv[vn.idx()].x(), result.uv[vn.idx()].y()),
                C(p_tri.x(), p_tri.y())));
        }
    }

    if (normalise) normalise_hyperbolic(result, mesh);
    return result;
}

// ── Convenience: save layout as OFF ──────────────────────────────────────────

/// Write a 2-D layout to disk in OFF format with z = 0.  Convenience
/// helper for quickly inspecting the UV result in any OFF viewer.
inline void save_layout_off(
    const std::string& path, ConformalMesh& mesh, const Layout2D& layout)
{
    std::ofstream ofs(path);
    ofs << "OFF\n" << mesh.number_of_vertices() << " " << mesh.number_of_faces() << " 0\n";
    for (auto v : mesh.vertices()) { auto& p = layout.uv[v.idx()]; ofs << p.x() << " " << p.y() << " 0\n"; }
    for (auto f : mesh.faces()) {
        ofs << "3";
        for (auto h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) ofs << " " << mesh.target(h).idx();
        ofs << "\n";
    }
}

/// Write a 3-D (spherical) layout to disk in OFF format.
inline void save_layout_off(
    const std::string& path, ConformalMesh& mesh, const Layout3D& layout)
{
    std::ofstream ofs(path);
    ofs << "OFF\n" << mesh.number_of_vertices() << " " << mesh.number_of_faces() << " 0\n";
    for (auto v : mesh.vertices()) { auto& p = layout.pos[v.idx()]; ofs << p.x() << " " << p.y() << " " << p.z() << "\n"; }
    for (auto f : mesh.faces()) {
        ofs << "3";
        for (auto h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) ofs << " " << mesh.target(h).idx();
        ofs << "\n";
    }
}

} // namespace conformallab