376 lines
23 KiB
Python
376 lines
23 KiB
Python
#!/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
|
|
APPLICATION_AGREEMENT_TOLERANCE = 2.0e-12
|
|
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 applications(self):
|
|
fields = ("direction", "compact_operator", "compact_inverse",
|
|
"operator_norm", "operator_error", "operator_relative_error",
|
|
"compact_operator_seconds", "root_operator_seconds",
|
|
"inverse_norm", "inverse_error", "inverse_relative_error",
|
|
"compact_inverse_seconds", "root_inverse_seconds")
|
|
rows = read_rows(self.folder / "gravity_applications.csv", fields)
|
|
outcome = {"recorded_directions": len(rows), "measurements": [],
|
|
"relative_tolerance": APPLICATION_AGREEMENT_TOLERANCE,
|
|
"absolute_floor": APPLICATION_AGREEMENT_TOLERANCE,
|
|
"inverse_equivalence_assumption": "Default independent material/gravity recipe; not a universal requirement for custom compact inverses"}
|
|
self.outcomes["gravity_applications"] = outcome
|
|
if not rows:
|
|
self.incomplete("gravity_applications_available", reason="No application directions recorded")
|
|
return
|
|
directions = set()
|
|
paths = set()
|
|
for row in rows:
|
|
direction = int(row["direction"])
|
|
if direction < 0 or direction in directions:
|
|
raise ValueError("Gravity application direction indices must be unique and nonnegative")
|
|
directions.add(direction)
|
|
compact_operator = flag(row, "compact_operator")
|
|
compact_inverse = flag(row, "compact_inverse")
|
|
paths.add((compact_operator, compact_inverse))
|
|
measurement = {"direction": direction,
|
|
"operator_path": "compact" if compact_operator else "root_restriction",
|
|
"inverse_path": "compact" if compact_inverse else "root_restriction"}
|
|
for kind in ("operator", "inverse"):
|
|
prefix = f"application_{direction}_{kind}"
|
|
scale, error, relative = (number(row, f"{kind}_{suffix}") for suffix in ("norm", "error", "relative_error"))
|
|
actual_seconds, reference_seconds = (number(row, f"{path}_{kind}_seconds") for path in ("compact", "root"))
|
|
self.check(f"{prefix}_finite_nonnegative", all(nonnegative(value) for value in
|
|
(scale, error, relative, actual_seconds, reference_seconds)))
|
|
# Match the real-FEM fixture's error/max(1, reference norm)
|
|
# comparison. Near a null action this is an absolute floor;
|
|
# a large reported relative error must not replace this gate.
|
|
tolerance = APPLICATION_AGREEMENT_TOLERANCE * max(1.0, scale)
|
|
self.check(f"{prefix}_same_action", nonnegative(scale) and nonnegative(error) and error <= tolerance,
|
|
error=error, reference_norm=scale, allowed=tolerance)
|
|
expected_relative = error / scale if scale > 0.0 else (0.0 if error == 0.0 else math.inf)
|
|
self.check(f"{prefix}_relative_error_arithmetic", nonnegative(relative) and
|
|
math.isclose(relative, expected_relative, rel_tol=1e-12, abs_tol=1e-30))
|
|
measurement[kind] = {"reference_norm": scale, "error": error, "relative_error": relative,
|
|
"allowed": tolerance, "selected_seconds": actual_seconds,
|
|
"root_seconds": reference_seconds}
|
|
outcome["measurements"].append(measurement)
|
|
self.check("gravity_application_paths_consistent", len(paths) == 1)
|
|
outcome["compact_operator_directions"] = sum(flag(row, "compact_operator") for row in rows)
|
|
outcome["compact_inverse_directions"] = sum(flag(row, "compact_inverse") for row in rows)
|
|
|
|
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, floor)) and
|
|
(nonnegative(relative) or abs(action) <= roundoff(action)))
|
|
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", any(
|
|
(nonnegative(item[2]) and item[2] <= 1e-5) or item[1] <= 20.0 * item[3] for item in valid),
|
|
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 (folder / "gravity_applications.csv").exists():
|
|
sections.add("applications")
|
|
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)
|
|
elif section == "gravity":
|
|
assessment.gravity()
|
|
elif section == "applications":
|
|
assessment.applications()
|
|
else:
|
|
raise ValueError(f"Unknown required probe: {section}")
|
|
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", "applications"), 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())
|