#include <array>
#include <bit>
#include <cmath>
#include <cstdint>
#include <limits>
#include <string>
#include <tuple>
#include <type_traits>

#include <catch2/catch_test_macros.hpp>
#include <catch2/catch_approx.hpp>

#include "serif/utils/serialize/numerics.hpp"

namespace {
    namespace serialization = serif::utils::serialize;
    namespace blocks = serif::discretization::blocks;
    namespace dependencies = serif::utils::dependencies;
    namespace model = serif::model;
    namespace equations = serif::equations;
    namespace dimensions = serif::dimensions;
    namespace normalization = serif::solver::normalization;
    using EOS = serif::eos::models::Polytrope;
    using Coordinates = dependencies::KeyList<blocks::density::mass, blocks::gravity::potential,
        blocks::surface_deformation::coefficients, blocks::barotropic_constant::mass_normalization>;
    using Rows = dependencies::KeyList<blocks::density::mass, blocks::gravity::potential,
        blocks::surface_deformation::shape_equilibrium, blocks::barotropic_constant::mass_normalization>;
    using Layout = blocks::SystemLayout<Coordinates, Rows>;

    /** @brief Exercise the complete concrete barotropic equation description without a mesh. */
    [[nodiscard]] auto physical_description(const EOS& eos) {
        return model::EquilibriumEquations{
            .interior_equations = equations::EquationRegistry{equations::BarotropicDensityBalance{eos, dimensions::PressureValue{0.25}},
                equations::PrimalPoissonBalance{12.5}, equations::MassNormalization{4.0}},
            .shape_force = equations::SurfacePotentialBalance{},
            .surface_condition = serif::geometry::surface::Isobaric{dimensions::PressureValue{0.25}}
        };
    }
}

TEST_CASE("Serialized stellar physics reconstructs every concrete barotropic component", "[unit][serialize][serialize_components]") {
    auto original = model::make_stellar_model(EOS{1.5, 2.75}, physical_description(EOS{1.5, 2.75}),
        model::FixedTotalMass{dimensions::MassValue{4.0}},
        model::FixedAngularMomentum{{.angular_momentum = dimensions::AngularMomentumValue{0.3}, .axis = {0, 1, 0}, .center = {0.2, -0.4, 0.1}}});
    const auto node = serialization::encode(original);
    const auto restored = serialization::decode<std::remove_cvref_t<decltype(original)>>(node);
    REQUIRE(restored);
    CHECK(serialization::encode(*restored) == node);
    CHECK(restored->equation_of_state().polytropic_constant() == 2.75);
    CHECK(restored->template specification<model::FixedAngularMomentum>().center()[1] == -0.4);

    auto invalid = node;
    invalid.children.at("eos").children.at("index") = serialization::encode(-2.0);
    const auto rejected = serialization::decode<std::remove_cvref_t<decltype(original)>>(invalid);
    REQUIRE_FALSE(rejected);
    CHECK(rejected.error().path == "eos");
}

TEST_CASE("Scalar equilibrium archives retain the density balance, primal gravity and the isobar level", "[unit][serialize][serialize_components]") {
    const EOS eos{1.5, 2.75};
    const std::array<double, 3> axis{0.0, 1.0, 0.0};
    const std::array<double, 3> center{0.2, -0.4, 0.1};
    const auto original = model::make_stellar_model(eos, model::EquilibriumEquations{
        .interior_equations = equations::EquationRegistry{equations::BarotropicDensityBalance{eos, dimensions::PressureValue{0.25}, axis, center},
            equations::PrimalPoissonBalance{12.5}, equations::MassNormalization{4.0},
            equations::AngularMomentumConstraint{0.3, axis, center}},
        .shape_force = equations::SurfacePotentialBalance{axis, center},
        .surface_condition = serif::geometry::surface::Isobaric{dimensions::PressureValue{0.25}}},
        model::FixedTotalMass{dimensions::MassValue{4.0}},
        model::FixedAngularMomentum{{.angular_momentum = dimensions::AngularMomentumValue{0.3}, .axis = axis, .center = center}});
    const auto node = serialization::encode(original);
    const auto restored = serialization::decode<std::remove_cvref_t<decltype(original)>>(node);
    REQUIRE(restored);
    CHECK(serialization::encode(*restored) == node);
    CHECK(restored->equations().interior_equations.template get<equations::PrimalPoissonBalance>().four_pi_G() == 12.5);
    CHECK(restored->equations().shape_force.rotation().axis_origin() == center);
    CHECK(restored->template specification<model::FixedTotalMass>().target().value() == 4.0);
}

TEST_CASE("Barotropic density archives preserve surface pressure and a tilted displaced rotation axis", "[unit][serialize][serialize_components][barotropic_density]") {
    using Balance = equations::BarotropicDensityBalance<EOS>;
    const std::array<double, 3> axis{0.0, 3.0, 4.0}, center{0.2, -0.4, 0.1};
    const Balance original{EOS{2.0, 0.4}, dimensions::PressureValue{0.05}, axis, center};
    const auto node = serialization::encode(original);
    CHECK(node.type == "equation.barotropic_density");
    const auto restored = serialization::decode<Balance>(node);
    REQUIRE(restored);
    CHECK(restored->surface_pressure().value() == 0.05);
    CHECK(restored->equation_of_state().polytropic_index() == 2.0);
    CHECK(restored->equation_of_state().polytropic_constant() == 0.4);
    CHECK(restored->rotation().axis_origin() == center);

    equations::EquationPoint<Balance::Operands> point;
    point.fields.get<Balance::Density>() = 1.44;
    point.fields.get<Balance::Potential>() = -0.7;
    point.fields.get<Balance::SurfaceLevel>() = -0.9;
    point.fields.get<Balance::AngularVelocity>() = 0.3;
    point.position = {1.2, 0.7, -0.2};
    const auto expected = original.evaluate(point), actual = restored->evaluate(point);
    // For n=2, K=0.4, P_surface=0.05: rho_surface=0.25 and A(rho)-A_surface=0.84.
    CHECK(actual.residual == Catch::Approx(1.04 - 0.09 * original.rotation().centrifugal_potential(point.position)).margin(2.0e-15));
    CHECK(actual.residual == Catch::Approx(expected.residual).margin(2.0e-15));
    for (std::size_t index = 0; index < actual.field_derivatives.values.size(); ++index) {
        CHECK(actual.field_derivatives.values[index] == Catch::Approx(expected.field_derivatives.values[index]).margin(2.0e-15));
    }
    for (int component = 0; component < 3; ++component) {
        CHECK(restored->rotation().angular_velocity_vector()[component] == Catch::Approx(original.rotation().angular_velocity_vector()[component]).margin(2.0e-15));
        CHECK(actual.position_derivative[component] == Catch::Approx(expected.position_derivative[component]).margin(2.0e-15));
    }

    // Deserialization must own formerly borrowed physics and reconstruct its finite-pressure reference.
    const auto borrowed_record = [&]() {
        const EOS temporary{2.0, 0.4};
        return serialization::encode(equations::BarotropicDensityBalance{model::borrow_eos(temporary), dimensions::PressureValue{0.05}, axis, center});
    }();
    using Borrowed = equations::BarotropicDensityBalance<model::EosReference<EOS>>;
    const auto reowned = serialization::decode<Borrowed>(borrowed_record);
    REQUIRE(reowned);
    CHECK(reowned->equation_of_state().get().polytropic_constant() == 0.4);
    CHECK(reowned->evaluate(point).residual == Catch::Approx(expected.residual).margin(2.0e-15));
    CHECK_FALSE(serialization::decode<equations::SurfacePotentialBalance>(node));
    auto invalid = node;
    invalid.children.at("surface_pressure") = serialization::encode(dimensions::PressureValue{-0.1});
    CHECK_FALSE(serialization::decode<Balance>(invalid));
    invalid = node;
    invalid.children.at("axis") = serialization::encode(std::array<double, 3>{0.0, 0.0, 0.0});
    CHECK_FALSE(serialization::decode<Balance>(invalid));
}

TEST_CASE("Deferred model archives restore current physics without pretending to recreate a callback", "[unit][serialize][serialize_components]") {
    const double coefficient = 2.0;
    auto original = model::make_stellar_model(EOS{1.0, 1.0}, [coefficient](const auto& star) {
        return model::EquilibriumEquations{equations::EquationRegistry{equations::PrimalPoissonBalance{coefficient * star.equation_of_state().polytropic_index()}}};
    });
    const auto node = serialization::encode(original);
    CHECK(node.children.at("authored_definition").type == "unsupported");
    CHECK_FALSE(serialization::decode<std::remove_cvref_t<decltype(original)>>(node));
    const auto restored = serialization::restore_stellar_model<decltype(original)>(node);
    REQUIRE(restored);
    CHECK(restored->equations().interior_equations.template get<equations::PrimalPoissonBalance>().four_pi_G() == coefficient);
    CHECK(serialization::encode(restored->equations()) == node.children.at("resolved_equations"));
}

TEST_CASE("Physical state serialization preserves coefficients and rejects mismatched coordinate schemas", "[unit][serialize][serialize_components]") {
    const Layout layout{{2, 3, 2, 1}, {2, 3, 2, 1}};
    std::vector<double> coefficients{0.0, -0.0, 1.0e-250, -1.0e200, 2.0, 0.2, -0.3, 1.5};
    const model::StellarState original{layout, coefficients};
    auto node = serialization::encode(original);
    const auto restored = serialization::decode<std::remove_cvref_t<decltype(original)>>(node);
    REQUIRE(restored);
    REQUIRE(restored->values().size() == coefficients.size());
    for (std::size_t index = 0; index < coefficients.size(); ++index) {
        CHECK(std::bit_cast<std::uint64_t>(restored->values()[index]) == std::bit_cast<std::uint64_t>(coefficients[index]));
    }
    node.children.at("layout").children.at("values").children.at("0").children.at("id") = serialization::encode(std::string{"gravity.potential"});
    CHECK_FALSE(serialization::decode<std::remove_cvref_t<decltype(original)>>(node));
    node = serialization::encode(original);
    node.children.at("coefficients") = serialization::encode(std::vector<double>{1.0});
    CHECK_FALSE(serialization::decode<std::remove_cvref_t<decltype(original)>>(node));
}

TEST_CASE("Default numerical recipes retain resolved parameters rather than selecting fresh defaults", "[unit][serialize][serialize_components]") {
    using PrimalValues = dependencies::KeyList<blocks::density::mass, blocks::gravity::potential,
        blocks::barotropic_constant::mass_normalization, blocks::rotation::angular_velocity, blocks::surface_deformation::coefficients>;
    using PrimalRows = dependencies::KeyList<blocks::density::mass, blocks::gravity::potential,
        blocks::barotropic_constant::mass_normalization, blocks::rotation::angular_velocity, blocks::surface_deformation::shape_equilibrium>;
    const blocks::SystemLayout<PrimalValues, PrimalRows> primal_layout{{2, 3, 1, 1, 2}, {2, 3, 1, 1, 2}};
    const auto primal_inverse = serif::solver::preconditioner::make_primal_preconditioner(primal_layout);
    const auto restored_primal = serialization::decode<std::remove_cvref_t<decltype(primal_inverse)>>(serialization::encode(primal_inverse));
    REQUIRE(restored_primal);
    CHECK(serialization::encode(*restored_primal) == serialization::encode(primal_inverse));

    auto options = serif::solver::SolverOptions<>{};
    options.method.maximum_iterations = 123;
    options.linear_solver.restart_length = 31;
    options.line_search.contraction = 0.4;
    const auto restored_options = serialization::decode<decltype(options)>(serialization::encode(options));
    REQUIRE(restored_options);
    CHECK(serialization::encode(*restored_options) == serialization::encode(options));

    const auto numerical = normalization::make_normalization(normalization::NormalizationOptions{
        .scheme = normalization::LaneEmdenApproximation{.n = 2.5, .mass = dimensions::MassValue{4.}, .central_density = dimensions::DensityValue{2.}},
        .coordinates = normalization::Dimensionless{}, .gravitational_constant = dimensions::QuantityValue<dimensions::GravitationalConstant>{0.75}});
    const auto restored_numerical = serialization::decode<std::remove_cvref_t<decltype(numerical)>>(serialization::encode(numerical));
    REQUIRE(restored_numerical);
    CHECK(serialization::encode(*restored_numerical) == serialization::encode(numerical));
    const auto factors = normalization::make_diagonal_normalization({2., 4.}, {0.1, 0.3});
    REQUIRE(factors);
    auto saved_factors = serialization::encode(*factors);
    REQUIRE(serialization::decode<normalization::DiagonalNormalization>(saved_factors));
    saved_factors.children.at("state") = serialization::encode(std::vector<double>{-1., 4.});
    CHECK_FALSE(serialization::decode<normalization::DiagonalNormalization>(saved_factors));
}

TEST_CASE("Surface representations retain physical centers and harmonic selection", "[unit][serialize][serialize_components]") {
    using namespace serif::geometry::deformation;
    const DomainDeformation original{SphericalHarmonicSurface{.reference_center = {0.1, -0.2, 0.3}, .maximum_degree = 8, .maximum_azimuthal_order = 3},
        PowerLawRadialExtension{.radial_power = 1.75}, FixedInfinityRadialExtension{}};
    const auto restored = serialization::decode<std::remove_cvref_t<decltype(original)>>(serialization::encode(original));
    REQUIRE(restored);
    CHECK(*restored == original);
    const auto nodal = serialization::decode<NodalRadialSurface>(serialization::encode(NodalRadialSurface{{1., 2., 3.}}));
    REQUIRE(nodal);
    CHECK(nodal->reference_center[2] == 3.0);
}

TEST_CASE("Quantity serializers reject incorrect dimensions and unknown schema versions", "[unit][serialize][serialize_components]") {
    auto node = serialization::encode(dimensions::MassValue{7.0});
    CHECK_FALSE(serialization::decode<dimensions::DensityValue>(node));
    node.version = 17;
    CHECK_FALSE(serialization::decode<dimensions::MassValue>(node));
}

TEST_CASE("Quadrature role controls survive serialization without compiled caches", "[unit][serialize][serialize_components]") {
    namespace quadrature = serif::discretization::quadrature;
    quadrature::RuleSet rules;
    rules.set_role_control(quadrature::QuadratureRole::projection, {.fixed_order = 12, .boost = 2});
    rules.set_fallback_control({.boost = 3});
    const quadrature::Policy original{rules};
    const auto restored = serialization::decode<quadrature::Policy>(serialization::encode(original));
    REQUIRE(restored);
    CHECK(restored->rules().role_control(quadrature::QuadratureRole::projection).fixed_order == 12);
    CHECK(restored->rules().fallback_control().boost == 3);
}

TEST_CASE("Prescribed input serialization distinguishes owned values from borrowed runtime objects", "[unit][serialize][serialize_components]") {
    const auto inputs = dependencies::make_dependency_context(dependencies::own<blocks::rotation::angular_velocity>(0.3));
    const auto restored = serialization::decode<std::remove_cvref_t<decltype(inputs)>>(serialization::encode(inputs));
    REQUIRE(restored);
    CHECK(restored->template get<blocks::rotation::angular_velocity>() == 0.3);
    double prescribed = 0.4;
    const auto borrowed = dependencies::make_dependency_context(dependencies::borrow<blocks::rotation::angular_velocity>(prescribed));
    const auto node = serialization::encode(borrowed);
    CHECK(node.children.at("resources").children.at("0").children.at("reconstruction").type == "unsupported");
    CHECK_FALSE(serialization::decode<std::remove_cvref_t<decltype(borrowed)>>(node));
}

TEST_CASE("Decoded numerical options reject invalid budgets and approximation controls", "[unit][serialize][serialize_components]") {
    auto options = serialization::encode(serif::solver::SolverOptions<>{});
    options.children.at("method").children.at("maximum_iterations") = serialization::encode(-3);
    const auto invalid_solver = serialization::decode<serif::solver::SolverOptions<>>(options);
    REQUIRE_FALSE(invalid_solver);
    CHECK(invalid_solver.error().path == "method");

    using AMG = serif::solver::preconditioner::strategies::inverse::FixedAMG;
    CHECK_FALSE(serialization::decode<AMG>(serialization::encode(AMG{.cycles = 0})));
    using Surface = serif::geometry::deformation::SphericalHarmonicSurface;
    CHECK_FALSE(serialization::decode<Surface>(serialization::encode(Surface{.maximum_degree = 2, .maximum_azimuthal_order = 3})));
    using Integration = serif::utils::laneemden::IntegrationOptions;
    CHECK_FALSE(serialization::decode<Integration>(serialization::encode(Integration{.maximum_steps = 0})));
}

TEST_CASE("Custom quadrature form controls remain inspectable when their descriptor cannot be reowned", "[unit][serialize][serialize_components]") {
    namespace quadrature = serif::discretization::quadrature;
    quadrature::RuleSet rules;
    rules.set_form_control({.name = "custom_mass", .trial = {.field = {"density"}, .name = "mass"}, .test = {.field = {"density"}, .name = "mass"}}, {.boost = 7});
    const auto node = serialization::encode(rules);
    CHECK(node.children.at("forms").children.at("reconstruction").type == "unsupported");
    const auto boost = serialization::decode<int>(node.children.at("forms").children.at("0").children.at("control").children.at("boost"));
    REQUIRE(boost);
    CHECK(*boost == 7);
    CHECK_FALSE(serialization::decode<quadrature::RuleSet>(node));
}

TEST_CASE("Configured discretizations preserve owned normalization inputs and resolved choices", "[unit][serialize][serialize_components]") {
    auto star = model::make_stellar_model(EOS{1.5, 2.75}, physical_description(EOS{1.5, 2.75}),
        model::FixedTotalMass{dimensions::MassValue{4.0}});
    auto numerical = normalization::NormalizationOptions{
        .scheme = normalization::LaneEmdenApproximation<normalization::FromModel, dimensions::DensityValue>{.central_density = dimensions::DensityValue{2.0}},
        .sources = dependencies::make_dependency_context(dependencies::own<blocks::rotation::angular_velocity>(0.3))};
    auto original = model::make_discretization(model::DiscretizationOptions{
        .model = std::move(star), .finite_elements = std::vector<double>{0.0, 0.5, 1.0}, .normalization = std::move(numerical)});
    const auto node = serialization::encode(original);
    CHECK(node.type == "model.configured_discretization");
    const auto restored = serialization::decode<decltype(original)>(node);
    REQUIRE(restored);
    CHECK(serialization::encode(*restored) == node);
    CHECK(restored->normalization.scheme.mass == dimensions::MassValue{4.0});
    CHECK(restored->normalization.sources.template get<blocks::rotation::angular_velocity>() == 0.3);
}