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

// conformallab_cli.cpp
//
// ConformalLab++ command-line interface.
//
// Usage:
//   conformallab_core -i input.off [-o layout.off] [-g euclidean|spherical|hyper_ideal]
//                    [-j result.json] [-x result.xml] [-s] [-v]
//
// The tool:
//   1. Loads an OFF/OBJ/PLY mesh.
//   2. Sets up DOF maps + computes λ° from the input geometry.
//   3. Pins one vertex (Euclidean/Spherical) or uses all-free DOFs (HyperIdeal).
//   4. Runs Newton until convergence.
//   5. Computes a 2-D (Euclidean / HyperIdeal) or 3-D (Spherical) layout.
//   6. Saves the layout as an OFF file and optionally serialises the result
//      to JSON and/or XML.
//   7. Optionally shows the input mesh in a viewer (-s flag, requires WITH_VIEWER).

#include "conformal_mesh.hpp"
#include "mesh_io.hpp"
#include "euclidean_functional.hpp"
#include "spherical_functional.hpp"
#include "hyper_ideal_functional.hpp"
#include "newton_solver.hpp"
#include "layout.hpp"
#include "serialization.hpp"

#include <CLI11.hpp>
#include <iostream>
#include <iomanip>
#include <string>
#include <vector>
#include <cmath>
#include <algorithm>

// Viewer is always available when WITH_CGAL=ON (implies WITH_VIEWER=ON)
#include "viewer_utils.h"
#include <Eigen/Core>

namespace cl = conformallab;
using cl::ConformalMesh;
using cl::Vertex_index;
using cl::Edge_index;

// ─────────────────────────────────────────────────────────────────────────────
// Shared helpers (mirroring test_pipeline.cpp patterns)
// ─────────────────────────────────────────────────────────────────────────────

// Pin vertex 0, assign 0..n-1 to the rest
static int pin_first_vertex(ConformalMesh& mesh, cl::EuclideanMaps& maps)
{
    auto vit = mesh.vertices().begin();
    Vertex_index v0 = *vit++;
    maps.v_idx[v0] = -1;
    int idx = 0;
    for (; vit != mesh.vertices().end(); ++vit)
        maps.v_idx[*vit] = idx++;
    return idx;
}

// Natural theta for Euclidean: make x=0 the equilibrium
static void set_natural_euclidean_theta(ConformalMesh& mesh, cl::EuclideanMaps& maps, int n)
{
    std::vector<double> x0(static_cast<std::size_t>(n), 0.0);
    auto G = cl::euclidean_gradient(mesh, x0, maps);
    for (auto v : mesh.vertices()) {
        int iv = maps.v_idx[v];
        if (iv < 0) continue;
        maps.theta_v[v] -= G[static_cast<std::size_t>(iv)];
    }
}

// Natural theta for HyperIdeal at base point (b=1, a=0.5) to avoid x=0 singularity
static std::vector<double> set_natural_hyper_ideal_targets(
    ConformalMesh& mesh, cl::HyperIdealMaps& maps, int n)
{
    const auto sz = static_cast<std::size_t>(n);
    std::vector<double> xbase(sz, 0.0);
    for (auto v : mesh.vertices()) {
        int iv = maps.v_idx[v];
        if (iv >= 0) xbase[static_cast<std::size_t>(iv)] = 1.0;
    }
    for (auto e : mesh.edges()) {
        int ie = maps.e_idx[e];
        if (ie >= 0) xbase[static_cast<std::size_t>(ie)] = 0.5;
    }
    auto G = cl::evaluate_hyper_ideal(mesh, xbase, maps, false).gradient;
    for (auto v : mesh.vertices()) {
        int iv = maps.v_idx[v];
        if (iv < 0) continue;
        maps.theta_v[v] += G[static_cast<std::size_t>(iv)];
    }
    for (auto e : mesh.edges()) {
        int ie = maps.e_idx[e];
        if (ie < 0) continue;
        maps.theta_e[e] += G[static_cast<std::size_t>(ie)];
    }
    return xbase;
}

// ─────────────────────────────────────────────────────────────────────────────
// Euclidean pipeline
// ─────────────────────────────────────────────────────────────────────────────
static int run_euclidean(ConformalMesh& mesh,
                         const std::string& out_layout,
                         const std::string& out_json,
                         const std::string& out_xml,
                         bool verbose)
{
    // Setup
    auto maps = cl::setup_euclidean_maps(mesh);
    cl::compute_euclidean_lambda0_from_mesh(mesh, maps);

    // DOF assignment: pin vertex 0
    int n = pin_first_vertex(mesh, maps);
    if (n <= 0) { std::cerr << "Error: mesh has only one vertex.\n"; return 1; }

    // Natural target angles
    set_natural_euclidean_theta(mesh, maps, n);

    // Newton
    std::vector<double> x0(static_cast<std::size_t>(n), 0.0);
    auto res = cl::newton_euclidean(mesh, x0, maps);

    if (!res.converged && verbose)
        std::cerr << "[warn] Newton did not converge (|grad|=" << res.grad_inf_norm << ")\n";

    // Layout
    cl::Layout2D layout = cl::euclidean_layout(mesh, res.x, maps);

    // Output
    if (!out_layout.empty()) cl::save_layout_off(out_layout, mesh, layout);
    if (!out_json.empty())
        cl::save_result_json(out_json, res, "euclidean",
                             static_cast<int>(mesh.number_of_vertices()),
                             static_cast<int>(mesh.number_of_faces()),
                             &layout);
    if (!out_xml.empty())
        cl::save_result_xml(out_xml, res, "euclidean",
                            static_cast<int>(mesh.number_of_vertices()),
                            static_cast<int>(mesh.number_of_faces()),
                            &layout);

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "Euclidean: converged=" << (res.converged ? "yes" : "no")
              << "  iter=" << res.iterations
              << "  |grad|_inf=" << std::scientific << std::setprecision(3)
              << res.grad_inf_norm << "\n";
    if (!out_layout.empty()) std::cout << "  layout → " << out_layout << "\n";
    if (!out_json.empty())   std::cout << "  json   → " << out_json   << "\n";
    if (!out_xml.empty())    std::cout << "  xml    → " << out_xml    << "\n";
    return 0;
}

// ─────────────────────────────────────────────────────────────────────────────
// Spherical pipeline
// ─────────────────────────────────────────────────────────────────────────────
static int run_spherical(ConformalMesh& mesh,
                          const std::string& out_layout,
                          const std::string& out_json,
                          const std::string& out_xml,
                          bool verbose)
{
    auto maps = cl::setup_spherical_maps(mesh);
    cl::compute_lambda0_from_mesh(mesh, maps);
    int n = cl::assign_vertex_dof_indices(mesh, maps);

    std::vector<double> x0(static_cast<std::size_t>(n), 0.0);
    auto res = cl::newton_spherical(mesh, x0, maps);

    if (!res.converged && verbose)
        std::cerr << "[warn] Newton did not converge (|grad|=" << res.grad_inf_norm << ")\n";

    cl::Layout3D layout = cl::spherical_layout(mesh, res.x, maps);

    if (!out_layout.empty()) cl::save_layout_off(out_layout, mesh, layout);
    if (!out_json.empty())
        cl::save_result_json(out_json, res, "spherical",
                             static_cast<int>(mesh.number_of_vertices()),
                             static_cast<int>(mesh.number_of_faces()),
                             nullptr, &layout);
    if (!out_xml.empty())
        cl::save_result_xml(out_xml, res, "spherical",
                            static_cast<int>(mesh.number_of_vertices()),
                            static_cast<int>(mesh.number_of_faces()),
                            nullptr, &layout);

    std::cout << "Spherical: converged=" << (res.converged ? "yes" : "no")
              << "  iter=" << res.iterations
              << "  |grad|_inf=" << std::scientific << std::setprecision(3)
              << res.grad_inf_norm << "\n";
    if (!out_layout.empty()) std::cout << "  layout → " << out_layout << "\n";
    if (!out_json.empty())   std::cout << "  json   → " << out_json   << "\n";
    if (!out_xml.empty())    std::cout << "  xml    → " << out_xml    << "\n";
    return 0;
}

// ─────────────────────────────────────────────────────────────────────────────
// HyperIdeal pipeline
// ─────────────────────────────────────────────────────────────────────────────
static int run_hyper_ideal(ConformalMesh& mesh,
                            const std::string& out_layout,
                            const std::string& out_json,
                            const std::string& out_xml,
                            bool verbose)
{
    auto maps = cl::setup_hyper_ideal_maps(mesh);
    int n = cl::assign_all_dof_indices(mesh, maps);

    // Natural targets at base point (b=1, a=0.5)
    auto xbase = set_natural_hyper_ideal_targets(mesh, maps, n);

    // Perturb slightly from the base and solve back
    std::vector<double> x0 = xbase;
    for (auto& v : x0) v += 0.3;

    auto res = cl::newton_hyper_ideal(mesh, x0, maps);

    if (!res.converged && verbose)
        std::cerr << "[warn] Newton did not converge (|grad|=" << res.grad_inf_norm << ")\n";

    cl::Layout2D layout = cl::hyper_ideal_layout(mesh, res.x, maps);

    if (!out_layout.empty()) cl::save_layout_off(out_layout, mesh, layout);
    if (!out_json.empty())
        cl::save_result_json(out_json, res, "hyper_ideal",
                             static_cast<int>(mesh.number_of_vertices()),
                             static_cast<int>(mesh.number_of_faces()),
                             &layout);
    if (!out_xml.empty())
        cl::save_result_xml(out_xml, res, "hyper_ideal",
                            static_cast<int>(mesh.number_of_vertices()),
                            static_cast<int>(mesh.number_of_faces()),
                            &layout);

    std::cout << "HyperIdeal: converged=" << (res.converged ? "yes" : "no")
              << "  iter=" << res.iterations
              << "  |grad|_inf=" << std::scientific << std::setprecision(3)
              << res.grad_inf_norm << "\n";
    if (!out_layout.empty()) std::cout << "  layout → " << out_layout << "\n";
    if (!out_json.empty())   std::cout << "  json   → " << out_json   << "\n";
    if (!out_xml.empty())    std::cout << "  xml    → " << out_xml    << "\n";
    return 0;
}

// ─────────────────────────────────────────────────────────────────────────────
// main
// ─────────────────────────────────────────────────────────────────────────────
int main(int argc, char* argv[])
{
    CLI::App app{"ConformalLab++ — discrete conformal mapping"};

    std::string input_file;
    std::string out_layout;
    std::string out_json;
    std::string out_xml;
    std::string geometry = "euclidean";
    bool show    = false;
    bool verbose = false;

    app.add_option("-i,--input",    input_file, "Input mesh (OFF/OBJ/PLY)")->required();
    app.add_option("-o,--output",   out_layout, "Output layout OFF file");
    app.add_option("-j,--json",     out_json,   "Save result as JSON");
    app.add_option("-x,--xml",      out_xml,    "Save result as XML");
    app.add_option("-g,--geometry", geometry,   "Target geometry: euclidean|spherical|hyper_ideal")
        ->check(CLI::IsMember({"euclidean", "spherical", "hyper_ideal"}));
    app.add_flag("-s,--show",    show,    "Visualise input mesh (requires WITH_VIEWER)");
    app.add_flag("-v,--verbose", verbose, "Verbose output");

    CLI11_PARSE(app, argc, argv);

    // ── Load mesh ─────────────────────────────────────────────────────────────
    ConformalMesh mesh;
    if (!cl::read_mesh(input_file, mesh) || mesh.is_empty()) {
        std::cerr << "Error: cannot load mesh from '" << input_file << "'\n";
        return EXIT_FAILURE;
    }
    std::cout << "Loaded: " << input_file
              << "  (" << mesh.number_of_vertices() << "v, "
              << mesh.number_of_faces() << "f)\n";

    // ── Optional viewer ───────────────────────────────────────────────────────
    if (show) {
        // Convert ConformalMesh to Eigen matrices for the libigl viewer
        const std::size_t nv = mesh.number_of_vertices();
        const std::size_t nf = mesh.number_of_faces();
        Eigen::MatrixXd V(static_cast<Eigen::Index>(nv), 3);
        Eigen::MatrixXi F(static_cast<Eigen::Index>(nf), 3);

        for (auto v : mesh.vertices()) {
            auto p = mesh.point(v);
            V.row(v.idx()) << p.x(), p.y(), p.z();
        }
        Eigen::Index fi = 0;
        for (auto f : mesh.faces()) {
            int ci = 0;
            for (auto h : CGAL::halfedges_around_face(mesh.halfedge(f), mesh))
                F(fi, ci++) = static_cast<int>(mesh.target(h).idx());
            ++fi;
        }
        viewer_utils::simple_visualize(V, F);
    }

    // ── Dispatch ──────────────────────────────────────────────────────────────
    if (geometry == "euclidean")
        return run_euclidean(mesh, out_layout, out_json, out_xml, verbose);
    if (geometry == "spherical")
        return run_spherical(mesh, out_layout, out_json, out_xml, verbose);
    if (geometry == "hyper_ideal")
        return run_hyper_ideal(mesh, out_layout, out_json, out_xml, verbose);

    std::cerr << "Unknown geometry: " << geometry << "\n";
    return EXIT_FAILURE;
}