tboudreaux/MeanField
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Plot audited paired Newton convergence with actual thresholds

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-19 14:57:43 -04:00
parent d62487f2d1
commit a458af0a7a
2 changed files with 384 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

237
experiments/plot_coupled_convergence.py Normal file
View File

@@ -0,0 +1,237 @@
#!/usr/bin/env python3
"""Plot completed, artifact-passed coupled-convergence audit summaries.
Example (run rendering serially with numerical campaigns):
python plot_coupled_convergence.py analysis.json --output NEW_FIGURE_DIRECTORY
Use --pair n3 or --pair n2 to request exactly one complete pair. Default: both.
The script reads summary JSON only; it does not rerun the underlying FEM/artifact
audit. Matplotlib is imported only when figures are actually rendered.
"""
from __future__ import annotations
import argparse
from dataclasses import dataclass
import hashlib
import json
import math
from pathlib import Path
import sys
PAIRS = {
"n3": ("n3_J0_control", "n3_J0_tight"),
"n2": ("n2_restart_raw", "n2_restart_gravity_path"),
}
LABELS = {
"n3_J0_control": "Control linear tolerance",
"n3_J0_tight": "Tighter linear tolerance",
"n2_restart_raw": "Raw trial path",
"n2_restart_gravity_path": "Gravity-completed trial path",
}
def require(condition, message):
if not condition:
raise ValueError(message)
def finite(value, name):
number = float(value)
require(math.isfinite(number), f"Nonfinite {name}")
return number
def near(a, b, name):
require(math.isclose(a, b, rel_tol=3e-11, abs_tol=0.), f"Inconsistent {name}: {a} versus {b}")
def digest_bytes(data):
return hashlib.sha256(data).hexdigest()
@dataclass
class Curve:
identifier: str
campaign: str
audit_path: str
producer_commit: str
mesh_sha256: str
metadata: dict
input_state_hashes: dict
time: list[float]
residual: list[float]
accepted: list[bool]
threshold: float
outcome: str
failure: str | None
def curve(case, campaign, audit_path):
name = case["id"]
require(case.get("complete") is True and case.get("artifact_passed") is True,
f"Case {name} is incomplete or did not pass the artifact audit")
require(case.get("kind") == "newton" and case.get("outcome") in ("converged", "not_converged"),
f"Case {name} has no completed Newton trajectory")
initial = finite(case["initial_residual"], "initial residual")
require(initial > 0, "Log plots require strictly positive residuals; exact zeros are not silently floored")
atol, rtol = finite(case["atol"], "atol"), finite(case["rtol"], "rtol")
require(atol >= 0 and rtol >= 0, "Nonlinear tolerances must be nonnegative")
threshold = finite(case["actual_threshold"], "actual stopping threshold")
near(threshold, max(atol, rtol * initial), "actual stopping threshold")
require(threshold > 0, "A zero stopping threshold cannot be drawn on an unmodified log scale")
history = case["history"]
require(isinstance(history, list) and len(history) == int(case["attempted_updates"]), "Missing Newton history rows")
times, residuals, accepted, durations = [0.], [initial], [], []
for index, row in enumerate(history):
require(int(row["iteration"]) == index, "Missing, repeated or unordered Newton iteration")
require(str(row["accepted"]) in ("0", "1"), "Invalid acceptance flag")
seconds = finite(row["seconds"], "Newton duration")
residual = finite(row["F"], "Newton residual")
require(seconds >= 0 and residual > 0, "Negative time or nonpositive log residual")
is_accepted = str(row["accepted"]) == "1"
if not is_accepted:
near(residual, residuals[-1], "restored residual after failed Newton attempt")
durations.append(seconds)
times.append(math.fsum(durations))
residuals.append(residual)
accepted.append(is_accepted)
near(times[-1], finite(case["newton_seconds"], "reported Newton time"), "summed Newton time")
near(residuals[-1], finite(case["final_residual"], "final residual"), "final residual")
require(sum(accepted) == int(case["accepted_updates"]), "Accepted-update count differs from history")
require(case["outcome"] != "converged" or residuals[-1] <= threshold,
"Convergence outcome exceeds the actual stopping threshold")
return Curve(name, campaign["path"], audit_path, campaign["producer_commit"], campaign["mesh_sha256"],
case["metadata"], case.get("run", {}).get("input_state_hashes", {}),
times, residuals, accepted, threshold, case["outcome"], case.get("failure"))
def load_pairs(paths, selected="both"):
require(selected in ("both", "n3", "n2"), "Unknown requested curve pair")
requested = list(PAIRS) if selected == "both" else [selected]
wanted = {name for pair in requested for name in PAIRS[pair]}
curves, sources = {}, []
for path in paths:
path = Path(path).resolve()
raw = path.read_bytes()
audit = json.loads(raw)
require(audit.get("complete") is True and audit.get("artifact_passed") is True,
f"Only completed, artifact-passed audit summaries may be plotted: {path}")
sources.append({"path": str(path), "sha256": digest_bytes(raw),
"auditor_sha256": audit["auditor_sha256"],
"benchmark_checker_sha256": audit.get("benchmark_checker_sha256"),
"audit_scope": audit.get("scope"), "audit_limitations": audit.get("limitations", [])})
for campaign in audit["campaigns"]:
require(campaign.get("complete") is True and campaign.get("artifact_passed") is True,
"An audit contains an incomplete or failed campaign")
for case in campaign["cases"]:
if case.get("id") in wanted:
require(case["id"] not in curves, f"Ambiguous duplicate case {case['id']} across audits")
curves[case["id"]] = curve(case, campaign, str(path))
require(set(curves) == wanted, f"Missing requested curves: {sorted(wanted - set(curves))}")
pairs = {}
for name in requested:
left, right = (curves[key] for key in PAIRS[name])
require(left.mesh_sha256 == right.mesh_sha256, "A paired comparison uses different mesh hashes")
for key in ("n", "J", "K", "rho_c", "surface_prescription", "surface_degree", "surface_maximum_order",
"density_order", "enthalpy_order", "gravity_potential_order", "gravity_flux_order", "displacement_order"):
require(left.metadata.get(key) == right.metadata.get(key), f"Paired physical/discrete setting differs: {key}")
require(float(left.metadata["n"]) == (3. if name == "n3" else 2.) and float(left.metadata["J"]) == 0.,
"A named control pair has unexpected physical parameters")
if name == "n2":
require(bool(left.input_state_hashes) and left.input_state_hashes == right.input_state_hashes,
"Restart pair does not declare the same hash-bound input state files")
pairs[name] = (left, right)
return pairs, sources
def render_pair(name, curves):
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
figure, axis = plt.subplots(figsize=(8.5, 5.7), layout="constrained")
colors = ("#1769a8", "#c35b16")
markers = ("o", "s")
for candidate, color, marker in zip(curves, colors, markers):
label = LABELS[candidate.identifier]
tolerance = candidate.metadata.get("linear_tol")
if tolerance is not None:
label += f" (linear rtol {float(tolerance):g})"
label += "; " + ("converged" if candidate.outcome == "converged" else "not converged")
axis.plot(candidate.time, candidate.residual, color=color, marker=marker, markersize=4,
linewidth=1.6, label=label)
failed = [i + 1 for i, accepted in enumerate(candidate.accepted) if not accepted]
if failed:
axis.scatter([candidate.time[i] for i in failed], [candidate.residual[i] for i in failed],
marker="x", s=65, linewidths=1.6, color=color, zorder=4)
if curves[0].threshold == curves[1].threshold:
axis.axhline(curves[0].threshold, color="#555555", linestyle=":", linewidth=1.4,
label=f"Actual stopping threshold (both): {curves[0].threshold:.4g}")
else:
for candidate, color in zip(curves, colors):
axis.axhline(candidate.threshold, color=color, linestyle=":", linewidth=1.3,
label=f"{LABELS[candidate.identifier]} threshold: {candidate.threshold:.4g}")
axis.set_yscale("log")
axis.set_xlim(left=0., right=max(1., max(c.time[-1] for c in curves) * 1.04))
axis.set_xlabel("Cumulative reported Newton time (s)")
axis.set_ylabel("Normalized coupled residual")
axis.grid(True, which="major", alpha=.22)
axis.legend(fontsize=8, loc="best")
title = "n = 3, J = 0: linear-solve tolerance control" if name == "n3" else "n = 2, J = 0: common-snapshot trial-path control"
axis.set_title(title, pad=12)
figure.supxlabel("Includes the initial norm at t = 0 and every completed Newton attempt; x marks a failed attempt.\n"
"Time excludes setup/observers. Residual convergence does not establish unbordered equilibrium or virial balance.",
fontsize=8)
return figure, matplotlib.__version__
def sidecar(name, curves, sources):
return {
"pair": name, "source_audits": sources,
"plotter_sha256": digest_bytes(Path(__file__).read_bytes()),
"curves": [{"id": c.identifier, "campaign": c.campaign, "audit_path": c.audit_path,
"producer_commit": c.producer_commit, "mesh_sha256": c.mesh_sha256,
"metadata": c.metadata, "input_state_hashes": c.input_state_hashes,
"cumulative_newton_seconds": c.time, "normalized_residual": c.residual,
"accepted_updates_after_t0": c.accepted, "actual_stopping_threshold": c.threshold,
"outcome": c.outcome, "failure": c.failure} for c in curves],
"scope": ["Consumes the completed artifact-passed audit JSON; underlying file/FEM audits are not rerun.",
"Initial residual and stopping threshold are taken from the solver metric audited in that summary.",
"Cumulative time sums completed history rows; setup, observers and any pre-solver initialization are excluded.",
"Failed attempts retain their restored residual and elapsed time; no missing or incomplete case is drawn.",
"Zero residuals or thresholds are refused rather than silently floored on the logarithmic axis.",
"The n=3 norm may include the central-density border; neither residual decline nor its stopping threshold proves physical equilibrium or virial balance."],
}
def main():
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("audits", nargs="+", type=Path)
parser.add_argument("--output", required=True, type=Path, help="New output directory outside the source checkout")
parser.add_argument("--pair", choices=("both", "n3", "n2"), default="both")
args = parser.parse_args()
require(not args.output.exists(), "Refusing to overwrite an existing figure directory")
require(not args.output.resolve().is_relative_to(Path(__file__).resolve().parents[1]),
"Write plot artifacts outside the source checkout")
pairs, sources = load_pairs(args.audits, args.pair)
# Preflight every requested curve before importing Matplotlib or creating
# output, so missing/incomplete controls cannot silently produce half a set.
args.output.mkdir(parents=True)
for name, curves in pairs.items():
figure, version = render_pair(name, curves)
report = sidecar(name, curves, sources)
report["matplotlib_version"], report["python_version"] = version, sys.version
report["figures"] = []
for suffix in ("png", "pdf"):
path = args.output / f"{name}_paired_convergence.{suffix}"
figure.savefig(path, dpi=200)
report["figures"].append({"path": str(path.resolve()), "sha256": digest_bytes(path.read_bytes())})
from matplotlib import pyplot as plt
plt.close(figure)
path = args.output / f"{name}_paired_convergence.json"
path.write_text(json.dumps(report, indent=2, allow_nan=False) + "\n")
print(path)
if __name__ == "__main__":
main()

147
experiments/test_plot_coupled_convergence.py Normal file
View File

@@ -0,0 +1,147 @@
#!/usr/bin/env python3
"""Cheap summary-schema checks; opt in to rendering only during a serial interval.
python -m unittest -v test_plot_coupled_convergence.py
MEANFIELD_PLOT_RENDER_TESTS=1 python -m unittest -v test_plot_coupled_convergence.py
"""
import copy
import hashlib
import json
import os
from pathlib import Path
import tempfile
import unittest
import plot_coupled_convergence as plots
def case(name):
n = 3 if name.startswith("n3") else 2
values = [.5, .25, .25]
return {
"id": name, "complete": True, "artifact_passed": True, "kind": "newton", "outcome": "not_converged",
"failure": "Synthetic failed attempt", "initial_residual": 1., "final_residual": .25,
"rtol": 1e-8, "atol": 0., "actual_threshold": 1e-8,
"attempted_updates": 3, "accepted_updates": 2, "newton_seconds": 6.,
"metadata": {"n": str(n), "J": "0", "K": ".1", "rho_c": "1", "density_order": "3",
"enthalpy_order": "3", "surface_prescription": "harmonics", "surface_degree": "4", "linear_tol": ".03"},
"run": {"input_state_hashes": {"/synthetic/common_state.txt": "a" * 64} if n == 2 else {}},
"history": [{"iteration": str(i), "accepted": "1" if i < 2 else "0", "seconds": str(t), "F": str(f)}
for i, (t, f) in enumerate(zip((2., 3., 1.), values))],
}
def audit():
return {
"complete": True, "artifact_passed": True, "auditor_sha256": "b" * 64,
"benchmark_checker_sha256": "c" * 64, "scope": "synthetic schema control", "limitations": ["not real FEM"],
"campaigns": [{"path": "/synthetic/campaign", "producer_commit": "d" * 40, "mesh_sha256": "e" * 64,
"complete": True, "artifact_passed": True,
"cases": [case(name) for pair in plots.PAIRS.values() for name in pair]}],
}
class SchemaChecks(unittest.TestCase):
def setUp(self):
self.temporary = tempfile.TemporaryDirectory()
self.directory = Path(self.temporary.name)
self.path = self.directory / "audit.json"
self.data = audit()
self.write()
def tearDown(self):
self.temporary.cleanup()
def write(self):
self.path.write_text(json.dumps(self.data))
def test_t0_cumulative_time_threshold_and_failed_attempt_are_preserved(self):
pairs, sources = plots.load_pairs([self.path])
self.assertEqual(set(pairs), {"n3", "n2"})
curve = pairs["n3"][0]
self.assertEqual(curve.time, [0., 2., 5., 6.])
self.assertEqual(curve.residual, [1., .5, .25, .25])
self.assertEqual(curve.accepted, [True, True, False])
self.assertEqual(curve.threshold, 1e-8)
self.assertEqual(sources[0]["sha256"], hashlib.sha256(self.path.read_bytes()).hexdigest())
report = plots.sidecar("n3", pairs["n3"], sources)
self.assertEqual(report["curves"][0]["outcome"], "not_converged")
self.assertEqual(report["curves"][0]["actual_stopping_threshold"], 1e-8)
json.dumps(report, allow_nan=False)
def test_incomplete_audit_and_case_are_refused(self):
self.data["complete"] = False
self.write()
with self.assertRaisesRegex(ValueError, "completed"):
plots.load_pairs([self.path])
self.data = audit()
self.data["campaigns"][0]["cases"][0]["complete"] = False
self.write()
with self.assertRaisesRegex(ValueError, "incomplete"):
plots.load_pairs([self.path])
def test_missing_and_duplicate_requested_curves_are_refused(self):
self.data["campaigns"][0]["cases"].pop()
self.write()
with self.assertRaisesRegex(ValueError, "Missing requested"):
plots.load_pairs([self.path])
self.data = audit()
self.data["campaigns"][0]["cases"].append(copy.deepcopy(self.data["campaigns"][0]["cases"][0]))
self.write()
with self.assertRaisesRegex(ValueError, "duplicate"):
plots.load_pairs([self.path])
def test_explicit_single_pair_does_not_require_an_unrequested_pair(self):
self.data["campaigns"][0]["cases"] = self.data["campaigns"][0]["cases"][:2]
self.write()
pairs, _ = plots.load_pairs([self.path], "n3")
self.assertEqual(list(pairs), ["n3"])
def test_initial_dependent_actual_threshold_is_checked(self):
self.data["campaigns"][0]["cases"][0]["actual_threshold"] = 1e-6
self.write()
with self.assertRaisesRegex(ValueError, "stopping threshold"):
plots.load_pairs([self.path])
def test_restart_source_mismatch_is_refused(self):
self.data["campaigns"][0]["cases"][-1]["run"]["input_state_hashes"] = {"/different_state": "f" * 64}
self.write()
with self.assertRaisesRegex(ValueError, "hash-bound"):
plots.load_pairs([self.path])
def test_nonpositive_log_values_are_not_silently_floored(self):
self.data["campaigns"][0]["cases"][0]["history"][0]["F"] = "0"
self.write()
with self.assertRaisesRegex(ValueError, "nonpositive"):
plots.load_pairs([self.path])
def test_missing_history_and_falsely_converged_outcome_are_refused(self):
self.data["campaigns"][0]["cases"][0]["history"].pop()
self.write()
with self.assertRaisesRegex(ValueError, "Missing Newton"):
plots.load_pairs([self.path])
self.data = audit()
self.data["campaigns"][0]["cases"][0]["outcome"] = "converged"
self.write()
with self.assertRaisesRegex(ValueError, "exceeds"):
plots.load_pairs([self.path])
@unittest.skipUnless(os.environ.get("MEANFIELD_PLOT_RENDER_TESTS") == "1", "Rendering reserved for serial validation interval")
def test_render_has_log_axis_t0_and_actual_threshold(self):
pairs, _ = plots.load_pairs([self.path])
for name, curves in pairs.items():
figure, _ = plots.render_pair(name, curves)
axis = figure.axes[0]
self.assertEqual(axis.get_yscale(), "log")
self.assertEqual(axis.lines[0].get_xdata()[0], 0.)
self.assertEqual(list(axis.lines[-1].get_ydata()), [1e-8, 1e-8])
for extension in ("png", "pdf"):
path = self.directory / f"{name}.{extension}"
figure.savefig(path)
self.assertGreater(path.stat().st_size, 1000)
from matplotlib import pyplot as plt
plt.close(figure)
if __name__ == "__main__":
unittest.main()