From a458af0a7a14b03368b80a27ae208928e2e07250 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Sat, 19 Sep 2026 14:57:43 -0400 Subject: [PATCH] Plot audited paired Newton convergence with actual thresholds --- experiments/plot_coupled_convergence.py | 237 +++++++++++++++++++ experiments/test_plot_coupled_convergence.py | 147 ++++++++++++ 2 files changed, 384 insertions(+) create mode 100644 experiments/plot_coupled_convergence.py create mode 100644 experiments/test_plot_coupled_convergence.py diff --git a/experiments/plot_coupled_convergence.py b/experiments/plot_coupled_convergence.py new file mode 100644 index 0000000..da8b7b9 --- /dev/null +++ b/experiments/plot_coupled_convergence.py @@ -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() diff --git a/experiments/test_plot_coupled_convergence.py b/experiments/test_plot_coupled_convergence.py new file mode 100644 index 0000000..2e4af56 --- /dev/null +++ b/experiments/test_plot_coupled_convergence.py @@ -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()