800 lines
52 KiB
Python
800 lines
52 KiB
Python
#!/usr/bin/env python3
|
|
"""Audit frozen shape/field artifacts without reevaluating FEM.
|
|
|
|
Use the NumPy-enabled Python environment. The SVD is taken directly on the tall
|
|
matrix of saved block actions: squaring its condition number through a Gram
|
|
matrix would obscure precisely the weak response this experiment investigates.
|
|
Physical observations use broken scalar gradients, not RT transport or DG jump
|
|
distributions. Neither these observations nor a small restricted singular value
|
|
identify a physical gauge or establish the rank of the complete Jacobian.
|
|
"""
|
|
import argparse
|
|
import csv
|
|
import hashlib
|
|
import json
|
|
import math
|
|
import os
|
|
from pathlib import Path
|
|
import sys
|
|
|
|
# Avoid competing with a later numerical campaign through a threaded BLAS pool.
|
|
for _name in ("OPENBLAS_NUM_THREADS", "OMP_NUM_THREADS", "MKL_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"):
|
|
os.environ[_name] = "1"
|
|
import numpy as np
|
|
|
|
P = "weak_shape_"
|
|
PP = P + "physical_"
|
|
EPS = np.finfo(float).eps
|
|
RTOL = 5e-10
|
|
FIELDS = ("density", "specific_enthalpy", "gravity_potential")
|
|
INTEGRALS = ("base_integral", "pullback_integral", "transport_integral", "eulerian_integral",
|
|
"base_squared_integral", "pullback_squared_integral", "transport_squared_integral",
|
|
"eulerian_squared_integral", "pullback_transport_cross_integral")
|
|
SQUARES = INTEGRALS[4:8]
|
|
|
|
|
|
def digest(path):
|
|
h = hashlib.sha256()
|
|
with Path(path).open("rb") as stream:
|
|
for chunk in iter(lambda: stream.read(1024 * 1024), b""):
|
|
h.update(chunk)
|
|
return h.hexdigest()
|
|
|
|
|
|
def rows(path):
|
|
with Path(path).open(newline="") as stream:
|
|
reader = csv.DictReader(stream)
|
|
if not reader.fieldnames or len(set(reader.fieldnames)) != len(reader.fieldnames):
|
|
raise ValueError(f"Missing or duplicate CSV columns: {path}")
|
|
result = list(reader)
|
|
if any(None in row or any(v is None for v in row.values()) for row in result):
|
|
raise ValueError(f"Malformed CSV row: {path}")
|
|
return result
|
|
|
|
|
|
def finite(value):
|
|
number = float(value)
|
|
if not math.isfinite(number):
|
|
raise ValueError(f"Nonfinite available value: {value}")
|
|
return number
|
|
|
|
|
|
def integer(value):
|
|
number = finite(value)
|
|
if number != int(number):
|
|
raise ValueError(f"Nonintegral count/index: {value}")
|
|
return int(number)
|
|
|
|
|
|
def vector(path):
|
|
with Path(path).open() as stream:
|
|
size = int(stream.readline().strip())
|
|
data = np.loadtxt(stream, dtype=float, ndmin=1)
|
|
if size <= 0 or data.ndim != 1 or data.size != size or not np.isfinite(data).all():
|
|
raise ValueError(f"Invalid finite vector: {path}")
|
|
return data
|
|
|
|
|
|
def norm(value):
|
|
return float(np.linalg.norm(value))
|
|
|
|
|
|
def checked_einsum(subscripts, *operands):
|
|
# NumPy 2.2.6 / Accelerate can flag invalid/divide/overflow for a bounded
|
|
# 257x6 matrix product whose answer agrees with scalar fsum to 4e-15.
|
|
# Direct contractions avoid that BLAS path; finite and reconstruction
|
|
# checks remain mandatory, with no warning suppression or relaxed bounds.
|
|
if any(not np.isfinite(operand).all() for operand in operands):
|
|
raise ValueError("Nonfinite contraction input")
|
|
result = np.einsum(subscripts, *operands, optimize=False)
|
|
if not np.isfinite(result).all():
|
|
raise ValueError("Nonfinite contraction result")
|
|
return result
|
|
|
|
|
|
def near(a, b, scale=0.):
|
|
return math.isfinite(a) and math.isfinite(b) and abs(a - b) <= (
|
|
RTOL * max(abs(a), abs(b)) + 2048 * EPS * abs(scale) + 1e-28)
|
|
|
|
|
|
def vector_near(a, b, scale=0.):
|
|
return a.shape == b.shape and near(norm(a - b), 0., scale + norm(a) + norm(b))
|
|
|
|
|
|
def ratio(a, b):
|
|
return finite(a / b) if b > 0 else None
|
|
|
|
|
|
class Audit:
|
|
def __init__(self):
|
|
self.count = 0
|
|
self.failures = []
|
|
|
|
def check(self, name, condition):
|
|
self.count += 1
|
|
if not bool(condition):
|
|
self.failures.append(name)
|
|
|
|
def close(self, name, actual, expected, scale=0.):
|
|
self.check(name, near(actual, expected, scale))
|
|
|
|
|
|
def table(path, key):
|
|
result = {}
|
|
for row in rows(path):
|
|
index = key(row)
|
|
if index in result:
|
|
raise ValueError(f"Duplicate row {index}: {path}")
|
|
result[index] = row
|
|
return result
|
|
|
|
|
|
def layout(path, size):
|
|
result = {}
|
|
all_rows = rows(path)
|
|
if any(r["kind"] not in ("value", "residual") for r in all_rows):
|
|
raise ValueError("Unknown layout kind")
|
|
for kind in ("value", "residual"):
|
|
entries = [r for r in all_rows if r["kind"] == kind]
|
|
offset, names, blocks = 0, set(), []
|
|
for row in entries:
|
|
start, count = integer(row["offset"]), integer(row["size"])
|
|
if start != offset or count <= 0 or row["block"] in names:
|
|
raise ValueError(f"Invalid {kind} layout partition")
|
|
names.add(row["block"])
|
|
blocks.append((row["block"], start, count))
|
|
offset += count
|
|
if offset != size:
|
|
raise ValueError(f"Incomplete {kind} layout partition")
|
|
result[kind] = blocks
|
|
return result
|
|
|
|
|
|
def restricted_svd(direction, actions, blocks, audit):
|
|
norms = np.array([norm(direction[o:o+n]) for _, o, n in blocks])
|
|
active = [i for i, x in enumerate(norms) if x > 0.]
|
|
for i, x in enumerate(norms):
|
|
if x == 0.:
|
|
audit.check(f"zero_direction_block_action:{i}", norm(actions[i]) == 0.)
|
|
if not active:
|
|
raise ValueError("No nonzero block directions")
|
|
matrix = np.column_stack([actions[i] / norms[i] for i in active])
|
|
if not np.isfinite(matrix).all():
|
|
raise ValueError("Nonfinite scaled block action")
|
|
u, singular, vh = np.linalg.svd(matrix, full_matrices=False)
|
|
weights = vh[-1].copy()
|
|
# Fix sign for reproducible reporting; repeated singular values still permit
|
|
# different equally valid vectors, which the JSON explicitly acknowledges.
|
|
if weights[np.argmax(np.abs(weights))] < 0:
|
|
weights *= -1
|
|
combination = np.zeros_like(direction)
|
|
action = np.zeros_like(direction)
|
|
coefficients = {}
|
|
for weight, i in zip(weights, active):
|
|
name, start, count = blocks[i]
|
|
coefficient = float(weight / norms[i])
|
|
combination[start:start+count] = coefficient * direction[start:start+count]
|
|
action += coefficient * actions[i]
|
|
coefficients[name] = {"orthonormal_basis_weight": float(weight),
|
|
"saved_block_direction_multiplier": coefficient}
|
|
scale = float(singular[0])
|
|
audit.close("restricted_svd:unit_direction", norm(combination), 1.)
|
|
audit.close("restricted_svd:smallest_action_norm", norm(action), float(singular[-1]), scale)
|
|
audit.check("restricted_svd:action_reconstruction", vector_near(action, checked_einsum("ij,j->i",matrix,weights), scale))
|
|
audit.check("restricted_svd:factorization", vector_near(matrix, checked_einsum("ij,jk->ik",u*singular,vh), norm(matrix)))
|
|
audit.check("restricted_svd:finite", np.isfinite(singular).all() and np.isfinite(vh).all())
|
|
return {
|
|
"scope": "Jacobian restricted to nonzero normalized saved block directions; not full Jacobian rank or nullspace",
|
|
"basis": "disjoint block directions divided by their normalized coefficient two-norm",
|
|
"active_blocks": [blocks[i][0] for i in active], "singular_values": singular.tolist(),
|
|
"condition_ratio_in_restricted_span": ratio(float(singular[0]), float(singular[-1])),
|
|
"smallest_direction_norm": norm(combination), "smallest_action_norm": norm(action),
|
|
"coefficients": coefficients,
|
|
"reconstructed_direction_float64_sha256": hashlib.sha256(combination.astype("<f8").tobytes()).hexdigest(),
|
|
"reconstructed_action_float64_sha256": hashlib.sha256(action.astype("<f8").tobytes()).hexdigest(),
|
|
"smallest_singular_value_gap": float(singular[-2] - singular[-1]) if len(singular) > 1 else None,
|
|
"qualification": "basis metric is coefficient normalization; vector sign fixed, degenerate subspace basis not unique"
|
|
}, action
|
|
|
|
|
|
def surface_fit(surface, physical, physical_direction, blocks, audit):
|
|
shape = next(b for b in blocks if b[0] == "surface_deformation")
|
|
_, offset, size = shape
|
|
audit.check("surface:complete_parameter_indices", sorted(surface) == list(range(size)))
|
|
ordered = [surface[i] for i in range(size)]
|
|
audit.check("surface:unique_volume_dofs", len({integer(r["volume_scalar_true_dof"]) for r in ordered}) == size)
|
|
audit.check("surface:nonnegative_volume_dofs", all(integer(r["volume_scalar_true_dof"]) >= 0 for r in ordered))
|
|
xyz = np.array([[finite(r["unit_radial_" + axis]) for axis in "xyz"] for r in ordered])
|
|
data = np.array([finite(r["physical_direction_parameter"]) for r in ordered])
|
|
for i, row in enumerate(ordered):
|
|
r0, r = finite(row["reference_radius"]), finite(row["physical_radius"])
|
|
b = finite(row["base_parameter"])
|
|
audit.check(f"surface:radii:{i}", r0 > 0 and r > 0)
|
|
audit.close(f"surface:unit_radial:{i}", norm(xyz[i]), 1.)
|
|
audit.close(f"surface:radius_identity:{i}", r, r0 + b, abs(r0) + abs(b))
|
|
audit.close(f"surface:base:{i}", b, float(physical[offset+i]))
|
|
audit.close(f"surface:direction:{i}", data[i], float(physical_direction[offset+i]))
|
|
x, y, z = xyz.T
|
|
columns = [np.ones(size), x, y, z, x*y, x*z, y*z, x*x-y*y, 3*z*z-1]
|
|
fits = []
|
|
for degree, count in ((0, 1), (1, 4), (2, 9)):
|
|
matrix = np.column_stack(columns[:count])
|
|
coef, _, rank, singular = np.linalg.lstsq(matrix, data, rcond=None)
|
|
residual = data - checked_einsum("ij,j->i",matrix,coef)
|
|
audit.check(f"surface:fit_finite:{degree}", np.isfinite(coef).all() and np.isfinite(singular).all())
|
|
fits.append({"maximum_degree": degree, "rank": int(rank), "basis_coefficients": coef.tolist(),
|
|
"relative_nodal_l2_remainder": ratio(norm(residual), norm(data)),
|
|
"nodal_l2_remainder": norm(residual), "singular_values": singular.tolist()})
|
|
return {"scope": "unweighted nodal least squares in real l=0,1,2 polynomial spans; not area-weighted harmonic coefficients or gauge identification",
|
|
"basis": ["1", "nx", "ny", "nz", "nx*ny", "nx*nz", "ny*nz", "nx^2-ny^2", "3*nz^2-1"],
|
|
"direction_nodal_l2": norm(data), "direction_maximum_absolute": float(np.max(np.abs(data))), "fits": fits}
|
|
|
|
|
|
def physical_analysis(folder, physical, physical_direction, blocks, audit, required):
|
|
def read(suffix, key):
|
|
name = PP + suffix + ".csv"
|
|
required.add(name)
|
|
return table(folder / name, key)
|
|
metadata = read("metadata", lambda r: r["key"])
|
|
meta = {k: r["value"] for k, r in metadata.items()}
|
|
audit.check("physical:schema", meta["schema"] == "weak_shape_physical_v1")
|
|
audit.check("physical:declared_scalar_scope", "broken" in meta["dg_derivatives"] and "not evaluated" in meta["gravity_vector"])
|
|
elem = read("elements", lambda r: (integer(r["quadrature"]), integer(r["element"]), r["field"]))
|
|
geometry = read("geometry", lambda r: (integer(r["quadrature"]), integer(r["element"])))
|
|
affine = read("affine", lambda r: (integer(r["quadrature"]), integer(r["element"])))
|
|
surface = read("surface_nodes", lambda r: integer(r["parameter"]))
|
|
summary = read("summary", lambda r: (integer(r["quadrature"]), r["field"], r["metric"]))
|
|
observed = {key: finite(r["value"]) for key, r in summary.items()}
|
|
expected_summary = set()
|
|
def compare(q, field, metric, expected, scale=0.):
|
|
key = (q, field, metric)
|
|
expected_summary.add(key)
|
|
audit.close("physical:summary:" + str(key), observed[key], expected, scale)
|
|
audit.check("physical:quadrature_coverage", {q for q, _ in geometry} == {10, 14})
|
|
ids = {e for q, e in geometry if q == 10}
|
|
audit.check("physical:same_elements_at_both_rules", ids and ids == {e for q, e in geometry if q == 14})
|
|
audit.check("physical:affine_coverage", set(affine) == set(geometry))
|
|
audit.check("physical:field_coverage", set(elem) == {(q, e, f) for q, e in geometry for f in FIELDS})
|
|
outputs = []
|
|
for q in (10, 14):
|
|
geo = [geometry[(q, e)] for e in sorted(ids)]
|
|
aff = [affine[(q, e)] for e in sorted(ids)]
|
|
for e, g, a in zip(sorted(ids), geo, aff):
|
|
for key in g:
|
|
if key not in ("quadrature", "element"):
|
|
finite(g[key])
|
|
volume = finite(g["physical_volume"])
|
|
audit.check(f"physical:geometry_positive:{q}:{e}", volume > 0 and integer(g["points"]) > 0)
|
|
for prefix in ("reference", "physical"):
|
|
lo, hi = finite(g[prefix+"_radius_min"]), finite(g[prefix+"_radius_max"])
|
|
audit.check(f"physical:radius_range:{q}:{e}:{prefix}", 0 <= lo <= hi)
|
|
for key in ("base_displacement_squared_integral", "direction_displacement_squared_integral"):
|
|
audit.check(f"physical:geometry_square:{q}:{e}:{key}", finite(g[key]) >= 0)
|
|
for key, family in (("density_element_order", "density_family_order"),
|
|
("enthalpy_element_order", "enthalpy_family_order"),
|
|
("potential_element_order", "potential_family_order"),
|
|
("displacement_element_order", "displacement_family_order")):
|
|
audit.check(f"physical:order:{q}:{e}:{key}", integer(g[key]) == integer(meta[family]))
|
|
for key in ("reference_volume", "value_error_squared_integral", "reference_value_squared_integral",
|
|
"gradient_error_squared_integral", "reference_gradient_squared_integral"):
|
|
audit.check(f"physical:affine_nonnegative:{q}:{e}:{key}", finite(a[key]) >= 0)
|
|
audit.check(f"physical:affine_volume:{q}:{e}", finite(a["reference_volume"]) > 0)
|
|
audit.check(f"physical:affine_points:{q}:{e}", integer(a["points"]) == integer(g["points"]))
|
|
audit.close(f"physical:affine_identity:{q}:{e}", finite(a["reference_gradient_squared_integral"]), 3*finite(a["reference_volume"]))
|
|
totals_g = {k: math.fsum(finite(g[k]) for g in geo) for k in
|
|
("physical_volume", "base_displacement_squared_integral", "direction_displacement_squared_integral")}
|
|
for key, value in totals_g.items():
|
|
compare(q, "geometry", key, value)
|
|
compare(q, "geometry", "stellar_elements", len(ids))
|
|
totals_a = {k: math.fsum(finite(a[k]) for a in aff) for k in
|
|
("reference_volume", "value_error_squared_integral", "reference_value_squared_integral",
|
|
"gradient_error_squared_integral", "reference_gradient_squared_integral")}
|
|
for key, value in totals_a.items():
|
|
compare(q, "affine", key, value)
|
|
for label, num, den in (("relative_value_l2_error", "value_error_squared_integral", "reference_value_squared_integral"),
|
|
("relative_gradient_l2_error", "gradient_error_squared_integral", "reference_gradient_squared_integral")):
|
|
compare(q, "affine", label, math.sqrt(totals_a[num] / totals_a[den]))
|
|
field_outputs = {}
|
|
for field in FIELDS:
|
|
records = [elem[(q, e, field)] for e in sorted(ids)]
|
|
for e, row, g in zip(sorted(ids), records, geo):
|
|
values = {k: finite(row[k]) for k in INTEGRALS}
|
|
audit.check(f"physical:gradient_kind:{q}:{e}:{field}", row["gradient_kind"] == ("continuous_H1" if field == "specific_enthalpy" else "broken_DG"))
|
|
for key in ("attribute", "geometry", "points"):
|
|
audit.check(f"physical:element_binding:{q}:{e}:{field}:{key}", integer(row[key]) == integer(g[key]))
|
|
audit.close(f"physical:volume_binding:{q}:{e}:{field}", finite(row["physical_volume"]), finite(g["physical_volume"]))
|
|
audit.check(f"physical:squares:{q}:{e}:{field}", all(values[k] >= 0 for k in SQUARES))
|
|
for prefix in ("base","pullback","transport","eulerian"):
|
|
first,second=values[prefix+"_integral"],values[prefix+"_squared_integral"]
|
|
product=finite(row["physical_volume"])*second
|
|
audit.check(f"physical:moment_cauchy:{q}:{e}:{field}:{prefix}", first*first <= product+2048*EPS*(first*first+product)+1e-28)
|
|
# Signed element moments may themselves cancel. Cauchy bounds
|
|
# their absolute quadrature sums using the saved positive moments.
|
|
scale = math.sqrt(finite(row["physical_volume"])) * (math.sqrt(values["pullback_squared_integral"])+math.sqrt(values["transport_squared_integral"]))
|
|
audit.close(f"physical:signed_identity:{q}:{e}:{field}", values["eulerian_integral"], values["pullback_integral"]-values["transport_integral"], scale)
|
|
d2, t2, e2, cross = (values[k] for k in ("pullback_squared_integral", "transport_squared_integral", "eulerian_squared_integral", "pullback_transport_cross_integral"))
|
|
audit.close(f"physical:squared_identity:{q}:{e}:{field}", e2, d2+t2-2*cross, d2+t2+2*abs(cross))
|
|
audit.check(f"physical:cauchy:{q}:{e}:{field}", abs(cross) <= math.sqrt(d2)*math.sqrt(t2) + 2048*EPS*(d2+t2) + 1e-28)
|
|
totals = {k: math.fsum(finite(r[k]) for r in records) for k in INTEGRALS}
|
|
for key in (*SQUARES, "pullback_transport_cross_integral"):
|
|
compare(q, field, key, totals[key])
|
|
for prefix in ("pullback", "transport", "eulerian"):
|
|
compare(q, field, prefix+"_l2_norm", math.sqrt(totals[prefix+"_squared_integral"]))
|
|
d2, t2, e2, cross = (totals[k] for k in ("pullback_squared_integral", "transport_squared_integral", "eulerian_squared_integral", "pullback_transport_cross_integral"))
|
|
compare(q, field, "squared_norm_identity_error", e2-d2-t2+2*cross, d2+t2+2*abs(cross))
|
|
audit.close(f"physical:total_squared_identity:{q}:{field}", e2, d2+t2-2*cross, d2+t2+2*abs(cross))
|
|
# Each element spans a radius interval. Midpoint bins are explicitly
|
|
# element-group summaries, not integrals over exact spherical shells.
|
|
bins = []
|
|
outer = max(finite(g["physical_radius_max"]) for g in geo)
|
|
for low, high in ((0., .25), (.25, .5), (.5, .75), (.75, .9), (.9, 1.0000001)):
|
|
chosen = [i for i, g in enumerate(geo) if low <= .5*(finite(g["physical_radius_min"])+finite(g["physical_radius_max"]))/outer < high]
|
|
bins.append({"normalized_element_midradius_interval": [low, min(high, 1.)], "elements": len(chosen),
|
|
"eulerian_squared_integral": math.fsum(finite(records[i]["eulerian_squared_integral"]) for i in chosen),
|
|
"pullback_squared_integral": math.fsum(finite(records[i]["pullback_squared_integral"]) for i in chosen)})
|
|
audit.check(f"physical:radial_bin_coverage:{q}:{field}", sum(b["elements"] for b in bins) == len(ids))
|
|
top = sorted(zip(sorted(ids), records), key=lambda pair: finite(pair[1]["eulerian_squared_integral"]), reverse=True)[:10]
|
|
field_outputs[field] = {"integrals": totals, "eulerian_over_pullback_l2": ratio(math.sqrt(e2), math.sqrt(d2)),
|
|
"eulerian_over_pullback_plus_transport_l2": ratio(math.sqrt(e2), math.sqrt(d2)+math.sqrt(t2)),
|
|
"pullback_l2_per_unit_displacement_l2": ratio(math.sqrt(d2), math.sqrt(totals_g["direction_displacement_squared_integral"])),
|
|
"transport_l2_per_unit_displacement_l2": ratio(math.sqrt(t2), math.sqrt(totals_g["direction_displacement_squared_integral"])),
|
|
"eulerian_l2_per_unit_displacement_l2": ratio(math.sqrt(e2), math.sqrt(totals_g["direction_displacement_squared_integral"])),
|
|
"eulerian_rms": math.sqrt(e2/totals_g["physical_volume"]),
|
|
"pullback_transport_cosine": ratio(cross, math.sqrt(d2)*math.sqrt(t2)),
|
|
"radial_element_groups": bins,
|
|
"largest_eulerian_elements": [{"element": e, "squared_integral": finite(r["eulerian_squared_integral"]),
|
|
"fraction": ratio(finite(r["eulerian_squared_integral"]), e2)} for e, r in top]}
|
|
totals_g["direction_displacement_l2"] = math.sqrt(totals_g["direction_displacement_squared_integral"])
|
|
totals_g["direction_displacement_rms"] = math.sqrt(totals_g["direction_displacement_squared_integral"]/totals_g["physical_volume"])
|
|
outputs.append({"quadrature": q, "geometry": totals_g, "affine": totals_a, "fields": field_outputs})
|
|
for metric, expected in (("accepted_trial_prepared_normalization_unchanged", 1.), ("production_mapping_cache_unchanged", 1.),
|
|
("trial_preparations", 0.), ("volume_direction_actions", 1.)):
|
|
compare(0, "invariant", metric, expected)
|
|
for metric in ("maximum_scatter_roundtrip_error", "maximum_surface_trace_error"):
|
|
key = (0, "reconstruction", metric)
|
|
expected_summary.add(key)
|
|
audit.check("physical:"+metric, 0 <= observed[key] <= 512*EPS*max(1., float(np.max(np.abs(physical))), float(np.max(np.abs(physical_direction)))))
|
|
audit.check("physical:summary_coverage", set(observed) == expected_summary)
|
|
sensitivity = {}
|
|
for field in FIELDS:
|
|
sensitivity[field] = {}
|
|
for key in (*SQUARES, "pullback_transport_cross_integral"):
|
|
a, b = (out["fields"][field]["integrals"][key] for out in outputs)
|
|
sensitivity[field][key] = {"q10": a, "q14": b, "difference": b-a,
|
|
"relative_to_q14": ratio(abs(b-a), abs(b))}
|
|
return {"metadata": meta, "quadratures": outputs, "quadrature_sensitivity": sensitivity,
|
|
"surface": surface_fit(surface, physical, physical_direction, blocks, audit),
|
|
"scope": "recorded stellar broken-scalar integrals; no interface jumps, RT transport, exact shell integrals or continuum gauge test"}
|
|
|
|
|
|
def producer_binding(direction_path, folder, audit):
|
|
producer = direction_path.parent
|
|
mapping = {"nonlinear_step_physical_base.txt":"physical_base", "nonlinear_step_normalized_base.txt":"base",
|
|
"nonlinear_step_state_factors.txt":"state_factors", "nonlinear_step_residual_factors.txt":"residual_factors",
|
|
"nonlinear_step_residual.txt":"residual", "nonlinear_step_action_full.txt":"action"}
|
|
present = [name for name in mapping if (producer/name).exists()]
|
|
if not present:
|
|
return {"status":"producer_neighbors_unavailable", "bound":False}
|
|
if len(present) != len(mapping) or not (producer/"nonlinear_step_layout.csv").exists():
|
|
raise ValueError("Incomplete nonlinear-step direction producer artifacts")
|
|
saved_hashes, equal = {}, {}
|
|
for name, current in mapping.items():
|
|
path=producer/name
|
|
saved_hashes[name]=digest(path)
|
|
a,b=vector(path),vector(folder/(P+current+".txt"))
|
|
exact=np.array_equal(a,b)
|
|
equal[name]=exact
|
|
audit.check("producer:vector:"+name, exact if "factors" in name else vector_near(a,b))
|
|
path=producer/"nonlinear_step_layout.csv"
|
|
saved_hashes[path.name]=digest(path)
|
|
audit.check("producer:layout", rows(path)==rows(folder/(P+"layout.csv")))
|
|
run_path=producer.parent/(producer.name+"_run.json")
|
|
run_bound=False
|
|
if run_path.exists():
|
|
run=json.loads(run_path.read_text())
|
|
for name, actual in saved_hashes.items():
|
|
audit.check("producer:recorded_output_hash:"+name, run["output_hashes"].get(name)==actual)
|
|
audit.check("producer:direction_output_hash", run["output_hashes"].get(direction_path.name)==digest(direction_path))
|
|
saved_hashes[run_path.name]=digest(run_path)
|
|
run_bound=True
|
|
return {"status":"producer_neighbors_compared", "bound":run_bound, "producer":str(producer),
|
|
"hashes":saved_hashes,"vectors_bitwise_equal":equal,
|
|
"qualification":"neighbor hashes checked against producer run when present; no FEM rerun"}
|
|
|
|
|
|
def provenance(folder, required, audit):
|
|
root = folder.parent
|
|
source_file, run_file = root / "provenance.json", root / (folder.name + "_run.json")
|
|
if not source_file.exists() and not run_file.exists():
|
|
return {"status": "unavailable", "bound": False}
|
|
if not source_file.exists() or not run_file.exists():
|
|
raise ValueError("Partially present campaign provenance")
|
|
source, run = (json.loads(p.read_text()) for p in (source_file, run_file))
|
|
audit.check("provenance:case_id", run["id"] == folder.name)
|
|
command = run["command"]
|
|
final_mode = "--weak-shape-final" in command
|
|
audit.check("provenance:completed_probe_exit", run["exit_code"] == 0 or (final_mode and run["exit_code"] == 3))
|
|
audit.check("provenance:required_output_hashes", required <= set(run["output_hashes"]))
|
|
for name, expected in run["output_hashes"].items():
|
|
path = folder / name
|
|
audit.check("provenance:contained_output:"+name, path.resolve().is_relative_to(folder.resolve()))
|
|
audit.check("provenance:output_hash:"+name, digest(path) == expected)
|
|
for name, expected in run.get("input_state_hashes", {}).items():
|
|
audit.check("provenance:input_hash:"+name, digest(Path(name)) == expected)
|
|
case = next(c for c in source["cases"] if c["id"] == folder.name)
|
|
expected = [str(root / "plateau_diagnostic"), *case["arguments"], "--mesh", source["mesh"], "--output", str(folder)]
|
|
audit.check("provenance:command", command == expected)
|
|
audit.check("provenance:binary_hash", digest(Path(command[0])) == source["binary_sha256"])
|
|
audit.check("provenance:mesh_hash", digest(Path(source["mesh"])) == source["mesh_sha256"])
|
|
manifest = root / "cases.json"
|
|
if manifest.exists():
|
|
audit.check("provenance:case_manifest", json.loads(manifest.read_text()) == source["cases"])
|
|
meta_path=folder/"metadata.txt"
|
|
if meta_path.exists():
|
|
meta=dict(line.split("=",1) for line in meta_path.read_text().splitlines() if "=" in line)
|
|
audit.check("provenance:mode", ("--weak-shape-probe" in command) != final_mode)
|
|
audit.check("provenance:metadata_mode", meta["weak_shape_probe"]==str(int(not final_mode)) and meta["weak_shape_final"]==str(int(final_mode)))
|
|
audit.check("provenance:metadata_scope", meta["weak_shape_direction_space"]=="normalized_current_factors" and
|
|
meta["weak_shape_operator"]=="full_monolithic" and meta["weak_shape_trial_policy"]=="raw_prepare_only")
|
|
for option,key in (("--linear-tol","linear_tol"),("--restart","restart"),("--max-linear","max_linear"),
|
|
("--max-newton","max_newton"),("--index","n"),("--angular-momentum","J")):
|
|
if option in command:
|
|
audit.check("provenance:metadata_argument:"+option, finite(command[command.index(option)+1])==finite(meta[key]))
|
|
if final_mode:
|
|
audit.check("provenance:reported_convergence", integer(meta["converged"])==int(run["exit_code"]==0))
|
|
audit.check("provenance:final_accepted_state", np.array_equal(vector(folder/"accepted_state.txt"),vector(folder/(P+"physical_base.txt"))))
|
|
producer = None
|
|
for option, saved in (("--probe-state", "physical_base"), ("--probe-direction", "direction")):
|
|
if option in command:
|
|
input_path = Path(command[command.index(option)+1])
|
|
if not input_path.is_absolute():
|
|
input_path = Path(source["source"]) / input_path
|
|
input_path = input_path.resolve()
|
|
audit.check("provenance:input_bound:"+option, str(input_path) in run.get("input_state_hashes", {}))
|
|
audit.check("provenance:input_vector:"+option, np.array_equal(vector(input_path), vector(folder/(P+saved+".txt"))))
|
|
if option == "--probe-direction":
|
|
producer = producer_binding(input_path,folder,audit)
|
|
return {"status": "runner_hashes_checked", "bound": True, "source_commit": source["commit"],
|
|
"provenance_sha256": digest(source_file), "run_sha256": digest(run_file), "binary_sha256": source["binary_sha256"],
|
|
"direction_producer":producer}
|
|
|
|
|
|
def inspect(folder):
|
|
folder = Path(folder).resolve()
|
|
audit, required, observations = Audit(), set(), {}
|
|
try:
|
|
def get(name):
|
|
filename = P + name + ".txt"
|
|
required.add(filename)
|
|
return vector(folder / filename)
|
|
required.update(P+n+".csv" for n in ("summary", "layout", "blocks", "trials"))
|
|
summary = {k: finite(r["value"]) for k, r in table(folder/(P+"summary.csv"), lambda r:r["metric"]).items()}
|
|
for key, value in {"schema_version": 1, "accepted_commits": 0, "probe_exception": 0,
|
|
"accepted_buffers_bitwise_identical": 1, "normalization_factors_bitwise_identical": 1,
|
|
"direction_buffers_bitwise_identical": 1, "maximum_restore_error": 0,
|
|
"restored_action_error": 0}.items():
|
|
audit.check("summary:"+key, summary[key] == value)
|
|
x, xp, f, p, pp, jp, defect, factors, rf = [get(n) for n in
|
|
("base", "physical_base", "residual", "direction", "physical_direction", "action", "defect", "state_factors", "residual_factors")]
|
|
n = x.size
|
|
audit.check("vector_sizes", all(v.size == n for v in (xp,f,p,pp,jp,defect,factors,rf)))
|
|
audit.check("summary:state_size", summary["state_size"] == n)
|
|
audit.check("positive_factors", np.all(factors > 0) and np.all(rf > 0))
|
|
audit.check("denormalized_base", vector_near(xp, x/factors))
|
|
audit.check("denormalized_direction", vector_near(pp, p/factors))
|
|
audit.check("defect_identity", vector_near(defect, f+jp, norm(f)+norm(jp)))
|
|
audit.check("restored_residual_exact", np.array_equal(f, get("restored_residual")))
|
|
audit.check("restored_action_exact", np.array_equal(jp, get("restored_action")))
|
|
fn, pn, jn, dn = map(norm, (f,p,jp,defect))
|
|
if fn <= 0 or pn <= 0:
|
|
raise ValueError("Probe requires nonzero base residual and direction")
|
|
threshold = max(summary["outer_absolute_tolerance"], summary["outer_tolerance"]*fn)
|
|
audit.check("valid_linear_controls", 0 < summary["outer_tolerance"] < 1 and summary["outer_absolute_tolerance"] >= 0 and
|
|
integer(summary["maximum_iterations"]) > 0 and integer(summary["restart"]) > 0)
|
|
for key, value in {"base_norm":fn, "direction_norm":pn, "action_norm":jn, "verified_defect_norm":dn,
|
|
"verified_relative_defect":dn/fn, "verified_threshold":threshold}.items():
|
|
audit.close("summary:"+key, summary[key], value)
|
|
audit.check("verified_gate_truth", summary["verified_gate"] in (0,1) and bool(summary["verified_gate"]) == (dn <= threshold))
|
|
audit.check("supplied_direction_flag", summary["supplied_direction"] in (0,1))
|
|
if not summary["supplied_direction"]:
|
|
audit.close("linear_reported_true_norm", summary["linear_reported_true_norm"], dn, fn+jn)
|
|
audit.check("linear_iterations", 0 <= integer(summary["linear_iterations"]) <= integer(summary["maximum_iterations"]))
|
|
blocks = layout(folder/(P+"layout.csv"), n)
|
|
actions = [get("action_b"+str(i)) for i in range(len(blocks["value"]))]
|
|
audit.check("split_column_count", summary["split_columns"] == len(actions))
|
|
audit.check("direct_action_count", summary["direct_jacobian_applications"] == len(actions)+2)
|
|
audit.check("preparation_count", summary["preparations"] == 6 + int(not summary["supplied_direction"]))
|
|
audit.check("finite_nonnegative_probe_time", summary["total_probe_seconds"] >= 0)
|
|
total = np.sum(actions, axis=0)
|
|
action_scale = math.fsum(norm(a) for a in actions)
|
|
audit.check("split_action_sum", vector_near(total, jp, action_scale))
|
|
audit.close("reported_split_error", summary["split_action_sum_error"], norm(total-jp), action_scale)
|
|
reported = table(folder/(P+"blocks.csv"), lambda r:(integer(r["column"]),r["residual_block"]))
|
|
audit.check("block_table_coverage", set(reported) == {(i,r[0]) for i in range(len(actions)) for r in blocks["residual"]})
|
|
for i,(name,start,count) in enumerate(blocks["value"]):
|
|
for rname,offset,size in blocks["residual"]:
|
|
row = reported[(i,rname)]
|
|
audit.check(f"block_metadata:{i}:{rname}", row["value_block"]==name and integer(row["value_offset"])==start and
|
|
integer(row["value_size"])==count and integer(row["residual_offset"])==offset and integer(row["residual_size"])==size)
|
|
audit.close(f"block_direction_norm:{i}:{rname}", finite(row["direction_norm"]), norm(p[start:start+count]))
|
|
audit.close(f"block_action_norm:{i}:{rname}", finite(row["action_norm"]), norm(actions[i][offset:offset+size]))
|
|
svd, smallest_action = restricted_svd(p, actions, blocks["value"], audit)
|
|
svd["smallest_action_residual_blocks"] = {name:norm(smallest_action[o:o+s]) for name,o,s in blocks["residual"]}
|
|
trial_rows = table(folder/(P+"trials.csv"), lambda r:integer(r["trial"]))
|
|
audit.check("trial_coverage", set(trial_rows)=={0,1,2})
|
|
trials = []
|
|
for i,alpha in enumerate((1/128,1/512,1/2048)):
|
|
row = trial_rows[i]
|
|
audit.check(f"trial_alpha:{i}", finite(row["alpha"])==alpha)
|
|
audit.check(f"trial_prepared_flag:{i}", row["prepared"] in ("0","1"))
|
|
audit.check(f"trial_candidate:{i}", vector_near(get("candidate_a"+str(i)), x+alpha*p, norm(x)+alpha*pn))
|
|
prediction = f+alpha*jp
|
|
audit.close(f"trial_prediction_norm:{i}", finite(row["linear_prediction_norm"]), norm(prediction))
|
|
prepared = row["prepared"]=="1"
|
|
entry = {"trial":i, "alpha":alpha, "prepared":prepared, "reason":integer(row["reason"]), "stage":integer(row["stage"])}
|
|
if prepared:
|
|
audit.check(f"trial_prepared_reasons:{i}", entry["reason"] == entry["stage"] == -1)
|
|
ft = get("residual_a"+str(i)); remainder = ft-prediction
|
|
values = {"residual_norm":norm(ft), "remainder_norm":norm(remainder), "remainder_over_alpha":norm(remainder)/alpha}
|
|
for key,value in values.items():
|
|
audit.close(f"trial:{i}:{key}", finite(row[key]), value, (fn+alpha*jn+norm(ft)) / (alpha if key=="remainder_over_alpha" else 1))
|
|
entry.update(values)
|
|
entry.update(residual_over_base=norm(ft)/fn, remainder_over_alpha_squared=norm(remainder)/alpha**2,
|
|
linear_dot_remainder=float(checked_einsum("i,i->",prediction,remainder)))
|
|
else:
|
|
audit.check(f"trial_failed_reasons:{i}", entry["reason"]>=0 and entry["stage"]>=0)
|
|
audit.check(f"trial_failed_vector_absent:{i}", not (folder/(P+"residual_a"+str(i)+".txt")).exists())
|
|
audit.check(f"trial_failed_metrics:{i}", all(math.isnan(float(row[k])) for k in ("residual_norm","remainder_norm","remainder_over_alpha")))
|
|
trials.append(entry)
|
|
physical_data = physical_analysis(folder,xp,pp,blocks["value"],audit,required)
|
|
if (folder/"metadata.txt").exists():
|
|
required.add("metadata.txt")
|
|
meta=dict(line.split("=",1) for line in (folder/"metadata.txt").read_text().splitlines() if "=" in line)
|
|
for key,source in (("outer_tolerance","linear_tol"),("restart","restart"),("maximum_iterations","max_linear")):
|
|
audit.check("metadata:linear_control:"+key,summary[key]==finite(meta[source]))
|
|
if meta.get("weak_shape_final")=="1":
|
|
required.add("accepted_state.txt")
|
|
binding = provenance(folder,required,audit)
|
|
observations = {"summary":summary, "provenance":binding,
|
|
"linear_classification":"verified_newton_direction" if dn<=threshold else "finite_direction_failed_true_gate",
|
|
"verified_linear_gate":dn<=threshold, "restricted_block_span_svd":svd, "trials":trials,
|
|
"block_action_cancellation_factor":ratio(action_scale,jn), "physical":physical_data,
|
|
"input_vector_hashes":{name:digest(folder/(P+name+".txt")) for name in ("physical_base","base","direction","state_factors","residual_factors")},
|
|
"layout_sha256":digest(folder/(P+"layout.csv"))}
|
|
except (OSError, ValueError, KeyError, TypeError, IndexError, StopIteration, OverflowError, ZeroDivisionError, np.linalg.LinAlgError) as error:
|
|
audit.check("complete_consistent_artifacts: " + str(error),False)
|
|
passed = not audit.failures
|
|
return {"directory":str(folder), "artifact_pass":passed, "observations_valid":passed,
|
|
"checks_count":audit.count, "failures":audit.failures,
|
|
"required_artifacts":sorted(required), "observations":observations if passed else {}}
|
|
|
|
|
|
def comparisons(cases):
|
|
output=[]
|
|
for i,left in enumerate(cases):
|
|
for right in cases[i+1:]:
|
|
if not left["artifact_pass"] or not right["artifact_pass"]:
|
|
continue
|
|
a,b=left["observations"],right["observations"]
|
|
same_layout=a["layout_sha256"]==b["layout_sha256"]
|
|
same={k:a["input_vector_hashes"][k]==b["input_vector_hashes"][k] for k in a["input_vector_hashes"]}
|
|
output.append({"left":left["directory"],"right":right["directory"],"same_layout":same_layout,"same_saved_vectors":same,
|
|
"same_state_metric_comparison":same_layout and all(same[k] for k in ("physical_base","base","state_factors","residual_factors")),
|
|
"qualification":"different discretizations or factors require physical comparisons; raw coefficient norms are not comparable"})
|
|
return output
|
|
|
|
|
|
def tiny_fixture(folder, failed_gate=False):
|
|
"""Construct a complete eight-coordinate synthetic producer, with no FEM."""
|
|
folder.mkdir(parents=True,exist_ok=True)
|
|
def save(name,v):
|
|
(folder/(P+name+".txt")).write_text(str(len(v))+"\n"+"\n".join(format(float(x),".17g") for x in v)+"\n")
|
|
def write(name,records):
|
|
with (folder/name).open("w",newline="") as stream:
|
|
out=csv.DictWriter(stream,fieldnames=list(records[0]));out.writeheader();out.writerows(records)
|
|
sf=np.arange(1.,9.);xp=np.array([2.,0.,0.,0.,0.,0.,0.,1.]);x=xp*sf
|
|
p=np.arange(1.,9.)*1e-4;f=-p;jp=p*(.9 if failed_gate else 1.)
|
|
for name,v in {"base":x,"physical_base":xp,"residual":f,"direction":p,"physical_direction":p/sf,
|
|
"action":jp,"defect":f+jp,"state_factors":sf,"residual_factors":np.ones(8),
|
|
"restored_residual":f,"restored_action":jp}.items():save(name,v)
|
|
blocks=[("density",0,1),("surface_deformation",1,6),("constant",7,1)]
|
|
write(P+"layout.csv",[{"kind":kind,"block":name,"offset":o,"size":n} for kind in ("value","residual") for name,o,n in blocks])
|
|
br=[]
|
|
for i,(name,o,n) in enumerate(blocks):
|
|
a=np.zeros(8);a[o:o+n]=jp[o:o+n];save("action_b"+str(i),a)
|
|
for rn,ro,rs in blocks:
|
|
br.append(dict(column=i,value_block=name,value_offset=o,value_size=n,residual_block=rn,residual_offset=ro,
|
|
residual_size=rs,direction_norm=norm(p[o:o+n]),action_norm=norm(a[ro:ro+rs])))
|
|
write(P+"blocks.csv",br)
|
|
summary=dict(schema_version=1,state_size=8,supplied_direction=1,accepted_commits=0,outer_tolerance=.03,
|
|
outer_absolute_tolerance=0,maximum_iterations=200,restart=160,base_norm=norm(f),direction_norm=norm(p),
|
|
action_norm=norm(jp),verified_defect_norm=norm(f+jp),verified_relative_defect=norm(f+jp)/norm(f),
|
|
verified_threshold=.03*norm(f),verified_gate=int(not failed_gate),split_columns=3,split_action_sum_error=0,
|
|
restored_action_error=0,probe_exception=0,accepted_buffers_bitwise_identical=1,normalization_factors_bitwise_identical=1,
|
|
direction_buffers_bitwise_identical=1,maximum_restore_error=0,preparations=6,direct_jacobian_applications=5,total_probe_seconds=1.)
|
|
write(P+"summary.csv",[dict(metric=k,value=v) for k,v in summary.items()])
|
|
tr=[]
|
|
for i,a in enumerate((1/128,1/512,1/2048)):
|
|
save("candidate_a"+str(i),x+a*p);pred=f+a*jp;ft=pred+a*a*np.ones(8)*1e-5;save("residual_a"+str(i),ft)
|
|
tr.append(dict(trial=i,alpha=a,prepared=1,reason=-1,stage=-1,residual_norm=norm(ft),remainder_norm=norm(ft-pred),
|
|
remainder_over_alpha=norm(ft-pred)/a,linear_prediction_norm=norm(pred)))
|
|
write(P+"trials.csv",tr)
|
|
meta=dict(schema="weak_shape_physical_v1",dg_derivatives="broken DG derivatives",gravity_vector="not evaluated",
|
|
density_family_order="2",enthalpy_family_order="3",potential_family_order="2",displacement_family_order="3")
|
|
write(PP+"metadata.csv",[dict(key=k,value=v) for k,v in meta.items()])
|
|
elements=[];geometry=[];affine=[];ps=[]
|
|
for q in (10,14):
|
|
g=dict(quadrature=q,element=0,attribute=1,geometry=5,points=8,reference_mesh_geometry_order=4,density_element_order=2,
|
|
enthalpy_element_order=3,potential_element_order=2,displacement_element_order=3,reference_radius_min=.1,
|
|
reference_radius_max=.9,physical_radius_min=.1,physical_radius_max=.9,physical_volume=1.,
|
|
base_displacement_squared_integral=0.,direction_displacement_squared_integral=.2)
|
|
geometry.append(g)
|
|
ar=dict(quadrature=q,element=0,points=8,reference_volume=1.,value_error_squared_integral=.01,
|
|
reference_value_squared_integral=1.,gradient_error_squared_integral=.03,reference_gradient_squared_integral=3.)
|
|
affine.append(ar)
|
|
for field in FIELDS:
|
|
vals=dict(zip(INTEGRALS,(2.,3.,2.,1.,4.,9.,4.,1.,6.)))
|
|
elements.append(dict(quadrature=q,element=0,attribute=1,geometry=5,field=field,
|
|
gradient_kind="continuous_H1" if field=="specific_enthalpy" else "broken_DG",points=8,physical_volume=1.,**vals))
|
|
measured={k:vals[k] for k in (*SQUARES,"pullback_transport_cross_integral")}
|
|
measured.update(pullback_l2_norm=3.,transport_l2_norm=2.,eulerian_l2_norm=1.,squared_norm_identity_error=0.)
|
|
ps.extend(dict(quadrature=q,field=field,metric=k,value=v) for k,v in measured.items())
|
|
ps.extend(dict(quadrature=q,field="geometry",metric=k,value=v) for k,v in
|
|
dict(stellar_elements=1,physical_volume=1.,base_displacement_squared_integral=0.,direction_displacement_squared_integral=.2).items())
|
|
am={k:v for k,v in ar.items() if k not in ("quadrature","element","points")}
|
|
am.update(relative_value_l2_error=.1,relative_gradient_l2_error=.1)
|
|
ps.extend(dict(quadrature=q,field="affine",metric=k,value=v) for k,v in am.items())
|
|
ps.extend(dict(quadrature=0,field="invariant",metric=k,value=v) for k,v in
|
|
dict(accepted_trial_prepared_normalization_unchanged=1,production_mapping_cache_unchanged=1,trial_preparations=0,volume_direction_actions=1).items())
|
|
ps.extend(dict(quadrature=0,field="reconstruction",metric=k,value=0.) for k in
|
|
("maximum_scatter_roundtrip_error","maximum_surface_trace_error"))
|
|
for name,data in (("elements",elements),("geometry",geometry),("affine",affine),("summary",ps)):write(PP+name+".csv",data)
|
|
xyz=[[1.,0.,0.],[-1.,0.,0.],[0.,1.,0.],[0.,-1.,0.],[0.,0.,1.],[0.,0.,-1.]]
|
|
write(PP+"surface_nodes.csv",[dict(parameter=i,volume_scalar_true_dof=i,unit_radial_x=v[0],unit_radial_y=v[1],
|
|
unit_radial_z=v[2],reference_radius=1.,physical_radius=1.,base_parameter=0.,physical_direction_parameter=p[i+1]/sf[i+1]) for i,v in enumerate(xyz)])
|
|
|
|
|
|
def self_test():
|
|
"""Tiny algebraic controls only; never opens campaign artifacts."""
|
|
import tempfile
|
|
count=0
|
|
def require(condition):
|
|
nonlocal count
|
|
assert condition
|
|
count+=1
|
|
audit=Audit()
|
|
p=np.array([2.,3.,0.])
|
|
actions=[np.array([2.,0.,0.]),np.array([-3.,3e-8,0.]),np.zeros(3)]
|
|
out,action=restricted_svd(p,actions,[("a",0,1),("b",1,1),("zero",2,1)],audit)
|
|
require(not audit.failures)
|
|
require(near(out["singular_values"][-1],1e-8/math.sqrt(2),1e-8))
|
|
require(out["active_blocks"]==["a","b"])
|
|
require(norm(action)<1e-7)
|
|
bad=Audit();restricted_svd(p,actions[:2]+[np.ones(3)],[("a",0,1),("b",1,1),("zero",2,1)],bad)
|
|
require(bool(bad.failures))
|
|
require(not vector_near(np.array([1.,2.]),np.array([1.,2.001])))
|
|
require(near(1e-20,0.,1.))
|
|
# This size triggers spurious Accelerate matmul flags on the investigation
|
|
# host. Use independent scalar sums, not another BLAS product, as the oracle.
|
|
a=np.linspace(-1.,1.,257*6).reshape(257,6);v=np.arange(1.,7.)
|
|
expected=np.array([math.fsum(float(x*y) for x,y in zip(row,v)) for row in a])
|
|
with np.errstate(divide="raise",over="raise",invalid="raise"):
|
|
require(vector_near(checked_einsum("ij,j->i",a,v),expected,norm(expected)))
|
|
b=np.arange(18.).reshape(6,3)/18.
|
|
expected_matrix=np.array([[math.fsum(float(a[i,k]*b[k,j]) for k in range(6)) for j in range(3)] for i in range(257)])
|
|
require(vector_near(checked_einsum("ij,jk->ik",a,b),expected_matrix,norm(expected_matrix)))
|
|
try:checked_einsum("i,i->",np.array([math.inf]),np.ones(1));raise AssertionError("nonfinite contraction accepted")
|
|
except ValueError:require(True)
|
|
with tempfile.TemporaryDirectory() as directory:
|
|
d=Path(directory)
|
|
(d/"v").write_text("3\n1\n2\n3\n")
|
|
require(np.array_equal(vector(d/"v"),np.array([1.,2.,3.])))
|
|
for text in ("2\n1\n", "2\n1\nnan\n", "2\n1\ninf\n"):
|
|
(d/"v").write_text(text)
|
|
try:vector(d/"v");raise AssertionError("invalid vector accepted")
|
|
except ValueError:require(True)
|
|
(d/"layout").write_text("kind,block,offset,size\nvalue,a,0,2\nresidual,b,0,2\n")
|
|
require(layout(d/"layout",2)["value"]==[("a",0,2)])
|
|
(d/"layout").write_text("kind,block,offset,size\nvalue,a,1,2\nresidual,b,0,2\n")
|
|
try:layout(d/"layout",2);raise AssertionError("bad layout accepted")
|
|
except ValueError:require(True)
|
|
xyz=np.array([[1.,0.,0.],[-1.,0.,0.],[0.,1.,0.],[0.,-1.,0.],[0.,0.,1.],[0.,0.,-1.]])
|
|
values=2.+3.*xyz[:,0]
|
|
surf={i:{"volume_scalar_true_dof":str(i),"unit_radial_x":str(v[0]),"unit_radial_y":str(v[1]),"unit_radial_z":str(v[2]),
|
|
"reference_radius":"1","physical_radius":"1","base_parameter":"0","physical_direction_parameter":str(values[i])} for i,v in enumerate(xyz)}
|
|
audit=Audit();fits=surface_fit(surf,np.zeros(6),values,[("surface_deformation",0,6)],audit)
|
|
require(not audit.failures)
|
|
require(fits["fits"][1]["relative_nodal_l2_remainder"]<1e-14)
|
|
require(fits["fits"][0]["relative_nodal_l2_remainder"]>.1)
|
|
with tempfile.TemporaryDirectory() as directory:
|
|
d=Path(directory).resolve()/"case"
|
|
tiny_fixture(d)
|
|
result=inspect(d)
|
|
assert result["artifact_pass"],result["failures"]
|
|
require(result["observations"]["verified_linear_gate"])
|
|
require(result["observations"]["provenance"]["status"]=="unavailable")
|
|
require(result["observations"]["physical"]["quadratures"][0]["fields"]["density"]["eulerian_over_pullback_plus_transport_l2"]==.2)
|
|
tiny_fixture(d,failed_gate=True)
|
|
result=inspect(d)
|
|
require(result["artifact_pass"] and result["observations"]["linear_classification"]=="finite_direction_failed_true_gate")
|
|
for name in ("defect.txt","restored_action.txt","candidate_a1.txt","physical_base.txt","action_b1.txt"):
|
|
tiny_fixture(d)
|
|
path=d/(P+name);lines=path.read_text().splitlines();lines[1]=str(float(lines[1])+.001);path.write_text("\n".join(lines)+"\n")
|
|
require(not inspect(d)["artifact_pass"])
|
|
for filename,key in (("elements","pullback_transport_cross_integral"),("summary","value"),("geometry","physical_volume")):
|
|
tiny_fixture(d);path=d/(PP+filename+".csv");data=rows(path);data[0][key]=str(float(data[0][key])+.1)
|
|
with path.open("w",newline="") as stream:
|
|
writer=csv.DictWriter(stream,fieldnames=list(data[0]));writer.writeheader();writer.writerows(data)
|
|
require(not inspect(d)["artifact_pass"])
|
|
for missing in (False,True):
|
|
tiny_fixture(d);path=d/(PP+"summary.csv");data=rows(path)
|
|
if missing:data=[r for r in data if r["metric"]!="production_mapping_cache_unchanged"]
|
|
else:
|
|
next(r for r in data if r["metric"]=="production_mapping_cache_unchanged")["value"]="0"
|
|
with path.open("w",newline="") as stream:
|
|
writer=csv.DictWriter(stream,fieldnames=list(data[0]));writer.writeheader();writer.writerows(data)
|
|
require(not inspect(d)["artifact_pass"])
|
|
tiny_fixture(d)
|
|
path=d/(P+"trials.csv");data=rows(path);data[0].update(prepared="0",reason="2",stage="3",
|
|
residual_norm="nan",remainder_norm="nan",remainder_over_alpha="nan")
|
|
with path.open("w",newline="") as stream:
|
|
writer=csv.DictWriter(stream,fieldnames=list(data[0]));writer.writeheader();writer.writerows(data)
|
|
(d/(P+"residual_a0.txt")).unlink()
|
|
require(inspect(d)["artifact_pass"])
|
|
tiny_fixture(d)
|
|
root=d.parent
|
|
executable=root/"plateau_diagnostic";executable.write_text("synthetic binary fixture")
|
|
mesh=root/"mesh";mesh.write_text("synthetic mesh fixture")
|
|
case={"id":"case","arguments":["--weak-shape-probe"]}
|
|
prov={"commit":"0"*40,"source":str(root),"cases":[case],"mesh":str(mesh),"mesh_sha256":digest(mesh),"binary_sha256":digest(executable)}
|
|
(root/"provenance.json").write_text(json.dumps(prov))
|
|
run={"id":"case","exit_code":0,"input_state_hashes":{},"command":[str(executable),"--weak-shape-probe","--mesh",str(mesh),"--output",str(d)],
|
|
"output_hashes":{p.name:digest(p) for p in d.iterdir()}}
|
|
(root/"case_run.json").write_text(json.dumps(run))
|
|
result=inspect(d)
|
|
assert result["artifact_pass"],result["failures"]
|
|
require(result["observations"]["provenance"]["bound"])
|
|
run["output_hashes"][P+"action.txt"]="0"*64
|
|
(root/"case_run.json").write_text(json.dumps(run))
|
|
require(not inspect(d)["artifact_pass"])
|
|
producer=root/"producer";producer.mkdir()
|
|
for old,new in {"physical_base":"physical_base","normalized_base":"base","state_factors":"state_factors",
|
|
"residual_factors":"residual_factors","residual":"residual","action_full":"action","direction_full":"direction"}.items():
|
|
(producer/("nonlinear_step_"+old+".txt")).write_bytes((d/(P+new+".txt")).read_bytes())
|
|
(producer/"nonlinear_step_layout.csv").write_bytes((d/(P+"layout.csv")).read_bytes())
|
|
(root/"producer_run.json").write_text(json.dumps({"output_hashes":{p.name:digest(p) for p in producer.iterdir()}}))
|
|
audit=Audit();bound=producer_binding(producer/"nonlinear_step_direction_full.txt",d,audit)
|
|
require(bound["bound"] and not audit.failures)
|
|
(producer/"nonlinear_step_state_factors.txt").write_text("8\n"+"2\n"*8)
|
|
audit=Audit();producer_binding(producer/"nonlinear_step_direction_full.txt",d,audit)
|
|
require(bool(audit.failures))
|
|
print(f"{count} tiny synthetic checks passed")
|
|
|
|
|
|
def main():
|
|
parser=argparse.ArgumentParser(description=__doc__)
|
|
parser.add_argument("directories",nargs="*",type=Path)
|
|
parser.add_argument("--output",type=Path)
|
|
parser.add_argument("--self-test",action="store_true")
|
|
args=parser.parse_args()
|
|
if args.self_test:
|
|
self_test();return 0
|
|
if not args.directories or args.output is None:
|
|
parser.error("Provide one or more probe directories and --output")
|
|
destination=args.output.resolve()
|
|
if destination.exists() or any(destination.is_relative_to(d.resolve()) for d in args.directories):
|
|
parser.error("Use a fresh output JSON outside the immutable probe directories")
|
|
cases=[inspect(d) for d in args.directories]
|
|
result={"schema_version":1,"passed":all(c["artifact_pass"] for c in cases),"checker_sha256":digest(Path(__file__)),
|
|
"numpy_version":np.__version__,"tolerances":{"relative":RTOL,"cancellation_roundoff_multiplier":2048,"absolute_floor":1e-28},
|
|
"scope":"Independent recorded-vector algebra/hash and physical-integral identity audit; no FEM rerun, quadrature rerun, full-rank, gauge, convergence or physical-model certificate",
|
|
"cases":cases,"comparisons":comparisons(cases)}
|
|
destination.parent.mkdir(parents=True,exist_ok=True)
|
|
destination.write_text(json.dumps(result,indent=2,allow_nan=False)+"\n")
|
|
print(json.dumps({"passed":result["passed"],"cases":len(cases),"output":str(destination)}))
|
|
return 0 if result["passed"] else 1
|
|
|
|
|
|
if __name__=="__main__":
|
|
sys.exit(main())
|