#pragma once

// Include after import mean_field. Root splitting removes the discontinuous
// sign indicator from the integration coordinate; tangential and mapped-volume
// integration still require the quadrature convergence study recorded here.
#include "polytrope_physical_state.hpp"

#include <algorithm>
#include <array>
#include <cmath>
#include <cstring>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <limits>
#include <string>
#include <vector>

namespace experiment::plateau_diagnostics {
namespace negative_density_regions_detail {

using Real = long double;

struct Polynomial final {
    Real c, b, a;
    [[nodiscard]] Real Value(Real t) const { return std::fma(std::fma(a, t, b), t, c); }
    [[nodiscard]] Real Integral(Real left, Real right) const {
        const Real width = right - left;
        // Subtracting antiderivatives loses a thin pocket's entire mass when
        // both antiderivatives are large. This centered quadratic identity does
        // not subtract nearly equal endpoint values.
        return width * (Value(left + width / 2) + a * width * width / 12);
    }
    [[nodiscard]] Real Scale() const { return std::max({std::abs(a), std::abs(b), std::abs(c)}); }
};

struct Interval final { Real left, right; bool negative; };

inline Polynomial Interpolate(Real left, Real middle, Real right) {
    return {left, 4 * middle - 3 * left - right, 2 * (left + right - 2 * middle)};
}

inline std::vector<Interval> Partition(const Polynomial& p) {
    if (!std::isfinite(p.a) || !std::isfinite(p.b) || !std::isfinite(p.c))
        throw std::domain_error("Nonfinite density slice polynomial.");
    std::vector<Real> points{0, 1};
    const auto add = [&](Real root) {
        if (std::isfinite(root) && root > 0 && root < 1) points.push_back(root);
    };
    const Real scale = p.Scale();
    if (scale > 0) {
        const Real a = p.a / scale, b = p.b / scale, c = p.c / scale;
        // No relative cutoff turns a small quadratic into a line: a thin
        // negative interval can be the very effect this diagnostic must retain.
        if (a == 0) {
            if (b != 0) add(-c / b);
        } else {
            const Real discriminant = std::fma(-4 * a, c, b * b);
            if (discriminant == 0) {
                add(-b / (2 * a));
            } else if (discriminant > 0) {
                const Real q = -Real{0.5} * (b + std::copysign(std::sqrt(discriminant), b));
                // The product-of-roots form avoids cancellation at a root near
                // the face. q cannot vanish for a positive discriminant.
                add(q / a);
                add(c / q);
            }
        }
    }
    std::sort(points.begin(), points.end());
    points.erase(std::unique(points.begin(), points.end()), points.end());
    std::vector<Interval> result;
    for (std::size_t i = 1; i < points.size(); ++i) {
        const Real left = points[i - 1], right = points[i];
        if (!(right > left)) throw std::logic_error("Invalid density root partition.");
        const Real value = p.Value(left + (right - left) / 2);
        const bool negative = value < 0;
        // Two extra interior signs catch missed roots without assigning a
        // finite-width pocket to an exactly tangent double root.
        for (const Real fraction : {Real{0.25}, Real{0.75}}) {
            const Real check = p.Value(left + fraction * (right - left));
            if ((negative && check > 0) || (!negative && check < 0))
                throw std::runtime_error("Density slice root partition has inconsistent interval signs.");
        }
        result.push_back({left, right, negative});
    }
    return result;
}

inline std::ofstream Open(const std::filesystem::path& path) {
    if (std::filesystem::exists(path)) throw std::invalid_argument("Density region output already exists: " + path.string());
    std::ofstream out(path);
    out.exceptions(std::ios::failbit | std::ios::badbit);
    out << std::setprecision(20);
    return out;
}

// These are analytic, uniform-measure one-dimensional controls, independent of
// FEM and geometry. The fractional/cubic data are interpolated at three Gauss
// nodes, reproducing the underresolved positive-profile mechanism of interest.
inline void SelfChecks(const std::filesystem::path& path) {
    auto out = Open(path);
    out << "case,negative_length,expected_negative_length,negative_mass,expected_negative_mass,passed\n";
    const auto check = [&](const char* name, const Polynomial& p, Real expectedLength, Real expectedMass) {
        Real length = 0, mass = 0;
        for (const auto& interval : Partition(p)) if (interval.negative) {
            length += interval.right - interval.left;
            mass -= p.Integral(interval.left, interval.right);
        }
        const bool passed = std::abs(length - expectedLength) <= 2e-14L * std::max(Real{1e-12}, expectedLength) &&
            std::abs(mass - expectedMass) <= 2e-13L * std::max(Real{1e-12}, expectedMass);
        out << name << ',' << length << ',' << expectedLength << ',' << mass << ',' << expectedMass << ',' << passed << '\n';
        if (!passed) throw std::runtime_error(std::string("Density region analytic self-check failed: ") + name);
    };
    check("positive_quadratic", {1, 0, 1}, 0, 0);
    check("negative_constant", {-2, 0, 0}, 1, 2);
    check("linear_crossing", {-0.25L, 1, 0}, 0.25L, 0.03125L);
    check("double_root_nonnegative", {0.25L, -1, 1}, 0, 0);
    check("near_linear_positive", {1, 1, 1e-28L}, 0, 0);
    // Exact binary coefficients isolate the root algorithm from the separate
    // conditioning of subtracting a tiny decimal offset from a large value.
    check("near_double_negative", {0.25L - std::ldexp(Real{1}, -40), -1, 1},
          std::ldexp(Real{1}, -19), 4 * std::ldexp(Real{1}, -60) / 3);
    const Real d = std::sqrt(Real{15}) / 10;
    const auto fractional = [](Real t) { return t * std::sqrt(t); };
    const Real y0 = fractional(Real{0.5} - d), y1 = fractional(Real{0.5}), y2 = fractional(Real{0.5} + d);
    const Real a = (y0 + y2 - 2 * y1) / (2 * d * d);
    const Real b = (y2 - y0) / (2 * d) - a;
    const Polynomial fractionalInterpolant{y1 - b / 2 - a / 4, b, a};
    check("gauss_q2_t_power_1p5", fractionalInterpolant,
          0.04548403175047501373888L, 0.0005262295634284655L);
    check("gauss_q2_t_power_3", {0.05L, -0.6L, 1.5L}, std::sqrt(Real{6}) / 15, std::sqrt(Real{6}) / 2250);
}

struct Axis final {
    int coordinate{0};
    int orientation{1};
    int face{-1};
    int faceCount{0};
};

inline std::vector<Axis> ChooseAxes(polytrope_validation::PhysicalState& state) {
    auto& mesh = *state.finiteElements.mesh;
    std::vector<Axis> axes(static_cast<std::size_t>(mesh.GetNE()));
    for (int face = 0; face < mesh.GetNumFaces(); ++face) {
        int left = -1, right = -1;
        mesh.GetFaceElements(face, &left, &right);
        const bool leftStar = left >= 0 && state.isStellar(left);
        const bool rightStar = right >= 0 && state.isStellar(right);
        if (leftStar == rightStar) continue;
        const int element = leftStar ? left : right;
        auto& axis = axes[static_cast<std::size_t>(element)];
        ++axis.faceCount;
        if (axis.face >= 0) continue;
        auto* transformation = mesh.GetFaceElementTransformations(face);
        if (!transformation) throw std::runtime_error("Missing density region face transformation.");
        mfem::IntegrationPoint facePoint, point;
        facePoint.Set2(0.5, 0.5);
        if (leftStar) transformation->Loc1.Transform(facePoint, point);
        else transformation->Loc2.Transform(facePoint, point);
        const std::array<double, 3> coordinates{point.x, point.y, point.z};
        int selected = -1;
        for (int d = 0; d < 3; ++d) {
            if (std::abs(coordinates[d]) < 1e-12 || std::abs(coordinates[d] - 1) < 1e-12) {
                if (selected >= 0) throw std::runtime_error("Ambiguous density region reference face.");
                selected = d;
            }
        }
        if (selected < 0) throw std::runtime_error("Density region face is not a tensor-product cell face.");
        axis.coordinate = selected;
        axis.orientation = coordinates[selected] < 0.5 ? 1 : -1;
        axis.face = face;
    }
    return axes;
}

inline mfem::IntegrationPoint Point(const Axis& axis, double u, double v, Real t) {
    std::array<double, 3> coordinate{};
    coordinate[axis.coordinate] = static_cast<double>(axis.orientation > 0 ? t : 1 - t);
    coordinate[(axis.coordinate + 1) % 3] = u;
    coordinate[(axis.coordinate + 2) % 3] = v;
    mfem::IntegrationPoint result;
    result.Set3(coordinate[0], coordinate[1], coordinate[2]);
    return result;
}

struct Measure final {
    Real volume{0}, mass{0}, negativeVolume{0}, negativeMass{0};
    std::size_t lines{0}, negativeIntervals{0}, points{0};
    double minimumRho{std::numeric_limits<double>::infinity()};
    double witnessH{std::numeric_limits<double>::quiet_NaN()};
    double witnessT{std::numeric_limits<double>::quiet_NaN()}, witnessU{witnessT}, witnessV{witnessT};
    std::array<double, 3> witnessPosition{};
    int witnessElement{-1};
    double negativeRadiusMin{std::numeric_limits<double>::infinity()}, negativeRadiusMax{-negativeRadiusMin};
    Real negativeTMin{std::numeric_limits<Real>::infinity()}, negativeTMax{-negativeTMin};
    double maximumReconstructionError{0}, maximumReconstructionTolerance{0};

    void Add(const Measure& other) {
        volume += other.volume; mass += other.mass;
        negativeVolume += other.negativeVolume; negativeMass += other.negativeMass;
        lines += other.lines; negativeIntervals += other.negativeIntervals; points += other.points;
        if (other.minimumRho < minimumRho) {
            minimumRho = other.minimumRho; witnessH = other.witnessH;
            witnessT = other.witnessT; witnessU = other.witnessU; witnessV = other.witnessV;
            witnessPosition = other.witnessPosition; witnessElement = other.witnessElement;
        }
        negativeRadiusMin = std::min(negativeRadiusMin, other.negativeRadiusMin);
        negativeRadiusMax = std::max(negativeRadiusMax, other.negativeRadiusMax);
        negativeTMin = std::min(negativeTMin, other.negativeTMin);
        negativeTMax = std::max(negativeTMax, other.negativeTMax);
        maximumReconstructionError = std::max(maximumReconstructionError, other.maximumReconstructionError);
        maximumReconstructionTolerance = std::max(maximumReconstructionTolerance, other.maximumReconstructionTolerance);
    }
};

inline Real Observed(Real value) { return std::isfinite(value) ? value : std::numeric_limits<Real>::quiet_NaN(); }

inline void Header(std::ostream& out) {
    out << "tangent_lines,negative_intervals,integration_points,volume,mass,negative_volume,negative_mass,"
           "negative_volume_fraction,negative_mass_fraction,minimum_density,witness_h,witness_element,witness_t,witness_u,witness_v,"
           "witness_x,witness_y,witness_z,negative_radius_minimum,negative_radius_maximum,negative_interval_t_minimum,"
           "negative_interval_t_maximum,maximum_reconstruction_error,maximum_reconstruction_tolerance\n";
}

inline void Write(std::ostream& out, const Measure& m) {
    out << m.lines << ',' << m.negativeIntervals << ',' << m.points << ',' << m.volume << ',' << m.mass << ','
        << m.negativeVolume << ',' << m.negativeMass << ',' << m.negativeVolume / m.volume << ','
        << (m.mass != 0 ? m.negativeMass / std::abs(m.mass) : std::numeric_limits<Real>::quiet_NaN()) << ','
        << m.minimumRho << ',' << m.witnessH << ',' << m.witnessElement << ',' << m.witnessT << ',' << m.witnessU << ',' << m.witnessV;
    for (double x : m.witnessPosition) out << ',' << x;
    out << ',' << Observed(m.negativeRadiusMin) << ',' << Observed(m.negativeRadiusMax) << ','
        << Observed(m.negativeTMin) << ',' << Observed(m.negativeTMax) << ','
        << m.maximumReconstructionError << ',' << m.maximumReconstructionTolerance << '\n';
}

inline bool Identical(const mfem::Vector& a, const mfem::Vector& b) {
    return a.Size() == b.Size() && (a.Size() == 0 || std::memcmp(a.GetData(), b.GetData(), sizeof(mfem::real_t) * a.Size()) == 0);
}

} // namespace negative_density_regions_detail

inline void RecordDensityNegativeRegions(polytrope_validation::PhysicalState& state, double n, double K,
                                         const std::filesystem::path& output) {
    namespace d = negative_density_regions_detail;
    using d::Real;
    if (!std::isfinite(n) || n < 1 || !std::isfinite(K) || K <= 0)
        throw std::invalid_argument("Density regions require finite n>=1 and K>0.");
    std::filesystem::create_directories(output);
    // Fail the one-dimensional arithmetic controls before any FEM evaluation.
    d::SelfChecks(output / "density_negative_regions_selfchecks.csv");
    const auto& fem = state.finiteElements;
    if (!fem.mesh || !fem.densityFes || !fem.enthalpyFes || !fem.compactificationCoordinate ||
        fem.mesh->Dimension() != 3 || !fem.mesh->Conforming() || fem.densityFes->GetVDim() != 1 ||
        dynamic_cast<const mfem::L2_FECollection*>(fem.densityFes->FEColl()) == nullptr ||
        state.density.ParFESpace() != fem.densityFes.get() || state.enthalpy.ParFESpace() != fem.enthalpyFes.get())
        throw std::invalid_argument("Density regions require the current scalar L2 density space on a conforming 3D mesh.");
    int ranks = 0;
    MPI_Comm_size(fem.mesh->GetComm(), &ranks);
    if (ranks != 1) throw std::invalid_argument("Density regions currently require one MPI rank.");
    for (int element = 0; element < fem.mesh->GetNE(); ++element) if (state.isStellar(element)) {
        const auto* fe = fem.densityFes->GetFE(element);
        if (!fe || fe->GetOrder() != 2 || fe->GetMapType() != mfem::FiniteElement::VALUE ||
            fe->GetGeomType() != mfem::Geometry::CUBE || fe->GetDof() != 27)
            throw std::invalid_argument("Density root splitting requires scalar VALUE Q2 hexahedral density elements.");
    }
    const auto axes = d::ChooseAxes(state);
    const mfem::Vector savedDensity(state.density), savedEnthalpy(state.enthalpy), savedPotential(state.potential),
        savedGravity(state.gravityGradientReference), savedDisplacement(state.displacement);
    mean_field::mapping::GridFunctionMappingEvaluator mapper(state.domainMapper, state.displacement, *fem.compactificationCoordinate);
    mean_field::mapping::MappingPointContext mapped;
    auto elements = d::Open(output / "density_negative_regions_elements.csv");
    elements << "tangential_order,interval_order,element,attribute,axis,orientation,surface_face,surface_face_count,";
    d::Header(elements);
    auto summary = d::Open(output / "density_negative_regions_summary.csv");
    summary << "tangential_order,interval_order,";
    d::Header(summary);
    auto intervals = d::Open(output / "density_negative_regions_intervals.csv");
    intervals << "tangential_order,element,attribute,axis,orientation,surface_face,u,v,tangential_weight,t_left,t_right,"
                 "coefficient_c,coefficient_b,coefficient_a,negative_volume,negative_mass,minimum_density,witness_h,"
                 "radius_minimum,radius_maximum,left_x,left_y,left_z,right_x,right_y,right_z\n";
    constexpr int intervalOrder = 14;
    const auto& intervalRule = mfem::IntRules.Get(mfem::Geometry::SEGMENT, intervalOrder);
    if (intervalRule.GetNPoints() != 8) throw std::logic_error("Unexpected density interval Gauss rule.");
    for (int tangentialOrder : {14, 22, 30}) {
        const auto& tangentRule = mfem::IntRules.Get(mfem::Geometry::SEGMENT, tangentialOrder);
        if (tangentRule.GetNPoints() != tangentialOrder / 2 + 1)
            throw std::logic_error("Unexpected density tangential Gauss rule.");
        d::Measure global;
        for (int element = 0; element < fem.mesh->GetNE(); ++element) {
            if (!state.isStellar(element)) continue;
            const auto& axis = axes[static_cast<std::size_t>(element)];
            auto* transformation = fem.mesh->GetElementTransformation(element);
            d::Measure local;
            for (int iu = 0; iu < tangentRule.GetNPoints(); ++iu) for (int iv = 0; iv < tangentRule.GetNPoints(); ++iv) {
                const auto& up = tangentRule.IntPoint(iu);
                const auto& vp = tangentRule.IntPoint(iv);
                const double u = up.x, v = vp.x;
                const Real tangentWeight = static_cast<Real>(up.weight) * vp.weight;
                const auto value = [&](Real t) {
                    const auto point = d::Point(axis, u, v, t);
                    const double rho = state.density.GetValue(element, point);
                    if (!std::isfinite(rho)) throw std::domain_error("Nonfinite density slice value.");
                    return rho;
                };
                const auto polynomial = d::Interpolate(value(0), value(Real{0.5}), value(1));
                const Real tolerance = 512 * std::numeric_limits<double>::epsilon() * polynomial.Scale();
                const auto verify = [&](Real t, double actual) {
                    const Real error = std::abs(polynomial.Value(t) - actual);
                    local.maximumReconstructionError = std::max(local.maximumReconstructionError, static_cast<double>(error));
                    local.maximumReconstructionTolerance = std::max(local.maximumReconstructionTolerance, static_cast<double>(tolerance));
                    if (!std::isfinite(error) || error > tolerance)
                        throw std::runtime_error("Frozen density slice fails independent Q2 reconstruction.");
                };
                for (Real t : {Real{0.125}, Real{0.25}, Real{0.375}, Real{0.625}, Real{0.75}, Real{0.875}})
                    verify(t, value(t));
                const auto parts = d::Partition(polynomial);
                ++local.lines;
                const auto evaluate = [&](Real t, d::Measure& measure, bool negativeRegion) {
                    const auto point = d::Point(axis, u, v, t);
                    if (mapper.EvaluatePoint(*transformation, point, mapped) != mean_field::mapping::MappingStatus::valid || mapped.compactified)
                        throw std::domain_error("Invalid stellar density-region mapping.");
                    const double rho = state.density.GetValue(element, point), h = state.enthalpy.GetValue(element, point);
                    if (!std::isfinite(rho) || !std::isfinite(h)) throw std::domain_error("Nonfinite density-region witness.");
                    verify(t, rho);
                    if (rho < measure.minimumRho) {
                        measure.minimumRho = rho; measure.witnessH = h; measure.witnessElement = element;
                        measure.witnessT = static_cast<double>(t); measure.witnessU = u; measure.witnessV = v;
                        for (int component = 0; component < 3; ++component) measure.witnessPosition[component] = mapped.physical_position(component);
                    }
                    if (negativeRegion) {
                        const double radius = mapped.physical_position.Norml2();
                        if (!std::isfinite(radius)) throw std::domain_error("Nonfinite density-region radius.");
                        measure.negativeRadiusMin = std::min(measure.negativeRadiusMin, radius);
                        measure.negativeRadiusMax = std::max(measure.negativeRadiusMax, radius);
                    }
                    transformation->SetIntPoint(&point);
                    const double referenceDeterminant = transformation->Jacobian().Det();
                    const Real determinant = static_cast<Real>(referenceDeterminant) * mapped.mapping_determinant;
                    if (!(referenceDeterminant > 0) || !std::isfinite(referenceDeterminant) ||
                        !(mapped.mapping_determinant > 0) || !std::isfinite(mapped.mapping_determinant) ||
                        !(determinant > 0) || !std::isfinite(determinant))
                        throw std::domain_error("Nonpositive or nonfinite density-region physical measure.");
                    return determinant;
                };
                for (const auto& part : parts) {
                    d::Measure piece;
                    std::array<double, 3> leftPosition{}, rightPosition{};
                    (void)evaluate(part.left, piece, part.negative);
                    for (int component = 0; component < 3; ++component) leftPosition[component] = mapped.physical_position(component);
                    (void)evaluate(part.right, piece, part.negative);
                    for (int component = 0; component < 3; ++component) rightPosition[component] = mapped.physical_position(component);
                    if (polynomial.a != 0) {
                        const Real stationary = -polynomial.b / (2 * polynomial.a);
                        if (stationary > part.left && stationary < part.right) (void)evaluate(stationary, piece, part.negative);
                    }
                    for (int it = 0; it < intervalRule.GetNPoints(); ++it) {
                        const auto& qp = intervalRule.IntPoint(it);
                        const Real t = part.left + (part.right - part.left) * qp.x;
                        const Real determinant = evaluate(t, piece, part.negative);
                        const Real weight = tangentWeight * (part.right - part.left) * qp.weight * determinant;
                        const Real rho = polynomial.Value(t);
                        if (!(weight > 0) || !std::isfinite(weight) || !std::isfinite(rho) ||
                            (part.negative && rho > 0) || (!part.negative && rho < 0))
                            throw std::runtime_error("Density region quadrature disagrees with its root partition.");
                        piece.volume += weight;
                        piece.mass += weight * rho;
                        ++piece.points;
                    }
                    if (part.negative) {
                        piece.negativeIntervals = 1;
                        piece.negativeVolume = piece.volume;
                        piece.negativeMass = -piece.mass;
                        piece.negativeTMin = part.left;
                        piece.negativeTMax = part.right;
                        if (tangentialOrder == 30) {
                            intervals << tangentialOrder << ',' << element << ',' << fem.mesh->GetAttribute(element) << ','
                                << axis.coordinate << ',' << axis.orientation << ',' << axis.face << ',' << u << ',' << v << ','
                                << tangentWeight << ',' << part.left << ',' << part.right << ',' << polynomial.c << ',' << polynomial.b << ','
                                << polynomial.a << ',' << piece.negativeVolume << ',' << piece.negativeMass << ',' << piece.minimumRho << ','
                                << piece.witnessH << ',' << piece.negativeRadiusMin << ',' << piece.negativeRadiusMax;
                            for (double x : leftPosition) intervals << ',' << x;
                            for (double x : rightPosition) intervals << ',' << x;
                            intervals << '\n';
                        }
                    }
                    local.Add(piece);
                }
            }
            if (!(local.volume > 0) || !std::isfinite(local.volume) || !std::isfinite(local.mass) ||
                !std::isfinite(local.negativeVolume) || !std::isfinite(local.negativeMass) || local.negativeMass < 0)
                throw std::runtime_error("Invalid integrated density-region measure.");
            elements << tangentialOrder << ',' << intervalOrder << ',' << element << ',' << fem.mesh->GetAttribute(element) << ','
                << axis.coordinate << ',' << axis.orientation << ',' << axis.face << ',' << axis.faceCount << ',';
            d::Write(elements, local);
            global.Add(local);
        }
        if (!(global.volume > 0)) throw std::invalid_argument("Density region diagnostic found no stellar volume.");
        summary << tangentialOrder << ',' << intervalOrder << ',';
        d::Write(summary, global);
        summary.flush(); elements.flush(); intervals.flush();
    }
    const bool unchanged = d::Identical(savedDensity, state.density) && d::Identical(savedEnthalpy, state.enthalpy) &&
        d::Identical(savedPotential, state.potential) && d::Identical(savedGravity, state.gravityGradientReference) &&
        d::Identical(savedDisplacement, state.displacement);
    if (!unchanged) throw std::logic_error("Density region observation changed a source field.");
    auto scope = d::Open(output / "density_negative_regions_scope.txt");
    scope << "source_fields_bitwise_identical=1\ntrial_preparations=0\nmpi_ranks=1\n"
          << "n=" << n << "\nK=" << K << "\ndensity_basis=" << fem.densityFes->FEColl()->Name() << '\n'
          << "Mass uses the reconstructed Q2 polynomial, verified against direct FE values; no density or EOS clipping.\n"
          << "Volume weights contain one reference determinant and one deformation determinant.\n"
          << "All stellar cells are included; surface axes point inward from the lowest-index stellar/nonstellar face.\n"
          << "Other cells use reference axis0. Tangential orders14/22/30 use8/12/16 points per axis; interval order14 uses8.\n"
          << "Whole-cell denominators integrate every root-partition interval with the same sliced quadrature.\n"
          << "Negative intervals are exact for each reconstructed Q2 line up to floating arithmetic; tangential integration is sampled.\n"
          << "Witnesses and radial bounds are samples, not global extrema or polynomial sign certificates.\n"
          << "Near-double-root regions below reconstruction uncertainty are not certified; compare quadrature levels and errors.\n";
}

} // namespace experiment::plateau_diagnostics