#include <algorithm>
#include <chrono>
#include <cmath>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <limits>
#include <map>
#include <mfem.hpp>
#include <numbers>
#include <span>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>
import mean_field;
#include "plateau_enthalpy_step_diagnostics.hpp"
#include "plateau_physical_diagnostics.hpp"
#include "plateau_preconditioner_diagnostics.hpp"
#include "plateau_preconditioner_benchmark.hpp"
#include "plateau_newton_preconditioner.hpp"
#include "plateau_pressure_obstruction_diagnostics.hpp"
#include "plateau_shape_curvature_diagnostics.hpp"
#include "plateau_numerical_strategy_diagnostics.hpp"
#include "plateau_nonlinear_step_diagnostics.hpp"
#include "plateau_coupled_trial_diagnostics.hpp"
#include "plateau_coupled_trial_physical.hpp"
#include "plateau_damped_coupled_diagnostics.hpp"
#include "plateau_gravity_application_diagnostics.hpp"
#include "plateau_positivity_strategy_diagnostics.hpp"
#include "plateau_polynomial_certificate.hpp"
#include "shape_field_weak_diagnostics.hpp"
#include "reduced_surface_diagnostics.hpp"
#include "spherical_harmonic_surface_diagnostics.hpp"
#include "negative_density_diagnostics.hpp"
#include "negative_density_regions.hpp"
#include "density_order_projection.hpp"
#include "continuation_snapshot.hpp"
#include "harmonic_linear_trials.hpp"

using namespace mean_field;
namespace {
using Clock = std::chrono::steady_clock;
struct Options {
    double n = 1., angularMomentum = 0., atol = 0., rtol = 1e-8, linearTolerance = .03;
    int maxNewton = 12, maxLinear = 200, restart = 40, probeRestart = 40;
    bool probes = false, noCentral = false;
    bool strategyProbe = false, positivityProbe = false, positivityCertificate = false;
    bool preconditionerProbe = false, nonlinearStepProbe = false;
    bool weakShapeProbe = false, weakShapeFinal = false;
    bool reducedSurfaceProbe = false, reducedSurfaceAlgebraCheck = false;
    bool densityDiagnostic = false;
    bool densityRegions = false;
    bool densityOrderProjection = false, densityOrderFinal = false;
    int surfaceDegree = 4, surfaceMaximumOrder = -1;
    std::string surfacePrescription = "nodal";
    bool coupledTrialProbe = false, coupledSeedControl = false;
    bool dampedCoupledProbe = false;
    std::string coupledSearch = "near";
    int coupledBetaSteps = 4;
    double coupledBetaReduction = .5, coupledExcessRatio = 1.;
    double coupledTolerance = .03, coupledMaximumRatio = 1.;
    int coupledIterations = 160, coupledRestart = 160, coupledAttempts = 1;
    std::string preconditionerVariant = "independent";
    int schurCorrectionSweeps = 1;
    double schurRelaxation = 0.5;
    std::string strategy = "baseline";
    double gravityTolerance = 1e-10;
    int gravityIterations = 300, gravityRestart = 80;
    bool pressureProbe = false, shapeCurvature = false, surfaceTrace = false;
    std::vector<int> captureIterations;
    std::string mesh = "sandbox.smesh";
    std::filesystem::path output, probeState, probeDirection, replayCandidate;
    std::filesystem::path densitySnapshot, restartSnapshot;
    bool continuation = false, allowUnconvergedRestart = false, benchmarkTrials = false;
};
Options Parse(int argc, char **argv) {
    Options o;
    for (int i = 1; i < argc; ++i) {
        std::string a = argv[i];
        auto v = [&]() {
            if (++i >= argc)
                throw std::invalid_argument("Missing value");
            return std::string(argv[i]);
        };
        if (a == "--index")
            o.n = std::stod(v());
        else if (a == "--angular-momentum")
            o.angularMomentum = std::stod(v());
        else if (a == "--output")
            o.output = v();
        else if (a == "--mesh")
            o.mesh = v();
        else if (a == "--max-newton")
            o.maxNewton = std::stoi(v());
        else if (a == "--max-linear")
            o.maxLinear = std::stoi(v());
        else if (a == "--restart")
            o.restart = std::stoi(v());
        else if (a == "--probe-restart")
            o.probeRestart = std::stoi(v());
        else if (a == "--linear-tol")
            o.linearTolerance = std::stod(v());
        else if (a == "--atol")
            o.atol = std::stod(v());
        else if (a == "--rtol")
            o.rtol = std::stod(v());
        else if (a == "--probe-state")
            o.probeState = v();
        else if (a == "--probe-direction")
            o.probeDirection = v();
        else if (a == "--replay-candidate")
            o.replayCandidate = v();
        else if (a == "--strategy")
            o.strategy = v();
        else if (a == "--gravity-tol")
            o.gravityTolerance = std::stod(v());
        else if (a == "--gravity-iterations")
            o.gravityIterations = std::stoi(v());
        else if (a == "--gravity-restart")
            o.gravityRestart = std::stoi(v());
        else if (a == "--strategy-probe")
            o.strategyProbe = true;
        else if (a == "--positivity-probe")
            o.positivityProbe = true;
        else if (a == "--positivity-certificate")
            o.positivityCertificate = true;
        else if (a == "--preconditioner-probe")
            o.preconditionerProbe = true;
        else if (a == "--nonlinear-step-probe")
            o.nonlinearStepProbe = true;
        else if (a == "--reduced-surface-probe")
            o.reducedSurfaceProbe = true;
        else if (a == "--reduced-surface-algebra-check")
            o.reducedSurfaceAlgebraCheck = true;
        else if (a == "--density-diagnostic")
            o.densityDiagnostic = true;
        else if (a == "--density-regions") {
            o.densityDiagnostic = true;
            o.densityRegions = true;
        }
        else if (a == "--density-order-projection") {
            o.densityDiagnostic = true;
            o.densityOrderProjection = true;
        }
        else if (a == "--density-order-final")
            o.densityOrderFinal = true;
        else if (a == "--restart-snapshot")
            o.restartSnapshot = v();
        else if (a == "--continuation")
            o.continuation = true;
        else if (a == "--allow-unconverged-restart")
            o.allowUnconvergedRestart = true;
        else if (a == "--benchmark-trials")
            o.benchmarkTrials = true;
        else if (a == "--density-snapshot")
            o.densitySnapshot = v();
        else if (a == "--surface-prescription")
            o.surfacePrescription = v();
        else if (a == "--surface-maximum-order")
            o.surfaceMaximumOrder = std::stoi(v());
        else if (a == "--surface-degree")
            o.surfaceDegree = std::stoi(v());
        else if (a == "--weak-shape-probe")
            o.weakShapeProbe = true;
        else if (a == "--weak-shape-final")
            o.weakShapeFinal = true;
        else if (a == "--damped-coupled-probe")
            o.dampedCoupledProbe = true;
        else if (a == "--coupled-search")
            o.coupledSearch = v();
        else if (a == "--coupled-beta-steps")
            o.coupledBetaSteps = std::stoi(v());
        else if (a == "--coupled-beta-reduction")
            o.coupledBetaReduction = std::stod(v());
        else if (a == "--coupled-excess-ratio")
            o.coupledExcessRatio = std::stod(v());
        else if (a == "--coupled-trial-probe")
            o.coupledTrialProbe = true;
        else if (a == "--coupled-seed-control")
            o.coupledSeedControl = true;
        else if (a == "--coupled-tol")
            o.coupledTolerance = std::stod(v());
        else if (a == "--coupled-iterations")
            o.coupledIterations = std::stoi(v());
        else if (a == "--coupled-restart")
            o.coupledRestart = std::stoi(v());
        else if (a == "--coupled-max-ratio")
            o.coupledMaximumRatio = std::stod(v());
        else if (a == "--coupled-attempts")
            o.coupledAttempts = std::stoi(v());
        else if (a == "--preconditioner-variant")
            o.preconditionerVariant = v();
        else if (a == "--schur-sweeps")
            o.schurCorrectionSweeps = std::stoi(v());
        else if (a == "--schur-relaxation")
            o.schurRelaxation = std::stod(v());
        else if (a == "--pressure-probe")
            o.pressureProbe = true;
        else if (a == "--shape-curvature")
            o.shapeCurvature = true;
        else if (a == "--surface-trace")
            o.surfaceTrace = true;
        else if (a == "--probes")
            o.probes = true;
        else if (a == "--no-central")
            o.noCentral = true;
        else if (a == "--capture-iterations") {
            std::istringstream entries(v());
            std::string entry;
            while (std::getline(entries, entry, ',')) {
                size_t used = 0;
                const int iteration = std::stoi(entry, &used);
                if (used != entry.size() || iteration < 0)
                    throw std::invalid_argument("Invalid capture iteration");
                o.captureIterations.push_back(iteration);
            }
            if (o.captureIterations.empty())
                throw std::invalid_argument("Empty capture iteration list");
        }
        else
            throw std::invalid_argument("Unknown option: " + a);
    }
    if (o.output.empty() || o.n < 1 || o.n >= 5 || o.maxNewton < 0 || o.maxLinear < 1 || o.restart < 1)
        throw std::invalid_argument("Invalid options");
    if ((o.continuation || o.allowUnconvergedRestart) && o.restartSnapshot.empty())
        throw std::invalid_argument("Continuation/restart permission requires --restart-snapshot");
    if (!o.restartSnapshot.empty() && (o.densityDiagnostic || !o.probeState.empty() ||
        !o.probeDirection.empty() || !o.replayCandidate.empty() || o.surfacePrescription != "harmonics"))
        throw std::invalid_argument("Snapshot restart requires native harmonic coordinates and no legacy replay");
    if (o.benchmarkTrials && !o.preconditionerProbe)
        throw std::invalid_argument("--benchmark-trials requires --preconditioner-probe");
    if (o.surfacePrescription != "nodal" && o.surfacePrescription != "harmonics")
        throw std::invalid_argument("Surface prescriptions are nodal or harmonics");
    if (!o.densitySnapshot.empty() && !o.densityDiagnostic)
        throw std::invalid_argument("--density-snapshot requires --density-diagnostic");
    if ((o.densityOrderProjection && o.densityRegions) || (o.densityDiagnostic && o.densityOrderFinal))
        throw std::invalid_argument("Choose one frozen density observation or the coupled final observation");
    if (o.densityDiagnostic && (o.maxNewton != 0 || o.strategy != "baseline" || o.probes ||
        o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || o.positivityProbe ||
        o.positivityCertificate || o.preconditionerProbe || o.nonlinearStepProbe || o.reducedSurfaceProbe ||
        o.reducedSurfaceAlgebraCheck || o.weakShapeProbe || o.weakShapeFinal || o.coupledTrialProbe ||
        o.dampedCoupledProbe || !o.probeState.empty() || !o.probeDirection.empty() ||
        !o.replayCandidate.empty() || !o.captureIterations.empty()))
        throw std::invalid_argument("Density diagnostics require --max-newton 0 and no solver/probe strategies");
    if (o.surfacePrescription == "harmonics") {
        deformation::SphericalHarmonicSurfaceOptions surfaceOptions{.maximumDegree = o.surfaceDegree};
        if (o.surfaceMaximumOrder != -1) surfaceOptions.maximumAzimuthalOrder = o.surfaceMaximumOrder;
        surfaceOptions.validate();
        // Older frozen probes interpret snapshots as nodal radial values. A
        // size match is insufficient to establish modal snapshot provenance.
        if (o.reducedSurfaceProbe || o.reducedSurfaceAlgebraCheck || o.weakShapeProbe || o.weakShapeFinal ||
            o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe ||
            o.positivityProbe || o.positivityCertificate ||
            o.coupledTrialProbe || o.dampedCoupledProbe || !o.probeState.empty() || !o.probeDirection.empty() ||
            !o.replayCandidate.empty() || (o.strategy != "baseline" && o.strategy != "curvature"))
            throw std::invalid_argument("Native harmonic diagnostics exclude legacy nodal probes and unsupported strategies");
    } else if (o.surfaceMaximumOrder != -1) {
        throw std::invalid_argument("--surface-maximum-order requires --surface-prescription harmonics");
    }
    if (o.reducedSurfaceProbe || o.reducedSurfaceAlgebraCheck) {
        if ((o.reducedSurfaceProbe && o.reducedSurfaceAlgebraCheck) || o.maxNewton != 0 ||
            o.strategy != "baseline" || o.probes || o.pressureProbe || o.shapeCurvature ||
            o.surfaceTrace || o.strategyProbe || o.positivityProbe || o.positivityCertificate ||
            o.preconditionerProbe || o.nonlinearStepProbe || o.coupledTrialProbe || o.dampedCoupledProbe ||
            o.weakShapeProbe || o.weakShapeFinal || !o.probeDirection.empty() ||
            !o.replayCandidate.empty() || !o.captureIterations.empty() || !o.noCentral ||
            (o.surfaceDegree != 0 && o.surfaceDegree != 2 && o.surfaceDegree != 4 &&
             o.surfaceDegree != 6 && o.surfaceDegree != 8) ||
            !std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.)
            throw std::invalid_argument("Reduced-surface study requires a standalone unbordered frozen probe; degree 0 is the full control");
    }
    if (o.weakShapeProbe || o.weakShapeFinal) {
        if ((o.weakShapeProbe && o.weakShapeFinal) || o.probes || o.pressureProbe || o.shapeCurvature ||
            o.surfaceTrace || o.strategyProbe || o.positivityProbe || o.positivityCertificate ||
            o.preconditionerProbe || o.nonlinearStepProbe || o.coupledTrialProbe || o.dampedCoupledProbe ||
            !o.replayCandidate.empty() || !o.captureIterations.empty() ||
            (o.weakShapeProbe && (o.maxNewton != 0 || o.strategy != "baseline")) ||
            (o.weakShapeFinal && (o.maxNewton == 0 || !o.probeState.empty() || !o.probeDirection.empty())) ||
            !std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.)
            throw std::invalid_argument("Weak-shape diagnostics require a standalone frozen probe or a final live endpoint probe");
    }
    if (o.probeRestart != 40 && o.probeRestart != 160)
        throw std::invalid_argument("--probe-restart must be 40 or 160");
    if ((o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || o.positivityProbe || o.coupledTrialProbe || o.dampedCoupledProbe) && (o.probeState.empty() || o.probeDirection.empty()))
        throw std::invalid_argument("Focused probes require --probe-state and --probe-direction");
    if (o.positivityCertificate && o.probeState.empty())
        throw std::invalid_argument("--positivity-certificate requires --probe-state");
    if (o.positivityCertificate && (o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace ||
        o.strategyProbe || o.positivityProbe || o.preconditionerProbe || o.nonlinearStepProbe || o.coupledTrialProbe || o.dampedCoupledProbe || !o.probeDirection.empty() || !o.replayCandidate.empty()))
        throw std::invalid_argument("--positivity-certificate is a certificate-only mode; omit other probe options");
    if (o.preconditionerProbe) {
        experiment::plateau_diagnostics::ParsePreconditionerBenchmarkVariant(o.preconditionerVariant);
        if (o.maxNewton != 0 || o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace ||
            o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.nonlinearStepProbe || o.coupledTrialProbe || o.dampedCoupledProbe || !o.probeDirection.empty() ||
            !o.replayCandidate.empty() || !o.captureIterations.empty() || o.strategy != "baseline")
            throw std::invalid_argument("--preconditioner-probe requires --max-newton 0 and no nonlinear/probe strategies");
        if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.0)
            throw std::invalid_argument("A preconditioner benchmark requires a finite positive --linear-tol");
    }
    if (o.nonlinearStepProbe) {
        if (o.maxNewton != 0 || o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace ||
            o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.preconditionerProbe || o.coupledTrialProbe || o.dampedCoupledProbe ||
            !o.probeDirection.empty() || !o.replayCandidate.empty() || !o.captureIterations.empty() ||
            o.strategy != "baseline")
            throw std::invalid_argument("--nonlinear-step-probe requires --max-newton 0 and no other probes or strategy; it compares raw and curvature paths itself");
        if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.0)
            throw std::invalid_argument("A nonlinear-step probe requires a finite positive --linear-tol");
    }
    if (o.coupledTrialProbe || o.dampedCoupledProbe) {
        if (o.maxNewton != 0 || o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace ||
            o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.preconditionerProbe ||
            o.nonlinearStepProbe || (o.coupledTrialProbe && o.dampedCoupledProbe) || !o.replayCandidate.empty() || !o.captureIterations.empty() ||
            o.coupledAttempts != 1 || (o.strategy != "coupled" && o.strategy != "curvature-coupled" && o.strategy != "backtracked-coupled" && o.strategy != "curvature-backtracked-coupled"))
            throw std::invalid_argument("--coupled-trial-probe requires --max-newton 0 and strategy coupled or curvature-coupled, without other probes");
        if (o.coupledTrialProbe && o.strategy != "coupled" && o.strategy != "curvature-coupled")
            throw std::invalid_argument("The independent coupled probe requires its original full-step strategy");
        if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.)
            throw std::invalid_argument("A coupled trial probe requires a finite positive outer --linear-tol");
    }
    if (o.coupledSeedControl && !o.coupledTrialProbe) {
        throw std::invalid_argument("--coupled-seed-control requires --coupled-trial-probe");
    }
    if (!o.preconditionerProbe) {
        if (o.preconditionerVariant != "independent" && o.preconditionerVariant != "gravity-schur" &&
            o.preconditionerVariant != "stellar-ldu" && o.preconditionerVariant != "stellar-ldu-gravity-schur")
            throw std::invalid_argument("Newton preconditioner variants are independent, gravity-schur, stellar-ldu, or stellar-ldu-gravity-schur");
        if (o.preconditionerVariant != "independent") {
            if (o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe ||
                o.positivityProbe || o.positivityCertificate || (!o.probeState.empty() && !o.nonlinearStepProbe && !o.coupledTrialProbe && !o.dampedCoupledProbe && !o.weakShapeProbe && !o.reducedSurfaceProbe) ||
                (!o.probeDirection.empty() && !o.coupledTrialProbe && !o.dampedCoupledProbe && !o.weakShapeProbe && !o.reducedSurfaceProbe) ||
                !o.replayCandidate.empty())
                throw std::invalid_argument("A custom Newton preconditioner cannot be combined with frozen or certificate probes");
            if (o.strategy != "baseline" && o.strategy != "positivity" && o.strategy != "curvature" &&
                o.strategy != "positivity-curvature" && o.strategy != "coupled" && o.strategy != "curvature-coupled" && o.strategy != "backtracked-coupled" && o.strategy != "curvature-backtracked-coupled")
                throw std::invalid_argument("Unsupported custom Newton preconditioner strategy");
            if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.0 || o.schurCorrectionSweeps < 0 ||
                !std::isfinite(o.schurRelaxation) || o.schurRelaxation <= 0.0)
                throw std::invalid_argument("Invalid custom Newton preconditioner controls");
        }
    }
    if ((!o.probeDirection.empty() || !o.replayCandidate.empty()) && o.probeState.empty())
        throw std::invalid_argument("A supplied direction or candidate requires --probe-state");
    if (o.strategy != "baseline" && o.strategy != "positivity" && o.strategy != "curvature" &&
        o.strategy != "reduced" && o.strategy != "positivity-curvature" && o.strategy != "positivity-reduced" &&
        o.strategy != "coupled" && o.strategy != "curvature-coupled" && o.strategy != "backtracked-coupled" && o.strategy != "curvature-backtracked-coupled")
        throw std::invalid_argument("Unknown numerical strategy");
    if (!std::isfinite(o.gravityTolerance) || o.gravityTolerance <= 0.0 || o.gravityIterations < 1 || o.gravityRestart < 1)
        throw std::invalid_argument("Invalid gravity solve controls");
    if (!std::isfinite(o.coupledTolerance) || o.coupledTolerance <= 0. || o.coupledIterations < 1 ||
        o.coupledRestart < 1 || !std::isfinite(o.coupledMaximumRatio) || o.coupledMaximumRatio <= 0. || o.coupledAttempts < 1)
        throw std::invalid_argument("Invalid coupled correction controls");
    if (o.coupledSearch != "none" && o.coupledSearch != "full" && o.coupledSearch != "damped" && o.coupledSearch != "near")
        throw std::invalid_argument("Unknown coupled search mode");
    if (!o.dampedCoupledProbe && o.coupledSearch != "near")
        throw std::invalid_argument("--coupled-search selects controls only in --damped-coupled-probe");
    solver::BacktrackedCoupledStep{.maximumSteps=o.coupledBetaSteps, .reduction=o.coupledBetaReduction}.Validate();
    solver::NearAcceptanceCoupledAttempt{.maximumExcessRatio=o.coupledExcessRatio}.Validate();
    if (std::filesystem::exists(o.output))
        throw std::runtime_error("Output directory already exists");
    std::filesystem::create_directories(o.output);
    return o;
}
std::ofstream File(const std::filesystem::path &p) {
    std::ofstream s(p);
    s.exceptions(std::ios::failbit | std::ios::badbit);
    s << std::setprecision(17);
    return s;
}
double Norm(std::span<const mfem::real_t> v, int offset, int size) {
    long double sum = 0;
    for (int i = offset; i < offset + size; ++i)
        sum += (long double)v[i] * v[i];
    return std::sqrt(sum);
}
double Norm(const mfem::Vector &v, int offset, int size) {
    return Norm(std::span<const mfem::real_t>{v.GetData(), (size_t)v.Size()}, offset, size);
}
struct Scales {
    double K, rho, xi, mu;
};
Scales Derive(double n) {
    if (n == 1.)
        return {2. * utils::G / std::numbers::pi, std::numbers::pi / 4., std::numbers::pi, std::numbers::pi};
    auto le = seed::integrateLaneEmden(n, 2000., 1e-3);
    if (!le.firstZeroCoordinate)
        throw std::runtime_error("No Lane Emden root");
    const double xi = *le.firstZeroCoordinate, mu = -xi * xi * le.thetaDerivative(le.thetaDerivative.Size() - 1);
    const double rho = xi * xi * xi / (4 * std::numbers::pi * mu);
    return {4 * std::numbers::pi * utils::G * std::pow(rho, 1. - 1. / n) / ((n + 1) * xi * xi), rho, xi, mu};
}
void SaveVector(const std::filesystem::path &p, const mfem::Vector &v) {
    auto s = File(p);
    s << v.Size() << '\n';
    v.Print(s, 1);
}
mfem::Vector LoadVector(const std::filesystem::path &path, int expectedSize) {
    // Legacy nodal probes may share a total vector size with another shape
    // space. Honor known coordinate provenance before reading its numbers.
    for (const auto &directory : {path.parent_path(), path.parent_path().parent_path()}) {
        std::ifstream metadata(directory / "metadata.txt");
        std::string line;
        while (std::getline(metadata, line))
            if (line == "surface_prescription=harmonics")
                throw std::invalid_argument("Legacy nodal probes cannot load a harmonic snapshot: " + path.string());
    }
    std::ifstream input(path);
    int count = -1;
    input >> count;
    if (!input || count != expectedSize)
        throw std::runtime_error("Snapshot size mismatch: " + path.string());
    mfem::Vector result(count);
    for (int i = 0; i < count; ++i) {
        input >> result(i);
        if (!input || !std::isfinite(result(i)))
            throw std::runtime_error("Invalid snapshot: " + path.string());
    }
    std::string extra;
    if (input >> extra)
        throw std::runtime_error("Unexpected trailing snapshot data: " + path.string());
    return result;
}

template <class Problem>
mfem::Vector LoadDensitySnapshot(const Problem &problem, const Options &o, const Scales &scales) {
    std::ifstream metadata(o.densitySnapshot / "metadata.txt");
    if (!metadata) throw std::invalid_argument("Density replay requires snapshot metadata.");
    std::map<std::string, std::string> saved;
    std::string line;
    while (std::getline(metadata, line)) {
        const auto separator = line.find('=');
        if (separator != std::string::npos) saved.emplace(line.substr(0, separator), line.substr(separator + 1));
    }
    const auto number = [&](const std::string &key, double expected) {
        const double actual = std::stod(saved.at(key));
        if (!std::isfinite(actual) || std::abs(actual - expected) > 1e-13 * std::max(1., std::abs(expected)))
            throw std::invalid_argument("Density snapshot model mismatch: " + key);
    };
    number("n", o.n); number("K", scales.K); number("J", o.angularMomentum);
    number("M", 1.); number("R_reference", 1.); number("central_phase", !o.noCentral);
    if (saved.at("surface_prescription") != o.surfacePrescription)
        throw std::invalid_argument("Density snapshot surface prescription mismatch.");
    if (o.surfacePrescription == "harmonics") {
        number("surface_degree", o.surfaceDegree);
        number("surface_maximum_order", o.surfaceMaximumOrder);
    }
    // A matching vector size does not establish matching coordinates. Require
    // every named state and residual block to match the current manifest.
    std::map<std::string, std::pair<int, int>> expected;
    for (const auto &block : problem.GetManifest().valueBlocks())
        expected.emplace("value," + std::string(block.stableId), std::pair{block.offset, block.size});
    for (const auto &block : problem.GetManifest().residualBlocks())
        expected.emplace("residual," + std::string(block.stableId), std::pair{block.offset, block.size});
    std::ifstream layout(o.densitySnapshot / "layout.csv");
    if (!std::getline(layout, line) || line != "kind,block,offset,size,scale")
        throw std::invalid_argument("Density replay requires the saved block layout.");
    while (std::getline(layout, line)) {
        std::istringstream row(line);
        std::array<std::string, 5> fields;
        for (auto &field : fields) if (!std::getline(row, field, ','))
            throw std::invalid_argument("Malformed saved density layout.");
        const auto key = fields[0] + "," + fields[1];
        const auto found = expected.find(key);
        if (found == expected.end() || found->second != std::pair{std::stoi(fields[2]), std::stoi(fields[3])})
            throw std::invalid_argument("Density snapshot coordinate layout mismatch: " + key);
        expected.erase(found);
    }
    if (!expected.empty() || layout.bad() || metadata.bad())
        throw std::invalid_argument("Incomplete density snapshot provenance.");
    return experiment::plateau_diagnostics::coupled_physical_detail::Load(
        o.densitySnapshot / "accepted_state.txt", problem.StateSize());
}

template <class State>
void FocusedProbes(State &s, const fem::FEM &f, const Options &o, const mfem::Vector &physicalBase) {
    mfem::Vector base;
    s.normalizedOperator->NormalizeState(physicalBase, base);
    auto restore = [&]() {
        s.trialNormalizedState = base;
        if (!s.PrepareTrial())
            throw std::runtime_error("Focused probe base is inadmissible");
    };
    restore();
    const mfem::Vector residual(s.trialNormalizedResidual);
    SaveVector(o.output / "frozen_physical_base.txt", physicalBase);
    SaveVector(o.output / "frozen_normalized_base.txt", base);
    SaveVector(o.output / "frozen_normalized_residual.txt", residual);
    auto layout = File(o.output / "layout.csv");
    layout << "kind,block,offset,size,scale\n";
    for (const auto &b : s.Problem().GetManifest().residualBlocks())
        layout << "residual," << b.stableId << ',' << b.offset << ',' << b.size << ',' << b.scale << '\n';
    for (const auto &b : s.Problem().GetManifest().valueBlocks())
        layout << "value," << b.stableId << ',' << b.offset << ',' << b.size << ',' << b.scale << '\n';
    auto checks = File(o.output / "frozen_input_checks.csv");
    checks << "metric,value\nbase_residual_norm," << residual.Norml2() << '\n';
    const auto capturedBase = o.probeState.parent_path() / "normalized_base.txt";
    if (std::filesystem::is_regular_file(capturedBase)) {
        mfem::Vector difference = LoadVector(capturedBase, base.Size());
        difference -= base;
        checks << "captured_normalized_base_difference," << difference.Norml2() << '\n';
        if (difference.Norml2() > 64. * std::numeric_limits<double>::epsilon() * std::max(1., base.Norml2()))
            throw std::runtime_error("Frozen normalization disagrees with the captured base");
    }
    const auto capturedResidual = o.probeState.parent_path() / "normalized_residual.txt";
    if (std::filesystem::is_regular_file(capturedResidual)) {
        mfem::Vector difference = LoadVector(capturedResidual, residual.Size());
        difference -= residual;
        checks << "captured_residual_difference," << difference.Norml2() << '\n';
    }
    if (!o.replayCandidate.empty()) {
        s.trialNormalizedState = LoadVector(o.replayCandidate, base.Size());
        const auto prepared = s.PrepareTrial();
        auto replay = File(o.output / "replayed_candidate.csv");
        replay << "admissible,reason,stage,eos_code,minimum_jacobian,residual_norm\n";
        if (prepared)
            replay << "1,-1,-1,-1," << prepared->minimumJacobianDeterminant << ','
                   << s.trialNormalizedResidual.Norml2() << '\n';
        else
            replay << "0," << int(prepared.error().reason) << ',' << int(prepared.error().stage) << ','
                   << int(prepared.error().thermodynamicErrorCode) << ','
                   << prepared.error().minimumJacobianDeterminant << ",nan\n";
        restore();
    }
    if (!o.probeDirection.empty()) {
        const auto direction = LoadVector(o.probeDirection, base.Size());
        mfem::Vector action(base.Size());
        restore();
        s.normalizedOperator->Mult(direction, action);
        SaveVector(o.output / "probe_normalized_correction.txt", direction);
        checks << "direction_norm," << direction.Norml2() << "\naction_norm," << action.Norml2() << '\n';
        action += residual;
        checks << "linear_defect_norm," << action.Norml2() << "\nlinear_defect_relative,"
               << action.Norml2() / std::max(residual.Norml2(), 1e-300) << '\n';
        checks.flush();
        s.normalizedCorrection = direction;
        auto bounds = File(o.output / "frozen_enthalpy_bound.csv");
        bounds << "metric,value\n";
        for (const auto &[key, value] : experiment::plateau_diagnostics::MeasureEnthalpyStep(s, f, o.n, &physicalBase))
            bounds << key << ',' << value << '\n';
        bounds.flush();
        if (o.strategyProbe) {
            s.physicalInverse->Refresh();
            using Form = typename std::remove_cvref_t<decltype(s.Problem())>::FormType;
            solver::PreparedGravitySubsystem gravity(*s.normalizedOperator, *s.normalizedInverse,
                solver::MakeGravityBlockPartition<Form>(s.Problem().GetManifest().layout()),
                s.Problem().GetCommunicator());
            plateau_diagnostics::ProbeGravityApplications(gravity, o.output);
            experiment::plateau_diagnostics::ProbeNumericalStrategies(s, f, base, direction, o.output);
        }
        if (o.positivityProbe)
            experiment::plateau_diagnostics::ProbePolynomialEnthalpy(s, f, base, direction, o.output);
        if (o.surfaceTrace)
            experiment::plateau_diagnostics::ProbeSurfaceTrace(s, f, o.n, base, direction, o.output);
        if (o.pressureProbe)
            experiment::plateau_diagnostics::ProbePressureObstruction(s, f, o.n, base, direction, o.output);
        if (o.shapeCurvature)
            experiment::plateau_diagnostics::ProbeShapeCurvature(s, f, o.n, base, direction, o.output);
    }
    restore();
    mfem::Vector difference(s.trialNormalizedResidual);
    difference -= residual;
    checks << "final_restoration_residual_difference," << difference.Norml2() << '\n';
}

template <class State>
void FrozenProbes(State &s, const fem::FEM &f, const Options &o, const mfem::Vector *physicalBase = nullptr) {
    auto rows = s.Problem().GetManifest().residualBlocks();
    mfem::Vector frozenState(s.acceptedNormalizedState);
    if (physicalBase)
        s.normalizedOperator->NormalizeState(*physicalBase, frozenState);
    auto restoreFrozen = [&]() {
        s.trialNormalizedState = frozenState;
        auto p = s.PrepareTrial();
        if (!p)
            throw std::runtime_error("Frozen diagnostic state is inadmissible");
    };
    restoreFrozen();
    s.physicalInverse->Refresh();
    const mfem::Vector base(s.trialNormalizedResidual);
    if (base.Norml2() == 0.)
        return;
    auto results = File(o.output / "frozen_linear.csv");
    results << "restart,tolerance,iterations,converged,true_relative,seconds\n";
    auto histories = File(o.output / "frozen_history.csv");
    histories << "restart,tolerance,iteration,reported_norm,final\n";
    auto blocks = File(o.output / "frozen_blocks.csv");
    blocks << "restart,tolerance,block,F_norm,Jp_plus_F_norm,global_correction_norm\n";
    mfem::Vector selectedCorrection, baselineAction(base.Size());
    // These are two cold direct MFEM solves, not the configured production
    // backend or any recycled Krylov history.
    for (int restart : {40, 160}) {
        mfem::FGMRESSolver linear(MPI_COMM_WORLD);
        solver::ResidualHistoryMonitor monitor;
        linear.SetPreconditioner(*s.normalizedInverse);
        linear.SetOperator(*s.normalizedOperator);
        linear.SetKDim(restart);
        linear.SetRelTol(0.);
        linear.SetAbsTol(o.linearTolerance * base.Norml2());
        linear.SetMaxIter(o.maxLinear);
        linear.SetPrintLevel(-1);
        linear.SetMonitor(monitor);
        linear.iterative_mode = false;
        mfem::Vector b(base), p(base.Size()), action(base.Size());
        b *= -1.;
        p = 0.;
        const auto start = Clock::now();
        linear.Mult(b, p);
        double seconds = std::chrono::duration<double>(Clock::now() - start).count();
        s.normalizedOperator->Mult(p, action);
        action += base;
        const double rel = action.Norml2() / base.Norml2();
        results << restart << ',' << o.linearTolerance << ',' << linear.GetNumIterations() << ','
                << linear.GetConverged() << ',' << rel << ',' << seconds << '\n';
        results.flush();
        for (const auto &h : monitor.GetHistory())
            histories << restart << ',' << o.linearTolerance << ',' << h.iteration << ',' << h.reportedNorm << ','
                      << h.final << '\n';
        histories.flush();
        for (const auto &r : rows)
            blocks << restart << ',' << o.linearTolerance << ',' << r.stableId << ',' << Norm(base, r.offset, r.size)
                   << ',' << Norm(action, r.offset, r.size) << ',' << p.Norml2() << '\n';
        blocks.flush();
        std::cout << "Frozen restart=" << restart << " true=" << rel << " seconds=" << seconds << std::endl;
        if (restart == o.probeRestart)
            selectedCorrection = p;
    }
    std::cout << "Frozen probe direction restart=" << o.probeRestart << std::endl;
    s.normalizedCorrection = selectedCorrection;
    s.normalizedOperator->Mult(selectedCorrection, baselineAction);
    auto positivity = File(o.output / "enthalpy_step_bound.csv");
    positivity << "metric,value\n";
    for (const auto &[key, value] : experiment::plateau_diagnostics::MeasureEnthalpyStep(s, f, o.n, physicalBase))
        positivity << key << ',' << value << '\n';
    positivity.flush();
    SaveVector(o.output / "probe_normalized_correction.txt", selectedCorrection);
    auto fd = File(o.output / "finite_differences.csv");
    fd << "epsilon,admissible,block,absolute_error,relative_to_Jp,action_norm\n";
    const auto safe = s.EstimateLargestSafeStepSize(1., .9);
    if (!std::isfinite(safe.stepSize) || safe.stepSize <= 0.)
        return;
    auto trials = File(o.output / "frozen_trials.csv");
    trials << "alpha,admissible,rejection_reason,F_norm,merit_ratio\n";
    // Diagnostic trials are never committed. Each samples the same frozen Newton direction.
    for (double a : {1., .5, .1, .01, .001, .0001, .00001}) {
        const double epsilon = std::min(1., safe.stepSize) * a;
        s.trialNormalizedState = frozenState;
        s.trialNormalizedState.Add(epsilon, selectedCorrection);
        const auto prep = s.PrepareTrial();
        if (!prep) {
            trials << epsilon << ",0," << int(prep.error().reason) << ",nan,nan\n";
            fd << epsilon << ",0,all,nan,nan,nan\n";
        } else {
            trials << epsilon << ",1,-1," << s.trialNormalizedResidual.Norml2() << ','
                   << std::pow(s.trialNormalizedResidual.Norml2() / base.Norml2(), 2) << '\n';
            mfem::Vector error(s.trialNormalizedResidual);
            error -= base;
            error /= epsilon;
            error -= baselineAction;
            for (const auto &r : rows) {
                double e = Norm(error, r.offset, r.size), j = Norm(baselineAction, r.offset, r.size);
                fd << epsilon << ",1," << r.stableId << ',' << e << ',' << e / std::max(j, 1e-300) << ',' << j << '\n';
            }
        }
        trials.flush();
        fd.flush();
        restoreFrozen();
    }
}

template <bool FullDefaultStack, class Context>
int RunContext(Context &context, const Options &o, const Scales &scales) {
    std::vector<operators::RootBlockDescriptor> rows, values;
    mfem::Vector residualFactors, stateFactors;
    std::vector<int> centralRows;
    int centralOffset = -1;
    if (!o.restartSnapshot.empty()) {
        mfem::Vector restartPhysical;
        solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &) {
            const auto directory = std::filesystem::is_directory(o.restartSnapshot) ? o.restartSnapshot :
                (std::filesystem::exists(o.restartSnapshot.parent_path() / "metadata.txt") ?
                 o.restartSnapshot.parent_path() : o.restartSnapshot.parent_path().parent_path());
            const auto saved = experiment::continuation::ParseSnapshotMetadata(directory);
            const auto number = [&](const char *key, double expected) {
                const double actual = std::stod(saved.at(key));
                if (!std::isfinite(actual) || std::abs(actual - expected) >
                    1e-13 * std::max(1., std::abs(expected)))
                    throw std::invalid_argument(std::string("Restart metadata mismatch: ") + key);
            };
            if (saved.at("surface_prescription") != o.surfacePrescription)
                throw std::invalid_argument("Restart surface coordinate prescription differs");
            number("surface_degree", o.surfaceDegree); number("surface_maximum_order", o.surfaceMaximumOrder);
            number("density_order", field::Density::scalarOrder); number("enthalpy_order", field::Enthalpy::scalarOrder);
            number("gravity_potential_order", field::Gravity::potentialOrder);
            number("gravity_flux_order", field::Gravity::fluxOrder);
            number("displacement_order", field::Displacement::Vector::familyOrder);
            number("M", 1.); number("R_reference", 1.);
            if (!o.continuation) {
                number("n", o.n); number("K", scales.K); number("J", o.angularMomentum);
                number("central_phase", !o.noCentral);
            }
            restartPhysical = experiment::continuation::LoadMappedSnapshot(s.Problem(),
                o.restartSnapshot, o.continuation, o.allowUnconvergedRestart);
        });
        solver::detail::StellarEquilibriumContextDiagnostics::InitializePhysicalState(context, restartPhysical);
    }
    const bool contextUsesGravitySchur = !o.preconditionerProbe &&
        (o.preconditionerVariant == "gravity-schur" || o.preconditionerVariant == "stellar-ldu-gravity-schur");
    auto metadata = File(o.output / "metadata.txt");
    metadata << "n=" << o.n << "\nK=" << scales.K << "\nrho_c=" << scales.rho << "\nJ=" << o.angularMomentum
             << "\nM=1\nR_reference=1\nxi1=" << scales.xi << "\nmu1=" << scales.mu << "\ncentral_phase=" << !o.noCentral
             << "\nrestart=" << o.restart << "\nprobe_restart=" << o.probeRestart
             << "\nfrozen_probe_backend=" << (o.preconditionerProbe ? "direct_mfem_fgmres_preconditioner_benchmark" :
                 (o.coupledTrialProbe || o.dampedCoupledProbe) ? "supplied_outer_direction_verified_with_full_J_and_auxiliary_fgmres" :
                 (o.reducedSurfaceProbe || o.nonlinearStepProbe || ((o.weakShapeProbe || o.weakShapeFinal) && o.probeDirection.empty()))
                     ? "prepared_fgmres_cold_with_independent_true_gate" :
                 o.positivityCertificate ? "polynomial_certificate_no_linear_solve" :
                 o.probeDirection.empty() ? "direct_mfem_fgmres_cold" : "supplied_direction_no_linear_solve")
             << "\nlinear_tol=" << o.linearTolerance << "\nmax_linear=" << o.maxLinear
             << "\nmax_newton=" << o.maxNewton << "\natol=" << o.atol << "\nrtol=" << o.rtol
             << "\nstrategy=" << o.strategy << "\ngravity_tolerance=" << o.gravityTolerance
             << "\ngravity_iterations=" << o.gravityIterations << "\ngravity_restart=" << o.gravityRestart
             << "\npreconditioner_probe=" << o.preconditionerProbe
             << "\nnonlinear_step_probe=" << o.nonlinearStepProbe
             << "\nsurface_prescription=" << o.surfacePrescription
             << "\nsurface_maximum_order=" << o.surfaceMaximumOrder
             << "\nreduced_surface_probe=" << o.reducedSurfaceProbe << "\nsurface_degree=" << o.surfaceDegree
             << "\nweak_shape_probe=" << o.weakShapeProbe << "\nweak_shape_final=" << o.weakShapeFinal
             << "\ndamped_coupled_probe=" << o.dampedCoupledProbe
             << "\ncoupled_search=" << o.coupledSearch
             << "\ncoupled_beta_steps=" << o.coupledBetaSteps
             << "\ncoupled_beta_reduction=" << o.coupledBetaReduction
             << "\ncoupled_excess_ratio=" << o.coupledExcessRatio
             << "\ncoupled_trial_probe=" << o.coupledTrialProbe
             << "\ncoupled_seed_control=" << o.coupledSeedControl
             << "\ncoupled_tolerance=" << o.coupledTolerance
             << "\ncoupled_iterations=" << o.coupledIterations
             << "\ncoupled_restart=" << o.coupledRestart
             << "\ncoupled_maximum_ratio=" << o.coupledMaximumRatio
             << "\ncoupled_attempts_per_iteration=" << o.coupledAttempts
             << "\npreconditioner_variant=" << o.preconditionerVariant
             << "\ncontext_preconditioner_variant=" << (o.preconditionerProbe ? "independent" : o.preconditionerVariant)
             << "\ncontext_gravity_schur_sweeps="
             << (contextUsesGravitySchur ? o.schurCorrectionSweeps : 0)
             << "\ncontext_gravity_schur_relaxation=" << (contextUsesGravitySchur ? o.schurRelaxation : 0.0)
             << "\nschur_correction_sweeps=" << o.schurCorrectionSweeps << "\nschur_relaxation=" << o.schurRelaxation
             << "\ncompiler=" << __VERSION__ << "\ncompiled=" << __DATE__ << ' ' << __TIME__ << '\n';
    metadata << "restart_snapshot=" << o.restartSnapshot.string() << "\ncontinuation=" << o.continuation
             << "\nallow_unconverged_restart=" << o.allowUnconvergedRestart
             << "\nbenchmark_trials=" << o.benchmarkTrials << '\n';
    metadata << "density_order=" << field::Density::scalarOrder
             << "\nenthalpy_order=" << field::Enthalpy::scalarOrder
             << "\ngravity_potential_order=" << field::Gravity::potentialOrder
             << "\ngravity_flux_order=" << field::Gravity::fluxOrder
             << "\ndisplacement_order=" << field::Displacement::Vector::familyOrder
             << "\ndensity_order_projection=" << o.densityOrderProjection
             << "\ndensity_order_final=" << o.densityOrderFinal << '\n';
    if (o.densityDiagnostic) {
        metadata << "density_diagnostic=1\ndensity_region_integral=" << o.densityRegions << "\ndensity_source="
                 << (o.densitySnapshot.empty() ? "lane_emden_seed" : o.densitySnapshot.string()) << '\n';
        metadata.flush();
        solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &f) {
            namespace d = experiment::plateau_diagnostics::coupled_physical_detail;
            const mfem::Vector accepted(s.AcceptedPhysicalState()), normalized(s.acceptedNormalizedState),
                residual(s.acceptedNormalizedResidual);
            const auto snapshot = o.densitySnapshot.empty() ? mfem::Vector(accepted) : LoadDensitySnapshot(s.Problem(), o, scales);
            s.normalizedOperator->NormalizeState(snapshot, s.trialNormalizedState);
            if (!s.PrepareTrial()) throw std::runtime_error("Density snapshot cannot be prepared by the physical operator.");
            d::PreparedSnapshotView view(s);
            experiment::polytrope_validation::PhysicalState physical(view, f);
            std::cout << "Frozen density diagnostic F=" << s.trialNormalizedResidual.Norml2() << std::endl;
            SaveVector(o.output / "observed_state.txt", s.normalizedOperator->GetPhysicalState());
            auto layout = File(o.output / "layout.csv");
            layout << "kind,block,offset,size,scale\n";
            for (const auto &block : s.Problem().GetManifest().valueBlocks())
                layout << "value," << block.stableId << ',' << block.offset << ',' << block.size << ',' << block.scale << '\n';
            for (const auto &block : s.Problem().GetManifest().residualBlocks())
                layout << "residual," << block.stableId << ',' << block.offset << ',' << block.size << ',' << block.scale << '\n';
            auto physicalOutput = File(o.output / "density_physical.csv");
            physicalOutput << "quadrature,metric,value\n";
            for (const auto &[key, value] : experiment::plateau_diagnostics::MeasurePhysical(physical, o.n, scales.K, 18))
                physicalOutput << "18," << key << ',' << value << '\n';
            physicalOutput.flush();
            const auto &closure = s.Problem().GetPhysicalOperator().GetBarotropicClosureOperator();
            auto ruleOutput = File(o.output / "production_closure_rules.csv");
            ruleOutput << "element,quadrature_order,quadrature_points\n";
            closure.VisitMappedGeometryRules([&](int element, const mfem::IntegrationRule &rule) {
                ruleOutput << element << ',' << rule.GetOrder() << ',' << rule.GetNPoints() << '\n';
            });
            mfem::Vector closureResidual;
            closure.BuildResidual(closureResidual);
            mfem::ParGridFunction closureFull(f.densityFes.get());
            const auto adapter = field::make_field_dof_grid_function_adapter<field::Density,
                utils::domain::CoreEnvelopeVacuumDomainSchema>(*f.densityFes);
            adapter.scatter(closureResidual, closureFull);
            SaveVector(o.output / "production_closure_residual_full.txt", closureFull);
            SaveVector(o.output / "observed_density_coefficients.txt", physical.density);
            if (o.densityOrderProjection)
                experiment::plateau_diagnostics::RecordDensityOrderProjection(physical, o.n, scales.K, o.output);
            else if (o.densityRegions)
                experiment::plateau_diagnostics::RecordDensityNegativeRegions(physical, o.n, scales.K, o.output);
            else
                experiment::plateau_diagnostics::RecordNegativeDensity(physical, o.n, scales.K, o.output);
            if (!d::Identical(accepted, s.AcceptedPhysicalState()) || !d::Identical(normalized, s.acceptedNormalizedState) ||
                !d::Identical(residual, s.acceptedNormalizedResidual))
                throw std::runtime_error("Density diagnostic changed the accepted state.");
            s.RestoreAccepted();
            auto invariants = File(o.output / "density_invariants.csv");
            invariants << "metric,value\naccepted_physical_unchanged," << d::Identical(accepted, s.AcceptedPhysicalState())
                << "\naccepted_normalized_unchanged," << d::Identical(normalized, s.acceptedNormalizedState)
                << "\naccepted_residual_unchanged," << d::Identical(residual, s.acceptedNormalizedResidual)
                << "\nprepared_residual_norm," << physical.normalizedBorderedResidualNorm << '\n';
            if (!d::Identical(accepted, s.AcceptedPhysicalState())) throw std::runtime_error("Density restoration failed.");
        });
        return 0;
    }
    if (o.weakShapeProbe || o.weakShapeFinal)
        metadata << "weak_shape_direction_space=normalized_current_factors\n"
                 << "weak_shape_operator=full_monolithic\n"
                 << "weak_shape_trial_policy=raw_prepare_only\n";
    metadata.flush();
    if (o.reducedSurfaceProbe) {
        metadata << "reduced_surface_input=" << (o.probeState.empty() ? "lane_emden_seed" : "physical_snapshot") << '\n'
                 << "reduced_surface_trial_policy=raw_prepare_only_no_commit\n"
                 << "probe_state=" << o.probeState << '\n';
        metadata.flush();
        solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &f) {
            mfem::Vector base(s.acceptedNormalizedState);
            if (!o.probeState.empty()) {
                const auto physical = LoadVector(o.probeState, s.AcceptedPhysicalState().Size());
                s.normalizedOperator->NormalizeState(physical, base);
            }
            experiment::plateau_diagnostics::ProbeReducedSurface(s, f, base, o.output, o.surfaceDegree, o.n, scales.K, {
                .linear = {.relativeTolerance = o.linearTolerance, .absoluteTolerance = 0.,
                           .maximumIterations = o.maxLinear}, .restartLength = o.restart});
        });
        return 0;
    }
    if (o.weakShapeProbe) {
        metadata << "weak_shape_input=" << (o.probeState.empty() ? "lane_emden_seed" : "physical_snapshot") << '\n'
                 << "probe_state=" << o.probeState << "\nprobe_direction=" << o.probeDirection << '\n';
        metadata.flush();
        solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &f) {
            mfem::Vector base(s.acceptedNormalizedState), direction;
            if (!o.probeState.empty()) {
                const auto physical = LoadVector(o.probeState, s.AcceptedPhysicalState().Size());
                s.normalizedOperator->NormalizeState(physical, base);
            }
            if (!o.probeDirection.empty()) direction = LoadVector(o.probeDirection, base.Size());
            experiment::plateau_diagnostics::ProbeWeakShapeField(s, f, base,
                o.probeDirection.empty() ? nullptr : &direction, o.output, {
                    .linear = {.relativeTolerance = o.linearTolerance, .absoluteTolerance = 0.,
                               .maximumIterations = o.maxLinear}, .restartLength = o.restart});
        });
        return 0;
    }
    if (o.dampedCoupledProbe) {
        metadata << "coupled_trial_input=physical_snapshot_and_supplied_normalized_direction\n"
                 << "coupled_trial_admissibility=sampled_enthalpy_without_repair\n"
                 << "coupled_trial_alphas_are_independent=0\n"
                 << "probe_state=" << o.probeState << "\nprobe_direction=" << o.probeDirection << '\n';
        metadata.flush();
        solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &f) {
            const auto physical = LoadVector(o.probeState, s.AcceptedPhysicalState().Size());
            mfem::Vector base;
            s.normalizedOperator->NormalizeState(physical, base);
            const auto direction = LoadVector(o.probeDirection, base.Size());
            experiment::plateau_diagnostics::ProbeDampedCoupled(s, base, direction, o.output, {
                .common = {.outerLinearTolerance = o.linearTolerance,
                    .gravity = {.control = {.relativeTolerance = o.gravityTolerance, .absoluteTolerance = 1e-14,
                                            .maximumIterations = o.gravityIterations},
                                .backend = {.restartLength = o.gravityRestart, .printLevel = -1}},
                    .coupled = {.control = {.relativeTolerance = o.coupledTolerance, .absoluteTolerance = 1e-14,
                                            .maximumIterations = o.coupledIterations},
                                .backend = {.restartLength = o.coupledRestart, .printLevel = -1},
                                .maximumCorrectionRatio = o.coupledMaximumRatio},
                    .gravityPreconditionPath = o.strategy.starts_with("curvature-")},
                .mode = o.coupledSearch,
                .search = {.maximumSteps = o.coupledBetaSteps, .reduction = o.coupledBetaReduction},
                .attempt = {.maximumExcessRatio = o.coupledExcessRatio}});
            experiment::plateau_diagnostics::MeasureCoupledTrialSnapshots(s, f, o.output, o.n, scales.K);
        });
        return 0;
    }
    if (o.coupledTrialProbe) {
        metadata << "coupled_trial_input=physical_snapshot_and_supplied_normalized_direction\n"
                 << "coupled_trial_admissibility=sampled_enthalpy_without_repair\n"
                 << "coupled_trial_alphas_are_independent=1\n"
                 << "probe_state=" << o.probeState << "\nprobe_direction=" << o.probeDirection << '\n';
        metadata.flush();
        solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &f) {
            const auto physical = LoadVector(o.probeState, s.AcceptedPhysicalState().Size());
            mfem::Vector base;
            s.normalizedOperator->NormalizeState(physical, base);
            const auto direction = LoadVector(o.probeDirection, base.Size());
            experiment::plateau_diagnostics::ProbeCoupledTrial(s, base, direction, o.output, {
                .outerLinearTolerance = o.linearTolerance,
                .gravity = {.control = {.relativeTolerance = o.gravityTolerance, .absoluteTolerance = 1e-14,
                                        .maximumIterations = o.gravityIterations},
                            .backend = {.restartLength = o.gravityRestart, .printLevel = -1}},
                .coupled = {.control = {.relativeTolerance = o.coupledTolerance, .absoluteTolerance = 1e-14,
                                        .maximumIterations = o.coupledIterations},
                            .backend = {.restartLength = o.coupledRestart, .printLevel = -1},
                            .maximumCorrectionRatio = o.coupledMaximumRatio},
                .gravityPreconditionPath = o.strategy == "curvature-coupled",
                .seedControl = o.coupledSeedControl});
            experiment::plateau_diagnostics::MeasureCoupledTrialSnapshots(s, f, o.output, o.n, scales.K);
        });
        return 0;
    }
    if (o.nonlinearStepProbe) {
        metadata << "nonlinear_step_input=" << (o.probeState.empty() ? "lane_emden_seed" : "physical_snapshot")
                 << "\nnonlinear_step_paths=raw,curvature\nnonlinear_step_positivity_repair=0\n";
        if (!o.probeState.empty()) metadata << "probe_state=" << o.probeState << '\n';
        metadata.flush();
        solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &) {
            mfem::Vector base(s.acceptedNormalizedState);
            if (!o.probeState.empty()) {
                const auto snapshot = LoadVector(o.probeState, s.AcceptedPhysicalState().Size());
                s.normalizedOperator->NormalizeState(snapshot, base);
            }
            experiment::plateau_diagnostics::ProbeNonlinearStep(s, base, o.output, {
                .linear = {.relativeTolerance = o.linearTolerance, .absoluteTolerance = 0.,
                           .maximumIterations = o.maxLinear},
                .restartLength = o.restart,
                .gravity = {.control = {.relativeTolerance = o.gravityTolerance, .absoluteTolerance = 1e-14,
                                        .maximumIterations = o.gravityIterations},
                            .backend = {.restartLength = o.gravityRestart, .printLevel = -1}}});
        });
        return 0;
    }
    // Custom Newton recipes do not instantiate the nine-variant frozen suite.
    if constexpr (FullDefaultStack) if (o.preconditionerProbe) {
        metadata << "benchmark_input=" << (!o.restartSnapshot.empty() ? "checked_restart_snapshot" :
            o.probeState.empty() ? "lane_emden_seed" : "physical_snapshot") << '\n';
        if (!o.probeState.empty()) metadata << "probe_state=" << o.probeState << '\n';
        solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &f) {
            if (o.probeState.empty()) {
                s.trialNormalizedState = s.acceptedNormalizedState;
            } else {
                const auto snapshot = LoadVector(o.probeState, s.AcceptedPhysicalState().Size());
                s.normalizedOperator->NormalizeState(snapshot, s.trialNormalizedState);
            }
            const auto preparationStart = Clock::now();
            if (!s.PrepareTrial())
                throw std::runtime_error("Preconditioner benchmark state is inadmissible to the physical operator");
            metadata << "benchmark_frozen_preparation_seconds="
                     << std::chrono::duration<double>(Clock::now() - preparationStart).count() << '\n';
            metadata.flush();
            experiment::plateau_diagnostics::BenchmarkPreconditioners(s, o.output,
                {.variant = o.preconditionerVariant, .restart = o.restart,
                 .maximumIterations = o.maxLinear, .relativeTolerance = o.linearTolerance,
                 .schurCorrectionSweeps = o.schurCorrectionSweeps, .schurRelaxation = o.schurRelaxation});
            if (o.benchmarkTrials)
                experiment::plateau_diagnostics::RecordHarmonicLinearTrials(s, f, o.output, o.n);
        });
        // A missed linear gate is a measured outcome, not failure to execute
        // the experiment. The CSV retains both recurrence and true gate results.
        return 0;
    }
    if constexpr (FullDefaultStack) if (!o.probeState.empty()) {
        if (o.maxNewton != 0)
            throw std::invalid_argument("--probe-state requires --max-newton 0");
        metadata << "probe_state=" << o.probeState << '\n';
        if (!o.probeDirection.empty())
            metadata << "probe_direction=" << o.probeDirection << '\n';
        if (!o.replayCandidate.empty())
            metadata << "replay_candidate=" << o.replayCandidate << '\n';
        metadata << "pressure_probe=" << o.pressureProbe << "\nshape_curvature=" << o.shapeCurvature
                 << "\nsurface_trace=" << o.surfaceTrace << "\npositivity_certificate=" << o.positivityCertificate << '\n';
        metadata.flush();
        solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &f) {
            const auto snapshot = LoadVector(o.probeState, s.AcceptedPhysicalState().Size());
            if (o.positivityCertificate) {
                // A legacy endpoint may already violate admissibility. Its
                // polynomial must be inspectable without preparing that state.
                mfem::Vector normalizedSnapshot;
                s.normalizedOperator->NormalizeState(snapshot,normalizedSnapshot);
                experiment::plateau_diagnostics::CertifyPolynomialEnthalpy(s,normalizedSnapshot,o.output);
                return;
            }
            s.normalizedOperator->NormalizeState(snapshot, s.trialNormalizedState);
            auto prep = s.PrepareTrial();
            if (!prep)
                throw std::runtime_error("Snapshot is inadmissible");
            std::cout << "Snapshot F=" << s.trialNormalizedResidual.Norml2() << std::endl;
            if (!o.probeDirection.empty() || !o.replayCandidate.empty()) {
                FocusedProbes(s, f, o, snapshot);
                return;
            }
            experiment::plateau_diagnostics::InspectPreconditioner(s, o.output);
            FrozenProbes(s, f, o, &snapshot);
        });
        return 0;
    }
    auto physicalMetrics = [&](auto &s, const fem::FEM &f, const std::string &label) {
        std::cout << label << " physical diagnostics" << std::endl;
        experiment::polytrope_validation::PhysicalState physical(s, f);
        auto out = File(o.output / (label + "_physical.csv"));
        out << "quadrature,metric,value\n";
        if (o.surfacePrescription == "harmonics") {
            deformation::SphericalHarmonicSurfaceOptions surfaceOptions{.maximumDegree = o.surfaceDegree};
            if (o.surfaceMaximumOrder != -1) surfaceOptions.maximumAzimuthalOrder = o.surfaceMaximumOrder;
            experiment::plateau_diagnostics::RecordHarmonicSurface(s, f, surfaceOptions, o.output, label);
        }
        for (int q : {14, 18}) {
            auto metrics = experiment::plateau_diagnostics::MeasurePhysical(physical, o.n, scales.K, q);
            metrics["normalized_bordered_residual"] = physical.normalizedBorderedResidualNorm;
            metrics["normalized_unbordered_residual"] = physical.normalizedUnborderedResidualNorm;
            metrics["central_border"] = physical.centralBorder;
            metrics["normalized_border_action"] = physical.normalizedCentralBorderActionNorm;
            metrics["minimum_sampled_jacobian"] = s.acceptedMinimumJacobianDeterminant;
            metrics["mass_absolute_error"] = std::abs(metrics["mass"] - 1.);
            metrics["angular_momentum_absolute_error"] = std::abs(metrics["angular_momentum_z"] - o.angularMomentum);
            std::cout << "  q=" << q << " virial_h=" << metrics["enthalpy_virial_error"] << " M=" << metrics["mass"]
                      << std::endl;
            for (auto &[k, v] : metrics)
                out << q << ',' << k << ',' << v << '\n';
            out.flush();
        }
    };
    solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &f) {
        auto r = s.Problem().GetManifest().residualBlocks(), v = s.Problem().GetManifest().valueBlocks();
        rows.assign(r.begin(), r.end());
        values.assign(v.begin(), v.end());
        residualFactors = s.normalizedOperator->GetNormalization().ResidualFactors();
        stateFactors = s.normalizedOperator->GetNormalization().StateFactors();
        for (const auto &b : values)
            if (b.stableId == "fixed_central_density.border")
                centralOffset = b.offset;
        experiment::polytrope_validation::PhysicalState physical(s, f);
        centralRows = physical.centralHydrostaticRows;
        auto layout = File(o.output / "layout.csv");
        layout << "kind,block,offset,size,scale\n";
        for (auto &b : rows)
            layout << "residual," << b.stableId << ',' << b.offset << ',' << b.size << ',' << b.scale << '\n';
        for (auto &b : values)
            layout << "value," << b.stableId << ',' << b.offset << ',' << b.size << ',' << b.scale << '\n';
        physicalMetrics(s, f, "initial");
    });
    auto history = File(o.output / "newton.csv");
    history << "iteration,accepted,F,alpha,trials,linear_status,linear_iterations,linear_true,seconds,linear_seconds,"
               "prepare_seconds,refresh_seconds,geometry_limited,geometry_safe_step,central_border,unbordered_F,"
               "linear_initial_residual,linear_rhs_norm,linear_initial_over_rhs\n";
    auto blocklog = File(o.output / "blocks.csv");
    blocklog << "iteration,stage,block,normalized_norm,physical_norm,max_abs_normalized\n";
    auto triallog = File(o.output / "trials.csv");
    triallog << "iteration,trial,alpha,disposition,source,F,min_jacobian\n";
    auto coupledLog = File(o.output / "coupled_trials.csv");
    coupledLog << "iteration,trial,alpha,diagnostics_available,status,attempts_total,kept_total,fallbacks_total,"
                  "iterations_total,preparations_total,restorations_total,solve_seconds_total,verification_seconds_total,"
                  "preparation_seconds_total,refresh_seconds_total,original_merit,corrected_merit,has_solve,"
                  "linear_status,linear_iterations,rhs_norm,reported_true_norm,verified_true_norm,threshold,verified_gate,"
                  "correction_norm,maximum_correction_norm\n";
    auto coupledSearchLog = File(o.output / "coupled_search.csv");
    coupledSearchLog << "iteration,trial,alpha,status,deferred_total,selected_beta,candidate_count,base_merit,predicted_merit,predicted_decrease,excess_ratio\n";
    auto coupledCandidatesLog = File(o.output / "coupled_candidates.csv");
    coupledCandidatesLog << "iteration,trial,candidate,beta,status,merit,accepted,preparation_seconds\n";
    auto stateConsistency = File(o.output / "trial_state_consistency.csv");
    stateConsistency << "iteration,trial,normalized_physical_roundtrip_error,relative_error\n";
    mfem::Vector captureBase;
    auto captureDirectory = [&](int iteration) {
        return o.output / ("capture_" + std::to_string(iteration));
    };
    auto shouldCapture = [&](int iteration) {
        return std::find(o.captureIterations.begin(), o.captureIterations.end(), iteration) != o.captureIterations.end();
    };
    auto copySpan = [](std::span<const mfem::real_t> values) {
        mfem::Vector result(static_cast<int>(values.size()));
        for (int i = 0; i < result.Size(); ++i)
            result(i) = values[i];
        return result;
    };
    auto record = [&](int it, const char *stage, std::span<const mfem::real_t> residual) {
        for (const auto &b : rows) {
            long double ps = 0;
            double mx = 0;
            for (int i = b.offset; i < b.offset + b.size; ++i) {
                const double p = residual[i] / residualFactors(i);
                ps += (long double)p * p;
                mx = std::max(mx, std::abs(residual[i]));
            }
            blocklog << it << ',' << stage << ',' << b.stableId << ',' << Norm(residual, b.offset, b.size) << ','
                     << std::sqrt(ps) << ',' << mx << '\n';
        }
        blocklog.flush();
    };
    auto observer = solver::nonlinear::makeObserver(
        [&](const solver::nonlinear::BeforeIteration &e) {
            record(e.iteration, "before", e.normalizedResidual);
            if (shouldCapture(e.iteration)) {
                const auto directory = captureDirectory(e.iteration);
                std::filesystem::create_directory(directory);
                captureBase = copySpan(e.normalizedState);
                SaveVector(directory / "normalized_base.txt", captureBase);
                SaveVector(directory / "physical_base.txt", copySpan(e.physicalState));
                SaveVector(directory / "normalized_residual.txt", copySpan(e.normalizedResidual));
            }
            std::cout << "Newton " << e.iteration << " F=" << e.residualNorm << std::endl;
        },
        [&](const solver::nonlinear::AfterLineSearchTrial &e) {
            coupledLog << e.iteration << ',' << e.trial << ',' << e.stepLength << ',' << (e.numericalStrategies != nullptr);
            if (e.numericalStrategies) {
                const auto &c = e.numericalStrategies->coupledTrial;
                coupledSearchLog << e.iteration << ',' << e.trial << ',' << e.stepLength << ',' << int(c.status) << ','
                    << c.deferredAttempts << ',' << c.selectedBeta << ',' << c.candidateCount << ',' << c.baseMerit << ','
                    << c.linearPredictionMerit << ',' << c.predictedDecrease << ',' << c.predictionExcessRatio << '\n';
                int candidateIndex = 0;
                for (const auto &candidate : c.candidates)
                    coupledCandidatesLog << e.iteration << ',' << e.trial << ',' << candidateIndex++ << ','
                        << candidate.beta << ',' << int(candidate.status) << ',' << candidate.merit << ','
                        << candidate.accepted << ',' << candidate.preparationSeconds << '\n';
                coupledSearchLog.flush();
                coupledCandidatesLog.flush();
                coupledLog << ',' << static_cast<int>(c.status) << ',' << c.attempts << ',' << c.kept << ','
                    << c.fallbacks << ',' << c.iterations << ',' << c.preparations << ',' << c.restorations << ','
                    << c.solveSeconds << ',' << c.verificationSeconds << ',' << c.preparationSeconds << ','
                    << c.refreshSeconds << ',' << c.originalMerit << ',' << c.correctedMerit << ',' << c.currentSolve.has_value();
                if (c.currentSolve) {
                    const auto &r = *c.currentSolve;
                    coupledLog << ',' << static_cast<int>(r.linearSolve.status) << ',' << r.linearSolve.iterations << ','
                        << r.linearSolve.rightHandSideNorm << ',' << r.linearSolve.trueResidualNorm << ','
                        << r.verifiedTrueResidualNorm << ',' << r.convergenceThreshold << ',' << r.verifiedGate << ','
                        << r.correctionNorm << ',' << r.maximumCorrectionNorm;
                } else coupledLog << ",-1,0,nan,nan,nan,nan,0,nan,nan";
            } else coupledLog << ",-1,0,0,0,0,0,0,0,0,0,0,nan,nan,0,-1,0,nan,nan,nan,nan,0,nan,nan";
            coupledLog << '\n';
            coupledLog.flush();
            if (e.candidatePhysicalState.size() != e.candidateNormalizedState.size() ||
                e.candidateNormalizedState.size() != static_cast<size_t>(stateFactors.Size()))
                throw std::runtime_error("Trial observer exposed incompatible state representations");
            long double stateError = 0.0L, stateNorm = 0.0L;
            for (int i = 0; i < stateFactors.Size(); ++i) {
                const long double difference = e.candidatePhysicalState[i] * stateFactors(i) - e.candidateNormalizedState[i];
                stateError += difference*difference;
                stateNorm += static_cast<long double>(e.candidateNormalizedState[i])*e.candidateNormalizedState[i];
            }
            const double relativeStateError = std::sqrt(stateError)/std::max(1.0L,std::sqrt(stateNorm));
            stateConsistency << e.iteration << ',' << e.trial << ',' << std::sqrt(stateError) << ',' << relativeStateError << '\n';
            stateConsistency.flush();
            if (std::isfinite(relativeStateError) && relativeStateError > 64.0*std::numeric_limits<double>::epsilon())
                throw std::runtime_error("Trial observer physical state disagrees with its normalized state");
            triallog << e.iteration << ',' << e.trial << ',' << e.stepLength << ',' << int(e.disposition) << ','
                     << e.rejectionSource << ',' << (e.metric ? e.metric->residualNorm : NAN) << ','
                     << e.minimumJacobianDeterminant.value_or(NAN) << '\n';
            triallog.flush();
            if (e.trial == 0 && shouldCapture(e.iteration)) {
                if (!(e.stepLength > 0.) || !std::isfinite(e.stepLength) ||
                    e.candidateNormalizedState.size() != static_cast<size_t>(captureBase.Size()))
                    throw std::runtime_error("Invalid candidate for diagnostic capture");
                const auto directory = captureDirectory(e.iteration);
                const auto candidate = copySpan(e.candidateNormalizedState);
                SaveVector(directory / "normalized_candidate.txt", candidate);
                SaveVector(directory / "physical_candidate.txt", copySpan(e.candidatePhysicalState));
                mfem::Vector direction(candidate), reconstructed(captureBase);
                direction -= captureBase;
                direction /= e.stepLength;
                reconstructed.Add(e.stepLength, direction);
                reconstructed -= candidate;
                SaveVector(directory / "normalized_effective_trial_direction.txt", direction);
                SaveVector(directory / "normalized_direction.txt", copySpan(e.normalizedDirection));
                auto info = File(directory / "capture.txt");
                info << "iteration=" << e.iteration << "\nfirst_trial_alpha=" << e.stepLength
                     << "\ndisposition=" << int(e.disposition) << "\nsource=" << e.rejectionSource
                     << "\nbase_norm=" << captureBase.Norml2() << "\ndirection_norm=" << direction.Norml2()
                     << "\ncandidate_reconstruction_error=" << reconstructed.Norml2()
                     << "\ndirection_roundoff_scale_estimate="
                     << 8. * std::numeric_limits<double>::epsilon() *
                            (captureBase.Norml2() + candidate.Norml2()) / e.stepLength
                     << "\ndirection_source=exact_observer_tangent\n"
                     << "effective_trial_direction_source=candidate_difference_divided_by_alpha\n"
                     << "roundoff_scale_applies_to_effective_trial_direction_only=1\n"
                     << "roundoff_scale_is_estimate_not_certificate=1\n";
            }
        },
        [&](const solver::nonlinear::AfterIteration &e) {
            record(e.iteration, "after", e.normalizedResidual);
            double border = centralOffset < 0 ? 0 : e.physicalState[centralOffset];
            mfem::Vector unbordered((int)e.normalizedResidual.size());
            for (int i = 0; i < unbordered.Size(); ++i)
                unbordered(i) = e.normalizedResidual[i];
            for (int i : centralRows)
                unbordered(i) -= residualFactors(i) * border;
            const double linearInitial = e.linearSolve ? e.linearSolve->initialResidualNorm : NAN;
            const double linearRhs = e.linearSolve ? e.linearSolve->rightHandSideNorm : NAN;
            double linearInitialRelative = NAN;
            if (linearRhs > 0.)
                linearInitialRelative = linearInitial / linearRhs;
            else if (linearRhs == 0. && std::isfinite(linearInitial))
                linearInitialRelative = linearInitial == 0. ? 0. : std::numeric_limits<double>::infinity();
            history << e.iteration << ',' << e.stepAccepted << ',' << e.residualNorm << ',' << e.acceptedStepLength
                    << ',' << e.lineSearchTrials << ',' << (e.linearSolve ? int(e.linearSolve->status) : -1) << ','
                    << (e.linearSolve ? e.linearSolve->iterations : 0) << ','
                    << (e.linearSolve ? e.linearSolve->relativeTrueResidualNorm : NAN) << ',' << e.iterationSeconds
                    << ',' << (e.linearSolve ? e.linearSolve->solveSeconds : 0) << ',' << e.trialPreparationSeconds
                    << ',' << e.preconditionerRefreshSeconds << ','
                    << (e.geometryPreflight ? e.geometryPreflight->limitedByGeometry : false) << ','
                    << (e.geometryPreflight ? e.geometryPreflight->stepSize : NAN) << ',' << border << ','
                    << unbordered.Norml2() << ',' << linearInitial << ',' << linearRhs << ','
                    << linearInitialRelative << '\n';
            history.flush();
            std::cout << "  alpha=" << e.acceptedStepLength << " F=" << e.residualNorm
                      << " linear=" << (e.linearSolve ? e.linearSolve->relativeTrueResidualNorm : NAN)
                      << " lambda=" << border << " unbordered=" << unbordered.Norml2()
                      << " seconds=" << e.iterationSeconds << std::endl;
        });
    bool converged = false;
    if (o.maxNewton > 0) {
        namespace nl = solver::nonlinear;
        const nl::NewtonOptions options{
            .relativeTolerance = o.rtol, .absoluteTolerance = o.atol, .maximumIterations = o.maxNewton,
            .linearSolve = {.relativeTolerance = o.linearTolerance, .absoluteTolerance = 0., .maximumIterations = o.maxLinear}};
        const nl::GravitySolveOptions gravity{
            .control = {.relativeTolerance = o.gravityTolerance, .absoluteTolerance = 1e-14,
                        .maximumIterations = o.gravityIterations},
            .backend = {.restartLength = o.gravityRestart, .printLevel = -1}};
        const nl::CoupledCurvatureCorrection coupled{
            .correction = {.control = {.relativeTolerance = o.coupledTolerance, .absoluteTolerance = 1e-14,
                                       .maximumIterations = o.coupledIterations},
                           .backend = {.restartLength = o.coupledRestart, .printLevel = -1},
                           .maximumCorrectionRatio = o.coupledMaximumRatio},
            .maximumAttemptsPerIteration = o.coupledAttempts};
        auto run = [&](auto numerics) {
            auto solve = solver::make(context, nl::Newton(options, nl::NormalizedResidualMetric{}, numerics), observer);
            return solve.evaluate();
        };
        auto report = [&] {
            if (o.strategy == "positivity") return run(nl::NewtonNumerics{nl::PolynomialEnthalpy{}, nl::MonolithicGravity{}});
            if (o.strategy == "curvature") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::GravityCurvatureCorrection{gravity}});
            const nl::BacktrackedCoupledCurvatureCorrection backtracked{
                .correction = coupled.correction, .maximumAttemptsPerIteration = o.coupledAttempts,
                .search = {.maximumSteps = o.coupledBetaSteps, .reduction = o.coupledBetaReduction},
                .attempt = {.maximumExcessRatio = o.coupledExcessRatio}};
            if (o.strategy == "backtracked-coupled") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::MonolithicGravity{}, backtracked});
            if (o.strategy == "curvature-backtracked-coupled") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::GravityCurvatureCorrection{gravity}, backtracked});
            if (o.strategy == "coupled") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::MonolithicGravity{}, coupled});
            if (o.strategy == "curvature-coupled") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::GravityCurvatureCorrection{gravity}, coupled});
            if constexpr (FullDefaultStack)
                if (o.strategy == "reduced") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::EliminatedGravity{gravity}});
            if (o.strategy == "positivity-curvature") return run(nl::NewtonNumerics{nl::PolynomialEnthalpy{}, nl::GravityCurvatureCorrection{gravity}});
            if constexpr (FullDefaultStack)
                if (o.strategy == "positivity-reduced") return run(nl::NewtonNumerics{nl::PolynomialEnthalpy{}, nl::EliminatedGravity{gravity}});
            return run(nl::NewtonNumerics{});
        }();
        const auto &strategyDiagnostics = report.diagnostics().numericalStrategies;
        auto strategyMetrics = File(o.output / "numerical_strategies.csv");
        strategyMetrics << "metric,value\n";
        const auto metric = [&](const char *name, auto value) { strategyMetrics << name << ',' << value << '\n'; };
        metric("polynomial_checks", strategyDiagnostics.polynomialChecks);
        metric("direction_repairs", strategyDiagnostics.directionRepairs);
        metric("last_repair_fraction", strategyDiagnostics.lastRepairFraction);
        metric("last_merit_directional_derivative", strategyDiagnostics.lastMeritDirectionalDerivative);
        metric("gravity_solves", strategyDiagnostics.gravitySolves);
        metric("gravity_iterations", strategyDiagnostics.gravityIterations);
        metric("gravity_failures", strategyDiagnostics.gravityFailures);
        metric("gravity_prepared_failures", strategyDiagnostics.gravityPreparedFailures);
        metric("compact_gravity_operator", strategyDiagnostics.compactGravityOperator);
        metric("compact_gravity_inverse", strategyDiagnostics.compactGravityInverse);
        metric("linear_report_space", static_cast<int>(strategyDiagnostics.lastLinearReportSpace));
        metric("gravity_solve_seconds", strategyDiagnostics.gravitySolveSeconds);
        metric("last_gravity_true_residual", strategyDiagnostics.lastGravityTrueResidual);
        metric("last_gravity_correction_norm", strategyDiagnostics.lastGravityCorrectionNorm);
        metric("last_gravity_prepared_defect", strategyDiagnostics.lastGravityPreparedDefect);
        metric("last_gravity_affine_disagreement", strategyDiagnostics.lastGravityAffineDisagreement);
        metric("last_gravity_arithmetic_allowance", strategyDiagnostics.lastGravityArithmeticAllowance);
        metric("last_gravity_prepared_target_threshold", strategyDiagnostics.lastGravityPreparedTargetThreshold);
        metric("initial_gravity_residual", strategyDiagnostics.initialGravityResidual);
        metric("completed_initial_gravity_residual", strategyDiagnostics.completedInitialGravityResidual);
        metric("last_lifted_relative_true_residual", strategyDiagnostics.lastLiftedRelativeTrueResidual);
        const auto &coupledDiagnostics = strategyDiagnostics.coupledTrial;
        metric("coupled_attempts", coupledDiagnostics.attempts);
        metric("coupled_deferred_attempts", coupledDiagnostics.deferredAttempts);
        metric("coupled_kept", coupledDiagnostics.kept);
        metric("coupled_fallbacks", coupledDiagnostics.fallbacks);
        metric("coupled_iterations", coupledDiagnostics.iterations);
        metric("coupled_preparations", coupledDiagnostics.preparations);
        metric("coupled_restorations", coupledDiagnostics.restorations);
        metric("coupled_solve_seconds", coupledDiagnostics.solveSeconds);
        metric("coupled_verification_seconds", coupledDiagnostics.verificationSeconds);
        metric("coupled_preparation_seconds", coupledDiagnostics.preparationSeconds);
        metric("coupled_refresh_seconds", coupledDiagnostics.refreshSeconds);
        converged = report.converged();
        metadata << "converged=" << converged << "\naccepted_steps=" << report.completedNonlinearIterations() << '\n';
        if (!converged)
            metadata << "failure=" << report.failure().message << '\n';
        metadata.flush();
    }
    solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &f) {
        SaveVector(o.output / "accepted_state.txt", s.AcceptedPhysicalState());
        physicalMetrics(s, f, "final");
        if (o.densityOrderFinal) {
            experiment::polytrope_validation::PhysicalState physical(s, f);
            experiment::plateau_diagnostics::RecordDensityOrderState(physical, o.n, scales.K, o.output);
            auto rules = File(o.output / "final_closure_rules.csv");
            rules << "element,quadrature_order,quadrature_points\n";
            s.Problem().GetPhysicalOperator().GetBarotropicClosureOperator().VisitMappedGeometryRules(
                [&](int element, const mfem::IntegrationRule &rule) {
                    rules << element << ',' << rule.GetOrder() << ',' << rule.GetNPoints() << '\n';
                });
        }
        const auto &schur = s.physicalInverse->GetGroupedPreconditioner().GetFactorization().GetSchurComplement();
        auto sh = File(o.output / "border_schur.csv");
        for (int i = 0; i < schur.Height(); ++i) {
            for (int j = 0; j < schur.Width(); ++j)
                sh << (j ? "," : "") << schur(i, j);
            sh << '\n';
        }
        if constexpr (FullDefaultStack) if (o.probes) {
            experiment::plateau_diagnostics::InspectPreconditioner(s, o.output);
            FrozenProbes(s, f, o);
        }
        if (o.weakShapeFinal)
            experiment::plateau_diagnostics::ProbeWeakShapeField(s, f, s.acceptedNormalizedState, nullptr, o.output, {
                .linear = {.relativeTolerance = o.linearTolerance, .absoluteTolerance = 0.,
                           .maximumIterations = o.maxLinear}, .restartLength = o.restart});
    });
    return converged ? 0 : 3;
}
template <class Model, class Surface>
int SetupWithSurface(Model model, fem::FEM f, Surface surface, const Options &o, const Scales &scales) {
    auto d = equilibrium::makeStellarDiscretization(
        std::move(f), normalization::PhysicalRieszDiagonal{dimensions::LengthValue{1.}, utils::G},
        std::move(surface));
    const auto withRecipe = [&](auto recipe) {
        using Recipe = decltype(recipe);
        using Problem = equilibrium::StellarEquilibriumProblem<Model, decltype(d)>;
        static_assert(preconditioning::StellarPreconditionerRuntimeAvailableFor<Recipe, Problem>);
        using Inverse = preconditioning::PreparedStellarInverseType<Recipe, Problem>;
        static_assert(solver::LinearPreconditionerApplicationContractAvailable<Inverse>);
        static_assert(solver::linearPreconditionerApplicationContract<Inverse> ==
                      preconditioning::ApplicationContract::stationary_linear);
        auto context = solver::makeContext(std::move(model), std::move(d), std::move(recipe),
                                           solver::linear::FGMRES({.restartLength = o.restart, .printLevel = -1}),
                                           seed::LaneEmden({.centralDensity = dimensions::DensityValue{scales.rho}}));
        constexpr bool fullDefaultStack = std::same_as<Recipe, preconditioning::DefaultStellarPreconditioner>;
        return RunContext<fullDefaultStack>(context, o, scales);
    };
    if (o.preconditionerProbe || o.preconditionerVariant == "independent")
        return withRecipe(preconditioning::makePreconditioner());
    if (o.preconditionerVariant == "gravity-schur")
        return withRecipe(experiment::plateau_diagnostics::GravitySchurNewtonPreconditioner{
            {.correctionSweeps = o.schurCorrectionSweeps, .relaxation = o.schurRelaxation}});
    if (o.preconditionerVariant == "stellar-ldu")
        return withRecipe(experiment::plateau_diagnostics::StellarLDUNewtonPreconditioner{
            {.correctionSweeps = 0, .relaxation = 1.0}});
    if (o.preconditionerVariant == "stellar-ldu-gravity-schur")
        return withRecipe(experiment::plateau_diagnostics::StellarLDUGravitySchurNewtonPreconditioner{
            {.correctionSweeps = o.schurCorrectionSweeps, .relaxation = o.schurRelaxation}});
    throw std::logic_error("An unsupported Newton preconditioner bypassed option validation.");
}
template <class Model> int Setup(Model model, fem::FEM f, const Options &o, const Scales &scales) {
    if (o.surfacePrescription == "harmonics") {
        deformation::SphericalHarmonicSurfaceOptions options{.maximumDegree = o.surfaceDegree};
        if (o.surfaceMaximumOrder != -1) options.maximumAzimuthalOrder = o.surfaceMaximumOrder;
        return SetupWithSurface(std::move(model), std::move(f),
                                deformation::SphericalHarmonicSurface{options}, o, scales);
    }
    mfem::Vector center(3);
    center = 0.0;
    return SetupWithSurface(std::move(model), std::move(f),
                            deformation::NodalRadialSurface{std::move(center)}, o, scales);
}
int Run(const Options &o) {
    const auto scales = Derive(o.n);
    std::cout << "n=" << o.n << " K=" << scales.K << " rho_c=" << scales.rho << " J=" << o.angularMomentum << std::endl;
    utils::Args args;
    args.mesh_file = o.mesh;
    args.p.rtol = 1e-12;
    args.p.atol = 1e-12;
    auto f = fem::setup_fem(o.mesh, args, 0);
    if (!f.okay())
        throw std::runtime_error("FEM setup failed");
    auto eos = eos::Polytrope({.n = o.n, .K = scales.K});
    auto surf = surface::Isobaric({.Psurf = dimensions::PressureValue{0.}});
    auto mass = integral::FixedTotalMass({.Mtotal = dimensions::MassValue{1.}});
    auto angular = integral::FixedAngularMomentum(
        {.Jtotal = dimensions::AngularMomentumValue{o.angularMomentum}, .axis = {0., 0., 1.}, .center = {0., 0., 0.}});
    if (o.noCentral)
        return Setup(model::StellarModel(eos, surf, mass, angular), std::move(f), o, scales);
    return Setup(model::StellarModel(eos, surf, mass, angular,
                                     constraint::FixedCentralDensity({.RhoC = dimensions::DensityValue{scales.rho}})),
                 std::move(f), o, scales);
}
} // namespace
int main(int argc, char **argv) {
    mfem::Mpi::Init(argc, argv);
    int result = 0;
    try {
        int count;
        MPI_Comm_size(MPI_COMM_WORLD, &count);
        if (count != 1)
            throw std::runtime_error("Single rank only");
        mfem::Device device("cpu");
        std::cout << std::setprecision(12);
        const auto options = Parse(argc, argv);
        if (options.reducedSurfaceAlgebraCheck) {
            const auto checks = experiment::plateau_diagnostics::RunReducedSurfaceAlgebraSelfCheck(MPI_COMM_WORLD);
            std::cout << "Reduced surface algebra self-check passed: " << checks.checks
                      << " checks, maximum error " << checks.maximumError << std::endl;
        } else {
            result = Run(options);
        }
    } catch (const std::exception &e) {
        std::cerr << "Diagnostic failure: " << e.what() << std::endl;
        result = 2;
    }
    mfem::Mpi::Finalize();
    return result;
}