# Derivative Backend Policy — Design Document

**Status:** DESIGN-FROZEN in v0.30a; **`finite_difference` v1 IMPLEMENTED in v0.30c** along with the `compute_derivatives(backend="auto")` dispatcher. Boundary-aware readiness warnings and the `residual_domain_policy` reporting field also land in v0.30c. **Residual evaluator auto-dispatch and interior-only diagnostics landed in v0.30d**: Heat, Burgers, advection-diffusion, and reaction-diffusion evaluators now route through `compute_derivatives(backend="auto")` and consume `recommended_residual_domain_policy` / `recommended_boundary_trim_width` from the `DerivativeBatch.config`. KdV and the weak evaluators remain periodic-only per v0.30 scope.

This document defines the policies governing PDELie's derivative backends after v0.30. The runtime in v0.30a is unchanged.

## Backends covered

After v0.30, two derivative backends exist:

- `spectral_fd` (existing): FFT-based spectral spatial derivatives + 2nd-order centered finite differences in time. Periodic-x only.
- `finite_difference` (new in v0.30 proper): centered finite differences with one-sided stencils at boundaries. Nonperiodic-x only.

A third name, `compute_derivatives(..., backend="auto")`, is a dispatcher, not a backend.

## `spectral_fd` policy

**Unchanged from current behavior.**

- Periodic-x is required. The check at `src/pdelie/derivatives/spectral_fd.py:32-33` stays in place.
- Spatial derivatives are computed via FFT phase factors `(ik)^n`. This is mathematically correct only for periodic data; loosening the check would silently produce garbage at the boundaries due to Gibbs ringing.
- Output is unchanged: `DerivativeBatch(backend="spectral_fd", ...)`.

After v0.30 the BC check reads `metadata["boundary_conditions"]["x"]["type"] == "periodic"` against the new structured spec, not the legacy string. The semantics are identical for callers.

## `finite_difference` policy (v0.30 proper)

A new module `src/pdelie/derivatives/finite_difference.py` will implement `compute_finite_difference_derivatives(field, *, max_spatial_order=2)`.

### Accepted BCs

- `"dirichlet"`
- `"neumann"`
- `"open_unknown"`

Periodic data passed to `finite_difference` raises `ScopeValidationError`. Periodic users must use `spectral_fd`. The dispatcher (below) handles the auto-routing.

### Order limits

The stable v0.30 surface supports:

- `u_t` via `np.gradient(values, dt, axis=t_axis, edge_order=2)` (unchanged from the time-derivative computation in `spectral_fd`)
- `u_x` via `np.gradient(values, dx, axis=x_axis, edge_order=2)`
- `u_xx` via `np.gradient(u_x, dx, axis=x_axis, edge_order=2)`

`u_xxx` and `u_xxxx` are **not** part of the stable v0.30 surface on nonperiodic data. Calling `compute_finite_difference_derivatives(field, max_spatial_order=3)` or higher raises `ScopeValidationError` unless an `experimental_high_order=True` flag is set, and even then the output is documented as advisory and may change.

#### Why no high-order FD on nonperiodic data

Repeated `np.gradient` propagates boundary error. The one-sided edge stencil at order 2 has the same asymptotic O(h²) interior accuracy as centered FD, but a larger constant. Each derivative step pulls the boundary error one stencil width into the interior. By `u_xxxx` on N≈64-point grids, the boundary-corrupted region dominates the interior region; the residual statistics become unreliable.

For periodic data the same chain works because FFT-based derivatives have no boundary problem at all. For nonperiodic data the right answer is a higher-order one-sided stencil family (Fornberg-style), or weak/Galerkin integration to bypass derivative estimation entirely. Neither is in scope for v0.30; both are candidates for later releases.

### `DerivativeBatch.config` requirements

The backend must populate:

```
{
    "spatial_method": "finite_difference_centered",
    "temporal_method": "finite_difference",
    "temporal_edge_order": 2,
    "spatial_edge_order": 2,
    "spatial_max_order": <int>,
    "boundary_handling": "dirichlet" | "neumann" | "open_unknown",
    "boundary_left_specified": <bool>,
    "boundary_right_specified": <bool>,
}
```

`boundary_left_specified` and `boundary_right_specified` read directly from the structured `BoundaryConditionSpec` faces. They allow downstream consumers (residual evaluators, readiness reports) to flag derivative results that depended on an unspecified boundary.

`DerivativeBatch.boundary_assumptions` records `"<bc_type> in x; finite differences in time"` for human inspection.

## Dispatcher: `compute_derivatives(field, *, backend="auto", max_spatial_order=2)`

A new helper, exported from `pdelie.derivatives`, selects the right backend based on the field's BC.

### Dispatch table

- `field.metadata["boundary_conditions"]["x"]["type"] == "periodic"` → `spectral_fd`
- `field.metadata["boundary_conditions"]["x"]["type"] in {"dirichlet", "neumann", "open_unknown"}` → `finite_difference`

### Explicit backend selection

- `backend="spectral_fd"` on nonperiodic data → `ScopeValidationError`
- `backend="finite_difference"` on periodic data → `ScopeValidationError`
- `backend="auto"` always succeeds for any supported BC type

### Never silent

The dispatcher must record its selection in `DerivativeBatch.config`:

```
{
    ...,
    "backend_selected_by_boundary_condition": True,
    "backend_selection_reason": "periodic_x_boundary" | "nonperiodic_x_boundary",
}
```

This is non-negotiable. A user reading a `DerivativeBatch` from disk must be able to tell which backend produced it without re-deriving the choice from `boundary_conditions`.

### Order capping under dispatch

When `backend="auto"` resolves to `finite_difference`, requesting `max_spatial_order > 2` raises `ScopeValidationError` with the same message as the direct path. The dispatcher does not silently downgrade orders.

## Residual evaluation domain policy

The existing residual evaluators (`HeatResidualEvaluator`, `BurgersResidualEvaluator`, `AdvectionDiffusionResidualEvaluator`, etc.) currently compute residuals over the full grid. On nonperiodic FD-derived data, near-boundary residuals are dominated by one-sided-stencil error and are not representative of interior residual quality.

### New summary field

`summarize_residual_batch` (in `src/pdelie/reporting/summaries.py`) gains a new field `"residual_domain_policy"` with vocabulary:

- `"full_grid"` — RMS / max over all grid points (the historical behavior for periodic data)
- `"interior_only"` — RMS / max over the interior excluding a boundary band
- `"drop_boundary_width_k"` — explicit boundary width `k` dropped from both ends

### Default policy

- Periodic data: `"full_grid"`. No boundary problem exists for spectral derivatives, so no need to trim.
- Nonperiodic FD-derived data: `"interior_only"` with default boundary width `k = 4`. The choice of `k = 4` matches the effective stencil width of `u_xx` computed as `np.gradient(np.gradient(u))`: centered second derivative consumes a 5-point stencil, edge-order-2 handling at each boundary consumes up to 3 extra points of error propagation; `k = 4` provides a safe margin.

### Reporting both

When a residual is computed on nonperiodic data, the summary SHOULD include both `"residual_domain_policy": "interior_only"` (the primary statistic) and a nested `"full_grid_diagnostic"` block with the same metrics computed over the full grid, so users can see the boundary contamination explicitly. This is recommended, not required, in v0.30 proper.

## Translation finite-transform policy on nonperiodic data

`verify_translation_generator` at `src/pdelie/verification/finite_transform.py:64` and `_apply_uniform_translation` at `:21` remain **periodic-only**. The FFT phase-shift method requires periodicity; there is no correct alternative on a finite domain without a domain policy.

### What is explicitly rejected

A `_apply_uniform_translation_open` that uses `np.roll` + extrapolation is **not approved**. `np.roll` does not translate a non-periodic function — it wraps values that have no physical meaning at the wrap point, and any extrapolation scheme (constant, linear, reflective) introduces O(1) artifacts in the residual.

### What is deferred

The correct approach for nonperiodic translation is overlap-crop:

1. Apply translation `ε`: shift the coordinate axis by `ε`.
2. Crop the field to the overlap interval `[a + |ε|, b − |ε|]` where `[a, b]` is the original domain.
3. Compute the residual on the cropped interior only.

This requires a domain policy (how to handle multiple overlapping translations, what minimum overlap is acceptable, how to compare orbits of different lengths). The full design is deferred to **v0.31.5 — Nonperiodic Orbit/Action Scope Decision** in the roadmap.

## Numerical sanity tests required for v0.30 proper

These are not implemented in v0.30a (it ships no runtime code). They are specified here so v0.30 proper has a clear acceptance contract.

### Manufactured analytic functions

For each backend × order combination, test against closed-form derivatives of:

- `u(t, x) = sin(k·x − ω·t)` on `[0, 2π]` (periodic check for spectral; reflected to `[0, π]` for FD where periodicity does not apply)
- `u(t, x) = tanh(α·x − v·t)` on `[−5, 5]` (smooth, non-periodic, Dirichlet boundaries known)
- `u(t, x) = x³ − 3xt²` on `[−1, 1]` (polynomial, exact closed-form derivatives at all orders)

For each test:
- Compare `compute_finite_difference_derivatives` outputs against the closed-form `u_t`, `u_x`, `u_xx` evaluated at the grid points.
- Assert max-absolute error in the interior decreases at rate O(h²) as the grid is refined.
- Assert max-absolute error in the **boundary band** (first/last `k` points for `k = 4`) is bounded but does **not** assert O(h²) — boundary stencils have a larger error constant.

### What not to test

**Do not test against `np.gradient`.** The new backend is built on `np.gradient` internally; testing it against itself is meaningless. Tests must use closed-form analytic derivatives.

### Convergence study

For at least one manufactured solution, run the backend at `N = 32, 64, 128, 256` and assert the interior-only max error halves to within 10% per grid doubling (the expected O(h²) order). The boundary-band max error must not blow up — a coarse bound like `< 10× interior_max` is sufficient.

## Implementation pointer (for v0.30 proper)

The most surgical path:

1. Add `src/pdelie/derivatives/finite_difference.py` with `compute_finite_difference_derivatives` mirroring the structure of `src/pdelie/derivatives/spectral_fd.py`.
2. Add `compute_derivatives` to `src/pdelie/derivatives/__init__.py` as a thin dispatcher.
3. Update `src/pdelie/residuals/heat_1d.py`, `burgers_1d.py`, `advection_diffusion_1d.py`, `reaction_diffusion_1d.py` to call `compute_derivatives(..., backend="auto")` instead of `compute_spectral_fd_derivatives(...)` when no derivatives are supplied. KdV stays periodic-only (its evaluator already requires `equation == "kdv_normalized"`, which only the periodic generator produces).
4. Update `summarize_residual_batch` to surface `residual_domain_policy`.
5. Add `compute_finite_difference_derivatives` to `tests/test_public_api.py::test_runtime_package_api_is_importable` and `tests/test_api_stability_audit.py::_ROOT_RUNTIME_NAMES`-style negative checks (it must not be a root export).
6. Update `docs/specs/API_STABILITY.md` with a "Runtime public API for the frozen `v0.30` Milestone 2 slice" subsection describing `compute_finite_difference_derivatives` and `compute_derivatives`.

The lazy-import pattern at `src/pdelie/discovery/pysindy_adapter.py:12-20` does not apply: `finite_difference` uses only `numpy`, which is already a core dependency.