#!/usr/bin/env python3 """Read-only, two-direction inexact-Newton selection; no new FEM evaluation. Usage: python analyze_nonlinear_direction_pair.py LOOSE_CASE TIGHT_CASE NEW_OUTPUT_DIRECTORY The reports remain preliminary until the separate paired campaign audit passes. """ import argparse import csv import hashlib import json import math from pathlib import Path import platform import sys import numpy as np def sha(path): return hashlib.sha256(path.read_bytes()).hexdigest() def vector(path): values = np.loadtxt(path, dtype=float) assert values.ndim == 1 and int(values[0]) == values.size - 1 values = values[1:] assert np.all(np.isfinite(values)) return values def norm(value): # A distinct summation path checks the small dense model against all saved # coefficients, without asking the FEM operator for any new action. wide = np.asarray(value, dtype=np.longdouble) return float(np.sqrt(np.sum(wide * wide, dtype=np.longdouble))) def write_vector(path, value): with path.open("w") as stream: stream.write(f"{value.size}\n") np.savetxt(stream, value, fmt="%.17g") def combine(columns, weights): # Explicit two-column arithmetic also provides an independent reconstruction # path from the tall QR factors; no large dense product is needed here. return columns[:, 0] * weights[0] + columns[:, 1] * weights[1] def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("loose", type=Path) parser.add_argument("tight", type=Path) parser.add_argument("output", type=Path) args = parser.parse_args() np.seterr(divide="raise", over="raise", invalid="raise") output = args.output.resolve() for campaign in (args.loose.resolve().parent, args.tight.resolve().parent): if output == campaign or campaign in output.parents: parser.error("Keep analysis outside the immutable campaigns") if output.exists(): parser.error("Choose a new output directory") prefix = "nonlinear_step_" shared = ["physical_base.txt", "normalized_base.txt", "residual.txt", "rhs.txt", "state_factors.txt", "residual_factors.txt", "layout.csv"] shared_hashes = {} for suffix in shared: a, b = args.loose / (prefix + suffix), args.tight / (prefix + suffix) assert sha(a) == sha(b), f"Same-state comparison failed: {suffix}" shared_hashes[suffix] = sha(a) inputs = {} for label, directory in [("loose", args.loose), ("tight", args.tight)]: names = [prefix + name for name in shared + ["direction_full.txt", "action_full.txt"]] names += ["metadata.txt"] inputs[label] = {name: {"path": str((directory / name).resolve()), "sha256": sha(directory / name)} for name in names} p = np.column_stack([vector(d / (prefix + "direction_full.txt")) for d in [args.loose, args.tight]]) a = np.column_stack([vector(d / (prefix + "action_full.txt")) for d in [args.loose, args.tight]]) f = vector(args.loose / (prefix + "residual.txt")) f_norm = norm(f) assert p.shape == a.shape == (f.size, 2) and f_norm > 0 eta = 0.03 # A negligible inward margin prevents formatting/roundoff at the active # inequality from turning an algebraically feasible saved vector into a fail. safe_eta = eta * (1.0 - 1.0e-10) q, r = np.linalg.qr(p, mode="reduced") for i in range(2): if r[i, i] < 0: q[:, i] *= -1 r[i, :] *= -1 r_condition = float(np.linalg.cond(r)) assert np.isfinite(r_condition) and r_condition < 1.0e10 q_action = np.empty_like(a) q_action[:, 0] = a[:, 0] / r[0, 0] q_action[:, 1] = (a[:, 1] - r[0, 1] * q_action[:, 0]) / r[1, 1] b = q_action / f_norm f_unit = f / f_norm # Work with QR/SVD of the action itself, never normal equations: the two # input directions can be almost parallel and squared conditioning obscures # whether a small norm reduction is actually resolved. u, t = np.linalg.qr(b, mode="reduced") left, sigma, vt = np.linalg.svd(t, full_matrices=False) action_basis = np.column_stack([combine(u, left[:, i]) for i in range(2)]) beta = action_basis.T @ f_unit perpendicular = f_unit - combine(action_basis, beta) perpendicular_norm = norm(perpendicular) assert sigma[-1] > 0 and perpendicular_norm < safe_eta sigma_scaled = sigma / sigma[0] squared = sigma_scaled * sigma_scaled def solve(target): def residual(mu): retained = beta * mu / (squared + mu) return math.hypot(perpendicular_norm, norm(retained)) lo, hi = 0.0, 1.0 while residual(hi) < target: hi *= 2 for _ in range(120): mid = (lo + hi) / 2 if residual(mid) <= target: lo = mid else: hi = mid mu = lo y = -(sigma_scaled / (squared + mu)) * (beta / sigma[0]) z = vt.T @ y coefficients = np.linalg.solve(r, z) return z, coefficients, mu, residual(mu) z, coefficients, mu, projected_defect = solve(safe_eta) exact_z, _, _, _ = solve(eta) candidate = combine(p, coefficients) candidate_action = combine(a, coefficients) candidate_defect = f + candidate_action candidate_norm = norm(candidate) verified_defect = norm(candidate_defect) / f_norm assert verified_defect <= eta reconstruction_error = norm(candidate - combine(q, z)) action_reconstruction_error = norm(candidate_action - combine(q_action, z)) # Stationarity and the active boundary provide a convex KKT certificate # within this two-dimensional subspace; they say nothing about global rank. mu_physical = mu * sigma[0] ** 2 reduced_residual = t @ z + u.T @ f_unit kkt = t.T @ reduced_residual + mu_physical * z kkt_scale = norm(t.T @ reduced_residual) + mu_physical * norm(z) originals = [] for i, label in enumerate(["loose", "tight"]): ai = a[:, i] / f_norm norm_ai = norm(ai) center = -float(np.dot(f_unit, ai)) / norm_ai ** 2 ray_minimum_defect = norm(f_unit + center * ai) assert ray_minimum_defect < safe_eta radius = math.sqrt(safe_eta ** 2 - ray_minimum_defect ** 2) / norm_ai ray_scale = max(0.0, center - radius) original_norm = norm(p[:, i]) original_defect = norm(f + a[:, i]) / f_norm assert original_defect <= (0.03 if i == 0 else 0.003) ray_norm = ray_scale * original_norm originals.append({"label": label, "norm": original_norm, "relative_saved_action_defect": original_defect, "candidate_norm_reduction_fraction": 1 - candidate_norm / original_norm, "best_feasible_positive_ray_scale": ray_scale, "best_feasible_positive_ray_norm": ray_norm, "extra_fractional_norm_reduction_over_ray": 1 - candidate_norm / ray_norm}) assert candidate_norm <= min(x["best_feasible_positive_ray_norm"] for x in originals) * (1 + 1e-10) assert norm(kkt) <= 1e-9 * kkt_scale assert abs(projected_defect - verified_defect) < 1e-10 assert reconstruction_error < 1e-12 * candidate_norm pair_difference = p[:, 1] - p[:, 0] difference_action = a[:, 1] - a[:, 0] output_vectors = {"candidate_normalized.txt": candidate, "candidate_saved_action.txt": candidate_action, "candidate_saved_defect.txt": candidate_defect} output_hashes = {} args.output.mkdir(exist_ok=False, parents=True) for name, value in output_vectors.items(): write_vector(args.output / name, value) output_hashes[name] = sha(args.output / name) assert np.array_equal(vector(args.output / name), value) report = { "status": "preliminary_pending_independent_paired_campaign_audit", "scope": "Algebraic selection in span of two saved same-state normalized directions. No new FEM action or nonlinear trial.", "python": sys.executable, "python_version": platform.python_version(), "numpy_version": np.__version__, "script_sha256": sha(Path(__file__)), "inputs": inputs, "shared_hashes": shared_hashes, "size": int(f.size), "base_residual_norm": f_norm, "eta": eta, "safe_eta": safe_eta, "originals": originals, "pair": {"cosine": float(np.dot(p[:, 0], p[:, 1])) / (norm(p[:, 0]) * norm(p[:, 1])), "difference_norm": norm(pair_difference), "difference_relative_to_loose": norm(pair_difference) / norm(p[:, 0]), "difference_action_norm": norm(difference_action), "difference_directional_gain": norm(difference_action) / norm(pair_difference)}, "reduced_problem": {"direction_R_condition_2": r_condition, "action_T_singular_values_relative_to_F": sigma.tolist(), "action_T_condition_2": float(sigma[0] / sigma[-1]), "irreducible_relative_defect_in_span": perpendicular_norm, "orthogonality_error": norm((q.T @ q - np.eye(2)).ravel()), "scaled_ridge_parameter": mu, "KKT_relative_residual": norm(kkt) / kkt_scale, "safety_margin_norm_increment": candidate_norm - norm(exact_z)}, "candidate": {"original_direction_coefficients": coefficients.tolist(), "orthonormal_coordinates": z.tolist(), "norm": candidate_norm, "relative_defect_projected": projected_defect, "relative_defect_full_saved_vectors": verified_defect, "direction_reconstruction_error": reconstruction_error, "action_reconstruction_error": action_reconstruction_error, "output_sha256": output_hashes}, } (args.output / "assessment.json").write_text(json.dumps(report, indent=2, allow_nan=False) + "\n") print(json.dumps({k: report[k] for k in ["status", "originals", "pair", "reduced_problem", "candidate"]}, indent=2)) if __name__ == "__main__": main()