# Evaluation — Mixed-Derivative ADI Benchmark

Run:

```bash
python grade.py <agent_output_dir>                 # full grading
python grade.py <agent_output_dir> --budget 900    # allow a slower solver
python grade.py <agent_output_dir> --strict        # refuse to certify anything inconclusive
python grade.py <agent_output_dir> --no-exec       # diagnostics only -> UNGRADED, no score
```

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**. It does not trust `results.json`, and it does not trust
`reference_arrays.npz` either — both must agree with each other *and* with what
`adi_solver.py` actually produces when executed. In particular, the regenerated
finest numerical field must match `U_finest` pointwise; matching only its scalar
maximum error is insufficient.

**Every acceptance threshold below is published verbatim in `problem_spec.md`.**
The agent is told the bar in advance, so no scheme is privileged by a hidden band.

## Score behavior

| Score | Condition |
|-------|-----------|
| **1.0** | Tier 1, Tier 2, and Tier 3 all pass — including regeneration of **all five** meshes from the submitted solver. |
| **0.5** | Tier 1, Tier 2 and the coarse-mesh provenance checks pass, but Tier 3 is missing, fails, or is unverified. |
| **0.0** | Required files missing (including `adi_solver.py`), Tier 1/Tier 2 array schema violations, `results.json` inconsistent with the arrays, wrong mesh sequence or storage order, non-finite errors, a disallowed dependency, a linear solve too large or too wide to be a tridiagonal 1-D sweep, a solver that cannot reproduce its own submitted convergence errors, or Tier 1 / Tier 2 fail. |
| *(none)* | **`--no-exec`**: no provenance established, so the run reports `RESULT: UNGRADED` with no `SCORE:` line and exit code 2. Diagnostics only. |

Only `h_values`, `max_norm_errors` and `eoc` are mandatory arrays. The five Tier 3
arrays may be **omitted honestly** for partial credit — their absence degrades to
a Tier 3 failure (0.5), never to a schema failure.

## Visible benchmark thresholds

Absolute criteria, not a band around any reference implementation:

- **Tier 1** (`h = 0.1`): `e_h(0.1)` finite and `< 1.0`. A nonsense guard only.
  Magnitude on the coarsest mesh genuinely cannot discriminate — a solver that
  drops the mixed term scores `2.40e-3` there, *better* than the perfectly valid
  MCS at `θ = 0.8` (`2.57e-3`). All discrimination lives in Tier 2.
- **Tier 2** (convergence, the core check), all of:
  - errors strictly decreasing under refinement;
  - total error reduction `e(0.1) / e(0.00625) ≥ 100×`;
  - least-squares fitted order over all five meshes in `[1.70, 2.30]`;
  - both of the two finest experimental orders in `[1.60, 2.60]`;
  - finest-mesh error `< 1e-2` (again a nonsense guard, not a discriminator).
- **Tier 3** (finest mesh): the grid arrays equal to the true uniform 161×161
  grid (`indexing='ij'`); `U_exact_finest` equal to the analytic exact solution;
  `error_finest == U_finest − U_exact_finest`; `U_finest` carrying exact
  Dirichlet data on the boundary; `max_norm_errors[4]` equal to the error
  recomputed from `U_finest` (so the table's last row and the Tier 3 surfaces
  provably come from the same run); the finest meshes regenerated from the
  submitted solver; **and `U_finest` equal pointwise to the regenerated finest
  field** (`rtol=1e-10`, `atol=1e-12`).
- **Provenance**: `adi_solver.py` present, importable, exposing
  `solve(h, T=1.0, problem=None)`,
  importing nothing outside the standard library + NumPy + SciPy, and
  regenerating `e_h` to within `1e-6` relative of the submitted values.
  *Scoring note:* a regenerated **error** that does not match the submitted
  convergence table is a provenance failure → **0.0**. A regenerated finest
  **field** that does not match `U_finest` pointwise is a Tier 3 failure →
  **0.5**, exactly like every other Tier 3 defect: the table itself is still
  genuine, only the artifact is inconsistent. The pointwise requirement is
  listed once, under Tier 3 above.
- **ADI structure**: on an `m × m` grid, no linear system of dimension greater
  than `2(m+1)` and no system of half-bandwidth greater than 1.
- **Hidden instance** (the check that the solver *solves*): `solve()` is called on
  an undisclosed problem — different SPD coefficients and a smooth non-separable
  initial state with no closed-form solution — and must land within `3e-2` of a
  reference the grader computes with its own independent MCS solver, with the
  error falling by at least `2×` between `h = 0.1` and `h = 0.05`.

## What the evaluator checks

1. **Array schema** — `h_values`, `max_norm_errors`, `eoc` present, shaped `(5,)`,
   finite, errors positive; `h_values` exactly `0.1·2^-p` **stored in descending
   order**; `eoc[0]` NaN and `eoc[1:]` recomputed from `max_norm_errors`. Tier 3
   arrays, when present, shaped `(161,161)` and finite.
2. **Sidecar consistency** — `results.json` `convergence_table` rows (`p` indices
   `0..4`, `h`, `max_norm_error`, `eoc_order`) against the arrays;
   `diffusion_matrix` (`a=1.0`, `b=0.8`, `c=0.5`), `final_time`,
   `time_step_rule` (must actually state `k = h`, matched by regex, not merely
   contain the letters), `scheme`, `finest_mesh`, and
   `observed_order_of_convergence` supported by the submitted errors.
3. **Numerical tiers** — as tabulated above, recomputing the exact solution and
   every order analytically rather than reading the agent's claims.
4. **Dependency allowlist** — an AST pass over `adi_solver.py` rejecting any
   top-level import outside stdlib + NumPy + SciPy.
5. **Reproducibility** — `adi_solver.py` is imported in an isolated subprocess and
   `solve()` is called to regenerate the study. `h = 0.1, 0.05, 0.025` are
   regenerated before any Tier 1/Tier 2 credit; `h = 0.0125, 0.00625` are
   regenerated before any Tier 3 credit. Regenerated errors must match the
   submitted ones, the returned `(X, Y)` must be the correct grid, and the
   regenerated finest `U` must match submitted `U_finest` pointwise.
6. **The hidden instance** — 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²)` error: it reproduces its own arrays, passes the
   dimension and bandwidth audits, and shows a textbook order of 2. The grader
   therefore calls `solve(h, T, problem)` on a problem chosen by the grader and
   compares against its own independently computed reference. Calibration: the
   gold solver scores `2.0e-3 → 6.3e-4` (ratio 3.19); a solver returning zeros,
   the initial state, or the benchmark's manufactured field is off by `~2.2e-1`.
7. **Causal sweep audit** — accuracy on hidden data rules out specialization but
   does not alone prove that the observed solves produced the returned field. On
   one hidden mesh the grader replaces every intercepted solve result with NaNs
   and reruns `solve()`. A finite field unchanged to pointwise `atol=1e-12`
   (`rtol=0`) fails with `0.0`; a genuine ADI data path changes, becomes
   non-finite, or raises. If a hand-written solver exposes no interceptable call,
   the existing manual-review rule applies.
8. **ADI structure at run time** — the subprocess wraps every linear-solve entry
   point in `scipy.linalg`, `scipy.sparse.linalg`, `numpy.linalg` **and
   `scipy.linalg.lapack`** (`gtsv`, `gesv`, `getrf`, `gbsv`, …), recording both
   the dimension *and* the half-bandwidth of each system. Genuine ADI on an
   `m × m` grid solves tridiagonal systems of dimension `m−1`; the forbidden
   monolithic 2-D solve is `(m−1)²`. At `h = 0.1` that is 9 versus 81, against a
   limit of 22. Bandwidth is checked separately so that a *small but dense*
   solve is caught regardless of its dimension.

### Known limits of the automated checks

Stated plainly so reviewers know where human judgement is still required. Each
of these prints a machine-readable `REVIEW_REQUIRED:` line, and `--strict`
converts any of them into a 0.0.

- A **hand-written Thomas algorithm is explicitly permitted** by the spec and
  makes no library call, so it records zero solves and the sweep audit is
  `INCONCLUSIVE`. Default mode does not fail it (that would reject a legitimate
  submission); `--strict` does. Use `--strict` when no human will read the flags.
  The flag is raised **per mesh**, not per batch: a solver making observable
  calls on the coarse meshes and switching to an uninstrumented path on the fine
  ones is flagged for exactly those meshes, and named as the hybrid-evasion
  pattern it resembles. Any mesh with zero observable solves counts.
- **Malformed `results.json` values score 0.0, they do not crash.** Every
  agent-supplied numeric field goes through guarded accessors, so a string, a
  `null` or a non-finite value where a number belongs produces a clean
  `FAIL … SCORE: 0.0` rather than a traceback.
- The dependency allowlist is **static**. A submission using `importlib` or
  `__import__` is flagged.
- The instrumentation covers the documented NumPy/SciPy/LAPACK entry points. A
  submission calling into a compiled extension or a private code path by other
  means would not be observed.
- `--no-exec` leaves provenance entirely unverified, so it produces **no score at
  all**: the run prints `RESULT: UNGRADED` and exits with code 2. A diagnostic run
  cannot be mistaken for a graded one.

## Threshold calibration

`scheme_survey.py` runs the in 't Hout–Foulon ADI family across a wide parameter
range and applies the thresholds above. Reproduce with `python scheme_survey.py`.

**The precise coverage claim — accepted:** MCS for `θ ∈ [1/3, 8]`, Craig–Sneyd at
`θ = 1/2`, Hundsdorfer–Verwer for `θ ∈ [1/2, 8]`.

| scheme | e(0.1) | e(finest) | reduction | fitted | last two orders | verdict |
|---|---|---|---|---|---|---|
| MCS θ=1/3 | 5.353e-4 | 2.565e-6 | 209× | 1.929 | 1.954, 1.966 | accept |
| MCS θ=1/2 | 7.725e-4 | 3.417e-6 | 226× | 1.957 | 1.979, 1.990 | accept |
| MCS θ=1 | 4.628e-3 | 1.983e-5 | 233× | 1.960 | 1.871, 1.905 | accept |
| MCS θ=1.5 | 1.354e-2 | 5.818e-5 | 233× | 1.956 | 1.802, 1.847 | accept |
| MCS θ=3 | 7.760e-2 | 2.757e-4 | 282× | 2.033 | 1.896, 1.746 | accept |
| MCS θ=5 | 1.490e-1 | 7.418e-4 | 201× | 1.989 | 2.071, 1.737 | accept |
| CS θ=1/2 | 7.725e-4 | 3.417e-6 | 226× | 1.957 | 1.979, 1.990 | accept |
| HV θ=1/2 | 1.612e-4 | 9.567e-7 | 169× | 1.846 | 1.690, 1.813 | accept |
| HV θ=1/2+√3/6 | 8.964e-4 | 3.292e-6 | 272× | 2.021 | 2.007, 2.000 | accept |
| HV θ=2 | 1.369e-2 | 8.325e-5 | 164× | 1.822 | 1.738, 1.818 | accept |
| HV θ=8 | 1.586e-1 | 1.405e-3 | 113× | 1.772 | 1.848, 1.691 | accept |
| MCS θ=8 | 1.953e-1 | 1.744e-3 | 112× | 1.788 | 2.446, 1.924 | accept |
| MCS θ=20 | 2.299e-1 | 1.075e-2 | 21× | 1.058 | 1.329, 2.569 | reject |
| CS θ=1/3 | 7.5e-3 | 3.1e+40 | — | — | — | reject (unstable) |
| CS θ=0.8 | 4.401e-3 | 1.521e-4 | 29× | 1.205 | 1.113, 1.057 | reject (order 1) |
| Douglas θ=1/2 | 6.375e-4 | 3.445e-5 | 18.5× | 1.006 | 0.888, 0.947 | reject (order 1) |
| Douglas θ=1 | 7.247e-3 | 2.228e-4 | 32.5× | 1.247 | 1.138, 1.068 | reject (order 1) |
| mixed term dropped | 2.400e-3 | 1.946e-3 | 1.2× | 0.068 | 0.015, 0.004 | reject (order 0) |

**Why the accepted band ends where it does.** MCS and HV are second order
*asymptotically* for any `θ`, but at large `θ` the asymptotic regime lies beyond
`h = 0.00625`, so the order observed on the prescribed meshes stops being ≈ 2.
The band ends at **MCS θ = 10**, which achieves only a 79.9× total reduction and a
fitted order of 1.652 — short of the published 100× and `[1.70, 2.30]`. MCS at
`θ = 8` is *accepted*: its pointwise orders run `[0.62, 1.82, 2.45, 1.92]`, an
untidy pre-asymptotic transient, but its total reduction (112×) and fitted order
(1.788) clear the published bar, and the upper pointwise limit of `2.60` is a
sanity guard rather than a discriminator — measured across this whole family it
excludes nothing except large-`θ` MCS, while every genuine failure is caught by
the *lower* limit, the reduction test, or monotonicity. The instruction requires
the observed order to approach 2 on exactly `h = 0.1·2^-p`, `p = 0..4`, so
rejecting `θ ≥ 10` is correct on the task's own stated terms, not a hidden
preference for one method. Craig–Sneyd is second order **only** at `θ = 1/2`; away
from it CS degrades to order ≈ 1.2, which is why it appears on both sides of the
table.

*For information only, not a threshold:* the bundled reference implementation
(MCS, θ = 1/3) produces `5.353e-4, 1.478e-04, 3.885e-05, 1.002e-05, 2.565e-06`.
No acceptance criterion is derived from these numbers.

## What a correct run looks like

A stable ADI scheme that incorporates the mixed term and demonstrates
second-order convergence on the prescribed mesh sequence — any member of the
CS / MCS / HV family in the accepted parameter range; all three required files
written; `adi_solver.py` exposing `solve(h, T=1.0, problem=None)` and reproducing the submitted
errors on all five meshes when re-executed; only tridiagonal 1-D sweeps solved
implicitly; arrays truthful and internally consistent; Tier 3 either passing for
full credit or omitted / failing honestly for partial credit.

## What a wrong run looks like

Missing `adi_solver.py` or the Tier 1/Tier 2 outputs; `results.json` inconsistent
with the arrays; `h_values` stored in the wrong order or with wrong `p` indices; a
hand-edited `eoc` array; Tier 3 surfaces spliced from a different run than the
convergence table; a solver that cannot reproduce the numbers submitted alongside
it; a scheme that ignores or naively lags the mixed derivative (order → 0, error
plateau near `2e-3`); a first-order scheme such as Douglas or Craig–Sneyd off
`θ = 1/2`; an unstable parameter choice such as CS at θ = 1/3 (overflow); a
disallowed PDE library; or a directly-assembled 2-D sparse solve in place of
tridiagonal ADI sweeps — detected even when it is *more* accurate than the
reference, because the recorded system dimension and bandwidth give it away.