"""Full-residual 17/33/65 audit for the AP1-R2c lambda=1 DAE branch.""" from __future__ import annotations from dataclasses import dataclass import json from pathlib import Path import numpy as np from numpy.polynomial.chebyshev import chebfit, chebval from scipy.interpolate import BarycentricInterpolator, make_interp_spline from scipy.optimize import minimize_scalar from ap1_r2c_implicit_dae_preflight import positive_friedmann_root from ap1_r2c_implicit_dae_solver import ( ImplicitDAEResidual, Q_KEYS, X_KEYS, _Dt, ) from ap1_r2c_self_consistent_candidate import ( classical_terms, portal_terms, standard_stress, ) G22_THRESHOLDS = { "max_abs_friedmann_residual": 1.0e-6, "ward_normalized": 2.0e-6, "validation_source_relative_change": 1.0e-4, "validation_wronskian_relative_error": 1.0e-11, } DAE_RESIDUAL_GATE = 1.0e-4 @dataclass class LiftCandidate: label: str vector: np.ndarray norm: float def load_lambda1_checkpoint(report_path: Path) -> np.ndarray: report = json.loads(report_path.read_text(encoding="utf-8")) stages = [stage for stage in report.get("stages", []) if stage.get("converged") and float(stage.get("lambda", -1.0)) == 1.0] if not stages: raise ValueError("converged lambda=1 checkpoint required") vector = np.asarray(stages[-1].get("solution_vector", []), dtype=float) if vector.shape != (120,): raise ValueError("complete 17-node lambda=1 checkpoint required") return vector def _physical_diagnostics(system: ImplicitDAEResidual, vector: np.ndarray) -> dict: residual = system.evaluate(vector) X, Q, N_seed = system.unpack(vector) states_without_H = np.column_stack([X[key] for key in X_KEYS]) rho_cl = np.array([ classical_terms(np.r_[row, system.seed[4]])["rho"] for row in states_without_H ]) rho_std, _ = standard_stress(N_seed + system.Nrel, system.pilot) total = rho_cl + rho_std + Q["rho"] H = positive_friedmann_root(total) Hdot = _Dt(H, H, system.D_N) fresh, wronskian = system._sources(X, H, Hdot) friedmann = (3.0 * H**2 - total) / max( float(np.max(np.abs(total))), 1.0e-300 ) mass2 = np.asarray(portal_terms( X["sigma"], X["theta"], system.chi_p )[0]) rho_dot = _Dt(fresh["rho"], H, system.D_N) mass2_dot = _Dt(mass2, H, system.D_N) ward_terms = ( rho_dot, 3.0 * H * (fresh["rho"] + fresh["pressure"]), -0.5 * mass2_dot * fresh["chi2"], ) ward = sum(ward_terms) ward_scale = max( *(float(np.max(np.abs(term))) for term in ward_terms), 1.0e-300 ) source_change = max( float(np.max(np.abs(fresh[key] - Q[key])) / max(np.max(np.abs(fresh[key])), 1.0e-300)) for key in Q_KEYS ) values = { "dae_max_abs_scaled_residual": float(np.max(np.abs(residual))), "max_abs_friedmann_residual": float(np.max(np.abs(friedmann))), "ward_normalized": float(np.max(np.abs(ward)) / ward_scale), "validation_source_relative_change": source_change, "validation_wronskian_relative_error": wronskian, "finite": bool(np.all(np.isfinite(vector)) and np.all(np.isfinite(residual))), "positive_H": bool(system.last_blocks["min_H"] > 0.0), "min_H": system.last_blocks["min_H"], "min_PX": system.last_blocks["min_PX"], "min_K": system.last_blocks["min_K"], } values["gate_pass"] = { "dae_residual": values["dae_max_abs_scaled_residual"] <= DAE_RESIDUAL_GATE, **{ key: values[key] <= threshold for key, threshold in G22_THRESHOLDS.items() }, "structural": bool( values["finite"] and values["positive_H"] and values["min_PX"] > 0.0 and values["min_K"] > 0.0 ), } return values def _lift_arrays(source: ImplicitDAEResidual, vector: np.ndarray, target: ImplicitDAEResidual, method: str, degree: int | None = None) -> np.ndarray: X, Q, N_seed = source.unpack(vector) if method == "barycentric": lift = lambda values: BarycentricInterpolator( source.Nrel, values )(target.Nrel) elif method == "chebyshev": x_source = 2.0 * source.Nrel / source.cfg.delta_N - 1.0 x_target = 2.0 * target.Nrel / target.cfg.delta_N - 1.0 lift = lambda values: chebval( x_target, chebfit(x_source, values, int(degree)) ) elif method == "bspline": if degree is None or not 1 <= int(degree) <= 5: raise ValueError("bounded B-spline degree required") lift = lambda values: make_interp_spline( source.Nrel, values, k=int(degree) )(target.Nrel) else: raise ValueError("recognized spectral lift required") return target.pack( {key: lift(X[key]) for key in X_KEYS}, {key: lift(Q[key]) for key in Q_KEYS}, N_seed, ) def endpoint_bspline_lift(source: ImplicitDAEResidual, vector: np.ndarray, target: ImplicitDAEResidual, degree: int = 2) -> np.ndarray: """Lift one validated parent while preserving both interval endpoints.""" if target.cfg.candidate_nodes != source.cfg.candidate_nodes + 1: raise ValueError("single-level p continuation required") lifted = _lift_arrays(source, vector, target, "bspline", degree) source_X, source_Q, source_N = source.unpack(vector) target_X, target_Q, target_N = target.unpack(lifted) for key in X_KEYS: if not np.allclose(target_X[key][[0, -1]], source_X[key][[0, -1]], rtol=0.0, atol=1.0e-14): raise RuntimeError("B-spline lift did not preserve X endpoints") for key in Q_KEYS: scale = max(float(np.max(np.abs(source_Q[key]))), 1.0e-300) if not np.allclose(target_Q[key][[0, -1]], source_Q[key][[0, -1]], rtol=0.0, atol=1.0e-13 * scale): raise RuntimeError("B-spline lift did not preserve source endpoints") if target_N != source_N: raise RuntimeError("B-spline lift changed the seed position") return lifted def richardson_bspline_lift(source: ImplicitDAEResidual, vector: np.ndarray, target: ImplicitDAEResidual, beta: float) -> np.ndarray: """Endpoint-preserving extrapolation of the degree-4/5 B-spline lifts.""" if not np.isfinite(beta) or not -3.0 <= beta <= 5.0: raise ValueError("bounded finite Richardson factor required") degree4 = endpoint_bspline_lift(source, vector, target, degree=4) degree5 = endpoint_bspline_lift(source, vector, target, degree=5) return degree5 + float(beta) * (degree5 - degree4) def degree_pair_bspline_lift(source: ImplicitDAEResidual, vector: np.ndarray, target: ImplicitDAEResidual, lower_degree: int, upper_degree: int, beta: float) -> np.ndarray: """Endpoint-preserving bounded extrapolation of two spline degrees.""" if not (1 <= lower_degree < upper_degree <= 5): raise ValueError("ordered B-spline degrees in 1..5 required") if not np.isfinite(beta) or not -3.0 <= beta <= 5.0: raise ValueError("bounded finite degree-pair factor required") lower = endpoint_bspline_lift( source, vector, target, degree=lower_degree ) upper = endpoint_bspline_lift( source, vector, target, degree=upper_degree ) return upper + float(beta) * (upper - lower) def select_p_lift_projection_bank( source: ImplicitDAEResidual, vector: np.ndarray, target: ImplicitDAEResidual, betas: tuple[float, ...] = ( 1.5, 1.65, 1.75, 1.85, 2.0, 2.15, 2.2, ), ) -> tuple[np.ndarray, dict]: """Select only after full residual checks of bounded p-lift candidates.""" candidates = [] def register(label: str, lifted: np.ndarray) -> None: projected, projection = balanced_algebraic_projection(target, lifted) residual = target.evaluate(projected) candidates.append(( float(np.max(np.abs(residual))), label, projected, projection, )) register( "endpoint_quadratic_bspline", endpoint_bspline_lift(source, vector, target, degree=2), ) register( "endpoint_quintic_bspline", endpoint_bspline_lift(source, vector, target, degree=5), ) for beta in betas: register( f"endpoint_richardson_degree4_5_beta_{beta:.6g}", richardson_bspline_lift(source, vector, target, beta), ) # A bounded cross-degree audit at the validated 37-node parent isolated # the adjacent degree-3/4 family. Keep only its PASS-bracketing factors; # this is a fixed method extension, not an adaptive gate search. for beta in (3.5, 4.5, 5.0): register( f"endpoint_degree3_4_beta_{beta:.6g}", degree_pair_bspline_lift( source, vector, target, 3, 4, beta ), ) norm, label, chosen, projection = min(candidates, key=lambda item: item[0]) return chosen.copy(), { "method": "bounded_full_residual_checked_endpoint_p_lift_bank", "selected": label, "selected_norm": norm, "candidate_norms": {item[1]: item[0] for item in candidates}, "selected_projection": projection, } def balanced_algebraic_projection(system: ImplicitDAEResidual, vector: np.ndarray) -> tuple[np.ndarray, dict]: """Balance the algebraic chi2 closure against its coupled background block. Pressure is exactly algebraic and does not enter any other residual block. At fixed X, rho and N_seed, chi2 enters the background equations linearly while its freshly evaluated source is independent of Q_chi2. The complete residual along the bounded chi2 projection segment is therefore affine. We minimize its infinity norm on alpha in [0,1] and independently verify both the affine prediction and the full physical residual. """ vector = np.asarray(vector, dtype=float) n = system.cfg.candidate_nodes pressure = slice(5 * n, 6 * n) chi2 = slice(6 * n, 7 * n) raw_residual = system.evaluate(vector) pressure_projected = vector.copy() pressure_projected[pressure] -= raw_residual[pressure] residual0 = system.evaluate(pressure_projected) fully_projected = pressure_projected.copy() fully_projected[chi2] -= residual0[chi2] residual1 = system.evaluate(fully_projected) def predicted_norm(alpha: float) -> float: residual = (1.0 - alpha) * residual0 + alpha * residual1 return float(np.max(np.abs(residual))) optimum = minimize_scalar( predicted_norm, bounds=(0.0, 1.0), method="bounded", options={"xatol": 1.0e-12}, ) alpha = float(optimum.x) balanced = pressure_projected.copy() balanced[chi2] -= alpha * residual0[chi2] verified = system.evaluate(balanced) predicted = (1.0 - alpha) * residual0 + alpha * residual1 affine_error = float(np.max(np.abs(verified - predicted))) if affine_error > 1.0e-12: raise RuntimeError("chi2 projection segment failed affine verification") candidates = [ ("pressure_only", pressure_projected, residual0), ("pressure_chi2_full", fully_projected, residual1), ("pressure_chi2_balanced", balanced, verified), ] label, chosen, residual = min( candidates, key=lambda item: float(np.max(np.abs(item[2]))) ) return chosen.copy(), { "method": "full_residual_verified_minimax_partial_chi2_projection", "alpha_bounds": [0.0, 1.0], "selected_alpha": alpha, "selected": label, "selected_norm": float(np.max(np.abs(residual))), "pressure_only_norm": float(np.max(np.abs(residual0))), "full_chi2_projection_norm": float(np.max(np.abs(residual1))), "balanced_projection_norm": float(np.max(np.abs(verified))), "affine_verification_error": affine_error, } def _register_projected_candidates(system: ImplicitDAEResidual, label: str, vector: np.ndarray, candidates: list[LiftCandidate]) -> None: try: residual = system.evaluate(vector) except (ValueError, RuntimeError, FloatingPointError): candidates.append(LiftCandidate(label, vector, float("inf"))) return n = system.cfg.candidate_nodes variants = [("raw", vector.copy())] pressure = vector.copy() pressure[5*n:6*n] -= residual[5*n:6*n] variants.append(("pressure", pressure)) pressure_chi2 = pressure.copy() pressure_chi2[6*n:7*n] -= residual[6*n:7*n] variants.append(("pressure_chi2", pressure_chi2)) for suffix, trial in variants: try: norm = float(np.max(np.abs(system.evaluate(trial)))) except (ValueError, RuntimeError, FloatingPointError): norm = float("inf") candidates.append(LiftCandidate(f"{label}_{suffix}", trial, norm)) def select_spectral_lift(target: ImplicitDAEResidual, parents: list[tuple[int, ImplicitDAEResidual, np.ndarray]]) -> tuple[np.ndarray, dict]: candidates: list[LiftCandidate] = [] for nodes, source, vector in parents: _register_projected_candidates( target, f"from_{nodes}_barycentric", _lift_arrays(source, vector, target, "barycentric"), candidates, ) max_degree = min(16, nodes - 1) # The complete bounded low-degree bank is cheap compared with a # nonlinear 65-node solve and avoids assuming that only even or # hand-picked degrees are numerically favorable. Every lift is still # checked against the full target residual before selection. degrees = range(1, max_degree + 1) for degree in degrees: _register_projected_candidates( target, f"from_{nodes}_chebyshev_degree_{degree}", _lift_arrays(source, vector, target, "chebyshev", degree), candidates, ) chosen = min(candidates, key=lambda item: item.norm) return chosen.vector, { "method": "bounded_full_residual_checked_spectral_lift_bank", "selected": chosen.label, "selected_norm": chosen.norm, "candidate_norms": {item.label: item.norm for item in candidates}, } def classify(levels: list[dict]) -> dict: complete = len(levels) == 3 all_gates = complete and all( all(level["diagnostics"]["gate_pass"].values()) for level in levels ) monotonic = {} for key in G22_THRESHOLDS: values = [level["diagnostics"][key] for level in levels] monotonic[key] = complete and all( values[index + 1] <= values[index] for index in range(len(values) - 1) ) passed = bool(all_gates and all(monotonic.values())) return { "pass": passed, "classification": "G22_PASS" if passed else "G22_ORANGE_NOT_PASS", "monotonic_refinement": monotonic, "all_level_gates_pass": bool(all_gates), } def run_suite(state_path: Path, homotopy_report: Path) -> dict: vector17 = load_lambda1_checkpoint(homotopy_report) system17 = ImplicitDAEResidual(state_path, nodes=17, uv_homotopy=1.0) levels = [{ "candidate_nodes": 17, "source": "validated_lambda1_checkpoint", "diagnostics": _physical_diagnostics(system17, vector17), "solution_vector": vector17.tolist(), }] parents = [(17, system17, vector17)] for nodes in (33, 65): system = ImplicitDAEResidual(state_path, nodes=nodes, uv_homotopy=1.0) vector, lift = select_spectral_lift(system, parents) diagnostics = _physical_diagnostics(system, vector) level = { "candidate_nodes": nodes, "source": "spectral_lift_from_validated_lower_resolution_only", "lift": lift, "diagnostics": diagnostics, "accepted_as_checkpoint": bool( diagnostics["gate_pass"]["dae_residual"] and diagnostics["gate_pass"]["structural"] ), } if level["accepted_as_checkpoint"]: level["solution_vector"] = vector.tolist() parents.append((nodes, system, vector)) else: level["nonpass_vector_stored_or_used_as_seed"] = False levels.append(level) return { "schema": "apeiron-ap1-r2c-implicit-dae-multiresolution-v1.0", "classification": classify(levels)["classification"], "thresholds_frozen": G22_THRESHOLDS, "dae_residual_gate": DAE_RESIDUAL_GATE, "levels": levels, "verdict": classify(levels), "old_solver_or_physical_map_called": False, "equations_physics_or_gates_changed": False, "claim_boundary": "17/33/65 numerical and physical gate audit only; no production observable or empirical claim released", } def main() -> None: import argparse parser = argparse.ArgumentParser() parser.add_argument("state", type=Path) parser.add_argument("homotopy_report", type=Path) parser.add_argument("--output", type=Path) args = parser.parse_args() rendered = json.dumps( run_suite(args.state, args.homotopy_report), indent=2 ) + "\n" if args.output: args.output.write_text(rendered, encoding="utf-8") else: print(rendered, end="") if __name__ == "__main__": main()