#pragma once
// spherical_hessian.hpp
//
// Analytical Hessian of the spherical discrete conformal energy —
// the spherical cotangent-Laplace operator.
//
// Ported from de.varylab.discreteconformal.functional.SphericalFunctional
// (the hessian() method, vertex DOFs).
//
// ┌──────────────────────────────────────────────────────────────────────────┐
// │  Hessian formula (vertex DOFs only)                                      │
// │                                                                          │
// │  For a spherical face (v1, v2, v3) with vertex angles α1, α2, α3:       │
// │                                                                          │
// │    Spherical cotangent weight for edge (vi, vj) with opposite vk:        │
// │      β_k = (π − αi − αj + αk) / 2                                       │
// │      w_k  = cot(β_k)  =  1/tan(β_k)                                     │
// │                                                                          │
// │    Euclidean limit: α1+α2+α3 → π, β_k → αk, w_k → cot(αk). ✓          │
// │                                                                          │
// │  Hessian contributions per face:                                         │
// │    H[vi, vi] += w_ij + w_ik      (diagonal: weights of incident edges)  │
// │    H[vi, vj] -= w_ij             (off-diagonal: weight of edge ij)      │
// │                                                                          │
// │  where w_ij is the weight of the edge between vi and vj (opposite vk):  │
// │    w_ij = cot(β_k)  with  β_k = (π − αi − αj + αk) / 2.               │
// └──────────────────────────────────────────────────────────────────────────┘
//
// Requires Eigen (header-only).  Returns Eigen::SparseMatrix<double>.

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

namespace conformallab {

// ── Spherical cotangent weight helper ────────────────────────────────────────
//
// Given the three face angles α1, α2, α3 of a spherical triangle, return the
// three edge cotangent weights w1 (edge opp v1), w2 (edge opp v2), w3 (edge opp v3).
//
// w_k = cot(β_k)  where β_k = (π − α_adj1 − α_adj2 + α_opp) / 2
//                             = (π − αi − αj + αk) / 2  for edge (vi,vj), opposite vk
//
// Mapping in our CGAL halfedge convention:
//   h0 = halfedge(f): edge v1-v2 → opposite v3 → w = cot(β3), β3=(π-α1-α2+α3)/2
//   h1: edge v2-v3 → opposite v1 → w = cot(β1), β1=(π-α2-α3+α1)/2
//   h2: edge v3-v1 → opposite v2 → w = cot(β2), β2=(π-α3-α1+α2)/2
//
// Returns valid=false if any β_k is out of range (degenerate face).
struct SpherCotWeights { double w12, w23, w31; bool valid; };

inline SpherCotWeights spherical_cot_weights(double alpha1, double alpha2, double alpha3)
{
    // β for each edge:
    //   β3 = (π - α1 - α2 + α3)/2  — weight for edge v1-v2 (opposite v3)
    //   β1 = (π - α2 - α3 + α1)/2  — weight for edge v2-v3 (opposite v1)
    //   β2 = (π - α3 - α1 + α2)/2  — weight for edge v3-v1 (opposite v2)
    const double beta3 = (PI - alpha1 - alpha2 + alpha3) * 0.5;
    const double beta1 = (PI - alpha2 - alpha3 + alpha1) * 0.5;
    const double beta2 = (PI - alpha3 - alpha1 + alpha2) * 0.5;

    // Each β_k must be in (0, π/2] for the weight to be positive and well-defined.
    // For degenerate or very flat triangles some β may be ≤ 0 or ≥ π/2.
    if (beta1 <= 0.0 || beta2 <= 0.0 || beta3 <= 0.0) return {0.0, 0.0, 0.0, false};

    const double tb1 = std::tan(beta1);
    const double tb2 = std::tan(beta2);
    const double tb3 = std::tan(beta3);

    if (std::abs(tb1) < 1e-15 || std::abs(tb2) < 1e-15 || std::abs(tb3) < 1e-15)
        return {0.0, 0.0, 0.0, false};

    // w_ij = cot(β_k) where β_k is for the edge opposite vk.
    // w12 is for edge v1-v2 (opposite v3): cot(β3)
    // w23 is for edge v2-v3 (opposite v1): cot(β1)
    // w31 is for edge v3-v1 (opposite v2): cot(β2)
    return {1.0 / tb3, 1.0 / tb1, 1.0 / tb2, true};
}

// ── Analytical Hessian ────────────────────────────────────────────────────────
//
// Returns the n×n sparse Hessian matrix H where n = spherical_dimension(mesh, m).
// x – current DOF vector.
//
// Derivation: G_v = θ_v − Σ_f α_v^f  →  H[i,j] = −Σ_f ∂α_i^f/∂u_j
//
// For a face (v1,v2,v3) with arc-lengths l12,l23,l31 and angles α1,α2,α3,
// differentiating the spherical law of cosines
//   cos(l_opp) = cos(l_a)cos(l_b) + sin(l_a)sin(l_b)cos(α)
// gives:
//   ∂α1/∂l12 = [cot(l12)cos(α1) − cot(l31)] / sin(α1)     (adjacent side)
//   ∂α1/∂l31 = [cot(l31)cos(α1) − cot(l12)] / sin(α1)     (adjacent side)
//   ∂α1/∂l23 = sin(l23) / [sin(l12)sin(l31)sin(α1)]        (opposite side)
//
// Chain rule with ∂l_ij/∂u_k = tan(l_ij/2)  (from l = 2·asin(exp(λ/2))):
//   ∂α1/∂u1 = ∂α1/∂l12·t12 + ∂α1/∂l31·t31
//   ∂α1/∂u2 = ∂α1/∂l12·t12 + ∂α1/∂l23·t23
//   ∂α1/∂u3 = ∂α1/∂l23·t23 + ∂α1/∂l31·t31
inline Eigen::SparseMatrix<double> spherical_hessian(
    ConformalMesh&             mesh,
    const std::vector<double>& x,
    const SphericalMaps&       m)
{
    const int n = spherical_dimension(mesh, m);

    std::vector<Eigen::Triplet<double>> trips;
    trips.reserve(static_cast<std::size_t>(n) * 9);

    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 = spher_dof_val(m.v_idx[v1], x);
        double u2 = spher_dof_val(m.v_idx[v2], x);
        double u3 = spher_dof_val(m.v_idx[v3], x);

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

        const double l12 = spherical_l(lam12);
        const double l23 = spherical_l(lam23);
        const double l31 = spherical_l(lam31);

        SphericalFaceAngles fa = spherical_angles(l12, l23, l31);
        if (!fa.valid) continue;

        const double sinl12 = std::sin(l12), cosl12 = std::cos(l12);
        const double sinl23 = std::sin(l23), cosl23 = std::cos(l23);
        const double sinl31 = std::sin(l31), cosl31 = std::cos(l31);

        if (sinl12 < 1e-15 || sinl23 < 1e-15 || sinl31 < 1e-15) continue;

        const double cot12 = cosl12 / sinl12;
        const double cot23 = cosl23 / sinl23;
        const double cot31 = cosl31 / sinl31;

        // ∂l_ij/∂λ_ij = tan(l_ij/2)
        const double t12 = std::tan(l12 * 0.5);
        const double t23 = std::tan(l23 * 0.5);
        const double t31 = std::tan(l31 * 0.5);

        const double sinA1 = std::sin(fa.alpha1), cosA1 = std::cos(fa.alpha1);
        const double sinA2 = std::sin(fa.alpha2), cosA2 = std::cos(fa.alpha2);
        const double sinA3 = std::sin(fa.alpha3), cosA3 = std::cos(fa.alpha3);

        if (sinA1 < 1e-15 || sinA2 < 1e-15 || sinA3 < 1e-15) continue;

        // ∂α1/∂l_jk  (α1 at v1; opposite l23, adjacent l12,l31)
        const double dA1_dl12 = (cot12 * cosA1 - cot31) / sinA1;
        const double dA1_dl31 = (cot31 * cosA1 - cot12) / sinA1;
        const double dA1_dl23 = sinl23 / (sinl12 * sinl31 * sinA1);

        // ∂α2/∂l_jk  (α2 at v2; opposite l31, adjacent l12,l23)
        const double dA2_dl12 = (cot12 * cosA2 - cot23) / sinA2;
        const double dA2_dl23 = (cot23 * cosA2 - cot12) / sinA2;
        const double dA2_dl31 = sinl31 / (sinl12 * sinl23 * sinA2);

        // ∂α3/∂l_jk  (α3 at v3; opposite l12, adjacent l23,l31)
        const double dA3_dl23 = (cot23 * cosA3 - cot31) / sinA3;
        const double dA3_dl31 = (cot31 * cosA3 - cot23) / sinA3;
        const double dA3_dl12 = sinl12 / (sinl23 * sinl31 * sinA3);

        // Chain rule: ∂α_i/∂u_j   (u1 affects l12,l31; u2 affects l12,l23; u3 affects l23,l31)
        const double dA1_du1 = dA1_dl12 * t12 + dA1_dl31 * t31;
        const double dA1_du2 = dA1_dl12 * t12 + dA1_dl23 * t23;
        const double dA1_du3 = dA1_dl23 * t23 + dA1_dl31 * t31;

        const double dA2_du1 = dA2_dl12 * t12 + dA2_dl31 * t31;
        const double dA2_du2 = dA2_dl12 * t12 + dA2_dl23 * t23;
        const double dA2_du3 = dA2_dl23 * t23 + dA2_dl31 * t31;

        const double dA3_du1 = dA3_dl12 * t12 + dA3_dl31 * t31;
        const double dA3_du2 = dA3_dl12 * t12 + dA3_dl23 * t23;
        const double dA3_du3 = dA3_dl23 * t23 + dA3_dl31 * t31;

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

        // H[vi, vj] -= ∂α_i/∂u_j   (G_v = θ_v − Σ α_v, so ∂G_i/∂u_j = −∂α_i/∂u_j)
        if (i1 >= 0) trips.emplace_back(i1, i1, -dA1_du1);
        if (i2 >= 0) trips.emplace_back(i2, i2, -dA2_du2);
        if (i3 >= 0) trips.emplace_back(i3, i3, -dA3_du3);

        if (i1 >= 0 && i2 >= 0) {
            trips.emplace_back(i1, i2, -dA1_du2);
            trips.emplace_back(i2, i1, -dA2_du1);
        }
        if (i2 >= 0 && i3 >= 0) {
            trips.emplace_back(i2, i3, -dA2_du3);
            trips.emplace_back(i3, i2, -dA3_du2);
        }
        if (i3 >= 0 && i1 >= 0) {
            trips.emplace_back(i3, i1, -dA3_du1);
            trips.emplace_back(i1, i3, -dA1_du3);
        }
    }

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

// ── Finite-difference Hessian check ──────────────────────────────────────────
//
// Compares the analytical Hessian column-by-column against
//   H_fd[:, j] = (G(x + ε·eⱼ) − G(x − ε·eⱼ)) / (2ε).
inline bool hessian_check_spherical(
    ConformalMesh&             mesh,
    const std::vector<double>& x0,
    const SphericalMaps&       m,
    double                     eps = 1e-5,
    double                     tol = 1e-4)
{
    const int n = static_cast<int>(x0.size());
    auto H = spherical_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 = spherical_gradient(mesh, xp, m);
        auto Gm = spherical_gradient(mesh, xm, m);

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

        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