From 7e1ab69fcd191717bdf88c918c35ae6ddc725f23 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 17 Sep 2026 10:51:07 -0400 Subject: [PATCH] Add read-only polynomial endpoint certification mode --- experiments/plateau_diagnostic.cpp | 23 ++++- .../plateau_polynomial_certificate.hpp | 92 +++++++++++++++++++ 2 files changed, 112 insertions(+), 3 deletions(-) create mode 100644 experiments/plateau_polynomial_certificate.hpp diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp index fb27fa8..3ea38a4 100644 --- a/experiments/plateau_diagnostic.cpp +++ b/experiments/plateau_diagnostic.cpp @@ -23,6 +23,7 @@ import mean_field; #include "plateau_numerical_strategy_diagnostics.hpp" #include "plateau_gravity_application_diagnostics.hpp" #include "plateau_positivity_strategy_diagnostics.hpp" +#include "plateau_polynomial_certificate.hpp" using namespace mean_field; namespace { @@ -31,7 +32,7 @@ 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; + bool strategyProbe = false, positivityProbe = false, positivityCertificate = false; std::string strategy = "baseline"; double gravityTolerance = 1e-10; int gravityIterations = 300, gravityRestart = 80; @@ -89,6 +90,8 @@ Options Parse(int argc, char **argv) { o.strategyProbe = true; else if (a == "--positivity-probe") o.positivityProbe = true; + else if (a == "--positivity-certificate") + o.positivityCertificate = true; else if (a == "--pressure-probe") o.pressureProbe = true; else if (a == "--shape-curvature") @@ -121,6 +124,11 @@ Options Parse(int argc, char **argv) { throw std::invalid_argument("--probe-restart must be 40 or 160"); if ((o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || o.positivityProbe) && (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.probeDirection.empty() || !o.replayCandidate.empty())) + throw std::invalid_argument("--positivity-certificate is a certificate-only mode; omit other probe options"); 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" && @@ -389,7 +397,8 @@ template int RunContext(Context &context, const Options &o, cons 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.probeDirection.empty() ? "direct_mfem_fgmres_cold" : "supplied_direction_no_linear_solve") + << "\nfrozen_probe_backend=" << (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 @@ -405,10 +414,18 @@ template int RunContext(Context &context, const Options &o, cons if (!o.replayCandidate.empty()) metadata << "replay_candidate=" << o.replayCandidate << '\n'; metadata << "pressure_probe=" << o.pressureProbe << "\nshape_curvature=" << o.shapeCurvature - << "\nsurface_trace=" << o.surfaceTrace << '\n'; + << "\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) diff --git a/experiments/plateau_polynomial_certificate.hpp b/experiments/plateau_polynomial_certificate.hpp new file mode 100644 index 0000000..4c54e08 --- /dev/null +++ b/experiments/plateau_polynomial_certificate.hpp @@ -0,0 +1,92 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +// Include after import mean_field. This callback only reads the supplied state; +// its WithState owner performs the usual accepted-state restoration afterward. +namespace experiment::plateau_diagnostics { +namespace polynomial_certificate_detail { +inline const char *Status(mean_field::solver::PolynomialSignStatus status) { + using Status = mean_field::solver::PolynomialSignStatus; + switch (status) { + case Status::certified_nonnegative: return "certified_nonnegative"; + case Status::unresolved: return "unresolved"; + case Status::negative: return "negative"; + case Status::non_finite: return "non_finite"; + } + return "unknown"; +} +inline std::ofstream Open(const std::filesystem::path &path) { + if (std::filesystem::exists(path)) + throw std::runtime_error("Refusing to overwrite a polynomial endpoint certificate."); + std::ofstream stream(path); + stream.exceptions(std::ios::badbit | std::ios::failbit); + stream << std::setprecision(17); + return stream; +} +inline bool Identical(const mfem::Vector &left,const mfem::Vector &right) { + return left.Size() == right.Size() && (left.Size() == 0 || + std::memcmp(left.GetData(),right.GetData(),static_cast(left.Size())*sizeof(double)) == 0); +} +} + +template +void CertifyPolynomialEnthalpy(State &state,const mfem::Vector &fullNormalized, + const std::filesystem::path &output) { + namespace detail = polynomial_certificate_detail; + using Clock = std::chrono::steady_clock; + int ranks = 0; + if (MPI_Comm_size(state.Problem().GetCommunicator(),&ranks) != MPI_SUCCESS || ranks != 1) + throw std::invalid_argument("Polynomial endpoint output requires exactly one MPI rank."); + const mfem::Vector acceptedPhysical(state.AcceptedPhysicalState()), + acceptedNormalized(state.acceptedNormalizedState), acceptedResidual(state.acceptedNormalizedResidual), + correction(state.normalizedCorrection), trialState(state.trialNormalizedState), + trialResidual(state.trialNormalizedResidual), candidatePhysical(state.candidatePhysicalState); + const auto generation = state.Problem().GetPreparationGeneration(); + const auto start = Clock::now(); + mean_field::solver::PreparedStellarEnthalpyStrategy prepared(state); + const auto certificateStart = Clock::now(); + const auto report = prepared.CheckState(fullNormalized); + const double checkSeconds = std::chrono::duration(Clock::now()-certificateStart).count(); + const double workspaceSeconds = std::chrono::duration(certificateStart-start).count(); + auto certificate = detail::Open(output / "polynomial_endpoint_certificate.csv"); + certificate << "status,lower_bound,upper_bound,visited_subcells,maximum_depth,witness_element,x,y,z,certificate_seconds,workspace_seconds\n" + << detail::Status(report.status) << ',' << report.lowerBound << ',' << report.upperBound << ',' + << report.visitedSubcells << ',' << report.maximumDepth << ',' << report.witnessElement << ',' + << report.witnessReferencePoint[0] << ',' << report.witnessReferencePoint[1] << ',' + << report.witnessReferencePoint[2] << ',' << checkSeconds << ',' << workspaceSeconds << '\n'; + const bool sameAcceptedPhysical = detail::Identical(acceptedPhysical,state.AcceptedPhysicalState()); + const bool sameAcceptedNormalized = detail::Identical(acceptedNormalized,state.acceptedNormalizedState); + const bool sameAcceptedResidual = detail::Identical(acceptedResidual,state.acceptedNormalizedResidual); + const bool sameCorrection = detail::Identical(correction,state.normalizedCorrection); + const bool sameTrialState = detail::Identical(trialState,state.trialNormalizedState); + const bool sameTrialResidual = detail::Identical(trialResidual,state.trialNormalizedResidual); + const bool sameCandidate = detail::Identical(candidatePhysical,state.candidatePhysicalState); + const bool samePreparation = generation == state.Problem().GetPreparationGeneration(); + auto invariants = detail::Open(output / "polynomial_endpoint_invariants.csv"); + invariants << "metric,value\naccepted_physical_bitwise_identical," << sameAcceptedPhysical + << "\naccepted_normalized_bitwise_identical," << sameAcceptedNormalized + << "\naccepted_residual_bitwise_identical," << sameAcceptedResidual + << "\ncorrection_bitwise_identical," << sameCorrection + << "\ntrial_state_bitwise_identical," << sameTrialState + << "\ntrial_residual_bitwise_identical," << sameTrialResidual + << "\ncandidate_physical_bitwise_identical," << sameCandidate + << "\npreparation_generation_unchanged," << samePreparation + << "\nhelper_trial_preparations,0\nhelper_jacobian_applications,0\nhelper_linear_solves,0\n"; + auto scope = detail::Open(output / "polynomial_endpoint_scope.txt"); + scope << "The helper constructs a positivity workspace and certifies only the supplied normalized state.\n" + << "It does not prepare that state, apply a Jacobian, repair a direction, or solve any equations.\n" + << "The existing WithState wrapper restores the accepted context after this callback returns;\n" + << "that standard restoration is outside the helper counters and timings reported here.\n" + << "Bounds are polynomial enclosures, not measured minimum/maximum values. Witness coordinates\n" + << "identify a proved negative point when status is negative; unresolved is never a certificate.\n"; + if (!(sameAcceptedPhysical && sameAcceptedNormalized && sameAcceptedResidual && sameCorrection && + sameTrialState && sameTrialResidual && sameCandidate && samePreparation)) + throw std::runtime_error("Polynomial endpoint certification changed a context buffer or preparation."); +} +}