// test_pipeline.cpp
//
// Phase 4c — End-to-end pipeline tests and library-user examples.
//
// These tests exercise the full conformallab++ pipeline as a user would:
//
//   1. Build (or load) a mesh
//   2. Set up maps and assign DOFs
//   3. Configure target angles / targets
//   4. Solve with Newton
//   5. Inspect / export the result
//
// Each test mirrors a realistic usage scenario documented in the README.
//
// Tests:
//   1. Pipeline_Euclidean_TriangleToEquilibrium
//        Read mesh → setup Euclidean maps → solve → verify convergence
//   2. Pipeline_Spherical_TetrahedronToEquilibrium
//        Setup spherical tetrahedron → solve → verify angles sum to 4π
//   3. Pipeline_HyperIdeal_TriangleRoundTrip
//        Build triangle → setup HyperIdeal → solve → verify G ≈ 0
//   4. Pipeline_MeshIO_SolveAndExport
//        Build mesh → solve → write OFF → reload → verify vertex count intact
//   5. Pipeline_AllThreeGeometries_SameTopology
//        Same quad-strip mesh solved under all three geometries: all converge

#include "conformal_mesh.hpp"
#include "mesh_builder.hpp"
#include "mesh_io.hpp"
#include "euclidean_functional.hpp"
#include "spherical_functional.hpp"
#include "hyper_ideal_functional.hpp"
#include "newton_solver.hpp"
#include <gtest/gtest.h>
#include <cmath>
#include <vector>
#include <filesystem>

using namespace conformallab;

// ────────────────────────────────────────────────────────────────────────────
// Shared helpers
// ────────────────────────────────────────────────────────────────────────────

static void pin_first_vertex_euclidean(ConformalMesh& mesh, EuclideanMaps& maps, int& n)
{
    auto vit = mesh.vertices().begin();
    Vertex_index v0 = *vit++;
    maps.v_idx[v0] = -1;
    int idx = 0;
    for (; vit != mesh.vertices().end(); ++vit)
        maps.v_idx[*vit] = idx++;
    n = idx;
}

static void set_natural_euclidean_theta(ConformalMesh& mesh, EuclideanMaps& maps, int n)
{
    std::vector<double> x0(static_cast<std::size_t>(n), 0.0);
    auto G = euclidean_gradient(mesh, x0, maps);
    for (auto v : mesh.vertices()) {
        int iv = maps.v_idx[v];
        if (iv < 0) continue;
        maps.theta_v[v] -= G[static_cast<std::size_t>(iv)];
    }
}

static std::vector<double> set_natural_hyper_ideal_targets(
    ConformalMesh& mesh, HyperIdealMaps& maps, int n,
    double b_base = 1.0, double a_base = 0.5)
{
    const auto sz = static_cast<std::size_t>(n);
    std::vector<double> xbase(sz, 0.0);
    for (auto v : mesh.vertices()) {
        int iv = maps.v_idx[v];
        if (iv >= 0) xbase[static_cast<std::size_t>(iv)] = b_base;
    }
    for (auto e : mesh.edges()) {
        int ie = maps.e_idx[e];
        if (ie >= 0) xbase[static_cast<std::size_t>(ie)] = a_base;
    }
    auto G = evaluate_hyper_ideal(mesh, xbase, maps, false).gradient;
    for (auto v : mesh.vertices()) {
        int iv = maps.v_idx[v];
        if (iv < 0) continue;
        maps.theta_v[v] += G[static_cast<std::size_t>(iv)];
    }
    for (auto e : mesh.edges()) {
        int ie = maps.e_idx[e];
        if (ie < 0) continue;
        maps.theta_e[e] += G[static_cast<std::size_t>(ie)];
    }
    return xbase;
}

// ════════════════════════════════════════════════════════════════════════════
// Test 1 — Euclidean full pipeline: triangle → equilibrium
//
// Simulates a user doing:
//   auto mesh = make_triangle();
//   auto maps = setup_euclidean_maps(mesh);
//   compute_euclidean_lambda0_from_mesh(mesh, maps);
//   // … set target angles and DOF indices …
//   auto result = newton_euclidean(mesh, x0, maps);
//   assert(result.converged);
// ════════════════════════════════════════════════════════════════════════════

TEST(Pipeline, Euclidean_TriangleToEquilibrium)
{
    // ── Step 1: build mesh ────────────────────────────────────────────────
    auto mesh = make_triangle();

    // ── Step 2: set up maps ───────────────────────────────────────────────
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);

    // ── Step 3: assign DOFs (pin v0) ──────────────────────────────────────
    int n = 0;
    pin_first_vertex_euclidean(mesh, maps, n);
    ASSERT_EQ(n, 2);

    // ── Step 4: choose natural target angles → x* = 0 ────────────────────
    set_natural_euclidean_theta(mesh, maps, n);

    // ── Step 5: solve ─────────────────────────────────────────────────────
    std::vector<double> x0(static_cast<std::size_t>(n), -0.1);
    auto result = newton_euclidean(mesh, x0, maps);

    // ── Step 6: verify ────────────────────────────────────────────────────
    EXPECT_TRUE(result.converged)
        << "Euclidean pipeline: triangle should converge; "
           "grad_inf_norm = " << result.grad_inf_norm;
    EXPECT_LT(result.grad_inf_norm, 1e-8);
    EXPECT_EQ(static_cast<int>(result.x.size()), n);
}

// ════════════════════════════════════════════════════════════════════════════
// Test 2 — Spherical full pipeline: tetrahedron → equilibrium
//
// The spherical tetrahedron equilibrium x* = 0 is built into the maps.
// After solving, the total angle defect Σ(Θ_v − Σα_v) should be ≈ 0.
// ════════════════════════════════════════════════════════════════════════════

TEST(Pipeline, Spherical_TetrahedronToEquilibrium)
{
    // ── Steps 1–3 ─────────────────────────────────────────────────────────
    auto mesh = make_spherical_tetrahedron();
    auto maps = setup_spherical_maps(mesh);
    compute_lambda0_from_mesh(mesh, maps);
    int n = assign_vertex_dof_indices(mesh, maps);

    // ── Step 4: solve ─────────────────────────────────────────────────────
    std::vector<double> x0(static_cast<std::size_t>(n), -0.2);
    auto result = newton_spherical(mesh, x0, maps);

    // ── Step 5: verify convergence ────────────────────────────────────────
    EXPECT_TRUE(result.converged)
        << "Spherical pipeline: tetrahedron should converge; "
           "grad_inf_norm = " << result.grad_inf_norm;
    EXPECT_LT(result.grad_inf_norm, 1e-8);

    // ── Step 6: verify geometric invariant — total angle defect ≈ 0 ──────
    auto G_final = spherical_gradient(mesh, result.x, maps);
    double total_defect = 0.0;
    for (double gv : G_final) total_defect += gv;
    EXPECT_NEAR(total_defect, 0.0, 1e-7)
        << "Spherical: total angle defect should vanish at equilibrium";
}

// ════════════════════════════════════════════════════════════════════════════
// Test 3 — HyperIdeal full pipeline: triangle → equilibrium
// ════════════════════════════════════════════════════════════════════════════

TEST(Pipeline, HyperIdeal_TriangleRoundTrip)
{
    // ── Steps 1–3 ─────────────────────────────────────────────────────────
    auto mesh = make_triangle();
    auto maps = setup_hyper_ideal_maps(mesh);
    int n = assign_all_dof_indices(mesh, maps);

    // ── Step 4: natural targets (equilibrium at b=1.0, a=0.5) ────────────
    auto xbase = set_natural_hyper_ideal_targets(mesh, maps, n);

    // Verify: gradient at xbase must be ≈ 0 before solving
    auto G_at_base = evaluate_hyper_ideal(mesh, xbase, maps, false).gradient;
    double max_g = 0.0;
    for (double v : G_at_base) max_g = std::max(max_g, std::abs(v));
    ASSERT_LT(max_g, 1e-10) << "Natural target setup: gradient at base should be ~0";

    // ── Step 5: perturb and solve ─────────────────────────────────────────
    std::vector<double> x0 = xbase;
    for (auto& v : x0) v += 0.3;

    auto result = newton_hyper_ideal(mesh, x0, maps);

    // ── Step 6: verify ────────────────────────────────────────────────────
    EXPECT_TRUE(result.converged)
        << "HyperIdeal pipeline: triangle should converge; "
           "grad_inf_norm = " << result.grad_inf_norm;
    EXPECT_LT(result.grad_inf_norm, 1e-8);

    // Solution should be close to xbase (same equilibrium)
    for (int i = 0; i < n; ++i) {
        EXPECT_NEAR(result.x[static_cast<std::size_t>(i)],
                    xbase[static_cast<std::size_t>(i)], 1e-6)
            << "DOF " << i << " should recover the equilibrium value";
    }
}

// ════════════════════════════════════════════════════════════════════════════
// Test 4 — Mesh I/O in the pipeline: solve → write → reload → check
//
// Demonstrates: compute a conformal factor on a mesh, write it to OFF, reload
// and check that the mesh topology is preserved.
// ════════════════════════════════════════════════════════════════════════════

TEST(Pipeline, MeshIO_SolveAndExport)
{
    // ── Build and solve ───────────────────────────────────────────────────
    auto mesh = make_quad_strip();
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);

    int n = 0;
    pin_first_vertex_euclidean(mesh, maps, n);
    set_natural_euclidean_theta(mesh, maps, n);

    std::vector<double> x0(static_cast<std::size_t>(n), -0.1);
    auto result = newton_euclidean(mesh, x0, maps);
    ASSERT_TRUE(result.converged) << "Solver must converge before export test";

    // ── Write mesh ────────────────────────────────────────────────────────
    const std::string tmp_path = "/tmp/conformallab_pipeline_test.off";
    ASSERT_NO_THROW(save_mesh(tmp_path, mesh));
    ASSERT_TRUE(std::filesystem::exists(tmp_path));

    // ── Reload and verify topology ────────────────────────────────────────
    ConformalMesh mesh2;
    ASSERT_NO_THROW(mesh2 = load_mesh(tmp_path));

    EXPECT_EQ(mesh2.number_of_vertices(), mesh.number_of_vertices())
        << "Vertex count must survive OFF round-trip";
    EXPECT_EQ(mesh2.number_of_faces(), mesh.number_of_faces())
        << "Face count must survive OFF round-trip";

    std::filesystem::remove(tmp_path);
}

// ════════════════════════════════════════════════════════════════════════════
// Test 5 — All three geometries, same quad-strip topology
//
// Validates that the solver infrastructure works uniformly: the same mesh
// topology is solvable under Euclidean, Spherical, and HyperIdeal geometries.
// ════════════════════════════════════════════════════════════════════════════

TEST(Pipeline, AllThreeGeometries_QuadStrip)
{
    // ── Euclidean ──────────────────────────────────────────────────────────
    {
        auto mesh = make_quad_strip();
        auto maps = setup_euclidean_maps(mesh);
        compute_euclidean_lambda0_from_mesh(mesh, maps);
        int n = 0;
        pin_first_vertex_euclidean(mesh, maps, n);
        set_natural_euclidean_theta(mesh, maps, n);
        std::vector<double> x0(static_cast<std::size_t>(n), -0.1);
        auto res = newton_euclidean(mesh, x0, maps);
        EXPECT_TRUE(res.converged) << "Euclidean: quad strip should converge";
    }

    // ── HyperIdeal ────────────────────────────────────────────────────────
    {
        auto mesh = make_quad_strip();
        auto maps = setup_hyper_ideal_maps(mesh);
        int n = assign_all_dof_indices(mesh, maps);
        auto xbase = set_natural_hyper_ideal_targets(mesh, maps, n);
        std::vector<double> x0 = xbase;
        for (auto& v : x0) v += 0.1;
        auto res = newton_hyper_ideal(mesh, x0, maps);
        EXPECT_TRUE(res.converged) << "HyperIdeal: quad strip should converge";
    }

    // ── Spherical (tetrahedron: smallest closed mesh with all vertices free) ─
    {
        auto mesh = make_spherical_tetrahedron();
        auto maps = setup_spherical_maps(mesh);
        compute_lambda0_from_mesh(mesh, maps);
        int n = assign_vertex_dof_indices(mesh, maps);
        std::vector<double> x0(static_cast<std::size_t>(n), -0.1);
        auto res = newton_spherical(mesh, x0, maps);
        EXPECT_TRUE(res.converged) << "Spherical: tetrahedron should converge";
    }
}