Further optimzations of flow tracking

development/flow/PORTING_CAMPAIGNS.md

@@ -0,0 +1,113 @@
+# Flow tracer: porting campaigns for Apple Silicon and MSVC
+
+Two self-contained prompts to start an optimization session on a machine that
+is not covered by the Linux/GCC measurements in `PERFORMANCE_NOTES.md`. Each
+prompt assumes a checkout of branch `flow-optim-fable` (commit `54fd4a9`) and
+the harnesses in this directory.
+
+State at the time of writing:
+
+- **arm64 macOS** gets no specialized tracer. The CMake block
+  `BIOIMAGE_FLOW_FMA_DISPATCH` matches only x86, so the portable scalar header
+  kernel runs. The packed kernel uses 26 distinct 128-bit intrinsics that all
+  map 1:1 onto NEON, and NEON+FMA is baseline on arm64, so no runtime dispatch
+  is needed.
+- **Windows/MSVC** already builds and dispatches the packed kernel
+  (`/arch:AVX2`, `__cpuid`/`_xgetbv`), but it has never been measured there.
+  The open items are inlining under `__forceinline`, COMDAT/ODR leakage, and
+  the `/fp:precise` contraction question for the scalar fallback.
+
+## Prompt: Apple Silicon (arm64 clang)
+
+```text
+Context: bioimage-cpp, branch flow-optim-fable (commit 54fd4a9). On x86 the flow
+tracer (`bic.flow.compute_flow_density`) has a packed-channel 128-bit SSE+FMA kernel in
+src/cpp/flow/flow_density_fma.cxx, dispatched from `try_trace_all_fma` in
+include/bioimage_cpp/flow/flow_density.hxx, giving -30..-40 % on Zen 3 and Tiger Lake.
+On arm64 nothing is specialized: the CMake block `BIOIMAGE_FLOW_FMA_DISPATCH` matches
+only x86, so this Mac runs the portable scalar header kernel. Read
+development/flow/PERFORMANCE_NOTES.md, section "Packed-channel FMA tracer (2026-09-02)",
+for the design, the per-step diagnosis, and the rejected experiments before you start.
+
+Goal: a NEON port of the packed kernel for arm64 macOS, gated on the same
+correctness bar, with paired numbers recorded in PERFORMANCE_NOTES.md.
+
+Steps:
+1. Build (pip install -e . --no-build-isolation) and confirm
+   `_core._flow_trace_backend()` reports "scalar". Run
+   development/flow/paired_bench.py --a X.so --b X.so --cases fixture3d,fixture2d
+   without --cpu (no taskset/sched_setaffinity here); every row must read "noise".
+   If single-thread spread exceeds ~3 %, address P/E-core placement or use more repeats.
+2. Baseline the scalar kernel at 1T and at the physical P-core count. Write the
+   numbers down before changing anything.
+3. Port src/cpp/flow/flow_density_fma.cxx to <arm_neon.h> in a new TU
+   (e.g. src/cpp/flow/flow_density_neon.cxx). All 26 intrinsics used are 128-bit and
+   map 1:1 (vld1q_f32, vfmaq_f32, vcvtq_s32_f32, vgetq_lane_s32, vminq/vmaxq,
+   vcgtq/vcltq + vbslq, vextq/vzipq/vdupq_laneq, a movemask idiom via vshrn/vaddvq).
+   Check the truncating-convert difference: NEON saturates and maps NaN to 0, SSE
+   returns INT_MIN; verify the clip/bounds/mask tests do not depend on SSE behavior.
+   No runtime dispatch and no CPUID are needed (NEON+FMA is baseline on arm64);
+   keep BIOIMAGE_CPP_FLOW_FORCE_SCALAR as the opt-out and make trace_backend()
+   return "neon". Keep BIOIMAGE_FORCE_INLINE on the per-step helpers and confirm
+   with `otool -tV` that nothing in the hot loop is called out of line.
+4. Re-sweep the lockstep lane count K (4, 6, 8): the x86 choice of 4 came from an
+   FP-register-file bottleneck that M-series cores may not have.
+5. Gate: differential_check.py --a base.so --b neon.so (400 cases, sums equal,
+   diff fraction below 1 %) and --env-b BIOIMAGE_CPP_FLOW_FORCE_SCALAR=1;
+   check_flow_density.py --dim both; full pytest; the parity test
+   tests/test_flow.py::test_forced_scalar_matches_fma_backend must still pass or be
+   generalized to the "neon" name.
+6. Decide, based on measured parity, whether to unify the SSE and NEON kernels behind
+   a small detail/simd4.hxx wrapper or keep two TUs. Do not unify before you have
+   numbers for both.
+7. Add a "Reproduction on Apple Silicon" subsection to PERFORMANCE_NOTES.md with
+   the paired table, K sweep, and anything rejected. Do not touch MIGRATION_GUIDE.md
+   unless the public API changes.
+```
+
+## Prompt: Windows x86-64 (MSVC)
+
+```text
+Context: bioimage-cpp, branch flow-optim-fable (commit 54fd4a9). The flow tracer
+(`bic.flow.compute_flow_density`) has a packed-channel SSE+FMA kernel in
+src/cpp/flow/flow_density_fma.cxx that CMake already builds under MSVC with
+/arch:AVX2 and dispatches at runtime via __cpuid/_xgetbv in
+include/bioimage_cpp/flow/flow_density.hxx (`runtime_fma_supported`,
+`try_trace_all_fma`). It has only been measured with GCC on Linux (-30..-40 %).
+Read development/flow/PERFORMANCE_NOTES.md, section "Packed-channel FMA tracer
+(2026-09-02)", especially "Codegen trap" and the residual MSVC caveats.
+
+Goal: validate and, where needed, tune the kernel under MSVC on Windows x86-64, and
+record paired numbers in PERFORMANCE_NOTES.md.
+
+Steps:
+1. Build with the MSVC toolchain cibuildwheel uses (pip install -e . --no-build-isolation
+   from a VS developer shell). Confirm `_core._flow_trace_backend()` reports "fma".
+   Note that CMakeLists.txt passes -O3, which cl ignores; confirm /O2 /Ob2 are in
+   effect for the Release config.
+2. Make development/flow/paired_bench.py and differential_check.py run on Windows:
+   they call `taskset` when --cpu is given, so add a Windows pinning branch
+   (psutil cpu_affinity or `start /affinity`) or run without --cpu. Set the power
+   plan to High performance. Calibrate with --a X.pyd --b X.pyd; every row must read
+   "noise" before any A/B is trusted.
+3. Static codegen check with `dumpbin /disasm` on the object files: (a) every
+   BIOIMAGE_FORCE_INLINE helper (`sample_flow`, `round_to_flat_index`,
+   `position_is_in_mask`, `trace_particle`, `trace_particle_block`, the TU's
+   `step_packed`/`sample_packed`) must be inlined into the hot loop; (b) no VEX-encoded
+   instructions may appear outside the FMA TU's own functions, since /OPT:ICF and
+   COMDAT folding can otherwise pick an AVX copy of a shared helper for the
+   portable path. If (a) fails, that is the first thing to fix (on GCC it cost 1.5-2x).
+4. Paired A/B: build the baseline commit 492ff9b and the branch into separate dirs,
+   swap the two _core .pyd files, run paired_bench.py on fixture3d/fixture2d and the
+   adversarial cases at 1T and at the physical core count. Also measure
+   BIOIMAGE_CPP_FLOW_FORCE_SCALAR=1 against the branch: on MSVC /fp:precise the
+   scalar-FMA fallback is not expected to contract, so it may equal the portable kernel.
+5. Gate: differential_check.py base vs branch and fma vs FORCE_SCALAR (400 cases,
+   sums equal, diff fraction below 1 %); check_flow_density.py --dim both; full
+   pytest. Report whether densities are bitwise identical to the GCC build on the
+   fixtures, and if not, whether the difference is confined to contraction order.
+6. Only if step 3 or 4 shows a real gap: tune (inlining pragmas, /Ob3, lane count K),
+   each change measured paired and gated as in step 5.
+7. Add a "Reproduction on Windows/MSVC" subsection to PERFORMANCE_NOTES.md with the
+   paired table, the dumpbin findings, and anything rejected.
+```

development/flow/_flow_cases.py

@@ -0,0 +1,176 @@
+"""Shared inputs and helpers for the flow-density development harnesses.
+
+Used by ``perf_kernel.py``, ``paired_bench.py`` and ``differential_check.py``.
+Not part of the package or the test suite.
+
+``preload_core`` lets a harness load ``bioimage_cpp._core`` from an explicit
+``.so`` path before ``import bioimage_cpp`` runs. This is how two prebuilt
+kernels are compared without rebuilding between A and B (rebuild rounds drift
+thermally and the editable install otherwise loads ``_core`` from
+site-packages).
+"""
+
+from __future__ import annotations
+
+import importlib.util
+import os
+import sys
+from pathlib import Path
+from typing import Callable
+
+import numpy as np
+
+DEFAULTS: dict = {
+    "n_iter": 50,
+    "dt": 0.2,
+    "tol": 0.005,
+    "method": "rk2",
+    "restrict_to_mask": True,
+}
+
+
+def preload_core(so_path: str | os.PathLike | None) -> None:
+    """Pre-seed ``sys.modules['bioimage_cpp._core']`` from ``so_path``.
+
+    Must be called before ``bioimage_cpp`` is imported. A ``None`` path is a
+    no-op (the installed extension is used).
+    """
+    if so_path is None:
+        return
+    if "bioimage_cpp" in sys.modules or "bioimage_cpp._core" in sys.modules:
+        raise RuntimeError("preload_core must run before bioimage_cpp is imported")
+    resolved = str(Path(so_path).resolve())
+    spec = importlib.util.spec_from_file_location("bioimage_cpp._core", resolved)
+    if spec is None or spec.loader is None:
+        raise RuntimeError(f"cannot create an import spec for {resolved}")
+    module = importlib.util.module_from_spec(spec)
+    sys.modules["bioimage_cpp._core"] = module
+    spec.loader.exec_module(module)
+    import bioimage_cpp
+
+    # A pre-seeded sys.modules entry is not bound as a package attribute by
+    # the import system, so do it explicitly for `bioimage_cpp._core` users.
+    bioimage_cpp._core = module
+    loaded = sys.modules["bioimage_cpp._core"].__file__
+    if loaded != resolved:
+        raise RuntimeError(f"expected _core from {resolved}, got {loaded}")
+
+
+def pin_cpus(cpus: str | None) -> list[int]:
+    """Pin the current process to ``cpus`` ('36,37') before allocating data."""
+    if not cpus:
+        return sorted(os.sched_getaffinity(0))
+    selected = {int(c) for c in cpus.split(",") if c}
+    os.sched_setaffinity(0, selected)
+    return sorted(selected)
+
+
+def load_fixture(ndim: int, timeout: float = 60.0) -> tuple[np.ndarray, np.ndarray]:
+    """Registered fixture as ``(flow float32, mask bool)`` (see check_flow_density.py)."""
+    from bioimage_cpp._data import load_flow_data
+
+    dist, fg, _ = load_flow_data(ndim, timeout=timeout)
+    flow = np.ascontiguousarray(-dist, dtype=np.float32)
+    mask = np.ascontiguousarray(fg > 0.5)
+    return flow, mask
+
+
+def _stripes(n_orbiter_columns_of_4: int) -> tuple[np.ndarray, np.ndarray]:
+    from benchmark_interleave import SHAPE_3D, _stripe_flow
+
+    return _stripe_flow(n_orbiter_columns_of_4), np.ones(SHAPE_3D, dtype=bool)
+
+
+def _random3d(scale: float) -> tuple[np.ndarray, np.ndarray]:
+    from benchmark_interleave import SHAPE_3D
+
+    rng = np.random.default_rng(0)
+    flow = rng.normal(scale=scale, size=(3,) + SHAPE_3D).astype(np.float32)
+    return flow, np.ones(SHAPE_3D, dtype=bool)
+
+
+def _random2d(scale: float) -> tuple[np.ndarray, np.ndarray]:
+    from benchmark_midpoint_reuse import SHAPE_2D
+
+    rng = np.random.default_rng(0)
+    flow = rng.normal(scale=scale, size=(2,) + SHAPE_2D).astype(np.float32)
+    return flow, np.ones(SHAPE_2D, dtype=bool)
+
+
+# name -> zero-argument builder returning (flow, mask). The random cases use a
+# fresh default_rng(0) each, so they are deterministic but not bit-identical to
+# the sequentially drawn arrays inside benchmark_interleave/midpoint_reuse.
+CASES: dict[str, Callable[[], tuple[np.ndarray, np.ndarray]]] = {
+    "fixture3d": lambda: load_fixture(3),
+    "fixture2d": lambda: load_fixture(2),
+    "stripes1": lambda: _stripes(1),
+    "stripes3": lambda: _stripes(3),
+    "random1": lambda: _random3d(1.0),
+    "random10": lambda: _random3d(10.0),
+    "sweep5": lambda: _random3d(5.0),
+    "sweep20": lambda: _random3d(20.0),
+    "sweep40": lambda: _random3d(40.0),
+    "sweep2d10": lambda: _random2d(10.0),
+}
+
+FIXTURE_CASES = ("fixture3d", "fixture2d")
+ADVERSARIAL_CASES = (
+    "stripes1", "stripes3", "random1", "random10",
+    "sweep5", "sweep20", "sweep40", "sweep2d10",
+)
+
+
+def build_case(name: str) -> tuple[np.ndarray, np.ndarray]:
+    try:
+        builder = CASES[name]
+    except KeyError as error:
+        raise SystemExit(f"unknown case {name!r}; choose from {', '.join(CASES)}") from error
+    return builder()
+
+
+def make_runner(bare: bool, ndim: int) -> Callable:
+    """Return ``run(flow, mask_or_u8, threads, **params) -> density``.
+
+    ``bare`` calls the binding directly (mask must be uint8), bypassing the
+    Python wrapper's conversions; otherwise the public wrapper is used.
+    """
+    import bioimage_cpp as bic
+
+    if not bare:
+        def run(flow, mask, threads, **params):
+            return bic.flow.compute_flow_density(
+                flow, mask, sigma=None, number_of_threads=threads, **params
+            )
+        return run
+
+    core = sys.modules["bioimage_cpp._core"]
+    fn = getattr(core, f"_compute_flow_density_{ndim}d_float32")
+
+    def run_bare(flow, mask_u8, threads, **params):
+        p = {**DEFAULTS, **params}
+        return fn(
+            flow, mask_u8, p["n_iter"], p["dt"], p["tol"], p["method"],
+            p["restrict_to_mask"], threads,
+        )
+    return run_bare
+
+
+def add_param_args(parser) -> None:
+    parser.add_argument("--n-iter", type=int, default=None)
+    parser.add_argument("--dt", type=float, default=None)
+    parser.add_argument("--tol", type=float, default=None)
+    parser.add_argument("--method", choices=("euler", "rk2"), default=None)
+    import argparse
+
+    parser.add_argument(
+        "--restrict-to-mask", action=argparse.BooleanOptionalAction, default=None
+    )
+
+
+def params_from_args(args) -> dict:
+    params = {}
+    for name in ("n_iter", "dt", "tol", "method", "restrict_to_mask"):
+        value = getattr(args, name, None)
+        if value is not None:
+            params[name] = value
+    return params