"""Automated grader for the mixed-derivative ADI benchmark. Usage: python grade.py [--budget SECONDS] [--no-exec] Required files in the output dir: adi_solver.py submitted solver (REQUIRED; imported and executed) reference_arrays.npz numerical arrays results.json summary sidecar Scoring: 1.0 Tier 1 + Tier 2 + Tier 3 all pass 0.5 Tier 1 + Tier 2 (and all provenance checks) pass, Tier 3 missing/fails 0.0 required files missing, schema/consistency violations, wrong mesh sequence, instability, disallowed dependencies, a non-ADI linear solve, a solver that cannot reproduce its own submitted convergence errors, or Tier 1 / Tier 2 fail Design notes ------------ The grader is self-contained: it knows the diffusion matrix and the manufactured exact solution, recomputes the exact solution analytically, recomputes every order of convergence from the submitted arrays, and re-runs the submitted solver. *Scheme-agnostic thresholds.* The task permits ANY ADI scheme that demonstrates stable second-order convergence on the prescribed mesh sequence. The acceptance criteria are therefore stated as absolute accuracy and convergence requirements, never as a band around one privileged reference method, and every one of them is published verbatim in problem_spec.md. They were calibrated against the Craig-Sneyd / Modified Craig-Sneyd / Hundsdorfer-Verwer family; the measured coverage is MCS for theta in [1/3, 8], CS at theta = 1/2, and HV for theta in [1/2, 8] (see scheme_survey.py, which exits non-zero if this claim ever diverges from the thresholds below). Schemes that are second order only ASYMPTOTICALLY, whose asymptotic regime lies beyond h = 0.00625 (MCS at theta >= 10), are rejected because the task requires the order to approach 2 on exactly this mesh sequence. The classic failure modes -- mixed term dropped, mixed term lagged at first order (Douglas), Craig-Sneyd away from theta = 1/2, and unstable parameter choices -- fail on several criteria at once. *Provenance.* Submitted arrays are worthless unless they came from the submitted code, so the solver is imported in an isolated subprocess and used to regenerate the study. Two distinct checks, with two distinct consequences: * the regenerated max-norm ERRORS must match the submitted convergence table. A mismatch means the table was not produced by this solver -> 0.0. * the regenerated finest FIELD must match the submitted ``U_finest`` pointwise, not merely by its scalar max-norm. A mismatch means the Tier 3 artifact was not produced by the same run -> Tier 3 fails -> 0.5, which is how every other Tier 3 defect is scored. That subprocess also instruments every linear-algebra entry point and records the dimension and bandwidth of each system solved, which is how the "1-D tridiagonal sweeps only" rule is enforced: genuine ADI on an m x m grid solves tridiagonal systems of dimension O(m), while the forbidden monolithic 2-D solve is O(m^2). It additionally poisons intercepted solve results on a hidden general-data run; an unchanged finite field proves that the calls were discarded dummy work. """ import argparse import ast import json import os import re import subprocess import sys import tempfile import numpy as np # ---- problem ground truth (must match adi_solver.py) ------------------------ A_COEF, B_COEF, C_COEF = 1.0, 0.8, 0.5 MODES = [(1.0, 0.7, 1.3, 0.0), (0.7, 1.4, -0.6, 0.3), (0.5, 1.0, 0.8, 0.7)] H_VALUES = [0.1 * 2.0**(-p) for p in range(5)] T_FINAL = 1.0 FINEST_POINTS = 161 # ---- acceptance thresholds -------------------------------------------------- # ALL of these are published verbatim in problem_spec.md, so the agent knows the # bar and no scheme is privileged by a hidden band. See scheme_survey.py for the # calibration and EVALUATION.md for the exact coverage claim. # # The two magnitude caps below are NONSENSE GUARDS, not discriminators. They are # set far above every scheme that converges at second order on this mesh # sequence (worst surveyed: e(0.1) = 2.09e-1, e(finest) = 2.62e-3 for MCS at # theta = 10). Discrimination comes entirely from the convergence criteria. TIER1_MAX_ERR = 1.0 # finite and not absurd FINEST_MAX_ERR = 1.0e-2 # finite and not absurd # Total error reduction across the 16x refinement. Second order implies ~256x. # Accepted schemes surveyed: 112..281. First-order Douglas: 18.5..32.5. # Craig-Sneyd off its second-order parameter: 28.9..51.5. Dropped term: 1.2. MIN_TOTAL_REDUCTION = 100.0 # Least-squares fitted order over all five meshes. # Accepted: 1.772..2.033. Douglas: 1.006..1.247. CS off-parameter: 1.205..1.410. LS_ORDER_LO, LS_ORDER_HI = 1.70, 2.30 # The two finest pointwise experimental orders. # The LOWER bound does the discriminating work: accepted worst case 1.680 (HV at # theta = 5) vs Douglas 0.888..1.068, Craig-Sneyd off-parameter 1.057..1.254 and # the dropped mixed term 0.004..0.015. # The UPPER bound is only a sanity guard -- measured over the whole surveyed # family it rejects nothing except large-theta MCS, whose pre-asymptotic # transient can overshoot (theta = 8 peaks at 2.446). It is therefore set wide # enough to admit that overshoot rather than to exclude a stable, second-order # scheme on a technicality; such schemes are still caught, when they deserve to # be, by the total-reduction and fitted-order criteria. PT_ORDER_LO, PT_ORDER_HI = 1.60, 2.60 # Provenance: regenerated errors must match the submitted ones. REPRO_RTOL = 1.0e-6 # Pointwise provenance tolerance for the finest submitted numerical field. # This is deliberately much tighter than the convergence thresholds: the same # submitted solver is being rerun in the same pinned environment. SURFACE_RTOL, SURFACE_ATOL = 1.0e-10, 1.0e-12 # A genuine 1-D sweep on an m x m grid has dimension m-1 (or m+1 if the boundary # nodes are carried); the forbidden monolithic 2-D system has dimension (m-1)^2. # The factor 2 leaves room for legitimate variants while still excluding the # monolithic solve by a wide margin (at m = 10: limit 22, forbidden 81). SWEEP_DIM_LIMIT = lambda m: 2 * (m + 1) # noqa: E731 # Any implicit solve must be tridiagonal, per the specification. MAX_BANDWIDTH = 1 DEFAULT_BUDGET = 600.0 # seconds per regeneration phase COARSE_H = H_VALUES[:3] # regenerated before Tier 1/Tier 2 credit FINE_H = H_VALUES[3:] # regenerated before Tier 3 credit ALLOWED_TOP_LEVEL = set(sys.stdlib_module_names) | {"numpy", "scipy"} # ---- hidden instance used to prove the solver actually solves the PDE ------- # Observing a correctly sized tridiagonal call does NOT establish that the call # produced the returned field: a solver can make one legal sweep, discard it, and # return the manufactured analytic solution plus a synthetic O(h^2) error. The # only robust answer is to call solve() on a problem whose solution the agent # cannot write down, and check that it converges to an independently computed # reference. These coefficients and this initial state appear nowhere in the # solver-visible specification. # Weak enough that the solution is still O(0.1) at T = 1 -- with strong # diffusion it decays to numerical zero and the comparison becomes vacuous # (a solver returning zeros would "converge"). Correlation c/sqrt(ab) = # -0.51, so the mixed term remains as influential as in the benchmark. HIDDEN = {"a": 0.12, "b": 0.08, "c": -0.05, "p1": 0.7, "p2": -0.4, "p3": 0.9} HIDDEN_H = [0.1, 0.05] # meshes the agent is asked to solve HIDDEN_H_REF = 0.00625 # the grader's own reference mesh (16x / 8x finer) # Calibrated: the gold solver scores 2.0e-3 -> 6.3e-4 (ratio 3.19) on this # instance, while a solver returning zeros -- or the benchmark's manufactured # field, or the initial state -- is off by ~2.2e-1. The cap sits ~48x above the # gold error and ~7x below the cheapest cheat, and the ratio test separates a # genuine solver (~3-4) from any h-independent field (~1). HIDDEN_MAX_ERR = 3.0e-2 HIDDEN_MIN_RATIO = 2.0 def hidden_problem(prm): """A smooth, non-separable initial state vanishing on the boundary, with homogeneous Dirichlet data. It is NOT an eigenfunction of the operator, so no closed-form solution exists for the agent to substitute.""" def u0(X, Y): return (np.sin(np.pi * X) * np.sin(np.pi * Y) * (1.0 + prm["p1"] * X + prm["p2"] * Y + prm["p3"] * X * Y)) def g(X, Y, t): return np.zeros_like(X * Y * 1.0) return {"a": prm["a"], "b": prm["b"], "c": prm["c"], "u0": u0, "g": g} def grader_mcs(h, T, prob): """The grader's OWN MCS(theta=1/3) solver for arbitrary problem data. Independent of the submitted module; used only to produce a reference for the hidden instance. Any stable second-order ADI scheme converges to the same limit, so comparing against this reference stays scheme-neutral. """ from scipy.linalg import solve_banded th = 1.0 / 3.0 a_c, b_c, c_c = prob["a"], prob["b"], prob["c"] m = int(round(1.0 / h)); k = h; N = int(round(T / k)); h2 = h * h ax = np.linspace(0.0, 1.0, m + 1) X, Y = np.meshgrid(ax, ax, indexing="ij") def fill(V, t): G = np.array(prob["g"](X, Y, t), dtype=float, copy=True) G[1:m, 1:m] = V return G def Dxx(G): return (G[2:m + 1, 1:m] - 2 * G[1:m, 1:m] + G[0:m - 1, 1:m]) / h2 def Dyy(G): return (G[1:m, 2:m + 1] - 2 * G[1:m, 1:m] + G[1:m, 0:m - 1]) / h2 def Dxy(G): return (G[2:m + 1, 2:m + 1] - G[2:m + 1, 0:m - 1] - G[0:m - 1, 2:m + 1] + G[0:m - 1, 0:m - 1]) / (4 * h2) def F0(V, t): return 2 * c_c * Dxy(fill(V, t)) def F1(V, t): return a_c * Dxx(fill(V, t)) def F2(V, t): return b_c * Dyy(fill(V, t)) def F(V, t): G = fill(V, t) return 2 * c_c * Dxy(G) + a_c * Dxx(G) + b_c * Dyy(G) rx, ry = th * k * a_c / h2, th * k * b_c / h2 abx = np.zeros((3, m - 1)); abx[0, 1:] = -rx; abx[1, :] = 1 + 2 * rx; abx[2, :-1] = -rx aby = np.zeros((3, m - 1)); aby[0, 1:] = -ry; aby[1, :] = 1 + 2 * ry; aby[2, :-1] = -ry def sx(R, t): gb = np.asarray(prob["g"](X, Y, t), dtype=float); Rb = R.copy() Rb[0, :] += rx * gb[0, 1:m]; Rb[-1, :] += rx * gb[m, 1:m] return solve_banded((1, 1), abx, Rb) def sy(R, t): gb = np.asarray(prob["g"](X, Y, t), dtype=float); Rb = R.copy() Rb[:, 0] += ry * gb[1:m, 0]; Rb[:, -1] += ry * gb[1:m, m] return solve_banded((1, 1), aby, Rb.T).T U = np.asarray(prob["u0"](X, Y), dtype=float)[1:m, 1:m].copy() for n in range(N): tn, tn1 = n * k, (n + 1) * k F0n, F1n, F2n = F0(U, tn), F1(U, tn), F2(U, tn) Fn = F0n + F1n + F2n Y0 = U + k * Fn Y1 = sx(Y0 - th * k * F1n, tn1) Y2 = sy(Y1 - th * k * F2n, tn1) Yh0 = (Y0 + th * k * (F0(Y2, tn1) - F0n) + (0.5 - th) * k * (F(Y2, tn1) - Fn)) Yh1 = sx(Yh0 - th * k * F1n, tn1) U = sy(Yh1 - th * k * F2n, tn1) Ufull = np.array(prob["g"](X, Y, T), dtype=float, copy=True) Ufull[1:m, 1:m] = U return X, Y, Ufull def u_exact(x, y, t): s = np.zeros_like(x * y * 1.0) for (amp, al, be, ph) in MODES: lam = A_COEF * al * al + 2 * C_COEF * al * be + B_COEF * be * be s = s + amp * np.exp(-lam * t) * np.sin(al * x + be * y + ph) return s def fail(msg): print(f"FAIL: {msg}") print("SCORE: 0.0") sys.exit(0) def note(msg): print(f"NOTE: {msg}") def req_float(value, what): """Coerce a submission-supplied value to a finite float, or fail cleanly. Every numeric field in results.json is agent-controlled, so a bare float() can raise on a string or None and crash the grader with no verdict. A malformed submission must score 0.0, not produce a traceback. """ # Reject non-numeric JSON types outright. bool is checked first because it # is a subclass of int in Python -- float(True) == 1.0 would otherwise let # "a": true satisfy a numeric field. Strings are rejected too: "1" is not a # JSON number, and the published schema asks for numbers. if isinstance(value, bool) or not isinstance( value, (int, float, np.integer, np.floating)): fail(f"{what} must be a JSON number; got {value!r}") try: fv = float(value) except (TypeError, ValueError, OverflowError): # A JSON integer may be arbitrarily long; float() raises OverflowError # on one too large to represent. fail(f"{what} must be a finite number; got {type(value).__name__} " f"out of floating-point range") if not np.isfinite(fv): fail(f"{what} must be finite; got {value!r}") return fv def npz_array(arr, key, what=None): """Fetch one NPZ member safely and require a real numeric dtype. np.load(..., allow_pickle=False) raises on an object-typed member, which would otherwise escape as an uncaught ValueError and exit with no verdict. String, byte-string, boolean and complex members are rejected too: they would be silently coerced by np.asarray(..., float) and pass as valid data. """ what = what or f"reference_arrays.npz['{key}']" if key not in arr: fail(f"reference_arrays.npz missing array '{key}'") try: a = arr[key] except Exception as e: fail(f"{what} could not be read ({type(e).__name__}: {e}); object-typed " f"or pickled members are not accepted") a = np.asanyarray(a) if a.dtype == object or a.dtype.kind in ("U", "S", "b", "V", "O"): fail(f"{what} has dtype {a.dtype!r}; a real numeric array is required " f"(string, byte, boolean and object arrays are rejected)") if not np.issubdtype(a.dtype, np.number): fail(f"{what} has non-numeric dtype {a.dtype!r}") if np.iscomplexobj(a): fail(f"{what} is complex; a real array is required") return a def req_int(value, what): """Coerce to an EXACT integer, or fail cleanly (int() would truncate).""" fv = req_float(value, what) if fv != np.floor(fv): fail(f"{what} must be an exact integer; got {value!r}") return int(fv) # ============================================================================= # Subprocess runner: instruments linear solves, then regenerates the study. # ============================================================================= RUNNER = r''' import json, os, sys import numpy as np OUT_DIR, H_LIST_JSON, CAPTURE_DIR = sys.argv[1], sys.argv[2], sys.argv[3] HIDDEN_JSON = sys.argv[4] if len(sys.argv) > 4 else "null" POISON_SOLVES = json.loads(sys.argv[5]) if len(sys.argv) > 5 else False H_LIST = json.loads(H_LIST_JSON) HIDDEN = json.loads(HIDDEN_JSON) A_COEF, B_COEF, C_COEF = 1.0, 0.8, 0.5 MODES = [(1.0, 0.7, 1.3, 0.0), (0.7, 1.4, -0.6, 0.3), (0.5, 1.0, 0.8, 0.7)] def u_exact(x, y, t): s = np.zeros_like(x * y * 1.0) for (amp, al, be, ph) in MODES: lam = A_COEF * al * al + 2 * C_COEF * al * be + B_COEF * be * be s = s + amp * np.exp(-lam * t) * np.sin(al * x + be * y + ph) return s STATS = {"max_dim": 0, "n_solves": 0, "bandwidths": set(), "max_band": 0, "dense_calls": 0} def _record(dim, band=None): STATS["n_solves"] += 1 if int(dim) > STATS["max_dim"]: STATS["max_dim"] = int(dim) if band is not None: STATS["bandwidths"].add(int(band)) if int(band) > STATS["max_band"]: STATS["max_band"] = int(band) def _dense_bandwidth(a): """Half-bandwidth of a dense square matrix (0 = diagonal, 1 = tridiagonal). A dense call whose matrix is not tridiagonal is not a 1-D sweep, whatever its dimension -- this closes the 'small dense solve' gap. """ a = np.asarray(a) if a.ndim != 2 or a.shape[0] != a.shape[1]: return None nz = np.nonzero(np.abs(a) > 0) if nz[0].size == 0: return 0 return int(np.max(np.abs(nz[0] - nz[1]))) def _patch(): """Wrap every linear-solve entry point a submission could plausibly use. Patching happens BEFORE adi_solver is imported, so `from scipy.linalg import solve_banded`-style bindings pick up the wrapper. """ import numpy.linalg as nla import scipy.linalg as sla def wrap(mod, name, probe): """probe(*a, **kw) -> (dimension, half_bandwidth_or_None).""" orig = getattr(mod, name, None) if orig is None: return def poison(value): """Replace numerical solve results by NaNs for the causal audit.""" if callable(value): def poisoned_callable(*a, **kw): return poison(value(*a, **kw)) return poisoned_callable if isinstance(value, tuple): return tuple(poison(v) for v in value) if isinstance(value, list): return [poison(v) for v in value] if isinstance(value, np.ndarray) and value.dtype.kind in "fc": out = np.array(value, copy=True) out[...] = np.nan return out return value def wrapper(*a, **kw): try: dim, band = probe(*a, **kw) _record(dim, band) except Exception: pass result = orig(*a, **kw) return poison(result) if POISON_SOLVES else result wrapper.__name__ = name setattr(mod, name, wrapper) def _banded(l_and_u, ab, *a, **kw): lu = tuple(int(v) for v in l_and_u) return np.asarray(ab).shape[-1], max(lu) def _dense(a, *r, **kw): STATS["dense_calls"] += 1 return np.asarray(a).shape[-1], _dense_bandwidth(a) # banded / tridiagonal entry points -- the legitimate ADI route wrap(sla, "solve_banded", _banded) wrap(sla, "solveh_banded", lambda ab, *a, **kw: (np.asarray(ab).shape[-1], np.asarray(ab).shape[0] - 1)) wrap(sla, "cholesky_banded", lambda ab, *a, **kw: (np.asarray(ab).shape[-1], np.asarray(ab).shape[0] - 1)) wrap(sla, "cho_solve_banded", lambda cb_lower, *a, **kw: (np.asarray(cb_lower[0]).shape[-1], np.asarray(cb_lower[0]).shape[0] - 1)) # dense entry points -- allowed only if the matrix is actually tridiagonal wrap(sla, "solve", _dense) wrap(sla, "lu_factor", _dense) wrap(sla, "solve_triangular", _dense) wrap(sla, "lu_solve", lambda lu_piv, *r, **kw: (np.asarray(lu_piv[0]).shape[-1], None)) wrap(nla, "solve", _dense) wrap(nla, "inv", _dense) wrap(nla, "lstsq", _dense) # low-level LAPACK entry points (evasion route named by the audit) try: import scipy.linalg.lapack as lap for nm in ("dgtsv", "sgtsv", "zgtsv", "cgtsv", # tridiagonal "dptsv", "sptsv"): # sym. tridiagonal wrap(lap, nm, lambda dl, *r, **kw: (np.asarray(dl).size + 1, 1)) for nm in ("dgesv", "sgesv", "zgesv", "cgesv", # dense general "dgetrf", "sgetrf", "dposv", "dpbsv"): wrap(lap, nm, _dense) for nm in ("dgbsv", "sgbsv"): # banded general wrap(lap, nm, lambda kl, ku, ab, *r, **kw: (np.asarray(ab).shape[-1], max(int(kl), int(ku)))) except Exception: pass try: import scipy.sparse.linalg as ssl def _sparse(A, *r, **kw): n = int(np.asarray(A.shape)[-1]) try: # a sparse tridiagonal sweep is legitimate coo = A.tocoo() band = int(np.max(np.abs(coo.row - coo.col))) if coo.nnz else 0 except Exception: band = None return n, band for nm in ("spsolve", "splu", "factorized", "cg", "gmres", "bicgstab", "spsolve_triangular", "minres", "lgmres", "use_solver"): wrap(ssl, nm, _sparse) except Exception: pass _patch() sys.path.insert(0, OUT_DIR) try: import adi_solver except Exception as e: print(json.dumps({"error": f"import failed: {type(e).__name__}: {e}"})) sys.exit(0) if not hasattr(adi_solver, "solve"): print(json.dumps({"error": "adi_solver.py does not define " "solve(h, T=1.0, problem=None)"})) sys.exit(0) results = [] for run_index, h in enumerate(H_LIST): STATS["max_dim"] = 0 STATS["n_solves"] = 0 STATS["max_band"] = 0 STATS["dense_calls"] = 0 m = int(round(1.0 / h)) try: X, Y, U = adi_solver.solve(h, 1.0) except TypeError: try: X, Y, U = adi_solver.solve(h) except Exception as e: print(json.dumps({"error": f"solve({h}) raised {type(e).__name__}: {e}"})) sys.exit(0) except Exception as e: print(json.dumps({"error": f"solve({h}) raised {type(e).__name__}: {e}"})) sys.exit(0) X, Y, U = np.asarray(X, float), np.asarray(Y, float), np.asarray(U, float) if U.shape != (m + 1, m + 1) or X.shape != U.shape or Y.shape != U.shape: print(json.dumps({"error": f"solve({h}) returned shapes " f"{X.shape}/{Y.shape}/{U.shape}, expected " f"{(m + 1, m + 1)} for all three"})) sys.exit(0) ax = np.linspace(0.0, 1.0, m + 1) gx, gy = np.meshgrid(ax, ax, indexing="ij") grid_ok = bool(np.allclose(X, gx, rtol=0.0, atol=1e-12) and np.allclose(Y, gy, rtol=0.0, atol=1e-12)) err = float(np.max(np.abs(U - u_exact(X, Y, 1.0)))) # Capture the regenerated field outside the submitted output directory. # The parent grader loads it before this private temporary directory is # removed and compares the finest field pointwise with U_finest. solution_file = f"regenerated_U_{run_index}.npy" np.save(os.path.join(CAPTURE_DIR, solution_file), U, allow_pickle=False) results.append({ "h": h, "error": err, "finite": bool(np.all(np.isfinite(U))), "grid_ok": grid_ok, "max_dim": STATS["max_dim"], "n_solves": STATS["n_solves"], "max_band": STATS["max_band"], "dense_calls": STATS["dense_calls"], "dim_limit": 2 * (m + 1), "solution_file": solution_file, }) hidden_out = [] if HIDDEN is not None: prm = HIDDEN def _u0(X, Y): return (np.sin(np.pi * X) * np.sin(np.pi * Y) * (1.0 + prm["p1"] * X + prm["p2"] * Y + prm["p3"] * X * Y)) def _g(X, Y, t): return np.zeros_like(X * Y * 1.0) prob = {"a": prm["a"], "b": prm["b"], "c": prm["c"], "u0": _u0, "g": _g} for h in prm["h_list"]: STATS["max_dim"] = 0; STATS["n_solves"] = 0 STATS["max_band"] = 0; STATS["dense_calls"] = 0 m = int(round(1.0 / h)) try: X, Y, U = adi_solver.solve(h, 1.0, prob) except TypeError as e: print(json.dumps({"error": f"solve() does not accept a `problem` " f"argument ({e}); the specification " f"requires solve(h, T=1.0, problem=None)"})) sys.exit(0) except Exception as e: print(json.dumps({"error": f"solve({h}, problem=hidden) raised " f"{type(e).__name__}: {e}"})) sys.exit(0) U = np.asarray(U, float) if U.shape != (m + 1, m + 1): print(json.dumps({"error": f"solve({h}, problem=hidden) returned " f"shape {U.shape}, expected {(m+1, m+1)}"})) sys.exit(0) np.save(os.path.join(CAPTURE_DIR, f"hidden_{m}.npy"), U, allow_pickle=False) hidden_out.append({"h": h, "m": m, "finite": bool(np.all(np.isfinite(U))), "n_solves": STATS["n_solves"], "max_dim": STATS["max_dim"], "max_band": STATS["max_band"], "dim_limit": 2 * (m + 1)}) print(json.dumps({ "results": results, "hidden": hidden_out, "bandwidths": sorted(STATS["bandwidths"]), })) ''' def run_solver(out_dir, h_list, budget, hidden=None, poison_solves=False): """Execute the submitted solver in an isolated subprocess. Returns (payload_dict, error_string). Exactly one is None. """ with tempfile.TemporaryDirectory() as td: runner_path = os.path.join(td, "_repro_runner.py") with open(runner_path, "w") as fh: fh.write(RUNNER) env = dict(os.environ) env["PYTHONDONTWRITEBYTECODE"] = "1" try: proc = subprocess.run( [sys.executable, runner_path, os.path.abspath(out_dir), json.dumps(list(h_list)), td, json.dumps(hidden), json.dumps(bool(poison_solves))], capture_output=True, text=True, timeout=budget, env=env, cwd=td, ) except subprocess.TimeoutExpired: return None, f"timed out after {budget:.0f}s" if proc.returncode != 0: tail = (proc.stderr or "").strip().splitlines()[-3:] return None, f"runner exited {proc.returncode}: {' | '.join(tail)}" out = (proc.stdout or "").strip().splitlines() if not out: return None, "runner produced no output" try: payload = json.loads(out[-1]) except Exception as e: return None, f"unparseable runner output ({e})" if "error" in payload: return None, payload["error"] try: payload["solutions"] = [ np.load(os.path.join(td, row.pop("solution_file")), allow_pickle=False) for row in payload["results"] ] payload["hidden_fields"] = [ np.load(os.path.join(td, f"hidden_{row['m']}.npy"), allow_pickle=False) for row in payload.get("hidden", []) ] except Exception as e: return None, f"could not load regenerated solution capture ({e})" return payload, None def check_imports(solver_path): """Static allowlist over top-level imports. Returns a manual-review flag.""" try: with open(solver_path) as fh: src = fh.read() except Exception as e: fail(f"adi_solver.py unreadable ({e})") try: tree = ast.parse(src) except SyntaxError as e: fail(f"adi_solver.py does not parse ({e})") offenders, dynamic = set(), False for node in ast.walk(tree): if isinstance(node, ast.Import): for a in node.names: top = a.name.split(".")[0] if top not in ALLOWED_TOP_LEVEL: offenders.add(top) elif isinstance(node, ast.ImportFrom): if node.level == 0 and node.module: top = node.module.split(".")[0] if top not in ALLOWED_TOP_LEVEL: offenders.add(top) elif isinstance(node, ast.Name) and node.id == "__import__": dynamic = True elif isinstance(node, ast.Attribute) and node.attr in ( "import_module", "exec_module", "spec_from_file_location"): dynamic = True if offenders: fail("adi_solver.py imports disallowed module(s): " + ", ".join(sorted(offenders)) + " (Python + NumPy + SciPy only)") return dynamic # ============================================================================= def main(out_dir, budget, do_exec, strict=False): solver_path = os.path.join(out_dir, "adi_solver.py") if not os.path.isfile(solver_path): fail("adi_solver.py missing (the solver is a required deliverable)") try: with open(os.path.join(out_dir, "results.json")) as fh: res = json.load(fh) except Exception as e: fail(f"results.json missing or invalid JSON ({e})") if not isinstance(res, dict): fail("results.json top level must be a JSON object") try: arr = np.load(os.path.join(out_dir, "reference_arrays.npz")) except Exception as e: fail(f"reference_arrays.npz missing or unreadable ({e})") # ---------------------------------------------------------------- schema # Tier 1/Tier 2 arrays are mandatory. The Tier 3 arrays are NOT: the rubric # promises 0.5 for an honest Tier 3 omission, so their absence must degrade # to a Tier 3 failure, never to a schema failure. scalar_keys = {"h_values": (5,), "max_norm_errors": (5,), "eoc": (5,)} grid_keys = ("X_finest", "Y_finest", "U_finest", "U_exact_finest", "error_finest") for key in scalar_keys: if key not in arr: fail(f"reference_arrays.npz missing required array '{key}'") for key, shape in scalar_keys.items(): got = npz_array(arr, key).shape if got != shape: fail(f"'{key}' has shape {got}, expected {shape}") missing_tier3 = [k for k in grid_keys if k not in arr] h_values = np.asarray(npz_array(arr, "h_values"), float) errors = np.asarray(npz_array(arr, "max_norm_errors"), float) eoc = np.asarray(npz_array(arr, "eoc"), float) if not np.all(np.isfinite(errors)): fail("non-finite entries in max_norm_errors (NaN/Inf) -> unstable scheme") if not np.all(errors > 0): fail("max_norm_errors must be strictly positive") if not np.all(np.isfinite(h_values)): fail("non-finite entries in h_values") # the schema requires descending storage order -- enforce it as stored, # rather than sorting first and losing the requirement if not np.allclose(h_values, H_VALUES, rtol=1e-12, atol=0.0): if np.allclose(np.sort(h_values)[::-1], H_VALUES, rtol=1e-12, atol=0.0): fail("h_values must be stored in DESCENDING order 0.1*2^-p, p=0..4 " f"(got {np.array2string(h_values, precision=5)})") fail("h_values do not match the required mesh sequence 0.1*2^-p, p=0..4") # eoc: entry 0 has no predecessor; entries 1..4 recomputed from the errors. # NaN specifically -- +/-Inf is non-finite too but is not what the schema # asks for, and would otherwise slip past a mere isfinite() test. if not np.isnan(eoc[0]): fail(f"eoc[0] must be NaN (the coarsest mesh has no predecessor); " f"got {eoc[0]!r}") if not np.all(np.isfinite(eoc[1:])): fail(f"eoc[1:] must all be finite; got " f"{np.array2string(eoc[1:], precision=4)}") recomputed = np.log2(errors[:-1] / errors[1:]) if not np.allclose(eoc[1:], recomputed, rtol=1e-6, atol=1e-12): fail(f"eoc inconsistent with max_norm_errors: submitted " f"{np.array2string(eoc[1:], precision=4)} vs recomputed " f"{np.array2string(recomputed, precision=4)}") tier3_ok, reason = True, "" if missing_tier3: tier3_ok = False reason = f"missing Tier 3 array(s): {', '.join(missing_tier3)}" else: for key in grid_keys: a = np.asarray(npz_array(arr, key), float) if a.shape != (FINEST_POINTS, FINEST_POINTS): tier3_ok, reason = False, ( f"'{key}' has shape {a.shape}, expected " f"({FINEST_POINTS}, {FINEST_POINTS})") break if not np.all(np.isfinite(a)): tier3_ok, reason = False, f"'{key}' contains non-finite entries" break n_present = 3 + len(grid_keys) - len(missing_tier3) print(f"PASS schema: {n_present}/8 arrays present, correctly shaped, finite; " f"mesh sequence (descending) and eoc consistent") # --------------------------------------------------- results.json sidecar table = res.get("convergence_table") if not isinstance(table, list) or len(table) != 5: fail("results.json convergence_table must be a list of 5 rows") for i, r in enumerate(table): if not isinstance(r, dict): fail(f"results.json convergence_table row {i} is not an object") for key in ("p", "h", "max_norm_error"): if key not in r: fail(f"results.json convergence_table row {i} is missing '{key}'") if "eoc_order" not in r: fail(f"results.json convergence_table row {i} is missing 'eoc_order'") rows = sorted(table, key=lambda r: -req_float(r["h"], "convergence_table 'h'")) j_h = np.array([req_float(r["h"], "convergence_table 'h'") for r in rows]) j_err = np.array([req_float(r["max_norm_error"], "convergence_table 'max_norm_error'") for r in rows]) j_eoc = [r["eoc_order"] for r in rows] # Validate the RAW p before coercion -- int() truncates, so 0.25 would # otherwise pass as 0, and a string would crash a bare float(). j_p = [req_int(r["p"], "convergence_table 'p'") for r in rows] if not np.allclose(j_h, H_VALUES, rtol=1e-12, atol=0.0): fail("results.json convergence_table h values do not match the mesh sequence") if j_p != list(range(5)): fail(f"results.json convergence_table 'p' values are {j_p}; each row must " f"carry p = 0..4 matching its h = 0.1*2^-p") if not np.allclose(j_err, errors, rtol=1e-6, atol=1e-18): fail("results.json max_norm_error values inconsistent with reference_arrays.npz") if j_eoc[0] is not None: fail(f"results.json eoc_order for the coarsest mesh must be null; " f"got {j_eoc[0]!r}") j_eoc_tail = np.array([req_float(v, "convergence_table 'eoc_order'") for v in j_eoc[1:]]) if not np.allclose(j_eoc_tail, recomputed, rtol=1e-6, atol=1e-12): fail("results.json eoc_order values inconsistent with the submitted errors") dm = res.get("diffusion_matrix") if not isinstance(dm, dict): fail("results.json missing diffusion_matrix") for name, want in (("a", A_COEF), ("b", B_COEF), ("c", C_COEF)): if name not in dm: fail(f"results.json diffusion_matrix is missing '{name}'") got = req_float(dm[name], f"diffusion_matrix['{name}']") if not np.isclose(got, want, rtol=1e-12): fail(f"results.json diffusion_matrix['{name}'] = {dm[name]!r}, " f"expected {want}") if "final_time" not in res: fail("results.json missing final_time") if not np.isclose(req_float(res["final_time"], "final_time"), T_FINAL, rtol=1e-12): fail(f"results.json final_time must be {T_FINAL}") tsr = str(res.get("time_step_rule", "")) if not re.search(r"\bk\b\s*=\s*\bh\b", tsr): fail(f"results.json time_step_rule = {tsr!r} must actually state 'k = h'") if not isinstance(res.get("scheme"), str) or not res["scheme"].strip(): fail("results.json 'scheme' must be a non-empty string naming the " "scheme used") # Optional extras: not required by problem_spec.md, but if a submission # reports them they must be sane rather than arbitrary. 'theta' is # scheme-specific so only its type is constrained; 'asymptotic_order' is a # claim about this PDE and must be 2. if "asymptotic_order" in res: if req_int(res["asymptotic_order"], "results.json asymptotic_order") != 2: fail("results.json asymptotic_order must be 2 (the scheme is " "second order in h with k = h)") if "theta" in res: req_float(res["theta"], "results.json theta") fm = res.get("finest_mesh") if not isinstance(fm, dict): fail("results.json missing finest_mesh") for key in ("h", "grid_points_per_side", "max_norm_error"): if key not in fm: fail(f"results.json finest_mesh is missing '{key}'") if not np.isclose(req_float(fm["h"], "finest_mesh.h"), H_VALUES[-1], rtol=1e-12): fail(f"results.json finest_mesh.h must be {H_VALUES[-1]}") if req_int(fm["grid_points_per_side"], "finest_mesh.grid_points_per_side") != FINEST_POINTS: fail(f"results.json finest_mesh.grid_points_per_side must be {FINEST_POINTS}") if not np.isclose(req_float(fm["max_norm_error"], "finest_mesh.max_norm_error"), errors[-1], rtol=1e-6): fail("results.json finest_mesh.max_norm_error inconsistent with the arrays") ls_order = float(-np.polyfit(np.arange(5), np.log2(errors), 1)[0]) if "observed_order_of_convergence" not in res: fail("results.json missing observed_order_of_convergence") ooc = req_float(res["observed_order_of_convergence"], "observed_order_of_convergence") if abs(ooc - ls_order) > 0.25: fail(f"results.json observed_order_of_convergence = {ooc:.3f} is not " f"supported by the submitted errors (fitted order {ls_order:.3f})") print("PASS sidecar: results.json fully consistent with the arrays " "(table, orders, diffusion matrix, finest mesh, metadata)") # ============================ Tier 1 ===================================== if not (errors[0] < TIER1_MAX_ERR): fail(f"Tier 1: e_h(0.1) = {errors[0]:.3e} >= {TIER1_MAX_ERR:.1e} " f"(no usable solution on the coarsest mesh)") print(f"PASS Tier 1: e_h(0.1) = {errors[0]:.6e}, finite and stable") # ============================ Tier 2 ===================================== if not np.all(np.diff(errors) < 0): fail(f"Tier 2: errors are not strictly decreasing under refinement " f"({np.array2string(errors, precision=3)})") reduction = errors[0] / errors[-1] if reduction < MIN_TOTAL_REDUCTION: fail(f"Tier 2: total error reduction over the 16x refinement is only " f"{reduction:.1f}x (need >= {MIN_TOTAL_REDUCTION:.0f}x for second " f"order; a dropped or first-order-lagged mixed term lands here)") if not (LS_ORDER_LO <= ls_order <= LS_ORDER_HI): fail(f"Tier 2: fitted order of convergence {ls_order:.3f} outside " f"[{LS_ORDER_LO}, {LS_ORDER_HI}]") if not all(PT_ORDER_LO <= v <= PT_ORDER_HI for v in recomputed[-2:]): fail(f"Tier 2: the two finest experimental orders " f"{np.array2string(recomputed[-2:], precision=3)} are not both in " f"[{PT_ORDER_LO}, {PT_ORDER_HI}]") if not (errors[-1] < FINEST_MAX_ERR): fail(f"Tier 2: finest-mesh error {errors[-1]:.3e} >= {FINEST_MAX_ERR:.1e}") print(f"PASS Tier 2: orders = {np.array2string(recomputed, precision=3)}, " f"fitted {ls_order:.3f}, total reduction {reduction:.1f}x, " f"finest {errors[-1]:.3e}") # ======================= Provenance / ADI enforcement ==================== review = [] # machine-readable manual-review reasons def flag(reason): review.append(reason) print(f"REVIEW_REQUIRED: {reason}") dynamic_imports = check_imports(solver_path) print("PASS dependencies: adi_solver.py imports only stdlib + NumPy + SciPy") if dynamic_imports: flag("adi_solver.py performs dynamic imports; the static dependency " "check is not exhaustive") fine_verified = False if not do_exec: flag("--no-exec given: solver execution, reproducibility and the " "sweep audit were SKIPPED; provenance is unverified") else: hidden_spec = dict(HIDDEN); hidden_spec["h_list"] = list(HIDDEN_H) payload, err = run_solver(out_dir, COARSE_H, budget, hidden=hidden_spec) if err is not None: fail(f"could not regenerate results from adi_solver.py: {err}") rows = payload["results"] for r in rows: if not r["finite"]: fail(f"regenerated solution at h={r['h']:.5f} is non-finite") if not r["grid_ok"]: fail(f"solve({r['h']:.5f}) returned coordinate arrays that are not " f"the uniform grid with indexing='ij'") regen = np.array([r["error"] for r in rows]) want = errors[:len(regen)] if not np.allclose(regen, want, rtol=REPRO_RTOL, atol=0.0): fail("adi_solver.py does not reproduce the submitted errors: " f"regenerated {np.array2string(regen, precision=6)} vs submitted " f"{np.array2string(want, precision=6)}") print(f"PASS provenance: adi_solver.py regenerated e_h at " f"h = {', '.join(f'{h:g}' for h in COARSE_H)} and matched the " f"submitted values to {REPRO_RTOL:.0e} relative") unobserved = [] # meshes where no library solve was seen # ---- the hidden instance: does the solver actually solve the PDE? ---- # Everything above proves the submitted numbers came from this code. It # does NOT prove the code solves anything: a solver can make one legal # tridiagonal sweep, discard it, and return the manufactured analytic # solution plus a synthetic O(h^2) error. So we hand it a problem whose # solution it cannot write down and compare against our own reference. hid = payload.get("hidden") or [] if len(hid) != len(HIDDEN_H): fail("adi_solver.py did not return solutions for the hidden instance; " "solve(h, T, problem) must accept and solve an arbitrary problem") prob = hidden_problem(HIDDEN) _, _, Uref = grader_mcs(HIDDEN_H_REF, T_FINAL, prob) hid_err = [] for row, Uag in zip(hid, payload["hidden_fields"]): if not row["finite"]: fail(f"hidden instance at h={row['h']:.5f}: non-finite solution") step = int(round(row["h"] / HIDDEN_H_REF)) e = float(np.max(np.abs(np.asarray(Uag, float) - Uref[::step, ::step]))) hid_err.append(e) if hid_err[-1] >= HIDDEN_MAX_ERR: fail(f"the solver does not solve the hidden instance: at " f"h={HIDDEN_H[-1]:.5f} it differs from an independently computed " f"reference by {hid_err[-1]:.3e} (limit {HIDDEN_MAX_ERR:.1e}). " f"Reproducing the benchmark's own numbers is not enough -- the " f"solver must actually integrate the PDE for data it has not " f"seen. A field returned analytically, or copied from the " f"manufactured solution, fails here.") ratio = hid_err[0] / hid_err[-1] if hid_err[-1] > 0 else float("inf") if ratio < HIDDEN_MIN_RATIO: fail(f"the solver does not converge on the hidden instance: errors " f"{np.array2string(np.array(hid_err), precision=3)} give a " f"refinement ratio of {ratio:.2f} (need >= {HIDDEN_MIN_RATIO}). " f"An h-independent or analytically constructed field lands here.") print(f"PASS hidden instance: solve() integrates an undisclosed problem " f"(a={HIDDEN['a']}, b={HIDDEN['b']}, c={HIDDEN['c']}, non-separable " f"initial state) to within {hid_err[-1]:.3e} of an independently " f"computed reference, converging by {ratio:.2f}x under refinement") # Accuracy on general data rules out a solver specialised to the public # manufactured solution, but it does not by itself prove that the # observed tridiagonal calls produced the field: an explicit substepper # plus a discarded dummy solve could still pass. Re-run one hidden mesh # with every intercepted numerical solve result replaced by NaNs. A # genuine ADI data path must then change, become non-finite, or raise. # Hand-written Thomas solvers have no interceptable call and remain a # documented manual-review case handled by audit_sweeps below. observed_index = next((i for i, row in enumerate(hid) if row["n_solves"]), None) if observed_index is not None: poison_h = hid[observed_index]["h"] poison_spec = dict(HIDDEN); poison_spec["h_list"] = [poison_h] poisoned, poison_err = run_solver( out_dir, [], budget, hidden=poison_spec, poison_solves=True) if poison_err is not None: # The normal run above succeeded. An exception specifically # while solving the poisoned hidden instance therefore proves # that an intercepted solve result entered the computational # path. Import/runner failures are inconclusive, not evidence. if "problem=hidden) raised" in poison_err: print("PASS causal sweep audit: poisoning intercepted solve " "results caused the hidden solve to raise") else: flag("causal sweep audit could not complete after solve-result " f"poisoning ({poison_err})") else: poison_rows = poisoned.get("hidden") or [] poison_fields = poisoned.get("hidden_fields") or [] if len(poison_rows) != 1 or len(poison_fields) != 1: fail("causal sweep audit did not return exactly one hidden field") changed = (not poison_rows[0]["finite"] or not np.allclose(poison_fields[0], payload["hidden_fields"][observed_index], rtol=0.0, atol=1e-12, equal_nan=True)) if not changed: fail("ADI authenticity: poisoning every intercepted linear-" "solve result left the hidden solution unchanged. The " "observed tridiagonal calls do not influence the returned " "field (discarded/dummy solves are not an ADI method).") print("PASS causal sweep audit: poisoning intercepted solve " "results materially changed the hidden solution") def audit_sweeps(rs, label): """Dimension AND bandwidth audit over one batch of regenerated runs. Every mesh with zero observable solves is recorded, not just a batch that is entirely unobserved: otherwise a hybrid solver could make legal instrumented calls on the coarse meshes and switch to an uninstrumented path on the fine ones, drawing no flag at all. """ for r in rs: if not r["n_solves"]: unobserved.append(r["h"]) continue m = int(round(1.0 / r["h"])) if r["max_dim"] > r["dim_limit"]: fail(f"non-ADI linear solve ({label}): at h={r['h']:.5f} (m={m}) " f"a system of dimension {r['max_dim']} was solved, exceeding " f"the 1-D sweep limit {r['dim_limit']}. Only tridiagonal " f"one-dimensional sweeps may be solved implicitly; the full " f"2-D operator must not be assembled and inverted.") if r["max_band"] > MAX_BANDWIDTH: fail(f"non-tridiagonal implicit solve ({label}): at " f"h={r['h']:.5f} a system of half-bandwidth " f"{r['max_band']} was solved (dimension {r['max_dim']}). " f"The specification requires tridiagonal 1-D sweeps.") audit_sweeps(rows, "coarse meshes") audit_sweeps(hid, "hidden instance") total_solves = sum(r["n_solves"] for r in rows) if total_solves: # report the mesh where the solve came closest to the per-mesh limit tight = max((r for r in rows if r["n_solves"]), key=lambda r: r["max_dim"] / r["dim_limit"]) worst_bw = max(r["max_band"] for r in rows) dense = sum(r["dense_calls"] for r in rows) print(f"PASS ADI structure: {total_solves} linear solves observed on " f"the coarse meshes; tightest was dimension {tight['max_dim']} " f"against a limit of {tight['dim_limit']} at h={tight['h']:g} " f"(monolithic would be {(int(round(1/tight['h'])) - 1)**2}); " f"max half-bandwidth {worst_bw} (tridiagonal)" + (f"; {dense} dense calls, all tridiagonal" if dense else "")) # Tier 3 credit requires the finest meshes to be regenerated too -- # unverified artifacts do not earn full marks. payload2, err2 = run_solver(out_dir, FINE_H, budget) if err2 is not None: flag(f"finest-mesh regeneration did not complete ({err2}); " f"max_norm_errors[3:] and the Tier 3 surfaces are NOT proven to " f"come from adi_solver.py") else: regen2 = np.array([r["error"] for r in payload2["results"]]) for r in payload2["results"]: if not r["finite"]: fail(f"regenerated solution at h={r['h']:.5f} is non-finite") if not r["grid_ok"]: fail(f"solve({r['h']:.5f}) returned coordinate arrays that are " f"not the uniform grid with indexing='ij'") if not np.allclose(regen2, errors[3:], rtol=REPRO_RTOL, atol=0.0): fail("adi_solver.py does not reproduce the submitted errors on the " f"finest meshes: regenerated {np.array2string(regen2, precision=6)} " f"vs submitted {np.array2string(errors[3:], precision=6)}") if tier3_ok: regenerated_finest = np.asarray(payload2["solutions"][-1], float) submitted_finest = np.asarray(npz_array(arr, "U_finest"), float) if not np.allclose(regenerated_finest, submitted_finest, rtol=SURFACE_RTOL, atol=SURFACE_ATOL): max_delta = float(np.max(np.abs( regenerated_finest - submitted_finest))) tier3_ok, reason = False, ( "U_finest does not match the pointwise field regenerated " "from adi_solver.py " f"(max difference {max_delta:.3e}; tolerances " f"rtol={SURFACE_RTOL:.0e}, atol={SURFACE_ATOL:.0e})") else: print("PASS Tier 3 provenance: U_finest matches the pointwise " "field regenerated from adi_solver.py") audit_sweeps(payload2["results"], "finest meshes") fine_verified = True print("PASS provenance: all five meshes regenerated from adi_solver.py") # One flag covering every mesh whose solves were not observable. A # hand-written Thomas algorithm is permitted, so this is a review flag # rather than a failure -- but --strict refuses to certify it, and a # hybrid solver that hides only its fine-mesh path is now caught too. if unobserved: hs = ", ".join(f"h={v:g}" for v in sorted(unobserved, reverse=True)) flag(f"no library linear solves were observed on {len(unobserved)} of " f"the regenerated meshes ({hs}). A hand-written Thomas algorithm " f"is permitted by the specification, so this is not by itself a " f"violation -- but the sweep audit is INCONCLUSIVE for those " f"meshes and the ADI structure must be confirmed by reading " f"adi_solver.py. Note that observable solves on some meshes and " f"none on others is exactly what a hybrid evasion looks like.") # ============================ Tier 3 ===================================== # Full credit requires the Tier 3 surfaces to be BOTH analytically correct # and provably generated by the submitted solver. if tier3_ok and do_exec and not fine_verified: tier3_ok, reason = False, ( "the finest-mesh results were not regenerated from adi_solver.py, so " "the Tier 3 surfaces are unverified (raise --budget to allow a slower " "solver to complete)") if not tier3_ok: print(f"PARTIAL Tier 3 failed: {reason}") return finish(review, strict, 0.5, do_exec) X = np.asarray(npz_array(arr, "X_finest"), float) Y = np.asarray(npz_array(arr, "Y_finest"), float) U = np.asarray(npz_array(arr, "U_finest"), float) Ue_sub = np.asarray(npz_array(arr, "U_exact_finest"), float) Esurf = np.asarray(npz_array(arr, "error_finest"), float) ax = np.linspace(0.0, 1.0, FINEST_POINTS) gx, gy = np.meshgrid(ax, ax, indexing="ij") Ue = u_exact(gx, gy, T_FINAL) recomputed_err = float(np.max(np.abs(U - Ue))) # rtol=0.0 is essential: NumPy's default rtol=1e-5 would otherwise accept a # coordinate shifted by ~1e-6 on an O(1) grid, which is not 'exact'. if (not np.allclose(X, gx, rtol=0.0, atol=1e-12) or not np.allclose(Y, gy, rtol=0.0, atol=1e-12)): tier3_ok, reason = False, ("X_finest/Y_finest are not the uniform 161x161 " "grid on [0,1]^2 with indexing='ij'") elif not np.allclose(Ue_sub, Ue, rtol=1e-10, atol=1e-12): tier3_ok, reason = False, ("U_exact_finest does not match the analytic exact " "solution at T = 1") elif recomputed_err >= FINEST_MAX_ERR: tier3_ok, reason = False, (f"finest max-norm error {recomputed_err:.3e} " f">= {FINEST_MAX_ERR:.1e}") elif not np.allclose(Esurf, U - Ue, rtol=0.0, atol=1e-12): tier3_ok, reason = False, "error_finest != U_finest - U_exact_finest" elif not np.isclose(recomputed_err, errors[-1], rtol=1e-6): tier3_ok, reason = False, ( f"max_norm_errors[4] = {errors[-1]:.6e} but the error recomputed from " f"U_finest is {recomputed_err:.6e}: the convergence table and the " f"Tier 3 surfaces did not come from the same run") else: bnd = max(np.max(np.abs(U[0, :] - Ue[0, :])), np.max(np.abs(U[-1, :] - Ue[-1, :])), np.max(np.abs(U[:, 0] - Ue[:, 0])), np.max(np.abs(U[:, -1] - Ue[:, -1]))) if bnd > 1e-12: tier3_ok, reason = False, (f"U_finest does not carry the exact Dirichlet " f"data on the boundary (max deviation {bnd:.3e})") if tier3_ok: print(f"PASS Tier 3: 161x161 grid, exact solution, error surface and " f"boundary data all consistent; max-norm {recomputed_err:.6e}") else: print(f"PARTIAL Tier 3 failed: {reason}") return finish(review, strict, 1.0 if tier3_ok else 0.5, do_exec) def finish(review, strict, score, do_exec=True): """Emit the review summary, then exactly one verdict line. The score is resolved BEFORE anything is printed, so a parser that takes the first (or the last) 'SCORE:' match reads the same value either way. Without solver execution nothing about provenance has been established, so no official score is produced at all -- a diagnostic run must never be mistakable for a graded one. Exit code 2 distinguishes 'ungraded' from a graded 0.0. """ if review and strict: score = 0.0 if not review: print("REVIEW: none required - every check was conclusive") else: print(f"REVIEW: {len(review)} item(s) need human confirmation:") for r in review: print(f" - {r}") if strict: print("STRICT: an inconclusive check cannot be certified " "automatically -> score forced to 0.0") if not do_exec: print("RESULT: UNGRADED - provenance was not verified (--no-exec). " "No official score is produced; re-run without --no-exec to grade.") sys.exit(2) print(f"SCORE: {score}") if __name__ == "__main__": ap = argparse.ArgumentParser(description="Grade a mixed-derivative ADI submission.") ap.add_argument("output_dir") ap.add_argument("--budget", type=float, default=DEFAULT_BUDGET, help="seconds allowed per solver-regeneration phase") ap.add_argument("--no-exec", action="store_true", help="skip executing the submitted solver (diagnostics only; " "provenance is then unverified)") ap.add_argument("--strict", action="store_true", help="treat any REVIEW_REQUIRED flag as a failure. Use when no " "human will read the manual-review items; note this also " "rejects the spec-permitted hand-written Thomas algorithm.") a = ap.parse_args() main(a.output_dir, a.budget, not a.no_exec, a.strict)