#pragma once

// Include after import mean_field. These observations use the already prepared
// base and private grid functions; they do not prepare a trial or alter a solver
// buffer. Scalar Eulerian variations are local, broken derivatives for DG
// fields. They do not contain jump distributions and cannot establish a gauge.
#include <algorithm>
#include <array>
#include <cmath>
#include <cstring>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <limits>
#include <memory>
#include <stdexcept>
#include <string>
#include <string_view>
#include <utility>

#include <mfem.hpp>
#include <mpi.h>

namespace experiment::plateau_diagnostics {
namespace shape_field_physical_detail {
using DomainSchema = mean_field::utils::domain::CoreEnvelopeVacuumDomainSchema;

// DomainMapper caches reference geometry even through its const interface.
// A private instance prevents the observation from warming the solver's cache.
// Its exterior policy deliberately has no numerical behavior: all observations
// below must be interior under both this map and the prepared production map.
class RejectExteriorMap final : public mean_field::mapping::compactification::ExteriorDomainMap {
public:
    mean_field::mapping::MappingStatus Evaluate(
        const mean_field::mapping::compactification::ExteriorMapInput &,
        mean_field::mapping::compactification::ExteriorMapResult &) const override {
        throw std::logic_error("Shape-field physical diagnostic cannot evaluate the exterior domain.");
    }
    mean_field::mapping::MappingStatus EvaluateVariation(
        const mean_field::mapping::compactification::ExteriorMapInput &,
        const mean_field::mapping::compactification::ExteriorMapResult &,
        const mean_field::mapping::compactification::ExteriorMapDirection &,
        mean_field::mapping::compactification::ExteriorMapVariation &) const override {
        throw std::logic_error("Shape-field physical diagnostic cannot vary the exterior domain.");
    }
    std::string_view GetName() const noexcept override { return "diagnostic_interior_only"; }
};

inline void Require(bool condition, const std::string &message) {
    if (!condition) throw std::runtime_error("Shape-field physical diagnostic: " + message);
}

inline bool Finite(const mfem::Vector &vector) {
    for (int i = 0; i < vector.Size(); ++i)
        if (!std::isfinite(vector(i))) return false;
    return true;
}

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

inline std::string Quote(std::string_view text) {
    std::string result = "\"";
    for (const char ch : text) {
        if (ch == '"') result += '"';
        result += ch;
    }
    return result + '"';
}

inline std::ofstream Output(const std::filesystem::path &directory, std::string_view suffix) {
    const auto path = directory / ("weak_shape_physical_" + std::string(suffix) + ".csv");
    Require(!std::filesystem::exists(path), "output already exists: " + path.string());
    std::ofstream output(path);
    output.exceptions(std::ios::failbit | std::ios::badbit);
    output << std::setprecision(std::numeric_limits<long double>::max_digits10);
    return output;
}

struct Block { int offset; int size; };

template <class Blocks>
Block RequireBlock(const Blocks &blocks, std::string_view name, int size) {
    for (const auto &block : blocks) {
        if (block.stableId != name) continue;
        Require(block.offset >= 0 && block.size >= 0 && block.offset <= size - block.size,
                "invalid root block " + std::string(name));
        return {block.offset, block.size};
    }
    throw std::invalid_argument("Shape-field physical diagnostic requires root block " + std::string(name));
}

template <class Field, class Blocks>
double Scatter(const Blocks &blocks, std::string_view name, const mfem::Vector &values,
               mfem::ParGridFunction &field) {
    const auto block = RequireBlock(blocks, name, values.Size());
    const auto adapter = mean_field::field::make_field_dof_grid_function_adapter<Field, DomainSchema>(
        *field.ParFESpace());
    Require(block.size == adapter.dof_map().reduced_size(), "field adapter size disagrees with " + std::string(name));
    mfem::Vector reduced(block.size);
    for (int i = 0; i < block.size; ++i) reduced(i) = values(block.offset + i);
    adapter.scatter(reduced, field);
    const auto gathered = adapter.gather(field);
    double maximumError = 0.;
    for (int i = 0; i < block.size; ++i)
        maximumError = std::max(maximumError, std::abs(gathered(i) - reduced(i)));
    Require(Finite(field) && std::isfinite(maximumError), "nonfinite field reconstruction");
    return maximumError;
}

// Every entry is a raw additive physical-volume integral. Keeping the signed
// cross term lets an independent reader verify E^2 = D^2 + T^2 - 2 D.T.
struct ScalarIntegrals {
    long double base{}, pullback{}, transport{}, eulerian{};
    long double baseSquared{}, pullbackSquared{}, transportSquared{}, eulerianSquared{}, cross{};

    void Add(double weight, double f, double delta, double advection) {
        const long double w = weight, b = f, d = delta, t = advection, e = d - t;
        base += w * b;
        pullback += w * d;
        transport += w * t;
        eulerian += w * e;
        baseSquared += w * b * b;
        pullbackSquared += w * d * d;
        transportSquared += w * t * t;
        eulerianSquared += w * e * e;
        cross += w * d * t;
    }

    ScalarIntegrals &operator+=(const ScalarIntegrals &other) {
        base += other.base; pullback += other.pullback; transport += other.transport; eulerian += other.eulerian;
        baseSquared += other.baseSquared; pullbackSquared += other.pullbackSquared;
        transportSquared += other.transportSquared; eulerianSquared += other.eulerianSquared; cross += other.cross;
        return *this;
    }

    void Write(std::ostream &output) const {
        for (const auto value : {base, pullback, transport, eulerian, baseSquared, pullbackSquared,
                                 transportSquared, eulerianSquared, cross})
            Require(std::isfinite(value), "nonfinite integrated scalar quantity");
        output << ',' << base << ',' << pullback << ',' << transport << ',' << eulerian << ',' << baseSquared
               << ',' << pullbackSquared << ',' << transportSquared << ',' << eulerianSquared << ',' << cross;
    }
};

inline constexpr std::string_view scalarColumns =
    ",base_integral,pullback_integral,transport_integral,eulerian_integral,base_squared_integral,"
    "pullback_squared_integral,transport_squared_integral,eulerian_squared_integral,pullback_transport_cross_integral";
} // namespace shape_field_physical_detail

template <class State>
void MeasureShapeFieldPhysical(State &state, const mean_field::fem::FEM &fem,
                               const mfem::Vector &physicalBase, const mfem::Vector &physicalDirection,
                               const std::filesystem::path &output) {
    namespace d = shape_field_physical_detail;
    using namespace mean_field;
    const auto &problem = state.Problem();
    int ranks = 0;
    MPI_Comm_size(problem.GetCommunicator(), &ranks);
    d::Require(ranks == 1, "exactly one MPI rank is required");
    d::Require(fem.mesh && fem.displacementFes && fem.surfaceDeformationFes && fem.densityFes &&
               fem.enthalpyFes && fem.gravityPotentialFes && fem.gravityFluxFes && fem.compactificationCoordinate,
               "incomplete finite-element spaces");
    d::Require(fem.mesh->SpaceDimension() == 3 && fem.mesh->Dimension() == 3,
               "three-dimensional volume elements are required");
    d::Require(state.normalizedOperator->IsPrepared() && problem.GetPhysicalOperator().IsPrepared(),
               "the base operator is not prepared");
    d::Require(physicalBase.Size() == problem.StateSize() && physicalDirection.Size() == physicalBase.Size() &&
               d::Finite(physicalBase) && d::Finite(physicalDirection) &&
               d::Identical(physicalBase, state.normalizedOperator->GetPhysicalState()),
               "the finite supplied base must exactly equal the prepared physical state");

    // Capture what a caller could otherwise mistake for a silent preparation.
    // The only production action below is the immutable deformation Jacobian.
    const auto generation = problem.GetPreparationGeneration();
    const mfem::Vector accepted(state.acceptedNormalizedState), acceptedResidual(state.acceptedNormalizedResidual),
        acceptedPhysical(state.AcceptedPhysicalState()), trial(state.trialNormalizedState),
        trialResidual(state.trialNormalizedResidual), candidatePhysical(state.candidatePhysicalState),
        preparedResidual(state.normalizedOperator->GetPhysicalResidual()),
        stateFactors(state.normalizedOperator->GetNormalization().StateFactors()),
        residualFactors(state.normalizedOperator->GetNormalization().ResidualFactors());
    const auto &physical = problem.GetPhysicalOperator();
    const mfem::Vector baseDisplacement(physical.GetGeneratedVolumeDisplacement());
    mfem::Vector volumeDirection(fem.displacementFes->GetTrueVSize());
    problem.BuildVolumeDisplacementDirection(physicalDirection, volumeDirection);
    d::Require(volumeDirection.Size() == fem.displacementFes->GetTrueVSize() && d::Finite(volumeDirection),
               "invalid generated volume direction");

    mfem::ParGridFunction rho(fem.densityFes.get()), drho(fem.densityFes.get()),
        h(fem.enthalpyFes.get()), dh(fem.enthalpyFes.get()),
        phi(fem.gravityPotentialFes.get()), dphi(fem.gravityPotentialFes.get()),
        displacement(fem.displacementFes.get()), velocity(fem.displacementFes.get()),
        affine(fem.displacementFes.get());
    const auto blocks = problem.GetManifest().valueBlocks();
    const std::array<double, 6> scatterErrors{
        d::Scatter<field::Density>(blocks, "density", physicalBase, rho),
        d::Scatter<field::Density>(blocks, "density", physicalDirection, drho),
        d::Scatter<field::Enthalpy>(blocks, "specific_enthalpy", physicalBase, h),
        d::Scatter<field::Enthalpy>(blocks, "specific_enthalpy", physicalDirection, dh),
        d::Scatter<field::Gravity>(blocks, "gravity_potential", physicalBase, phi),
        d::Scatter<field::Gravity>(blocks, "gravity_potential", physicalDirection, dphi)};
    displacement.SetFromTrueDofs(baseDisplacement);
    velocity.SetFromTrueDofs(volumeDirection);
    mfem::VectorFunctionCoefficient identity(3, [](const mfem::Vector &x, mfem::Vector &value) { value = x; });
    affine.ProjectCoefficient(identity);
    d::Require(d::Finite(displacement) && d::Finite(velocity) && d::Finite(affine), "nonfinite private displacement field");
    const auto &productionMapper = problem.GetDiscretization().domainMapper();
    const auto productionCacheBefore = productionMapper.GetGeometryCacheStatistics();
    mapping::DomainMapper privateMapper(
        utils::DomainMapperOptions{.dimension = 3,
            .vacuum_element_attribute = d::DomainSchema::template material_attribute<utils::domain::Vacuum>()},
        std::make_unique<d::RejectExteriorMap>());
    privateMapper.SetGeometryCacheEnabled(false);
    mapping::GridFunctionMappingEvaluator mappingEvaluator(privateMapper, displacement, *fem.compactificationCoordinate);
    mappingEvaluator.InvalidateCache();

    // The default surface descriptor erases its center. Checking both actual
    // displacement traces prevents a different radial prescription from being
    // silently interpreted in the zero-center basis used for this CSV.
    const auto &domain = physical.GetDomainDeformation();
    mfem::Vector center(3); center = 0.;
    const deformation::NodalRadialSurface prescription(std::move(center));
    d::Require(domain.descriptor().surfaceDeformation == prescription.descriptor() &&
               fem.surfaceDeformationFes->GetVDim() == 1 && fem.displacementFes->GetVDim() == 3 &&
               fem.displacementFes->GetOrdering() == mfem::Ordering::byNODES &&
               fem.displacementFes->FEColl() == fem.surfaceDeformationFes->FEColl() &&
               fem.displacementFes->GetTrueVSize() == 3 * fem.surfaceDeformationFes->GetTrueVSize(),
               "the default nodal radial surface and its vector companion are required");
    const auto radial = deformation::compileSurfaceDeformationPrescription(prescription,
        deformation::SurfaceDeformationCompilationContext{
            *fem.surfaceDeformationFes, field::make_stellar_surface_scalar_dof_map(*fem.surfaceDeformationFes)});
    const auto shape = d::RequireBlock(blocks, "surface_deformation", physicalBase.Size());
    d::Require(radial.parameterCount() == shape.size && shape.size == domain.parameterCount(),
               "surface parameter counts disagree");
    const int scalarSize = fem.surfaceDeformationFes->GetTrueVSize();
    double traceError = 0., traceScale = 1.;
    for (int i = 0; i < shape.size; ++i)
        for (int component = 0; component < 3; ++component) {
            const int dof = radial.surfaceDofMap().volume_true_dof(i) + component * scalarSize;
            const double direction = radial.radialDirection(i, component);
            const double expectedBase = physicalBase(shape.offset + i) * direction;
            const double expectedDirection = physicalDirection(shape.offset + i) * direction;
            traceError = std::max({traceError, std::abs(baseDisplacement(dof) - expectedBase),
                                  std::abs(volumeDirection(dof) - expectedDirection)});
            traceScale = std::max({traceScale, std::abs(expectedBase), std::abs(expectedDirection)});
        }
    d::Require(std::isfinite(traceError) && traceError <= 512. * std::numeric_limits<double>::epsilon() * traceScale,
               "actual generated surface trace disagrees with the reconstructed radial basis");

    std::filesystem::create_directories(output);
    auto metadata = d::Output(output, "metadata");
    auto elements = d::Output(output, "elements");
    auto geometry = d::Output(output, "geometry");
    auto affineOutput = d::Output(output, "affine");
    auto surface = d::Output(output, "surface_nodes");
    auto summary = d::Output(output, "summary");
    metadata << "key,value\n";
    const auto meta = [&](std::string_view key, std::string_view value) {
        metadata << d::Quote(key) << ',' << d::Quote(value) << '\n';
    };
    meta("schema", "weak_shape_physical_v1");
    meta("scope", "stellar interior only; prepared base; physical direction supplied; no trial preparation");
    meta("scalar_formula", "delta_f_E = delta_f_at_fixed_reference - grad_physical_f dot volume_displacement_direction");
    meta("dg_derivatives", "density and potential gradients are element-local broken derivatives; jumps excluded");
    meta("interpretation", "local shape-field compensation diagnostic; no nullspace or pure-gauge certification");
    meta("quadrature", "fixed integration-rule orders 10 and 14, independent of field order; compare for integration sensitivity");
    meta("radius", "Euclidean distance from zero; per-element min/max are quadrature-sample extrema");
    meta("surface_samples", "nodal values without area weights; not harmonic coefficients or surface integrals");
    meta("affine_test", "projection of reference physical coordinate x into displacement space; independent of the actual domain lift");
    meta("affine_weights", "reference physical mesh volume; gradient is with respect to reference physical coordinates");
    meta("physical_weights", "integration_point.weight * element_transformation.Weight * mapping_determinant");
    meta("mapping", "private uncached DomainMapper; same interior x+u implementation; exterior evaluation forbidden");
    meta("scalar_integrals", "raw additive integrals; transport = grad_physical_f dot v; no distributional interface contribution");
    meta("gravity_vector", "not evaluated; RT Piola variation is outside this scalar diagnostic");
    meta("reference_mesh_geometry_collection", fem.mesh->GetNodes()
        ? fem.mesh->GetNodes()->FESpace()->FEColl()->Name() : "linear_mesh_without_nodal_grid_function");
    const std::array<std::pair<std::string_view, const mfem::ParFiniteElementSpace *>, 6> spaces{{
        {"density", fem.densityFes.get()}, {"specific_enthalpy", fem.enthalpyFes.get()},
        {"gravity_potential", fem.gravityPotentialFes.get()}, {"displacement", fem.displacementFes.get()},
        {"surface_deformation", fem.surfaceDeformationFes.get()}, {"gravity_gradient", fem.gravityFluxFes.get()}}};
    for (const auto &[name, space] : spaces) {
        meta(std::string(name) + "_collection", space->FEColl()->Name());
        meta(std::string(name) + "_true_size", std::to_string(space->GetTrueVSize()));
    }
    meta("density_family_order", std::to_string(field::Density::Scalar::familyOrder));
    meta("enthalpy_family_order", std::to_string(field::Enthalpy::Scalar::familyOrder));
    meta("potential_family_order", std::to_string(field::Gravity::Potential::familyOrder));
    meta("gravity_flux_family_order", std::to_string(field::Gravity::Flux::familyOrder));
    meta("displacement_family_order", std::to_string(field::Displacement::Vector::familyOrder));

    elements << "quadrature,element,attribute,geometry,field,gradient_kind,points,physical_volume"
             << d::scalarColumns << '\n';
    geometry << "quadrature,element,attribute,geometry,points,reference_mesh_geometry_order,density_element_order,enthalpy_element_order,"
                "potential_element_order,displacement_element_order,reference_radius_min,reference_radius_max,"
                "physical_radius_min,physical_radius_max,physical_volume,base_displacement_squared_integral,"
                "direction_displacement_squared_integral\n";
    affineOutput << "quadrature,element,points,reference_volume,value_error_squared_integral,"
                    "reference_value_squared_integral,gradient_error_squared_integral,reference_gradient_squared_integral\n";
    surface << "parameter,volume_scalar_true_dof,unit_radial_x,unit_radial_y,unit_radial_z,reference_radius,"
               "physical_radius,base_parameter,physical_direction_parameter\n";
    for (int i = 0; i < shape.size; ++i) {
        const double radius = radial.referenceRadius(i), base = physicalBase(shape.offset + i);
        d::Require(std::isfinite(radius) && radius > 0. && radius + base > 0., "invalid radial surface radius");
        surface << i << ',' << radial.surfaceDofMap().volume_true_dof(i);
        for (int component = 0; component < 3; ++component) surface << ',' << radial.radialDirection(i, component);
        surface << ',' << radius << ',' << radius + base << ',' << base << ',' << physicalDirection(shape.offset + i) << '\n';
    }
    summary << "quadrature,field,metric,value\n";
    const auto measure = [&](int q, std::string_view fieldName, std::string_view metric, long double value) {
        d::Require(std::isfinite(value), "nonfinite summary " + std::string(metric));
        summary << q << ',' << fieldName << ',' << metric << ',' << value << '\n';
    };
    const std::array<std::string_view, 3> names{"density", "specific_enthalpy", "gravity_potential"};
    const std::array<mfem::ParGridFunction *, 3> fields{&rho, &h, &phi}, directions{&drho, &dh, &dphi};
    mfem::Vector referenceGradient(3), physicalGradient(3), u(3), v(3), x(3), affineValue(3);
    mfem::DenseMatrix affineGradient(3);
    mapping::MappingPointContext mapped;
    for (const int quadrature : {10, 14}) {
        std::array<d::ScalarIntegrals, 3> totals;
        long double totalVolume = 0., totalU2 = 0., totalV2 = 0., affineVolume = 0., affineError2 = 0.,
            affineReference2 = 0., affineGradientError2 = 0.;
        int stellarElements = 0;
        for (int element = 0; element < fem.mesh->GetNE(); ++element) {
            const int attribute = fem.mesh->GetAttribute(element);
            if (!d::DomainSchema::template attribute_belongs_to<utils::domain::Stellar>(attribute)) continue;
            ++stellarElements;
            auto *transformation = fem.mesh->GetElementTransformation(element);
            d::Require(!productionMapper.IsCompactifiedElement(*transformation),
                       "production mapper classifies a stellar element as compactified");
            const auto &rule = mfem::IntRules.Get(transformation->GetGeometryType(), quadrature);
            std::array<d::ScalarIntegrals, 3> integrals;
            long double volume = 0., u2 = 0., v2 = 0., referenceVolume = 0., affineE2 = 0., affineR2 = 0., affineGE2 = 0.;
            double rmin = std::numeric_limits<double>::infinity(), rmax = 0.,
                physicalRmin = std::numeric_limits<double>::infinity(), physicalRmax = 0.;
            for (int point = 0; point < rule.GetNPoints(); ++point) {
                const auto &ip = rule.IntPoint(point);
                d::Require(mappingEvaluator.EvaluatePoint(*transformation, ip, mapped) == mapping::MappingStatus::valid &&
                           !mapped.compactified && std::isfinite(mapped.mapping_determinant) && mapped.mapping_determinant > 0.,
                           "invalid or compactified stellar mapping");
                transformation->SetIntPoint(&ip);
                const double referenceWeight = ip.weight * transformation->Weight();
                const double weight = referenceWeight * mapped.mapping_determinant;
                d::Require(std::isfinite(referenceWeight) && referenceWeight > 0. && std::isfinite(weight) && weight > 0.,
                           "invalid quadrature weight");
                displacement.GetVectorValue(element, ip, u);
                velocity.GetVectorValue(element, ip, v);
                d::Require(d::Finite(u) && d::Finite(v) && d::Finite(mapped.reference_position) && d::Finite(mapped.physical_position),
                           "nonfinite sampled displacement or position");
                const double r = mapped.reference_position.Norml2(), physicalR = mapped.physical_position.Norml2();
                d::Require(std::isfinite(r) && std::isfinite(physicalR), "overflowing sampled radius");
                rmin = std::min(rmin, r); rmax = std::max(rmax, r);
                physicalRmin = std::min(physicalRmin, physicalR); physicalRmax = std::max(physicalRmax, physicalR);
                volume += weight;
                for (int component = 0; component < 3; ++component) {
                    u2 += static_cast<long double>(weight) * u(component) * u(component);
                    v2 += static_cast<long double>(weight) * v(component) * v(component);
                }
                for (std::size_t fieldIndex = 0; fieldIndex < fields.size(); ++fieldIndex) {
                    const double base = fields[fieldIndex]->GetValue(element, ip);
                    const double delta = directions[fieldIndex]->GetValue(element, ip);
                    transformation->SetIntPoint(&ip);
                    fields[fieldIndex]->GetGradient(*transformation, referenceGradient);
                    mapped.inverse_mapping_jacobian.MultTranspose(referenceGradient, physicalGradient);
                    const double transport = physicalGradient * v;
                    d::Require(std::isfinite(base) && std::isfinite(delta) && std::isfinite(transport) && d::Finite(physicalGradient),
                               "nonfinite scalar or broken gradient");
                    integrals[fieldIndex].Add(weight, base, delta, transport);
                }
                transformation->SetIntPoint(&ip);
                transformation->Transform(ip, x);
                affine.GetVectorValue(element, ip, affineValue);
                affine.GetVectorGradient(*transformation, affineGradient);
                d::Require(d::Finite(x) && d::Finite(affineValue), "nonfinite affine projection sample");
                referenceVolume += referenceWeight;
                for (int row = 0; row < 3; ++row) {
                    const long double error = static_cast<long double>(affineValue(row)) - x(row);
                    affineE2 += referenceWeight * error * error;
                    affineR2 += static_cast<long double>(referenceWeight) * x(row) * x(row);
                    for (int column = 0; column < 3; ++column) {
                        const long double gradientError = affineGradient(row, column) - (row == column ? 1. : 0.);
                        d::Require(std::isfinite(gradientError), "nonfinite affine projection gradient");
                        affineGE2 += referenceWeight * gradientError * gradientError;
                    }
                }
            }
            for (std::size_t fieldIndex = 0; fieldIndex < fields.size(); ++fieldIndex) {
                elements << quadrature << ',' << element << ',' << attribute << ',' << int(transformation->GetGeometryType())
                         << ',' << names[fieldIndex] << ',' << (fieldIndex == 1 ? "continuous_H1" : "broken_DG")
                         << ',' << rule.GetNPoints() << ',' << volume;
                integrals[fieldIndex].Write(elements);
                elements << '\n';
                totals[fieldIndex] += integrals[fieldIndex];
            }
            geometry << quadrature << ',' << element << ',' << attribute << ',' << int(transformation->GetGeometryType())
                     << ',' << rule.GetNPoints() << ',' << (fem.mesh->GetNodes()
                         ? fem.mesh->GetNodes()->FESpace()->GetElementOrder(element) : 1)
                     << ',' << fem.densityFes->GetElementOrder(element)
                     << ',' << fem.enthalpyFes->GetElementOrder(element) << ',' << fem.gravityPotentialFes->GetElementOrder(element)
                     << ',' << fem.displacementFes->GetElementOrder(element) << ',' << rmin << ',' << rmax << ',' << physicalRmin
                     << ',' << physicalRmax << ',' << volume << ',' << u2 << ',' << v2 << '\n';
            affineOutput << quadrature << ',' << element << ',' << rule.GetNPoints() << ',' << referenceVolume << ','
                         << affineE2 << ',' << affineR2 << ',' << affineGE2 << ',' << 3. * referenceVolume << '\n';
            totalVolume += volume; totalU2 += u2; totalV2 += v2;
            affineVolume += referenceVolume; affineError2 += affineE2; affineReference2 += affineR2; affineGradientError2 += affineGE2;
        }
        d::Require(stellarElements > 0 && totalVolume > 0. && affineVolume > 0. && affineReference2 > 0., "empty stellar integrals");
        for (std::size_t i = 0; i < totals.size(); ++i) {
            const auto &integral = totals[i];
            measure(quadrature, names[i], "base_squared_integral", integral.baseSquared);
            measure(quadrature, names[i], "pullback_squared_integral", integral.pullbackSquared);
            measure(quadrature, names[i], "transport_squared_integral", integral.transportSquared);
            measure(quadrature, names[i], "eulerian_squared_integral", integral.eulerianSquared);
            measure(quadrature, names[i], "pullback_transport_cross_integral", integral.cross);
            measure(quadrature, names[i], "pullback_l2_norm", std::sqrt(integral.pullbackSquared));
            measure(quadrature, names[i], "transport_l2_norm", std::sqrt(integral.transportSquared));
            measure(quadrature, names[i], "eulerian_l2_norm", std::sqrt(integral.eulerianSquared));
            measure(quadrature, names[i], "squared_norm_identity_error",
                    integral.eulerianSquared - integral.pullbackSquared - integral.transportSquared + 2. * integral.cross);
        }
        measure(quadrature, "geometry", "stellar_elements", stellarElements);
        measure(quadrature, "geometry", "physical_volume", totalVolume);
        measure(quadrature, "geometry", "base_displacement_squared_integral", totalU2);
        measure(quadrature, "geometry", "direction_displacement_squared_integral", totalV2);
        measure(quadrature, "affine", "reference_volume", affineVolume);
        measure(quadrature, "affine", "value_error_squared_integral", affineError2);
        measure(quadrature, "affine", "reference_value_squared_integral", affineReference2);
        measure(quadrature, "affine", "gradient_error_squared_integral", affineGradientError2);
        measure(quadrature, "affine", "reference_gradient_squared_integral", 3. * affineVolume);
        measure(quadrature, "affine", "relative_value_l2_error", std::sqrt(affineError2 / affineReference2));
        measure(quadrature, "affine", "relative_gradient_l2_error", std::sqrt(affineGradientError2 / (3. * affineVolume)));
    }
    measure(0, "reconstruction", "maximum_scatter_roundtrip_error", *std::max_element(scatterErrors.begin(), scatterErrors.end()));
    measure(0, "reconstruction", "maximum_surface_trace_error", traceError);
    const auto productionCacheAfter = productionMapper.GetGeometryCacheStatistics();
    const bool cacheUnchanged = productionCacheBefore.hits == productionCacheAfter.hits &&
        productionCacheBefore.misses == productionCacheAfter.misses && productionCacheBefore.points == productionCacheAfter.points;
    const bool unchanged = cacheUnchanged &&
        generation == problem.GetPreparationGeneration() && state.normalizedOperator->IsPrepared() &&
        d::Identical(physicalBase, state.normalizedOperator->GetPhysicalState()) &&
        d::Identical(preparedResidual, state.normalizedOperator->GetPhysicalResidual()) &&
        d::Identical(accepted, state.acceptedNormalizedState) && d::Identical(acceptedResidual, state.acceptedNormalizedResidual) &&
        d::Identical(acceptedPhysical, state.AcceptedPhysicalState()) && d::Identical(trial, state.trialNormalizedState) &&
        d::Identical(trialResidual, state.trialNormalizedResidual) && d::Identical(candidatePhysical, state.candidatePhysicalState) &&
        d::Identical(baseDisplacement, physical.GetGeneratedVolumeDisplacement()) &&
        d::Identical(stateFactors, state.normalizedOperator->GetNormalization().StateFactors()) &&
        d::Identical(residualFactors, state.normalizedOperator->GetNormalization().ResidualFactors());
    measure(0, "invariant", "accepted_trial_prepared_normalization_unchanged", unchanged);
    measure(0, "invariant", "production_mapping_cache_unchanged", cacheUnchanged);
    measure(0, "invariant", "trial_preparations", 0);
    measure(0, "invariant", "volume_direction_actions", 1);
    d::Require(unchanged, "solver state or normalization changed during physical observation");
}
} // namespace experiment::plateau_diagnostics