From 0b0ab9cd1836da75553673149fd15da4eb33d36f Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Fri, 18 Sep 2026 08:47:19 -0400 Subject: [PATCH] Plot audited nonlinear paths and matched direction controls --- experiments/plot_nonlinear_step.py | 485 +++++++++++++++++++++++++++++ 1 file changed, 485 insertions(+) create mode 100644 experiments/plot_nonlinear_step.py diff --git a/experiments/plot_nonlinear_step.py b/experiments/plot_nonlinear_step.py new file mode 100644 index 0000000..2e3a7ab --- /dev/null +++ b/experiments/plot_nonlinear_step.py @@ -0,0 +1,485 @@ +#!/usr/bin/env python3 +"""Plot independently audited frozen nonlinear-step campaigns. + +Example (no FEM solves or checker reruns): + /opt/homebrew/anaconda3/envs/general/bin/python experiments/plot_nonlinear_step.py \ + --audited-campaign /path/late_v1 /path/late_review.json \ + --audited-campaign /path/controls_v1 /path/controls_review.json \ + --output /path/figures + +The review must come from check_nonlinear_step.py and pass for every supplied +case. Figures and extracted JSON are written outside the immutable campaigns. +Raw denotes the affine trial; curvature denotes its gravity-only correction. +Neither a passing artifact audit nor an improving sampled step validates a +physical stellar model. +""" + +import argparse +import csv +from dataclasses import dataclass +import hashlib +import json +import math +import os +from pathlib import Path +import tempfile + +import numpy as np + + +PREFIX = "nonlinear_step_" +PATHS = ("raw", "curvature") +BASE_IDENTITY_FILES = ( + "physical_base.txt", "normalized_base.txt", "residual.txt", "rhs.txt", + "state_factors.txt", "residual_factors.txt", "layout.csv", +) +COLORS = ("#2563a6", "#ca6c24", "#278064", "#9a438c", "#6d65a7") +STYLES = {"raw": "--", "curvature": "-"} +PATH_LABELS = {"raw": "Affine", "curvature": "Gravity corrected"} +STATUS_LABELS = { + "usable": ("o", "#54856b", "Usable"), + "raw_preparation_failed": ("x", "#b42c3b", "Preparation failed"), + "raw_preparation_blocked": ("x", "#b42c3b", "Preparation blocked correction"), + "inner_solve_failed": ("D", "#c07c11", "Inner solve failed"), + "corrected_preparation_failed": ("X", "#b42c3b", "Corrected preparation failed"), + "prepared_gravity_gate_failed": ("s", "#81549c", "Prepared gravity check failed"), +} + + +def sha256(path): + result = hashlib.sha256() + with Path(path).open("rb") as stream: + for chunk in iter(lambda: stream.read(1024 * 1024), b""): + result.update(chunk) + return result.hexdigest() + + +def read_rows(path): + with path.open(newline="") as stream: + return list(csv.DictReader(stream)) + + +def read_vector(path): + with path.open() as stream: + size = int(stream.readline()) + result = np.loadtxt(stream, ndmin=1) + if result.shape != (size,) or size <= 0 or not np.isfinite(result).all(): + raise ValueError(f"Invalid saved vector: {path}") + return result + + +def finite(value, name): + result = float(value) + if not math.isfinite(result): + raise ValueError(f"Expected a finite {name}: {value}") + return result + + +def close(left, right): + return math.isclose(left, right, rel_tol=2e-11, abs_tol=1e-28) + + +def clean(value): + """Unavailable observations are JSON null, never zero or nonstandard NaN.""" + if isinstance(value, dict): + return {key: clean(item) for key, item in value.items()} + if isinstance(value, (list, tuple)): + return [clean(item) for item in value] + if isinstance(value, (np.integer, np.floating)): + value = value.item() + if isinstance(value, float) and not math.isfinite(value): + return None + return value + + +@dataclass +class Case: + identifier: str + campaign: Path + folder: Path + review_path: Path + review: dict + metadata: dict + summary: dict + paths: list + hashes: dict + linear: dict + orders: dict + + @property + def eta(self): + return finite(self.summary["linear_relative_tolerance"], "forcing tolerance") + + @property + def is_seed(self): + return self.metadata.get("nonlinear_step_input") == "lane_emden_seed" + + @property + def identity(self): + return tuple(self.hashes[PREFIX + name] for name in BASE_IDENTITY_FILES) + + @property + def base_norm(self): + return finite(self.summary["base_residual_norm"], "base residual norm") + + @property + def label(self): + suffix = "" if self.review["direction_quality"]["passed"] else " (linear gate failed)" + return rf"$\eta={self.eta:g}$" + suffix + + +def load_campaign(campaign, review_path): + campaign, review_path = campaign.resolve(), review_path.resolve() + review = json.loads(review_path.read_text()) + provenance_path = campaign / "provenance.json" + provenance = json.loads(provenance_path.read_text()) + if review.get("passed") is not True: + raise ValueError(f"Audit did not pass: {review_path}") + if review.get("provenance_sha256") != sha256(provenance_path): + raise ValueError(f"Audit belongs to different provenance: {campaign}") + if review.get("source_commit") != provenance.get("commit"): + raise ValueError(f"Source commit does not match audit: {campaign}") + case_ids = [case["id"] for case in provenance["cases"]] + if not case_ids or len(set(case_ids)) != len(case_ids) or set(review.get("cases", {})) != set(case_ids): + raise ValueError(f"Campaign and audit case lists disagree: {campaign}") + if json.loads((campaign / "cases.json").read_text()) != provenance["cases"]: + raise ValueError(f"Case snapshot differs from provenance: {campaign}") + cases = [] + for identifier in case_ids: + audited = review["cases"][identifier] + if audited.get("passed") is not True or audited.get("observations_valid") is not True: + raise ValueError(f"Case lacks a passing vector audit: {identifier}") + folder = campaign / identifier + record = json.loads((campaign / f"{identifier}_run.json").read_text()) + if record.get("exit_code") != 0 or record.get("timed_out") or not record.get("finished_utc"): + raise ValueError(f"Case is not complete: {identifier}") + recorded_hashes = record["output_hashes"] + # Validate every artifact required by the passed audit, not just the + # small CSV subset plotted below. The campaign itself remains read-only. + for name in audited["required_artifacts"]: + artifact = folder / name + if artifact.resolve().parent != folder.resolve(): + raise ValueError(f"Unexpected artifact path: {name}") + if name not in recorded_hashes or sha256(artifact) != recorded_hashes[name]: + raise ValueError(f"Saved artifact hash changed: {artifact}") + metadata = dict(line.split("=", 1) for line in (folder / "metadata.txt").read_text().splitlines() if "=" in line) + if metadata != audited["metadata"]: + raise ValueError(f"Metadata no longer matches the audit: {identifier}") + summary = {row["metric"]: float(row["value"]) for row in read_rows(folder / (PREFIX + "summary.csv"))} + if clean(summary) != audited["summary"]: + raise ValueError(f"Summary no longer matches the audit: {identifier}") + if metadata.get("nonlinear_step_paths") != "raw,curvature" or metadata.get("nonlinear_step_positivity_repair") != "0": + raise ValueError(f"Unsupported trial-path semantics: {identifier}") + audited_paths = {(row["sample"], row["path"]): row for row in audited["paths"] if row["direction"] == "full"} + paths = [] + maximum_drift = finite(summary["maximum_original_restoration_error"], "restoration drift") + epsilon = finite(summary["machine_epsilon"], "machine epsilon") + for row in read_rows(folder / (PREFIX + "paths.csv")): + if row["direction"] != "full": + continue + sample, path = int(row["sample"]), row["path"] + audited_row = audited_paths[(sample, path)] + flags = {key: row[key] == "1" for key in ("prepared", "usable", "increment_resolved", "inner_converged", "prepared_gate")} + if any(flags[key] != audited_row[key] for key in flags): + raise ValueError(f"Path status differs from audit: {identifier}/{sample}/{path}") + alpha = finite(row["alpha"], "alpha") + if alpha != audited_row["alpha"] or path not in PATHS: + raise ValueError(f"Path identity differs from audit: {identifier}") + point = {"sample": sample, "alpha": alpha, "path": path, **flags, "status": audited_row["status"]} + if flags["prepared"]: + for name in ("base_norm", "action_norm", "residual_norm", "linear_prediction_norm", "remainder_norm", + "remainder_over_alpha", "relative_derivative_error", "slope", "heuristic_roundoff_scale", + "linear_dot_remainder"): + point[name] = finite(row[name], name) + if not close(point["residual_norm"], audited_row["actual_norm"]) or not close(point["linear_prediction_norm"], audited_row["predicted_norm"]): + raise ValueError(f"Residual norm differs from audit: {identifier}/{sample}/{path}") + # Match the independent checker's final restoration envelope; + # corrected-path preparation has a separate arithmetic estimate. + floor = max(8 * maximum_drift, 128 * epsilon * ( + 1 + point["base_norm"] + alpha * point["action_norm"] + point["residual_norm"])) + if path == "curvature": + floor += finite(row["gravity_arithmetic_allowance"], "gravity arithmetic allowance") + for name in ("gravity_correction_norm", "gravity_target_defect", "gravity_affine_disagreement", + "gravity_arithmetic_allowance", "gravity_target_threshold"): + point[name] = finite(row[name], name) + point["floor"] = floor + point["resolved_for_order"] = flags["usable"] and flags["increment_resolved"] and point["remainder_norm"] > 10 * floor + point["residual_ratio"] = point["residual_norm"] / point["base_norm"] + point["residual_change_percent"] = 100 * (point["residual_ratio"] - 1) + point["merit_change_ratio"] = point["residual_ratio"] ** 2 - 1 + direction_norm = finite(summary["full_direction_norm"], "direction norm") + point["remainder_over_alpha_squared_direction_norm_squared"] = ( + point["remainder_norm"] / (alpha * direction_norm) ** 2 if direction_norm > 0 else None) + point["signed_remainder_merit_term_over_base_merit"] = 2 * point["linear_dot_remainder"] / point["base_norm"] ** 2 + point["squared_remainder_merit_term_over_base_merit"] = (point["remainder_norm"] / point["base_norm"]) ** 2 + point["observed_decrease"] = flags["usable"] and point["residual_ratio"] < 1 + point["satisfies_default_unrepaired_merit_test"] = flags["usable"] and point["residual_ratio"] ** 2 <= 1 - 1e-4 * alpha + else: + point.update(residual_ratio=None, residual_change_percent=None, remainder_over_alpha=None, + relative_derivative_error=None, resolved_for_order=False, observed_decrease=False) + paths.append(point) + if len(paths) != 16: + raise ValueError(f"Expected eight full-direction samples for both paths: {identifier}") + orders = {row["path"]: row["resolved_intervals"] for row in audited["vector_order_observations"] + if row["direction"] == "full" and row["block"] == "all"} + cases.append(Case(identifier, campaign, folder, review_path, audited, metadata, summary, paths, + recorded_hashes, read_rows(folder / (PREFIX + "linear.csv"))[0], orders)) + return cases + + +def compare_cases(left, right): + """Only equal frozen equations, coordinates and RHS support an eta control.""" + identical = {name: left.hashes[PREFIX + name] == right.hashes[PREFIX + name] for name in BASE_IDENTITY_FILES} + result = {"left": left.identifier, "right": right.identifier, "same_input_hashes": identical, + "available": all(identical.values())} + if not result["available"]: + return result + p, q = [read_vector(case.folder / (PREFIX + "direction_full.txt")) for case in (left, right)] + jp, jq = [read_vector(case.folder / (PREFIX + "action_full.txt")) for case in (left, right)] + difference, image = q - p, jq - jp + pn, qn, wn, jwn = (float(np.linalg.norm(vector)) for vector in (p, q, difference, image)) + result.update(direction_difference_norm=wn, action_difference_norm=jwn, + relative_direction_difference=wn / pn if pn else None, + relative_action_difference_to_base_residual=jwn / left.base_norm, + sampled_action_over_direction_difference=jwn / wn if wn else None, + direction_cosine=float(np.dot(p, q)) / (pn * qn) if pn * qn else None, + both_newton_linear_gates_pass=all(case.review["direction_quality"]["passed"] for case in (left, right)), + scope="A sampled weak response is not a nullspace proof, gauge identification or condition estimate.") + result["difference_by_block"] = {"value": {}, "residual": {}} + for row in read_rows(left.folder / (PREFIX + "layout.csv")): + values = difference if row["kind"] == "value" else image + start, size = int(row["offset"]), int(row["size"]) + result["difference_by_block"][row["kind"]][row["block"]] = float(np.linalg.norm(values[start:start + size])) + left_points = {(point["alpha"], point["path"]): point for point in left.paths} + result["paired_full_paths"] = [] + for point in right.paths: + earlier = left_points.get((point["alpha"], point["path"])) + if earlier is None: + continue + result["paired_full_paths"].append({ + "alpha": point["alpha"], "path": point["path"], "left_status": earlier["status"], "right_status": point["status"], + "left_residual_ratio": earlier["residual_ratio"], "right_residual_ratio": point["residual_ratio"], + "both_usable": earlier["usable"] and point["usable"], + }) + return result + + +def extract_case(case): + best, largest_improving = {}, {} + for path in PATHS: + usable = [point for point in case.paths if point["path"] == path and point["usable"]] + best[path] = min(usable, key=lambda point: point["residual_ratio"]) if usable else None + improving = [point["alpha"] for point in usable if point["observed_decrease"]] + largest_improving[path] = max(improving) if improving else None + by_sample = {(point["sample"], point["path"]): point for point in case.paths} + correction_comparisons = [] + for sample in sorted({point["sample"] for point in case.paths}): + raw, corrected = (by_sample[(sample, path)] for path in PATHS) + if not (raw["usable"] and corrected["usable"]): + continue + correction_comparisons.append({ + "sample": sample, "alpha": raw["alpha"], + "raw_remainder_over_corrected_remainder": raw["remainder_norm"] / corrected["remainder_norm"] if corrected["remainder_norm"] else None, + "raw_signed_prediction_dot_remainder": raw["linear_dot_remainder"], + "corrected_signed_prediction_dot_remainder": corrected["linear_dot_remainder"], + "corrected_minus_raw_residual_ratio": corrected["residual_ratio"] - raw["residual_ratio"], + }) + return { + "campaign": str(case.campaign), "review": str(case.review_path), "review_sha256": sha256(case.review_path), + "state_kind": "seed" if case.is_seed else "saved_state", "forcing_tolerance": case.eta, + "preconditioner": case.metadata["context_preconditioner_variant"], "base_residual_norm": case.base_norm, + "direction_quality": case.review["direction_quality"], "linear_iterations": int(case.linear["iterations"]), + "linear_seconds": finite(case.linear["seconds"], "linear seconds"), + "direction_norm": finite(case.summary["full_direction_norm"], "direction norm"), + "merit_slope_over_squared_base_norm": float(case.summary["full_merit_slope"]) / case.base_norm ** 2, + "total_probe_seconds": finite(case.summary["total_seconds"], "probe seconds"), + "full_paths": case.paths, "best_sampled_usable_trial": best, + "corrected_versus_raw": correction_comparisons, + "largest_sampled_improving_alpha": largest_improving, + "audited_vector_orders": case.orders, + "shape_fields_action_cancellation": case.review["shape_fields_action_cancellation"], + "scopes": [ + "Every trial uses a fixed sampled alpha; no line search or model update was performed.", + "Merit-test labels use the default unmodified-direction criterion, not a claim of production acceptance.", + "Vector-order intervals come only from the passing independent audit; heuristic-floor points are excluded.", + "Directional Taylor consistency does not require a passing outer solve; Newton forcing comparisons do.", + "Remainder/(alpha squared times direction norm squared) describes this sampled trial path; it is not a universal Hessian norm.", + "A smaller remainder can produce worse merit when its signed projection against the linear prediction is more adverse.", + ], + } + + +def group_cases(cases): + groups = {} + for case in cases: + groups.setdefault((case.is_seed, case.identity), []).append(case) + # Seed and late state deliberately have independent panels and axis limits. + ordered = sorted(groups.values(), key=lambda group: (group[0].is_seed, group[0].identifier)) + for group in ordered: + group.sort(key=lambda case: (-case.eta, case.identifier)) + return ordered + + +def draw_figure(cases, output, small_step_maximum): + import matplotlib + matplotlib.use("Agg") + from matplotlib import pyplot as plt + from matplotlib.lines import Line2D + + plt.rcParams.update({"font.size": 10, "axes.titlesize": 12, "axes.labelsize": 10, + "pdf.fonttype": 42, "ps.fonttype": 42, "savefig.dpi": 220}) + groups = group_cases(cases) + figure, axes = plt.subplots(4, len(groups), figsize=(6.3 * len(groups), 11.3), squeeze=False, + gridspec_kw={"height_ratios": [1.7, 1.45, 1.7, .7]}) + status_seen = set() + for column, group in enumerate(groups): + first = group[0] + ratio_ax, change_ax, derivative_ax, status_ax = axes[:, column] + state_label = "Seed" if first.is_seed else "Late saved state" + ratio_ax.set_title(f"{state_label} | $\\|F_0\\|={first.base_norm:.5g}$") + ratio_ax.axhline(1., color="#555555", lw=.9, zorder=0) + change_ax.axhline(0., color="#555555", lw=.9, zorder=0) + statuses = [] + for index, case in enumerate(group): + color = COLORS[index % len(COLORS)] + for path in PATHS: + points = sorted((point for point in case.paths if point["path"] == path), key=lambda point: point["alpha"]) + label = f"{case.label} · {PATH_LABELS[path]}" + xs = np.array([point["alpha"] for point in points]) + ys = np.array([point["residual_ratio"] if point["usable"] else np.nan for point in points], dtype=float) + ratio_ax.plot(xs, ys, STYLES[path], color=color, marker="o", ms=4, lw=1.5, label=label) + selected = xs <= small_step_maximum * (1 + 1e-12) + changes = np.array([point["residual_change_percent"] if point["usable"] else np.nan for point in points], dtype=float) + change_ax.plot(xs[selected], changes[selected], STYLES[path], color=color, marker="o", ms=4, lw=1.5) + positive = {point["alpha"]: point for point in points if point["usable"] and point["remainder_over_alpha"] > 0} + # Draw interval segments only when the independently audited + # vector evidence supports an order observation. No fitted + # asymptotic slope is extrapolated through failed or noisy data. + for interval in case.orders.get(path, []): + pair = [positive.get(interval[key]) for key in ("alpha_small", "alpha_large")] + if all(point is not None and point["resolved_for_order"] for point in pair): + derivative_ax.plot([point["alpha"] for point in pair], [point["remainder_over_alpha"] for point in pair], + STYLES[path], color=color, lw=1.5) + for point in positive.values(): + derivative_ax.plot(point["alpha"], point["remainder_over_alpha"], "o" if path == "raw" else "s", ms=4.7, + color=color, markerfacecolor=color if point["resolved_for_order"] else "none") + row_index = len(statuses) + statuses.append(label) + for point in points: + marker, status_color, _ = STATUS_LABELS[point["status"]] + status_seen.add(point["status"]) + status_ax.scatter(point["alpha"], row_index, marker=marker, c=status_color, + s=21 if point["usable"] else 34, linewidths=1.1) + failed_prepared = [point for point in points if point["prepared"] and not point["usable"]] + if failed_prepared: + ratio_ax.scatter([point["alpha"] for point in failed_prepared], + [point["residual_ratio"] for point in failed_prepared], + marker="s", facecolors="none", edgecolors="#81549c", s=38) + ratio_ax.set_ylabel(r"$\|F(\mathrm{trial})\|/\|F_0\|$") + ratio_ax.set_yscale("log") + ratio_ax.legend(loc="best", fontsize=8, framealpha=.92) + change_ax.set_ylabel("Residual change (%)") + change_ax.set_title(rf"Small steps only: $\alpha\leq {small_step_maximum:g}$", fontsize=10) + change_ax.ticklabel_format(axis="y", style="plain", useOffset=False) + derivative_ax.set_ylabel(r"$\|F(\mathrm{trial})-F_0-\alpha Jp\|/\alpha$") + derivative_ax.set_yscale("log") + derivative_ax.set_title("Vector derivative error; filled points resolve the heuristic floor", fontsize=9) + status_ax.set_yticks(range(len(statuses)), statuses, fontsize=8) + status_ax.set_ylim(-.6, len(statuses) - .4) + status_ax.set_ylabel("Trial status", fontsize=9) + for axis in (ratio_ax, change_ax, derivative_ax, status_ax): + axis.set_xscale("log") + axis.set_xlabel(r"Step length $\alpha$") + axis.grid(True, which="major", color="#dddddd", lw=.65) + axis.set_axisbelow(True) + all_alphas = [point["alpha"] for case in group for point in case.paths] + for axis in (ratio_ax, derivative_ax, status_ax): + axis.set_xlim(min(all_alphas) / 1.3, max(all_alphas) * 1.3) + small = [alpha for alpha in all_alphas if alpha <= small_step_maximum] + if small: + change_ax.set_xlim(min(small) / 1.2, max(small) * 1.2) + if not any(point["usable"] and point["alpha"] <= small_step_maximum for case in group for point in case.paths): + change_ax.text(.5, .5, "No usable small-step observations", transform=change_ax.transAxes, + ha="center", va="center", fontsize=9) + status_handles = [Line2D([0], [0], marker=STATUS_LABELS[key][0], color="none", + markeredgecolor=STATUS_LABELS[key][1], markerfacecolor=STATUS_LABELS[key][1], + label=STATUS_LABELS[key][2], markersize=5) + for key in STATUS_LABELS if key in status_seen] + figure.suptitle("Frozen nonlinear steps\nAffine and gravity-corrected trial paths", fontsize=14, y=.997) + figure.legend(handles=status_handles, loc="lower center", bbox_to_anchor=(.5, .061), ncol=min(4, len(status_handles)), + frameon=False, fontsize=8) + figure.text(.5, .019, + "All norms use fixed diagonal normalization. Missing trials appear only in the status strip; they are not zero residuals.\n" + "Derivative-error segments use audited, usable, resolved vector observations above 10× the heuristic floor.\n" + "A smaller vector remainder need not improve merit: its signed projection against the prediction also matters.\n" + "Cold linear directions and sampled trials do not establish nonlinear convergence, virial balance or pointwise physical validity.", + ha="center", va="bottom", fontsize=8, color="#404040") + figure.tight_layout(rect=(0., .115, 1., .965), h_pad=1.1, w_pad=2.2) + for extension in ("png", "pdf"): + figure.savefig(output / f"nonlinear_step_comparison.{extension}", bbox_inches="tight") + plt.close(figure) + + +def main(argv=None): + parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument("--audited-campaign", nargs=2, action="append", required=True, metavar=("CAMPAIGN", "REVIEW_JSON")) + parser.add_argument("--output", type=Path, required=True) + parser.add_argument("--small-step-maximum", type=float, default=1 / 1024, + help="Largest alpha in the separate percent-change panel (default 1/1024).") + args = parser.parse_args(argv) + if not math.isfinite(args.small_step_maximum) or args.small_step_maximum <= 0: + parser.error("--small-step-maximum must be finite and positive") + output = args.output.resolve() + campaigns = [Path(campaign).resolve() for campaign, _ in args.audited_campaign] + if any(output == campaign or campaign in output.parents for campaign in campaigns): + parser.error("The output directory must be outside every immutable campaign") + if output.exists(): + parser.error("The output directory already exists; choose a new directory") + try: + cases = [] + for campaign, audit in args.audited_campaign: + cases.extend(load_campaign(Path(campaign), Path(audit))) + if len({case.identifier for case in cases}) != len(cases): + raise ValueError("Case IDs must be unique across the supplied campaigns") + if any(case.base_norm <= 0 for case in cases): + raise ValueError("Residual ratios require a positive base norm") + comparisons = [compare_cases(left, right) for i, left in enumerate(cases) for right in cases[i + 1:] + if left.identity == right.identity] + report = { + "schema_version": 1, "plot_script_sha256": sha256(Path(__file__)), + "audits_required_to_pass": True, "small_step_maximum": args.small_step_maximum, + "campaigns": [{"path": str(Path(campaign).resolve()), "review": str(Path(audit).resolve()), + "review_sha256": sha256(Path(audit))} for campaign, audit in args.audited_campaign], + "cases": {case.identifier: extract_case(case) for case in cases}, + "same_state_comparisons": comparisons, + "scope": "Scientific extraction of independently audited frozen data; no model convergence or physical-validation certificate.", + "order_scope": "Audited intervals require usable paths, resolved additions and remainder above10 times the heuristic floor. No order-two requirement is imposed; a correct n=1.5 EOS tangent can have nonquadratic remainder near zero enthalpy.", + "comparison_scope": "Same-state direction comparisons require byte-identical saved base/RHS/factors/layout. Their sampled J response is not a condition estimate.", + } + output.mkdir(parents=True) + # Font discovery needs a cache, but never writes into the source tree, + # immutable campaigns or the user's default plotting configuration. + with tempfile.TemporaryDirectory(prefix="meanfield_plot_config_") as cache: + cache_variables = ("MPLCONFIGDIR", "XDG_CACHE_HOME") + old_caches = {name: os.environ.get(name) for name in cache_variables} + for name in cache_variables: + os.environ[name] = cache + try: + draw_figure(cases, output, args.small_step_maximum) + finally: + for name, previous in old_caches.items(): + if previous is None: + os.environ.pop(name, None) + else: + os.environ[name] = previous + (output / "nonlinear_step_comparison.json").write_text(json.dumps(clean(report), indent=2, allow_nan=False) + "\n") + except (OSError, KeyError, TypeError, ValueError) as error: + parser.exit(1, f"Nonlinear-step plot failed: {error}\n") + print(output / "nonlinear_step_comparison.png") + print(output / "nonlinear_step_comparison.pdf") + print(output / "nonlinear_step_comparison.json") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main())