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

// euclidean_hessian.hpp
//
// Analytical Hessian of the Euclidean discrete conformal energy —
// the cotangent-Laplace operator.
//
// Ported from de.varylab.discreteconformal.functional.EuclideanCyclicFunctional
// (the hessian() method).
//
// ┌──────────────────────────────────────────────────────────────────────────┐
// │  Hessian formula (vertex DOFs only)                                      │
// │                                                                          │
// │  For a face (v1, v2, v3) with effective log-lengths Λ̃ij and             │
// │  side lengths lij = exp(Λ̃ij/2):                                         │
// │                                                                          │
// │    t12 = −l12+l23+l31,  t23 = l12−l23+l31,  t31 = l12+l23−l31         │
// │    denom2 = 2·sqrt(t12·t23·t31·l123) = 8·Area                           │
// │                                                                          │
// │    cot_k = (t_adj1·l123 − t_adj2·t_opp) / denom2                        │
// │          = cotangent of the angle αk at vertex k                         │
// │                                                                          │
// │  Hessian contributions per face:                                         │
// │    H[vi, vi] += cot_vj + cot_vk         (diagonal, both non-opp angles) │
// │    H[vi, vj] -= cot_vk                  (off-diagonal, for variable vi,vj│
// │                                                                          │
// │  This is exactly the cotangent-Laplace operator from Pinkall–Polthier.   │
// │                                                                          │
// │  Pinned vertices (v_idx = −1) contribute to diagonal of neighbours but   │
// │  do not create a column/row in H themselves.                             │
// └──────────────────────────────────────────────────────────────────────────┘
//
// Requires Eigen (header-only).  The Hessian is returned as an
// Eigen::SparseMatrix<double> for direct use in the Phase-4 Newton solver
// (Eigen::SimplicialLDLT).
//
// The Hessian is symmetric positive semi-definite for any valid mesh with
// no degenerate faces.  The null space is spanned by the uniform-shift
// vector 1 on closed surfaces (Euler characteristic = 0).

#include "euclidean_functional.hpp"
#include <Eigen/Sparse>
#include <vector>
#include <cmath>
#include <stdexcept>

namespace conformallab {

// ── Cotangent weight helper ───────────────────────────────────────────────────
//
// Given three Euclidean SIDE LENGTHS l12, l23, l31 (already exp(Λ̃/2)),
// return the three cotangent weights (cot1, cot2, cot3).
//
// cot_k = (t_adj·l123 − t_opp·t_other) / (8·Area)
//
// Returns {0,0,0} for degenerate faces (triangle inequality violated or Area=0).
/// Three Euclidean cotangent weights `(cot1, cot2, cot3)` for the
/// vertices opposite to edges (l₂₃, l₃₁, l₁₂) of a triangle, plus a
/// `valid` flag that is `false` when the triangle is degenerate.
struct EuclCotWeights {
    double cot1; ///< Cotangent at vertex 1 (opposite to l₂₃).
    double cot2; ///< Cotangent at vertex 2 (opposite to l₃₁).
    double cot3; ///< Cotangent at vertex 3 (opposite to l₁₂).
    bool   valid;///< `false` when the triangle is degenerate (triangle inequality violated or area = 0).
};

/// Compute the three Euclidean cotangent weights from edge lengths.
/// Returns `{0,0,0,false}` for degenerate triangles.
inline EuclCotWeights euclidean_cot_weights(double l12, double l23, double l31)
{
    const double t12 = -l12 + l23 + l31;
    const double t23 = +l12 - l23 + l31;
    const double t31 = +l12 + l23 - l31;

    if (t12 <= 0.0 || t23 <= 0.0 || t31 <= 0.0)
        return {0.0, 0.0, 0.0, false};

    const double l123   = l12 + l23 + l31;
    const double denom2_sq = t12 * t23 * t31 * l123;
    if (denom2_sq <= 0.0) return {0.0, 0.0, 0.0, false};

    // denom2 = 2·sqrt(t12·t23·t31·l123) = 8·Area
    const double denom2 = 2.0 * std::sqrt(denom2_sq);

    // cot at v1 (opposite l23): adjacent t-values are t12 and t31.
    // cot at v2 (opposite l31): adjacent t-values are t12 and t23.
    // cot at v3 (opposite l12): adjacent t-values are t23 and t31.
    return {
        (t23 * l123 - t31 * t12) / denom2,  // cot1
        (t31 * l123 - t12 * t23) / denom2,  // cot2
        (t12 * l123 - t23 * t31) / denom2,  // cot3
        true
    };
}

/// Analytical Euclidean Hessian (cotangent Laplacian), sparse.
/// Only vertex DOFs are supported — the function asserts that no edge
/// DOF is variable.  `x` is used to compute effective log-lengths Λ̃ᵢⱼ.
inline Eigen::SparseMatrix<double> euclidean_hessian(
    ConformalMesh&             mesh,
    const std::vector<double>& x,
    const EuclideanMaps&       m)
{
    const int n = euclidean_dimension(mesh, m);

    // Only the vertex block of the cyclic Hessian is implemented here.  If any
    // edge DOF is variable (assign_euclidean_all_dof_indices), the edge-edge and
    // vertex-edge blocks present in the Java reference (conformalHessian) are
    // missing, which would leave singular zero rows/cols.  Fail loudly instead
    // of silently returning a rank-deficient matrix (matches the doc contract).
    for (auto e : mesh.edges()) {
        if (m.e_idx[e] >= 0)
            throw std::logic_error(
                "euclidean_hessian: edge DOFs are not supported "
                "(only the vertex-block cotangent Laplacian is implemented)");
    }

    // Collect triplets (row, col, value) — setFromTriplets sums duplicates.
    std::vector<Eigen::Triplet<double>> trips;
    trips.reserve(static_cast<std::size_t>(n) * 7);   // rough estimate

    for (auto f : mesh.faces()) {
        Halfedge_index h0 = mesh.halfedge(f);
        Halfedge_index h1 = mesh.next(h0);
        Halfedge_index h2 = mesh.next(h1);

        Vertex_index v1 = mesh.source(h0);
        Vertex_index v2 = mesh.source(h1);
        Vertex_index v3 = mesh.source(h2);

        Edge_index e12 = mesh.edge(h0);
        Edge_index e23 = mesh.edge(h1);
        Edge_index e31 = mesh.edge(h2);

        // Effective log-lengths.
        double u1 = eucl_dof_val(m.v_idx[v1], x);
        double u2 = eucl_dof_val(m.v_idx[v2], x);
        double u3 = eucl_dof_val(m.v_idx[v3], x);

        double lam12 = m.lambda0[e12] + u1 + u2 + eucl_dof_val(m.e_idx[e12], x);
        double lam23 = m.lambda0[e23] + u2 + u3 + eucl_dof_val(m.e_idx[e23], x);
        double lam31 = m.lambda0[e31] + u3 + u1 + eucl_dof_val(m.e_idx[e31], x);

        // Side lengths (centered to avoid overflow, same as euclidean_angles).
        const double mu  = (lam12 + lam23 + lam31) / 6.0;
        const double l12 = std::exp((lam12 - 2.0 * mu) * 0.5);
        const double l23 = std::exp((lam23 - 2.0 * mu) * 0.5);
        const double l31 = std::exp((lam31 - 2.0 * mu) * 0.5);

        auto [cot1, cot2, cot3, valid] = euclidean_cot_weights(l12, l23, l31);
        if (!valid) continue;

        const int i1 = m.v_idx[v1];
        const int i2 = m.v_idx[v2];
        const int i3 = m.v_idx[v3];

        // ── Diagonal contributions ──────────────────────────────────────────
        // H[v1,v1] += (cot2 + cot3)/2  (Pinkall–Polthier factor of 1/2)
        // H[v2,v2] += (cot3 + cot1)/2
        // H[v3,v3] += (cot1 + cot2)/2
        if (i1 >= 0) trips.emplace_back(i1, i1, (cot2 + cot3) * 0.5);
        if (i2 >= 0) trips.emplace_back(i2, i2, (cot3 + cot1) * 0.5);
        if (i3 >= 0) trips.emplace_back(i3, i3, (cot1 + cot2) * 0.5);

        // ── Off-diagonal contributions (only for variable pairs) ────────────
        // Edge v1-v2 opposite α3: H[v1,v2] -= cot3/2
        if (i1 >= 0 && i2 >= 0) {
            trips.emplace_back(i1, i2, -cot3 * 0.5);
            trips.emplace_back(i2, i1, -cot3 * 0.5);
        }
        // Edge v2-v3 opposite α1: H[v2,v3] -= cot1/2
        if (i2 >= 0 && i3 >= 0) {
            trips.emplace_back(i2, i3, -cot1 * 0.5);
            trips.emplace_back(i3, i2, -cot1 * 0.5);
        }
        // Edge v3-v1 opposite α2: H[v3,v1] -= cot2/2
        if (i3 >= 0 && i1 >= 0) {
            trips.emplace_back(i3, i1, -cot2 * 0.5);
            trips.emplace_back(i1, i3, -cot2 * 0.5);
        }
    }

    Eigen::SparseMatrix<double> H(n, n);
    H.setFromTriplets(trips.begin(), trips.end());
    return H;
}

// ── Finite-difference Hessian check ──────────────────────────────────────────
/// FD Hessian check for the Euclidean functional.  Compares analytic
/// `H` column-by-column to `(G(x+εeⱼ) − G(x−εeⱼ)) / (2ε)`; returns
/// `true` iff max relative error is below `tol`.
inline bool hessian_check_euclidean(
    ConformalMesh&             mesh,
    const std::vector<double>& x0,
    const EuclideanMaps&       m,
    double                     eps = 1e-5,
    double                     tol = 1e-4)
{
    const int n = static_cast<int>(x0.size());
    auto H = euclidean_hessian(mesh, x0, m);

    std::vector<double> xp = x0, xm = x0;
    bool ok = true;

    for (int j = 0; j < n; ++j) {
        const std::size_t sj = static_cast<std::size_t>(j);
        xp[sj] = x0[sj] + eps;
        xm[sj] = x0[sj] - eps;

        auto Gp = euclidean_gradient(mesh, xp, m);
        auto Gm = euclidean_gradient(mesh, xm, m);

        xp[sj] = xm[sj] = x0[sj];   // restore

        for (int i = 0; i < n; ++i) {
            double fd_ij  = (Gp[static_cast<std::size_t>(i)]
                           - Gm[static_cast<std::size_t>(i)]) / (2.0 * eps);
            double H_ij   = H.coeff(i, j);
            double err    = std::abs(H_ij - fd_ij);
            double scale  = std::max(1.0, std::abs(H_ij));
            if (err / scale > tol) ok = false;
        }
    }
    return ok;
}

} // namespace conformallab