#pragma once

// Frozen-field order controls, not solver recipes. All target spaces borrow the
// same mesh and retain the source enthalpy and deformation without modification.
#include "negative_density_projection.hpp"
#include "density_polynomial_regions.hpp"

#include <iostream>
#include <chrono>
#include <map>
#include <span>

namespace experiment::plateau_diagnostics {
namespace density_order_detail {
namespace projection = negative_density_projection_detail;
using ObserverClock = std::chrono::steady_clock;

inline double Seconds(ObserverClock::time_point start, ObserverClock::time_point finish) {
    return std::chrono::duration<double>(finish - start).count();
}

inline void WriteTimings(const std::filesystem::path& path, ObserverClock::time_point start,
                         ObserverClock::time_point prepared, ObserverClock::time_point errors,
                         ObserverClock::time_point finished) {
    auto out = region::Open(path);
    out << "phase,seconds\nsetup_and_projection," << Seconds(start, prepared)
        << "\nindependent_error_quadrature," << Seconds(prepared, errors)
        << "\nregion_observation_and_output," << Seconds(errors, finished)
        << "\nobserver_total," << Seconds(start, finished) << '\n';
}

inline int ValidateState(polytrope_validation::PhysicalState& state, double n, double K) {
    const auto& fem = state.finiteElements;
    if (!std::isfinite(n) || n < 1 || !std::isfinite(K) || K <= 0 || !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 order observations require the current conforming 3D scalar L2 state and finite n>=1,K>0.");
    int ranks = 0;
    MPI_Comm_size(fem.mesh->GetComm(), &ranks);
    if (ranks != 1) throw std::invalid_argument("Density order observations require one MPI rank.");
    int degree = -1;
    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->GetOrder() > 4 || fe->GetMapType() != mfem::FiniteElement::VALUE ||
            fe->GetGeomType() != mfem::Geometry::CUBE || fe->GetDof() != (fe->GetOrder() + 1) * (fe->GetOrder() + 1) * (fe->GetOrder() + 1))
            throw std::invalid_argument("Density order observations require scalar VALUE tensor Q2/Q3/Q4 elements.");
        if (degree >= 0 && degree != fe->GetOrder()) throw std::invalid_argument("Variable density order is outside this frozen study.");
        degree = fe->GetOrder();
    }
    if (degree < 0) throw std::invalid_argument("Density order observation found no stellar elements.");
    return degree;
}

struct SourceCheckpoint final {
    mfem::Vector density, enthalpy, potential, gravity, displacement;
    explicit SourceCheckpoint(const polytrope_validation::PhysicalState& state)
        : density(state.density), enthalpy(state.enthalpy), potential(state.potential),
          gravity(state.gravityGradientReference), displacement(state.displacement) {}
    void Verify(const polytrope_validation::PhysicalState& state) const {
        if (!region::Identical(density, state.density) || !region::Identical(enthalpy, state.enthalpy) ||
            !region::Identical(potential, state.potential) || !region::Identical(gravity, state.gravityGradientReference) ||
            !region::Identical(displacement, state.displacement))
            throw std::logic_error("Density order diagnostic changed a source field.");
    }
};

inline void WriteScope(const std::filesystem::path& path, bool projections, int sourceDegree, double n, double K) {
    auto out = region::Open(path);
    out << "source_fields_bitwise_identical=1\ntrial_preparations=0\nmpi_ranks=1\nsource_density_degree=" << sourceDegree
        << "\nn=" << n << "\nK=" << K << "\nprojection_controls=" << projections << '\n'
        << "All measures use unchanged source enthalpy and physical geometry.\n"
        << "eos_nodal_interpolation is interpolation of rho_EOS(h_h), NOT original Lane-Emden initialization.\n"
        << "eos_l2_q26 uses one fixed order26 rule for both mass and target source at each density order.\n"
        << "Surface records are one-sided samples, not area fractions.\n"
        << "Density coefficients are in different FE spaces; their array differences across orders are not physical norms.\n"
        << "EOS errors and moments are weighted physical integrals; full and negative measures share root-sliced quadrature.\n"
        << "density_order_error_quadrature.csv independently compares unsliced order26/30 EOS error; it should be used for error-norm sensitivity.\n"
        << "Complete EOS integrals are undefined if negative h is sampled; density mass and negative-region measures remain reported.\n"
        << "Projection/interpolation construction rejects negative sampled enthalpy; it does not repair the source.\n"
        << "Root-sliced measures are numerical estimates, not global positivity certificates.\n"
        << "negative_mass_uncertainty_scale is the integral of a coefficient-scaled reconstruction tolerance, not a rigorous bound.\n"
        << "Resolved-negative intervals have a direct FE witness below their reconstruction tolerance; their full volume is still an estimate.\n"
        << "Uncertain shallow intervals and surface values are retained and separately labeled; no density field is clipped.\n"
        << "Near repeated roots and narrow unsampled tangential pockets remain unresolved limitations.\n";
}

struct DegreeFields final {
    int degree;
    std::unique_ptr<mfem::L2_FECollection> collection;
    std::unique_ptr<mfem::ParFiniteElementSpace> space;
    std::unique_ptr<mfem::ParGridFunction> interpolation, projection, constant;
    mfem::DenseMatrix shapes;
    DegreeFields(mfem::ParMesh& mesh, int order, const mfem::IntegrationRule& rule)
        : degree(order), collection(std::make_unique<mfem::L2_FECollection>(order, 3, mfem::BasisType::GaussLegendre, mfem::FiniteElement::VALUE)),
          space(std::make_unique<mfem::ParFiniteElementSpace>(&mesh, collection.get())),
          interpolation(std::make_unique<mfem::ParGridFunction>(space.get())),
          projection(std::make_unique<mfem::ParGridFunction>(space.get())),
          constant(std::make_unique<mfem::ParGridFunction>(space.get())) {
        *interpolation = 0.; *projection = 0.;
        mfem::ConstantCoefficient one(1.);
        constant->ProjectCoefficient(one);
        const auto* fe = collection->FiniteElementForGeometry(mfem::Geometry::CUBE);
        shapes.SetSize(fe->GetDof(), rule.GetNPoints());
        mfem::Vector shape(fe->GetDof());
        for (int q = 0; q < rule.GetNPoints(); ++q) {
            fe->CalcShape(rule.IntPoint(q), shape);
            for (int i = 0; i < shape.Size(); ++i) shapes(i, q) = shape(i);
        }
    }
};

struct ErrorField final {
    std::string_view label;
    int degree;
    const mfem::ParGridFunction* density;
};

struct ErrorMeasure final {
    Real volume{0}, mass{0}, targetMass{0}, eosSquared{0}, signedEos{0};
    Real minimumRho{std::numeric_limits<Real>::infinity()};
    double minimumH{std::numeric_limits<double>::infinity()};
    std::size_t negativeH{0};
};

inline void ErrorHeader(std::ostream& out) {
    out << "control,degree,element,attribute,quadrature_request,quadrature_achieved,quadrature_points,volume,mass,target_mass,"
           "eos_squared_integral,signed_eos_integral,minimum_density,minimum_enthalpy,negative_enthalpy_samples,eos_defined_everywhere_sampled\n";
}

inline void WriteError(std::ostream& out, std::string_view label, int degree, int element, int attribute,
                       int order, const mfem::IntegrationRule& rule, const ErrorMeasure& measure) {
    const Real nan = std::numeric_limits<Real>::quiet_NaN();
    if (!(measure.volume > 0) || !std::isfinite(measure.volume) || !std::isfinite(measure.mass) ||
        !std::isfinite(measure.minimumRho) || !std::isfinite(measure.minimumH) ||
        (measure.negativeH == 0 && (!std::isfinite(measure.targetMass) || !std::isfinite(measure.eosSquared) || !std::isfinite(measure.signedEos))))
        throw std::runtime_error("Nonfinite independent density error measure.");
    out << label << ',' << degree << ',' << element << ',' << attribute << ',' << order << ',' << rule.GetOrder() << ','
        << rule.GetNPoints() << ',' << measure.volume << ',' << measure.mass << ','
        << (measure.negativeH == 0 ? measure.targetMass : nan) << ',' << (measure.negativeH == 0 ? measure.eosSquared : nan) << ','
        << (measure.negativeH == 0 ? measure.signedEos : nan) << ',' << measure.minimumRho << ',' << measure.minimumH << ','
        << measure.negativeH << ',' << (measure.negativeH == 0) << '\n';
}

// A squared EOS error needs more quadrature than either projection moment.
// This pass reuses each mapped sample and reference basis across all supplied
// fields; it does not assemble another matrix or solve another projection.
inline void ObserveErrorRule(polytrope_validation::PhysicalState& state, std::span<const ErrorField> fields,
                             int order, double n, double K, std::ostream& out) {
    const auto& fem = state.finiteElements;
    const auto& rule = mfem::IntRules.Get(mfem::Geometry::CUBE, order);
    std::map<const mfem::FiniteElement*, std::unique_ptr<mfem::DenseMatrix>> shapeCache;
    std::vector<const mfem::DenseMatrix*> shapes(fields.size());
    std::vector<mfem::Vector> coefficients(fields.size());
    mean_field::mapping::GridFunctionMappingEvaluator mapper(state.domainMapper, state.displacement, *fem.compactificationCoordinate);
    mean_field::mapping::MappingPointContext mapped;
    for (int element = 0; element < fem.mesh->GetNE(); ++element) {
        if (!state.isStellar(element)) continue;
        auto* transformation = fem.mesh->GetElementTransformation(element);
        for (std::size_t field = 0; field < fields.size(); ++field) {
            const auto& density = *fields[field].density;
            const auto* fe = density.ParFESpace()->GetFE(element);
            if (fe->GetOrder() != fields[field].degree || fe->GetMapType() != mfem::FiniteElement::VALUE || fe->GetGeomType() != mfem::Geometry::CUBE)
                throw std::logic_error("Independent density error quadrature received an incompatible scalar field.");
            auto& cached = shapeCache[fe];
            if (!cached) {
                cached = std::make_unique<mfem::DenseMatrix>(fe->GetDof(), rule.GetNPoints());
                mfem::Vector shape(fe->GetDof());
                for (int q = 0; q < rule.GetNPoints(); ++q) {
                    fe->CalcShape(rule.IntPoint(q), shape);
                    for (int i = 0; i < shape.Size(); ++i) (*cached)(i, q) = shape(i);
                }
            }
            shapes[field] = cached.get();
            mfem::Array<int> dofs;
            auto* dofTransform = density.ParFESpace()->GetElementDofs(element, dofs);
            density.GetSubVector(dofs, coefficients[field]);
            if (dofTransform) dofTransform->InvTransformPrimal(coefficients[field]);
        }
        std::vector<ErrorMeasure> measures(fields.size());
        for (int q = 0; q < rule.GetNPoints(); ++q) {
            const auto& point = rule.IntPoint(q);
            if (mapper.EvaluatePoint(*transformation, point, mapped) != mean_field::mapping::MappingStatus::valid || mapped.compactified)
                throw std::domain_error("Invalid independent density error mapping.");
            const double h = state.enthalpy.GetValue(element, point);
            const Real target = EosTarget(h, n, K);
            transformation->SetIntPoint(&point);
            const double reference = transformation->Jacobian().Det();
            const Real weight = static_cast<Real>(point.weight) * reference * mapped.mapping_determinant;
            if (!(reference > 0) || !std::isfinite(reference) || !(mapped.mapping_determinant > 0) ||
                !std::isfinite(mapped.mapping_determinant) || !(weight > 0) || !std::isfinite(weight))
                throw std::domain_error("Invalid independent density error weight.");
            for (std::size_t field = 0; field < fields.size(); ++field) {
                auto& measure = measures[field];
                Real rho = 0;
                const auto* shape = shapes[field]->GetColumn(q);
                for (int i = 0; i < coefficients[field].Size(); ++i) rho += static_cast<Real>(shape[i]) * coefficients[field](i);
                if (!std::isfinite(rho)) throw std::domain_error("Nonfinite independently observed density.");
                measure.volume += weight; measure.mass += weight * rho; measure.targetMass += weight * target;
                measure.eosSquared += weight * (rho - target) * (rho - target); measure.signedEos += weight * (rho - target);
                measure.minimumRho = std::min(measure.minimumRho, rho); measure.minimumH = std::min(measure.minimumH, h);
                measure.negativeH += h < 0;
            }
        }
        for (std::size_t field = 0; field < fields.size(); ++field)
            WriteError(out, fields[field].label, fields[field].degree, element, fem.mesh->GetAttribute(element), order, rule, measures[field]);
    }
    out.flush();
}
} // namespace density_order_detail

inline void RecordDensityOrderProjection(polytrope_validation::PhysicalState& state, double n, double K,
                                          const std::filesystem::path& output) {
    namespace d = density_order_detail;
    const auto observerStart = d::ObserverClock::now();
    const int sourceDegree = d::ValidateState(state, n, K);
    std::filesystem::create_directories(output);
    d::RootSelfChecks(output / "density_order_root_selfchecks.csv");
    d::SourceCheckpoint checkpoint(state);
    const auto& fem = state.finiteElements;
    const auto& rule = mfem::IntRules.Get(mfem::Geometry::CUBE, 26);
    std::vector<std::unique_ptr<d::DegreeFields>> controls;
    for (int degree : {2, 3, 4}) controls.push_back(std::make_unique<d::DegreeFields>(*fem.mesh, degree, rule));
    auto elements = d::region::Open(output / "density_order_projection_elements.csv");
    elements << "degree,element,attribute,quadrature_request,quadrature_achieved,quadrature_points,dofs,basis_family,volume,target_mass,"
                "interpolation_mass,projection_mass,interpolation_eos_squared,projection_eos_squared,rhs_norm,normal_defect_norm,"
                "normal_relative_defect,constant_coefficient_error,constant_sample_error,interpolation_coefficient_minimum,"
                "projection_coefficient_minimum,projection_coefficient_maximum\n";
    auto nodes = d::region::Open(output / "density_order_nodes.csv");
    nodes << "degree,element,local_dof,global_dof,xi,eta,zeta,h,target_rho,interpolation_coefficient,projection_coefficient\n";
    auto errorQuadrature = d::region::Open(output / "density_order_error_quadrature.csv");
    d::ErrorHeader(errorQuadrature);
    mean_field::mapping::GridFunctionMappingEvaluator mapper(state.domainMapper, state.displacement, *fem.compactificationCoordinate);
    mean_field::mapping::MappingPointContext mapped;
    mfem::Vector weights(rule.GetNPoints()), targets(rule.GetNPoints());
    for (int element = 0; element < fem.mesh->GetNE(); ++element) {
        if (!state.isStellar(element)) continue;
        auto* transformation = fem.mesh->GetElementTransformation(element);
        // Geometry and EOS evaluation dominate observer cost; share these
        // samples across target orders without sharing their independent M,b.
        long double volume = 0, targetMass = 0;
        double minimumH = std::numeric_limits<double>::infinity();
        for (int q = 0; q < rule.GetNPoints(); ++q) {
            const auto& point = rule.IntPoint(q);
            if (mapper.EvaluatePoint(*transformation, point, mapped) != mean_field::mapping::MappingStatus::valid || mapped.compactified)
                throw std::domain_error("Invalid density order projection mapping.");
            const double h = state.enthalpy.GetValue(element, point);
            if (!std::isfinite(h) || h < 0) throw std::domain_error("Density order projection encountered inadmissible sampled enthalpy.");
            targets(q) = d::EosTarget(h, n, K);
            minimumH = std::min(minimumH, h);
            transformation->SetIntPoint(&point);
            const double reference = transformation->Jacobian().Det();
            weights(q) = point.weight * reference * mapped.mapping_determinant;
            if (!(reference > 0) || !std::isfinite(reference) || !(mapped.mapping_determinant > 0) ||
                !std::isfinite(mapped.mapping_determinant) || !(weights(q) > 0) || !std::isfinite(weights(q)))
                throw std::domain_error("Invalid density order projection weight.");
            volume += weights(q); targetMass += static_cast<long double>(weights(q)) * targets(q);
        }
        for (auto& control : controls) {
            const auto* fe = control->space->GetFE(element);
            const int size = fe->GetDof();
            mfem::Array<int> dofs;
            if (control->space->GetElementDofs(element, dofs) != nullptr || dofs.Size() != size)
                throw std::invalid_argument("Density order projection requires untransformed scalar DG element DOFs.");
            for (int i = 0; i < size; ++i) if (dofs[i] < 0) throw std::logic_error("Unexpected signed DG density DOF.");
            mfem::Vector interpolated(size), nodalH(size), rhs(size), constantRhs(size), constantExpected;
            const auto& nodalRule = fe->GetNodes();
            if (nodalRule.GetNPoints() != size) throw std::logic_error("Density target collection is not nodal.");
            for (int i = 0; i < size; ++i) {
                nodalH(i) = state.enthalpy.GetValue(element, nodalRule.IntPoint(i));
                if (!std::isfinite(nodalH(i)) || nodalH(i) < 0)
                    throw std::domain_error("Density order nodal interpolation encountered inadmissible enthalpy.");
                interpolated(i) = d::EosTarget(nodalH(i), n, K);
            }
            mfem::DenseMatrix mass(size);
            mass = 0.; rhs = 0.; constantRhs = 0.;
            for (int q = 0; q < rule.GetNPoints(); ++q) {
                const auto* shape = control->shapes.GetColumn(q);
                for (int j = 0; j < size; ++j) {
                    const double weighted = weights(q) * shape[j];
                    rhs(j) += weighted * targets(q); constantRhs(j) += weighted;
                    auto* column = mass.GetColumn(j);
                    for (int i = 0; i <= j; ++i) column[i] += weighted * shape[i];
                }
            }
            for (int j = 0; j < size; ++j) for (int i = 0; i < j; ++i) mass(j, i) = mass(i, j);
            const auto lower = d::projection::Factor(mass);
            const auto projected = d::projection::Solve(lower, rhs), constant = d::projection::Solve(lower, constantRhs);
            const auto defect = d::projection::Defect(mass, projected, rhs);
            d::projection::VerifySolve(mass, projected, rhs, defect);
            d::projection::VerifySolve(mass, constant, constantRhs, d::projection::Defect(mass, constant, constantRhs));
            control->constant->GetSubVector(dofs, constantExpected);
            double constantError = 0, constantSampleError = 0;
            for (int i = 0; i < size; ++i) constantError = std::max(constantError, std::abs(constant(i) - constantExpected(i)));
            long double interpolationMass = 0, projectionMass = 0, interpolationError = 0, projectionError = 0;
            d::ErrorMeasure interpolationMeasure, projectionMeasure;
            interpolationMeasure.volume = projectionMeasure.volume = volume;
            interpolationMeasure.targetMass = projectionMeasure.targetMass = targetMass;
            interpolationMeasure.minimumH = projectionMeasure.minimumH = minimumH;
            for (int q = 0; q < rule.GetNPoints(); ++q) {
                long double interpolatedValue = 0, projectedValue = 0, constantValue = 0;
                for (int i = 0; i < size; ++i) {
                    const double shape = control->shapes(i, q);
                    interpolatedValue += static_cast<long double>(shape) * interpolated(i);
                    projectedValue += static_cast<long double>(shape) * projected(i);
                    constantValue += static_cast<long double>(shape) * constant(i);
                }
                interpolationMass += weights(q) * interpolatedValue; projectionMass += weights(q) * projectedValue;
                interpolationError += weights(q) * (interpolatedValue - targets(q)) * (interpolatedValue - targets(q));
                projectionError += weights(q) * (projectedValue - targets(q)) * (projectedValue - targets(q));
                interpolationMeasure.signedEos += weights(q) * (interpolatedValue - targets(q));
                projectionMeasure.signedEos += weights(q) * (projectedValue - targets(q));
                interpolationMeasure.minimumRho = std::min(interpolationMeasure.minimumRho, interpolatedValue);
                projectionMeasure.minimumRho = std::min(projectionMeasure.minimumRho, projectedValue);
                constantSampleError = std::max(constantSampleError, static_cast<double>(std::abs(constantValue - 1)));
            }
            if (!std::isfinite(constantError) || !std::isfinite(constantSampleError) || constantError > 1e-10 || constantSampleError > 1e-10)
                throw std::runtime_error("Density order projection failed independent constant reconstruction.");
            control->interpolation->SetSubVector(dofs, interpolated); control->projection->SetSubVector(dofs, projected);
            interpolationMeasure.mass = interpolationMass; interpolationMeasure.eosSquared = interpolationError;
            projectionMeasure.mass = projectionMass; projectionMeasure.eosSquared = projectionError;
            d::WriteError(errorQuadrature, "eos_nodal_interpolation", control->degree, element, fem.mesh->GetAttribute(element), 26, rule, interpolationMeasure);
            d::WriteError(errorQuadrature, "eos_l2_q26", control->degree, element, fem.mesh->GetAttribute(element), 26, rule, projectionMeasure);
            const double rhsNorm = d::projection::Norm(rhs), defectNorm = d::projection::Norm(defect);
            elements << control->degree << ',' << element << ',' << fem.mesh->GetAttribute(element) << ",26," << rule.GetOrder() << ','
                << rule.GetNPoints() << ',' << size << ',' << std::quoted(control->collection->Name()) << ',' << volume << ',' << targetMass << ','
                << interpolationMass << ',' << projectionMass << ',' << interpolationError << ',' << projectionError << ','
                << rhsNorm << ',' << defectNorm << ',' << d::projection::Relative(defectNorm, rhsNorm) << ',' << constantError << ','
                << constantSampleError << ',' << interpolated.Min() << ',' << projected.Min() << ',' << projected.Max() << '\n';
            for (int i = 0; i < size; ++i) {
                const auto& point = nodalRule.IntPoint(i);
                nodes << control->degree << ',' << element << ',' << i << ',' << dofs[i] << ',' << point.x << ',' << point.y << ',' << point.z
                    << ',' << nodalH(i) << ',' << interpolated(i) << ',' << interpolated(i) << ',' << projected(i) << '\n';
            }
        }
    }
    elements.flush(); nodes.flush(); errorQuadrature.flush();
    std::vector<d::ErrorField> errorFields;
    for (const auto& control : controls) {
        errorFields.push_back({"eos_nodal_interpolation", control->degree, control->interpolation.get()});
        errorFields.push_back({"eos_l2_q26", control->degree, control->projection.get()});
    }
    const auto projectionsFinished = d::ObserverClock::now();
    d::ObserveErrorRule(state, errorFields, 30, n, K, errorQuadrature);
    const auto errorObservationFinished = d::ObserverClock::now();
    d::Outputs observations(output);
    for (const auto& control : controls) {
        std::cout << "Frozen density order Q" << control->degree << " interpolation/projection observations" << std::endl;
        const auto degree = std::to_string(control->degree);
        d::projection::WriteVector(output / ("density_order_Q" + degree + "_eos_nodal_interpolation_coefficients.txt"), *control->interpolation);
        d::projection::WriteVector(output / ("density_order_Q" + degree + "_eos_l2_q26_coefficients.txt"), *control->projection);
        d::Observe(state, *control->interpolation, control->degree, "eos_nodal_interpolation", n, K, observations);
        d::Observe(state, *control->projection, control->degree, "eos_l2_q26", n, K, observations);
    }
    checkpoint.Verify(state);
    d::WriteScope(output / "density_order_scope.txt", true, sourceDegree, n, K);
    d::WriteTimings(output / "density_order_timings.csv", observerStart, projectionsFinished,
                    errorObservationFinished, d::ObserverClock::now());
}

inline void RecordDensityOrderState(polytrope_validation::PhysicalState& state, double n, double K,
                                    const std::filesystem::path& output) {
    namespace d = density_order_detail;
    const auto observerStart = d::ObserverClock::now();
    const int degree = d::ValidateState(state, n, K);
    std::filesystem::create_directories(output);
    d::RootSelfChecks(output / "density_order_root_selfchecks.csv");
    d::SourceCheckpoint checkpoint(state);
    auto errorQuadrature = d::region::Open(output / "density_order_error_quadrature.csv");
    d::ErrorHeader(errorQuadrature);
    const std::array<d::ErrorField, 1> errorFields{{{"stored_density", degree, &state.density}}};
    const auto setupFinished = d::ObserverClock::now();
    for (int order : {26, 30}) d::ObserveErrorRule(state, errorFields, order, n, K, errorQuadrature);
    const auto errorObservationFinished = d::ObserverClock::now();
    d::Outputs observations(output);
    d::Observe(state, state.density, degree, "stored_density", n, K, observations);
    d::projection::WriteVector(output / "density_order_stored_density_coefficients.txt", state.density);
    checkpoint.Verify(state);
    d::WriteScope(output / "density_order_scope.txt", false, degree, n, K);
    d::WriteTimings(output / "density_order_timings.csv", observerStart, setupFinished,
                    errorObservationFinished, d::ObserverClock::now());
}

} // namespace experiment::plateau_diagnostics