367 lines
20 KiB
Python
367 lines
20 KiB
Python
#!/usr/bin/env python3
|
|
"""Describe saved radial surface directions by unweighted nodal angular fits.
|
|
|
|
These are projections of the sampled coefficient values, not surface-integral
|
|
harmonic coefficients. Cubic symmetry can suggest a mesh-related pattern, but
|
|
does not establish its cause, physical accuracy, or a nullspace of the solver.
|
|
No FEM operations or full-state vector reads are performed.
|
|
"""
|
|
import argparse
|
|
import csv
|
|
import hashlib
|
|
import itertools
|
|
import json
|
|
import math
|
|
import os
|
|
from pathlib import Path
|
|
import sys
|
|
|
|
for _name in ("OPENBLAS_NUM_THREADS", "OMP_NUM_THREADS", "MKL_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"):
|
|
os.environ[_name] = "1"
|
|
import numpy as np
|
|
|
|
SURFACE_FILE = "weak_shape_physical_surface_nodes.csv"
|
|
DEGREES = (0, 2, 4, 6, 8)
|
|
EPS = np.finfo(float).eps
|
|
NODE_TOLERANCE = 2e-10
|
|
|
|
|
|
def digest(path):
|
|
return hashlib.sha256(Path(path).read_bytes()).hexdigest()
|
|
|
|
|
|
def contract(subscripts, *arrays):
|
|
# Direct contraction avoids the Accelerate matmul flag defect documented
|
|
# in audit_shape_field_weak.py; finite results remain an explicit check.
|
|
if any(not np.isfinite(a).all() for a in arrays):
|
|
raise ValueError("Nonfinite contraction input")
|
|
result = np.einsum(subscripts, *arrays, optimize=False)
|
|
if not np.isfinite(result).all():
|
|
raise ValueError("Nonfinite contraction output")
|
|
return result
|
|
|
|
|
|
def norm(array):
|
|
result = math.sqrt(float(contract("i,i->", array, array)))
|
|
if not math.isfinite(result):
|
|
raise ValueError("Nonfinite norm")
|
|
return result
|
|
|
|
|
|
def harmonics(points, maximum_degree=8):
|
|
"""Continuous-L2-normalized real Y_lm, including Condon--Shortley phase.
|
|
|
|
Recurring on normalized associated Legendre functions avoids separately
|
|
generating large factorials and then cancelling them in normalization.
|
|
Negative m denotes sine, positive m cosine, with sqrt(2) for nonzero m.
|
|
"""
|
|
if maximum_degree < 0 or maximum_degree > 8:
|
|
raise ValueError("This bounded diagnostic supports degrees zero through eight")
|
|
z, phi = points[:, 2], np.arctan2(points[:, 1], points[:, 0])
|
|
radial = np.hypot(points[:, 0], points[:, 1])
|
|
legendre = {(0, 0): np.full(len(points), 1 / math.sqrt(4 * math.pi))}
|
|
for m in range(maximum_degree + 1):
|
|
if m:
|
|
legendre[m, m] = -math.sqrt((2 * m + 1) / (2 * m)) * radial * legendre[m - 1, m - 1]
|
|
if m < maximum_degree:
|
|
legendre[m + 1, m] = math.sqrt(2 * m + 3) * z * legendre[m, m]
|
|
for degree in range(m + 2, maximum_degree + 1):
|
|
a = math.sqrt((4 * degree * degree - 1) / (degree * degree - m * m))
|
|
b = math.sqrt(((2 * degree + 1) * ((degree - 1) ** 2 - m * m)) /
|
|
((2 * degree - 3) * (degree * degree - m * m)))
|
|
legendre[degree, m] = a * z * legendre[degree - 1, m] - b * legendre[degree - 2, m]
|
|
labels, columns = [], []
|
|
for degree in range(maximum_degree + 1):
|
|
for m in range(-degree, degree + 1):
|
|
value = legendre[degree, abs(m)]
|
|
if m:
|
|
value = math.sqrt(2) * value * (np.sin(-m * phi) if m < 0 else np.cos(m * phi))
|
|
labels.append({"degree": degree, "m": m})
|
|
columns.append(value)
|
|
matrix = np.column_stack(columns)
|
|
if not np.isfinite(matrix).all():
|
|
raise ValueError("Nonfinite harmonic basis")
|
|
return matrix, labels
|
|
|
|
|
|
def fit(matrix, values, column_scale=False):
|
|
"""Tall-matrix SVD with explicit reconstruction/orthogonality checks."""
|
|
if not np.isfinite(matrix).all() or not np.isfinite(values).all():
|
|
raise ValueError("Nonfinite fit input")
|
|
scales = np.sqrt(contract("ij,ij->j", matrix, matrix)) if column_scale else np.ones(matrix.shape[1])
|
|
scales = np.where(scales > 0, scales, 1.)
|
|
scaled = matrix / scales
|
|
u, singular, vh = np.linalg.svd(scaled, full_matrices=False)
|
|
if not all(np.isfinite(a).all() for a in (u, singular, vh)):
|
|
raise ValueError("Nonfinite SVD")
|
|
cutoff = EPS * max(matrix.shape) * singular[0]
|
|
rank = int(np.count_nonzero(singular > cutoff))
|
|
scores = contract("ij,i->j", u[:, :rank], values)
|
|
coefficients = contract("ji,j->i", vh[:rank], scores / singular[:rank]) / scales
|
|
projected = contract("ij,j->i", matrix, coefficients)
|
|
remainder = values - projected
|
|
reconstruction = contract("ij,jk->ik", u * singular, vh)
|
|
factorization_error = norm((reconstruction - scaled).ravel())
|
|
factorization_limit = 4096 * EPS * max(norm(scaled.ravel()), 1.)
|
|
if factorization_error > factorization_limit:
|
|
raise ValueError("SVD reconstruction failed")
|
|
# The retained singular directions define the actual numerical projection;
|
|
# omitted singular directions need not be resolved by the node set.
|
|
orthogonality = norm(contract("ij,i->j", u[:, :rank], remainder))
|
|
magnitude = norm(values)
|
|
limit = 8192 * EPS * max(magnitude, norm(projected), 1e-300)
|
|
if orthogonality > limit:
|
|
raise ValueError("Projection orthogonality failed")
|
|
error = norm(remainder)
|
|
return {"sample_count": len(values), "column_count": matrix.shape[1], "rank": rank,
|
|
"full_column_rank": rank == matrix.shape[1], "rank_cutoff": float(cutoff),
|
|
"singular_values": singular.tolist(), "column_scaling": "unit nodal L2" if column_scale else "none",
|
|
"column_scales": scales.tolist(),
|
|
"retained_condition_ratio": float(singular[0] / singular[rank - 1]) if rank else None,
|
|
"basis_coefficients": coefficients.tolist(), "nodal_l2_remainder": error,
|
|
"relative_nodal_l2_remainder": error / magnitude if magnitude else None,
|
|
"explained_nodal_squared_norm_fraction": 1 - (error / magnitude) ** 2 if magnitude else None,
|
|
"svd_reconstruction_error": factorization_error, "retained_projection_orthogonality_error": orthogonality}
|
|
|
|
|
|
def cubic_columns(points):
|
|
squared = points * points
|
|
c4 = np.sum(squared * squared, axis=1) - 3 / 5
|
|
# On the unit sphere, the reflection/permutation-invariant polynomial ring
|
|
# through degree eight is spanned by 1, e2, e3, e2^2. C4 is an affine
|
|
# replacement for e2; e3's spherical mean is 1/105. These columns are not
|
|
# individually harmonic or orthogonal, except the pure l=4 C4 column.
|
|
e3 = np.prod(squared, axis=1) - 1 / 105
|
|
return np.column_stack((np.ones(len(points)), c4, e3, c4 * c4))
|
|
|
|
|
|
def symmetry_projection(points, values):
|
|
"""Project onto all signed-permutation invariants if nodes form orbits.
|
|
|
|
A permutation of the entire node set preserves the unweighted metric.
|
|
Without that closure, a nodal group average is not a justified orthogonal
|
|
projection; invariant polynomial fits remain separately available.
|
|
"""
|
|
buckets = {}
|
|
for i, point in enumerate(points):
|
|
key = tuple(np.rint(point / NODE_TOLERANCE).astype(np.int64))
|
|
buckets.setdefault(key, []).append(i)
|
|
shifts = tuple(itertools.product((-1, 0, 1), repeat=3))
|
|
transforms = []
|
|
maximum_error = 0.
|
|
for permutation in itertools.permutations(range(3)):
|
|
for signs in itertools.product((-1, 1), repeat=3):
|
|
image_indices = []
|
|
for point in points[:, permutation] * np.asarray(signs):
|
|
key = tuple(np.rint(point / NODE_TOLERANCE).astype(np.int64))
|
|
candidates = []
|
|
for shift in shifts:
|
|
for index in buckets.get(tuple(key[k] + shift[k] for k in range(3)), ()):
|
|
distance = float(np.max(np.abs(point - points[index])))
|
|
if distance <= NODE_TOLERANCE:
|
|
candidates.append((distance, index))
|
|
if len(candidates) != 1:
|
|
return {"status": "unavailable_node_cloud_not_uniquely_closed", "node_tolerance": NODE_TOLERANCE,
|
|
"matching_node_count_at_failure": len(candidates),
|
|
"scope": "No full cubic projection claimed; polynomial invariant fits remain valid."}
|
|
distance, index = candidates[0]
|
|
maximum_error = max(maximum_error, distance)
|
|
image_indices.append(index)
|
|
if len(set(image_indices)) != len(points):
|
|
raise ValueError("Cubic transformation did not define a nodal permutation")
|
|
transforms.append(image_indices)
|
|
indices = np.asarray(transforms, dtype=int)
|
|
projected = np.mean(values[indices], axis=0)
|
|
error = norm(values - projected)
|
|
magnitude = norm(values)
|
|
idempotency = norm(np.mean(projected[indices], axis=0) - projected)
|
|
cross = abs(float(contract("i,i->", projected, values - projected)))
|
|
if idempotency > 4096 * EPS * max(magnitude, 1e-300) or cross > 4096 * EPS * max(magnitude ** 2, 1e-300):
|
|
raise ValueError("Cubic nodal projection identities failed")
|
|
return {"status": "closed_48_element_group", "node_tolerance": NODE_TOLERANCE,
|
|
"maximum_node_matching_error": maximum_error, "nodal_l2_remainder": error,
|
|
"relative_nodal_l2_remainder": error / magnitude if magnitude else None,
|
|
"invariant_nodal_squared_norm_fraction": 1 - (error / magnitude) ** 2 if magnitude else None,
|
|
"projection_idempotency_error": idempotency, "projection_cross_term": cross,
|
|
"scope": "All signed coordinate permutations in the saved Cartesian axes; unweighted nodal group average."}
|
|
|
|
|
|
def read_surface(path):
|
|
with path.open(newline="") as stream:
|
|
reader = csv.DictReader(stream)
|
|
required = {"parameter", "volume_scalar_true_dof", "unit_radial_x", "unit_radial_y", "unit_radial_z",
|
|
"reference_radius", "physical_radius", "base_parameter", "physical_direction_parameter"}
|
|
if not reader.fieldnames or len(set(reader.fieldnames)) != len(reader.fieldnames) or not required <= set(reader.fieldnames):
|
|
raise ValueError("Missing or duplicate surface CSV columns")
|
|
rows = list(reader)
|
|
if not rows or any(None in r or any(v is None for v in r.values()) for r in rows):
|
|
raise ValueError("Empty or malformed surface CSV")
|
|
data = {k: np.asarray([float(r[k]) for r in rows]) for k in required}
|
|
if any(not np.isfinite(a).all() for a in data.values()):
|
|
raise ValueError("Nonfinite surface data")
|
|
if not np.array_equal(data["parameter"], np.arange(len(rows))):
|
|
raise ValueError("Noncontiguous surface parameter indices")
|
|
identifiers = data["volume_scalar_true_dof"]
|
|
if np.any(identifiers < 0) or np.any(identifiers != np.floor(identifiers)) or len(set(identifiers)) != len(rows):
|
|
raise ValueError("Invalid or repeated surface degree of freedom")
|
|
if np.any(data["reference_radius"] <= 0) or np.any(data["physical_radius"] <= 0):
|
|
raise ValueError("Nonpositive surface radius")
|
|
discrepancy = np.abs(data["physical_radius"] - data["reference_radius"] - data["base_parameter"])
|
|
if np.any(discrepancy > 2048 * EPS * (data["physical_radius"] + data["reference_radius"] + np.abs(data["base_parameter"]))):
|
|
raise ValueError("Physical radius identity failed")
|
|
points = np.column_stack([data["unit_radial_" + axis] for axis in "xyz"])
|
|
lengths = np.sqrt(contract("ij,ij->i", points, points))
|
|
unit_error = float(np.max(np.abs(lengths - 1)))
|
|
if unit_error > 5e-10:
|
|
raise ValueError("Surface radial directions are not unit length")
|
|
# Remove decimal-rounding length error before the spherical recurrence.
|
|
return points / lengths[:, None], data["physical_direction_parameter"], unit_error, data
|
|
|
|
|
|
def analyze(path):
|
|
path = Path(path).resolve()
|
|
if path.is_dir():
|
|
path /= SURFACE_FILE
|
|
points, values, unit_error, surface_data = read_surface(path)
|
|
source_hash = digest(path)
|
|
binding = {"status": "input_hash_only"}
|
|
run_path = path.parent.parent / (path.parent.name + "_run.json")
|
|
if run_path.exists():
|
|
run = json.loads(run_path.read_text())
|
|
if run.get("output_hashes", {}).get(path.name) != source_hash:
|
|
raise ValueError("Surface file does not match runner output hash")
|
|
binding = {"status": "runner_surface_hash_checked", "run_path": str(run_path), "run_sha256": digest(run_path),
|
|
"scope": "Only this surface artifact binding; no independent audit of the complete FEM campaign."}
|
|
matrix, labels = harmonics(points)
|
|
harmonic_fits = [{"maximum_degree": degree, **fit(matrix[:, :(degree + 1) ** 2], values)} for degree in DEGREES]
|
|
cubic = cubic_columns(points)
|
|
invariant_fits = [{"maximum_degree": degree, **fit(cubic[:, :count], values, True)}
|
|
for degree, count in ((0, 1), (2, 1), (4, 2), (6, 3), (8, 4))]
|
|
return {"surface_file": str(path), "surface_sha256": source_hash, "binding": binding,
|
|
"sample_count": len(points), "maximum_original_unit_length_error": unit_error,
|
|
"base_surface_radius": {
|
|
"scope": "unweighted nodal statistics; extrema bound observed asphericity, not full-surface extrema",
|
|
"reference_minimum": float(np.min(surface_data["reference_radius"])),
|
|
"reference_maximum": float(np.max(surface_data["reference_radius"])),
|
|
"physical_minimum": float(np.min(surface_data["physical_radius"])),
|
|
"physical_maximum": float(np.max(surface_data["physical_radius"])),
|
|
"physical_mean": float(np.mean(surface_data["physical_radius"])),
|
|
"physical_span_over_mean": float(np.ptp(surface_data["physical_radius"]) / np.mean(surface_data["physical_radius"])),
|
|
"maximum_absolute_base_parameter": float(np.max(np.abs(surface_data["base_parameter"])))},
|
|
"direction_nodal_l2": norm(values), "direction_maximum_absolute": float(np.max(np.abs(values))),
|
|
"harmonic_basis": {"normalization": "continuous unit-sphere L2 orthonormal real Y_lm with Condon--Shortley phase",
|
|
"ordering": labels, "negative_m": "sine", "positive_m": "cosine"},
|
|
"harmonic_fits": harmonic_fits,
|
|
"cubic_invariant_basis": ["1", "x^4+y^4+z^4-3/5", "x^2*y^2*z^2-1/105", "(x^4+y^4+z^4-3/5)^2"],
|
|
"cubic_invariant_fits": invariant_fits, "full_cubic_nodal_projection": symmetry_projection(points, values)}
|
|
|
|
|
|
def self_test():
|
|
# Product quadrature integrates products through l=8 exactly, independently
|
|
# checking the recurrence's normalization, phase and degree organization.
|
|
z, weights = np.polynomial.legendre.leggauss(9)
|
|
phi = np.arange(18) * (2 * math.pi / 18)
|
|
points = np.asarray([(math.sqrt(1 - q * q) * math.cos(t), math.sqrt(1 - q * q) * math.sin(t), q) for q in z for t in phi])
|
|
quadrature = np.repeat(weights, len(phi)) * (2 * math.pi / len(phi))
|
|
matrix, labels = harmonics(points)
|
|
gram = contract("ki,k,kj->ij", matrix, quadrature, matrix)
|
|
assert np.max(np.abs(gram - np.eye(81))) < 8e-14
|
|
assert np.max(np.abs(matrix[:, 2] - math.sqrt(3 / (4 * math.pi)) * points[:, 2])) < 1e-15
|
|
assert np.max(np.abs(matrix[:, 3] + math.sqrt(3 / (4 * math.pi)) * points[:, 0])) < 1e-15
|
|
for degree in range(9):
|
|
columns = matrix[:, degree * degree:(degree + 1) ** 2]
|
|
assert np.max(np.abs(np.sum(columns * columns, axis=1) - (2 * degree + 1) / (4 * math.pi))) < 2e-14
|
|
cubic = cubic_columns(points)
|
|
assert fit(matrix[:, :9], cubic[:, 1])["relative_nodal_l2_remainder"] > .8
|
|
assert fit(matrix[:, :25], cubic[:, 1])["relative_nodal_l2_remainder"] < 1e-13
|
|
for degree in DEGREES:
|
|
values = matrix[:, degree * degree]
|
|
result = fit(matrix[:, :(degree + 1) ** 2], values)
|
|
assert result["rank"] == (degree + 1) ** 2 and result["relative_nodal_l2_remainder"] < 1e-13
|
|
assert fit(matrix[:, :1], np.zeros(len(points)))["relative_nodal_l2_remainder"] is None
|
|
assert fit(np.column_stack((matrix[:, 0], matrix[:, 0])), matrix[:, 0])["rank"] == 1
|
|
cloud = []
|
|
for seed in ((1, 2, 3), (1, 3, 5), (2, 3, 7), (1, 4, 9)):
|
|
for permutation in itertools.permutations(seed):
|
|
for signs in itertools.product((-1, 1), repeat=3):
|
|
point = np.asarray(permutation, dtype=float) * signs
|
|
cloud.append(point / np.linalg.norm(point))
|
|
cloud = np.asarray(cloud)
|
|
cubic = cubic_columns(cloud)
|
|
invariant = contract("ij,j->i", cubic, np.asarray((.4, 2., -7., .6)))
|
|
assert fit(cubic, invariant, True)["relative_nodal_l2_remainder"] < 1e-13
|
|
assert symmetry_projection(cloud, invariant)["relative_nodal_l2_remainder"] < 1e-13
|
|
antisymmetric = cloud[:, 0]
|
|
result = symmetry_projection(cloud, invariant + antisymmetric)
|
|
assert abs(result["nodal_l2_remainder"] - norm(antisymmetric)) < 1e-12
|
|
assert symmetry_projection(cloud[:-1], invariant[:-1])["status"] == "unavailable_node_cloud_not_uniquely_closed"
|
|
# Bind the analysis to a tiny real CSV and reject corruption in both parser
|
|
# and runner certificate, without reading any campaign artifacts.
|
|
import tempfile
|
|
with tempfile.TemporaryDirectory() as tmp:
|
|
root = Path(tmp)
|
|
folder = root / "case"
|
|
folder.mkdir()
|
|
path = folder / SURFACE_FILE
|
|
columns = ["parameter", "volume_scalar_true_dof", "unit_radial_x", "unit_radial_y", "unit_radial_z",
|
|
"reference_radius", "physical_radius", "base_parameter", "physical_direction_parameter"]
|
|
with path.open("w", newline="") as stream:
|
|
writer = csv.writer(stream)
|
|
writer.writerow(columns)
|
|
for i, (point, value) in enumerate(zip(cloud, invariant)):
|
|
writer.writerow([i, i, *point, 1., 1., 0., value])
|
|
run = root / "case_run.json"
|
|
run.write_text(json.dumps({"output_hashes": {SURFACE_FILE: digest(path)}}))
|
|
assert analyze(folder)["binding"]["status"] == "runner_surface_hash_checked"
|
|
run.write_text(json.dumps({"output_hashes": {SURFACE_FILE: "0" * 64}}))
|
|
try:
|
|
analyze(folder)
|
|
except ValueError as exc:
|
|
assert "runner output hash" in str(exc)
|
|
else:
|
|
raise AssertionError("Corrupted source hash passed")
|
|
path.write_text(path.read_text().replace(",1.0,1.0,0.0,", ",nan,1.0,0.0,", 1))
|
|
try:
|
|
read_surface(path)
|
|
except ValueError as exc:
|
|
assert "Nonfinite" in str(exc)
|
|
else:
|
|
raise AssertionError("Nonfinite source passed")
|
|
print("Tiny harmonic, invariant, symmetry, rank and artifact-binding fixtures passed")
|
|
|
|
|
|
def main():
|
|
parser = argparse.ArgumentParser(description=__doc__)
|
|
parser.add_argument("inputs", nargs="*", help="Probe directories or saved surface-node CSV files")
|
|
parser.add_argument("--output", type=Path)
|
|
parser.add_argument("--self-test", action="store_true")
|
|
args = parser.parse_args()
|
|
if args.self_test:
|
|
if args.inputs or args.output:
|
|
parser.error("Self-test does not consume campaign inputs or create an output")
|
|
self_test()
|
|
return 0
|
|
if not args.inputs or not args.output:
|
|
parser.error("Provide surface inputs and a fresh output JSON")
|
|
output = args.output.resolve()
|
|
if output.exists() or any(output.is_relative_to(Path(p).resolve() if Path(p).is_dir() else Path(p).resolve().parent) for p in args.inputs):
|
|
parser.error("Use a fresh output outside the immutable input directories")
|
|
cases = []
|
|
for path in args.inputs:
|
|
try:
|
|
cases.append({"analysis_pass": True, **analyze(path)})
|
|
except (ValueError, OSError, KeyError, np.linalg.LinAlgError) as exc:
|
|
cases.append({"input": str(Path(path).resolve()), "analysis_pass": False, "error": str(exc)})
|
|
result = {"schema_version": 1, "passed": all(c["analysis_pass"] for c in cases),
|
|
"analyzer_sha256": digest(__file__), "numpy_version": np.__version__,
|
|
"scope": "Unweighted nodal projections of the physical radial direction. Not area-weighted harmonic power, not a continuum symmetry or solver-nullspace certificate.",
|
|
"cases": cases}
|
|
output.parent.mkdir(parents=True, exist_ok=True)
|
|
output.write_text(json.dumps(result, indent=2, allow_nan=False) + "\n")
|
|
print(f"{'PASS' if result['passed'] else 'FAIL'}: {len(cases)} surface analyses; {output}")
|
|
return 0 if result["passed"] else 1
|
|
|
|
|
|
if __name__ == "__main__":
|
|
sys.exit(main())
|