486 lines
28 KiB
Python
486 lines
28 KiB
Python
#!/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())
|