// example_layout.cpp
//
// Phase 5 example — Euclidean conformal layout with JSON/XML serialisation.
//
// Demonstrates the full pipeline that a library user would run:
//
//   1. Build (or load) a mesh.
//   2. Set up Euclidean conformal maps + compute λ° from vertex positions.
//   3. Choose natural target angles (→ x* = 0 is the equilibrium).
//   4. Run Newton's method.
//   5. Unfold the mesh in the plane (euclidean_layout).
//   6. Save the layout as an OFF file for inspection in MeshLab/Blender.
//   7. Serialise the solver result and UV coordinates to JSON and XML.
//   8. Reload from JSON and verify the DOF vector is recovered.
//
// Build (from code/):
//   cmake -S . -B build -DWITH_CGAL=ON && cmake --build build -t example_layout
//
// Run:
//   ./build/examples/example_layout
//   ./build/examples/example_layout input.off layout.off result.json result.xml

#include "conformal_mesh.hpp"
#include "mesh_builder.hpp"
#include "mesh_io.hpp"
#include "euclidean_functional.hpp"
#include "newton_solver.hpp"
#include "layout.hpp"
#include "serialization.hpp"

#include <iostream>
#include <iomanip>
#include <vector>
#include <string>
#include <cmath>

using namespace conformallab;

// ── Helper: natural target angles so x* = 0 is the equilibrium ───────────────
//
// ⚠  TESTING CONVENTION — NOT A REAL CONFORMAL MAP
// Natural theta sets Θ_v = actual angle sum at x=0, making x* = 0 trivially
// the equilibrium (u_v ≈ 0, no deformation).  For real conformal flattening,
// see example_flatten.cpp which uses Θ_v = 2π (flat interior target).
static void set_natural_theta(ConformalMesh& mesh, EuclideanMaps& maps, int n)
{
    std::vector<double> x0(static_cast<std::size_t>(n), 0.0);
    auto G = 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)];
    }
}

// ── Helper: pin vertex 0, assign DOF indices 0..n-1 to the rest ──────────────
// Uses the gauge-vertex overload introduced in external-audit Finding-D:
// assign_euclidean_vertex_dof_indices(mesh, maps, gauge) pins the chosen
// vertex and assigns sequential indices in a single pass.
static int pin_first(ConformalMesh& mesh, EuclideanMaps& maps)
{
    return assign_euclidean_vertex_dof_indices(
        mesh, maps, *mesh.vertices().begin());
}

// ── Main ─────────────────────────────────────────────────────────────────────

int main(int argc, char* argv[])
{
    // ── Step 1: mesh ──────────────────────────────────────────────────────────
    ConformalMesh mesh;
    std::string out_off  = "layout.off";
    std::string out_json = "result.json";
    std::string out_xml  = "result.xml";

    if (argc >= 2) {
        // User provided an input mesh
        if (!read_mesh(argv[1], mesh) || mesh.is_empty()) {
            std::cerr << "Cannot load " << argv[1] << "\n";
            return 1;
        }
        if (argc >= 3) out_off  = argv[2];
        if (argc >= 4) out_json = argv[3];
        if (argc >= 5) out_xml  = argv[4];
    } else {
        // Use built-in quad-strip (6 vertices, 4 triangles)
        mesh = make_quad_strip();
        std::cout << "Using built-in quad-strip mesh (no arguments given).\n";
    }

    std::cout << "Mesh: " << mesh.number_of_vertices() << " vertices, "
              << mesh.number_of_faces() << " faces\n";

    // ── Step 2: maps + λ° ────────────────────────────────────────────────────
    auto maps = setup_euclidean_maps(mesh);
    compute_euclidean_lambda0_from_mesh(mesh, maps);

    // ── Step 3: DOF assignment + natural targets ──────────────────────────────
    int n = pin_first(mesh, maps);
    std::cout << "DOFs: " << n << " (vertex 0 pinned)\n";

    set_natural_theta(mesh, maps, n);

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

    std::cout << std::boolalpha
              << "Newton converged: " << res.converged
              << "  iterations: "     << res.iterations
              << "  |grad|_inf: "
              << std::scientific << std::setprecision(3) << res.grad_inf_norm << "\n";

    // Print scale factors u_v
    std::cout << "Conformal factors u_v:\n";
    for (auto v : mesh.vertices()) {
        int iv = maps.v_idx[v];
        double u = (iv >= 0) ? res.x[static_cast<std::size_t>(iv)] : 0.0;
        std::cout << "  v" << v.idx() << ": u = " << std::fixed << std::setprecision(8) << u << "\n";
    }

    // ── Step 5: Layout ────────────────────────────────────────────────────────
    auto layout = euclidean_layout(mesh, res.x, maps);

    std::cout << "Layout success: " << layout.success
              << "  has_seam: " << layout.has_seam << "\n";

    if (layout.success) {
        std::cout << "UV coordinates:\n";
        for (auto v : mesh.vertices()) {
            auto& p = layout.uv[v.idx()];
            std::cout << "  v" << v.idx()
                      << ": (" << std::fixed << std::setprecision(6)
                      << p.x() << ", " << p.y() << ")\n";
        }
    }

    // ── Step 6: Save layout OFF ───────────────────────────────────────────────
    save_layout_off(out_off, mesh, layout);
    std::cout << "Layout saved → " << out_off << "\n";

    // ── Step 7: Serialise JSON + XML ──────────────────────────────────────────
    save_result_json(out_json, res, "euclidean",
                     static_cast<int>(mesh.number_of_vertices()),
                     static_cast<int>(mesh.number_of_faces()),
                     &layout);
    std::cout << "JSON   saved → " << out_json << "\n";

    save_result_xml(out_xml, res, "euclidean",
                    static_cast<int>(mesh.number_of_vertices()),
                    static_cast<int>(mesh.number_of_faces()),
                    &layout);
    std::cout << "XML    saved → " << out_xml << "\n";

    // ── Step 8: JSON round-trip verification ──────────────────────────────────
    NewtonResult res2;
    std::string  geom;
    Layout2D     layout2;
    auto x2 = load_result_json(out_json, &res2, &geom, &layout2);

    std::cout << "Round-trip: geometry=\"" << geom << "\"  "
              << "DOFs=" << x2.size() << "  "
              << "converged=" << res2.converged << "\n";

    // Verify the DOF vector survived the round-trip
    bool ok = (x2.size() == res.x.size());
    for (std::size_t i = 0; i < x2.size() && ok; ++i)
        ok = (std::abs(x2[i] - res.x[i]) < 1e-10);
    std::cout << "DOF vector round-trip: " << (ok ? "OK ✓" : "MISMATCH ✗") << "\n";

    return res.converged ? 0 : 1;
}