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

// test_scalability_smoke.cpp
//
// Scalability smoke tests — convergence on large real-world meshes.
//
// PURPOSE
//   These tests verify that the Newton solver converges correctly on meshes
//   significantly larger than the unit tests (which use tiny synthetic meshes).
//   They do NOT assert on wall-clock time — timing is printed for information
//   only, so the tests remain stable on slow CI hardware (Raspberry Pi ARM64).
//
//   If Newton fails to converge here, it is a correctness regression, not a
//   performance regression.  See doc/math/complexity.md for timing context.
//
// MESHES USED
//   cathead.obj   V=131,  F=248,  genus=0, open  — small, sanity check
//   brezel.obj    V=6910, F=13824, genus=2, closed — large genus-2 mesh (χ=−2)
//   brezel2.obj   V=2622, F=5248,  genus=2, closed — smaller genus-2 mesh (χ=−2)
//
// NOTE: both brezel meshes are genus-2. The naming follows the Java original
// where "brezel2" is a different triangulation, not a different genus.
//
// EXPECTED RESULTS
//   Newton converges in < 30 iterations for all Euclidean meshes (strictly
//   convex energy, quadratic convergence from u=0).
//   Cut graph produces 2g seam edges: 2 for brezel, 4 for brezel2.
//
// Tests:
//   1. SmokeEuclidean.CatHead_SmallOpen
//   2. SmokeEuclidean.Brezel_LargeGenus2
//   3. SmokeEuclidean.Brezel2_Genus2_CutGraph

#include "conformal_mesh.hpp"
#include "mesh_io.hpp"
#include "euclidean_functional.hpp"
#include "gauss_bonnet.hpp"
#include "newton_solver.hpp"
#include "cut_graph.hpp"
#include <gtest/gtest.h>
#include <chrono>
#include <iostream>
#include <vector>
#include <cmath>
#include <string>

using namespace conformallab;
using Clock = std::chrono::steady_clock;
using Ms    = std::chrono::milliseconds;

// ── Helpers ──────────────────────────────────────────────────────────────────

static int setup_open_mesh_dofs(ConformalMesh& mesh, EuclideanMaps& maps)
{
    int idx = 0;
    for (auto v : mesh.vertices())
        maps.v_idx[v] = mesh.is_border(v) ? -1 : idx++;
    return idx;
}

static int setup_closed_mesh_dofs(ConformalMesh& mesh, EuclideanMaps& maps)
{
    auto vit = mesh.vertices().begin();
    maps.v_idx[*vit++] = -1;
    int idx = 0;
    for (; vit != mesh.vertices().end(); ++vit)
        maps.v_idx[*vit] = idx++;
    return idx;
}

static void apply_natural_theta(ConformalMesh& mesh, EuclideanMaps& maps, int n)
{
    std::vector<double> x0(static_cast<std::size_t>(n), 0.0);
    auto G0 = euclidean_gradient(mesh, x0, maps);
    for (auto v : mesh.vertices()) {
        int iv = maps.v_idx[v];
        if (iv >= 0) maps.theta_v[v] -= G0[static_cast<std::size_t>(iv)];
    }
}

// ── Test 1 — cathead.obj (V=131, F=248, open) ────────────────────────────────

TEST(SmokeEuclidean, CatHead_SmallOpen)
{
    const std::string path = std::string(CONFORMALLAB_DATA_DIR) + "/obj/cathead.obj";
    ConformalMesh mesh;
    ASSERT_NO_THROW(mesh = load_mesh(path)) << "cathead.obj not found: " << path;

    EXPECT_EQ(131u, mesh.number_of_vertices());
    EXPECT_EQ(248u, mesh.number_of_faces());

    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);
    const int n = setup_open_mesh_dofs(mesh, maps);
    apply_natural_theta(mesh, maps, n);

    // Start from a small perturbation so Newton actually iterates.
    std::vector<double> x0(static_cast<std::size_t>(n), -0.05);

    auto t0 = Clock::now();
    auto res = newton_euclidean(mesh, x0, maps, 1e-9, 200);
    auto dt  = std::chrono::duration_cast<Ms>(Clock::now() - t0).count();

    std::cout << "[SmokeEuclidean.CatHead] V=" << mesh.number_of_vertices()
              << " F=" << mesh.number_of_faces()
              << " iter=" << res.iterations
              << " ||G||=" << res.grad_inf_norm
              << " time=" << dt << "ms\n";

    EXPECT_TRUE(res.converged)   << "Newton did not converge on cathead.obj";
    EXPECT_LT(res.iterations, 30) << "Newton took ≥ 30 iterations — unexpected";
    EXPECT_LT(res.grad_inf_norm, 1e-8);
}

// ── Test 2 — brezel.obj (V=6910, F=13824, genus=2) ───────────────────────────
// Primary scalability target: largest mesh in the test suite.
// Newton is started from a small perturbation (x0 = −0.05) so it must
// actually iterate rather than exit immediately from the trivial equilibrium.

TEST(SmokeEuclidean, Brezel_LargeGenus2)
{
    const std::string path = std::string(CONFORMALLAB_DATA_DIR) + "/obj/brezel.obj";
    ConformalMesh mesh;
    ASSERT_NO_THROW(mesh = load_mesh(path)) << "brezel.obj not found: " << path;

    EXPECT_EQ(6910u,  mesh.number_of_vertices());
    EXPECT_EQ(13824u, mesh.number_of_faces());

    // Euler characteristic: V - E + F = −2 for genus-2 closed surface
    const int chi = static_cast<int>(mesh.number_of_vertices())
                  - static_cast<int>(mesh.number_of_edges())
                  + static_cast<int>(mesh.number_of_faces());
    EXPECT_EQ(-2, chi) << "brezel.obj must be genus-2 (χ=−2)";

    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);

    const int n = setup_closed_mesh_dofs(mesh, maps);
    enforce_gauss_bonnet(mesh, maps);
    apply_natural_theta(mesh, maps, n);

    // Start from a small perturbation so Newton actually iterates.
    std::vector<double> x0(static_cast<std::size_t>(n), -0.05);

    // Newton solve
    auto t0  = Clock::now();
    auto res = newton_euclidean(mesh, x0, maps, 1e-9, 200);
    auto dt_newton = std::chrono::duration_cast<Ms>(Clock::now() - t0).count();

    // Cut graph
    auto t1  = Clock::now();
    CutGraph cg = compute_cut_graph(mesh);
    auto dt_cut = std::chrono::duration_cast<Ms>(Clock::now() - t1).count();

    std::cout << "[SmokeEuclidean.Brezel] V=" << mesh.number_of_vertices()
              << " F=" << mesh.number_of_faces()
              << " iter=" << res.iterations
              << " ||G||=" << res.grad_inf_norm
              << " newton=" << dt_newton << "ms"
              << " cut=" << dt_cut << "ms\n";

    EXPECT_TRUE(res.converged)    << "Newton did not converge on brezel.obj";
    EXPECT_LT(res.iterations, 30) << "Newton took ≥ 30 iterations";
    EXPECT_LT(res.grad_inf_norm, 1e-8);

    // Genus-2: 2g = 4 seam edges
    EXPECT_EQ(4u, cg.cut_edge_indices.size())
        << "brezel.obj (genus 2) must yield 2g=4 cut edges";
    EXPECT_EQ(2, cg.genus);
}

// ── Test 3 — brezel2.obj (V=2622, F=5248, genus=2) ───────────────────────────

TEST(SmokeEuclidean, Brezel2_Genus2_CutGraph)
{
    const std::string path = std::string(CONFORMALLAB_DATA_DIR) + "/obj/brezel2.obj";
    ConformalMesh mesh;
    ASSERT_NO_THROW(mesh = load_mesh(path)) << "brezel2.obj not found: " << path;

    EXPECT_EQ(2622u, mesh.number_of_vertices());
    EXPECT_EQ(5248u, mesh.number_of_faces());

    const int chi = static_cast<int>(mesh.number_of_vertices())
                  - static_cast<int>(mesh.number_of_edges())
                  + static_cast<int>(mesh.number_of_faces());
    EXPECT_EQ(-2, chi) << "brezel2.obj must be genus-2 (χ=−2)";

    // Cut graph only — Newton on genus-2 requires full DOF setup
    // (tested separately in test_geometry_utils.cpp HomologyGenerators suite)
    auto t0  = Clock::now();
    CutGraph cg = compute_cut_graph(mesh);
    auto dt_cut = std::chrono::duration_cast<Ms>(Clock::now() - t0).count();

    std::cout << "[SmokeEuclidean.Brezel2] V=" << mesh.number_of_vertices()
              << " F=" << mesh.number_of_faces()
              << " cut=" << dt_cut << "ms"
              << " seams=" << cg.cut_edge_indices.size() << "\n";

    // Genus-2: 2g = 4 seam edges
    EXPECT_EQ(4u, cg.cut_edge_indices.size())
        << "brezel2.obj (genus 2) must yield 2g=4 cut edges";
    EXPECT_EQ(2, cg.genus);
}