WIP: SEP 44: read the transverse edge where the state lives, and what else a native run pays the host for

AGENTS.md

@@ -216,9 +216,9 @@ The long form is [SEP 0](docs/dev/seps/sep-0000-process.md); the index is
   starts and marked ready when the last commit sets `Final`; there are no per-phase pull
   requests. CI runs on every push to it and need not be green while it is a draft; a push
   that needs no run has `[skip ci]` in its head commit's message. The branch may start on
-  top of an open pull request's branch when pull requests are stacked behind CI. Within
-  each phase: `numpy_ref` first with its tests and manual page, then native CPU in
-  `native/`; `mlx` only where the SEP says so (the backend rules, below). A
+  top of another agent's open pull request, when the SEP needs that work first. Within
+  each phase: `numpy_ref` and native CPU in `native/`, in either order (ADR-068), with
+  tests and the manual page; `mlx` only where the SEP says so (the backend rules, below). A
   backend without the feature refuses it by name, never falls back silently; `Final`
   needs parity on `numpy_ref` and native CPU, or an exception written into the accepted
   SEP.
@@ -259,7 +259,8 @@ to it; a function that departs from it says so in the first line of its docstrin
   sentinel where raising is correct.
 - Respect the layers (`architecture.md §3`): physics code never contains performance
   tricks; backend code never contains physics knowledge; the kernel contract
-  (`kernel-contract.md`) is the only crossing.
+  (`kernel-contract.md`) and the protected driver surface (frozen in
+  `tests/backends/test_driver.py`) are the crossings.
 
 ## Environment and dependencies
 
@@ -317,13 +318,15 @@ no blanket `# noqa`). What is left to judgement:
 
 The long form is [`backends.md`](docs/dev/backends.md).
 
-- **Where features land** (ADR-061, ADR-064, `backends.md §10`). `numpy_ref` and `native`:
-  every new feature lands on both, and they come first for testing and development.
-  `mlx`: experimental; keep it runnable and correct, add nothing as part of other work. A
-  backend without a feature refuses it by name.
-- The **NumPy reference backend is the oracle**. Every other backend is validated against
-  it (one linear application, one nonlinear evaluation, one step, eight steps, each
-  statistic, a full run) with tolerances derived from summation bounds and asserted `> 0`.
+- **Where features land** (ADR-061, ADR-064, ADR-068, `backends.md §10`). `native` is the
+  primary backend. `numpy_ref` and `native`: every new feature lands on both, in either
+  order, and `native` may shape the contract. `mlx`: experimental; keep it runnable and
+  correct, add nothing as part of other work. A backend without a feature refuses it by
+  name.
+- The **NumPy reference backend is the oracle** for numbers. Every other backend is
+  validated against it (one linear application, one nonlinear evaluation, one step, eight
+  steps, each statistic, a full run) with tolerances derived from summation bounds and
+  asserted `> 0`.
 - Hot code lives only under `src/solprop/backends/`, `src/solprop/translate/` and
   `native/`.
 - **Native kernels** follow [`native/AGENTS.md`](native/AGENTS.md). The floating-point
@@ -353,11 +356,16 @@ The long form is [`backends.md`](docs/dev/backends.md).
     tolerance. Adaptive work counts may differ; accounting for actual work is exact.
 - **In-place and `out=`** only inside the backends, wrapped by the seam; every buffer is
   listed by `memory_ledger()`, which is tested.
-- Device arrays never cross the seam; backends take and return host arrays.
-- **`KERNEL_SCHEMA_VERSION` versions the whole seam** (the `KernelProblem` fields, the
-  protected `StepDriver` surface and the stage bundles of `translate/stage.py`, frozen in
-  `tests/backends/test_driver.py`). Changing any of them bumps the version, rewrites the
-  literal in the same commit and adds a release-note fragment.
+- Device arrays never cross the `Backend` protocol, whose methods take and return host
+  arrays. Between `StepDriver` and its subclass the state is an opaque handle (an MLX
+  array, a native state), and the shared code reads it only through `_observed_field`,
+  which copies it to the host.
+- **The seam has no version number** (ADR-068). A change to the `KernelProblem` fields,
+  the protected `StepDriver` surface or the stage bundles of `translate/stage.py` rewrites
+  the frozen tables in `tests/backends/test_driver.py` in the same commit and adds a
+  release-note fragment. solprop and solprop-native are released together and pinned to
+  one version, and an editable checkout's extension is compared with its `native/` sources
+  by digest (`check_build`).
 - Every hot-path change comes with a benchmark ledger entry
   ([`benchmarks.md`](docs/dev/benchmarks.md)), measured with `python -m benchmarks ab`
   on about ten steps or a short `ab --length` window; report the window measured and
@@ -471,6 +479,10 @@ The long form is [`testing.md`](docs/dev/testing.md).
   commits rather than one large one. Stage deliberately with `git add <paths>`.
 - Changes reach `main` as pull requests with green CI, or by a local merge when the
   maintainers ask for one.
+- **Few pull requests, few pushes.** CI runs its whole matrix on every push to every pull
+  request, so keep one agent's work in one pull request, or a few, rather than one per
+  step. Stack only on another agent's open pull request, when your work needs it before it
+  lands. Commit as often as you like and push in batches.
 - **Before opening an issue or a pull request,** read [`CONTRIBUTING.md`](CONTRIBUTING.md):
   search the issues, the pull requests and the SEPs, link what is related, and follow its
   sections and labels.

benchmarks/device_split.py

@@ -0,0 +1,292 @@
+# SPDX-FileCopyrightText: 2026 John C. Travers and contributors
+#
+# SPDX-License-Identifier: GPL-3.0-or-later
+
+"""Split a field-free point of a Kerr scan: the march, the check, the crossings.
+
+One input of a delay scan through :class:`solprop.driver.SimulationRunner`, with no
+fields kept: a fused-silica Kerr beam at 260 nm on fixed low-storage RK4IP steps, with
+and a stated input identity. It runs with ``boundary_check=None`` and
+with the default check, budgets relaxed so that nothing stops, and profiles one point
+of each. Then it times what a plane can ask of a device: a state's copy to the host, an
+observation and the common statistics.
+
+Usage::
+
+    uv run python -m benchmarks.device_split NX [--steps 8] [--planes 4] [--out FILE]
+
+It prints one JSON object per section (``record``, ``setup``, ``point`` twice,
+``calls``), and writes them to ``FILE`` too when given. ``point`` names each function's
+calls, own and cumulative seconds, and ``calls`` each crossing's minor page faults.
+The profiler inflates Python-heavy work, so compare ``walls_s``, which it doesn't run
+under. ``calls`` also splits a width: its native moments, and what the same azimuthal
+average takes on the host. The native backend is required, and a device build is what
+it's for; on a CPU build the crossings are plain copies.
+"""
+
+from __future__ import annotations
+
+import argparse
+import cProfile
+import json
+import platform
+import pstats
+import statistics
+import subprocess
+import sys
+import time
+import warnings
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+
+from solprop.constants import EPS0, MU0, C
+from solprop.driver import Simulation, SimulationSettings
+from solprop.grid import EnvelopeGrid, FreeGrid
+from solprop.monitor import BoundaryCheck
+from solprop.nonlinear.responses import KerrResponse
+from solprop.physdata.glasses import SellmeierGlass
+from solprop.pulses import GaussianPulse, SpatialPulse
+from solprop.spatial import SpatialMedium, cartesian_fluence_widths
+from solprop.stats import StatsSpec, default_spec
+from solprop.stepper import FixedSchedule, LowStorageRungeKutta, PIController
+
+__all__ = ["main"]
+
+STEP = 0.5e-6
+"""The fixed step, m."""
+
+NAMED = (
+    "_read_transverse",
+    "_edge_reading",
+    "node_power",
+    "column_power",
+    "_spacing",
+    "_observed_field",
+    "_boundary_energies",
+    "_take_step_stats",
+    "_monitor_plane",
+    "configuration_hash",
+    "_hashed_record",
+    "_check_replaced_input",
+    "attach_monitor",
+    "replay_step",
+    "c_einsum",
+    "field_to_host",
+)
+"""The functions a point's profile reports by name."""
+
+GLASS = SellmeierGlass(
+    name="malitson",
+    source="Malitson 1965",
+    wavelength_range=(0.21e-6, 3.71e-6),
+    citation="Malitson, JOSA 55, 1205 (1965), doi:10.1364/JOSA.55.001205",
+    poles=((0.6961663, 0.0684043), (0.4079426, 0.1162414), (0.8974794, 9.896161)),
+)
+
+
+def _shell(command: str) -> str:
+    return subprocess.run(
+        command, shell=True, capture_output=True, text=True, check=False
+    ).stdout.strip()
+
+
+def _simulation(
+    nx: int, steps: int, planes: int, check: BoundaryCheck | None
+) -> Simulation:
+    length = STEP * steps
+    half = 183e-6 * nx / 384 * 1.1
+    grid = EnvelopeGrid.from_limits(
+        length, 260e-9, (143e-9, 600e-9), 200e-15, thg=False, fft_size="smooth"
+    )
+    positions = [length * (index + 1) / planes for index in range(planes)]
+    settings = SimulationSettings(
+        spatial_absorber=False,
+        controller=PIController(rtol=1e-6, atol=1e-14),
+        save_planes=positions,
+        dz_max=length,
+        dz_init=STEP,
+        schedule=FixedSchedule.uniform(0.0, length, STEP, save_planes=positions),
+        integrator=LowStorageRungeKutta(),
+        backend="native",
+        boundary_check=check,
+        fields="none",
+        stats=StatsSpec(names=("energy",)),
+    )
+    return Simulation.from_free(
+        grid,
+        FreeGrid(half, nx, half, nx, collar_fraction=0.0),
+        SpatialMedium(GLASS),
+        [KerrResponse(chi3=2.2e-22 * MU0 * EPS0 * C**2, thg=False)],
+        SpatialPulse(GaussianPulse(wavelength=260e-9, fwhm=30e-15, energy=1e-7), 40e-6),
+        settings=settings,
+    )
+
+
+def _profile_summary(profile: cProfile.Profile, top: int = 30) -> dict[str, Any]:
+    rows, named = [], {}
+    # `Stats.stats` is pstats' documented table, untyped in its stubs.
+    table: dict[Any, Any] = pstats.Stats(profile).stats  # pyright: ignore[reportAttributeAccessIssue]
+    for (path, line, function), (_, calls, own, cumulative, _) in table.items():
+        where = f"{function}@{Path(path).name}:{line}"
+        rows.append((own, cumulative, calls, where))
+        if any(name in function for name in NAMED):
+            named[where] = [calls, round(own, 3), round(cumulative, 3)]
+    rows.sort(reverse=True)
+    return {
+        "named": named,
+        "top_own": [[round(a, 3), round(b, 3), c, d] for a, b, c, d in rows[:top]],
+    }
+
+
+def _minor_faults() -> int:
+    """Return this process's minor page faults so far: none are counted on Windows.
+
+    Returns
+    -------
+    int
+        The count.
+    """
+    if sys.platform == "win32":
+        return 0
+    import resource  # noqa: PLC0415 - not on Windows
+
+    return resource.getrusage(resource.RUSAGE_SELF).ru_minflt
+
+
+def _timed(call: Any, repeats: int) -> dict[str, float]:
+    times, faults = [], 0
+    for _ in range(repeats):
+        before = _minor_faults()
+        started = time.perf_counter()
+        value = call()
+        times.append(time.perf_counter() - started)
+        faults = _minor_faults() - before
+        del value
+    return {
+        "median_ms": round(1e3 * statistics.median(times), 2),
+        "first_ms": round(1e3 * times[0], 2),
+        "minor_faults": faults,
+    }
+
+
+def _calls(backend: Any) -> dict[str, Any]:
+    specs = {
+        name: StatsSpec(names=(name,))
+        for name in (
+            "energy",
+            "beam_spectrum",
+            "absorbed_power",
+            "on_axis_intensity",
+            "peak_intensity",
+            "beam_radius",
+            "fwhm_r",
+        )
+    }
+    specs["default"] = default_spec("free")
+    geometry = backend._problem.transform.geometry
+    x, y = geometry.axes
+    state = backend._native_state
+    return {
+        "crossings": {
+            "field_to_host": _timed(state.field_to_host, 3),
+            "observe_field": _timed(lambda: backend.observe(field=True), 3),
+            "observe_free": _timed(lambda: backend.observe(field=False), 3),
+        },
+        "statistics": {
+            name: _timed(lambda spec=spec: backend.stats(spec), 5)
+            for name, spec in specs.items()
+        },
+        # A width's parts: the fluence and its moments, one native call; and what the
+        # host's azimuthal average of a Gaussian's fluence takes on this machine.
+        "widths": {
+            "moments": _timed(
+                lambda: state.spatial_statistics(32, 0.5 * EPS0 * C, backend._z), 5
+            ),
+            "host_average": _timed(
+                lambda: cartesian_fluence_widths(
+                    x,
+                    y,
+                    np.exp(-2 * (x[None, :] ** 2 + y[:, None] ** 2) / 40e-6**2).ravel(),
+                ),
+                3,
+            ),
+        },
+    }
+
+
+def main(argv: list[str] | None = None) -> None:
+    """Run the split and print its sections.
+
+    Parameters
+    ----------
+    argv : list of str or None, optional
+        The arguments, by default the process's.
+    """
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("nx", type=int, help="transverse samples per axis")
+    parser.add_argument("--steps", type=int, default=8)
+    parser.add_argument("--planes", type=int, default=4)
+    parser.add_argument("--out", type=Path, default=None)
+    arguments = parser.parse_args(argv)
+    warnings.simplefilter("ignore")
+    sink = arguments.out.open("w", encoding="utf-8") if arguments.out else None
+
+    def emit(section: str, **values: Any) -> None:
+        line = json.dumps({"section": section, "nx": arguments.nx, **values})
+        sys.stdout.write(line + "\n")
+        sys.stdout.flush()
+        if sink is not None:
+            sink.write(line + "\n")
+            sink.flush()
+
+    import solprop_native  # noqa: PLC0415  # pyright: ignore[reportMissingImports]
+
+    emit(
+        "record",
+        cpu=_shell("lscpu | grep 'Model name' | cut -d: -f2 | xargs"),
+        cores=_shell("nproc"),
+        uptime=_shell("uptime"),
+        machine=platform.machine(),
+        transparent_hugepages=_shell("cat /sys/kernel/mm/transparent_hugepage/enabled"),
+        commit=_shell("git rev-parse --short HEAD"),
+        device=solprop_native.build_info().device_backend,
+    )
+    started = time.perf_counter()
+    base = _simulation(arguments.nx, arguments.steps, arguments.planes, None)
+    built = time.perf_counter()
+    problem = base.translate()
+    translated = time.perf_counter()
+    launch = np.array(problem.initial.field0)[0]
+    emit(
+        "setup",
+        from_free_s=round(built - started, 2),
+        translate_s=round(translated - built, 2),
+        field_GiB=round(launch.nbytes / 2**30, 3),
+        shape=list(launch.shape),
+    )
+    relaxed = BoundaryCheck(budget=1.0, transverse_budget=1.0)
+    for name, check in (("off", None), ("on", relaxed)):
+        simulation = _simulation(arguments.nx, arguments.steps, arguments.planes, check)
+        with simulation.runner() as runner:
+            walls = []
+            for point in range(2):
+                field = launch * np.exp(1j * point)
+                started = time.perf_counter()
+                runner.run(field, progress="none", input_identity=f"point-{point}")
+                walls.append(round(time.perf_counter() - started, 3))
+            field = launch * np.exp(0.5j)
+            profile = cProfile.Profile()
+            profile.enable()
+            runner.run(field, progress="none", input_identity="profiled")
+            profile.disable()
+            emit("point", check=name, walls_s=walls, **_profile_summary(profile))
+            if check is None:
+                emit("calls", **_calls(runner.backend))
+    if sink is not None:
+        sink.close()
+
+
+if __name__ == "__main__":
+    main()