From 516ac89bd87b066ba059f2de8845bcf3467f39e0 Mon Sep 17 00:00:00 2001 From: Tarik Moussa Date: Mon, 11 May 2026 23:10:23 +0200 Subject: [PATCH] feat(phase3b): port HyperIdealFunctional energy + gradient onto ConformalMesh MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Implements the hyper-ideal discrete conformal map functional on CGAL::Surface_mesh. The energy and analytic gradient are ported directly from HyperIdealFunctional.java; correctness is verified via a finite-difference gradient check (same eps=1E-5 / tol=1E-4 as Java). New files: include/hyper_ideal_geometry.hpp — ζ, ζ₁₃, ζ₁₄, ζ₁₅, lij, αij, σi, σij include/hyper_ideal_functional.hpp — HyperIdealMaps, evaluate_hyper_ideal, gradient_check tests/cgal/test_hyper_ideal_functional.cpp — 6 tests (1 skipped @Ignore) Test results (local, -DWITH_CGAL=ON): conformallab_cgal_tests: 21 registered | 20 passed | 1 skipped | 0 failed - GradientCheck_AllHyperIdealTriangle ✓ - GradientCheck_ExtendedDomain ✓ - GradientCheck_TetrahedronAllVariable ✓ - EnergyFiniteAtTestPoint ✓ - GradientCheck_MixedIdealHyperIdeal ✓ - GradientCheck_Fan6AllVariable ✓ Co-Authored-By: Claude Sonnet 4.6 --- code/include/hyper_ideal_functional.hpp | 344 ++++++++++++++++++ code/include/hyper_ideal_geometry.hpp | 138 +++++++ code/tests/cgal/CMakeLists.txt | 4 +- .../cgal/test_hyper_ideal_functional.cpp | 184 ++++++++++ 4 files changed, 668 insertions(+), 2 deletions(-) create mode 100644 code/include/hyper_ideal_functional.hpp create mode 100644 code/include/hyper_ideal_geometry.hpp create mode 100644 code/tests/cgal/test_hyper_ideal_functional.cpp diff --git a/code/include/hyper_ideal_functional.hpp b/code/include/hyper_ideal_functional.hpp new file mode 100644 index 0000000..e63a703 --- /dev/null +++ b/code/include/hyper_ideal_functional.hpp @@ -0,0 +1,344 @@ +#pragma once +// hyper_ideal_functional.hpp +// +// Energy and gradient of the hyper-ideal discrete conformal map functional +// evaluated on a ConformalMesh (CGAL::Surface_mesh). +// +// Ported from de.varylab.discreteconformal.functional.HyperIdealFunctional. +// +// ┌─────────────────────────────────────────────────────────────────────────┐ +// │ E(x) = Σ_faces U(f) − Σ_edges θ_e · a_e − Σ_vertices Θ_v · b_v │ +// │ │ +// │ ∂E/∂b_v = Σ_{faces adj. v} β_v(face) − Θ_v │ +// │ ∂E/∂a_e = α_e(face⁺) + α_e(face⁻) − θ_e │ +// └─────────────────────────────────────────────────────────────────────────┘ +// +// DOF vector layout (matches getDimension() ordering): +// x[v_idx[v]] = b_v for each variable vertex (log scale factor) +// x[e_idx[e]] = a_e for each variable edge (intersection angle) +// -1 in v_idx / e_idx means "pinned" (ideal / fixed at 0). +// +// Usage +// ───── +// auto mesh = make_tetrahedron(); +// auto maps = setup_hyper_ideal_maps(mesh); +// int n = assign_all_dof_indices(mesh, maps); // all variable +// std::vector x(n, 1.0); +// auto res = evaluate_hyper_ideal(mesh, x, maps); +// // res.energy, res.gradient + +#include "conformal_mesh.hpp" +#include "hyper_ideal_geometry.hpp" +#include "hyper_ideal_utility.hpp" +#include +#include +#include +#include + +namespace conformallab { + +// ── Property-map type aliases ───────────────────────────────────────────────── + +using VMapD = ConformalMesh::Property_map; +using VMapI = ConformalMesh::Property_map; +using EMapD = ConformalMesh::Property_map; +using EMapI = ConformalMesh::Property_map; + +// ── Persistent map bundle ───────────────────────────────────────────────────── + +struct HyperIdealMaps { + VMapI v_idx; // DOF index per vertex (-1 = pinned / ideal point) + EMapI e_idx; // DOF index per edge (-1 = fixed) + VMapD theta_v; // target cone angle Θ_v (parameter, not variable) + EMapD theta_e; // target intersection angle θ_e +}; + +// Add all needed persistent property maps and return handles. +// Defaults: theta_v = 2π (regular cone), theta_e = π (orthogonal circles). +inline HyperIdealMaps setup_hyper_ideal_maps(ConformalMesh& mesh) +{ + HyperIdealMaps m; + m.v_idx = mesh.add_property_map ("v:idx", -1 ).first; + m.e_idx = mesh.add_property_map ("e:idx", -1 ).first; + m.theta_v = mesh.add_property_map("v:theta", 2.0*PI ).first; + m.theta_e = mesh.add_property_map("e:theta", PI ).first; + return m; +} + +// Count variable DOFs: #variable_vertices + #variable_edges. +inline int hyper_ideal_dimension(const ConformalMesh& mesh, const HyperIdealMaps& m) +{ + int dim = 0; + for (auto v : mesh.vertices()) if (m.v_idx[v] >= 0) ++dim; + for (auto e : mesh.edges()) if (m.e_idx[e] >= 0) ++dim; + return dim; +} + +// Assign DOF indices 0..n-1: vertices first, then edges. +// All vertices and edges become variable. Returns total DOF count. +inline int assign_all_dof_indices(ConformalMesh& mesh, HyperIdealMaps& m) +{ + int idx = 0; + for (auto v : mesh.vertices()) m.v_idx[v] = idx++; + for (auto e : mesh.edges()) m.e_idx[e] = idx++; + return idx; +} + +// ── Evaluation result ───────────────────────────────────────────────────────── + +struct HyperIdealResult { + double energy = 0.0; + std::vector gradient; // empty when gradient was not requested +}; + +// ── Internal helpers ────────────────────────────────────────────────────────── + +// Get the DOF value from x, or 0.0 if pinned. +static inline double dof_val(int idx, const std::vector& x) +{ + return idx >= 0 ? x[static_cast(idx)] : 0.0; +} + +// Convert a CGAL halfedge index to a plain std::size_t (for vector indexing). +static inline std::size_t hidx(Halfedge_index h) +{ + return static_cast(static_cast(h)); +} + +// ── Per-face angle kernel ───────────────────────────────────────────────────── + +struct FaceAngles { + double alpha12, alpha23, alpha31; // dihedral angles at each edge + double beta1, beta2, beta3; // interior angles at each vertex + double a12, a23, a31; // edge DOF values (used in energy) + double b1, b2, b3; // vertex DOF values + bool v1b, v2b, v3b; // whether each vertex is variable +}; + +static FaceAngles compute_face_angles( + const ConformalMesh& mesh, + Face_index f, + const std::vector& x, + const HyperIdealMaps& m) +{ + 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); + + FaceAngles fa; + fa.v1b = m.v_idx[v1] >= 0; + fa.v2b = m.v_idx[v2] >= 0; + fa.v3b = m.v_idx[v3] >= 0; + + fa.a12 = dof_val(m.e_idx[e12], x); + fa.a23 = dof_val(m.e_idx[e23], x); + fa.a31 = dof_val(m.e_idx[e31], x); + fa.b1 = dof_val(m.v_idx[v1], x); + fa.b2 = dof_val(m.v_idx[v2], x); + fa.b3 = dof_val(m.v_idx[v3], x); + + // Clamp invalid inputs (mirrors Java log.warning + clamp) + if (fa.v1b && fa.v2b && fa.a12 < 0.0) fa.a12 = 0.0; + if (fa.v2b && fa.v3b && fa.a23 < 0.0) fa.a23 = 0.0; + if (fa.v3b && fa.v1b && fa.a31 < 0.0) fa.a31 = 0.0; + if (fa.v1b && fa.b1 < 0.0) fa.b1 = 0.01; + if (fa.v2b && fa.b2 < 0.0) fa.b2 = 0.01; + if (fa.v3b && fa.b3 < 0.0) fa.b3 = 0.01; + + double l12 = lij(fa.b1, fa.b2, fa.a12, fa.v1b, fa.v2b); + double l23 = lij(fa.b2, fa.b3, fa.a23, fa.v2b, fa.v3b); + double l31 = lij(fa.b3, fa.b1, fa.a31, fa.v3b, fa.v1b); + + // Guard degenerate lengths + if (l12 < 1E-12 && l23 < 1E-12 && l31 < 1E-12) + l12 = l23 = l31 = 1E-12; + + // Check triangle inequalities; degenerate cases get extreme angles + if (l12 > l23 + l31) { + fa.beta1 = 0.0; fa.beta2 = 0.0; fa.beta3 = PI; + fa.alpha12 = PI; fa.alpha23 = 0.0; fa.alpha31 = 0.0; + } else if (l23 > l12 + l31) { + fa.beta1 = PI; fa.beta2 = 0.0; fa.beta3 = 0.0; + fa.alpha12 = 0.0; fa.alpha23 = PI; fa.alpha31 = 0.0; + } else if (l31 > l12 + l23) { + fa.beta1 = 0.0; fa.beta2 = PI; fa.beta3 = 0.0; + fa.alpha12 = 0.0; fa.alpha23 = 0.0; fa.alpha31 = PI; + } else { + fa.beta1 = zeta(l12, l31, l23); + fa.beta2 = zeta(l23, l12, l31); + fa.beta3 = zeta(l31, l23, l12); + + fa.alpha12 = alpha_ij(fa.a12, fa.a23, fa.a31, + fa.b1, fa.b2, fa.b3, + fa.beta1, fa.beta2, fa.beta3, + fa.v1b, fa.v2b, fa.v3b); + fa.alpha23 = alpha_ij(fa.a23, fa.a31, fa.a12, + fa.b2, fa.b3, fa.b1, + fa.beta2, fa.beta3, fa.beta1, + fa.v2b, fa.v3b, fa.v1b); + fa.alpha31 = alpha_ij(fa.a31, fa.a12, fa.a23, + fa.b3, fa.b1, fa.b2, + fa.beta3, fa.beta1, fa.beta2, + fa.v3b, fa.v1b, fa.v2b); + } + return fa; +} + +// Per-face energy contribution U(f) (before subtracting θ·a and Θ·b terms). +static double face_energy(const FaceAngles& fa) +{ + double aa = fa.a12*fa.alpha12 + fa.a23*fa.alpha23 + fa.a31*fa.alpha31; + double bb = fa.b1 *fa.beta1 + fa.b2 *fa.beta2 + fa.b3 *fa.beta3; + + double V = 0.0; + if (fa.v1b && fa.v2b && fa.v3b) { + V = calculateTetrahedronVolume( + fa.beta1, fa.beta2, fa.beta3, + fa.alpha23, fa.alpha31, fa.alpha12); + } else if (!fa.v1b) { + V = calculateTetrahedronVolumeWithIdealVertexAtGamma( + fa.beta1, fa.alpha31, fa.alpha12, + fa.alpha23, fa.beta2, fa.beta3); + } else if (!fa.v2b) { + V = calculateTetrahedronVolumeWithIdealVertexAtGamma( + fa.beta2, fa.alpha12, fa.alpha23, + fa.alpha31, fa.beta3, fa.beta1); + } else { // !v3b + V = calculateTetrahedronVolumeWithIdealVertexAtGamma( + fa.beta3, fa.alpha23, fa.alpha31, + fa.alpha12, fa.beta1, fa.beta2); + } + return aa + bb + 2.0 * V; +} + +// ── Full evaluation ─────────────────────────────────────────────────────────── + +inline HyperIdealResult evaluate_hyper_ideal( + ConformalMesh& mesh, + const std::vector& x, + const HyperIdealMaps& m, + bool need_energy = true, + bool need_gradient = true) +{ + HyperIdealResult res; + + // Temporary per-halfedge storage for computed angles. + // Indexed by the integer value of Halfedge_index. + const std::size_t nh = mesh.number_of_halfedges(); + std::vector h_alpha(nh, 0.0); // α_ij stored on halfedge + std::vector h_beta (nh, 0.0); // β_i stored on opposite halfedge + + // ── Pass 1: angles + energy per face ───────────────────────────────────── + for (auto f : mesh.faces()) { + Halfedge_index h0 = mesh.halfedge(f); + Halfedge_index h1 = mesh.next(h0); + Halfedge_index h2 = mesh.next(h1); + + FaceAngles fa = compute_face_angles(mesh, f, x, m); + + // Store computed angles into temporary arrays. + // h_alpha[h] = α for the edge of h in this face. + h_alpha[hidx(h0)] = fa.alpha12; + h_alpha[hidx(h1)] = fa.alpha23; + h_alpha[hidx(h2)] = fa.alpha31; + + // h_beta[h] = β at the vertex OPPOSITE to h. + // β1 (at v1 = source(h0)) is opposite to h1 = e23. + h_beta[hidx(h1)] = fa.beta1; // h1 is across from v1 + h_beta[hidx(h2)] = fa.beta2; // h2 is across from v2 + h_beta[hidx(h0)] = fa.beta3; // h0 is across from v3 + + if (need_energy) + res.energy += face_energy(fa); + } + + // ── Pass 2: linear energy terms ────────────────────────────────────────── + if (need_energy) { + for (auto e : mesh.edges()) { + int ie = m.e_idx[e]; + if (ie >= 0) res.energy -= m.theta_e[e] * x[static_cast(ie)]; + } + for (auto v : mesh.vertices()) { + int iv = m.v_idx[v]; + if (iv >= 0) res.energy -= m.theta_v[v] * x[static_cast(iv)]; + } + } + + // ── Pass 3: gradient ───────────────────────────────────────────────────── + if (need_gradient) { + const int n = hyper_ideal_dimension(mesh, m); + res.gradient.assign(static_cast(n), 0.0); + + // ∂E/∂b_v = Σ_{faces adj. v} β_v(face) − Θ_v + // β_v(face) = h_beta[prev(h)] for the incoming halfedge h to v in that face. + for (auto v : mesh.vertices()) { + int iv = m.v_idx[v]; + if (iv < 0) continue; + for (auto h : CGAL::halfedges_around_target(v, mesh)) { + if (mesh.is_border(h)) continue; + res.gradient[static_cast(iv)] += h_beta[hidx(mesh.prev(h))]; + } + res.gradient[static_cast(iv)] -= m.theta_v[v]; + } + + // ∂E/∂a_e = α_e(face⁺) + α_e(face⁻) − θ_e + for (auto e : mesh.edges()) { + int ie = m.e_idx[e]; + if (ie < 0) continue; + auto h = mesh.halfedge(e); + auto ho = mesh.opposite(h); + if (!mesh.is_border(h)) res.gradient[static_cast(ie)] += h_alpha[hidx(h)]; + if (!mesh.is_border(ho)) res.gradient[static_cast(ie)] += h_alpha[hidx(ho)]; + res.gradient[static_cast(ie)] -= m.theta_e[e]; + } + } + + return res; +} + +// ── Finite-difference gradient check ───────────────────────────────────────── +// +// Returns true if |G[i] − fd[i]| / max(1, |G[i]|) < tol for all DOFs. +// eps = step size, tol = tolerance (same defaults as Java FunctionalTest). +inline bool gradient_check( + ConformalMesh& mesh, + const std::vector& x0, + const HyperIdealMaps& m, + double eps = 1E-5, + double tol = 1E-4) +{ + // Analytic gradient + auto res = evaluate_hyper_ideal(mesh, x0, m, false, true); + const auto& G = res.gradient; + const int n = static_cast(G.size()); + + std::vector xp = x0, xm = x0; + bool ok = true; + + for (int i = 0; i < n; ++i) { + std::size_t si = static_cast(i); + xp[si] = x0[si] + eps; + xm[si] = x0[si] - eps; + + double Ep = evaluate_hyper_ideal(mesh, xp, m, true, false).energy; + double Em = evaluate_hyper_ideal(mesh, xm, m, true, false).energy; + + xp[si] = xm[si] = x0[si]; + + double fd = (Ep - Em) / (2.0 * eps); + double err = std::abs(G[si] - fd); + double scale = std::max(1.0, std::abs(G[si])); + if (err / scale > tol) ok = false; + } + return ok; +} + +} // namespace conformallab diff --git a/code/include/hyper_ideal_geometry.hpp b/code/include/hyper_ideal_geometry.hpp new file mode 100644 index 0000000..72e0e6a --- /dev/null +++ b/code/include/hyper_ideal_geometry.hpp @@ -0,0 +1,138 @@ +#pragma once +// hyper_ideal_geometry.hpp +// +// Pure-math building blocks for the hyper-ideal discrete conformal map. +// Ported from de.varylab.discreteconformal.functional.HyperIdealUtility +// and the private helpers of HyperIdealFunctional (lij, αij, σi, σij). +// +// All functions are independent of the mesh type. +// +// Notation follows the original Java / paper: +// b_i, b_j – vertex variables (log scale factors, hyper-ideal vertices) +// a_ij – edge variable (intersection angle between horocycles) +// l_ij – effective hyperbolic edge length in the auxiliary triangle +// β_i – interior angle of the hyperbolic triangle at vertex i +// α_ij – dihedral angle of the tetrahedron at edge ij + +#include +#include + +namespace conformallab { + +constexpr double PI = 3.14159265358979323846264338328; + +// ── Length functions ───────────────────────────────────────────────────────── + +// ζ(x,y,z) — interior angle in a hyperbolic triangle with edge lengths +// x, y, z, opposite to the side of length z. +// Ports HyperIdealUtility.ζ(x, y, z). +inline double zeta(double x, double y, double z) +{ + double cx = std::cosh(x), cy = std::cosh(y), cz = std::cosh(z); + double sx = std::sinh(x), sy = std::sinh(y); + double nbd = (cx*cy - cz) / (sx*sy); + nbd = std::clamp(nbd, -1.0, 1.0); // guard floating-point rounding + return std::acos(nbd); +} + +// ζ₁₃(x,y,z) — third edge length in a right-angled hyperbolic hexagon. +// Ports HyperIdealUtility.ζ_13(x, y, z). +inline double zeta13(double x, double y, double z) +{ + double cx = std::cosh(x), cy = std::cosh(y), cz = std::cosh(z); + double sx = std::sinh(x), sy = std::sinh(y); + return std::acosh((cx*cy + cz) / (sx*sy)); +} + +// ζ₁₄(x,y) — edge length in a hyperbolic pentagon with one ideal vertex. +// Ports HyperIdealUtility.ζ_14(x, y). +inline double zeta14(double x, double y) +{ + double cy = std::cosh(y), sy = std::sinh(y); + return std::acosh((std::exp(x) + cy) / sy); +} + +// ζ₁₅(x) — length in a hyperbolic quadrilateral with two ideal vertices. +// Ports HyperIdealUtility.ζ_15(x). +inline double zeta15(double x) +{ + return 2.0 * std::asinh(std::exp(x / 2.0)); +} + +// ── Effective edge length ───────────────────────────────────────────────────── + +// l_ij: effective hyperbolic length of edge ij. +// b_i, b_j – vertex log scale factors (used only if vertex is hyper-ideal) +// a_ij – edge intersection-angle variable +// vi_var – true if vertex i is hyper-ideal (has a DOF b_i) +// vj_var – true if vertex j is hyper-ideal +// Ports HyperIdealFunctional.lij(). +inline double lij(double bi, double bj, double aij, bool vi_var, bool vj_var) +{ + if (vi_var && vj_var) return zeta13(bi, bj, aij); + if (vi_var) return zeta14(aij, bi); + if (vj_var) return zeta14(aij, bj); + return zeta15(aij); +} + +// ── Auxiliary angle functions ───────────────────────────────────────────────── + +// σᵢ(aᵢⱼ, aₖᵢ, aⱼₖ, vj_var, vk_var) — intermediate half-length at vertex i. +// Ports HyperIdealFunctional.σi(). +inline double sigma_i(double aij, double aki, double ajk, bool vj_var, bool vk_var) +{ + if (vj_var && vk_var) return zeta13(aij, aki, ajk); + if (vj_var) return zeta14(ajk - aki, aij); + if (vk_var) return zeta14(ajk - aij, aki); + return zeta15(ajk - aij - aki); +} + +// σᵢⱼ(aᵢⱼ, bᵢ, bⱼ, vj_var) — intermediate half-length for edge ij from vertex i. +// Ports HyperIdealFunctional.σij(). +inline double sigma_ij(double aij, double bi, double bj, bool vj_var) +{ + if (vj_var) return zeta13(aij, bi, bj); + return zeta14(-aij, bi); +} + +// α_ij: computed dihedral angle at edge ij in the face with vertices i, j, k. +// +// Arguments (cyclic role assignment): +// aij, ajk, aki – edge variables +// bi, bj, bk – vertex variables +// βi, βj, βk – interior angles of the auxiliary hyperbolic triangle +// vi_var, vj_var, vk_var – which vertices are hyper-ideal +// +// Ports HyperIdealFunctional.αij() (the private helper). +// Note: the vk_var case recurses once (never more than one level deep). +inline double alpha_ij( + double aij, double ajk, double aki, + double bi, double bj, double bk, + double beta_i, double beta_j, double beta_k, + bool vi_var, bool vj_var, bool vk_var) +{ + if (vi_var) { + double si = sigma_i (aij, aki, ajk, vj_var, vk_var); + double sij = sigma_ij(aij, bi, bj, vj_var); + double sik = sigma_ij(aki, bi, bk, vk_var); + return zeta(si, sij, sik); + } + if (vj_var) { + double sj = sigma_i (ajk, aij, aki, vk_var, vi_var); + double sjk = sigma_ij(ajk, bj, bk, vk_var); + double sji = sigma_ij(aij, bj, bi, vi_var); + return zeta(sj, sji, sjk); + } + if (vk_var) { + // Derive α_ij from α_jk (one level of recursion). + double a_jk = alpha_ij(ajk, aki, aij, + bj, bk, bi, + beta_j, beta_k, beta_i, + vj_var, vk_var, vi_var); + return PI - a_jk - beta_j; + } + // All ideal: closed-form formula. + return 0.5 * (PI + beta_k - beta_i - beta_j); +} + +} // namespace conformallab diff --git a/code/tests/cgal/CMakeLists.txt b/code/tests/cgal/CMakeLists.txt index ad1e0b3..229cbb8 100644 --- a/code/tests/cgal/CMakeLists.txt +++ b/code/tests/cgal/CMakeLists.txt @@ -12,8 +12,8 @@ add_executable(conformallab_cgal_tests # ── Phase 3a: mesh infrastructure ────────────────────────────────────── test_conformal_mesh.cpp - # ── Phase 3b: HyperIdealFunctional (to be added) ────────────────────── - # test_hyper_ideal_functional.cpp + # ── Phase 3b: HyperIdealFunctional ───────────────────────────────────── + test_hyper_ideal_functional.cpp # ── Phase 3c: SphericalFunctional (to be added) ────────────────────── # test_spherical_functional.cpp diff --git a/code/tests/cgal/test_hyper_ideal_functional.cpp b/code/tests/cgal/test_hyper_ideal_functional.cpp new file mode 100644 index 0000000..22da339 --- /dev/null +++ b/code/tests/cgal/test_hyper_ideal_functional.cpp @@ -0,0 +1,184 @@ +// test_hyper_ideal_functional.cpp +// +// Phase 3b — HyperIdealFunctional ported to ConformalMesh. +// +// Corresponds to de.varylab.discreteconformal.functional.HyperIdealFunctionalTest. +// +// The Java tests use CoHDS + HyperIdealGenerator to build complex meshes and +// then verify correctness via a finite-difference gradient check (FunctionalTest). +// Here we apply the same gradient-check strategy on simple hand-crafted meshes +// (triangle, tetrahedron) to validate the port independently of the generator. +// +// Test map (Java → C++) +// ────────────────────── +// testHessian (Ignored) → GradientCheck_Hessian (SKIPPED, not implemented) +// testGradientWithHyperIdeal… → GradientCheck_AllHyperIdealTriangle (ported) +// testGradientInExtendedDomain → GradientCheck_ExtendedDomain (ported) +// testGradientWithHyperelliptic → GradientCheck_TetrahedronAllVariable (ported) +// testFunctionalAtNaNValue → EnergyFiniteAtTestPoint (ported) + +#include "conformal_mesh.hpp" +#include "mesh_builder.hpp" +#include "hyper_ideal_functional.hpp" +#include +#include +#include + +using namespace conformallab; + +// ════════════════════════════════════════════════════════════════════════════ +// Helpers +// ════════════════════════════════════════════════════════════════════════════ + +// Build a mesh with all vertices and edges variable, and set reasonable +// DOF values: b_i = b_val, a_e = a_val. +static std::vector make_x_all_variable( + ConformalMesh& mesh, HyperIdealMaps& maps, + double b_val, double a_val) +{ + int n = assign_all_dof_indices(mesh, maps); + std::vector x(static_cast(n)); + + // Vertices first, then edges (matching assign_all_dof_indices order) + for (auto v : mesh.vertices()) + x[static_cast(maps.v_idx[v])] = b_val; + for (auto e : mesh.edges()) + x[static_cast(maps.e_idx[e])] = a_val; + + return x; +} + +// ════════════════════════════════════════════════════════════════════════════ +// @Ignore in Java: no Hessian implemented +// ════════════════════════════════════════════════════════════════════════════ + +TEST(HyperIdealFunctional, GradientCheck_Hessian) +{ + GTEST_SKIP() << "@Ignore in Java – Hessian not implemented in the functional"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Gradient check: single triangle, all vertices hyper-ideal +// +// Mirrors testGradientWithHyperIdealAndIdealPoints (simplified to single face). +// DOFs: 3 vertex b-values + 3 edge a-values = 6 total. +// Point: b_i = 1.0, a_e = 0.5. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(HyperIdealFunctional, GradientCheck_AllHyperIdealTriangle) +{ + auto mesh = make_triangle(); + auto maps = setup_hyper_ideal_maps(mesh); + auto x = make_x_all_variable(mesh, maps, /*b=*/1.0, /*a=*/0.5); + + EXPECT_TRUE(gradient_check(mesh, x, maps)) + << "Finite-difference gradient check failed on all-hyper-ideal triangle"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Gradient check in the "extended domain" +// +// Mirrors testGradientInTheExtendedDomain: larger DOF values +// (Java: x_i = 1.2 + |rnd|, so typically > 1.2). +// ════════════════════════════════════════════════════════════════════════════ + +TEST(HyperIdealFunctional, GradientCheck_ExtendedDomain) +{ + auto mesh = make_triangle(); + auto maps = setup_hyper_ideal_maps(mesh); + auto x = make_x_all_variable(mesh, maps, /*b=*/2.0, /*a=*/1.5); + + EXPECT_TRUE(gradient_check(mesh, x, maps)) + << "Finite-difference gradient check failed in extended domain"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Gradient check: tetrahedron (4 faces), all vertices and edges variable +// +// Mirrors testGradientWithHyperellipticCurve — larger mesh, closed surface. +// DOFs: 4 vertex b-values + 6 edge a-values = 10 total. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(HyperIdealFunctional, GradientCheck_TetrahedronAllVariable) +{ + auto mesh = make_tetrahedron(); + auto maps = setup_hyper_ideal_maps(mesh); + auto x = make_x_all_variable(mesh, maps, /*b=*/1.0, /*a=*/0.5); + + EXPECT_TRUE(gradient_check(mesh, x, maps)) + << "Finite-difference gradient check failed on all-variable tetrahedron"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Energy is finite (not NaN / Inf) +// +// Mirrors testFunctionalAtNaNValue: evaluates at a test point and checks +// the energy is a real number. Uses a quad-strip to exercise the interior +// edge path (adjacent faces sharing one non-border edge). +// ════════════════════════════════════════════════════════════════════════════ + +TEST(HyperIdealFunctional, EnergyFiniteAtTestPoint) +{ + auto mesh = make_quad_strip(); + auto maps = setup_hyper_ideal_maps(mesh); + int n = assign_all_dof_indices(mesh, maps); + + // Values taken from the Java testFunctionalAtNaNValue spirit: + // large-ish but positive DOF values that could expose degenerate paths. + std::vector x(static_cast(n), 1.5); + + auto res = evaluate_hyper_ideal(mesh, x, maps, true, false); + + EXPECT_FALSE(std::isnan(res.energy)) << "Energy must not be NaN"; + EXPECT_FALSE(std::isinf(res.energy)) << "Energy must not be Inf"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Gradient check: mixed vertices (some ideal = pinned, some hyper-ideal) +// +// One vertex is ideal (v_idx = -1, b = 0), the other two are hyper-ideal. +// This exercises the lij / αij branches for the ideal-vertex case. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(HyperIdealFunctional, GradientCheck_MixedIdealHyperIdeal) +{ + auto mesh = make_triangle(); + auto maps = setup_hyper_ideal_maps(mesh); + + // Make v0 ideal (pinned), v1 and v2 hyper-ideal. + auto vit = mesh.vertices().begin(); + Vertex_index v0 = *vit++; + Vertex_index v1 = *vit++; + Vertex_index v2 = *vit; + + maps.v_idx[v0] = -1; // ideal: b_0 = 0 (fixed) + maps.v_idx[v1] = 0; + maps.v_idx[v2] = 1; + + // All edges variable: indices 2, 3, 4 + int eidx = 2; + for (auto e : mesh.edges()) maps.e_idx[e] = eidx++; + + // DOF vector: [b1, b2, a_e0, a_e1, a_e2] + std::vector x = {1.0, 1.0, 0.5, 0.5, 0.5}; + + EXPECT_TRUE(gradient_check(mesh, x, maps)) + << "Finite-difference gradient check failed for mixed ideal / hyper-ideal"; +} + +// ════════════════════════════════════════════════════════════════════════════ +// Gradient check: fan mesh (n=6), all variable +// +// Exercises the full halfedge-around-vertex traversal for a high-valence +// interior vertex. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(HyperIdealFunctional, GradientCheck_Fan6AllVariable) +{ + auto mesh = make_fan(6); + auto maps = setup_hyper_ideal_maps(mesh); + auto x = make_x_all_variable(mesh, maps, /*b=*/1.0, /*a=*/0.5); + + EXPECT_TRUE(gradient_check(mesh, x, maps)) + << "Finite-difference gradient check failed on fan-6 mesh"; +}