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

// test_phase6.cpp
//
// Phase 6 — Tests for:
//   - gauss_bonnet.hpp  : euler_characteristic, genus, GB sum/rhs, check, enforce
//   - cut_graph.hpp     : compute_cut_graph (tree-cotree algorithm)
//   - layout.hpp Phase6 : exact hyperbolic trilateration math
//                         normalise_euclidean (centroid + PCA)
//                         normalise_hyperbolic (Möbius centering)

#include "conformal_mesh.hpp"
#include "mesh_builder.hpp"
#include "euclidean_functional.hpp"
#include "spherical_functional.hpp"
#include "hyper_ideal_functional.hpp"
#include "newton_solver.hpp"
#include "gauss_bonnet.hpp"
#include "cut_graph.hpp"
#include "layout.hpp"
#include <gtest/gtest.h>
#include <cmath>
#include <complex>
#include <vector>

using namespace conformallab;

// mesh_builder.hpp provides make_triangle(), make_tetrahedron(), make_quad_strip().

// ════════════════════════════════════════════════════════════════════════════
// GaussBonnet — topology helpers
// ════════════════════════════════════════════════════════════════════════════

TEST(GaussBonnet, EulerCharacteristic_Triangle)
{
    // V=3 E=3 F=1  →  χ = 1
    auto m = make_triangle();
    EXPECT_EQ(euler_characteristic(m), 1);
}

TEST(GaussBonnet, EulerCharacteristic_QuadStrip)
{
    // V=4 E=5 F=2  →  χ = 1
    auto m = make_quad_strip();
    EXPECT_EQ(euler_characteristic(m), 1);
}

TEST(GaussBonnet, EulerCharacteristic_Tetrahedron)
{
    // V=4 E=6 F=4  →  χ = 2
    auto m = make_tetrahedron();
    EXPECT_EQ(euler_characteristic(m), 2);
}

TEST(GaussBonnet, Genus_Tetrahedron_IsZero)
{
    auto m = make_tetrahedron();
    EXPECT_EQ(genus(m), 0);
}

TEST(GaussBonnet, RHS_Triangle)
{
    // χ=1  →  2π·χ = 2π
    auto m = make_triangle();
    EXPECT_NEAR(gauss_bonnet_rhs(m), 2.0 * M_PI, 1e-12);
}

TEST(GaussBonnet, RHS_Tetrahedron)
{
    // χ=2  →  2π·χ = 4π
    auto m = make_tetrahedron();
    EXPECT_NEAR(gauss_bonnet_rhs(m), 4.0 * M_PI, 1e-12);
}

// ════════════════════════════════════════════════════════════════════════════
// GaussBonnet — sum / deficit / check / enforce
// ════════════════════════════════════════════════════════════════════════════

TEST(GaussBonnet, DeficitIsNonZeroWithDefaultTheta)
{
    // Default theta_v = 2π at all V=4 vertices of quad-strip  →
    //   Σ(2π−2π) = 0,  rhs = 2π·1 = 2π  →  deficit = −2π ≠ 0.
    auto m = make_quad_strip();
    auto maps = setup_euclidean_maps(m);
    double def = gauss_bonnet_deficit(m, maps);
    EXPECT_GT(std::abs(def), 1e-6);
}

TEST(GaussBonnet, EnforceAdjustsTheta_Deficit_BecomesZero)
{
    auto m = make_quad_strip();
    auto maps = setup_euclidean_maps(m);
    // Set clearly wrong theta
    for (auto v : m.vertices()) maps.theta_v[v] = M_PI;
    EXPECT_GT(std::abs(gauss_bonnet_deficit(m, maps)), 1e-6);

    enforce_gauss_bonnet(m, maps);
    EXPECT_NEAR(gauss_bonnet_deficit(m, maps), 0.0, 1e-10);
}

TEST(GaussBonnet, CheckThrowsWhenViolated)
{
    auto m = make_triangle();
    auto maps = setup_euclidean_maps(m);
    // theta_v = 2π everywhere  →  sum = 0, rhs = 2π  →  |deficit| = 2π
    EXPECT_THROW(check_gauss_bonnet(m, maps), std::runtime_error);
}

TEST(GaussBonnet, CheckPassesAfterEnforce)
{
    auto m = make_triangle();
    auto maps = setup_euclidean_maps(m);
    enforce_gauss_bonnet(m, maps);
    EXPECT_NO_THROW(check_gauss_bonnet(m, maps));
}

TEST(GaussBonnet, SumAfterEnforce_EqualsRHS)
{
    auto m = make_tetrahedron();
    auto maps = setup_euclidean_maps(m);
    // theta_v starts at 2π everywhere; sum = 0, rhs = 4π
    enforce_gauss_bonnet(m, maps);
    double lhs = gauss_bonnet_sum(m, maps);
    double rhs = gauss_bonnet_rhs(m);
    EXPECT_NEAR(lhs, rhs, 1e-10);
}

TEST(GaussBonnet, ManuallySetAnalyticalTheta_PassesCheck)
{
    // Single triangle (χ=1): need Σ(2π−Θ_v) = 2π.
    // Set each of the 3 boundary vertices to Θ_v = 4π/3.
    // Then Σ(2π − 4π/3) = 3·(2π/3) = 2π. ✓
    auto m = make_triangle();
    auto maps = setup_euclidean_maps(m);
    for (auto v : m.vertices()) maps.theta_v[v] = 4.0 * M_PI / 3.0;
    EXPECT_NO_THROW(check_gauss_bonnet(m, maps));
}

// ════════════════════════════════════════════════════════════════════════════
// CutGraph — tree-cotree algorithm
// ════════════════════════════════════════════════════════════════════════════

TEST(CutGraph, Triangle_ZeroCutEdges)
{
    // Open disk → genus 0 → 0 cut edges.
    auto m = make_triangle();
    auto cg = compute_cut_graph(m);
    EXPECT_EQ(cg.cut_edge_indices.size(), 0u);
}

TEST(CutGraph, QuadStrip_ZeroCutEdges)
{
    // Open disk → genus 0 → 0 cut edges.
    auto m = make_quad_strip();
    auto cg = compute_cut_graph(m);
    EXPECT_EQ(cg.cut_edge_indices.size(), 0u);
}

TEST(CutGraph, Tetrahedron_ZeroCutEdges_ClosedSphere)
{
    // Closed genus-0 surface → 2g = 0 cut edges.
    auto m = make_tetrahedron();
    auto cg = compute_cut_graph(m);
    EXPECT_EQ(cg.cut_edge_indices.size(), 0u);
    EXPECT_EQ(cg.genus, 0);
}

TEST(CutGraph, FlagsVsIndicesConsistent)
{
    // Every index in cut_edge_indices has flag=true;
    // count of true flags == number of indices.
    auto m = make_tetrahedron();
    auto cg = compute_cut_graph(m);
    for (std::size_t idx : cg.cut_edge_indices)
        EXPECT_TRUE(cg.cut_edge_flags[idx]);

    std::size_t count = 0;
    for (bool b : cg.cut_edge_flags) count += b ? 1u : 0u;
    EXPECT_EQ(count, cg.cut_edge_indices.size());
}

TEST(CutGraph, IsCutMethodMatchesFlags)
{
    auto m = make_tetrahedron();
    auto cg = compute_cut_graph(m);
    for (auto e : m.edges()) {
        bool flag = cg.cut_edge_flags[static_cast<std::size_t>(e.idx())];
        EXPECT_EQ(cg.is_cut(e), flag);
    }
}

TEST(CutGraph, FlagsVectorHasCorrectSize)
{
    auto m = make_quad_strip();
    auto cg = compute_cut_graph(m);
    EXPECT_EQ(cg.cut_edge_flags.size(), m.number_of_edges());
}

// ════════════════════════════════════════════════════════════════════════════
// Exact hyperbolic trilateration — Möbius + law of cosines
// ════════════════════════════════════════════════════════════════════════════

namespace {

/// Poincaré-disk hyperbolic distance.
double hyp_dist_disk(Eigen::Vector2d a, Eigen::Vector2d b)
{
    double num = (a.x()-b.x())*(a.x()-b.x()) + (a.y()-b.y())*(a.y()-b.y());
    double da  = 1.0 - a.x()*a.x() - a.y()*a.y();
    double db_  = 1.0 - b.x()*b.x() - b.y()*b.y();
    double arg = 1.0 + 2.0*num/(da*db_);
    return std::acosh(std::max(1.0, arg));
}

/// Reproduce the exact trilateration formula from layout.hpp detail namespace.
Eigen::Vector2d trilaterate_hyp(Eigen::Vector2d pa, Eigen::Vector2d pb,
                                  double D, double da, double db)
{
    if (D < 1e-12 || da < 1e-12) return pa;
    if (std::sinh(da) * std::sinh(D) < 1e-14) return pa;
    using C = std::complex<double>;
    auto mobius_fwd = [](C z, C center) -> C {
        return (z - center) / (C(1.0) - std::conj(center) * z);
    };
    auto mobius_inv = [](C w, C center) -> C {
        return (w + center) / (C(1.0) + std::conj(center) * w);
    };
    C a(pa.x(), pa.y()), b(pb.x(), pb.y());
    C b_mapped = mobius_fwd(b, a);
    double theta_b = std::arg(b_mapped);
    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);
    C pc_origin = std::tanh(da * 0.5) * std::exp(C(0.0, theta_b + alpha));
    C pc_c = mobius_inv(pc_origin, a);
    return Eigen::Vector2d(pc_c.real(), pc_c.imag());
}

} // namespace

TEST(HyperbolicTrilateration, RecoveredDistances_Small)
{
    // pa at origin, pb at tanh(D/2) on x-axis.
    double D = 0.8, da = 0.5, db = 0.6;
    Eigen::Vector2d pa(0.0, 0.0);
    Eigen::Vector2d pb(std::tanh(D * 0.5), 0.0);

    Eigen::Vector2d pc = trilaterate_hyp(pa, pb, D, da, db);

    EXPECT_NEAR(hyp_dist_disk(pa, pb), D,  1e-10);
    EXPECT_NEAR(hyp_dist_disk(pa, pc), da, 1e-9);
    EXPECT_NEAR(hyp_dist_disk(pb, pc), db, 1e-9);
}

TEST(HyperbolicTrilateration, RecoveredDistances_Large)
{
    double D = 2.0, da = 1.5, db = 1.0;
    Eigen::Vector2d pa(0.0, 0.0);
    Eigen::Vector2d pb(std::tanh(D * 0.5), 0.0);

    Eigen::Vector2d pc = trilaterate_hyp(pa, pb, D, da, db);

    EXPECT_NEAR(hyp_dist_disk(pa, pc), da, 1e-8);
    EXPECT_NEAR(hyp_dist_disk(pb, pc), db, 1e-8);
}

TEST(HyperbolicTrilateration, ResultInsidePoincareDisk)
{
    double D = 1.2, da = 0.9, db = 0.7;
    Eigen::Vector2d pa(0.0, 0.0);
    Eigen::Vector2d pb(std::tanh(D * 0.5), 0.0);

    Eigen::Vector2d pc = trilaterate_hyp(pa, pb, D, da, db);
    EXPECT_LT(pc.norm(), 1.0);
}

TEST(HyperbolicTrilateration, OffOrigin_StillCorrect)
{
    // Place pa somewhere in the disk (not at origin) to test Möbius map.
    double D = 0.6, da = 0.4, db = 0.5;
    // pa at (0.3, 0.0) in Poincaré coords.
    Eigen::Vector2d pa(0.3, 0.0);
    // pb at distance D from pa — compute pb in Poincaré disk:
    // Möbius inverse: tanh(D/2) maps back from origin frame.
    using C = std::complex<double>;
    C a_c(pa.x(), pa.y());
    C pb_c = (std::tanh(D*0.5) + a_c) / (C(1.0) + std::conj(a_c) * std::tanh(D*0.5));
    Eigen::Vector2d pb(pb_c.real(), pb_c.imag());

    Eigen::Vector2d pc = trilaterate_hyp(pa, pb, D, da, db);
    EXPECT_LT(pc.norm(), 1.0);
    EXPECT_NEAR(hyp_dist_disk(pa, pc), da, 1e-8);
    EXPECT_NEAR(hyp_dist_disk(pb, pc), db, 1e-8);
}

// ════════════════════════════════════════════════════════════════════════════
// Normalisation — euclidean_layout (centroid) + hyper_ideal_layout (Möbius)
// ════════════════════════════════════════════════════════════════════════════

/// Build a solved euclidean layout for the quad-strip at the natural equilibrium.
static Layout2D make_euclidean_layout_normalised(bool normalise)
{
    auto mesh = make_quad_strip();
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);

    // Pin first vertex; assign free DOF indices to all others.
    auto vit = mesh.vertices().begin();
    maps.v_idx[*vit++] = -1;
    int idx = 0;
    for (; vit != mesh.vertices().end(); ++vit) maps.v_idx[*vit] = idx++;
    const int n = idx;

    std::vector<double> x0(static_cast<std::size_t>(n), 0.0);
    // Adjust theta so G(x=0) = 0 (natural equilibrium).
    auto G0 = euclidean_gradient(mesh, x0, maps);
    for (auto v : mesh.vertices()) {
        int iv = maps.v_idx[v];
        if (iv >= 0) maps.theta_v[v] -= G0[static_cast<std::size_t>(iv)];
    }

    auto res = newton_euclidean(mesh, x0, maps, 1e-10, 100);
    return euclidean_layout(mesh, res.x, maps, nullptr, nullptr, normalise);
}

TEST(Normalisation, Euclidean_CentroidAtOrigin)
{
    auto L = make_euclidean_layout_normalised(true);
    ASSERT_GT(L.uv.size(), 0u);

    double cx = 0.0, cy = 0.0;
    for (auto& p : L.uv) { cx += p.x(); cy += p.y(); }
    int n = static_cast<int>(L.uv.size());
    EXPECT_NEAR(cx / n, 0.0, 1e-9);
    EXPECT_NEAR(cy / n, 0.0, 1e-9);
}

TEST(Normalisation, Euclidean_NormalisePreservesEdgeLengthRatios)
{
    auto L_no  = make_euclidean_layout_normalised(false);
    auto L_yes = make_euclidean_layout_normalised(true);

    // Ratio of the lengths of the first two edges must be preserved.
    ASSERT_GE(L_no.uv.size(), 4u);
    // edges 0→1 and 0→2 (the two edges from vertex 0 in the quad-strip)
    double len0_no  = (L_no.uv[1]  - L_no.uv[0]).norm();
    double len1_no  = (L_no.uv[2]  - L_no.uv[0]).norm();
    double len0_yes = (L_yes.uv[1] - L_yes.uv[0]).norm();
    double len1_yes = (L_yes.uv[2] - L_yes.uv[0]).norm();

    if (len1_no > 1e-12 && len1_yes > 1e-12) {
        double ratio_no  = len0_no  / len1_no;
        double ratio_yes = len0_yes / len1_yes;
        EXPECT_NEAR(ratio_no, ratio_yes, 1e-9);
    }
}

/// Build a solved hyper-ideal layout for a triangle at natural equilibrium.
static Layout2D make_hyper_ideal_layout_normalised(bool normalise)
{
    auto mesh = make_triangle();
    auto maps = setup_hyper_ideal_maps(mesh);
    int n = assign_all_dof_indices(mesh, maps);

    std::vector<double> xbase(static_cast<std::size_t>(n));
    for (auto v : mesh.vertices()) {
        int iv = maps.v_idx[v]; if (iv >= 0) xbase[iv] = 1.0;
    }
    for (auto e : mesh.edges()) {
        int ie = maps.e_idx[e]; if (ie >= 0) xbase[ie] = 0.5;
    }
    auto G0 = evaluate_hyper_ideal(mesh, xbase, maps, false).gradient;
    for (auto v : mesh.vertices()) {
        int iv = maps.v_idx[v]; if (iv >= 0) maps.theta_v[v] += G0[iv];
    }
    for (auto e : mesh.edges()) {
        int ie = maps.e_idx[e]; if (ie >= 0) maps.theta_e[e] += G0[ie];
    }
    auto res = newton_hyper_ideal(mesh, xbase, maps, 1e-9, 100);
    return hyper_ideal_layout(mesh, res.x, maps, nullptr, nullptr, normalise);
}

TEST(Normalisation, HyperIdeal_PositionsInsidePoincareDisk)
{
    auto L = make_hyper_ideal_layout_normalised(true);
    ASSERT_GT(L.uv.size(), 0u);
    for (auto& p : L.uv) {
        EXPECT_LT(p.norm(), 1.0 + 1e-9)
            << "Vertex outside Poincaré disk: r=" << p.norm();
    }
}

TEST(Normalisation, HyperIdeal_CenteringReducesAverageRadius)
{
    // After Möbius centering the average Euclidean radius inside the disk
    // must be ≤ the unconstrained average.
    auto L_no  = make_hyper_ideal_layout_normalised(false);
    auto L_yes = make_hyper_ideal_layout_normalised(true);

    int n = static_cast<int>(L_no.uv.size());
    double r_no = 0.0, r_yes = 0.0;
    for (int i = 0; i < n; ++i) {
        r_no  += L_no.uv[i].norm();
        r_yes += L_yes.uv[i].norm();
    }
    EXPECT_LE(r_yes / n, r_no / n + 1e-9);
}