#pragma once

// Experiment-only: include after import mean_field. The caller holds the state
// lease. A physical snapshot uses the context's frozen normalization throughout.
#include <algorithm>
#include <array>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <cstring>
#include <exception>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <limits>
#include <stdexcept>
#include <string>
#include <string_view>
#include <type_traits>

namespace experiment::plateau_diagnostics {
struct NonlinearStepProbeOptions final {
    mean_field::solver::LinearSolveControl linear{.relativeTolerance = .03, .absoluteTolerance = 0.,
                                                 .maximumIterations = 200};
    int restartLength{40};
    mean_field::solver::nonlinear::GravitySolveOptions gravity{};
};

namespace nonlinear_step_detail {
using Clock = std::chrono::steady_clock;
inline double Seconds(Clock::time_point start) {
    return std::chrono::duration<double>(Clock::now() - start).count();
}
inline bool Finite(const mfem::Vector &v) {
    for (int i = 0; i < v.Size(); ++i) if (!std::isfinite(v(i))) return false;
    return true;
}
inline bool Identical(const mfem::Vector &a, const mfem::Vector &b) {
    return a.Size() == b.Size() &&
           (a.Size() == 0 || std::memcmp(a.GetData(), b.GetData(), sizeof(mfem::real_t) * a.Size()) == 0);
}
inline long double Dot(const mfem::Vector &a, const mfem::Vector &b, int offset, int size) {
    long double result = 0.;
    for (int i = offset; i < offset + size; ++i) result += static_cast<long double>(a(i)) * b(i);
    return result;
}
inline double Norm(const mfem::Vector &v, int offset, int size) {
    return std::sqrt(static_cast<double>(Dot(v, v, offset, size)));
}
inline double Norm(const mfem::Vector &v) { return Norm(v, 0, v.Size()); }
inline double Difference(const mfem::Vector &a, const mfem::Vector &b) {
    mfem::Vector difference(a);
    difference -= b;
    return Norm(difference);
}
inline double Ratio(double numerator, double denominator) {
    return std::isfinite(numerator) && std::isfinite(denominator) && denominator > 0.
               ? numerator / denominator : std::numeric_limits<double>::quiet_NaN();
}
inline std::ofstream Open(const std::filesystem::path &path) {
    if (std::filesystem::exists(path))
        throw std::invalid_argument("Nonlinear-step output already exists: " + path.string());
    std::ofstream file(path);
    file.exceptions(std::ios::failbit | std::ios::badbit);
    file << std::setprecision(17);
    return file;
}
inline void Save(const std::filesystem::path &path, const mfem::Vector &v) {
    auto file = Open(path);
    file << v.Size() << '\n';
    for (int i = 0; i < v.Size(); ++i) file << v(i) << '\n';
}
struct Observation final {
    bool prepared{false}, innerConverged{true}, preparedGate{true}, usable{false};
    mfem::Vector residual, remainder;
    double correctionNorm{std::numeric_limits<double>::quiet_NaN()};
    double targetDefect{std::numeric_limits<double>::quiet_NaN()};
    double affineDisagreement{std::numeric_limits<double>::quiet_NaN()};
    double arithmeticAllowance{std::numeric_limits<double>::quiet_NaN()};
    double targetThreshold{std::numeric_limits<double>::quiet_NaN()};
    double seconds{0.};
};
} // namespace nonlinear_step_detail

// The SOC target is G0 + alpha J_G p, not zero gravity residual. Its retained
// rows may have substantial curvature even when that target is met accurately.
// Full and split vectors let independent analysis test this trial map's tangent
// without mistaking scalar norm cancellation for a correct vector derivative.
template <class State>
void ProbeNonlinearStep(State &state, const mfem::Vector &normalizedBase,
                        const std::filesystem::path &output, NonlinearStepProbeOptions options = {}) {
    using namespace mean_field;
    namespace d = nonlinear_step_detail;
    using Form = typename std::remove_cvref_t<decltype(state.Problem())>::FormType;
    using namespace utils::blocks;
    constexpr value_block<type_index_v<surface_deformation::parameters::value, typename Form::value_blocks>> shape{};
    constexpr value_block<type_index_v<gravity::gradient::value, typename Form::value_blocks>> gradient{};
    constexpr value_block<type_index_v<gravity::poisson::value, typename Form::value_blocks>> potential{};
    options.linear.Validate();
    options.gravity.Validate();
    int ranks = 0;
    MPI_Comm_size(state.Problem().GetCommunicator(), &ranks);
    const int size = state.Problem().StateSize();
    if (ranks != 1 || size != state.Problem().EquationSize() || normalizedBase.Size() != size ||
        !d::Finite(normalizedBase) || options.restartLength < 1)
        throw std::invalid_argument("Nonlinear-step probes require one rank, compatible finite vectors and valid controls.");
    const auto start = d::Clock::now();
    const double nan = std::numeric_limits<double>::quiet_NaN();
    const double eps = std::numeric_limits<mfem::real_t>::epsilon();
    const mfem::Vector base(normalizedBase);
    const auto &normalization = state.normalizedOperator->GetNormalization();
    const mfem::Vector stateFactors(normalization.StateFactors()), residualFactors(normalization.ResidualFactors());
    const mfem::Vector savedAcceptedState(state.acceptedNormalizedState),
        savedAcceptedResidual(state.acceptedNormalizedResidual), savedAcceptedPhysical(state.AcceptedPhysicalState());
    const mfem::Vector savedCorrection(state.normalizedCorrection), savedPhysicalCorrection(state.physicalCorrection),
        savedVolumeDirection(state.volumeDisplacementDirection);
    const auto &layout = state.Problem().GetManifest().layout();
    const auto rows = state.Problem().GetManifest().residualBlocks();
    const int shapeOffset = layout.offset(shape), shapeSize = layout.size(shape);
    if (shapeSize <= 0) throw std::invalid_argument("The nonlinear-step probe requires a nonempty shape block.");
    constexpr std::array<const char *, 3> names{"full", "shape", "fields"};
    constexpr std::array<double, 8> alphas{1., .25, .0625, .015625, .00390625, .0009765625,
                                           .000244140625, .00006103515625};
    std::filesystem::create_directories(output);
    const auto file = [&](const std::string &name) { return output / ("nonlinear_step_" + name); };
    auto summary = d::Open(file("summary.csv"));
    auto linear = d::Open(file("linear.csv"));
    auto preparations = d::Open(file("preparations.csv"));
    auto paths = d::Open(file("paths.csv"));
    auto blocks = d::Open(file("blocks.csv"));
    auto mixed = d::Open(file("mixed.csv"));
    auto inner = d::Open(file("inner.csv"));
    auto restorations = d::Open(file("restorations.csv"));
    summary << "metric,value\n";
    linear << "status,converged,iterations,rhs_norm,reported_true_norm,verified_true_norm,verified_relative_true,"
              "threshold,verified_gate,operator_applications,inverse_applications,seconds\n";
    preparations << "sequence,kind,direction,alpha,prepared,rejection_reason,rejection_stage,seconds,residual_norm\n";
    paths << "sample,alpha,direction,path,prepared,inner_converged,prepared_gate,usable,base_norm,action_norm,"
             "residual_norm,merit,slope,linear_prediction_norm,remainder_norm,remainder_over_alpha,"
             "remainder_over_alpha_squared,relative_derivative_error,linear_dot_remainder,norm_identity_error,"
             "expected_step_norm,actual_step_norm,addition_error_norm,increment_resolved,gravity_correction_norm,"
             "gravity_target_defect,gravity_affine_disagreement,gravity_arithmetic_allowance,gravity_target_threshold,"
             "seconds,heuristic_roundoff_scale\n";
    blocks << "sample,alpha,direction,path,block,offset,size,base_norm,action_norm,linear_prediction_norm,residual_norm,"
              "remainder_norm,remainder_over_alpha,remainder_over_alpha_squared,linear_dot_remainder,"
              "norm_identity_error,heuristic_roundoff_scale\n";
    mixed << "sample,alpha,path,block,offset,size,usable,full_remainder_norm,shape_remainder_norm,fields_remainder_norm,"
             "mixed_remainder_norm,linear_dot_full,linear_dot_shape,linear_dot_fields,linear_dot_mixed,"
             "projection_identity_error,norm_identity_error\n";
    inner << "sample,alpha,direction,status,converged,iterations,rhs_norm,true_norm,relative_true,threshold,"
             "operator_applications,inverse_applications,seconds,compact_operator,compact_inverse\n";
    restorations << "sequence,kind,block,offset,size,base_norm,difference_norm\n";
    auto metric = [&](const char *name, auto value) { summary << name << ',' << value << '\n'; };
    metric("schema_version", 1);
    metric("linear_relative_tolerance", options.linear.relativeTolerance);
    metric("linear_absolute_tolerance", options.linear.absoluteTolerance);
    metric("linear_maximum_iterations", options.linear.maximumIterations);
    metric("linear_restart_length", options.restartLength);
    metric("gravity_relative_tolerance", options.gravity.control.relativeTolerance);
    metric("gravity_absolute_tolerance", options.gravity.control.absoluteTolerance);
    metric("gravity_maximum_iterations", options.gravity.control.maximumIterations);
    metric("gravity_restart_length", options.gravity.backend.restartLength);
    metric("shape_offset", shapeOffset);
    metric("shape_size", shapeSize);
    metric("machine_epsilon", eps);
    metric("polynomial_direction_repairs", 0);
    metric("accepted_commits", 0);
    metric("linear_gate_failure_skips_trials", 0);
    metric("path_seconds_include_shared_raw_preparation", 1);
    metric("roundoff_scale_is_subtraction_heuristic", 1);
    for (std::size_t sample = 0; sample < alphas.size(); ++sample)
        summary << "alpha_" << sample << ',' << alphas[sample] << '\n';
    mfem::Vector originalResidual(size);
    bool haveResidual = false;
    int preparationCount = 0, refreshCount = 0, directFullActions = 0, captureFullActions = 0;
    int assessmentGravityActions = 0, successfulPreparedPaths = 0;
    double preparationSeconds = 0., refreshSeconds = 0., directActionSeconds = 0., captureSeconds = 0.;
    double assessmentGravitySeconds = 0.;
    double restorationError = 0.;
    auto eachBlock = [&](auto &&operation) {
        operation("all", 0, size);
        for (const auto &row : rows) operation(row.stableId, row.offset, row.size);
    };
    auto prepare = [&](const mfem::Vector &candidate, const char *kind, const char *direction, double alpha) {
        state.trialNormalizedState = candidate;
        const auto preparationStart = d::Clock::now();
        const auto result = state.PrepareTrial();
        const double elapsed = d::Seconds(preparationStart);
        preparationSeconds += elapsed;
        ++preparationCount;
        preparations << preparationCount << ',' << kind << ',' << direction << ',' << alpha << ',' << bool(result)
                     << ',' << (result ? -1 : int(result.error().reason)) << ','
                     << (result ? -1 : int(result.error().stage)) << ',' << elapsed << ','
                     << (result ? d::Norm(state.trialNormalizedResidual) : nan) << '\n';
        preparations.flush();
        if (result && !d::Finite(state.trialNormalizedResidual))
            throw std::runtime_error("Nonlinear-step preparation produced a nonfinite residual.");
        return bool(result);
    };
    auto refresh = [&]() {
        const auto refreshStart = d::Clock::now();
        state.physicalInverse->Refresh();
        refreshSeconds += d::Seconds(refreshStart);
        ++refreshCount;
    };
    auto restore = [&](const char *kind) {
        if (!prepare(base, kind, "base", 0.))
            throw std::runtime_error("The nonlinear-step probe could not restore its base.");
        if (!haveResidual) { originalResidual = state.trialNormalizedResidual; haveResidual = true; }
        mfem::Vector difference(state.trialNormalizedResidual);
        difference -= originalResidual;
        restorationError = std::max(restorationError, d::Norm(difference));
        eachBlock([&](auto name, int offset, int count) {
            restorations << preparationCount << ',' << kind << ',' << name << ',' << offset << ',' << count << ','
                         << d::Norm(originalResidual, offset, count) << ',' << d::Norm(difference, offset, count) << '\n';
        });
        restorations.flush();
    };
    auto action = [&](const mfem::Vector &direction, mfem::Vector &image) {
        const auto actionStart = d::Clock::now();
        state.normalizedOperator->Mult(direction, image);
        directActionSeconds += d::Seconds(actionStart);
        ++directFullActions;
        if (!d::Finite(image)) throw std::runtime_error("Nonlinear-step J action is nonfinite.");
    };
    auto finish = [&]() {
        restore("final_restore");
        d::Save(file("restored_residual.txt"), state.trialNormalizedResidual);
        refresh();
        state.normalizedCorrection = savedCorrection;
        state.physicalCorrection = savedPhysicalCorrection;
        state.volumeDisplacementDirection = savedVolumeDirection;
        const bool acceptedUnchanged = d::Identical(savedAcceptedState, state.acceptedNormalizedState) &&
            d::Identical(savedAcceptedResidual, state.acceptedNormalizedResidual) &&
            d::Identical(savedAcceptedPhysical, state.AcceptedPhysicalState());
        const bool factorsUnchanged = d::Identical(stateFactors, normalization.StateFactors()) &&
            d::Identical(residualFactors, normalization.ResidualFactors());
        metric("maximum_original_restoration_error", restorationError);
        metric("accepted_state_change", d::Difference(savedAcceptedState, state.acceptedNormalizedState));
        metric("accepted_residual_change", d::Difference(savedAcceptedResidual, state.acceptedNormalizedResidual));
        metric("accepted_physical_change", d::Difference(savedAcceptedPhysical, state.AcceptedPhysicalState()));
        metric("accepted_buffers_bitwise_identical", acceptedUnchanged);
        metric("normalization_factors_bitwise_identical", factorsUnchanged);
        metric("preparations", preparationCount);
        metric("preparation_seconds", preparationSeconds);
        metric("inverse_refreshes", refreshCount);
        metric("inverse_refresh_seconds", refreshSeconds);
        metric("direct_full_operator_applications", directFullActions);
        metric("direct_full_operator_seconds", directActionSeconds);
        metric("capture_full_operator_applications", captureFullActions);
        metric("capture_seconds", captureSeconds);
        metric("assessment_gravity_operator_applications", assessmentGravityActions);
        metric("assessment_gravity_operator_seconds", assessmentGravitySeconds);
        metric("prepared_paths", successfulPreparedPaths);
        metric("total_seconds", d::Seconds(start));
        summary.flush();
        if (!acceptedUnchanged || !factorsUnchanged)
            throw std::runtime_error("The nonlinear-step probe changed accepted buffers or normalization.");
    };

    try {
        restore("initial_restore");
        refresh();
        mfem::Vector physicalBase;
        state.normalizedOperator->DenormalizeState(base, physicalBase);
        d::Save(file("physical_base.txt"), physicalBase);
        d::Save(file("normalized_base.txt"), base);
        d::Save(file("residual.txt"), originalResidual);
        d::Save(file("state_factors.txt"), stateFactors);
        d::Save(file("residual_factors.txt"), residualFactors);
        auto manifest = d::Open(file("layout.csv"));
        manifest << "kind,block,offset,size\n";
        for (const auto &value : state.Problem().GetManifest().valueBlocks())
            manifest << "value," << value.stableId << ',' << value.offset << ',' << value.size << '\n';
        for (const auto &row : rows)
            manifest << "residual," << row.stableId << ',' << row.offset << ',' << row.size << '\n';
        mfem::Vector rhs(originalResidual);
        rhs *= -1.;
        d::Save(file("rhs.txt"), rhs);
        std::array<mfem::Vector, 3> directions{mfem::Vector(size), mfem::Vector(size), mfem::Vector(size)};
        std::array<mfem::Vector, 3> actions{mfem::Vector(size), mfem::Vector(size), mfem::Vector(size)};
        directions[0] = 0.;
        auto backend = solver::linear::prepareLinearBackend(
            solver::linear::FGMRES({.restartLength = options.restartLength, .printLevel = -1}),
            *state.normalizedOperator, *state.normalizedInverse, state.Problem().GetCommunicator());
        const auto solve = backend.Solve(rhs, directions[0], options.linear);
        // Reprepare before the independent certificate: this also tests that a
        // solve has not changed the operator's promised frozen linearization.
        restore("verify_linear_base");
        action(directions[0], actions[0]);
        mfem::Vector defect(actions[0]);
        defect -= rhs;
        const double trueNorm = d::Norm(defect), rhsNorm = d::Norm(rhs);
        const double threshold = options.linear.ConvergenceThreshold(rhsNorm);
        const bool trueGate = d::Finite(directions[0]) && std::isfinite(trueNorm) && trueNorm <= threshold;
        linear << int(solve.status) << ',' << solve.Converged() << ',' << solve.iterations << ',' << rhsNorm << ','
               << solve.trueResidualNorm << ',' << trueNorm << ',' << d::Ratio(trueNorm, rhsNorm) << ',' << threshold
               << ',' << trueGate << ',' << solve.operatorApplications << ','
               << solve.inversePreconditionerApplications << ',' << solve.solveSeconds << '\n';
        linear.flush();
        metric("linear_verified_gate", trueGate);
        metric("linear_solve_operator_applications", solve.operatorApplications);
        metric("linear_solve_inverse_applications", solve.inversePreconditionerApplications);
        metric("linear_solve_seconds", solve.solveSeconds);
        metric("linear_solve_operator_seconds", solve.operatorSeconds);
        metric("linear_solve_inverse_seconds", solve.inversePreconditionerSeconds);
        if (!d::Finite(directions[0])) throw std::runtime_error("A frozen direction is nonfinite.");
        directions[1] = 0.;
        directions[2] = directions[0];
        for (int i = shapeOffset; i < shapeOffset + shapeSize; ++i) {
            directions[1](i) = directions[0](i);
            directions[2](i) = 0.;
        }
        for (int k = 1; k < 3; ++k) { restore("split_action_base"); action(directions[k], actions[k]); }
        mfem::Vector splitError(actions[0]);
        splitError -= actions[1];
        splitError -= actions[2];
        metric("action_split_error_norm", d::Norm(splitError));
        metric("base_residual_norm", d::Norm(originalResidual));
        metric("full_direction_norm", d::Norm(directions[0]));
        metric("full_merit_slope", d::Dot(originalResidual, actions[0], 0, size));
        for (int k = 0; k < 3; ++k) {
            d::Save(file(std::string("direction_") + names[k] + ".txt"), directions[k]);
            d::Save(file(std::string("action_") + names[k] + ".txt"), actions[k]);
        }
        restore("gravity_setup_base");
        refresh();
        auto partition = solver::MakeGravityBlockPartition<Form>(layout);
        solver::PreparedGravitySubsystem gravitySubsystem(*state.normalizedOperator, *state.normalizedInverse,
            partition, state.Problem().GetCommunicator(), solver::linear::FGMRES(options.gravity.backend));
        solver::PreparedGravityCurvatureCorrection correction(gravitySubsystem);
        metric("compact_gravity_operator", gravitySubsystem.UsesCompactOperator());
        metric("compact_gravity_inverse", gravitySubsystem.UsesCompactPreconditioner());
        mfem::Vector compactCorrection(partition.GravityValueSize()), gravityAction(partition.GravityResidualSize());
        mfem::Vector target(partition.GravityResidualSize()), beforeGravity(partition.GravityResidualSize());
        mfem::Vector source(partition.GravityResidualSize()), compactValues(partition.GravityValueSize());
        mfem::Vector scaleAction(partition.GravityResidualSize()), rootValues(size), fullCorrection(size);
        std::uint64_t innerOperatorApplications = 0, innerInverseApplications = 0;
        double innerOperatorSeconds = 0., innerInverseSeconds = 0.;
        const auto assessmentAction = [&](const mfem::Vector &values, mfem::Vector &image) {
            const auto actionStart = d::Clock::now();
            gravitySubsystem.GetOperator().Mult(values, image);
            assessmentGravitySeconds += d::Seconds(actionStart);
            ++assessmentGravityActions;
        };

        for (std::size_t sample = 0; sample < alphas.size(); ++sample) {
            const double alpha = alphas[sample];
            std::array<std::array<d::Observation, 3>, 2> observations;
            for (int k = 0; k < 3; ++k) {
                restore("capture_base");
                refresh();
                gravitySubsystem.Refresh();
                const auto captureStart = d::Clock::now();
                correction.CaptureBase(originalResidual, directions[k]);
                captureSeconds += d::Seconds(captureStart);
                ++captureFullActions;
                mfem::Vector candidate(base), predicted(originalResidual), expectedIncrement(directions[k]);
                candidate.Add(alpha, directions[k]);
                predicted.Add(alpha, actions[k]);
                d::Save(file(std::string("prediction_") + names[k] + "_a" + std::to_string(sample) + ".txt"), predicted);
                expectedIncrement *= alpha;
                mfem::Vector actualIncrement(candidate);
                actualIncrement -= base;
                auto &raw = observations[0][k];
                auto &curved = observations[1][k];
                const auto rawStart = d::Clock::now();
                raw.prepared = prepare(candidate, "raw", names[k], alpha);
                raw.seconds = d::Seconds(rawStart);
                raw.usable = raw.prepared;
                if (raw.prepared) {
                    raw.residual = state.trialNormalizedResidual;
                    raw.remainder = raw.residual;
                    raw.remainder -= predicted;
                }
                const auto correctedStart = d::Clock::now();
                curved.innerConverged = curved.preparedGate = false;
                if (raw.prepared) {
                    refresh();
                    gravitySubsystem.Refresh();
                    correction.BuildLinearPrediction(alpha, target);
                    const auto report = correction.SolveCorrection(alpha, raw.residual, compactCorrection,
                                                                    options.gravity.control);
                    innerOperatorApplications += report.operatorApplications;
                    innerInverseApplications += report.inversePreconditionerApplications;
                    innerOperatorSeconds += report.operatorSeconds;
                    innerInverseSeconds += report.inversePreconditionerSeconds;
                    inner << sample << ',' << alpha << ',' << names[k] << ',' << int(report.status) << ','
                          << report.Converged() << ',' << report.iterations << ',' << report.rightHandSideNorm << ','
                          << report.trueResidualNorm << ',' << report.relativeTrueResidualNorm << ','
                          << (std::isfinite(report.rightHandSideNorm) && report.rightHandSideNorm >= 0.
                                  ? options.gravity.control.ConvergenceThreshold(report.rightHandSideNorm) : nan) << ','
                          << report.operatorApplications << ',' << report.inversePreconditionerApplications << ','
                          << report.solveSeconds << ',' << gravitySubsystem.UsesCompactOperator() << ','
                          << gravitySubsystem.UsesCompactPreconditioner() << '\n';
                    inner.flush();
                    curved.innerConverged = report.Converged();
                    if (report.Converged()) {
                        partition.InjectGravityValues(compactCorrection, fullCorrection);
                        d::Save(file(std::string("correction_curvature_") + names[k] + "_a" +
                                     std::to_string(sample) + ".txt"), fullCorrection);
                        assessmentAction(compactCorrection, gravityAction);
                        partition.GatherGravityResiduals(raw.residual, beforeGravity);
                        source = beforeGravity;
                        std::array<double, 8> contributions{};
                        // Match production ApplyCorrection's arithmetic guard:
                        // separate the large cancelling g/Phi terms, rather than
                        // estimating preparation error from a tiny gravity defect.
                        const auto measure = [&](const mfem::Vector &values, auto block, bool subtractSource) {
                            rootValues = 0.;
                            for (int i = layout.offset(block); i < layout.offset(block) + layout.size(block); ++i)
                                rootValues(i) = values(i);
                            partition.GatherGravityValues(rootValues, compactValues);
                            assessmentAction(compactValues, scaleAction);
                            if (subtractSource) source -= scaleAction;
                            return d::Norm(scaleAction);
                        };
                        contributions[0] = measure(candidate, gradient, true);
                        contributions[1] = measure(candidate, potential, true);
                        contributions[2] = d::Norm(source);
                        contributions[3] = measure(fullCorrection, gradient, false);
                        contributions[4] = measure(fullCorrection, potential, false);
                        contributions[5] = d::Norm(beforeGravity);
                        contributions[6] = d::Norm(target);
                        curved.correctionNorm = d::Norm(compactCorrection);
                        partition.AddGravityValues(compactCorrection, candidate);
                        // Recompute the full prepared residual after correction;
                        // algebraic A delta alone cannot certify this trial map.
                        state.normalizedOperator->DenormalizeState(candidate, state.candidatePhysicalState);
                        curved.prepared = prepare(candidate, "curvature", names[k], alpha);
                        if (curved.prepared) {
                            curved.residual = state.trialNormalizedResidual;
                            curved.remainder = curved.residual;
                            curved.remainder -= predicted;
                            mfem::Vector actualGravity(partition.GravityResidualSize());
                            partition.GatherGravityResiduals(curved.residual, actualGravity);
                            contributions[7] = d::Norm(actualGravity);
                            mfem::Vector affine(beforeGravity);
                            affine += gravityAction;
                            affine -= actualGravity;
                            actualGravity -= target;
                            curved.targetDefect = d::Norm(actualGravity);
                            curved.affineDisagreement = d::Norm(affine);
                            const auto assessment = solver::detail::AssessGravityPreparedResidual(
                                contributions, curved.affineDisagreement, curved.targetDefect, report.trueResidualNorm);
                            curved.preparedGate = assessment.admissible;
                            curved.arithmeticAllowance = assessment.arithmeticAllowance;
                            curved.targetThreshold = assessment.targetThreshold;
                            curved.usable = curved.preparedGate;
                        }
                    }
                }
                // Curvature cost includes its shared raw preparation. Summing
                // path seconds double-counts that work; total_seconds does not.
                curved.seconds = raw.seconds + d::Seconds(correctedStart);
                for (int path = 0; path < 2; ++path) {
                    const auto &observation = observations[path][k];
                    const char *pathName = path == 0 ? "raw" : "curvature";
                    mfem::Vector intendedIncrement(expectedIncrement), pathIncrement(actualIncrement);
                    if (path == 1 && observation.innerConverged) {
                        intendedIncrement += fullCorrection;
                        pathIncrement = candidate;
                        pathIncrement -= base;
                    }
                    mfem::Vector pathAdditionError(pathIncrement);
                    pathAdditionError -= intendedIncrement;
                    const bool haveIncrement = path == 0 || observation.innerConverged;
                    const double expectedPathStep = haveIncrement ? d::Norm(intendedIncrement) : nan;
                    const double actualPathStep = haveIncrement ? d::Norm(pathIncrement) : nan;
                    const bool pathResolved = haveIncrement && expectedPathStep > 64. * eps * (1. + d::Norm(base)) &&
                                              actualPathStep > 0. && d::Norm(pathAdditionError) <= .01 * expectedPathStep;
                    const double residualNorm = observation.prepared ? d::Norm(observation.residual) : nan;
                    const double remainderNorm = observation.prepared ? d::Norm(observation.remainder) : nan;
                    const double slope = static_cast<double>(d::Dot(originalResidual, actions[k], 0, size));
                    const double projection = observation.prepared
                        ? static_cast<double>(d::Dot(predicted, observation.remainder, 0, size)) : nan;
                    const auto identityError = [&](int offset, int count) {
                        return observation.prepared ? static_cast<double>(
                            d::Dot(observation.residual, observation.residual, offset, count) -
                            d::Dot(predicted, predicted, offset, count) -
                            2. * d::Dot(predicted, observation.remainder, offset, count) -
                            d::Dot(observation.remainder, observation.remainder, offset, count)) : nan;
                    };
                    const auto arithmeticFloor = [&](int offset, int count) {
                        // Only a subtraction heuristic; production's gravity
                        // gate separately measures assembled cancellation scales.
                        return observation.prepared ? std::max(8. * restorationError, 128. * eps * (
                            1. + d::Norm(originalResidual, offset, count) + alpha * d::Norm(actions[k], offset, count) +
                            d::Norm(observation.residual, offset, count))) : nan;
                    };
                    paths << sample << ',' << alpha << ',' << names[k] << ',' << pathName << ','
                          << observation.prepared << ',' << observation.innerConverged << ','
                          << observation.preparedGate << ',' << observation.usable << ',' << d::Norm(originalResidual)
                          << ',' << d::Norm(actions[k]) << ',' << residualNorm << ',' << .5 * residualNorm * residualNorm
                          << ',' << slope << ',' << d::Norm(predicted) << ',' << remainderNorm << ','
                          << remainderNorm / alpha << ',' << remainderNorm / (alpha * alpha) << ','
                          << d::Ratio(remainderNorm / alpha, d::Norm(actions[k])) << ',' << projection << ','
                          << identityError(0, size) << ',' << expectedPathStep << ',' << actualPathStep << ','
                          << (haveIncrement ? d::Norm(pathAdditionError) : nan) << ',' << pathResolved << ',' << observation.correctionNorm << ','
                          << observation.targetDefect << ',' << observation.affineDisagreement << ','
                          << observation.arithmeticAllowance << ',' << observation.targetThreshold << ','
                          << observation.seconds << ',' << arithmeticFloor(0, size) << '\n';
                    if (observation.prepared) {
                        ++successfulPreparedPaths;
                        const std::string suffix = std::string(pathName) + '_' + names[k] + "_a" + std::to_string(sample) + ".txt";
                        d::Save(file("residual_" + suffix), observation.residual);
                        // The raw state must not accidentally inherit the SOC
                        // update from the shared trial buffer.
                        mfem::Vector recordedState(base);
                        recordedState.Add(alpha, directions[k]);
                        if (path == 1) recordedState = candidate;
                        d::Save(file("state_" + suffix), recordedState);
                    }
                    eachBlock([&](auto name, int offset, int count) {
                        const double norm = observation.prepared ? d::Norm(observation.remainder, offset, count) : nan;
                        blocks << sample << ',' << alpha << ',' << names[k] << ',' << pathName << ',' << name << ','
                               << offset << ',' << count << ',' << d::Norm(originalResidual, offset, count) << ','
                               << d::Norm(actions[k], offset, count) << ',' << d::Norm(predicted, offset, count) << ','
                               << (observation.prepared ? d::Norm(observation.residual, offset, count) : nan) << ','
                               << norm << ',' << norm / alpha << ',' << norm / (alpha * alpha) << ','
                               << (observation.prepared ? static_cast<double>(d::Dot(predicted, observation.remainder, offset, count)) : nan)
                               << ',' << identityError(offset, count) << ',' << arithmeticFloor(offset, count) << '\n';
                    });
                }
                paths.flush(); blocks.flush();
            }
            // For corrected trials this is an interaction of the corrected
            // trial map, not solely the Hessian of the original affine path.
            mfem::Vector fullPrediction(originalResidual);
            fullPrediction.Add(alpha, actions[0]);
            for (int path = 0; path < 2; ++path) {
                const auto &o = observations[path];
                const bool allPrepared = o[0].prepared && o[1].prepared && o[2].prepared;
                const bool usable = o[0].usable && o[1].usable && o[2].usable;
                mfem::Vector interaction(size);
                if (allPrepared) {
                    interaction = o[0].remainder;
                    interaction -= o[1].remainder;
                    interaction -= o[2].remainder;
                }
                eachBlock([&](auto name, int offset, int count) {
                    std::array<double, 4> norms{nan, nan, nan, nan};
                    std::array<long double, 4> projections{nan, nan, nan, nan};
                    double identity = nan;
                    if (allPrepared) {
                        for (int k = 0; k < 3; ++k) {
                            norms[k] = d::Norm(o[k].remainder, offset, count);
                            projections[k] = d::Dot(fullPrediction, o[k].remainder, offset, count);
                        }
                        norms[3] = d::Norm(interaction, offset, count);
                        projections[3] = d::Dot(fullPrediction, interaction, offset, count);
                        identity = static_cast<double>(d::Dot(o[0].residual, o[0].residual, offset, count) -
                            d::Dot(fullPrediction, fullPrediction, offset, count) -
                            2. * (projections[1] + projections[2] + projections[3]) -
                            d::Dot(o[0].remainder, o[0].remainder, offset, count));
                    }
                    mixed << sample << ',' << alpha << ',' << (path == 0 ? "raw" : "curvature") << ',' << name << ','
                          << offset << ',' << count << ',' << usable;
                    for (double norm : norms) mixed << ',' << norm;
                    for (long double projection : projections) mixed << ',' << projection;
                    mixed << ',' << projections[0] - projections[1] - projections[2] - projections[3] << ',' << identity << '\n';
                });
            }
            mixed.flush(); summary.flush();
        }
        metric("gravity_solves", gravitySubsystem.GetSolveCount());
        metric("gravity_iterations", gravitySubsystem.GetIterationCount());
        metric("gravity_failures", gravitySubsystem.GetFailureCount());
        metric("gravity_solve_seconds", gravitySubsystem.GetSolveSeconds());
        metric("gravity_solve_operator_applications", innerOperatorApplications);
        metric("gravity_solve_inverse_applications", innerInverseApplications);
        metric("gravity_solve_operator_seconds", innerOperatorSeconds);
        metric("gravity_solve_inverse_seconds", innerInverseSeconds);
        double maximumActionChange = 0.;
        for (int k = 0; k < 3; ++k) {
            restore("final_split_action_base");
            mfem::Vector restoredAction(size);
            action(directions[k], restoredAction);
            maximumActionChange = std::max(maximumActionChange, d::Difference(actions[k], restoredAction));
            d::Save(file(std::string("restored_action_") + names[k] + ".txt"), restoredAction);
        }
        metric("maximum_restored_action_difference", maximumActionChange);
    } catch (...) {
        const auto failure = std::current_exception();
        finish();
        std::rethrow_exception(failure);
    }
    finish();
}
} // namespace experiment::plateau_diagnostics