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

// test_euclidean_hessian.cpp
//
// Phase 3f — Euclidean cotangent-Laplace Hessian.
//
// The Hessian of the Euclidean discrete conformal energy is the well-known
// cotangent-Laplace operator (Pinkall–Polthier 1993, Springborn 2008).
//
// Tests:
//   1. Cotangent weights are analytically correct for simple triangles.
//   2. Hessian is symmetric.
//   3. Hessian has the null-space property H·1 = 0 (uniform-shift mode).
//   4. Hessian is positive semi-definite (all eigenvalues ≥ 0).
//   5. Finite-difference check H[i,j] ≈ (G_i(x+ε·eⱼ)−G_i(x−ε·eⱼ))/(2ε).
//
// All tests use meshes and maps built with Phase-3d infrastructure.

#include "conformal_mesh.hpp"
#include "mesh_builder.hpp"
#include "euclidean_hessian.hpp"
#include <gtest/gtest.h>
#include <Eigen/Dense>   // for dense conversion and eigenvalue solver
#include <cmath>
#include <vector>

using namespace conformallab;

// ════════════════════════════════════════════════════════════════════════════
// Cotangent weight: equilateral triangle → all cots = 1/√3 = cot(60°)
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanHessian, CotWeights_EquilateralTriangle)
{
    // Equilateral triangle with l = 1 (all log-lengths = 0).
    auto cw = euclidean_cot_weights(1.0, 1.0, 1.0);

    ASSERT_TRUE(cw.valid);
    const double expected = 1.0 / std::sqrt(3.0);   // cot(60°)
    EXPECT_NEAR(cw.cot1, expected, 1e-12);
    EXPECT_NEAR(cw.cot2, expected, 1e-12);
    EXPECT_NEAR(cw.cot3, expected, 1e-12);
}

// ════════════════════════════════════════════════════════════════════════════
// Cotangent weight: right-isosceles triangle (legs 1, hypotenuse √2)
//
//   v1=(0,0): right angle → cot(90°) = 0
//   v2=(1,0), v3=(0,1): 45° angles → cot(45°) = 1
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanHessian, CotWeights_RightIsoscelesTriangle)
{
    // l12=1, l23=√2, l31=1
    auto cw = euclidean_cot_weights(1.0, std::sqrt(2.0), 1.0);

    ASSERT_TRUE(cw.valid);
    EXPECT_NEAR(cw.cot1, 0.0, 1e-12);   // right angle at v1
    EXPECT_NEAR(cw.cot2, 1.0, 1e-12);   // 45° at v2
    EXPECT_NEAR(cw.cot3, 1.0, 1e-12);   // 45° at v3
}

// ════════════════════════════════════════════════════════════════════════════
// N5: cotangent weights on a SLIVER triangle match a high-precision reference.
//
// The area is now computed by Kahan's stable side-length formula instead of
// the naive 2·√(t12·t23·t31·l123).  On a thin (sliver) triangle the naive
// product of near-equal-length differences loses precision, biasing every
// cotangent weight.  Here we cross-check against an independent law-of-cosines
// reference in long double (80-bit on x86 CI — a genuine high-precision oracle;
// equal to double on ARM64, where it still serves as a regression cross-check).
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanHessian, CotWeights_SliverMatchesHighPrecisionReference)
{
    // Thin isosceles: short base l12, two near-unit legs (apex angle ≈ 0.01 rad).
    const double l12 = 0.01, l23 = 1.0, l31 = 1.0;
    auto cw = euclidean_cot_weights(l12, l23, l31);
    ASSERT_TRUE(cw.valid);

    // cot at the vertex opposite side `opp`, with adjacent sides s1, s2:
    //   cos = (s1² + s2² − opp²) / (2·s1·s2),  sin = √((1−cos)(1+cos)).
    auto cot_ref = [](long double opp, long double s1, long double s2) {
        long double cosA = (s1 * s1 + s2 * s2 - opp * opp) / (2.0L * s1 * s2);
        long double sinA = std::sqrt((1.0L - cosA) * (1.0L + cosA));
        return static_cast<double>(cosA / sinA);
    };
    const double c1 = cot_ref(l23, l12, l31);  // v1 opposite l23
    const double c2 = cot_ref(l31, l12, l23);  // v2 opposite l31
    const double c3 = cot_ref(l12, l23, l31);  // v3 opposite l12 (needle angle)

    auto rel = [](double a, double b) {
        return std::abs(a - b) / std::max(1.0, std::abs(b));
    };
    EXPECT_LT(rel(cw.cot1, c1), 1e-9);
    EXPECT_LT(rel(cw.cot2, c2), 1e-9);
    EXPECT_LT(rel(cw.cot3, c3), 1e-9);

    EXPECT_TRUE(std::isfinite(cw.cot1) &&
                std::isfinite(cw.cot2) &&
                std::isfinite(cw.cot3));
}

// ════════════════════════════════════════════════════════════════════════════
// Hessian is symmetric: H[i,j] == H[j,i]
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanHessian, HessianIsSymmetric)
{
    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.1);
    auto H = euclidean_hessian(mesh, x, maps);

    Eigen::MatrixXd Hd = Eigen::MatrixXd(H);
    EXPECT_NEAR((Hd - Hd.transpose()).norm(), 0.0, 1e-12)
        << "Hessian must be symmetric";
}

// ════════════════════════════════════════════════════════════════════════════
// Null-space property: H·1 = 0 for a closed surface (regular tetrahedron)
//
// The cotangent Laplacian on a closed mesh has the constant vector in its
// null space (each row sums to zero).
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanHessian, NullSpaceIsConstantVector_ClosedMesh)
{
    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.0);
    auto H = euclidean_hessian(mesh, x, maps);

    // 1-vector
    Eigen::VectorXd ones = Eigen::VectorXd::Ones(n);
    Eigen::VectorXd Hones = H * ones;

    EXPECT_NEAR(Hones.norm(), 0.0, 1e-10)
        << "H·1 must be zero on a closed mesh (cotangent Laplacian null-space)";
}

// ════════════════════════════════════════════════════════════════════════════
// Hessian is positive semi-definite: all eigenvalues ≥ 0
//
// Checked on a small mesh (regular tetrahedron, 4 vertices) using dense
// self-adjoint eigenvalue decomposition (only feasible for small n).
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanHessian, HessianIsPositiveSemiDefinite)
{
    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.0);
    auto H = euclidean_hessian(mesh, x, maps);

    Eigen::MatrixXd Hd = Eigen::MatrixXd(H);
    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(Hd);
    double min_ev = es.eigenvalues().minCoeff();

    EXPECT_GE(min_ev, -1e-10)
        << "All eigenvalues of the cotangent Laplacian must be ≥ 0; "
           "smallest = " << min_ev;
}

// ════════════════════════════════════════════════════════════════════════════
// Finite-difference Hessian check: single right-isosceles triangle
//
// H[i,j] ≈ (G_i(x+ε·eⱼ) − G_i(x−ε·eⱼ)) / (2ε)
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanHessian, FDCheck_Triangle)
{
    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);

    EXPECT_TRUE(hessian_check_euclidean(mesh, x, maps))
        << "FD Hessian check failed on right-isosceles triangle";
}

// ════════════════════════════════════════════════════════════════════════════
// Finite-difference Hessian check: quad strip (2 triangles, 1 interior edge)
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanHessian, FDCheck_QuadStrip)
{
    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.1);

    EXPECT_TRUE(hessian_check_euclidean(mesh, x, maps))
        << "FD Hessian check failed on quad strip";
}

// ════════════════════════════════════════════════════════════════════════════
// Finite-difference Hessian check: regular tetrahedron (closed, 4 faces)
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanHessian, FDCheck_Tetrahedron)
{
    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(hessian_check_euclidean(mesh, x, maps))
        << "FD Hessian check failed on regular tetrahedron";
}

// ════════════════════════════════════════════════════════════════════════════
// Finite-difference Hessian check: with mixed pinned/variable vertices
//
// One vertex pinned: the corresponding row/column must be absent from H
// while the diagonal of neighbouring variable vertices still gets the full
// cotangent contribution.
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanHessian, FDCheck_MixedPinnedVertices)
{
    auto mesh = make_quad_strip();
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);

    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.2, -0.15};

    EXPECT_TRUE(hessian_check_euclidean(mesh, x, maps))
        << "FD Hessian check failed for mixed pinned/variable vertices";
}

// ════════════════════════════════════════════════════════════════════════════
// Edge-DOF guard (Finding-G, java-port-audit item 2)
//
// euclidean_hessian() (vertex-only cotangent Laplacian) must throw
// std::logic_error when any edge DOF is active.  Without this guard the
// function would silently return a Hessian with zero rows/cols for the
// edge DOFs, causing SimplicialLDLT to fail in a hard-to-diagnose way.
// ════════════════════════════════════════════════════════════════════════════

TEST(EuclideanHessian, EdgeDOFGuard_Throws)
{
    auto mesh = make_tetrahedron();
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);
    assign_euclidean_all_dof_indices(mesh, maps);   // assigns vertex + edge DOFs

    const int n = euclidean_dimension(mesh, maps);
    std::vector<double> x(static_cast<std::size_t>(n), 0.0);

    EXPECT_THROW(euclidean_hessian(mesh, x, maps), std::logic_error)
        << "euclidean_hessian must throw when edge DOFs are present";
}

TEST(EuclideanHessian, EdgeDOFGuard_VertexOnlyDoesNotThrow)
{
    // Vertex-only layout must NOT trigger the guard.
    auto mesh = make_tetrahedron();
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);
    auto gauge = *mesh.vertices().begin();
    assign_euclidean_vertex_dof_indices(mesh, maps, gauge);

    const int n = euclidean_dimension(mesh, maps);
    std::vector<double> x(static_cast<std::size_t>(n), 0.0);

    EXPECT_NO_THROW(euclidean_hessian(mesh, x, maps))
        << "euclidean_hessian must not throw for vertex-only DOF layout";
}