#pragma once
// layout.hpp
//
// Phase 5/6 — 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 a root face, place its three vertices analytically.           │
// │  3. BFS over the dual graph: for each adjacent face, the shared edge is  │
// │     already placed; trilaterate the third vertex.                        │
// │  4. (Phase 6) If a CutGraph is supplied, cut edges are treated as        │
// │     boundary; crossing them records a HolonomyData entry instead.        │
// │                                                                          │
// │  For open meshes the placement is globally consistent.                   │
// │  For closed meshes without a cut: first visit wins, has_seam = true.    │
// │  For closed meshes with a CutGraph: each vertex is placed exactly once;  │
// │  the holonomy (translation / Möbius) of each cut is recorded.           │
// └──────────────────────────────────────────────────────────────────────────┘
//
// Euclidean trilateration   — exact (analytic formula in ℝ²)
// Spherical trilateration   — exact (spherical law of cosines on S²)
// Hyperbolic trilateration  — exact (hyperbolic law of cosines + Möbius maps
//                             in the Poincaré disk model)
//
// Normalisation
//   normalise_euclidean(layout)   — centroid → origin, major axis → x-axis
//   normalise_hyperbolic(layout)  — Möbius map centering to origin
//   normalise_spherical(layout)   — rotate centroid to north pole
//
// Holonomy (Euclidean, genus-g surfaces with CutGraph)
//   HolonomyData.translations[i]  — translation vector ω_i for cut edge i
//   (for a flat torus: ω_1, ω_2 are the lattice generators)

#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 <fstream>
#include <string>

namespace conformallab {

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

struct Layout2D {
    std::vector<Eigen::Vector2d> uv;      ///< uv[v.idx()] = 2-D position
    bool success   = false;
    bool has_seam  = false;               ///< true if mesh is closed and no CutGraph given
};

struct Layout3D {
    std::vector<Eigen::Vector3d> pos;     ///< pos[v.idx()] = 3-D position on S²
    bool success  = false;
    bool has_seam = false;
};

/// Per-cut-edge holonomy for the Euclidean case.
/// translations[i] is the translation ω associated with cut_edge_indices[i]
/// of the CutGraph.  For a flat torus, ω_1 and ω_2 are the lattice generators.
struct HolonomyData {
    std::vector<Eigen::Vector2d> translations;   // one per cut edge
    std::vector<std::size_t>     cut_edge_indices;
};

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

namespace detail {

// ─────────────────────────────────────────────────────────────────────────────
// Euclidean trilateration
// Given placed p_a, p_b and distances d_a (from p_a), d_b (from p_b),
// return the point on the LEFT (CCW) side of the directed edge 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());   // CCW perpendicular
    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;               // s > 0 → left side
}

// ─────────────────────────────────────────────────────────────────────────────
// Spherical trilateration
// Given unit vectors pa, pb on S² and arc-lengths da, db to the new point,
// return the unit vector on the LEFT side of the geodesic arc pa → pb.
//
// Solution: p_c = α·pa + β·pb + γ·(pa × pb),  γ > 0.
// Constraints: pa·p_c = cos(da),  pb·p_c = cos(db),  |p_c| = 1.
// ─────────────────────────────────────────────────────────────────────────────
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          cross_n2 = cross.squaredNorm();
    double gamma2 = 1.0 - alpha*alpha - beta*beta - 2.0*alpha*beta*c;
    double gamma  = (gamma2 > 0.0 && cross_n2 > 1e-28)
                    ? std::sqrt(gamma2 / cross_n2) : 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 map + hyperbolic law of cosines.
//
// pa, pb: Poincaré disk coordinates of two already-placed vertices
// D:  hyperbolic distance pa → pb  (from the DOF vector)
// da: hyperbolic distance pa → pc  (new vertex)
// db: hyperbolic distance pb → pc  (new vertex)
//
// Returns the Poincaré disk coordinate of pc on the LEFT (CCW) side of pa→pb.
//
// Algorithm:
//   1. Map pa → 0 via Möbius T: z ↦ (z−pa)/(1−conj(pa)·z)
//   2. θ_b = arg(T(pb))   — direction from origin towards T(pb)
//   3. Angle α at pa from the hyperbolic law of cosines:
//        cosh(db) = cosh(da)·cosh(D) − sinh(da)·sinh(D)·cos(α)
//        → cos(α) = (cosh(da)·cosh(D) − cosh(db)) / (sinh(da)·sinh(D))
//   4. pc in origin frame: tanh(da/2)·exp(i·(θ_b + α))   [left = +α]
//   5. Map back: T⁻¹(w) = (w + pa)/(1 + conj(pa)·w)
// ─────────────────────────────────────────────────────────────────────────────
inline Eigen::Vector2d trilaterate_hyp(
    const Eigen::Vector2d& pa, const Eigen::Vector2d& pb,
    double D,   // hyperbolic distance pa → pb
    double da,  // hyperbolic distance pa → pc
    double db)  // hyperbolic distance pb → pc
{
    using C = std::complex<double>;

    // Degenerate guard: very short edges or coincident points
    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());
    C b(pb.x(), pb.y());

    // Möbius map: T_a(z) = (z − a) / (1 − conj(a)·z)
    auto mobius_fwd = [](C z, C center) -> C {
        return (z - center) / (C(1.0) - std::conj(center) * z);
    };
    // Inverse: T_a⁻¹(w) = (w + a) / (1 + conj(a)·w)
    auto mobius_inv = [](C w, C center) -> C {
        return (w + center) / (C(1.0) + std::conj(center) * w);
    };

    C b_mapped = mobius_fwd(b, a);
    double theta_b = std::arg(b_mapped);

    // Hyperbolic law of cosines for angle α at pa:
    // cosh(db) = cosh(da)·cosh(D) − sinh(da)·sinh(D)·cos(α)
    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));
    double alpha = std::acos(cos_alpha);   // ∈ [0, π], sin > 0 for left side

    // pc in the frame where pa = 0 and pb is on the positive real axis,
    // then rotate back by theta_b:
    C pc_origin = std::tanh(da * 0.5) * std::exp(C(0.0, theta_b + alpha));

    // Map back to original disk
    C pc_c = mobius_inv(pc_origin, a);
    return Eigen::Vector2d(pc_c.real(), pc_c.imag());
}

// ─────────────────────────────────────────────────────────────────────────────
// Möbius centering of a Poincaré-disk layout
// Applies T_c(z) = (z − c)/(1 − conj(c)·z) where c is the Euclidean centroid
// of all vertex positions.  Maps c → 0, i.e. centers the cloud in the disk.
// ─────────────────────────────────────────────────────────────────────────────
inline void center_poincare_disk(std::vector<Eigen::Vector2d>& uv)
{
    if (uv.empty()) return;
    using C = std::complex<double>;

    // Euclidean centroid (good enough for moderate displacements from origin)
    C centroid(0.0, 0.0);
    for (auto& p : uv) centroid += C(p.x(), p.y());
    centroid /= static_cast<double>(uv.size());

    // Clamp to strictly inside the disk
    double r = std::abs(centroid);
    if (r > 0.999) centroid *= (0.999 / r);
    if (r < 1e-12) return;  // already centered

    for (auto& p : uv) {
        C z(p.x(), p.y());
        C w = (z - centroid) / (C(1.0) - std::conj(centroid) * z);
        p = Eigen::Vector2d(w.real(), w.imag());
    }
}

} // namespace detail

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

/// Euclidean: translate centroid to origin, rotate major axis to x-axis (PCA).
inline void normalise_euclidean(Layout2D& layout)
{
    if (!layout.success || layout.uv.empty()) return;
    const std::size_t n = layout.uv.size();

    // Centroid
    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;

    // PCA: 2×2 covariance
    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);
    // Eigenvectors in ascending order — we want the largest (index 1)
    Eigen::Vector2d major = eig.eigenvectors().col(1);
    // Rotation that aligns major axis with x-axis
    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;
}

/// Hyperbolic: Möbius map centering the Euclidean centroid to the disk origin.
inline void normalise_hyperbolic(Layout2D& layout)
{
    if (!layout.success || layout.uv.empty()) return;
    detail::center_poincare_disk(layout.uv);
}

/// Spherical: rotate so the Euclidean centroid of vertex positions points
/// towards the north pole (0, 0, 1).
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;  // unit vector of mean direction

    Eigen::Vector3d north(0.0, 0.0, 1.0);
    Eigen::Vector3d axis = mean.cross(north);
    double sin_a = axis.norm();
    double cos_a = mean.dot(north);
    if (sin_a < 1e-12) return;  // already at north (or south) pole
    axis /= sin_a;

    // Rodrigues rotation matrix
    Eigen::Matrix3d K;
    K <<     0.0, -axis.z(),  axis.y(),
         axis.z(),      0.0, -axis.x(),
        -axis.y(),  axis.x(),      0.0;
    Eigen::Matrix3d R = Eigen::Matrix3d::Identity()
                      + sin_a * K + (1.0 - cos_a) * K * K;

    for (auto& p : layout.pos) p = R * p;
}

// ── Euclidean layout ──────────────────────────────────────────────────────────
//
//   l_ij = exp( (λ°_ij + u_i + u_j + λ_e) / 2 )
//
// Optional CutGraph: if non-null, cut edges are treated as boundary;
// holonomy translations are computed for each cut edge and returned via
// *holonomy (if non-null).
// ─────────────────────────────────────────────────────────────────────────────
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();
    Layout2D result;
    result.uv.assign(nv, Eigen::Vector2d::Zero());
    if (mesh.number_of_faces() == 0) return result;

    auto get_u = [&](Vertex_index v) -> double {
        int iv = maps.v_idx[v]; return (iv >= 0) ? x[static_cast<std::size_t>(iv)] : 0.0;
    };
    auto get_ue = [&](Edge_index e) -> double {
        int ie = maps.e_idx[e]; return (ie >= 0) ? x[static_cast<std::size_t>(ie)] : 0.0;
    };
    auto edge_len = [&](Halfedge_index h) -> double {
        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(mesh.number_of_faces(), false);

    // ── Place root face ───────────────────────────────────────────────────────
    Face_index     f0 = *mesh.faces().begin();
    Halfedge_index h0 = mesh.halfedge(f0);
    Halfedge_index h1 = mesh.next(h0);
    Halfedge_index h2 = mesh.next(h1);
    Vertex_index   vA = mesh.source(h0);
    Vertex_index   vB = mesh.source(h1);
    Vertex_index   vC = mesh.source(h2);
    double lAB = edge_len(h0);
    double lCA = edge_len(h2);
    double lBC = edge_len(h1);

    result.uv[vA.idx()] = Eigen::Vector2d(0.0, 0.0);
    result.uv[vB.idx()] = 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;
    face_placed[f0.idx()] = true;

    // ── BFS ───────────────────────────────────────────────────────────────────
    std::queue<Halfedge_index> q;
    auto enqueue = [&](Face_index f) {
        for (Halfedge_index h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) {
            Halfedge_index h_opp = mesh.opposite(h);
            if (mesh.is_border(h_opp)) continue;
            // Treat cut edges as boundary (don't cross them in BFS)
            if (cut && cut->is_cut(mesh.edge(h))) continue;
            Face_index f_adj = mesh.face(h_opp);
            if (!face_placed[f_adj.idx()])
                q.push(h_opp);
        }
    };
    enqueue(f0);

    while (!q.empty()) {
        Halfedge_index h = q.front(); q.pop();
        Face_index f = mesh.face(h);
        if (face_placed[f.idx()]) continue;

        Vertex_index v_src = mesh.source(h);
        Vertex_index 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;
        } else {
            result.has_seam = true;
        }
        face_placed[f.idx()] = true;
        enqueue(f);
    }

    result.success = true;

    // ── Holonomy: compute translation for each cut edge ───────────────────────
    // For each cut edge e = (u,v), look at the face on the far side (h_opp).
    // Trilaterate the third vertex w of that face from uv[u] and uv[v],
    // and compare to the actual position uv[w].
    // Translation ω = trilaterated_position(w) − uv[w].
    if (cut && holonomy) {
        holonomy->cut_edge_indices = cut->cut_edge_indices;
        holonomy->translations.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 h  = mesh.halfedge(e);
            Halfedge_index h_opp = mesh.opposite(h);

            // Choose the halfedge whose face was NOT placed during BFS
            // (the one across the cut from the main BFS sweep).
            // If both sides were placed, use h_opp; if neither, skip.
            Halfedge_index h_cross = h_opp;
            if (mesh.is_border(h_cross)) h_cross = h;
            if (mesh.is_border(h_cross)) {
                holonomy->translations.push_back(Eigen::Vector2d::Zero());
                continue;
            }

            Vertex_index v_src = mesh.source(h_cross);
            Vertex_index v_tgt = mesh.target(h_cross);
            Vertex_index v_new = mesh.target(mesh.next(h_cross));

            if (!vertex_placed[v_new.idx()]) {
                holonomy->translations.push_back(Eigen::Vector2d::Zero());
                continue;
            }

            Eigen::Vector2d p_tri = detail::trilaterate_2d(
                result.uv[v_src.idx()], result.uv[v_tgt.idx()],
                edge_len(mesh.prev(h_cross)),
                edge_len(mesh.next(h_cross)));

            holonomy->translations.push_back(p_tri - result.uv[v_new.idx()]);
        }
    }

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

// ── Spherical layout ──────────────────────────────────────────────────────────
//
//   l_ij = 2·arcsin( min(exp(Λ_ij / 2), 1) )
//
// ─────────────────────────────────────────────────────────────────────────────
inline Layout3D spherical_layout(
    ConformalMesh&             mesh,
    const std::vector<double>& x,
    const SphericalMaps&       maps,
    bool                       normalise = false)
{
    const std::size_t nv = mesh.number_of_vertices();
    Layout3D result;
    result.pos.assign(nv, Eigen::Vector3d::Zero());
    if (mesh.number_of_faces() == 0) return result;

    auto get_u = [&](Vertex_index v) -> double {
        int iv = maps.v_idx[v]; return (iv >= 0) ? x[static_cast<std::size_t>(iv)] : 0.0;
    };
    auto get_ue = [&](Edge_index e) -> double {
        int ie = maps.e_idx[e]; return (ie >= 0) ? x[static_cast<std::size_t>(ie)] : 0.0;
    };
    auto arc_len = [&](Halfedge_index h) -> double {
        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(mesh.number_of_faces(), false);

    Face_index     f0 = *mesh.faces().begin();
    Halfedge_index h0 = mesh.halfedge(f0);
    Halfedge_index h1 = mesh.next(h0);
    Halfedge_index h2 = mesh.next(h1);
    Vertex_index   vA = mesh.source(h0);
    Vertex_index   vB = mesh.source(h1);
    Vertex_index   vC = mesh.source(h2);
    double lAB = arc_len(h0);
    double lCA = arc_len(h2);
    double lBC = arc_len(h1);

    result.pos[vA.idx()] = Eigen::Vector3d(0.0, 0.0, 1.0);
    result.pos[vB.idx()] = Eigen::Vector3d(std::sin(lAB), 0.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;
    face_placed[f0.idx()] = true;

    std::queue<Halfedge_index> q;
    auto enqueue = [&](Face_index f) {
        for (Halfedge_index h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) {
            Halfedge_index h_opp = mesh.opposite(h);
            if (!mesh.is_border(h_opp) && !face_placed[mesh.face(h_opp).idx()])
                q.push(h_opp);
        }
    };
    enqueue(f0);

    while (!q.empty()) {
        Halfedge_index h = q.front(); q.pop();
        Face_index f = mesh.face(h);
        if (face_placed[f.idx()]) continue;

        Vertex_index v_src = mesh.source(h);
        Vertex_index v_tgt = mesh.target(h);
        Vertex_index v_new = mesh.target(mesh.next(h));

        Eigen::Vector3d p = detail::trilaterate_sph(
            result.pos[v_src.idx()], result.pos[v_tgt.idx()],
            arc_len(mesh.prev(h)), arc_len(mesh.next(h)));

        if (!vertex_placed[v_new.idx()]) {
            result.pos[v_new.idx()]    = p;
            vertex_placed[v_new.idx()] = true;
        } else {
            result.has_seam = true;
        }
        face_placed[f.idx()] = true;
        enqueue(f);
    }

    result.success = true;
    if (normalise) normalise_spherical(result);
    return result;
}

// ── HyperIdeal layout (Poincaré disk) — exact trilateration ──────────────────
//
// Hyperbolic edge distance:
//   cosh(d_ij) = cosh(b_i + a_ij/2)·cosh(b_j + a_ij/2) − sinh(b_i)·sinh(b_j)
//
// Trilateration via Möbius map + hyperbolic law of cosines (exact, Phase 6).
//
// Optional CutGraph and HolonomyData for closed meshes.
// ─────────────────────────────────────────────────────────────────────────────
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();
    Layout2D result;
    result.uv.assign(nv, Eigen::Vector2d::Zero());
    if (mesh.number_of_faces() == 0) return result;

    auto get_b = [&](Vertex_index v) -> double {
        int iv = maps.v_idx[v]; return (iv >= 0) ? x[static_cast<std::size_t>(iv)] : 0.0;
    };
    auto get_a = [&](Edge_index e) -> double {
        int ie = maps.e_idx[e]; return (ie >= 0) ? x[static_cast<std::size_t>(ie)] : 0.0;
    };

    // Hyperbolic distance between two adjacent vertices
    auto hyp_dist = [&](Halfedge_index h) -> double {
        Vertex_index vi = mesh.source(h), vj = mesh.target(h);
        Edge_index   e  = mesh.edge(h);
        double bi = get_b(vi), bj = get_b(vj), a = get_a(e);
        double half_a = a * 0.5;
        double ch = std::cosh(bi + half_a) * std::cosh(bj + half_a)
                  - std::sinh(bi) * std::sinh(bj);
        return std::acosh(std::max(1.0, ch));
    };

    std::vector<bool> vertex_placed(nv, false);
    std::vector<bool> face_placed(mesh.number_of_faces(), false);

    // ── Place root face: vA at origin, vB on positive real axis ──────────────
    Face_index     f0 = *mesh.faces().begin();
    Halfedge_index h0 = mesh.halfedge(f0);
    Halfedge_index h1 = mesh.next(h0);
    Halfedge_index h2 = mesh.next(h1);
    Vertex_index   vA = mesh.source(h0);
    Vertex_index   vB = mesh.source(h1);
    Vertex_index   vC = mesh.source(h2);
    double dAB = hyp_dist(h0);
    double dCA = hyp_dist(h2);
    double dBC = hyp_dist(h1);

    result.uv[vA.idx()] = Eigen::Vector2d(0.0, 0.0);
    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;
    face_placed[f0.idx()] = true;

    // ── BFS ───────────────────────────────────────────────────────────────────
    std::queue<Halfedge_index> q;
    auto enqueue = [&](Face_index f) {
        for (Halfedge_index h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh)) {
            Halfedge_index h_opp = mesh.opposite(h);
            if (mesh.is_border(h_opp)) continue;
            if (cut && cut->is_cut(mesh.edge(h))) continue;
            Face_index f_adj = mesh.face(h_opp);
            if (!face_placed[f_adj.idx()])
                q.push(h_opp);
        }
    };
    enqueue(f0);

    while (!q.empty()) {
        Halfedge_index h = q.front(); q.pop();
        Face_index f = mesh.face(h);
        if (face_placed[f.idx()]) continue;

        Vertex_index v_src = mesh.source(h);
        Vertex_index v_tgt = mesh.target(h);
        Vertex_index v_new = mesh.target(mesh.next(h));

        double D  = hyp_dist(h);               // distance v_src → v_tgt
        double da = hyp_dist(mesh.prev(h));     // distance v_new → v_src
        double db = hyp_dist(mesh.next(h));     // distance v_tgt → v_new

        Eigen::Vector2d p = detail::trilaterate_hyp(
            result.uv[v_src.idx()], result.uv[v_tgt.idx()], D, da, db);

        if (!vertex_placed[v_new.idx()]) {
            result.uv[v_new.idx()]    = p;
            vertex_placed[v_new.idx()] = true;
        } else {
            result.has_seam = true;
        }
        face_placed[f.idx()] = true;
        enqueue(f);
    }

    result.success = true;

    // ── Holonomy ──────────────────────────────────────────────────────────────
    if (cut && holonomy) {
        holonomy->cut_edge_indices = cut->cut_edge_indices;
        holonomy->translations.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 h_opp  = mesh.opposite(mesh.halfedge(e));
            Halfedge_index h_cross = mesh.is_border(h_opp)
                                   ? mesh.halfedge(e) : h_opp;
            if (mesh.is_border(h_cross)) {
                holonomy->translations.push_back(Eigen::Vector2d::Zero());
                continue;
            }
            Vertex_index v_src = mesh.source(h_cross);
            Vertex_index v_tgt = mesh.target(h_cross);
            Vertex_index v_new = mesh.target(mesh.next(h_cross));
            if (!vertex_placed[v_new.idx()]) {
                holonomy->translations.push_back(Eigen::Vector2d::Zero());
                continue;
            }
            double D  = hyp_dist(h_cross);
            double da = hyp_dist(mesh.prev(h_cross));
            double db = hyp_dist(mesh.next(h_cross));
            Eigen::Vector2d p_tri = detail::trilaterate_hyp(
                result.uv[v_src.idx()], result.uv[v_tgt.idx()], D, da, db);
            holonomy->translations.push_back(p_tri - result.uv[v_new.idx()]);
        }
    }

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

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

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";
    }
}

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