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

// test_euclidean_functional.cpp
//
// Phase 3d — EuclideanCyclicFunctional ported to ConformalMesh.
//
// Corresponds to de.varylab.discreteconformal.functional.EuclideanCyclicFunctionalTest.
//
// Test map (Java → C++)
// ──────────────────────
//   testHessian (Ignored)       → GradientCheck_Hessian           (ported)
//   testGradient…Triangle       → GradientCheck_TriangleVertex     (ported)
//   testGradient…QuadStrip      → GradientCheck_QuadStripVertex    (ported)
//   testGradient…Tetrahedron    → GradientCheck_TetrahedronVertex  (ported)
//   testGradient…AllDofs        → GradientCheck_TetrahedronAllDofs (ported)
//   testFunctionalAtNaNValue    → AnglesFiniteAtKnownPoint          (ported)
//
// Energy model
// ────────────
//   Uses the Schläfli path integral E(x) = ∫₀¹⟨G(tx),x⟩dt (10-point GL).
//   The gradient check verifies G is curl-free.

#include "conformal_mesh.hpp"
#include "mesh_builder.hpp"
#include "euclidean_geometry.hpp"
#include "euclidean_functional.hpp"
#include "euclidean_hessian.hpp"
#include "clausen.hpp"
#include <gtest/gtest.h>
#include <cmath>
#include <vector>

using namespace conformallab;

// ════════════════════════════════════════════════════════════════════════════
// Cross-module Hessian check: euclidean_gradient() ↔ euclidean_hessian()
//
// Java @Ignore reason: "no Hessian implemented yet" — the Java functional
// test was written before the Hessian existed.  In C++ the analytic
// cotangent-Laplace Hessian (euclidean_hessian.hpp, Phase 3f) is complete.
//
// This test verifies cross-module consistency:
//   H[i,j] ≈ (G_i(x+ε·eⱼ) − G_i(x−ε·eⱼ)) / (2ε)
// using the gradient from euclidean_functional.hpp and the Hessian from
// euclidean_hessian.hpp.  A bug in DOF-index mapping or sign convention
// that affects both modules independently would only be caught here.
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, GradientCheck_Hessian)
{
    auto mesh = make_triangle();
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);
    int n = assign_euclidean_vertex_dof_indices(mesh, maps);

    std::vector<double> x(static_cast<std::size_t>(n), -0.1);

    // hessian_check_euclidean: H[i,j] ≈ FD(G)[i,j] using euclidean_gradient()
    EXPECT_TRUE(hessian_check_euclidean(mesh, x, maps))
        << "Cross-module: euclidean_gradient() and euclidean_hessian() are inconsistent";
}

// ════════════════════════════════════════════════════════════════════════════
// Angle formula: equilateral triangle → all angles = π/3
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, EquilateralTriangleAnglesArePiOver3)
{
    // All sides equal: l = 1.0, log-length = 0.
    auto fa = euclidean_angles(0.0, 0.0, 0.0);

    ASSERT_TRUE(fa.valid) << "Equilateral triangle must be valid";
    constexpr double PI_3 = 3.14159265358979323846 / 3.0;
    EXPECT_NEAR(fa.alpha1, PI_3, 1e-12);
    EXPECT_NEAR(fa.alpha2, PI_3, 1e-12);
    EXPECT_NEAR(fa.alpha3, PI_3, 1e-12);
}

// ════════════════════════════════════════════════════════════════════════════
// Angle formula: right isosceles triangle (legs 1, hypotenuse √2)
//
// For the 45-45-90 triangle: angles are π/4, π/4, π/2.
// From make_triangle default (v0=(0,0), v1=(1,0), v2=(0,1)):
//   e01: l=1, λ°=0
//   e12: l=√2, λ°=log(2)
//   e02: l=1, λ°=0
// Angle at v0 (opposite e12) = π/2.
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, RightIsoscelesTriangleAnglesCorrect)
{
    const double log2 = std::log(2.0);
    // lam12 = 0 (v0-v1, length 1), lam23 = log(2) (v1-v2, length √2), lam31 = 0 (v2-v0, length 1)
    // v1 = v0 in our ordering → remap: l01=1, l12=√2, l20=1
    // Using euclidean_angles(lam_v1v2, lam_v2v3, lam_v3v1):
    //   v1=(0,0), v2=(1,0), v3=(0,1)
    //   lam12 = log(1²) = 0, lam23 = log(√2 ²) = log2, lam31 = log(1²) = 0
    auto fa = euclidean_angles(0.0, log2, 0.0);

    ASSERT_TRUE(fa.valid);
    constexpr double PI   = 3.14159265358979323846;
    // v1=(0,0) is at the right-angle corner (opposite the hypotenuse l23=√2) → α1 = 90°.
    // v2=(1,0) and v3=(0,1) are the 45° corners (each opposite a leg of length 1).
    EXPECT_NEAR(fa.alpha1, PI / 2.0, 1e-12);   // angle at v1 (opposite l23=√2): 90°
    EXPECT_NEAR(fa.alpha2, PI / 4.0, 1e-12);   // angle at v2 (opposite l31=1):  45°
    EXPECT_NEAR(fa.alpha3, PI / 4.0, 1e-12);   // angle at v3 (opposite l12=1):  45°
}

// ════════════════════════════════════════════════════════════════════════════
// Angle sum = π for any valid Euclidean triangle
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, AngleSumEqualsPi)
{
    // Scalene triangle with log-lengths (0, 0.5, -0.3).
    auto fa = euclidean_angles(0.0, 0.5, -0.3);

    ASSERT_TRUE(fa.valid);
    constexpr double PI = 3.14159265358979323846;
    EXPECT_NEAR(fa.alpha1 + fa.alpha2 + fa.alpha3, PI, 1e-12);
}

// ════════════════════════════════════════════════════════════════════════════
// Degenerate triangle → valid = false
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, DegenerateTriangleReturnsFalse)
{
    // l12 = l23 = 1, l31 = 3 → violates triangle inequality.
    auto fa = euclidean_angles_from_lengths(1.0, 1.0, 3.0);
    EXPECT_FALSE(fa.valid);
}

// ════════════════════════════════════════════════════════════════════════════
// Gradient check: default right-isosceles triangle, vertex DOFs only
//
// Mirrors Java testGradient…SingleTriangle.
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, GradientCheck_TriangleVertex)
{
    auto mesh = make_triangle();   // (0,0)–(1,0)–(0,1)
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);
    int n = assign_euclidean_vertex_dof_indices(mesh, maps);

    // Small uniform conformal perturbation.
    std::vector<double> x(static_cast<std::size_t>(n), -0.1);

    EXPECT_TRUE(gradient_check_euclidean(mesh, x, maps))
        << "Gradient check failed on right-isosceles triangle (vertex DOFs)";
}

// ════════════════════════════════════════════════════════════════════════════
// Gradient check: quad strip (2 triangles, 1 interior edge), vertex DOFs only
//
// Mirrors Java testGradient…QuadStrip / testGradientInExtendedDomain.
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, GradientCheck_QuadStripVertex)
{
    auto mesh = make_quad_strip();
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);
    int n = assign_euclidean_vertex_dof_indices(mesh, maps);

    std::vector<double> x(static_cast<std::size_t>(n), -0.2);

    EXPECT_TRUE(gradient_check_euclidean(mesh, x, maps))
        << "Gradient check failed on quad strip (vertex DOFs)";
}

// ════════════════════════════════════════════════════════════════════════════
// Gradient check: regular tetrahedron, vertex DOFs only
//
// Closed surface (4 faces, 4 vertices, 6 interior edges).
// Exercises per-vertex angle-sum accumulation on multiple faces.
// Mirrors Java testGradient…Tetrahedron / testGradientWithHyperIdeal…
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, GradientCheck_TetrahedronVertex)
{
    auto mesh = make_tetrahedron();
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);
    int n = assign_euclidean_vertex_dof_indices(mesh, maps);

    std::vector<double> x(static_cast<std::size_t>(n), -0.15);

    EXPECT_TRUE(gradient_check_euclidean(mesh, x, maps))
        << "Gradient check failed on regular tetrahedron (vertex DOFs)";
}

// ════════════════════════════════════════════════════════════════════════════
// Gradient check: tetrahedron, all DOFs (vertex + edge)
//
// Exercises the edge-gradient branch G_e = α_opp⁺ + α_opp⁻ − π.
// Mirrors Java testGradientWithHyperellipticCurve.
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, GradientCheck_TetrahedronAllDofs)
{
    auto mesh = make_tetrahedron();
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);
    int n = assign_euclidean_all_dof_indices(mesh, maps);

    // 4 vertex DOFs + 6 edge DOFs = 10 total.
    std::vector<double> x(static_cast<std::size_t>(n), 0.0);
    // Set vertex DOFs slightly negative to keep triangles non-degenerate.
    for (int i = 0; i < 4; ++i)
        x[static_cast<std::size_t>(i)] = -0.15;

    EXPECT_TRUE(gradient_check_euclidean(mesh, x, maps))
        << "Gradient check failed on regular tetrahedron (all DOFs)";
}

// ════════════════════════════════════════════════════════════════════════════
// Angles are finite at a known interior point
//
// Mirrors Java testFunctionalAtNaNValue: stress-test the angle formula with
// large negative conformal factors (compressed triangle) to ensure no NaN/Inf.
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, AnglesFiniteAtKnownPoint)
{
    auto mesh = make_tetrahedron();
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);
    int n = assign_euclidean_vertex_dof_indices(mesh, maps);

    // Very compressed: u_i = -3 (all sides shrunk by exp(-3) ≈ 0.05).
    // Triangle stays well-formed (equilateral shrinks uniformly).
    std::vector<double> x(static_cast<std::size_t>(n), -3.0);
    auto G = euclidean_gradient(mesh, x, maps);

    for (std::size_t i = 0; i < G.size(); ++i) {
        EXPECT_FALSE(std::isnan(G[i])) << "Gradient component " << i << " is NaN";
        EXPECT_FALSE(std::isinf(G[i])) << "Gradient component " << i << " is Inf";
    }
}

// ════════════════════════════════════════════════════════════════════════════
// Gradient check: fan of 5 flat triangles, vertex DOFs only
//
// High-valence central vertex: exercises per-vertex angle accumulation
// across 5 incident faces.
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, GradientCheck_Fan5Vertex)
{
    auto mesh = make_fan(5);
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);
    int n = assign_euclidean_vertex_dof_indices(mesh, maps);

    std::vector<double> x(static_cast<std::size_t>(n), -0.05);

    EXPECT_TRUE(gradient_check_euclidean(mesh, x, maps))
        << "Gradient check failed on flat fan-5 mesh";
}

// ════════════════════════════════════════════════════════════════════════════
// Gradient check: mixed pinned/variable vertices
//
// Pins the first vertex (u_v0 = 0 fixed), lets the rest be variable.
// Verifies that the gradient accumulator skips pinned vertices correctly.
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanFunctional, GradientCheck_MixedPinnedVertices)
{
    auto mesh = make_quad_strip();
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);

    // Manually pin v0; assign v1, v2, v3 as DOFs 0, 1, 2.
    auto vit = mesh.vertices().begin();
    Vertex_index v0 = *vit++;
    Vertex_index v1 = *vit++;
    Vertex_index v2 = *vit++;
    Vertex_index v3 = *vit;

    maps.v_idx[v0] = -1;   // pinned
    maps.v_idx[v1] = 0;
    maps.v_idx[v2] = 1;
    maps.v_idx[v3] = 2;

    std::vector<double> x = {-0.1, -0.3, -0.2};

    EXPECT_TRUE(gradient_check_euclidean(mesh, x, maps))
        << "Gradient check failed for mixed pinned/variable vertices";
}

// ─────────────────────────────────────────────────────────────────────────────
// Golden-value oracle — pin the Euclidean angle formula and the 2·Л(α) energy
// term bit-for-bit against the upstream Java reference (EuclideanCyclicFunctional
// .triangleEnergyAndAlphas, lines 341-361), captured by running the compiled Java
// library (openjdk 17) with the real de.varylab…Clausen.Л on these exact edge
// lengths.  Companion to HyperIdealGoldenJava and SphericalGoldenJava: locks the
// absolute angle/energy values against an independent implementation, catching
// silent index/sign drift the curl-free path-integral gradient check cannot see.
//
// To regenerate: /tmp/oracle/EucOracle.java.  Values are Java printf %.17g.
// ─────────────────────────────────────────────────────────────────────────────
TEST(EuclideanGoldenJava, AngleAndLobachevskyEnergyFromLengths)
{
    auto check = [](double l12, double l23, double l31,
                    double a1_g, double a2_g, double a3_g, double L_g) {
        auto fa = euclidean_angles_from_lengths(l12, l23, l31);
        EXPECT_TRUE(fa.valid);
        EXPECT_NEAR(fa.alpha1, a1_g, 1e-12);
        EXPECT_NEAR(fa.alpha2, a2_g, 1e-12);
        EXPECT_NEAR(fa.alpha3, a3_g, 1e-12);
        const double Lterm = 2.0 * Lobachevsky(fa.alpha1)
                           + 2.0 * Lobachevsky(fa.alpha2)
                           + 2.0 * Lobachevsky(fa.alpha3);
        EXPECT_NEAR(Lterm, L_g, 1e-12);
    };
    check(1.0, 1.2, 0.9,
          1.3637649752769678, 0.82416964552030680, 0.95365803279251860,
          1.9456273836230942);
    check(1.0, 1.0, 1.0,
          1.0471975511965979, 1.0471975511965979, 1.0471975511965979,
          2.0298832128193070);
}

// ─────────────────────────────────────────────────────────────────────────────
// FULL-MESH golden oracle — the strongest cross-check: drives the REAL upstream
// EuclideanCyclicFunctional (openjdk 17) on a tetrahedron loaded from a shared
// OBJ (identical topology + geometry to make_tetrahedron()), and pins BOTH the
// per-vertex gradient G_v = Θ_v − Σα AND the energy difference ΔE = E(x) − E(0)
// bit-for-bit against it.
//
// This closes audit missing-test item 5 (full-mesh energy + gradient at a known
// x).  It is genuinely independent of the C++ implementation in two ways:
//   • the gradient is the upstream library's own analytic gradient (not an FD
//     check, which only proves curl-freeness of the C++ self-consistent energy);
//   • the energy is Java's CLOSED-FORM functional value, whereas C++ computes it
//     as a Gauss-Legendre PATH INTEGRAL of its gradient — two different methods
//     that must agree on ΔE (the initialEnergy constant and the φ·λ⁰ term cancel
//     in the difference, and there are no edge DOFs here).
//
// Setup parity (verified against UnwrapUtility.prepareInvariantDataEuclidean):
//   closed mesh, ALL 4 vertices variable (no pin), Θ_v = 2π, no edge DOFs,
//   λ°_e = 2·log(|p_i − p_j|), per-vertex u(P) = 0.10·X − 0.07·Y + 0.13·Z.
//
// To regenerate: /tmp/oracle/{tet.obj,EucMeshOracle.java}.  Values are Java %.17g.
// ─────────────────────────────────────────────────────────────────────────────
TEST(EuclideanGoldenJava, FullMeshGradientAndEnergy_Tetrahedron)
{
    constexpr double TWO_PI = 2.0 * 3.14159265358979323846264338328;

    auto mesh = make_tetrahedron();   // same 4 vertices as /tmp/oracle/tet.obj
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);

    // All four vertices are DOFs (Java prepareInvariantData makes every interior
    // vertex variable on a closed mesh); Θ_v = 2π; no edge DOFs.
    int idx = 0;
    for (auto v : mesh.vertices()) {
        maps.v_idx[v]   = idx++;
        maps.theta_v[v] = TWO_PI;
    }

    auto u_of = [](const Point3& p) {
        return 0.10 * p.x() - 0.07 * p.y() + 0.13 * p.z();
    };

    std::vector<double> x(static_cast<std::size_t>(idx), 0.0);
    for (auto v : mesh.vertices())
        x[static_cast<std::size_t>(maps.v_idx[v])] = u_of(mesh.point(v));

    auto G = euclidean_gradient(mesh, x, maps);

    // Java golden gradients keyed by vertex position (order-independent lookup).
    struct GoldRow { double X, Y, Z, G; };
    const GoldRow gold[4] = {
        { 1,  1,  1, 3.5277511803984396},
        { 1, -1, -1, 3.2837611905358440},
        {-1,  1, -1, 2.3437485291237100},
        {-1, -1,  1, 3.4111097143011780},
    };
    for (auto v : mesh.vertices()) {
        const auto& p = mesh.point(v);
        const double g = G[static_cast<std::size_t>(maps.v_idx[v])];
        bool matched = false;
        for (const auto& row : gold) {
            if (std::abs(p.x() - row.X) < 1e-9 &&
                std::abs(p.y() - row.Y) < 1e-9 &&
                std::abs(p.z() - row.Z) < 1e-9) {
                EXPECT_NEAR(g, row.G, 1e-12)
                    << "gradient mismatch at (" << p.x() << "," << p.y()
                    << "," << p.z() << ")";
                matched = true;
                break;
            }
        }
        EXPECT_TRUE(matched) << "unexpected vertex position";
    }

    // ΔE = E(x) − E(0).  C++ path integral vs Java closed-form functional value.
    std::vector<double> x0(static_cast<std::size_t>(idx), 0.0);
    const double dE = euclidean_energy(mesh, x, maps)
                    - euclidean_energy(mesh, x0, maps);
    EXPECT_NEAR(dE, 0.15962619236187336, 1e-12);
}