"""Strict PASS-only single-node p continuation for AP1-R2c G22.""" from __future__ import annotations from datetime import datetime, timezone import json from pathlib import Path import numpy as np from ap1_r2c_implicit_dae_multiresolution import ( DAE_RESIDUAL_GATE, G22_THRESHOLDS, _lift_arrays, _physical_diagnostics, balanced_algebraic_projection, endpoint_bspline_lift, richardson_bspline_lift, select_p_lift_projection_bank, ) from ap1_r2c_implicit_dae_homotopy import solve_stage from ap1_r2c_implicit_dae_solver import ImplicitDAEResidual def _all_gates_pass(diagnostics: dict) -> bool: return bool(all(diagnostics["gate_pass"].values())) def repeated_diagnostics(system: ImplicitDAEResidual, vector: np.ndarray, repeats: int = 3) -> dict: """Require repeated full residual checks and report conservative extrema.""" if repeats < 2: raise ValueError("at least two independent full checks required") samples = [_physical_diagnostics(system, vector) for _ in range(repeats)] upper_keys = ( "dae_max_abs_scaled_residual", "max_abs_friedmann_residual", "ward_normalized", "validation_source_relative_change", "validation_wronskian_relative_error", ) conservative = {key: max(sample[key] for sample in samples) for key in upper_keys} conservative.update({ "finite": all(sample["finite"] for sample in samples), "positive_H": all(sample["positive_H"] for sample in samples), "min_H": min(sample["min_H"] for sample in samples), "min_PX": min(sample["min_PX"] for sample in samples), "min_K": min(sample["min_K"] for sample in samples), }) conservative["gate_pass"] = { "dae_residual": conservative["dae_max_abs_scaled_residual"] <= DAE_RESIDUAL_GATE, **{key: conservative[key] <= limit for key, limit in G22_THRESHOLDS.items()}, "structural": bool( conservative["finite"] and conservative["positive_H"] and conservative["min_PX"] > 0.0 and conservative["min_K"] > 0.0 ), } return {"repeats": repeats, "samples": samples, "conservative": conservative} def bounded_residual_correction(system: ImplicitDAEResidual, initial: np.ndarray) -> tuple[np.ndarray, dict]: """Try one fixed safeguarded Newton-Krylov budget without storing NONPASS.""" initial = np.asarray(initial, dtype=float) initial_norm = float(np.max(np.abs(system.evaluate(initial)))) if initial_norm <= DAE_RESIDUAL_GATE: return initial.copy(), { "attempted": False, "reason": "initial_candidate_already_below_DAE_gate", "initial_norm": initial_norm, } corrected, stage = solve_stage( system, initial, max_iterations=4, rdiff=1.0e-7 ) # solve_stage includes its candidate for homotopy checkpointing. Here a # NONPASS must never be serialized, so retain diagnostic scalars only. stage = dict(stage) stage.pop("solution_vector", None) corrected_norm = float(np.max(np.abs(system.evaluate(corrected)))) improved = corrected_norm < initial_norm return (corrected.copy() if improved else initial.copy()), { "attempted": True, "max_iterations": 4, "rdiff": 1.0e-7, "initial_norm": initial_norm, "corrected_norm": corrected_norm, "selected_corrected_candidate": improved, "solver": stage, } def bounded_rho_projection_bank( system: ImplicitDAEResidual, initial: np.ndarray, alphas: tuple[float, ...] = ( 1.0e-5, 3.0e-5, 5.0e-5, 5.25e-5, 5.5e-5, 5.75e-5, 6.0e-5, 6.25e-5, ), ) -> tuple[np.ndarray, dict]: """Select a tiny rho-source projection only after full residual checks.""" if any(not np.isfinite(alpha) or abs(alpha) > 3.0e-4 for alpha in alphas): raise ValueError("bounded finite rho projection factors required") initial = np.asarray(initial, dtype=float) n = system.cfg.candidate_nodes rho = slice(4 * n, 5 * n) residual0 = system.evaluate(initial) candidates = [( float(np.max(np.abs(residual0))), "rho_alpha_0", initial.copy(), None, )] for alpha in alphas: raw = initial.copy() raw[rho] -= float(alpha) * residual0[rho] projected, projection = balanced_algebraic_projection(system, raw) residual = system.evaluate(projected) candidates.append(( float(np.max(np.abs(residual))), f"rho_alpha_{alpha:.6g}", projected, projection, )) norm, label, candidate, projection = min( candidates, key=lambda item: item[0] ) return candidate.copy(), { "method": "bounded_full_residual_checked_partial_rho_projection", "alpha_bounds": [-3.0e-4, 3.0e-4], "selected": label, "selected_norm": norm, "candidate_norms": {item[1]: item[0] for item in candidates}, "selected_algebraic_projection": projection, } def _frozen_rhs_background_direction( system: ImplicitDAEResidual, vector: np.ndarray, ) -> np.ndarray: """Solve the four linear Lobatto defect equations at the current RHS.""" n = system.cfg.candidate_nodes residual = system.evaluate(vector) matrix = system.cfg.delta_N * system.D_N.copy() matrix[0, :] = 0.0 matrix[0, 0] = 1.0 direction = np.zeros_like(vector) for block in range(4): defect = residual[block*n:(block+1)*n].copy() direction[block*n:(block+1)*n] = np.linalg.solve(matrix, -defect) return direction def bounded_frozen_rhs_cross_projection( source: ImplicitDAEResidual, parent: np.ndarray, target: ImplicitDAEResidual, reference: np.ndarray, betas: tuple[float, ...] = (2.2,), background_alphas: tuple[float, ...] = (-1.0e-5,), rho_alphas: tuple[float, ...] = (4.25e-5,), ) -> tuple[np.ndarray, dict]: """Test a fixed bounded lift/background/rho cross projection.""" if any(not np.isfinite(beta) or not -3.0 <= beta <= 5.0 for beta in betas): raise ValueError("bounded finite Richardson factors required") if any(not np.isfinite(alpha) or abs(alpha) > 3.0e-4 for alpha in (*background_alphas, *rho_alphas)): raise ValueError("bounded finite cross-projection factors required") n = target.cfg.candidate_nodes rho = slice(4 * n, 5 * n) reference = np.asarray(reference, dtype=float) reference_norm = float(np.max(np.abs(target.evaluate(reference)))) candidates = [(reference_norm, "reference", reference.copy(), None)] norms = {"reference": reference_norm} for beta in betas: lifted = richardson_bspline_lift(source, parent, target, beta) base, _ = balanced_algebraic_projection(target, lifted) direction = _frozen_rhs_background_direction(target, base) for background_alpha in background_alphas: background = base + float(background_alpha) * direction residual = target.evaluate(background) for rho_alpha in rho_alphas: raw = background.copy() raw[rho] -= float(rho_alpha) * residual[rho] projected, projection = balanced_algebraic_projection( target, raw ) final = target.evaluate(projected) label = ( f"beta_{beta:.6g}_background_{background_alpha:.6g}" f"_rho_{rho_alpha:.6g}" ) norm = float(np.max(np.abs(final))) norms[label] = norm candidates.append((norm, label, projected, projection)) norm, label, candidate, projection = min( candidates, key=lambda item: item[0] ) return candidate.copy(), { "method": "bounded_frozen_RHS_Lobatto_background_rho_cross_projection", "beta_bounds": [-3.0, 5.0], "alpha_bounds": [-3.0e-4, 3.0e-4], "selected": label, "selected_norm": norm, "candidate_norms": norms, "selected_algebraic_projection": projection, } def load_parent33(report_path: Path) -> np.ndarray: report = json.loads(report_path.read_text(encoding="utf-8")) levels = [level for level in report.get("levels", []) if level.get("candidate_nodes") == 33 and level.get("accepted_as_checkpoint") and len(level.get("solution_vector", [])) == 232] if not levels: raise ValueError("complete validated 33-node parent required") return np.asarray(levels[-1]["solution_vector"], dtype=float) def reconstruct_parent34(state_path: Path, report_path: Path, max_attempts: int = 12) -> tuple[np.ndarray, dict]: """Rebuild the documented 34 PASS and retain only a repeated full PASS.""" vector33 = load_parent33(report_path) system33 = ImplicitDAEResidual(state_path, nodes=33, uv_homotopy=1.0) system34 = ImplicitDAEResidual(state_path, nodes=34, uv_homotopy=1.0) attempts = [] for attempt in range(1, max_attempts + 1): vector = _lift_arrays(system33, vector33, system34, "barycentric") residual = system34.evaluate(vector) n = 34 vector[5*n:6*n] -= residual[5*n:6*n] vector[6*n:7*n] -= residual[6*n:7*n] first = _physical_diagnostics(system34, vector) entry = { "attempt": attempt, "first_dae_norm": first["dae_max_abs_scaled_residual"], "first_all_gates_pass": _all_gates_pass(first), } # Recheck only candidates with useful margin. This avoids promoting a # numerically marginal one-shot value to the persistent parent chain. if _all_gates_pass(first) and first["dae_max_abs_scaled_residual"] <= 8.8e-5: repeated = repeated_diagnostics(system34, vector) entry["repeated"] = repeated attempts.append(entry) if _all_gates_pass(repeated["conservative"]): return vector, { "method": "validated_33_barycentric_pressure_chi2_reconstruction", "attempts": attempts, "diagnostics": repeated, } else: attempts.append(entry) raise RuntimeError("bounded 34-node repeated-PASS reconstruction failed") def _atomic_write(path: Path, report: dict) -> None: temporary = path.with_suffix(path.suffix + ".tmp") temporary.write_text(json.dumps(report, indent=2) + "\n", encoding="utf-8") temporary.replace(path) def _new_report() -> dict: return { "schema": "apeiron-ap1-r2c-pass-only-p-continuation-v1.0", "updated_utc": datetime.now(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ"), "classification": "G22_ORANGE_P_CONTINUATION_IN_PROGRESS", "frozen_thresholds": { "dae_max_abs_scaled_residual": DAE_RESIDUAL_GATE, **G22_THRESHOLDS, }, "method": { "lift": "endpoint_preserving_bspline_and_richardson_single_p_bank", "algebraic_projection": "full_residual_verified_minimax_partial_chi2", "equations_changed": False, "physics_changed": False, "gates_changed": False, }, "levels": [], "production_curves_released": False, "nonpass_stored_or_used_as_seed": False, } def run(state_path: Path, multiresolution_report: Path, output: Path, target_nodes: int = 65) -> dict: if not 35 <= target_nodes <= 65: raise ValueError("target must remain inside frozen 35..65 audit range") if output.exists(): report = json.loads(output.read_text(encoding="utf-8")) if report.get("schema") != "apeiron-ap1-r2c-pass-only-p-continuation-v1.0": raise ValueError("recognized p-continuation checkpoint required") accepted = [level for level in report.get("levels", []) if level.get("accepted_as_checkpoint") and level.get("solution_vector")] if not accepted: raise ValueError("resume checkpoint lacks an accepted vector") parent_level = max(accepted, key=lambda level: level["candidate_nodes"]) parent_nodes = int(parent_level["candidate_nodes"]) parent_vector = np.asarray(parent_level["solution_vector"], dtype=float) # Prior failed diagnostics remain useful in their version history, but # the active checkpoint must contain only the new result at the next # level. A NONPASS is never retained as a continuation parent. report["superseded_nonpass_diagnostics"] = [ level for level in report.get("levels", []) if not level.get("accepted_as_checkpoint") and int(level.get("candidate_nodes", -1)) > parent_nodes ] report["levels"] = [ level for level in report.get("levels", []) if level.get("accepted_as_checkpoint") ] else: report = _new_report() parent_vector, reconstruction = reconstruct_parent34( state_path, multiresolution_report ) parent_nodes = 34 report["levels"].append({ "candidate_nodes": 34, "source": "reconstructed_documented_parent_PASS", "accepted_as_checkpoint": True, "reconstruction": reconstruction, "solution_vector": parent_vector.tolist(), }) _atomic_write(output, report) for nodes in range(parent_nodes + 1, target_nodes + 1): source = ImplicitDAEResidual( state_path, nodes=parent_nodes, uv_homotopy=1.0 ) target = ImplicitDAEResidual(state_path, nodes=nodes, uv_homotopy=1.0) candidate, lift_bank = select_p_lift_projection_bank( source, parent_vector, target ) candidate, rho_projection = bounded_rho_projection_bank( target, candidate ) lift_bank["bounded_rho_projection"] = rho_projection candidate, cross_projection = bounded_frozen_rhs_cross_projection( source, parent_vector, target, candidate ) lift_bank["bounded_frozen_rhs_cross_projection"] = cross_projection candidate, correction = bounded_residual_correction(target, candidate) lift_bank["bounded_residual_correction"] = correction diagnostics = repeated_diagnostics(target, candidate) passed = _all_gates_pass(diagnostics["conservative"]) level = { "candidate_nodes": nodes, "source": f"validated_{parent_nodes}_PASS_only", "lift_bank": lift_bank, "diagnostics": diagnostics, "accepted_as_checkpoint": passed, } if passed: level["solution_vector"] = candidate.tolist() parent_nodes = nodes parent_vector = candidate else: level["nonpass_vector_stored_or_used_as_seed"] = False report["levels"].append(level) report["updated_utc"] = datetime.now(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ") report["highest_validated_p_level"] = parent_nodes report["classification"] = ( "G22_P_CONTINUATION_REACHED_65_PENDING_FINAL_AUDIT" if parent_nodes == 65 else f"G22_ORANGE_P_CONTINUATION_STOPPED_AT_{parent_nodes}_PASS" ) _atomic_write(output, report) if not passed: break return report def main() -> None: import argparse parser = argparse.ArgumentParser() parser.add_argument("state", type=Path) parser.add_argument("multiresolution_report", type=Path) parser.add_argument("--output", type=Path, required=True) parser.add_argument("--target-nodes", type=int, default=65) args = parser.parse_args() report = run( args.state, args.multiresolution_report, args.output, args.target_nodes ) print(json.dumps({ "classification": report["classification"], "highest_validated_p_level": report.get("highest_validated_p_level", 34), "levels_recorded": len(report["levels"]), }, indent=2)) if __name__ == "__main__": main()