#include <algorithm>
#include <array>
#include <chrono>
#include <cmath>
#include <cstdlib>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <limits>
#include <numbers>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>
#include <mfem.hpp>

import mean_field;

#include "geometry_quality_diagnostics.hpp"

namespace {
using namespace mean_field;
using Clock = std::chrono::steady_clock;

struct Options {
    std::string mesh = "sandbox.smesh";
    std::filesystem::path output = "geometry_quality_results";
    std::vector<double> tolerances{0.03, 0.003};
    int maximumLinearIterations = 200;
    int advance = 0;
    bool finiteDifferences = true;
    bool blockActions = true;
    bool geometryOnly = false;
    bool saveVectors = false;
    bool warm = false;
    std::filesystem::path replayVectors;
    bool diagonalOnly=false;
};

Options Parse(int argc, char **argv) {
    Options o;
    for (int i = 1; i < argc; ++i) {
        const std::string arg = argv[i];
        auto value = [&]() -> std::string {
            if (++i >= argc) throw std::invalid_argument("Missing value after " + arg);
            return argv[i];
        };
        if (arg == "--mesh") o.mesh = value();
        else if (arg == "--output") o.output = value();
        else if (arg == "--tolerances") {
            o.tolerances.clear();
            std::istringstream stream(value());
            for (std::string token; std::getline(stream, token, ',');) {
                const double tolerance = std::stod(token);
                if (!(tolerance > 0.0 && tolerance < 1.0)) throw std::invalid_argument("Invalid tolerance");
                o.tolerances.push_back(tolerance);
            }
            if (o.tolerances.empty()) throw std::invalid_argument("Empty tolerance list");
        } else if (arg == "--max-linear-iterations") o.maximumLinearIterations = std::stoi(value());
        else if (arg == "--advance") o.advance = std::stoi(value());
        else if (arg == "--no-fd") o.finiteDifferences = false;
        else if (arg == "--no-block-actions") o.blockActions = false;
        else if (arg == "--geometry-only") o.geometryOnly = true;
        else if (arg == "--save-vectors") o.saveVectors = true;
        else if (arg == "--warm") o.warm = true;
        else if (arg == "--replay-vectors") o.replayVectors=value();
        else if (arg == "--diagonal-only") o.diagonalOnly=true;
        else if (arg == "--help") {
            std::cout << "geometry_quality_experiment [--mesh FILE] [--output NEW_DIRECTORY]\n"
                         "  [--tolerances 0.03,0.003] [--max-linear-iterations 200] [--advance N]\n"
                         "  [--warm] [--geometry-only] [--no-fd] [--no-block-actions] [--save-vectors]\n"
                         "  [--replay-vectors CASE_vectors.csv] (verifies saved accepted state, no linear solve)\n"
                         "  [--diagonal-only] (seed replay only; skips per-element scans and residual checks)\n"
                         "Single MPI rank only. Defaults inspect a frozen seed; --advance uses production Newton.\n";
            std::exit(0);
        } else throw std::invalid_argument("Unknown option: " + arg);
    }
    if (o.advance < 0 || o.maximumLinearIterations < 1) throw std::invalid_argument("Invalid iteration count");
    if (!o.replayVectors.empty()) o.tolerances.resize(1);
    if (o.diagonalOnly) {
        if (o.replayVectors.empty() || o.advance!=0 || o.geometryOnly) throw std::invalid_argument("--diagonal-only requires seed correction replay");
        o.finiteDifferences=false; o.blockActions=false;
    }
    return o;
}

std::ofstream File(const std::filesystem::path &path) {
    std::ofstream stream(path);
    stream.exceptions(std::ios::failbit | std::ios::badbit);
    stream << std::setprecision(17);
    return stream;
}

double Norm(const mfem::Vector &v, int offset, int size) {
    long double sum = 0;
    for (int i = offset; i < offset + size; ++i) sum += static_cast<long double>(v(i)) * v(i);
    return std::sqrt(sum);
}

double MaxAbs(const mfem::Vector &v, int offset, int size) {
    double result = 0;
    for (int i = offset; i < offset + size; ++i) result = std::max(result, std::abs(v(i)));
    return result;
}

template <typename State>
void Inspect(State &s, const fem::FEM &fem, const Options &options) {
    auto &problem = s.Problem();
    const auto &deformation = problem.GetPhysicalOperator().GetDomainDeformation();
    // Recompile only the public surface-coordinate descriptor to inspect its
    // radii/directions; the actual extension always remains the context's.
    mfem::Vector center(3); center=0.0;
    const deformation::SurfaceDeformationCompilationContext surfaceContext{
        *fem.surfaceDeformationFes,
        field::make_stellar_surface_scalar_dof_map<utils::domain::CoreEnvelopeVacuumDomainSchema>(*fem.surfaceDeformationFes)};
    const auto surface=deformation::compileSurfaceDeformationPrescription(deformation::NodalRadialSurface{center},surfaceContext);
    const auto values = problem.GetManifest().valueBlocks();
    const auto rows = problem.GetManifest().residualBlocks();
    const auto surfaceIterator = std::find_if(values.begin(), values.end(), [](const auto &b) {
        return b.stableId == "surface_deformation";
    });
    if (surfaceIterator == values.end()) throw std::logic_error("Experiment requires a surface-deformation block");
    const auto surfaceBlock = *surfaceIterator;
    const mfem::Vector acceptedVolume(problem.GetPhysicalOperator().GetGeneratedVolumeDisplacement());
    const mfem::Vector physicalResidual(s.normalizedOperator->GetPhysicalResidual());
    const mfem::Vector warmCorrection(s.normalizedCorrection);
    auto residuals = File(options.output / "blocks.csv");
    residuals << "case,kind,block,size,physical_l2,normalized_l2,physical_linf,normalized_linf,linear_residual_over_block_F\n";
    auto record = [&](const std::string &label, const std::string &kind, auto blocks,
                      const mfem::Vector &physical, const mfem::Vector &normalized, bool relative) {
        for (const auto &b : blocks) {
            const double denominator = relative ? Norm(s.acceptedNormalizedResidual, b.offset, b.size) : 0.0;
            residuals << label << ',' << kind << ',' << b.stableId << ',' << b.size << ','
                      << Norm(physical, b.offset, b.size) << ',' << Norm(normalized, b.offset, b.size) << ','
                      << MaxAbs(physical, b.offset, b.size) << ',' << MaxAbs(normalized, b.offset, b.size) << ','
                      << (relative && denominator > 0 ? Norm(normalized, b.offset, b.size) / denominator
                                                     : std::numeric_limits<double>::quiet_NaN()) << '\n';
        }
        residuals.flush();
    };
    record("accepted", "residual", rows, physicalResidual, s.acceptedNormalizedResidual, false);
    record("accepted", "state", values, s.AcceptedPhysicalState(), s.acceptedNormalizedState, false);

    auto metadata = File(options.output / "metadata.txt");
    metadata << "mesh=" << options.mesh << "\nmesh_bytes=" << std::filesystem::file_size(options.mesh)
             << "\ncompiled=" << __DATE__ << ' ' << __TIME__ << "\ncompiler=" << __VERSION__
             << "\npolynomial_increment=" << MEAN_FIELD_UNIFORM_POLYNOMIAL_ORDER_INCREMENT
             << "\nmpi_ranks=1\nmodel=n1_K2G_R2_over_pi_M1_R1_J0_fixed_central_density"
             << "\nnormalization=production_frozen_physical_Riesz_diagonal"
             << "\npreconditioner=production_default\nFGMRES_restart=40\nadvanced_steps=" << options.advance
             << "\nmax_linear_iterations=" << options.maximumLinearIterations
             << "\nwarm=" << options.warm << "\nstate_size=" << problem.StateSize()
             << "\nelement_count=" << fem.mesh->GetNE() << "\naccepted_residual=" << s.acceptedNormalizedResidual.Norml2()
             << "\ncoefficient_norms_are_unweighted_single_rank=true"
             << "\nsurface_mean_and_rms_are_nodal_not_area_weighted=true\ntolerances=";
    for (double t : options.tolerances) metadata << t << ',';
    metadata << '\n';
    metadata << "replay_vectors=" << options.replayVectors.string() << '\n';
    metadata << "diagonal_only=" << options.diagonalOnly << '\n';
    metadata.flush();

    auto vertexGeometry = [&](const std::string &label,const mfem::Vector &direction) {
        // Deliberately independent sampling of the known sandbox core corners.
        // Uses the unmodified production estimator with a vertex-only rule.
        std::vector<deformation::NewtonStepGeometryRule> rules;
        for (int element : {0,9,18,27,36,45,54,63}) {
            if (element<fem.mesh->GetNE() && fem.mesh->GetAttribute(element)==1 &&
                fem.mesh->GetElementBaseGeometry(element)==mfem::Geometry::CUBE)
                rules.push_back({element,mfem::Geometries.GetVertices(mfem::Geometry::CUBE)});
        }
        if (!rules.empty()) experiment::InspectGeometry(problem.GetDiscretization().domainMapper(),*fem.displacementFes,
            *fem.compactificationCoordinate,acceptedVolume,direction,rules,options.output,label+"_vertices");
    };

    auto geometry = [&](const std::string &label, const mfem::Vector &physicalDirection) {
        mfem::Vector volumeDirection(acceptedVolume.Size());
        problem.BuildVolumeDisplacementDirection(physicalDirection, volumeDirection);
        std::cout << "Geometry: " << label << std::endl;
        if (options.advance==0 && (label=="newton_0" || label=="uniform_contraction")) {
            int corner=0;
            double minimumSum=std::numeric_limits<double>::infinity();
            for (int i=0;i<surfaceBlock.size;++i) {
                double sum=0;
                for (int d=0;d<3;++d) sum+=surface.radialDirection(i,d);
                if (sum<minimumSum) { minimumSum=sum; corner=i; }
            }
            experiment::InspectCoreDiagonal(fem,volumeDirection,physicalDirection(surfaceBlock.offset+corner),options.output,label);
        }
        if (options.advance==0 && !options.replayVectors.empty()) vertexGeometry(label,volumeDirection);
        if (options.diagonalOnly) return experiment::GeometryReport{};
        return experiment::InspectGeometry(problem.GetDiscretization().domainMapper(), *fem.displacementFes,
            *fem.compactificationCoordinate, acceptedVolume, volumeDirection, s.geometryPreflightRules,
            options.output, label);
    };

    mfem::Vector contraction(problem.StateSize());
    contraction = 0.0;
    for (int i = 0; i < surfaceBlock.size; ++i) contraction(surfaceBlock.offset + i) = -surface.referenceRadius(i);
    geometry("uniform_contraction", contraction);
    // Controls, not candidate production prescriptions: bypass the surface
    // extension and inspect only core (attribute 1) samples. P3 interpolation
    // cannot represent the P4 physical mesh coordinates exactly.
    std::vector<deformation::NewtonStepGeometryRule> coreRules;
    for (const auto &rule : s.geometryPreflightRules) {
        if (fem.mesh->GetAttribute(rule.element)==1) coreRules.push_back(rule);
    }
    if (!options.diagonalOnly) for (bool radial : {false,true}) {
        mfem::VectorFunctionCoefficient coefficient(3,[radial](const mfem::Vector &x,mfem::Vector &u) {
            u=x;
            u *= radial ? -x.Norml2()/utils::RADIUS : -1.0;
        });
        mfem::ParGridFunction projected(fem.displacementFes.get());
        projected.ProjectCoefficient(coefficient);
        mfem::Vector direction; projected.GetTrueDofs(direction);
        experiment::InspectGeometry(problem.GetDiscretization().domainMapper(), *fem.displacementFes,
            *fem.compactificationCoordinate, acceptedVolume, direction, coreRules, options.output,
            radial ? "core_projected_physical_radial_contraction" : "core_projected_affine_contraction");
    }
    if (options.geometryOnly) return;

    auto solves = File(options.output / "solves.csv");
    solves << "case,tolerance,status,iterations,relative_true_residual,wall_seconds,safe_step,boundary_step,limiting_element,correction_difference_from_baseline,surface_difference_from_baseline,action_difference_over_F\n";
    auto surfaceFile = File(options.output / "surface.csv");
    surfaceFile << "case,parameter,x,y,z,radius,accepted_fraction,correction_fraction\n";
    auto surfaceSummary = File(options.output / "surface_summary.csv");
    surfaceSummary << "case,mean_fraction,rms_fraction,rms_nonmean_fraction,min_fraction,max_fraction\n";
    auto columns = File(options.output / "block_actions.csv");
    columns << "case,column,row,normalized_action_l2,normalized_dot_with_F\n";
    auto fd = File(options.output / "finite_differences.csv");
    fd << "case,epsilon,row,relative_error,absolute_error,action_norm\n";
    auto extensionChecks = File(options.output / "extension_checks.csv");
    extensionChecks << "case,epsilon,relative_volume_direction_error,trial_residual_norm\n";
    mfem::Vector baseline, baselineAction;

    for (std::size_t caseIndex = 0; caseIndex < options.tolerances.size(); ++caseIndex) {
        const std::string label = "newton_" + std::to_string(caseIndex);
        const double tolerance = options.tolerances[caseIndex];
        s.normalizedCorrection = options.warm ? warmCorrection : mfem::Vector(problem.StateSize());
        if (!options.warm) s.normalizedCorrection = 0.0;
        s.linearRightHandSide = s.acceptedNormalizedResidual;
        s.linearRightHandSide *= -1;
        std::cout << "Solve: " << label << " tolerance=" << tolerance << std::endl;
        const auto start = Clock::now();
        solver::LinearSolveReport solve;
        if (options.replayVectors.empty()) {
            solve=s.linearBackend->Solve(s.linearRightHandSide, s.normalizedCorrection,
                {.relativeTolerance=tolerance, .absoluteTolerance=0.0, .maximumIterations=options.maximumLinearIterations});
        } else {
            std::ifstream saved(options.replayVectors);
            if (!saved) throw std::runtime_error("Cannot open replay vectors");
            std::string line; std::getline(saved,line);
            if (line!="index,physical_state,normalized_state,physical_correction,normalized_correction")
                throw std::runtime_error("Unrecognized replay vector header");
            int index=0;
            while (std::getline(saved,line)) {
                std::replace(line.begin(),line.end(),',',' ');
                std::istringstream row(line);
                int savedIndex; double physicalState, normalizedState, physicalCorrection, normalizedCorrection;
                if (!(row>>savedIndex>>physicalState>>normalizedState>>physicalCorrection>>normalizedCorrection) ||
                    savedIndex!=index || index>=problem.StateSize()) throw std::runtime_error("Invalid replay vector row");
                if (std::abs(physicalState-s.AcceptedPhysicalState()(index))>1e-12*(1+std::abs(physicalState)) ||
                    std::abs(normalizedState-s.acceptedNormalizedState(index))>1e-12*(1+std::abs(normalizedState)))
                    throw std::runtime_error("Replay state differs from this frozen context");
                if (!std::isfinite(normalizedCorrection)) throw std::runtime_error("Nonfinite replay direction");
                s.normalizedCorrection(index++)=normalizedCorrection;
            }
            if (index!=problem.StateSize()) throw std::runtime_error("Incomplete replay vector file");
            mfem::Vector check(problem.EquationSize());
            s.normalizedOperator->Mult(s.normalizedCorrection,check);
            check+=s.acceptedNormalizedResidual;
            solve.relativeTrueResidualNorm=check.Norml2()/s.acceptedNormalizedResidual.Norml2();
            solve.status=solve.relativeTrueResidualNorm<=tolerance ? solver::LinearSolveStatus::converged : solver::LinearSolveStatus::maximum_iterations;
            solve.iterations=0;
            std::cout << "Replayed saved correction; true residual independently checked" << std::endl;
        }
        const double elapsed = std::chrono::duration<double>(Clock::now()-start).count();
        std::cout << "Solved: iterations=" << solve.iterations << " true_relative=" << solve.relativeTrueResidualNorm
                  << " wall=" << elapsed << std::endl;
        s.normalizedOperator->DenormalizeState(s.normalizedCorrection, s.physicalCorrection);
        const auto safe = s.EstimateLargestSafeStepSize(1.0, 0.9);
        geometry(label, s.physicalCorrection);
        if (caseIndex==0 && options.advance==0 && !options.replayVectors.empty()) {
            mfem::Vector zeroSurface(surfaceBlock.size), surfaceDirection(surfaceBlock.size);
            zeroSurface=0.0;
            for (int i=0;i<surfaceBlock.size;++i) surfaceDirection(i)=s.physicalCorrection(surfaceBlock.offset+i);
            const auto extensionContext=deformation::makeRadialDeformationExtensionCompilationContext(
                *fem.surfaceDeformationFes,*fem.displacementFes,*fem.logicalReferenceMesh);
            for (double power : {3.0,4.0}) {
                // Existing production prescription, but a geometry-only
                // intervention: this is NOT a Newton direction for the new
                // parameterization until that system is rebuilt and solved.
                auto alternativeSurface=deformation::compileSurfaceDeformationPrescription(
                    deformation::NodalRadialSurface{center},surfaceContext);
                auto alternativeInterior=deformation::compileInteriorDeformationExtension(
                    deformation::PowerLawRadialInteriorExtension{power},extensionContext);
                auto alternativeVacuum=deformation::compileVacuumDeformationExtension(
                    deformation::FixedInfinityRadialVacuumExtension{},extensionContext);
                auto alternative=deformation::composePreparedDomainDeformation(
                    std::move(alternativeSurface),std::move(alternativeInterior),std::move(alternativeVacuum),
                    *fem.surfaceDeformationFes,*fem.displacementFes,*fem.logicalReferenceMesh);
                mfem::Vector direction(acceptedVolume.Size());
                alternative.applyJacobian(zeroSurface,surfaceDirection,direction);
                const std::string control="newton_0_radial_power_"+std::to_string(static_cast<int>(power));
                int corner=0; double smallest=std::numeric_limits<double>::infinity();
                for (int i=0;i<surfaceBlock.size;++i) {
                    double sum=0; for (int d=0;d<3;++d) sum+=surface.radialDirection(i,d);
                    if (sum<smallest) { smallest=sum; corner=i; }
                }
                experiment::InspectCoreDiagonal(fem,direction,surfaceDirection(corner),options.output,control,power);
                vertexGeometry(control,direction);
                if (!options.diagonalOnly) experiment::InspectGeometry(problem.GetDiscretization().domainMapper(),*fem.displacementFes,
                    *fem.compactificationCoordinate,acceptedVolume,direction,s.geometryPreflightRules,options.output,control);
            }
        }
        mfem::Vector action(problem.EquationSize()), linearResidual(problem.EquationSize()), physicalLinear(problem.EquationSize());
        s.normalizedOperator->Mult(s.normalizedCorrection, action);
        linearResidual = action;
        linearResidual += s.acceptedNormalizedResidual;
        s.normalizedOperator->DenormalizeResidual(linearResidual, physicalLinear);
        record(label, "linear_residual", rows, physicalLinear, linearResidual, true);
        record(label, "correction", values, s.physicalCorrection, s.normalizedCorrection, false);
        if (caseIndex == 0) { baseline = s.normalizedCorrection; baselineAction=action; }
        mfem::Vector difference(s.normalizedCorrection);
        difference -= baseline;
        mfem::Vector actionDifference(action); actionDifference-=baselineAction;
        solves << label << ',' << tolerance << ',' << static_cast<int>(solve.status) << ',' << solve.iterations << ','
               << solve.relativeTrueResidualNorm << ',' << elapsed << ',' << safe.stepSize << ',' << safe.boundaryStepSize
               << ',' << safe.limitingElement << ',' << difference.Norml2()/std::max(baseline.Norml2(),1e-300) << ','
               << Norm(difference,surfaceBlock.offset,surfaceBlock.size)/std::max(Norm(baseline,surfaceBlock.offset,surfaceBlock.size),1e-300) << ','
               << actionDifference.Norml2()/std::max(s.acceptedNormalizedResidual.Norml2(),1e-300) << '\n';
        solves.flush();
        double sum=0, squared=0, minimum=std::numeric_limits<double>::infinity(), maximum=-minimum;
        for (int i=0; i<surfaceBlock.size; ++i) {
            const double radius=surface.referenceRadius(i);
            const double fraction=s.physicalCorrection(surfaceBlock.offset+i)/radius;
            sum += fraction; squared += fraction*fraction; minimum=std::min(minimum,fraction); maximum=std::max(maximum,fraction);
            surfaceFile << label << ',' << i;
            for (int d=0; d<3; ++d) surfaceFile << ',' << surface.referenceCenter()(d)+radius*surface.radialDirection(i,d);
            surfaceFile << ',' << radius << ',' << s.AcceptedPhysicalState()(surfaceBlock.offset+i)/radius << ',' << fraction << '\n';
        }
        const double mean=sum/surfaceBlock.size, meanSquared=squared/surfaceBlock.size;
        surfaceSummary << label << ',' << mean << ',' << std::sqrt(meanSquared) << ','
                       << std::sqrt(std::max(0.0,meanSquared-mean*mean)) << ',' << minimum << ',' << maximum << '\n';
        surfaceFile.flush(); surfaceSummary.flush();

        if (caseIndex == 0) {
            // Split only the surface component; other physical fields do not
            // enter the production volume-extension Jacobian.
            mfem::Vector meanDirection(problem.StateSize()), nonmeanDirection(s.physicalCorrection);
            meanDirection = 0.0;
            for (int i=0;i<surfaceBlock.size;++i) {
                meanDirection(surfaceBlock.offset+i)=mean*surface.referenceRadius(i);
                nonmeanDirection(surfaceBlock.offset+i)-=meanDirection(surfaceBlock.offset+i);
            }
            geometry(label+"_mean",meanDirection);
            geometry(label+"_nonmean",nonmeanDirection);
        }

        if (options.saveVectors) {
            auto vectors=File(options.output/(label+"_vectors.csv"));
            vectors << "index,physical_state,normalized_state,physical_correction,normalized_correction\n";
            for (int i=0;i<problem.StateSize();++i) vectors << i << ',' << s.AcceptedPhysicalState()(i) << ','
                << s.acceptedNormalizedState(i) << ',' << s.physicalCorrection(i) << ',' << s.normalizedCorrection(i) << '\n';
        }
        if (options.blockActions && caseIndex == 0) {
            std::cout << "Block-column actions" << std::endl;
            for (const auto &column : values) {
                mfem::Vector direction(problem.StateSize()), blockAction(problem.EquationSize());
                direction = 0.0;
                for (int i=column.offset;i<column.offset+column.size;++i) direction(i)=s.normalizedCorrection(i);
                s.normalizedOperator->Mult(direction,blockAction);
                for (const auto &row : rows) {
                    double dot=0;
                    for (int i=row.offset;i<row.offset+row.size;++i) dot+=blockAction(i)*s.acceptedNormalizedResidual(i);
                    columns << label << ',' << column.stableId << ',' << row.stableId << ','
                            << Norm(blockAction,row.offset,row.size) << ',' << dot << '\n';
                }
                columns.flush();
            }
        }
        if (options.finiteDifferences && caseIndex == 0) {
            // Forward differences remain inside the verified positive-alpha interval.
            // Frozen normalization and fresh dependency revisions match production trial preparation.
            mfem::Vector expectedVolumeDirection(acceptedVolume.Size());
            problem.BuildVolumeDisplacementDirection(s.physicalCorrection,expectedVolumeDirection);
            for (double multiplier : {1e-2,1e-3,1e-4}) {
                const double epsilon = std::min(1.0,safe.stepSize)*multiplier;
                if (!(epsilon>0)) throw std::runtime_error("No admissible finite-difference step");
                std::cout << "Finite difference epsilon=" << epsilon << std::endl;
                s.trialNormalizedState=s.acceptedNormalizedState;
                s.trialNormalizedState.Add(epsilon,s.normalizedCorrection);
                const auto preparation=s.PrepareTrial();
                if (!preparation) throw std::runtime_error("Finite-difference trial preparation rejected");
                mfem::Vector volumeError(problem.GetPhysicalOperator().GetGeneratedVolumeDisplacement());
                volumeError-=acceptedVolume;
                volumeError/=epsilon;
                volumeError-=expectedVolumeDirection;
                extensionChecks << label << ',' << epsilon << ','
                    << volumeError.Norml2()/std::max(expectedVolumeDirection.Norml2(),1e-300) << ','
                    << s.trialNormalizedResidual.Norml2() << '\n';
                extensionChecks.flush();
                mfem::Vector approximation(s.trialNormalizedResidual);
                approximation-=s.acceptedNormalizedResidual;
                approximation/=epsilon;
                approximation-=action;
                for (const auto &row:rows) {
                    const double error=Norm(approximation,row.offset,row.size), magnitude=Norm(action,row.offset,row.size);
                    fd << label << ',' << epsilon << ',' << row.stableId << ','
                       << (magnitude>0 ? error/magnitude : std::numeric_limits<double>::quiet_NaN()) << ',' << error << ',' << magnitude << '\n';
                }
                fd.flush();
            }
            s.RestoreAccepted();
        }
    }
}
} // namespace

int main(int argc, char **argv) {
    mfem::Mpi::Init(argc,argv);
    try {
        const auto options=Parse(argc,argv);
        int ranks=0; MPI_Comm_size(MPI_COMM_WORLD,&ranks);
        if (ranks!=1) throw std::invalid_argument("This diagnostic executable currently requires exactly one MPI rank");
        if (std::filesystem::exists(options.output)) throw std::invalid_argument("Output directory already exists; use a new --output path");
        std::filesystem::create_directories(options.output);
        mfem::Device device("cpu");
        utils::Args args; args.mesh_file=options.mesh; args.p.rtol=1e-12; args.p.atol=1e-12;
        const auto start=Clock::now();
        auto fem=fem::setup_fem(args.mesh_file,args,0);
        if (!fem.okay()) throw std::runtime_error("Finite-element setup failed");
        constexpr double radius=utils::RADIUS, mass=utils::MASS, G=utils::G;
        auto model=model::StellarModel(eos::Polytrope({.n=1.0,.K=2*G*radius*radius/std::numbers::pi}),
            surface::Isobaric({.Psurf=dimensions::PressureValue{0}}),
            integral::FixedTotalMass({.Mtotal=dimensions::MassValue{mass}}),
            integral::FixedAngularMomentum({.Jtotal=dimensions::AngularMomentumValue{0},.axis={0,0,1},.center={0,0,0}}),
            constraint::FixedCentralDensity({.RhoC=dimensions::DensityValue{std::numbers::pi*mass/(4*radius*radius*radius)}}));
        auto discretization=equilibrium::makeStellarDiscretization(std::move(fem),
            normalization::PhysicalRieszDiagonal{dimensions::LengthValue{radius},G});
        std::cout << "Constructing production context" << std::endl;
        auto context=solver::makeContext(std::move(model),std::move(discretization),
            preconditioning::makePreconditioner(),solver::linear::FGMRES({.restartLength=40,.printLevel=-1}));
        std::cout << "Setup seconds=" << std::chrono::duration<double>(Clock::now()-start).count() << std::endl;
        if (options.advance>0) {
            auto trajectory=File(options.output/"trajectory.csv");
            trajectory << "iteration,residual,step,limiting_element,boundary\n";
            auto observer=solver::nonlinear::makeObserver([](const solver::nonlinear::BeforeIteration &) {},
                [&](const solver::nonlinear::AfterIteration &event) {
                trajectory << event.iteration << ',' << event.residualNorm << ',' << event.acceptedStepLength << ','
                           << (event.geometryPreflight ? event.geometryPreflight->limitingElement : -1) << ','
                           << (event.geometryPreflight ? event.geometryPreflight->boundaryStepSize : 0) << '\n';
                trajectory.flush();
                std::cout << "Advance iteration=" << event.iteration << " residual=" << event.residualNorm << " step=" << event.acceptedStepLength << std::endl;
            });
            auto method=solver::nonlinear::Newton(solver::nonlinear::NewtonOptions{
                .relativeTolerance=1e-8,.absoluteTolerance=0,.maximumIterations=options.advance,
                .linearSolve={.relativeTolerance=0.03,.absoluteTolerance=0,.maximumIterations=200},.backtracking={}});
            auto solve=solver::make(context,method,observer);
            auto report=solve.evaluate();
            std::cout << "Production accepted steps=" << report.completedNonlinearIterations() << std::endl;
        }
        solver::detail::StellarEquilibriumContextDiagnostics::WithState(context,[&](auto &state, auto &fem) { Inspect(state,fem,options); });
        if (!context.isReady() || context.hasActiveSolver()) throw std::logic_error("Diagnostic access did not restore the context");
        std::cout << "Experiment complete: " << options.output << " total_seconds="
                  << std::chrono::duration<double>(Clock::now()-start).count() << std::endl;
        return 0;
    } catch (const std::exception &error) {
        std::cerr << "geometry_quality_experiment: " << error.what() << '\n';
        return 1;
    }
}