Check frozen numerical strategies without conflating incomplete evidence with success

experiments/check_numerical_strategies.py

@@ -0,0 +1,315 @@
+#!/usr/bin/env python3
+"""Check frozen numerical-strategy probes, without certifying a stellar model.
+
+Exit 0: requested evidence passes. Exit 3: an invariant or expectation fails.
+Exit 2: evidence is incomplete (including a correctly reported failed solve).
+Use check_plateau_physics.py separately for an end-to-end equilibrium endpoint.
+"""
+
+import argparse
+import csv
+import json
+import math
+from pathlib import Path
+import sys
+
+
+EPSILON = sys.float_info.epsilon
+SIGN_STATUSES = {"certified_nonnegative", "unresolved", "negative", "non_finite"}
+
+
+def read_rows(path, required):
+    with path.open(newline="") as stream:
+        reader = csv.DictReader(stream)
+        if reader.fieldnames is None or len(set(reader.fieldnames)) != len(reader.fieldnames):
+            raise ValueError(f"Missing or duplicate columns in {path.name}")
+        if not set(required).issubset(reader.fieldnames):
+            raise ValueError(f"Missing required columns in {path.name}: {sorted(set(required) - set(reader.fieldnames))}")
+        result = []
+        for row in reader:
+            if None in row or any(value is None for value in row.values()):
+                raise ValueError(f"Malformed CSV row in {path.name}")
+            result.append(row)
+    return result
+
+
+def metrics(path):
+    result = {}
+    for row in read_rows(path, ("metric", "value")):
+        if not row["metric"] or row["metric"] in result:
+            raise ValueError(f"Empty or duplicate metric in {path.name}")
+        result[row["metric"]] = float(row["value"])
+    return result
+
+
+def number(row, key):
+    return float(row[key])
+
+
+def flag(row, key):
+    if row[key] not in ("0", "1"):
+        raise ValueError(f"Invalid boolean {key}: {row[key]!r}")
+    return row[key] == "1"
+
+
+def nonnegative(value):
+    return math.isfinite(value) and value >= 0.0
+
+
+def roundoff(*values):
+    return 512.0 * EPSILON * (1.0 + sum(abs(value) for value in values if math.isfinite(value)))
+
+
+def clean_json(value):
+    if isinstance(value, float) and not math.isfinite(value):
+        return None
+    if isinstance(value, dict):
+        return {key: clean_json(item) for key, item in value.items()}
+    if isinstance(value, list):
+        return [clean_json(item) for item in value]
+    return value
+
+
+class Assessment:
+    def __init__(self, folder, restoration_tolerance):
+        self.folder = folder
+        self.restoration_tolerance = restoration_tolerance
+        self.checks = []
+        self.outcomes = {}
+
+    def check(self, name, condition, **details):
+        self.checks.append({"name": name, "status": "pass" if condition else "fail", **details})
+
+    def incomplete(self, name, **details):
+        self.checks.append({"name": name, "status": "incomplete", **details})
+
+    def restore(self, values, accepted, restored):
+        for key in accepted:
+            self.check(key, values[key] == 0.0, value=values[key], required=0.0)
+        for key in restored:
+            self.check(key, nonnegative(values[key]) and values[key] <= self.restoration_tolerance,
+                       value=values[key], tolerance=self.restoration_tolerance)
+
+    def polynomial(self, expect_negative):
+        certificates = read_rows(self.folder / "polynomial_certificates.csv", ("state", "status", "witness_element", "x", "y", "z"))
+        by_state = {row["state"]: row for row in certificates}
+        if len(by_state) != len(certificates):
+            raise ValueError("Duplicate polynomial certificate state")
+        for name in ("frozen_base", "context_seed_anchor", "original_full_candidate"):
+            if name not in by_state:
+                raise ValueError(f"Missing polynomial certificate: {name}")
+        if any(row["status"] not in SIGN_STATUSES for row in certificates):
+            raise ValueError("Unknown polynomial certificate status")
+        repairs = read_rows(self.folder / "polynomial_repair.csv", ("status", "minimum_safe_fraction", "base_status"))
+        if len(repairs) != 1:
+            raise ValueError("Polynomial repair summary must contain one row")
+        repair = repairs[0]
+        trials = read_rows(self.folder / "polynomial_repair_trials.csv", ("fraction", "alpha", "certificate", "prepared", "merit_slope"))
+        base_status = by_state["frozen_base"]["status"]
+        self.outcomes["polynomial"] = {"base": base_status, "anchor": by_state["context_seed_anchor"]["status"],
+                                       "repair": repair["status"], "trials": len(trials)}
+        self.check("repair_base_matches_certificate", repair["base_status"] == base_status)
+        if expect_negative:
+            self.check("expected_negative_base_detected", base_status == "negative", observed=base_status)
+            self.check("negative_base_repair_rejected", repair["status"] == "invalid_base", observed=repair["status"])
+        if base_status != "certified_nonnegative":
+            self.check("uncertified_base_is_rejected", repair["status"] == "invalid_base")
+            self.check("uncertified_base_has_no_repaired_trials", all(number(row, "fraction") == 0.0 for row in trials))
+            if base_status == "unresolved":
+                self.incomplete("polynomial_base_classification", reason="Subdivision budget did not resolve the base sign")
+        if base_status == "negative":
+            witness = by_state["frozen_base"]
+            self.check("negative_base_has_local_witness", int(witness["witness_element"]) >= 0 and
+                       all(0.0 <= number(witness, axis) <= 1.0 for axis in ("x", "y", "z")))
+        if repair["status"] == "feasible":
+            fraction = number(repair, "minimum_safe_fraction")
+            self.check("repair_fraction_in_unit_interval", math.isfinite(fraction) and 0.0 <= fraction <= 1.0,
+                       fraction=fraction)
+            repaired = [row for row in trials if number(row, "fraction") > 0.0]
+            if fraction > 0.0 and not repaired:
+                self.incomplete("certified_repair_trials", reason="No repaired candidates were evaluated")
+            for index, row in enumerate(repaired):
+                self.check(f"repaired_polynomial_{index}_certified", row["certificate"] == "certified_nonnegative")
+            self.outcomes["polynomial"]["prepared_repaired_trials"] = sum(flag(row, "prepared") for row in repaired)
+            self.outcomes["polynomial"]["descent_repaired_trials"] = sum(number(row, "merit_slope") < 0.0 for row in repaired)
+        elif repair["status"] != "invalid_base":
+            self.incomplete("polynomial_repair_available", outcome=repair["status"])
+        restoration = metrics(self.folder / "polynomial_restoration.csv")
+        self.restore(restoration,
+                     ("accepted_physical_difference", "accepted_normalized_difference", "accepted_residual_difference", "restored_correction_difference"),
+                     ("restored_residual_difference",))
+
+    def gravity(self):
+        summary = metrics(self.folder / "numerical_strategy_summary.csv")
+        self.restore(summary, ("accepted_state_change", "accepted_residual_change", "accepted_physical_change"),
+                     ("maximum_original_restoration_error",))
+        inner = read_rows(self.folder / "numerical_strategy_inner.csv", ("kind", "alpha", "index", "status", "converged", "right_hand_side_norm", "true_residual_norm", "relative_true_residual_norm", "threshold"))
+        lookup = {}
+        for index, row in enumerate(inner):
+            key = (row["kind"], number(row, "alpha"), int(row["index"]))
+            if key in lookup:
+                raise ValueError(f"Duplicate inner solve record: {key}")
+            lookup[key] = row
+            converged = flag(row, "converged")
+            self.check(f"inner_{index}_status_consistent", (int(row["status"]) == 0) == converged)
+            if converged:
+                rhs, residual, threshold = (number(row, key) for key in ("right_hand_side_norm", "true_residual_norm", "threshold"))
+                self.check(f"inner_{index}_true_defect_gate", all(nonnegative(x) for x in (rhs, residual, threshold)) and residual <= math.nextafter(threshold, math.inf),
+                           true_residual=residual, threshold=threshold)
+                expected_relative = residual / rhs if rhs > 0 else (0.0 if residual == 0 else math.inf)
+                self.check(f"inner_{index}_relative_defect", math.isclose(number(row, "relative_true_residual_norm"), expected_relative, rel_tol=1e-12, abs_tol=1e-30))
+            else:
+                self.incomplete(f"inner_{index}_solve_available", kind=row["kind"], alpha=number(row, "alpha"),
+                                outcome="numerical_solve_failure", solve_status=int(row["status"]))
+        self.outcomes["gravity_inner"] = {"recorded_calls": len(inner), "failed_calls": sum(not flag(row, "converged") for row in inner)}
+        self.curvature(summary, lookup)
+        self.reduced(summary, lookup)
+
+    def curvature(self, summary, lookup):
+        rows = read_rows(self.folder / "numerical_strategy_curvature.csv", ("alpha", "prepared", "inner_converged", "corrected_prepared", "raw_gravity_remainder", "corrected_gravity_remainder", "gravity_correction_norm", "correction_over_alpha_squared", "retained_state_change"))
+        if len(rows) != 4 or sorted(number(row, "alpha") for row in rows) != [.001, .01, .1, 1.]:
+            self.incomplete("curvature_probe_coverage", expected_alphas=[1, .1, .01, .001], recorded_rows=len(rows))
+        reliable = []
+        evaluated = 0
+        for row in rows:
+            alpha = number(row, "alpha")
+            if not math.isfinite(alpha) or alpha <= 0:
+                raise ValueError("A curvature step must be finite and positive")
+            prepared, solved, corrected = (flag(row, key) for key in ("prepared", "inner_converged", "corrected_prepared"))
+            self.check(f"curvature_{alpha:g}_stages_consistent", (not solved or prepared) and (not corrected or solved))
+            if not (prepared and solved and corrected):
+                self.incomplete(f"curvature_{alpha:g}_corrected_trial_available", raw_prepared=prepared,
+                                inner_converged=solved, corrected_prepared=corrected)
+                continue
+            evaluated += 1
+            solve = lookup[("curvature", alpha, 0)]
+            self.check(f"curvature_{alpha:g}_inner_record_consistent", flag(solve, "converged"))
+            raw, actual, correction, coefficient = (number(row, key) for key in ("raw_gravity_remainder", "corrected_gravity_remainder", "gravity_correction_norm", "correction_over_alpha_squared"))
+            floor = roundoff(*(number(row, key) for key in ("base_residual_norm", "raw_residual_norm", "corrected_residual_norm", "linear_prediction_norm")))
+            self.check(f"curvature_{alpha:g}_finite_nonnegative", all(nonnegative(x) for x in (raw, actual, correction, coefficient)))
+            self.check(f"curvature_{alpha:g}_retained_unchanged", number(row, "retained_state_change") == 0.0)
+            self.check(f"curvature_{alpha:g}_coefficient_arithmetic", math.isclose(coefficient, correction / alpha**2, rel_tol=1e-12, abs_tol=1e-30))
+            tolerance = number(solve, "threshold") + floor
+            self.check(f"curvature_{alpha:g}_remainder_cancelled", nonnegative(actual) and actual <= tolerance,
+                       corrected_remainder=actual, allowed=tolerance, raw_remainder=raw)
+            if alpha <= .1 and raw > 100.0 * floor and correction > 100.0 * floor and number(solve, "threshold") < .01 * raw:
+                reliable.append((alpha, correction, coefficient))
+        self.outcomes["curvature"] = {"fully_evaluated_steps": evaluated, "unavailable_steps": len(rows) - evaluated,
+                                      "resolved_small_steps": reliable}
+        if evaluated == 0:
+            self.incomplete("curvature_cancellation_available", reason="No prepared, solved, corrected trial")
+        reliable.sort(reverse=True)
+        if len(reliable) >= 2:
+            coarse, fine = reliable[-2:]
+            order = math.log(fine[1] / coarse[1]) / math.log(fine[0] / coarse[0])
+            self.check("curvature_at_least_second_order_trend", order >= 1.7, observed_order=order,
+                       minimum_test_order=1.7, alphas=[coarse[0], fine[0]],
+                       interpretation="A finite-range trend check, not an asymptotic proof")
+        else:
+            self.incomplete("curvature_second_order_trend", reason="Fewer than two small steps above the heuristic arithmetic/solve floor")
+
+    def reduced(self, summary, lookup):
+        tolerance = summary["gravity_completion_tolerance"]
+        self.check("gravity_completion_tolerance_valid", math.isfinite(tolerance) and tolerance > 0.0)
+        completions = read_rows(self.folder / "numerical_strategy_completions.csv", ("kind", "alpha", "completed", "final_gravity_norm", "retained_state_change"))
+        for index, row in enumerate(completions):
+            self.check(f"completion_{index}_retained_unchanged", number(row, "retained_state_change") == 0.0)
+            if flag(row, "completed"):
+                residual = number(row, "final_gravity_norm")
+                self.check(f"completion_{index}_prepared_gravity_gate", nonnegative(residual) and residual <= tolerance,
+                           final_gravity_norm=residual, tolerance=tolerance)
+        self.outcomes["gravity_completion"] = {"attempts": len(completions), "completed": sum(flag(row, "completed") for row in completions)}
+        if summary.get("reduced_base_completed") != 1.0:
+            self.incomplete("reduced_base_available", reason="Base gravity completion did not reach its prepared-residual tolerance")
+            return
+        for kind in ("reduced_action", "reduced_lift"):
+            if (kind, 0.0, 0) not in lookup or not flag(lookup[kind, 0.0, 0], "converged"):
+                self.incomplete("reduced_action_available", reason=f"{kind} did not converge")
+                return
+        action = summary["reduced_action_norm"]
+        lifted_g = summary["lifted_gravity_action_norm"]
+        agreement = summary["lifted_retained_action_difference"]
+        lift = lookup["reduced_lift", 0.0, 0]
+        floor = roundoff(action, number(lift, "right_hand_side_norm"))
+        self.check("lifted_tangent_gravity_gate", nonnegative(lifted_g) and lifted_g <= number(lift, "threshold") + floor,
+                   lifted_gravity_action=lifted_g, allowed=number(lift, "threshold") + floor)
+        allowed_agreement = 1e-8 * abs(action) + floor
+        self.check("schur_matches_lifted_full_action", nonnegative(agreement) and agreement <= allowed_agreement,
+                   difference=agreement, allowed=allowed_agreement)
+        differences = read_rows(self.folder / "numerical_strategy_reduced_derivative.csv", ("alpha", "plus_completed", "minus_completed", "relative_error", "difference_quotient_error", "plus_gravity_norm", "minus_gravity_norm", "residual_subtraction_roundoff_scale"))
+        reliable = []
+        valid = []
+        for row in differences:
+            if not (flag(row, "plus_completed") and flag(row, "minus_completed")):
+                self.incomplete(f"reduced_fd_{row['alpha']}_completed_pair_available",
+                                plus_completed=flag(row, "plus_completed"), minus_completed=flag(row, "minus_completed"))
+                continue
+            alpha, error, relative, floor = (number(row, key) for key in ("alpha", "difference_quotient_error", "relative_error", "residual_subtraction_roundoff_scale"))
+            valid.append((alpha, error, relative, floor))
+            self.check(f"reduced_fd_{alpha:g}_finite", all(nonnegative(x) for x in (error, relative, floor)))
+            self.check(f"reduced_fd_{alpha:g}_completed_sides", all(nonnegative(number(row, key)) and number(row, key) <= tolerance for key in ("plus_gravity_norm", "minus_gravity_norm")))
+            if error > 20.0 * floor:
+                reliable.append((alpha, error))
+        self.outcomes["reduced_central_difference"] = {"completed_pairs": len(valid), "measurements": valid}
+        if not valid:
+            self.incomplete("reduced_central_difference_available", reason="No pair completed on both sides")
+            return
+        self.check("reduced_central_difference_accuracy", min(item[2] for item in valid) <= 1e-5,
+                   minimum_relative_error=min(item[2] for item in valid), tolerance=1e-5)
+        reliable.sort(reverse=True)
+        orders = [math.log(b[1] / a[1]) / math.log(b[0] / a[0]) for a, b in zip(reliable, reliable[1:]) if a[0] != b[0]]
+        if any(order >= 1.5 for order in orders):
+            self.check("reduced_central_difference_second_order_trend", True, observed_orders=orders)
+        elif all(item[1] <= 20.0 * item[3] for item in valid):
+            self.check("reduced_central_difference_floor_consistent", True,
+                       interpretation="Errors are at the recorded residual-subtraction scale; order is unresolved")
+        else:
+            self.incomplete("reduced_central_difference_trend", observed_orders=orders,
+                            reason="No resolved second-order interval; inner completion error has no measured retained-residual amplification bound")
+
+
+def assess(folder, required=(), expect_negative=False, restoration_tolerance=0.0):
+    assessment = Assessment(folder, restoration_tolerance)
+    sections = set(required)
+    if (folder / "polynomial_certificates.csv").exists() or expect_negative:
+        sections.add("polynomial")
+    if (folder / "numerical_strategy_summary.csv").exists():
+        sections.add("gravity")
+    if not sections:
+        assessment.incomplete("probe_artifacts_available", reason="No recognized frozen strategy outputs")
+    for section in sorted(sections):
+        try:
+            if section == "polynomial":
+                assessment.polynomial(expect_negative)
+            else:
+                assessment.gravity()
+        except FileNotFoundError as error:
+            assessment.incomplete(f"{section}_artifacts_complete", error=str(error))
+        except KeyError as error:
+            assessment.incomplete(f"{section}_records_complete", missing=str(error))
+        except (ValueError, OSError, csv.Error, ZeroDivisionError) as error:
+            assessment.check(f"{section}_artifacts_valid", False, error=str(error))
+    statuses = {item["status"] for item in assessment.checks}
+    status = "fail" if "fail" in statuses else "incomplete" if "incomplete" in statuses else "pass"
+    report = {"schema_version": 1, "result_folder": str(folder.resolve()), "status": status,
+              "checks": assessment.checks, "outcomes": assessment.outcomes,
+              "scope": "Frozen algebra, admissibility, and restoration checks; use check_plateau_physics.py for equilibrium/virial acceptance"}
+    return clean_json(report), {"pass": 0, "incomplete": 2, "fail": 3}[status]
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("result_folder", type=Path)
+    parser.add_argument("--require", choices=("polynomial", "gravity"), action="append", default=[])
+    parser.add_argument("--expect-negative-base", action="store_true")
+    parser.add_argument("--restoration-tolerance", type=float, default=0.0)
+    args = parser.parse_args()
+    if not nonnegative(args.restoration_tolerance):
+        parser.error("--restoration-tolerance must be finite and nonnegative")
+    report, code = assess(args.result_folder, args.require, args.expect_negative_base, args.restoration_tolerance)
+    print(json.dumps(report, indent=2, allow_nan=False))
+    return code
+
+
+if __name__ == "__main__":
+    sys.exit(main())

experiments/check_numerical_strategies_README.md

@@ -0,0 +1,72 @@
+# Checking numerical-strategy experiments
+
+Run the checker on each completed frozen-probe directory:
+
+```sh
+python3 experiments/check_numerical_strategies.py RESULT_FOLDER --require gravity
+python3 experiments/check_numerical_strategies.py RESULT_FOLDER --require polynomial
+python3 experiments/check_numerical_strategies.py LATE_ROTATING_FOLDER --require polynomial --expect-negative-base
+python3 experiments/test_check_numerical_strategies.py
+```
+
+Both sections are discovered automatically when their CSV files exist. Repeating
+`--require` demands both sections. Exit 0 means the requested checks passed;
+3 means an invariant or numerical expectation failed; 2 means evidence is
+incomplete. An honestly reported failed inner solve, unavailable trial, or
+unresolved certificate must not turn into a passing verification. Expected
+negative-base detection and rejection are successful guard tests, not successful
+equilibrium solves.
+
+Accepted state, residual, physical state, saved correction, and retained-coordinate
+changes must be exactly zero. Restored residual drift defaults to exactly zero;
+`--restoration-tolerance` relaxes only that comparison when explicitly justified.
+The JSON output records the chosen tolerance.
+
+For successful curvature probes, the checker compares the fresh gravity remainder
+with the selected inner threshold plus a documented arithmetic scale, verifies
+that only gravity changed, and checks correction/alpha² arithmetic. It requires
+two small steps resolved above that scale to assess a second-order trend. An
+observed order of at least 1.7 passes the finite-range trend check; this does not
+prove an asymptotic rate. Higher order is allowed because a quadratic coefficient
+can vanish.
+
+Reduced probes must satisfy the freshly prepared gravity completion tolerance,
+the lifted gravity tangent gate, and agreement of the Schur action with the
+retained part of the full lifted action. Central differences must reach relative
+error at most 1e-5 and show a resolved interval with order at least 1.5, or reach
+the recorded residual-subtraction scale. A plateau above that scale is incomplete
+evidence rather than an automatic claim of a bad Jacobian: the retained-residual
+amplification of the remaining gravity error has not been bounded. The recorded
+roundoff scales are heuristics, not rigorous forward-error bounds.
+
+This script deliberately does not duplicate physical acceptance. For each
+end-to-end model, run `check_plateau_physics.py` on its endpoint. That separate
+gate requires recorded strict nonlinear convergence and all eight defined virial
+measurements at or below 1e-5. Its local EOS, hydrostatic, and sign measurements
+remain part of the physical assessment.
+
+## Bounded validation matrix
+
+These are proposed operational limits, not measured performance predictions.
+Use one numerical process at a time and preserve all capped failures as results.
+
+| Layer | Cases | Bound / acceptance |
+|---|---|---|
+| Algebra and lifecycle | Gravity subsystem, Schur reduction, polynomial certificate tests; one and two MPI ranks | 2 minutes per test group; true-defect, stale-state, empty-local, rollback, type-contract coverage |
+| Python checker | Synthetic passing, false-convergence, mutation, missing-data, negative-base, and NaN records | 30 seconds; standard-library unit tests |
+| Positivity frozen states | n=1 J=0.1 captured iteration 2 and iteration 11; a certified seed/control | 10 minutes per probe; iteration 11 must be negative and reject repair with `invalid_base`; no claim that the earlier base is positive before checking it |
+| Gravity curvature | n=1.5 J=0 captured iteration 1 | 15 minutes; alpha = 1, 0.1, 0.01, 0.001; gate fresh remainder cancellation and small-step order |
+| Reduced gravity | Same n=1.5 state, plus constrained n=3 saved seed direction | 15 minutes per probe; completion and three central-difference pairs; n=3 direction remains an approximate diagnostic correction |
+| Regression equilibrium | n=1 J=0 and low-rotation J=0.05 baseline; opt-in strategies individually | Preserve existing nonlinear/linear caps, add a 20-minute wall limit per model; strict solver and all virial gates remain mandatory |
+| Failure-prone equilibrium | n=1 J=0.1 and n=1.5 J=0, strategies separately before combinations | 8 Newton steps / 20 minutes; report residual progress, auxiliary work, rejection causes, and physical endpoint without equating the cap with convergence |
+
+At the supplied frozen-gravity defaults, there are at most 20 explicitly requested
+inner solves, each capped at 160 iterations, and at most 36 preparations in the
+normal completed path. The maximum covers four curvature corrections, seven
+gravity completions with up to two corrections each, and two reduced actions.
+Rejected preparations can reduce those counts. Root integration controls and
+end-to-end outer solves are separately bounded; their nested work is not covered
+by this frozen-probe count.
+
+Treat outer reduced-solve time as inclusive of nested gravity solves. Report the
+inner counters separately, but do not add the two times as if they were disjoint.

experiments/test_check_numerical_strategies.py

@@ -0,0 +1,170 @@
+#!/usr/bin/env python3
+"""Exercise checker acceptance and failure paths with synthetic probe records."""
+
+import csv
+import json
+from pathlib import Path
+import tempfile
+import unittest
+
+from check_numerical_strategies import assess
+
+
+def write_rows(folder, name, rows, fields=None):
+    fields = fields or list(rows[0])
+    with (folder / name).open("w", newline="") as stream:
+        writer = csv.DictWriter(stream, fieldnames=fields)
+        writer.writeheader()
+        writer.writerows(rows)
+
+
+def write_metrics(folder, name, values):
+    write_rows(folder, name, [{"metric": key, "value": value} for key, value in values.items()])
+
+
+def gravity_fixture(folder):
+    summary = {"accepted_state_change": 0, "accepted_residual_change": 0, "accepted_physical_change": 0,
+               "maximum_original_restoration_error": 0, "gravity_completion_tolerance": 1e-11,
+               "reduced_base_completed": 1, "reduced_action_norm": 2,
+               "lifted_gravity_action_norm": 1e-13, "lifted_retained_action_difference": 1e-13}
+    inner = []
+    curvature = []
+    for alpha in (1., .1, .01, .001):
+        rhs = .1 * alpha**2
+        threshold = max(1e-13, 1e-9 * rhs)
+        residual = threshold * .1
+        inner.append({"kind": "curvature", "alpha": alpha, "index": 0, "status": 0, "converged": 1,
+                      "right_hand_side_norm": rhs, "true_residual_norm": residual,
+                      "relative_true_residual_norm": residual / rhs, "threshold": threshold})
+        curvature.append({"alpha": alpha, "prepared": 1, "inner_converged": 1, "corrected_prepared": 1,
+                          "base_residual_norm": 1, "raw_residual_norm": 1, "corrected_residual_norm": .99,
+                          "linear_prediction_norm": .99, "raw_gravity_remainder": rhs,
+                          "corrected_gravity_remainder": residual, "gravity_correction_norm": .2 * alpha**2,
+                          "correction_over_alpha_squared": .2, "retained_state_change": 0})
+    for kind in ("reduced_action", "reduced_lift"):
+        inner.append({"kind": kind, "alpha": 0, "index": 0, "status": 0, "converged": 1,
+                      "right_hand_side_norm": 1, "true_residual_norm": 1e-13,
+                      "relative_true_residual_norm": 1e-13, "threshold": 1e-9})
+    completions = [{"kind": "reduced_base", "alpha": 0, "completed": 1,
+                    "final_gravity_norm": 5e-14, "retained_state_change": 0}]
+    differences = []
+    for alpha in (.01, .001, .0001):
+        for kind in ("reduced_plus", "reduced_minus"):
+            completions.append({"kind": kind, "alpha": alpha, "completed": 1,
+                                "final_gravity_norm": 5e-14, "retained_state_change": 0})
+        differences.append({"alpha": alpha, "plus_completed": 1, "minus_completed": 1,
+                            "relative_error": .25 * alpha**2, "difference_quotient_error": .5 * alpha**2,
+                            "plus_gravity_norm": 5e-14, "minus_gravity_norm": 5e-14,
+                            "residual_subtraction_roundoff_scale": 1e-14 / alpha})
+    write_metrics(folder, "numerical_strategy_summary.csv", summary)
+    write_rows(folder, "numerical_strategy_inner.csv", inner)
+    write_rows(folder, "numerical_strategy_curvature.csv", curvature)
+    write_rows(folder, "numerical_strategy_completions.csv", completions)
+    write_rows(folder, "numerical_strategy_reduced_derivative.csv", differences)
+    return summary, inner, curvature, completions, differences
+
+
+def polynomial_fixture(folder):
+    certificates = [{"state": name, "status": status, "witness_element": 459 if status == "negative" else -1,
+                     "x": .5, "y": .5, "z": .5}
+                    for name, status in (("frozen_base", "negative"), ("context_seed_anchor", "certified_nonnegative"),
+                                         ("original_full_candidate", "negative"))]
+    repair = [{"status": "invalid_base", "minimum_safe_fraction": 1, "base_status": "negative"}]
+    trials = [{"fraction": 0, "alpha": 1, "certificate": "negative", "prepared": 1, "merit_slope": -.1}]
+    write_rows(folder, "polynomial_certificates.csv", certificates)
+    write_rows(folder, "polynomial_repair.csv", repair)
+    write_rows(folder, "polynomial_repair_trials.csv", trials)
+    write_metrics(folder, "polynomial_restoration.csv", {
+        "accepted_physical_difference": 0, "accepted_normalized_difference": 0, "accepted_residual_difference": 0,
+        "restored_correction_difference": 0, "restored_residual_difference": 0})
+    return certificates, repair, trials
+
+
+class NumericalStrategyChecks(unittest.TestCase):
+    def setUp(self):
+        self.temporary = tempfile.TemporaryDirectory()
+        self.folder = Path(self.temporary.name)
+
+    def tearDown(self):
+        self.temporary.cleanup()
+
+    def test_complete_gravity_evidence_passes(self):
+        gravity_fixture(self.folder)
+        report, code = assess(self.folder, required=("gravity",))
+        self.assertEqual(code, 0, report)
+
+    def test_negative_late_base_is_an_expected_guard_success(self):
+        polynomial_fixture(self.folder)
+        report, code = assess(self.folder, expect_negative=True)
+        self.assertEqual(code, 0, report)
+
+    def test_negative_base_cannot_be_repaired_silently(self):
+        _, repair, _ = polynomial_fixture(self.folder)
+        repair[0]["status"] = "feasible"
+        write_rows(self.folder, "polynomial_repair.csv", repair)
+        report, code = assess(self.folder, expect_negative=True)
+        self.assertEqual(code, 3, report)
+
+    def test_false_convergence_report_fails(self):
+        _, inner, _, _, _ = gravity_fixture(self.folder)
+        inner[0]["true_residual_norm"] = 1.0
+        write_rows(self.folder, "numerical_strategy_inner.csv", inner)
+        report, code = assess(self.folder)
+        self.assertEqual(code, 3, report)
+
+    def test_changed_retained_state_fails(self):
+        _, _, curvature, _, _ = gravity_fixture(self.folder)
+        curvature[0]["retained_state_change"] = 1e-8
+        write_rows(self.folder, "numerical_strategy_curvature.csv", curvature)
+        report, code = assess(self.folder)
+        self.assertEqual(code, 3, report)
+
+    def test_failed_numerical_probe_is_incomplete_not_a_false_pass(self):
+        summary, inner, curvature, completions, _ = gravity_fixture(self.folder)
+        summary["reduced_base_completed"] = 0
+        for row in inner:
+            row["status"], row["converged"] = 1, 0
+        for row in curvature:
+            row["inner_converged"], row["corrected_prepared"] = 0, 0
+        for row in completions:
+            row["completed"] = 0
+        write_metrics(self.folder, "numerical_strategy_summary.csv", summary)
+        write_rows(self.folder, "numerical_strategy_inner.csv", inner)
+        write_rows(self.folder, "numerical_strategy_curvature.csv", curvature)
+        write_rows(self.folder, "numerical_strategy_completions.csv", completions)
+        report, code = assess(self.folder)
+        self.assertEqual(code, 2, report)
+        self.assertEqual(report["outcomes"]["gravity_inner"]["failed_calls"], 6)
+
+    def test_missing_probe_is_incomplete(self):
+        report, code = assess(self.folder, required=("gravity",))
+        self.assertEqual(code, 2, report)
+
+    def test_unfinished_restoration_record_is_incomplete(self):
+        summary, _, _, _, _ = gravity_fixture(self.folder)
+        del summary["accepted_state_change"]
+        write_metrics(self.folder, "numerical_strategy_summary.csv", summary)
+        report, code = assess(self.folder)
+        self.assertEqual(code, 2, report)
+
+    def test_nonfinite_restoration_fails_and_serializes_standard_json(self):
+        summary, _, _, _, _ = gravity_fixture(self.folder)
+        summary["accepted_state_change"] = float("nan")
+        write_metrics(self.folder, "numerical_strategy_summary.csv", summary)
+        report, code = assess(self.folder)
+        self.assertEqual(code, 3, report)
+        json.dumps(report, allow_nan=False)
+
+    def test_only_restoration_not_accepted_state_may_use_optional_tolerance(self):
+        summary, _, _, _, _ = gravity_fixture(self.folder)
+        summary["maximum_original_restoration_error"] = 1e-13
+        write_metrics(self.folder, "numerical_strategy_summary.csv", summary)
+        self.assertEqual(assess(self.folder)[1], 3)
+        self.assertEqual(assess(self.folder, restoration_tolerance=1e-12)[1], 0)
+        summary["accepted_state_change"] = 1e-13
+        write_metrics(self.folder, "numerical_strategy_summary.csv", summary)
+        self.assertEqual(assess(self.folder, restoration_tolerance=1e-12)[1], 3)
+
+
+if __name__ == "__main__":
+    unittest.main()