diff --git a/experiments/analyze_nonlinear_direction_pair.py b/experiments/analyze_nonlinear_direction_pair.py new file mode 100644 index 0000000..7ebf9d0 --- /dev/null +++ b/experiments/analyze_nonlinear_direction_pair.py @@ -0,0 +1,222 @@ +#!/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()