diff --git a/CMakeLists.txt b/CMakeLists.txt index 4254996..2301996 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -395,6 +395,7 @@ add_executable(tests tests/solver/gravity_reduction.cpp tests/solver/coupled_trial_correction.cpp tests/solver/stellar_numerical_strategies.cpp + tests/solver/continuation_snapshot.cpp ) target_link_libraries(tests PRIVATE mean_field test_mod Catch2::Catch2 Boost::boost) diff --git a/experiments/COUPLED_CONVERGENCE_STUDY_2026-09-19.md b/experiments/COUPLED_CONVERGENCE_STUDY_2026-09-19.md index fbba293..a5ffd8a 100644 --- a/experiments/COUPLED_CONVERGENCE_STUDY_2026-09-19.md +++ b/experiments/COUPLED_CONVERGENCE_STUDY_2026-09-19.md @@ -113,5 +113,15 @@ The exact continuum statement must not be applied uncritically to the present discrete equations: independent rho and P(h), weak EOS consistency, negative density, gravity approximation and central forcing can bypass its assumptions. Numerical Lane–Emden integration also makes the driver's K only -approximately critical. A small discrete virial in such a case would require -careful interpretation. No mass/K change has been made in this campaign. +approximately critical. Independent Lane–Emden integration gives +M_LE=0.9999999918534539 for the recorded K, allowing only +T/|W| <= 2.7155e-9 under those continuum assumptions. Thus the actual rounded +parameters permit negligible, rather than literally zero, rotational support. +A small discrete virial in such a case would require careful interpretation. + +The initial screen keeps its original mass/K parameters. A separate physical +family control is specified in [N3_COMPATIBLE_FAMILY_2026-09-19.md](N3_COMPATIBLE_FAMILY_2026-09-19.md) +and isolated on `research/n3-compatible-family-20260919`: reduce K by one +percent, leave M=1, omit the central border and prescribe J=.035. This tests +a necessary physical compatibility condition, not a preconditioner speedup +or a guarantee of equilibrium or a particular fraction of critical rotation. diff --git a/experiments/check_n3_mass_compatibility.py b/experiments/check_n3_mass_compatibility.py new file mode 100644 index 0000000..b839e86 --- /dev/null +++ b/experiments/check_n3_mass_compatibility.py @@ -0,0 +1,85 @@ +#!/usr/bin/env python3 +"""Independent n=3 Lane–Emden integration; no FEM or production code changes.""" +from fractions import Fraction +from pathlib import Path +import json +import math +import sys + +import scipy +from scipy.integrate import solve_ivp + + +def central_series(x, terms=10): + # Recurrence from theta'' + 2 theta'/x + theta^3 = 0 avoids + # a singular initial-value evaluation and a hand-transcribed series. + coeff = [Fraction(1)] + for k in range(terms - 1): + cube = sum( + coeff[i] * coeff[j] * coeff[k-i-j] + for i in range(k+1) for j in range(k-i+1) + ) + coeff.append(-cube / ((2*k+2)*(2*k+3))) + theta = math.fsum(float(a)*x**(2*k) for k, a in enumerate(coeff)) + mu = -math.fsum(2*k*float(a)*x**(2*k+1) for k, a in enumerate(coeff) if k) + return theta, mu + + +def solve(rtol, atol, maximum_step, start): + def rhs(x, y): + theta, mu = y + return -mu / x**2, x**2 * theta**3 + + def zero(x, y): + return y[0] + zero.terminal = True + zero.direction = -1 + + result = solve_ivp(rhs, (start, 10.0), central_series(start), + events=zero, method="DOP853", rtol=rtol, atol=atol, + max_step=maximum_step) + if not result.success or len(result.t_events[0]) != 1: + raise RuntimeError(result.message) + xi = float(result.t_events[0][0]) + mu = float(result.y_events[0][0, 1]) + kcrit = math.pi / (4.0*math.pi*mu)**(2.0/3.0) + rho_for_unit_radius = xi**3/(4.0*math.pi*mu) + mass = 4.0*math.pi*(driver_k/math.pi)**1.5*mu + radius = xi*math.sqrt(driver_k/math.pi)*driver_rho**(-1.0/3.0) + return dict(rtol=rtol, atol=atol, maximum_step=maximum_step, start=start, + steps=len(result.t)-1, function_evaluations=result.nfev, + xi1=xi, mu1=mu, K_for_mass_one=kcrit, + rho_c_for_mass_and_radius_one=rho_for_unit_radius, + M_LE_at_driver_K=mass, M_target_minus_M_LE=1.0-mass, + R_LE_at_driver_K_and_rho_c=radius, + K_driver_relative_to_critical_minus_one=driver_k/kcrit-1.0, + virial_bound_T_over_absW_max=0.5*(1.0-mass**(2.0/3.0))) + + +metadata = Path(sys.argv[1]).resolve() +values = dict(line.split("=", 1) for line in metadata.read_text().splitlines() if "=" in line) +if float(values["n"]) != 3.0 or float(values["M"]) != 1.0: + raise ValueError("This diagnostic uses the campaign's n=3, M=1, G=1 convention") +driver_k = float(values["K"]) +driver_rho = float(values["rho_c"]) +checks = [solve(3e-12, 1e-14, 0.2, 1e-3), + solve(3e-13, 1e-15, 0.1, 1e-3), + solve(3e-14, 1e-16, 0.05, 1e-3), + solve(3e-14, 1e-16, 0.025, 1e-2), + solve(3e-14, 1e-16, 0.025, 1e-4)] +mu = checks[-1]["mu1"] +def mass_for_k(k): + return 4.0*math.pi*(k/math.pi)**1.5*mu + +out = dict(metadata=str(metadata), python=sys.version, scipy=scipy.__version__, + convention="G=1; theta'= -mu/xi^2; mu'=xi^2 theta^3; first theta zero", + driver={name: values[name] for name in ("n", "K", "rho_c", "M", "J", "xi1", "mu1")}, + checks=checks, + mass_check_spread=max(x["M_LE_at_driver_K"] for x in checks)-min(x["M_LE_at_driver_K"] for x in checks), + full_metadata_K_binary64_ulp=math.ulp(driver_k), + M_LE_for_K_adjacent_binary64=[mass_for_k(math.nextafter(driver_k, -math.inf)), + mass_for_k(math.nextafter(driver_k, math.inf))], + K_quoted_12_decimal_rounding_interval=[0.3639390599515, 0.3639390599525], + M_LE_for_K_quoted_rounding_interval=[mass_for_k(0.3639390599515), mass_for_k(0.3639390599525)], + limitation="Tolerance/step/start sensitivity check, not an interval-arithmetic error certificate.") +print(json.dumps(out, indent=2)) diff --git a/experiments/summarize_coupled_convergence.py b/experiments/summarize_coupled_convergence.py new file mode 100644 index 0000000..ba13bae --- /dev/null +++ b/experiments/summarize_coupled_convergence.py @@ -0,0 +1,684 @@ +#!/usr/bin/env python3 +"""Verify completed campaign evidence and summarize coupled/linear convergence. + +This reads saved artifacts, not FEM operators. A verified record may describe +nonconvergence. Pending records never acquire a final numerical classification. +Run serially with numerical campaigns: hashing and vector checks are real work. +""" +import argparse +from collections import Counter, defaultdict +import csv +import hashlib +import json +import math +from pathlib import Path +import sys +import tempfile +import unittest + +import check_preconditioner_benchmark as benchmark + + +def require(condition, message): + if not condition: + raise ValueError(message) + + +def finite(value): + result = float(value) + require(math.isfinite(result), f'Nonfinite required scalar: {value}') + return result + + +def near(left, right, label, rtol=3e-10, atol=1e-26): + require(math.isfinite(left) and math.isfinite(right) and + abs(left-right) <= atol + rtol * max(abs(left), abs(right)), f'{label}: {left} != {right}') + + +def clean(value): + if isinstance(value, float) and not math.isfinite(value): + return None + if isinstance(value, dict): + return {k: clean(v) for k, v in value.items()} + if isinstance(value, (list, tuple)): + return [clean(v) for v in value] + return value + + +def rows(path): + with path.open() as stream: + return list(csv.DictReader(stream)) + + +def metadata(path): + result = {} + for line in path.read_text().splitlines(): + if not line: + continue + key, separator, value = line.partition('=') + require(separator and key and key not in result, 'Malformed or repeated metadata key') + result[key] = value + return result + + +def norm(values): + return math.sqrt(math.fsum(v*v for v in values)) + + +class Hashes: + def __init__(self): + self.consumed = {} + + def digest(self, path): + path = path.resolve() + key = str(path) + if key not in self.consumed: + digest = hashlib.sha256() + with path.open('rb') as stream: + for chunk in iter(lambda: stream.read(1024*1024), b''): + digest.update(chunk) + self.consumed[key] = digest.hexdigest() + return self.consumed[key] + + def verify(self, path, expected): + require(self.digest(path) == expected, f'Hash mismatch: {path}') + + +def physical(path): + result = {} + for row in rows(path): + q = row['quadrature'] + measurements = result.setdefault(q, {}) + require(row['metric'] not in measurements, 'Repeated physical metric') + measurements[row['metric']] = float(row['value']) + for measurements in result.values(): + required = {'complete_volume_measurement', 'pressure_from_enthalpy_defined_everywhere', + 'invalid_stellar_samples', 'valid_stellar_samples', 'invalid_enthalpy_eos_samples', + 'volume', 'binding_energy', 'force_binding_energy', 'enthalpy_pressure_integral', + 'kinetic_energy', 'enthalpy_virial_error', 'enthalpy_force_virial_error'} + require(required <= measurements.keys(), 'Incomplete physical virial accounting') + for name in ('complete_volume_measurement', 'pressure_from_enthalpy_defined_everywhere'): + require(measurements[name] in (0., 1.), 'Invalid physical validity flag') + for name in ('invalid_stellar_samples', 'valid_stellar_samples', 'invalid_enthalpy_eos_samples'): + count = finite(measurements[name]) + require(count >= 0 and count.is_integer(), 'Invalid physical sample count') + complete = measurements['complete_volume_measurement'] == 1. + if complete: + require(measurements['invalid_stellar_samples'] == 0 and + measurements['valid_stellar_samples'] > 0 and finite(measurements['volume']) > 0, + 'Complete physical measurement contradicts sample validity') + valid_pressure = complete and measurements['invalid_enthalpy_eos_samples'] == 0 + require((measurements['pressure_from_enthalpy_defined_everywhere'] == 1.) == valid_pressure, + 'Enthalpy-pressure validity contradicts sample accounting') + # Independently reconstruct both h-pressure virial scalars. This checks + # their accounting, not the underlying FEM quadrature or physical model. + w, p, t = (measurements[name] for name in + ('binding_energy', 'enthalpy_pressure_integral', 'kinetic_energy')) + if complete: + finite(w) + finite(t) + finite(measurements['force_binding_energy']) + if valid_pressure: + require(finite(p) >= 0, 'Invalid complete enthalpy pressure') + valid_virial = valid_pressure and math.isfinite(w) and w < 0 + for binding, label in [('binding_energy', 'enthalpy_virial_error'), + ('force_binding_energy', 'enthalpy_force_virial_error')]: + if valid_virial: + require(finite(measurements[label]) >= 0, 'Invalid reported virial error') + expected = abs(2*t + measurements[binding] + 3*p) / abs(w) + # The producer accumulates in long double, but serializes each + # cancelling term as double. Scale rounding by those terms. + rounding = 64*sys.float_info.epsilon*(abs(2*t)+abs(measurements[binding])+3*abs(p))/abs(w) + near(measurements[label], expected, label, atol=rounding) + else: + # Invalid physics is a retained numerical outcome. A finite + # virial in this situation would falsely certify missing terms. + require(math.isnan(measurements[label]), 'Undefined physical virial reported as finite') + return result + + +def both_h_virials_below_threshold(measurements, threshold=1e-5): + required = ('volume', 'binding_energy', 'force_binding_energy', 'enthalpy_pressure_integral', + 'kinetic_energy', 'enthalpy_virial_error', 'enthalpy_force_virial_error') + if not all(name in measurements and math.isfinite(measurements[name]) for name in required): + return False + return (measurements.get('complete_volume_measurement') == 1. and + measurements.get('pressure_from_enthalpy_defined_everywhere') == 1. and + measurements.get('invalid_stellar_samples') == 0. and + measurements.get('invalid_enthalpy_eos_samples') == 0. and + measurements.get('valid_stellar_samples', 0.) > 0. and + measurements['volume'] > 0. and measurements['binding_energy'] < 0. and + measurements['enthalpy_pressure_integral'] >= 0. and measurements['kinetic_energy'] >= 0. and + all(0. <= measurements[name] <= threshold for name in + ('enthalpy_virial_error', 'enthalpy_force_virial_error'))) + + +DISPOSITIONS = {0: 'accepted', 1: 'inadmissible_state', 2: 'non_finite_state', + 3: 'non_finite_residual', 4: 'insufficient_decrease', 5: 'auxiliary_solve_failure'} + + +def trial_summary(table, history, initial, check_merit): + grouped = defaultdict(list) + counts = Counter() + for row in table: + iteration, trial, code = (int(row[k]) for k in ('iteration', 'trial', 'disposition')) + alpha = finite(row['alpha']) + require(iteration >= 0 and code in DISPOSITIONS and 0 < alpha <= 1, 'Invalid trial coordinates') + require(trial == len(grouped[iteration]), 'Missing/repeated trial') + if grouped[iteration]: + previous = grouped[iteration][-1] + require(int(previous['disposition']) != 0 and alpha < float(previous['alpha']), + 'Trial after acceptance or nondecreasing backtracking') + if code == 0: + require(row['source'] == '', 'Accepted trial has a rejection source') + if code == 4: + require(row['source'] == 'residual metric', 'Merit failure source mismatch') + grouped[iteration].append(row) + counts[(row['source'], code)] += 1 + require(set(grouped) <= {int(r['iteration']) for r in history}, 'Trial has no completed iteration') + previous = initial + accepted_alphas = [] + for row in history: + trials = grouped[int(row['iteration'])] + require(len(trials) == int(row['trials']), 'Trial count differs from iteration record') + accepted = [r for r in trials if int(r['disposition']) == 0] + require(len(accepted) == int(row['accepted']), 'Trial/iteration acceptance mismatch') + for trial in trials: + code = int(trial['disposition']) + if code in (0, 4): + value, alpha = finite(trial['F']), finite(trial['alpha']) + require(value >= 0 and finite(trial['min_jacobian']) > 0, 'Invalid evaluated trial') + if check_merit: + limit = previous * math.sqrt(1 - 1e-4*alpha) + allowance = 3e-12 * max(previous, value) + require(value <= limit + allowance if code == 0 else value >= limit - allowance, + 'Recorded trial contradicts squared-merit acceptance') + if accepted: + require(accepted[0] is trials[-1], 'Accepted trial must be last') + near(float(accepted[0]['F']), float(row['F']), 'Accepted trial residual') + near(float(accepted[0]['alpha']), float(row['alpha']), 'Accepted alpha') + accepted_alphas.append(float(row['alpha'])) + else: + near(float(row['F']), previous, 'Failed iteration restored residual') + previous = float(row['F']) + return {'total': len(table), 'minimum_accepted_alpha': min(accepted_alphas, default=None), + 'source_disposition_counts': [dict(source=source, disposition=DISPOSITIONS[code], + disposition_code=code, count=count) for (source, code), count in sorted(counts.items())]} + + +def nonlinear(folder, meta): + history, blocks, trials = (rows(folder / name) for name in ('newton.csv', 'blocks.csv', 'trials.csv')) + initial_physical, final_physical = physical(folder/'initial_physical.csv'), physical(folder/'final_physical.csv') + require('18' in initial_physical and '18' in final_physical, 'Missing q18 physical diagnostics') + layout = rows(folder/'layout.csv') + residual_layout = [r for r in layout if r['kind'] == 'residual'] + names = {r['block'] for r in residual_layout} + require(bool(names) and len(names) == len(residual_layout), 'Invalid residual block catalog') + offset = 0 + for block in residual_layout: + require(int(block['offset']) == offset and int(block['size']) > 0, 'Noncontiguous residual layout') + offset += int(block['size']) + block_groups = defaultdict(dict) + for row in blocks: + key = (int(row['iteration']), row['stage']) + require(row['block'] not in block_groups[key], 'Repeated iteration/block observation') + value = finite(row['normalized_norm']) + require(value >= 0, 'Negative block norm') + block_groups[key][row['block']] = value + require(set(block_groups) == {(i, stage) for i in range(len(history)) for stage in ('before','after')}, + 'Incomplete or extra iteration/block stages') + require(all(set(group) == names for group in block_groups.values()), 'Incomplete residual block observations') + initial = finite(initial_physical['18']['normalized_bordered_residual']) + if (0, 'before') in block_groups: + before = norm(block_groups[(0, 'before')].values()) + if meta.get('strategy') in ('reduced', 'positivity-reduced'): + initial = before # Gravity completion precedes the solver's initial metric. + else: + near(before, initial, 'Initial block/physical residual binding') + require(initial >= 0, 'Negative initial norm') + accepted, previous, rates = 0, initial, [] + for i, row in enumerate(history): + require(int(row['iteration']) == i and row['accepted'] in ('0', '1'), 'Invalid iteration sequence') + value = finite(row['F']) + require(value >= 0, 'Negative nonlinear norm') + if (i, 'before') in block_groups: + near(norm(block_groups[(i, 'before')].values()), previous, 'Before block norm') + if (i, 'after') in block_groups: + near(norm(block_groups[(i, 'after')].values()), value, 'After block norm') + for key in ('seconds', 'linear_seconds', 'prepare_seconds', 'refresh_seconds'): + require(finite(row[key]) >= 0, 'Invalid iteration timing') + require(int(row['linear_iterations']) >= 0, 'Invalid Krylov count') + if row['accepted'] == '1': + require(0 < finite(row['alpha']) <= 1, 'Invalid accepted step') + require(row['linear_status'] == '0', 'Accepted Newton step has failed linear solve') + require(0 <= finite(row['linear_true']) <= finite(meta['linear_tol']), + 'Accepted Newton step fails its reported linear gate') + accepted += 1 + rates.append(value/previous if previous > 0 else None) + previous = value + near(previous, finite(final_physical['18']['normalized_bordered_residual']), 'Final residual binding') + require(accepted == int(meta['accepted_steps']) and meta['converged'] in ('0', '1'), 'Invalid completion metadata') + rtol, atol = finite(meta['rtol']), finite(meta['atol']) + require(rtol >= 0 and atol >= 0, 'Invalid nonlinear tolerances') + threshold = max(atol, rtol*initial) + declared = meta['converged'] == '1' + require(not declared or previous <= threshold, 'Convergence flag exceeds actual initial-dependent threshold') + stages = trial_summary(trials, history, initial, meta.get('strategy') in ('baseline', 'curvature', 'reduced')) + final_blocks = block_groups.get((len(history)-1, 'after'), {}) + strategies = {r['metric']: float(r['value']) for r in rows(folder/'numerical_strategies.csv')} + return {'kind': 'newton', 'outcome': 'converged' if declared else 'not_converged', + 'failure': meta.get('failure'), 'initial_residual': initial, 'final_residual': previous, + 'final_over_initial': previous/initial if initial else None, + 'rtol': rtol, 'atol': atol, 'actual_threshold': threshold, + 'final_residual_within_threshold': previous <= threshold, + 'accepted_updates': accepted, 'attempted_updates': len(history), + 'krylov_iterations': sum(int(r['linear_iterations']) for r in history), + 'newton_seconds': math.fsum(float(r['seconds']) for r in history), + 'linear_seconds': math.fsum(float(r['linear_seconds']) for r in history), + 'prepare_seconds': math.fsum(float(r['prepare_seconds']) for r in history), + 'refresh_seconds': math.fsum(float(r['refresh_seconds']) for r in history), + 'accepted_residual_ratios': rates, 'trials': stages, + 'final_block_squared_fractions': {k: v*v/(previous*previous) if previous else None for k,v in final_blocks.items()}, + 'initial_physical': initial_physical, 'final_physical': final_physical, + 'numerical_strategies': strategies, 'history': history, + 'linear_report_scope': 'reduced_or_auxiliary; consult lifted full check' if 'reduced' in meta.get('strategy','') else 'full_system', + 'initial_physical_scope': 'before optional gravity completion; first before-row defines completed initial norm'} + + +def linear_trials(folder, summary): + table = rows(folder/'linear_direction_trials.csv') + rhs = benchmark.vector(folder/'benchmark_normalized_rhs.txt') + f = [-v for v in rhs] + f_norm = norm(f) + result = {} + for role in ('baseline', 'candidate'): + p = benchmark.vector(folder/f'{role}_normalized_correction.txt') + action = benchmark.vector(folder/f'{role}_trial_jacobian_action.txt') + saved_defect = benchmark.vector(folder/f'{role}_normalized_linear_defect.txt') + require(len(f) == len(p) == len(action) == len(saved_defect), 'Trial action vector dimensions differ') + defect = [a+b for a,b in zip(action,f)] + error = norm([a-b for a,b in zip(defect,saved_defect)]) + require(error <= 2e-12*(norm(action)+f_norm), 'Saved full-J defect disagrees with re-applied action') + relative = norm(defect)/f_norm + slope = math.fsum(a*b for a,b in zip(f,action)) + recorded = summary['comparison'][role] + near(relative, float(recorded['true_relative']), 'Linear gate relative defect', atol=error/f_norm) + gate = norm(defect) <= float(recorded['relative_tolerance'])*f_norm + require(gate == (recorded['true_gate_pass'] == '1'), 'Trial direction gate mismatch') + samples = [row for row in table if row['role'] == role] + require([float(r['alpha']) for r in samples] == [1., .25, .0625, .015625, .00390625, .0009765625, .000244140625], + 'Incomplete or changed frozen trial grid') + observations = [] + for row in samples: + alpha = finite(row['alpha']) + for key, expected in [('base_norm', f_norm), ('linear_defect_relative', relative), + ('direction_norm', norm(p)), ('action_norm', norm(action)), ('merit_slope', slope)]: + near(float(row[key]), expected, key) + require(row['prepared'] in ('0','1') and row['decrease_screen'] in ('0','1'), 'Invalid trial flags') + if row['prepared'] == '1': + require(row['reason'] == '-1' and row['stage'] == '-1', 'Prepared trial has rejection codes') + predicted = norm([x+alpha*y for x,y in zip(f,action)]) + trial, remainder, projection = (finite(row[k]) for k in + ('trial_norm','remainder_norm','signed_remainder_projection')) + require(trial >= 0 and remainder >= 0, 'Negative trial/remainder norm') + near(float(row['predicted_norm']), predicted, 'Affine predicted norm') + near(float(row['relative_norm']), trial/f_norm, 'Trial relative norm') + near(float(row['remainder_over_alpha']), remainder/alpha, 'Scaled remainder') + identity = trial*trial - predicted*predicted - 2*projection - remainder*remainder + scale = trial*trial + predicted*predicted + 2*abs(projection) + remainder*remainder + allowance = 3e-10*scale + 1e-28 + require(abs(identity) <= allowance and abs(identity-float(row['merit_identity_error'])) <= allowance, + 'Scalar squared-merit identity contradicts reported norms') + require(abs(projection) <= predicted*remainder + allowance, 'Impossible signed remainder projection') + screen = trial*trial <= (1 - 1e-4*alpha)*f_norm*f_norm + require(screen == (row['decrease_screen'] == '1'), 'Squared-merit screen mismatch') + else: + require(int(row['reason']) >= 0 and int(row['stage']) >= 0 and row['decrease_screen'] == '0', + 'Unprepared trial claims evaluation or merit success') + require(all(math.isnan(float(row[k])) for k in ('trial_norm','relative_norm','remainder_norm', + 'remainder_over_alpha','predicted_norm','signed_remainder_projection','merit_identity_error')), + 'Unprepared trial has fabricated finite measurements') + observations.append(row) + result[role] = {'verified_linear_gate': gate, 'relative_true_defect': relative, 'merit_slope': slope, + 'largest_sampled_decreasing_alpha': max((float(r['alpha']) for r in samples if r['decrease_screen']=='1'), default=None), + 'observations': observations} + require({r['role'] for r in table} == {'baseline','candidate'}, 'Unknown trial role') + invariants = {r['metric']: r['value'] for r in rows(folder/'linear_direction_trial_invariants.csv')} + require(invariants == {'accepted_unchanged':'1','base_residual_restored':'1'}, 'Trial restoration invariant failed') + return {'scope': 'Jp+F and slope reconstructed from saved vectors; trial remainder identities checked only from saved scalars. No trial F/remainder vectors or independent FEM evaluation. Sampled grid is not the live halving line search; failed linear gates remain failed.', + 'roles': result} + + +def linear_history(folder, assessment): + table = rows(folder/'preconditioner_history.csv') + require({r['role'] for r in table} == {'baseline','candidate'}, 'Missing or unknown linear history role') + result = {} + for role in ('baseline','candidate'): + samples = [r for r in table if r['role'] == role] + record = assessment['comparison'][role] + iterations = [int(r['iteration']) for r in samples] + require(iterations == sorted(iterations) and iterations[0] >= 0, 'Nonsequential recurrence history') + require(all(r['variant'] == record['variant'] and r['final'] in ('0','1') for r in samples), + 'Linear history variant/flag mismatch') + require(all(finite(r['reported_norm']) >= 0 for r in samples), 'Invalid recurrence norm') + result[role] = {'samples':samples, 'minimum_reported_over_rhs': + min(float(r['reported_norm']) for r in samples)/float(record['rhs_norm']), + 'scope':'Recorded Krylov recurrence norms; not independently recomputed true residuals or a proved floor.'} + return result + + +def frozen_audit_scope(folder, assessment): + variants = {'surface-border-schur', 'gravity-schur-surface-border-schur'} + roles = [role for role in ('baseline', 'candidate') + if assessment['comparison'][role]['variant'] in variants] + files = ('preconditioner_surface_border_schur.csv', + 'preconditioner_surface_border_coordinates.csv', + 'preconditioner_surface_border_matrices.csv') + if roles: + require(all((folder/name).is_file() for name in files), 'Missing surface/scalar Schur observations') + return { + 'endpoint_scope': 'Saved full-system linear defects, block norms, gates and reported timing sums are audited.', + 'surface_border_schur_roles': roles, + 'surface_border_schur_internal_audit': 'not_performed' if roles else 'not_applicable', + 'surface_border_schur_artifact_scope': + 'Presence and producer hashes only; coarse matrices, coordinates, counters and condition estimates are not independently checked.' + if roles else None, + 'operator_work_scope': + 'Outer Krylov J counts exclude one internal J action per surface/scalar Schur application; inverse time includes it. No total operator-work comparison is certified.' + if roles else 'Application counters have the scope stated by the producing benchmark.'} + + +def verify_run_binding(campaign, spec, record, provenance): + command = record['command'] + require(command[1:-4] == spec['arguments'], 'Command differs from declared case') + require(command[-4] == '--mesh' and command[-2] == '--output', 'Command has unexpected argument layout') + require(Path(command[0]).resolve() == campaign/'plateau_diagnostic', 'Case did not use the retained executable') + require(Path(command[-3]).resolve() == Path(provenance['mesh']).resolve(), 'Case mesh differs from provenance') + require(Path(command[-1]).resolve() == campaign/spec['id'], 'Command output is not the case directory') + inputs = {str(Path(p).resolve()) for p in record['input_state_hashes']} + arguments = spec['arguments'] + for i, word in enumerate(arguments): + if word not in ('--probe-state','--probe-direction','--replay-candidate','--density-snapshot','--restart-snapshot'): + continue + path = Path(arguments[i+1]) + path = (path if path.is_absolute() else Path(provenance['source'])/path).resolve() + if word in ('--density-snapshot','--restart-snapshot'): + snapshot = path/'accepted_state.txt' if path.is_dir() else path + prior = snapshot.parent + if not (prior/'metadata.txt').is_file(): + prior = prior.parent + expected = {snapshot, prior/'metadata.txt', prior/'layout.csv', + prior.parent/'provenance.json', prior.parent/(prior.name+'_run.json')} + require({str(p) for p in expected} <= inputs, 'Restart input metadata/provenance is not hash-bound') + prior_run = json.loads((prior.parent/(prior.name+'_run.json')).read_text()) + prior_provenance = json.loads((prior.parent/'provenance.json').read_text()) + require(prior_provenance['mesh_sha256'] == provenance['mesh_sha256'], 'Restart mesh hash differs') + for saved in (snapshot, prior/'metadata.txt', prior/'layout.csv'): + require(prior_run['output_hashes'][str(saved.relative_to(prior))] == record['input_state_hashes'][str(saved)], + 'Restart input differs from its producing record') + else: + require(str(path) in inputs, 'Supplied input vector has no recorded hash') + + +def audit_case(campaign, spec, hashes, provenance=None): + name, folder = spec['id'], campaign/spec['id'] + run_path = campaign/f'{name}_run.json' + result = {'id': name, 'campaign': str(campaign), 'complete': False, 'artifact_passed': None} + if not run_path.exists(): + return {**result, 'outcome': 'not_started'} + try: + record = json.loads(run_path.read_text()) + except json.JSONDecodeError: + return {**result, 'outcome': 'incomplete_run_record'} + if 'finished_utc' not in record: + return {**result, 'outcome': 'running_or_interrupted'} + result.update(complete=True, run=record) + try: + hashes.digest(run_path) + require(record['id'] == name, 'Run/case ID mismatch') + verify_run_binding(campaign, spec, record, provenance) + listed = record['output_hashes'] + actual = {str(p.relative_to(folder)) for p in folder.rglob('*') if p.is_file() and p.suffix in ('.txt','.csv','.json')} + require(set(listed) == actual, 'Retained output set differs from run record') + for relative, expected in listed.items(): + path = folder/relative + require(path.resolve().is_relative_to(folder), 'Output hash path escapes case directory') + hashes.verify(path, expected) + require(record['inputs_unchanged'] is True, 'Runner recorded mutated inputs') + for path, expected in record['input_state_hashes'].items(): + hashes.verify(Path(path), expected) + if record.get('timed_out') or record['exit_code'] not in (0,3): + return {**result, 'artifact_passed': True, 'kind':'execution_failure', 'outcome':'timed_out' if record.get('timed_out') else 'execution_failure'} + meta = metadata(folder/'metadata.txt') + result['metadata'] = meta + if (folder/'preconditioner_summary.csv').exists(): + assessment = benchmark.inspect_case(folder) + require(assessment['passed'], 'Frozen benchmark audit failed: ' + str([c for c in assessment['checks'] if not c['passed']])) + result.update(kind='benchmark', outcome=assessment['comparison']['outcome'], benchmark=assessment) + result['benchmark_audit_scope'] = frozen_audit_scope(folder, assessment) + result['linear_history'] = linear_history(folder, assessment) + if (folder/'linear_direction_trials.csv').exists(): + result['linear_trials'] = linear_trials(folder, assessment) + else: + result.update(nonlinear(folder, meta)) + require((record['exit_code'] == 0) == (result['outcome'] == 'converged'), 'Exit code/convergence mismatch') + result['artifact_passed'] = True + except (OSError, ValueError, KeyError, TypeError, IndexError, OverflowError, ZeroDivisionError) as error: + result.update(artifact_passed=False, outcome='invalid_artifacts', error=str(error)) + return result + + +def audit_campaign(path, hashes): + campaign = path.resolve() + provenance_path = campaign/'provenance.json' + provenance = json.loads(provenance_path.read_text()) + hashes.digest(provenance_path) + if 'binary_sha256' not in provenance: + return {'path':str(campaign), 'complete':False, 'artifact_passed':None, 'outcome':'build_incomplete', 'cases':[]} + hashes.verify(campaign/'plateau_diagnostic', provenance['binary_sha256']) + hashes.verify(Path(provenance['mesh']), provenance['mesh_sha256']) + hashes.digest(campaign/'cases.json') + require(json.loads((campaign/'cases.json').read_text()) == provenance['cases'], 'Cases file differs from provenance') + ids = [spec['id'] for spec in provenance['cases']] + require(bool(ids) and len(ids) == len(set(ids)) and all(Path(name).name == name and name not in ('.','..') for name in ids), + 'Invalid or repeated case IDs') + cases = [audit_case(campaign, spec, hashes, provenance) for spec in provenance['cases']] + complete = all(c['complete'] for c in cases) + return {'path':str(campaign), 'producer_commit':provenance['commit'], 'mesh_sha256':provenance['mesh_sha256'], + 'complete':complete, 'artifact_passed':complete and all(c['artifact_passed'] for c in cases), 'cases':cases} + + +def flat_rows(campaigns): + result = [] + for campaign in campaigns: + for case in campaign['cases']: + meta, run = case.get('metadata',{}), case.get('run',{}) + row = {k:case.get(k) for k in ('id','kind','outcome','complete','artifact_passed','failure','initial_residual', + 'final_residual','final_over_initial','rtol','atol','actual_threshold','final_residual_within_threshold', + 'accepted_updates','attempted_updates','krylov_iterations', + 'newton_seconds','linear_seconds','prepare_seconds','refresh_seconds')} + row.update(campaign=campaign['path'], **{k:meta.get(k) for k in + ('n','K','rho_c','M','J','n3_K_factor','seed_analytic_mass_estimate','central_phase', + 'density_order','surface_prescription','surface_maximum_order','surface_degree', + 'strategy','preconditioner_variant')}) + row.update({k:run.get(k) for k in ('seconds','utc_elapsed_seconds','utc_minus_monotonic_seconds','child_cpu_seconds','exit_code')}) + q18 = case.get('final_physical',{}).get('18',{}) + for key in ('enthalpy_virial_error','enthalpy_force_virial_error','normalized_unbordered_residual', + 'central_border','normalized_border_action','hydrostatic_gradient_rms','mass_absolute_error','angular_momentum_absolute_error'): + row[key] = q18.get(key) + row['both_h_virials_below_1e-5'] = both_h_virials_below_threshold(q18) if q18 else None + row['minimum_accepted_alpha'] = case.get('trials',{}).get('minimum_accepted_alpha') + for stage in ('pressure force','residual metric'): + row[stage.replace(' ','_')+'_rejections'] = sum(r['count'] for r in + case.get('trials',{}).get('source_disposition_counts',[]) if r['source'] == stage and r['disposition_code'] != 0) + if 'benchmark' in case: + row['surface_border_schur_internal_audit'] = case['benchmark_audit_scope']['surface_border_schur_internal_audit'] + for role in ('baseline','candidate'): + summary = case['benchmark']['comparison'][role] + for key in ('variant','true_relative','true_gate_pass','iterations','krylov_seconds','construction_seconds','setup_solve_and_true_check_seconds'): + row[role+'_'+key] = summary[key] + result.append(row) + return result + + +class TinyChecks(unittest.TestCase): + def test_pending_is_not_success(self): + with tempfile.TemporaryDirectory() as directory: + result = audit_case(Path(directory), {'id':'pending','arguments':[]}, Hashes()) + self.assertFalse(result['complete']) + self.assertIsNone(result['artifact_passed']) + + def test_squared_merit_and_failed_stage(self): + trials = [dict(iteration='0',trial='0',alpha='1',disposition='1',source='pressure force',F='nan',min_jacobian='1'), + dict(iteration='0',trial='1',alpha='.5',disposition='0',source='',F='.8',min_jacobian='1')] + history = [dict(iteration='0',trials='2',accepted='1',F='.8',alpha='.5')] + result = trial_summary(trials, history, 1., True) + self.assertEqual(result['minimum_accepted_alpha'], .5) + self.assertEqual(len(result['source_disposition_counts']), 2) + trials[-1]['F'] = '1.1' + with self.assertRaises(ValueError): + trial_summary(trials, history, 1., True) + + def test_duplicate_metadata_and_hash_corruption(self): + with tempfile.TemporaryDirectory() as directory: + path = Path(directory)/'metadata.txt' + path.write_text('n=1\nn=2\n') + with self.assertRaises(ValueError): metadata(path) + with self.assertRaises(ValueError): Hashes().verify(path, '0'*64) + + def test_command_must_use_retained_binary_and_mesh(self): + with tempfile.TemporaryDirectory() as directory: + campaign = Path(directory).resolve() + spec = dict(id='case', arguments=['--n','2']) + provenance = dict(source=directory, mesh=str(campaign/'mesh')) + record = dict(command=[str(campaign/'plateau_diagnostic'),'--n','2','--mesh', + str(campaign/'mesh'),'--output',str(campaign/'case')], input_state_hashes={}) + verify_run_binding(campaign, spec, record, provenance) + record['command'][0] = str(campaign/'other_binary') + with self.assertRaises(ValueError): verify_run_binding(campaign, spec, record, provenance) + + def test_virial_cancellation_and_corruption(self): + with tempfile.TemporaryDirectory() as directory: + path = Path(directory)/'physical.csv' + values = dict(complete_volume_measurement=1., pressure_from_enthalpy_defined_everywhere=1., + invalid_stellar_samples=0., valid_stellar_samples=8., invalid_enthalpy_eos_samples=0., volume=1., + binding_energy=-3., force_binding_energy=-3.+1e-8, + enthalpy_pressure_integral=1., kinetic_energy=0., enthalpy_virial_error=1e-17, + enthalpy_force_virial_error=1e-8/3.) + def save(): + path.write_text('quadrature,metric,value\n' + ''.join(f'18,{k},{v}\n' for k,v in values.items())) + save() + self.assertTrue(both_h_virials_below_threshold(physical(path)['18'])) + values['enthalpy_force_virial_error'] = 1e-5 + save() + with self.assertRaises(ValueError): physical(path) + + def test_invalid_physics_is_retained_without_virial_success(self): + with tempfile.TemporaryDirectory() as directory: + path = Path(directory)/'physical.csv' + values = dict(complete_volume_measurement=0., pressure_from_enthalpy_defined_everywhere=0., + invalid_stellar_samples=1., valid_stellar_samples=7., invalid_enthalpy_eos_samples=0., + volume=math.nan, binding_energy=math.nan, force_binding_energy=math.nan, + enthalpy_pressure_integral=math.nan, kinetic_energy=math.nan, + enthalpy_virial_error=math.nan, enthalpy_force_virial_error=math.nan) + def save(): + path.write_text('quadrature,metric,value\n' + ''.join(f'18,{k},{v}\n' for k,v in values.items())) + save() + self.assertFalse(both_h_virials_below_threshold(physical(path)['18'])) + values['enthalpy_virial_error'] = values['enthalpy_force_virial_error'] = 0. + self.assertFalse(both_h_virials_below_threshold(values)) + save() + with self.assertRaises(ValueError): physical(path) + del values['binding_energy'] + save() + with self.assertRaises(ValueError): physical(path) + + def test_flat_summary_retains_physical_family(self): + meta = dict(n='3', K='.3603', rho_c='12.74', M='1', J='.035', n3_K_factor='.99', + central_phase='0', surface_prescription='harmonics', surface_maximum_order='4') + result = flat_rows([dict(path='/example', cases=[dict(metadata=meta)])])[0] + for name, value in meta.items(): + self.assertEqual(result[name], value) + + def test_new_schur_scope_is_explicit_and_requires_retained_observations(self): + with tempfile.TemporaryDirectory() as directory: + folder = Path(directory) + assessment = {'comparison': {'baseline': {'variant':'independent'}, + 'candidate': {'variant':'surface-border-schur'}}} + with self.assertRaises(ValueError): frozen_audit_scope(folder, assessment) + for name in ('preconditioner_surface_border_schur.csv', + 'preconditioner_surface_border_coordinates.csv', + 'preconditioner_surface_border_matrices.csv'): + (folder/name).write_text('retained,observation\n') + scope = frozen_audit_scope(folder, assessment) + self.assertEqual(scope['surface_border_schur_internal_audit'], 'not_performed') + self.assertEqual(scope['surface_border_schur_roles'], ['candidate']) + assessment['comparison']['candidate']['variant'] = 'gravity-schur' + self.assertEqual(frozen_audit_scope(folder, assessment)['surface_border_schur_internal_audit'], + 'not_applicable') + + def test_completed_failure_requires_matching_retained_outputs(self): + with tempfile.TemporaryDirectory() as directory: + campaign = Path(directory).resolve() + folder = campaign/'case' + folder.mkdir() + spec = dict(id='case', arguments=[]) + provenance = dict(source=directory, mesh=str(campaign/'mesh')) + record = dict(id='case', finished_utc='2026-09-19T00:00:00+00:00', exit_code=1, + command=[str(campaign/'plateau_diagnostic'),'--mesh',str(campaign/'mesh'),'--output',str(folder)], + output_hashes={}, input_state_hashes={}, inputs_unchanged=True) + (campaign/'case_run.json').write_text(json.dumps(record)) + result = audit_case(campaign, spec, Hashes(), provenance) + self.assertTrue(result['artifact_passed']) + self.assertEqual(result['outcome'],'execution_failure') + (folder/'unrecorded.csv').write_text('value\n1\n') + self.assertFalse(audit_case(campaign, spec, Hashes(), provenance)['artifact_passed']) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument('campaigns', nargs='*', type=Path) + parser.add_argument('--output', type=Path) + parser.add_argument('--self-test', action='store_true') + args = parser.parse_args() + if args.self_test: + suite = unittest.defaultTestLoader.loadTestsFromTestCase(TinyChecks) + return 0 if unittest.TextTestRunner(verbosity=2).run(suite).wasSuccessful() else 1 + if not args.campaigns or args.output is None: + parser.error('Provide campaign directories and --output analysis.json') + hashes = Hashes() + campaigns = [] + for path in args.campaigns: + try: + campaigns.append(audit_campaign(path, hashes)) + except (OSError, ValueError, KeyError, TypeError, IndexError) as error: + campaigns.append({'path':str(path.resolve()), 'complete':False, 'artifact_passed':False, + 'outcome':'invalid_campaign_provenance', 'error':str(error), 'cases':[]}) + result = {'artifact_passed':all(c['artifact_passed'] for c in campaigns), + 'complete':all(c['complete'] for c in campaigns), 'campaigns':campaigns, + 'auditor_sha256':hashes.digest(Path(__file__)), + 'benchmark_checker_sha256':hashes.digest(Path(benchmark.__file__)), 'consumed_sha256':hashes.consumed, + 'scope':'Provenance and retained-file hashes; recorded Newton/trial/block bindings; saved-vector final linear defects and scalar trial identities. Numerical failure is an outcome, not an artifact failure.', + 'limitations':['No independent FEM reassembly or physical-validity certificate.', + 'Linear histories are recurrence residuals; only saved endpoints have independently checked true defects. No proved residual floor.', + 'New surface/scalar Schur wrapper artifacts are retained and hash-bound, but their matrices, coordinates, counters and conditioning are outside this audit; endpoint gates and timing sums are checked.', + 'Trial remainder vectors are not saved. Scalar checks cannot verify the underlying FEM evaluations.', + 'Virial outcomes refer to the reported quadrature and h-derived pressure; they do not establish nonlinear convergence.', + 'Normalized residuals and rtol thresholds belong to their own discrete systems; Krylov counts are not equal work units.', + 'Newton time excludes observers/setup; benchmark setup/solve costs have separate scope. Clock discrepancies are retained without rescaling.']} + args.output.parent.mkdir(parents=True, exist_ok=True) + args.output.write_text(json.dumps(clean(result), indent=2, allow_nan=False)+'\n') + flat = flat_rows(campaigns) + with args.output.with_suffix('.csv').open('w') as stream: + writer = csv.DictWriter(stream, fieldnames=list(dict.fromkeys(k for row in flat for k in row))) + writer.writeheader() + writer.writerows(clean(flat)) + print(json.dumps({'complete':result['complete'], 'artifact_passed':result['artifact_passed'], + 'cases':len(flat), 'output':str(args.output)})) + invalid = any(c.get('outcome') == 'invalid_campaign_provenance' or + any(case['artifact_passed'] is False for case in c['cases']) for c in campaigns) + return 1 if invalid else 0 if result['artifact_passed'] else 2 + + +if __name__ == '__main__': + sys.exit(main())