#!/usr/bin/env python3 """Audit root-split density artifacts against the original frozen campaign. Saved polynomial algebra and aggregation are checked independently. Mapped quadrature weights are not exported, so this is not an independent FEM integral or a global positivity certificate. Tangential convergence remains an observation. Run after the serial numerical campaign, never alongside a timed FEM run. """ from __future__ import annotations import argparse from collections import defaultdict import csv from decimal import Decimal, localcontext import hashlib import json import math from pathlib import Path import tempfile import unittest EPS = 2.220446049250313e-16 PREFIX = "density_negative_regions_" ORDERS = (14, 22, 30) SUM_FIELDS = ("tangent_lines", "negative_intervals", "integration_points", "volume", "mass", "negative_volume", "negative_mass") NEGATIVE_BOUNDS = ("negative_radius_minimum", "negative_radius_maximum", "negative_interval_t_minimum", "negative_interval_t_maximum") CONSUMED = ("metadata.txt", "layout.csv", "observed_state.txt", "observed_density_coefficients.txt", "density_invariants.csv", "density_physical.csv") def require(condition, message): if not condition: raise AssertionError(message) def close(a, b, label, rtol=3e-11, atol=1e-25): require(math.isfinite(a) and math.isfinite(b) and math.isclose(a, b, rel_tol=rtol, abs_tol=atol), f"{label}: {a:.17g} != {b:.17g}") def digest(path): with path.open("rb") as stream: return hashlib.file_digest(stream, "sha256").hexdigest() def rows(path): with path.open(newline="") as stream: reader = csv.DictReader(stream) require(reader.fieldnames is not None and len(set(reader.fieldnames)) == len(reader.fieldnames), f"Missing or duplicate CSV columns: {path}") result = list(reader) require(all(None not in row and None not in row.values() for row in result), f"Malformed CSV: {path}") return result def metadata(path): entries = [line.split("=", 1) for line in path.read_text().splitlines() if "=" in line] require(len({key for key, _ in entries}) == len(entries), f"Duplicate metadata: {path}") return dict(entries) def keyed(table, key): result = {key(row): row for row in table} require(len(result) == len(table), "Duplicate artifact rows") return result def numeric(row): return {key: float(value) for key, value in row.items()} def finite(row, allowed_nan=()): for key, value in row.items(): require(math.isfinite(value) or (key in allowed_nan and math.isnan(value)), f"Unexpected nonfinite value in {key}") class Campaign: def __init__(self, path): self.path = path.resolve() self.ledger = {} self.provenance = json.loads(self.record(self.path / "provenance.json").read_text()) p = self.provenance require(digest(self.record(self.path / "plateau_diagnostic")) == p["binary_sha256"], "Producer binary changed") require(digest(self.record(Path(p["mesh"]))) == p["mesh_sha256"], "Producer mesh changed") require(json.loads(self.record(self.path / "cases.json").read_text()) == p["cases"], "Case manifest changed") self.cases = keyed(p["cases"], lambda case: case["id"]) self.runs = {} def record(self, path): path = path.resolve() self.ledger[str(path)] = digest(path) return path def run(self, name): if name not in self.runs: r = json.loads(self.record(self.path / f"{name}_run.json").read_text()) require(r["id"] == name and r["exit_code"] == 0 and r["inputs_unchanged"] is True, f"Incomplete or changed case: {name}") expected = [str(self.path / "plateau_diagnostic"), *self.cases[name]["arguments"], "--mesh", self.provenance["mesh"], "--output", str(self.path / name)] require(r["command"] == expected, f"Recorded command differs from manifest: {name}") for path, expected_hash in r["input_state_hashes"].items(): require(self.ledger.get(str(self.record(Path(path)))) == expected_hash, f"Changed source input: {path}") self.runs[name] = r return self.runs[name] def artifact(self, name, relative): r = self.run(name) require(relative in r["output_hashes"], f"Consumed artifact lacks producer hash: {name}/{relative}") path = self.record(self.path / name / relative) require(self.ledger[str(path)] == r["output_hashes"][relative], f"Changed output: {name}/{relative}") return path def vector_size(path): with path.open() as stream: size = int(next(stream).strip()) require(size > 0, "Empty saved vector") count = 0 for line in stream: require(math.isfinite(float(line)), "Nonfinite saved vector") count += 1 require(size == count, "Saved vector size differs from its header") return size def gauss_rule(count): """Independent Legendre recurrence on [0,1], used only for node binding.""" result = [] for j in range(1, count + 1): z = math.cos(math.pi * (j - .25) / (count + .5)) for _ in range(30): previous, current = 1., z for degree in range(2, count + 1): previous, current = current, ((2 * degree - 1) * z * current - (degree - 1) * previous) / degree derivative = count * (z * current - previous) / (z * z - 1) update = current / derivative z -= update if abs(update) < 2e-16: break result.append(((1 + z) / 2, 1 / ((1 - z * z) * derivative * derivative))) return sorted(result) def polynomial_intervals(c, b, a): """Direct high-precision quadratic formula differs from producer q/a,c/q.""" with localcontext() as context: context.prec = 70 c, b, a = [Decimal(str(x)) for x in (c, b, a)] points = [Decimal(0), Decimal(1)] if a == 0: roots = [] if b == 0 else [-c / b] else: discriminant = b * b - 4 * a * c roots = [] if discriminant < 0 else [(-b - discriminant.sqrt()) / (2 * a), (-b + discriminant.sqrt()) / (2 * a)] points = sorted(set(points + [t for t in roots if 0 < t < 1])) value = lambda t: c + t * (b + t * a) return [(float(left), float(right)) for left, right in zip(points, points[1:]) if value((left + right) / 2) < 0] def interval_check(row, rule): r = numeric(row) finite(r) left, right = r["t_left"], r["t_right"] require(r["tangential_order"] == 30 and 0 <= left < right <= 1, "Invalid finest negative interval") require(r["axis"] in (0, 1, 2) and r["orientation"] in (-1, 1), "Invalid interval axis") a, b, c = (r[key] for key in ("coefficient_a", "coefficient_b", "coefficient_c")) scale = max(abs(a), abs(b), abs(c)) require(scale > 0, "Zero polynomial cannot have a negative interval") allowance = 4096 * EPS * scale evaluate = lambda t: math.fsum((c, b * t, a * t * t)) for endpoint in (left, right): if 0 < endpoint < 1: require(abs(evaluate(endpoint)) <= allowance, "Interior interval endpoint is not a root") for fraction in (.125, .25, .5, .75, .875): require(evaluate(left + fraction * (right - left)) <= allowance, "Negative interval contains a positive value") candidates = [left, right] if a and left < -b / (2 * a) < right: candidates.append(-b / (2 * a)) minimum = min(map(evaluate, candidates)) require(minimum <= allowance, "Claimed negative interval has positive minimum") close(r["minimum_density"], minimum, "Direct FE and polynomial interval minimum", rtol=0, atol=allowance) require(r["negative_volume"] > 0 and r["negative_mass"] > 0, "Nonpositive negative-region integral") require(r["negative_mass"] / r["negative_volume"] <= -minimum + allowance, "Negative mean density exceeds the polynomial's maximum magnitude") weights = [] for coordinate in ("u", "v"): nearest = min(rule, key=lambda pair: abs(pair[0] - r[coordinate])) close(nearest[0], r[coordinate], "Tangential Gauss node", rtol=0, atol=4e-15) weights.append(nearest[1]) close(weights[0] * weights[1], r["tangential_weight"], "Tangential Gauss weight", rtol=2e-12) radii = [math.sqrt(math.fsum(r[f"{side}_{component}"] ** 2 for component in "xyz")) for side in ("left", "right")] require(0 <= r["radius_minimum"] <= r["radius_maximum"], "Invalid radius bounds") tolerance = 2e-13 * max(1., r["radius_maximum"]) require(r["radius_minimum"] <= min(radii) + tolerance and r["radius_maximum"] >= max(radii) - tolerance, "Sampled radius bounds exclude a saved endpoint") return r, minimum < -allowance def measure_check(r): empty = r["negative_intervals"] == 0 finite(r, NEGATIVE_BOUNDS if empty else ()) require(r["volume"] > 0 and 0 <= r["negative_volume"] <= r["volume"], "Invalid measured volume") require(r["negative_mass"] >= 0, "Invalid negative mass") for key in SUM_FIELDS[:3]: require(r[key] >= 0 and r[key].is_integer(), f"Noninteger count: {key}") close(r["negative_volume_fraction"], r["negative_volume"] / r["volume"], "Negative volume fraction") require(r["mass"] != 0, "Zero total mass requires an explicitly unsupported fraction convention") close(r["negative_mass_fraction"], r["negative_mass"] / abs(r["mass"]), "Negative mass fraction") require(0 <= r["maximum_reconstruction_error"] <= r["maximum_reconstruction_tolerance"], "Q2 reconstruction failed") require(all(0 <= r[key] <= 1 for key in ("witness_t", "witness_u", "witness_v")), "Witness outside reference cell") if empty: require(r["negative_mass"] == 0 and r["negative_volume"] == 0 and all(math.isnan(r[key]) for key in NEGATIVE_BOUNDS), "Empty negative region has nonempty observations") else: require(r["negative_mass"] > 0 and r["negative_volume"] > 0 and 0 <= r["negative_interval_t_minimum"] < r["negative_interval_t_maximum"] <= 1 and 0 <= r["negative_radius_minimum"] <= r["negative_radius_maximum"], "Invalid nonempty negative bounds") def audit_selfchecks(table): expected = {"positive_quadratic": (0., 0.), "negative_constant": (1., 2.), "linear_crossing": (.25, .03125), "double_root_nonnegative": (0., 0.), "near_linear_positive": (0., 0.), "near_double_negative": (2. ** -19, 4 * 2. ** -60 / 3), "gauss_q2_t_power_1p5": (.04548403175047501373888, .0005262295634284655), "gauss_q2_t_power_3": (math.sqrt(6) / 15, math.sqrt(6) / 2250)} indexed = keyed(table, lambda r: r["case"]) require(indexed.keys() == expected.keys(), "Analytic selfcheck coverage changed") for name, (length, mass) in expected.items(): r = indexed[name] require(r["passed"] == "1", "Producer arithmetic selfcheck failed") for key, target in (("negative_length", length), ("negative_mass", mass)): close(float(r[key]), target, f"Independent selfcheck {name}/{key}", rtol=3e-12, atol=1e-27) close(float(r["expected_" + key]), target, "Producer selfcheck expected value", rtol=3e-12, atol=1e-27) def audit_case(current, original, name): new = lambda filename: current.artifact(name, filename) old = lambda filename: original.artifact(name, filename) meta, prior_meta = metadata(new("metadata.txt")), metadata(old("metadata.txt")) for key in ("n", "K", "J", "rho_c", "M", "R_reference", "density_source", "surface_prescription", "surface_degree", "surface_maximum_order", "central_phase"): require(key in meta and key in prior_meta and meta[key] == prior_meta[key], f"Changed model metadata: {key}") require(meta["density_diagnostic"] == "1" and meta["density_region_integral"] == "1", "Wrong diagnostic mode") require(current.run(name)["input_state_hashes"] == original.run(name)["input_state_hashes"], "Changed frozen inputs") for filename in ("observed_state.txt", "observed_density_coefficients.txt", "layout.csv"): require(digest(new(filename)) == digest(old(filename)), f"Source is not byte-identical: {filename}") state_size = vector_size(new("observed_state.txt")) density_size = vector_size(new("observed_density_coefficients.txt")) layout = rows(new("layout.csv")) for kind in ("value", "residual"): blocks = sorted((int(r["offset"]), int(r["size"])) for r in layout if r["kind"] == kind) end = 0 for offset, size in blocks: require(offset == end and size > 0, "Invalid saved block layout") end += size require(end == state_size, "Layout does not cover the observed state") require(density_size % 27 == 0, "Full density storage is not Q2 tensor sized") invariants = {r["metric"]: float(r["value"]) for r in rows(new("density_invariants.csv"))} previous_invariants = {r["metric"]: float(r["value"]) for r in rows(old("density_invariants.csv"))} finite(invariants) for key in ("accepted_physical_unchanged", "accepted_normalized_unchanged", "accepted_residual_unchanged"): require(invariants[key] == previous_invariants[key] == 1, "Accepted-state observation/restoration failed") close(invariants["prepared_residual_norm"], previous_invariants["prepared_residual_norm"], "Frozen residual", atol=1e-14) scope_file = new(PREFIX + "scope.txt") scope = metadata(scope_file) require(scope["source_fields_bitwise_identical"] == "1" and scope["trial_preparations"] == "0" and scope["mpi_ranks"] == "1", "Region observer invariants failed") for key in ("n", "K"): close(float(scope[key]), float(meta[key]), "Observer model binding") text = scope_file.read_text() for phrase in ("no density or EOS clipping", "All stellar cells are included", "tangential integration is sampled", "not global extrema or polynomial sign certificates", "Near-double-root"): require(phrase in text, "Missing observer scope limitation") audit_selfchecks(rows(new(PREFIX + "selfchecks.csv"))) totals = keyed([numeric(r) for r in rows(new(PREFIX + "summary.csv"))], lambda r: int(r["tangential_order"])) elements = keyed([numeric(r) for r in rows(new(PREFIX + "elements.csv"))], lambda r: (int(r["tangential_order"]), int(r["element"]))) topology = keyed(rows(old("density_element_topology.csv")), lambda r: int(r["element"])) require(totals.keys() == set(ORDERS), "Missing quadrature levels") require(elements.keys() == {(order, element) for order in ORDERS for element in topology}, "Stellar support differs") original_summary = keyed(rows(old("density_summary.csv")), lambda r: (int(r["quadrature"]), r["field"])) physical = {r["metric"]: float(r["value"]) for r in rows(new("density_physical.csv")) if r["quadrature"] == "18"} denominator_differences = {} for order, total in totals.items(): require(total["interval_order"] == 14, "Changed interval quadrature") measure_check(total) subset = [elements[(order, element)] for element in topology] for r in subset: measure_check(r) top = topology[int(r["element"])] require(r["interval_order"] == 14 and r["tangent_lines"] == (order // 2 + 1) ** 2 and r["attribute"] == int(top["attribute"]) and r["surface_face_count"] == int(top["surface_face_count"]), "Element integration layout or topology differs") require(r["witness_element"] == r["element"] and r["axis"] in (0, 1, 2) and r["orientation"] in (-1, 1) and r["integration_points"] % 8 == 0 and 8 * r["tangent_lines"] <= r["integration_points"] <= 24 * r["tangent_lines"], "Element count/axis contradiction") require((r["surface_face"] >= 0) == (r["surface_face_count"] > 0), "Surface axis metadata inconsistent") if r["surface_face_count"] == 0: require(r["axis"] == 0 and r["orientation"] == 1, "Non-surface slicing convention changed") for key in SUM_FIELDS: close(math.fsum(r[key] for r in subset), total[key], f"Global element sum {order}/{key}") close(min(r["minimum_density"] for r in subset), total["minimum_density"], "Global minimum") witness = elements[(order, int(total["witness_element"]))] for key in ("minimum_density", "witness_h", "witness_t", "witness_u", "witness_v", "witness_x", "witness_y", "witness_z"): close(witness[key], total[key], "Global witness binding", atol=1e-15) nonempty = [r for r in subset if r["negative_intervals"]] for key, function in (("negative_radius_minimum", min), ("negative_radius_maximum", max), ("negative_interval_t_minimum", min), ("negative_interval_t_maximum", max)): if nonempty: close(function(r[key] for r in nonempty), total[key], "Negative-region bounds") for key in ("maximum_reconstruction_error", "maximum_reconstruction_tolerance"): close(max(r[key] for r in subset), total[key], "Reconstruction maximum") denominator_differences[str(order)] = {} for key in ("mass", "volume"): for old_order in (18, 26): reference = float(original_summary[(old_order, "state")][key]) close(total[key], reference, f"Independent full {key}, q{order}/original q{old_order}", rtol=2e-9, atol=1e-13) denominator_differences[str(order)][f"{key}_relative_to_q{old_order}"] = (total[key] - reference) / abs(reference) close(physical[key], float(original_summary[(18, "state")][key]), "Repeated q18 physical measure") interval_rows = rows(new(PREFIX + "intervals.csv")) by_element, by_line = defaultdict(list), defaultdict(list) unresolved = 0 rule = gauss_rule(16) for row in interval_rows: r, resolved = interval_check(row, rule) unresolved += not resolved element = int(r["element"]) require((30, element) in elements, "Interval on an unsupported element") for key in ("attribute", "axis", "orientation", "surface_face"): require(r[key] == elements[(30, element)][key], "Interval geometry metadata mismatch") by_element[element].append(r) by_line[(element, r["u"], r["v"])].append((row, r, resolved)) for line in by_line.values(): first = line[0][0] for row, _, _ in line: require(all(row[key] == first[key] for key in ("coefficient_c", "coefficient_b", "coefficient_a")), "One line has multiple polynomials") expected = polynomial_intervals(first["coefficient_c"], first["coefficient_b"], first["coefficient_a"]) actual = sorted((r["t_left"], r["t_right"]) for _, r, _ in line) require(all(actual[i][1] <= actual[i + 1][0] for i in range(len(actual) - 1)), "Overlapping negative intervals") if all(resolved for _, _, resolved in line): require(len(expected) == len(actual), "Missing or extra resolved negative interval") a, b, c = [float(first[k]) for k in ("coefficient_a", "coefficient_b", "coefficient_c")] uncertainty = 4096 * EPS * max(abs(a), abs(b), abs(c)) for got, wanted in zip(actual, expected): for x, y in zip(got, wanted): if y in (0., 1.): require(x == y, "Negative interval misses the physical cell endpoint") continue # Near a double root, position uncertainty scales as the # square root of coefficient uncertainty, not machine eps. allowance = max(64 * EPS, uncertainty / max(abs(b + 2 * a * y), math.sqrt(abs(a) * uncertainty))) close(x, y, "Independent Decimal root", rtol=0, atol=allowance) for element in topology: r = elements[(30, element)] pieces = by_element[element] require(len(pieces) == r["negative_intervals"], "Finest interval export coverage mismatch") for key in ("negative_mass", "negative_volume"): close(math.fsum(piece[key] for piece in pieces), r[key], "Finest interval-to-element sum") if pieces: allowance = r["maximum_reconstruction_tolerance"] + 128 * EPS * max(abs(r["minimum_density"]), 1e-30) close(min(piece["minimum_density"] for piece in pieces), r["minimum_density"], "Finest negative minimum", rtol=0, atol=allowance) for interval_key, element_key, function in (("t_left", "negative_interval_t_minimum", min), ("t_right", "negative_interval_t_maximum", max), ("radius_minimum", "negative_radius_minimum", min), ("radius_maximum", "negative_radius_maximum", max)): close(function(p[interval_key] for p in pieces), r[element_key], "Finest interval bounds") comparison = {} for key in ("negative_mass", "negative_volume", "negative_mass_fraction", "negative_volume_fraction"): comparison[key] = {"q14": totals[14][key], "q22": totals[22][key], "q30": totals[30][key], "q30_minus_q22": totals[30][key] - totals[22][key], "relative_q30_minus_q22": ((totals[30][key] - totals[22][key]) / totals[30][key] if totals[30][key] else None)} return {"artifact_passed": True, "n": float(meta["n"]), "state_size": state_size, "density_storage_size": density_size, "stellar_elements": len(topology), "finest_exported_negative_intervals": len(interval_rows), "finest_negative_lines": len(by_line), "intervals_below_coefficient_uncertainty": unresolved, "finest_affected_elements": len([v for v in by_element.values() if v]), "finest_negative_elements_not_touching_surface": sum(bool(by_element[e]) and int(topology[e]["surface_face_count"]) == 0 for e in topology), "denominator_comparisons": denominator_differences, "quadrature_comparison": comparison, "summary": {str(k): v for k, v in totals.items()}} class SelfTests(unittest.TestCase): def test_decimal_roots(self): self.assertEqual(polynomial_intervals("-.25", "1", "0"), [(0., .25)]) self.assertEqual(polynomial_intervals(".25", "-1", "1"), []) roots = polynomial_intervals(".05", "-.6", "1.5") close(roots[0][0], .2 - math.sqrt(6) / 30, "Cubic root") close(roots[0][1], .2 + math.sqrt(6) / 30, "Cubic root") def fixture(self): rule = gauss_rule(16) u, weight = rule[3] r = dict(tangential_order=30, element=1, attribute=1, axis=0, orientation=1, surface_face=2, u=u, v=u, tangential_weight=weight * weight, t_left=0, t_right=.25, coefficient_c=-.25, coefficient_b=1, coefficient_a=0, negative_volume=.25 * weight * weight, negative_mass=.03125 * weight * weight, minimum_density=-.25, witness_h=1, radius_minimum=math.sqrt(2 * u * u), radius_maximum=math.sqrt(.25 ** 2 + 2 * u * u), left_x=0, left_y=u, left_z=u, right_x=.25, right_y=u, right_z=u) return {k: str(v) for k, v in r.items()}, rule def test_valid_interval(self): r, rule = self.fixture() self.assertTrue(interval_check(r, rule)[1]) def test_corrupt_interval(self): for key, value in (("t_right", ".3"), ("negative_mass", "1"), ("u", ".5"), ("minimum_density", "-.5"), ("witness_h", "nan")): with self.subTest(key=key): r, rule = self.fixture() r[key] = value with self.assertRaises(AssertionError): interval_check(r, rule) def test_hash_coverage(self): with tempfile.TemporaryDirectory() as directory: folder = Path(directory) (folder / "case").mkdir() (folder / "case/file.txt").write_text("original") campaign = object.__new__(Campaign) campaign.path, campaign.ledger = folder, {} campaign.runs = {"case": {"output_hashes": {"file.txt": digest(folder / "case/file.txt")}}} campaign.artifact("case", "file.txt") (folder / "case/file.txt").write_text("changed") with self.assertRaises(AssertionError): campaign.artifact("case", "file.txt") with self.assertRaises(AssertionError): campaign.artifact("case", "unrecorded.txt") def clean(item): if isinstance(item, dict): return {str(k): clean(v) for k, v in item.items()} if isinstance(item, list): return [clean(v) for v in item] return None if isinstance(item, float) and not math.isfinite(item) else item def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("campaign", type=Path, nargs="?") parser.add_argument("--original", type=Path, help="Original negative-density campaign; defaults to sibling campaign") parser.add_argument("--output", type=Path) parser.add_argument("--self-test", action="store_true") args = parser.parse_args() if args.self_test: result = unittest.TextTestRunner(verbosity=2).run(unittest.defaultTestLoader.loadTestsFromTestCase(SelfTests)) if not result.wasSuccessful(): raise SystemExit(1) if args.campaign is None: return if args.campaign is None or args.output is None: parser.error("campaign and --output are required for artifact auditing") current = Campaign(args.campaign) original = Campaign(args.original or args.campaign.resolve().parent / "campaign") require(current.path != original.path, "A distinct source campaign is required") require(current.provenance["mesh_sha256"] == original.provenance["mesh_sha256"], "Campaign meshes differ") require(current.cases.keys() <= original.cases.keys(), "Region case has no original source case") report = {"schema": "negative_density_regions_audit_v1", "artifact_passed": True, "auditor_sha256": digest(Path(__file__)), "producer_commit": current.provenance["commit"], "original_commit": original.provenance["commit"], "scope": "Independent saved quadratic algebra, aggregation and source binding; not independent mapped FEM integration or global positivity certification.", "tolerances": {"aggregation_relative": 3e-11, "aggregation_absolute": 1e-25, "root_coefficient_roundoff_factor": 4096, "whole_measure_relative": 2e-9, "whole_measure_absolute": 1e-13}, "cases": {name: audit_case(current, original, name) for name in current.cases}} report["consumed_sha256"] = {**original.ledger, **current.ledger} args.output.write_text(json.dumps(clean(report), indent=2, allow_nan=False) + "\n") print("PASS", args.output) if __name__ == "__main__": main()