223 lines
10 KiB
Python
223 lines
10 KiB
Python
#!/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()
|