From 54fd4a9040cf919a022406b421dffc36cfb3dbbf Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Wed, 2 Sep 2026 22:58:37 +0200 Subject: [PATCH 1/2] Further optimzations of flow tracking --- development/flow/FLOW_OPTIM.md | 644 +++++++++++++++++++ development/flow/PERFORMANCE_NOTES.md | 243 ++++++- development/flow/_flow_cases.py | 176 +++++ development/flow/differential_check.py | 198 ++++++ development/flow/paired_bench.py | 189 ++++++ development/flow/perf_kernel.py | 311 +++++++++ include/bioimage_cpp/detail/finite.hxx | 59 ++ include/bioimage_cpp/detail/force_inline.hxx | 16 + include/bioimage_cpp/flow/flow_density.hxx | 236 ++++++- src/bindings/flow.cxx | 22 +- src/bioimage_cpp/flow/_flow.py | 16 +- src/cpp/flow/flow_density_fma.cxx | 502 ++++++++++++++- tests/test_flow.py | 66 ++ 13 files changed, 2638 insertions(+), 40 deletions(-) create mode 100644 development/flow/FLOW_OPTIM.md create mode 100644 development/flow/_flow_cases.py create mode 100644 development/flow/differential_check.py create mode 100644 development/flow/paired_bench.py create mode 100644 development/flow/perf_kernel.py create mode 100644 include/bioimage_cpp/detail/finite.hxx create mode 100644 include/bioimage_cpp/detail/force_inline.hxx diff --git a/development/flow/FLOW_OPTIM.md b/development/flow/FLOW_OPTIM.md new file mode 100644 index 0000000..581b11a --- /dev/null +++ b/development/flow/FLOW_OPTIM.md @@ -0,0 +1,644 @@ +# Flow-density optimization investigation + +Date: 2026-07-12 + +## Executive summary + +This document records a new performance investigation of +`bioimage_cpp.flow.compute_flow_density`, building on the experiments already +documented in `development/flow/PERFORMANCE_NOTES.md`. + +The investigation found three complementary changes that are substantially +more promising than the SIMD, channel-layout, scatter, and active-list ideas +previously explored: + +1. Trace a particle through all of its integration steps inside one parallel + worker invocation, instead of launching and joining workers once per global + integration step. +2. Replace the generic `2^D` corner-table construction with explicit bilinear + and trilinear interpolation. +3. Convert already-clipped, nonnegative sampling coordinates to integers by + truncation instead of calling `std::floor`. + +In a temporary implementation, the combined change reduced the registered 3D +fixture runtime from 7.96 s to 3.87 s with one thread and from 2.05 s to 0.90 s +with eight threads on an AMD EPYC 7513. This is a reduction of approximately +51% and 56%, respectively. The registered 2D fixture improved from 1.09 s to +0.54 s with one thread, a reduction of approximately 51%. + +The best temporary implementation passed all 17 tests in `tests/test_flow.py`, +retained the same stored-reference accuracy metrics, and produced a density +array that was bitwise identical to the current implementation on the full 3D +fixture. No algorithm changes from this investigation have been applied to the +repository; the experiments were made in a copy under `/tmp`. + +## Scope and constraints + +The goals of this investigation were: + +- Take all existing flow benchmarks and optimization notes into account. +- Identify remaining performance opportunities in the current implementation. +- Validate ideas by measurement rather than relying on instruction-count or + memory-bandwidth intuition alone. +- Keep experiments outside the repository until the evidence justified a + concrete recommendation. +- Preserve the existing dependency-free, portable C++20 implementation and + deterministic Python API. + +The investigation did not change public parameters, defaults, numerical data +types, the channel-first flow layout, the integration methods, or the density +scatter semantics. + +## Current implementation + +The current implementation is in +`include/bioimage_cpp/flow/flow_density.hxx`. Its main phases are: + +1. Zero the output and collect the coordinates of foreground voxels into an + array of particle positions. +2. For every global integration iteration: + - call `detail::parallel_for_chunks` over all particle positions; + - skip particles whose `alive` byte is zero; + - clip the current position; + - build a shared table of corner offsets and weights; + - sample all flow channels from that table; + - optionally take a second sample for RK2; + - test convergence and the foreground mask; + - update the position or mark the particle dead; + - join all worker threads; + - scan the complete `alive` array to determine whether tracing can stop. +3. Scatter the final particle positions into the density array. +4. Zero density outside the foreground mask. + +The existing profiler shows that approximately 99% of the runtime is in the +integration loop. Initialization, scatter, and mask-zeroing are each well below +1%, so this investigation continued to focus exclusively on particle tracing. + +### Generic interpolation work + +For every flow sample, `compute_corners` currently: + +- calls `std::floor` once per spatial axis; +- computes a fractional coordinate once per axis; +- visits all `2^D` corners; +- for every corner, revisits all `D` axes; +- selects the lower or upper coordinate; +- multiplies the corresponding weight; +- clamps the selected coordinate; +- multiplies it by the grid stride and accumulates the flat offset. + +The resulting offsets and weights are shared across flow channels, which was a +useful improvement over recomputing them independently for every channel. +Nevertheless, the generic construction performs much more coordinate, +branching, and weight-product work than fixed 2D and 3D interpolation require. + +### Parallel-loop structure + +`detail::parallel_for_chunks` creates `n_threads - 1` `std::thread` objects and +joins them before returning. The current integration loop calls it once per +global integration step. With the default `n_iter=50`, this can create and join +the worker set up to 50 times per function call. + +The iteration-major order also requires a global byte per particle to retain +the alive state, scans dead particles in later iterations, and performs a +sequential complete-array alive reduction after every step when convergence or +mask restriction is enabled. + +Every trajectory is independent until the final sequential density scatter. +There is therefore no algorithmic requirement for particles to advance in +lockstep. + +## Prior optimization work considered + +The existing `development/flow/PERFORMANCE_NOTES.md` was treated as the starting +point, not as a list of experiments to repeat. The following successful changes +were retained in all experiments: + +- Multithreaded particle tracing through `detail::parallel_for_chunks`. +- Per-particle convergence detection via `tol`. +- Optional termination when a proposed endpoint leaves the foreground mask. +- RK2 integration, which enables useful accuracy/runtime defaults. +- Sharing a corner table across the flow channels. +- Hoisting channel pointers, strides, and grid bounds. +- Existing phase profiling through `BIOIMAGE_PROFILE`. + +The following previously rejected approaches were not reintroduced: + +- Parallel or per-thread density scatter buffers, because scatter is below 1% + of runtime. +- Fusing the initialization passes, because initialization is approximately + 0.4% of runtime. +- A globally compacted active-particle list. +- Compiler-generated SIMD through target clones. +- Handwritten AVX2 for the existing corner sum. +- Channel-last/interleaved flow storage, including padded SIMD-friendly + storage. +- Half-precision flow storage. +- A standalone structure-of-arrays position layout whose main purpose would be + SIMD. + +The prior AVX2 and interleaved-layout experiment is particularly important. It +produced the intended packed FMA instructions but did not improve runtime. This +showed that reducing only interpolation arithmetic, while retaining the same +corner-address and load behavior, was insufficient on the earlier Intel test +machine. + +## Experimental environment + +All new source experiments were performed in a copied repository at: + +```text +/tmp/bic-flowopt.5WuTmL/repo +``` + +An isolated virtual environment and persistent build directory were also kept +under that temporary directory. The production repository remained unchanged +during the experiments. + +Hardware and runtime environment: + +- CPU: AMD EPYC 7513 32-Core Processor, Zen 3 generation. +- Host topology: two sockets, 32 physical cores per socket, two hardware + threads per core. +- CPUs available to the process: `32-35,96-99`, corresponding to eight logical + CPUs. +- Compiler optimization: the project's normal non-Debug `-O3` build. +- Python: 3.14. +- Performance counters: Linux `perf` was available. +- Registered fixture shape: `(48, 512, 512)` for 3D, with 729,236 foreground + particles. +- Registered fixture shape: `(520, 704)` for 2D, with 151,019 foreground + particles. +- Flow defaults: `n_iter=50`, `dt=0.2`, `tol=0.005`, `method="rk2"`, + `restrict_to_mask=True`, and no density smoothing. + +This is a different machine from the Intel Tiger Lake laptop used for the +existing performance notes. Absolute timings should therefore not be compared +between the two documents. Relative timings within each experiment are the +relevant measurements. + +The baseline was built from an untouched temporary copy and reproduced the +timings of the already installed current implementation, which reduced the +risk that build-directory or compiler differences biased the comparison. + +## Baseline measurements + +Registered fixture timings for the current implementation were: + +| Fixture | Threads | Median | Minimum | Repeats | +|---|---:|---:|---:|---:| +| 3D | 1 | 7.9621 s | 7.9464 s | 3 | +| 3D | 8 | 2.0529 s | 2.0472 s | 3 | +| 2D | 1 | 1.0876 s | 1.0845 s | 5 | + +The initial hardware-counter run on the current implementation reported high +instruction throughput and a low branch-miss rate: + +| Counter | Current implementation | +|---|---:| +| Cycles | 16.394 billion | +| Instructions | 51.949 billion | +| Instructions per cycle | 3.17 | +| Branches | 8.945 billion | +| Branch misses | 44.45 million, 0.50% | +| Cache references | 405.3 million | +| Cache misses | 63.35 million | + +These were whole-process `perf stat -r 3` measurements of the validation script, +not counters around only the C++ kernel. They include data loading and accuracy +calculation, and the events were multiplexed at approximately 83%. They are +still useful for comparing large changes made under the same conditions. + +The low branch-miss rate argues against branch prediction hints as a useful +primary optimization. The high total instruction count suggested that removing +work could be valuable even though the existing roofline analysis correctly +showed that the kernel does not saturate streaming memory bandwidth. + +## Successful experiment 1: particle-major traversal + +### Idea + +Move the integration loop inside the per-particle loop: + +```text +parallel over contiguous particle chunks once: + for each particle in the chunk: + for up to n_iter steps: + sample, test, and update this particle + break when it converges or leaves the mask +``` + +This is valid because particles do not interact during tracing. Only the final +density scatter combines their results, and that scatter remains sequential and +deterministic. + +### Work removed + +- Up to 49 of the 50 worker-thread creation/join cycles. +- The global `alive` vector. +- Alive-byte loads and dead-particle branches in later global iterations. +- The complete sequential alive reduction after each global iteration. +- Repeated loads and stores of a particle position between global iterations; + the compiler has a better opportunity to keep one trajectory in registers. + +### Measured result + +| Fixture | Threads | Current | Particle-major | Improvement | +|---|---:|---:|---:|---:| +| 3D | 1 | 7.9621 s | 6.8670 s | 13.8% | +| 3D | 8 | 2.0529 s | 1.8783 s | 8.5% | + +The one-thread improvement proves that repeated worker creation is not the only +mechanism. Removing alive-state traffic and keeping a trajectory locally are +also important. + +### Risks + +The current iteration-major traversal gives each thread the same static range +of particles on every iteration. Particle-major traversal also assigns static +contiguous ranges, but each thread now owns the entire remaining lifetime of +its range. If slow-converging trajectories are spatially clustered, threads may +finish at different times and reduce parallel utilization. + +The registered 3D fixture still improved at eight threads, but an implementation +should add a deliberately imbalanced benchmark or test before assuming this is +universal. A dynamic scheduler should not be introduced casually: it would add +new infrastructure, scheduling overhead, and portability complexity. First +measure whether real workloads exhibit enough imbalance to need it. + +## Successful experiment 2: direct bilinear/trilinear interpolation + +### Idea + +Replace the generic corner table with a small sampling description containing +the lower coordinate, upper coordinate, and fractional coordinate for each +axis. Sample each channel through nested interpolation. + +Conceptually, 3D sampling becomes: + +1. Interpolate along x for each of the four `(z, y)` pairs. +2. Interpolate the resulting values along y for each z plane. +3. Interpolate the two plane values along z. + +The 2D case performs two x interpolations followed by one y interpolation. + +At an upper grid boundary, the lower and upper coordinate are identical, which +preserves the current nearest-boundary behavior. + +### Work removed + +- The generic loop over every corner and every axis. +- Repeated lower/upper coordinate selection for every corner. +- Repeated per-corner bounds checks. +- Construction and storage of all product weights. +- Construction and storage of a full corner table when only fixed 2D or 3D + interpolation is supported. + +The flow loads themselves remain channel-first and use the existing data. This +is fundamentally different from the rejected SIMD experiment: the new approach +reduces address-generation, coordinate, branch, and weight-construction work, +not just the arithmetic instructions used to sum already-constructed corners. + +### Measured result + +| Fixture | Threads | Current | Direct interpolation | Improvement | +|---|---:|---:|---:|---:| +| 3D | 1 | 7.9621 s | 6.7392 s | 15.4% | +| 3D | 8 | 2.0529 s | 1.5889 s | 22.6% | +| 2D | 1 | 1.0876 s | 0.8853 s | 18.6% | + +Direct interpolation and particle-major traversal composed well: + +| Fixture | Threads | Current | Combined | Improvement | +|---|---:|---:|---:|---:| +| 3D | 1 | 7.9621 s | 5.0713 s | 36.3% | +| 3D | 8 | 2.0529 s | 1.2550 s | 38.9% | +| 2D | 1 | 1.0876 s | 0.7447 s | 31.5% | + +### Numerical considerations + +Nested interpolation changes the association order of floating-point +operations relative to multiplying eight precomputed weights and summing them. +The mathematical interpolation is the same, but intermediate float32 rounding +is not guaranteed to be bitwise identical for every possible flow field. + +On the registered fixtures, all reported accuracy metrics were unchanged. The +final best 3D density was also bitwise identical to the current density. +Nevertheless, an implementation should include targeted differential tests of +the sampler itself and end-to-end randomized tests, rather than relying on the +integer-valued density scatter to hide small trajectory differences. + +## Successful experiment 3: truncate clipped coordinates + +### Idea + +Every position passed to interpolation is clipped to the closed grid domain: + +- the current particle position is clipped before its first sample; +- the RK2 midpoint is clipped before its second sample. + +Sampling coordinates are therefore finite and nonnegative. For a nonnegative +float, conversion to an integer by truncation has the same result as +`floor(position)`. Replacing: + +```cpp +static_cast(std::floor(position[axis])) +``` + +with: + +```cpp +static_cast(position[axis]) +``` + +is exact under this precondition. + +This optimization should remain local to a helper whose contract requires +clipped coordinates. It should not silently change a general-purpose sampling +helper that might later be called with negative positions. + +### Measured result + +Adding truncation to the interpolation-only implementation reduced the 3D +one-thread runtime from 6.7392 s to 5.5406 s, an additional 17.8% improvement. + +Adding it to the combined particle-major and direct-interpolation implementation +produced the best result: + +| Fixture | Threads | Current | Best temporary variant | Improvement | +|---|---:|---:|---:|---:| +| 3D | 1 | 7.9621 s | 3.8657 s | 51.5% | +| 3D | 8 | 2.0529 s | 0.9046 s | 55.9% | +| 2D | 1 | 1.0876 s | 0.5351 s | 50.8% | + +The size of this effect was unexpected. It may depend strongly on compiler and +microarchitecture, so the exact gain must be checked on Intel x86-64, arm64, +and the wheel build toolchains. The semantic equivalence does not depend on the +CPU as long as the clipped-coordinate precondition holds. + +## Counter evidence for the combined change + +The best temporary implementation was measured with the same whole-process +`perf stat -r 3` command as the current implementation: + +| Counter | Current | Best temporary variant | Change | +|---|---:|---:|---:| +| Cycles | 16.394 B | 10.269 B | -37.4% | +| Instructions | 51.949 B | 25.391 B | -51.1% | +| Instructions per cycle | 3.17 | 2.48 | lower, but much less total work | +| Branches | 8.945 B | 2.248 B | -74.9% | +| Branch misses | 44.45 M | 43.79 M | approximately unchanged absolute count | +| Cache references | 405.3 M | 275.5 M | -32.0% | +| Cache misses | 63.35 M | 27.64 M | -56.4% | + +Because these counters include Python data loading and result checking, they +understate the proportional reduction inside the C++ kernel. They nevertheless +confirm the main mechanism: the improvement comes from executing much less +work. It is not an artifact of a different flow layout, reduced precision, +non-deterministic scatter, or architecture-specific packed SIMD. + +The reduced IPC does not indicate a regression. The current implementation +has abundant independent generic corner arithmetic that can retire at high +throughput. The optimized version removes much of that arithmetic, leaving a +larger fraction of the remaining cycles exposed to the unavoidable dependent +flow sampling. Total cycles and wall time still fall substantially. + +## Correctness validation + +Every major successful variant was checked with: + +```bash +python -m pytest tests/test_flow.py -q +``` + +All 17 tests passed. These tests cover: + +- 2D and 3D tracing. +- Euler and RK2 paths. +- Mask restriction enabled and disabled. +- Convergence enabled and disabled. +- Single-threaded and multithreaded equality. +- Non-contiguous Python inputs. +- Degenerate iteration counts and zero flows. +- Invalid inputs. +- Density smoothing integration. + +The registered 2D and 3D validation scripts continued to pass their stored +reference gates. For the final 3D variant: + +- relative difference against the stored reference remained 0.0470; +- Pearson correlation remained 0.9690; +- the density sum remained exactly 729,236 particles; +- the full density array was bitwise equal to the current implementation; +- zero of 12,582,912 voxels differed. + +Bitwise equality on this fixture is strong evidence, but it is not a proof that +nested interpolation will produce identical trajectories for every input. + +## New experiments that did not help + +### Conditional removal of start-of-step clipping + +With `restrict_to_mask=True`, every committed endpoint is in bounds, so clipping +the current position at the start of the next iteration is logically redundant. +The experiment wrapped clipping in a runtime `if (!restrict_to_mask)` condition. + +This regressed the best 3D one-thread result from approximately 3.87 s to +4.97 s. The likely cause is worse compiler optimization or hot-loop layout from +the additional invariant runtime branch. The experiment was reverted. + +If this is revisited, it should only be through compile-time specialization of +the mask-restricted and unrestricted kernels, followed by careful code-size and +cross-platform measurements. The current evidence does not justify that added +complexity. + +### Explicit displacement reuse + +The current code calculates `dt * step[axis]` once for the convergence maximum +and again when forming the proposed endpoint. A temporary implementation stored +the displacement in an array and reused it. + +The result was neutral: 3.8871 s versus 3.8657 s for the best implementation, +within normal run variation. The compiler already appears capable of handling +this arithmetic efficiently, or the extra temporary storage cancels the saved +multiplications. + +### Small-block iteration-major traversal + +Particle-major traversal introduces a dependent chain: the next flow address +for one particle is unknown until the current sample updates its position. A +possible latency-hiding strategy is to process a small block of neighboring +particles in iteration-major order inside each persistent worker. + +A block size of 16 was tested. It retained one parallel invocation but used a +small local alive table and interleaved 16 trajectories. It regressed the best +implementation: + +| Threads | Particle-major best | Block size 16 | Regression | +|---:|---:|---:|---:| +| 1 | 3.8657 s | 5.1927 s | 34.3% | +| 8 | 0.9046 s | 1.2043 s | 33.1% | + +The benefit of keeping one trajectory in registers outweighed any additional +memory-level parallelism or local spatial reuse. Other block sizes were not +tested because the result was not marginal and the simpler particle-major +design already performed well. + +## Revised performance diagnosis + +The existing notes concluded that the hot loop was not limited by streaming +memory bandwidth and that flow-load latency was likely important. That remains +compatible with the new evidence, but the earlier conclusion was too narrow if +interpreted as meaning that only prefetching could help. + +The current generic interpolation loop executes a very large amount of +coordinate, bounds, offset, weight, alive-state, and loop-control work around +the actual loads. It can sustain high IPC, but it also retires approximately +twice as many instructions as the best temporary implementation in the +whole-process comparison. + +The rejected SIMD experiment reduced arithmetic for an already-constructed +corner set without removing most address-generation and control work. The new +direct interpolation changes the algorithmic formulation of that work. The +particle-major traversal independently removes synchronization and global +state traffic. The truncation change removes an expensive coordinate operation +under a valid local precondition. + +After these changes, load latency may become a larger fraction of the remaining +runtime. Software prefetching is therefore still a possible later experiment, +but it is lower priority than implementing and validating the demonstrated +instruction-count reductions. The failed 16-particle interleaving experiment +also shows that adding machinery solely to expose more independent loads can +easily lose more than it gains. + +## Recommended implementation plan + +### 1. Introduce a narrow direct sampling helper + +Add a small internal sampling description or explicit 2D/3D helpers in +`bioimage_cpp::flow::detail`. Keep the API focused on clipped coordinates and +C-contiguous grid strides. Avoid introducing a generic interpolation framework +or new dependency. + +Recommended properties: + +- Separate, readable 2D and 3D code paths selected with `if constexpr` or two + small overloads. +- Compute lower, upper, and fractional coordinates once per sample position. +- Use truncation only after documenting or asserting the nonnegative clipped + precondition. +- Load each corner value into a named local before interpolation, both for + readability and to make single-load intent obvious. +- Preserve the current duplicated-coordinate behavior at upper boundaries. +- Reuse the same sampling description across all flow channels. + +The temporary proof of concept prioritized experimental speed and should be +cleaned up before landing, particularly by shortening long interpolation +expressions and documenting boundary behavior. + +### 2. Change to particle-major tracing + +Move the per-particle iteration loop inside one call to +`detail::parallel_for_chunks`. Preserve contiguous static chunks and the +existing final sequential scatter. This remains compatible with the project's +single threading primitive and deterministic behavior. + +The implementation can remove: + +- the global `alive` vector; +- the outer global iteration loop; +- the sequential complete alive scan. + +For each particle, convergence or mask exit becomes a simple `break` from its +local integration loop. + +### 3. Keep existing profiling scopes and add temporary subphase evidence + +The existing `iter_loop` profile scope should remain. During implementation, +temporary profiling or benchmark variants can distinguish: + +- sampling-coordinate preparation; +- first flow sample; +- RK2 midpoint and second sample; +- convergence/mask/update logic. + +Use the existing profiling macros rather than ad hoc timers. Fine-grained +scopes inside every particle step may themselves be intrusive, so use them only +in a dedicated profiling build or compare separately compiled variants. + +### 4. Add differential tests + +In addition to the existing suite, add tests for: + +- direct sampling at integer coordinates; +- fractional coordinates in 2D and 3D; +- the last coordinate on every axis; +- shapes containing an axis of length one; +- random positions near voxel-rounding boundaries; +- random finite flow fields comparing old and new implementations during + development; +- both Euler and RK2; +- `restrict_to_mask` enabled and disabled; +- `tol=0` and positive convergence tolerance; +- masks designed so different contiguous particle ranges require very + different iteration counts; +- equality across thread counts. + +The old sampler can be kept only in a development comparison or temporary test +helper while validating the change; it should not remain as duplicate +production infrastructure after confidence is established. + +### 5. Benchmark multiple architectures + +At minimum, repeat the registered benchmark matrix on: + +- an Intel x86-64 machine representative of the earlier Tiger Lake results; +- the AMD Zen 3 host used here; +- macOS arm64 with AppleClang; +- one Linux arm64 wheel environment when available; +- Windows x86-64 with MSVC. + +The particle-major and direct-interpolation changes are portable C++20 and are +expected to generalize. The magnitude of the truncation gain is the most likely +to vary with compiler and CPU. + +Recommended measurements: + +```bash +python development/flow/check_flow_density.py --dim both --repeats 5 --threads 1 +python development/flow/check_flow_density.py --dim 3 --repeats 5 --threads 2 +python development/flow/check_flow_density.py --dim 3 --repeats 5 --threads 4 +python development/flow/check_flow_density.py --dim 3 --repeats 5 --threads 8 +python -m pytest tests/test_flow.py -q +``` + +Where available, collect `perf stat` counters around a script that loads the +fixture once and executes several warm kernel calls. This would isolate the C++ +kernel better than the whole-process counters used in this investigation. + +## Priority ranking + +1. **Implement direct 2D/3D interpolation plus clipped-coordinate truncation.** + This removes the largest clearly identified body of unnecessary work and + has no threading/load-balance risk. +2. **Implement particle-major traversal.** It provides a separate substantial + gain and composes exceptionally well with direct interpolation. Validate + load balance on heterogeneous trajectories. +3. **Run cross-architecture benchmarks and hardware counters.** The combined + gain is large enough that it should survive ordinary noise, but compiler and + architecture coverage are essential for a wheel-oriented project. +4. **Only then reconsider prefetching or compile-time mode specialization.** + These are more complex and currently lack positive experimental evidence. + +## Conclusion + +The strongest remaining opportunity is not more SIMD or a new memory layout. +It is simplifying the scalar algorithm and matching the loop structure to the +independence of particle trajectories. + +The combined temporary implementation approximately halved 2D and 3D runtime, +improved both single-threaded and eight-threaded execution, retained all +existing test behavior, and produced a bitwise-identical full 3D density result +for the registered fixture. The evidence is strong enough to justify an +implementation pass, provided it includes targeted numerical tests, +load-balance checks, and validation on the project's major wheel platforms. diff --git a/development/flow/PERFORMANCE_NOTES.md b/development/flow/PERFORMANCE_NOTES.md index 9781f27..4416608 100644 --- a/development/flow/PERFORMANCE_NOTES.md +++ b/development/flow/PERFORMANCE_NOTES.md @@ -12,7 +12,7 @@ Hardware on which these numbers were taken: **11th Gen Intel Core i7-1185G7 @ 3.00 GHz** (Tiger Lake, 4 physical cores / 8 SMT threads, 12 MB L3, AVX2 + AVX-512 capable). The 2026-07-12 pass below was cross-checked against GPT 5.6-Sol's independent investigation, which measured the same changes on an -AMD EPYC 7513 (Zen 3); see `FLOW_OPTIM.md` in the repo root. Absolute timings +AMD EPYC 7513 (Zen 3); see `development/flow/FLOW_OPTIM.md`. Absolute timings differ between the hosts; the relative gains agree closely across both. ## Headline numbers @@ -40,7 +40,14 @@ The rewrite produces a **byte-for-byte identical** density on both fixtures: pre-rewrite build. It was also verified bitwise identical across 1056 randomized small/edge-case inputs (see "Correctness" below). -Two eras of work are recorded here: +**2026-09-02 update (Zen 3, see "Packed-channel FMA tracer" below):** the +default path on AVX+FMA x86-64 CPUs now uses a packed-channel vectorized +tracer. On the same machine and flags, the registered fixtures run +**1.052 s -> 0.637 s (3D, 1T)** and **0.110 s -> 0.069 s (2D, 1T)**, +-39 % and -37 %, with identical densities; 2T and 4T gain 30-42 %, adversarial +inputs 17-44 %. + +Three eras of work are recorded here: 1. **Threading + early-exit era** (earlier). Threading was the big lever; the algorithmic changes capped work so the threaded path had something to scale. @@ -48,6 +55,11 @@ Two eras of work are recorded here: multi-threaded runtime by cutting per-step instruction count and removing the per-iteration synchronization, on the existing channel-first layout with no SIMD, no relayout, and no precision change. +3. **Packed-channel FMA tracer** (2026-09-02). Another -37..-44 % by repacking + the flow field channel-last once and tracing with 4-lane packed arithmetic + inside the runtime-dispatched FMA translation unit, after a per-step counter + budget showed the scalar kernel retired ~330 instructions per step and spent + a third of its cycles with the FP register file full. ## Changes that landed @@ -266,7 +278,8 @@ to `trace_particle`); the RK2 selectors use K=3, with a plain `trace_particle` remainder loop. Measured with tightly paired `.so`-swap ABA runs (K1 = per-particle loop, -`paired_bench.py`, min-of-3 per invocation, Tiger Lake). Two flag regimes, +a local predecessor of today's `development/flow/paired_bench.py`, min-of-3 per +invocation, Tiger Lake). Two flag regimes, because the conda toolchain exports default `CXXFLAGS` (`-march=nocona -mtune=haswell -ftree-vectorize ...`) while wheel builds get plain flags — see the pitfall below: @@ -309,6 +322,190 @@ conda) get the plain flags, so plain-flag numbers are the shipping-relevant ones. When A/B-benchmarking local builds: `echo $CXXFLAGS` first, and keep the regime identical on both sides. +## Packed-channel FMA tracer (2026-09-02) + +Third optimization round, on an **AMD EPYC 7513 (Zen 3)** node under a SLURM +allocation of 2 physical cores / 4 logical CPUs (CPUs 36,37 + SMT siblings +100,101; NUMA node 1), governor `performance`, THP `always`, `nmi_watchdog=1` +(five programmable counters per `perf` group). GCC 14.3 (conda-forge), Python +3.14.5, NumPy 2.4.6, **plain flags** (`CXXFLAGS= CFLAGS=`; the conda toolchain's +`-march=nocona ... -O2` defaults were cleared), commit `492ff9b` as baseline. +All A/B numbers come from `development/flow/paired_bench.py` (fresh pinned +subprocess per measurement, A B B A per repeat, prebuilt `.so` files swapped +through `sys.modules`, never rebuilt between sides) and from the kernel-only +counter harness `development/flow/perf_kernel.py`. This machine cannot measure +beyond two physical cores; 8T-class numbers remain to be taken elsewhere. + +### What changed + +The default path (RK2, `tol > 0`, `restrict_to_mask`) on CPUs with AVX+FMA now +runs a **packed-channel kernel** in `src/cpp/flow/flow_density_fma.cxx`: + +- The flow field is repacked once (parallel over planes/rows) into channel-last + storage, D floats per voxel, on a grid padded by one replicated voxel per + axis (+2.5 % voxels on the 3D fixture). Every trilinear corner is then a + single 16-byte load carrying all channels, and the padding makes the upper + corner index always valid, so the scalar kernel's nearest-boundary rule + needs no clamp. Cost: a transient buffer the size of the flow input + (151 MB on the 3D fixture) and a pack pass of ~16 ms at 1T / ~10 ms at 2T. +- The particle position is one 4-lane register `[z, y, x, 0]`. Truncation, + fraction, midpoint, clip, convergence test (one compare + `movemask`), + bounds test and the rounding for the mask lookup are packed; the flat + offsets are formed in general-purpose registers from three extracted lanes. +- 3D interpolation keeps the nested-lerp association order of the scalar + kernel (fused `lo + f*(hi-lo)`, fused midpoint, unfused position update), + so the traced positions match the scalar-FMA kernel bit for bit on the + fixtures. 2D uses a four-weight sum (`mul -> fma -> add`), which shortens + the post-load chain from 16 to 11 cycles; its association order differs, + and densities were still bitwise identical on both fixtures and all 400 + randomized cases. +- Lockstep interleave of **K = 4** trajectories per group (K = 3 for the + scalar kernel), lanes addressed with compile-time indices so the state + stays in registers. +- The scalar-FMA instantiation of the header kernel remains as fallback when + the packed buffer cannot be allocated or its offsets exceed int32, and the + portable SSE2 kernel is untouched for selectors 0-6 and non-AVX CPUs. +- New runtime opt-out `BIOIMAGE_CPP_FLOW_FORCE_SCALAR=1` (mirrors the filters + module) and `_core._flow_trace_backend()` returning `"fma"` or `"scalar"`, + with a parity test in `tests/test_flow.py`. +- The binding's finiteness scan over the flow buffer (previously a serial + `std::isfinite` loop holding the GIL, ~30 ms of the 3D fixture) now runs + inside the GIL release through `detail/finite.hxx::all_finite` (branch-free + per-block exponent test, parallel with the caller's thread count). +- Hot helpers (header `sample_flow`, `round_to_flat_index`, + `position_is_in_mask`, `trace_particle`, `trace_particle_block`; the TU's + samplers and step) carry `BIOIMAGE_FORCE_INLINE` + (`include/bioimage_cpp/detail/force_inline.hxx`). See the codegen trap below. + +### Why: per-step diagnosis of the previous kernel + +The 2026-07-13 gate had closed prefetching on the grounds that misses were rare; +it had not converted the counters into a per-step budget. A profile-build step +counter (`TraceStats`, printed as `[bioimage flow trace]`) now gives: + +``` +3D fixture: 33,959,237 executed steps for 729,236 particles (46.6 per particle) + exits: converged 4.4 %, left mask 5.6 %, hit n_iter 90.0 % +2D fixture: 7,230,607 executed steps for 151,019 particles (47.9 per particle) + exits: converged 6.1 %, left mask 0.3 %, hit n_iter 93.6 % +``` + +Baseline kernel, 3D fixture, 1T, per executed step (perf_kernel, 8 warm calls): + +``` +cycles 114 instructions 326 loads 95 (49 essential) IPC 2.87 +L1d misses 1.4/step (1.5 %) L2 misses 0.03/step DRAM fills 0.007/step +dTLB misses ~0.0007/step (THP effective) +FP-register-file dispatch stalls: 40 cycles/step (35 % of cycles) +``` + +So the kernel was neither memory- nor branch-bound; it retired ~330 +instructions per step, half of them address/weight bookkeeping in the FP +register domain, and spent a third of its cycles with the FP register file full. +Static disassembly of the FMA TU confirmed ~330-360 instructions and 42-47 +stack operands per lane-step. That is what the packed kernel attacks: fewer +instructions, fewer FP-domain micro-ops, shorter chain. + +### Results + +Paired A/B (`paired_bench.py`, ABBA, best-of-subprocess minima, `.so` swap), +previous kernel (`492ff9b`) vs packed kernel, plain flags, Zen 3: + +``` +case threads cpus previous packed change noise parity +fixture3d 1 36 1.0515 s 0.6372 s -39.4 % 0.6 % identical +fixture2d 1 36 0.1099 s 0.0689 s -37.3 % 1.8 % identical +fixture3d 2 phys 36,37 0.5520 s 0.3388 s -38.6 % 0.4 % identical +fixture2d 2 phys 36,37 0.0561 s 0.0360 s -35.9 % 1.3 % identical +fixture3d 2 SMT 36,100 0.8489 s 0.4914 s -42.1 % 0.5 % identical +fixture3d 4 36,37,100,101 0.4576 s 0.2675 s -41.6 % 1.9 % identical +fixture2d 4 36,37,100,101 0.0441 s 0.0309 s -29.9 % 1.0 % identical +stripes 1/4 1 36 0.5817 s 0.4821 s -17.1 % 1.3 % identical +stripes 3/4 1 36 1.2842 s 0.7707 s -40.0 % 2.8 % identical +random s=1 1 36 1.5414 s 0.8931 s -42.1 % 3.0 % identical +random s=10 1 36 1.6568 s 0.9259 s -44.1 % 0.9 % identical +random s=5 1 36 1.1499 s 0.6816 s -40.7 % 0.3 % identical +random s=20 1 36 1.5028 s 0.8957 s -40.4 % 2.7 % identical +random s=40 1 36 0.8152 s 0.5425 s -33.5 % 3.2 % identical +2D random 10 1 36 0.1727 s 0.1103 s -36.1 % 0.4 % differs (2D tree) +``` + +The base-vs-base calibration of the same harness reads `noise` on both fixtures +(|dMin| <= 0.12 %). The stripes 1/4 case (three of four lanes converge on the +first step) gains least, as expected for a lockstep group; K = 4 still beat +K = 3 there by 3.6 %. + +Accuracy gate (`check_flow_density.py --dim both --repeats 3`): unchanged, +3D `rel_diff = 0.0470`, `pearson = 0.9690`; 2D `rel_diff = 0.0201`, +`pearson = 0.9975`; particle sums 729,236 / 151,019; densities bitwise equal to +the previous kernel on both fixtures. `differential_check.py` (400 randomized +cases: length-1 axes, tiny grids, Euler/RK2, `tol`/`dt`/`n_iter` variants, mask +on/off, 1/4 threads) is 400/400 bitwise identical previous-vs-packed, and +399/400 packed-vs-forced-scalar (one scale-50 random flow: 0.63 % of voxels +differ by <= 2 counts, sums equal; contraction noise on a chaotic field). Full +suite: 1457 passed, 7 skipped. + +Per executed step, 3D fixture, 1T (2D in parentheses): + +``` + previous packed +cycles/step 114 (56) 67 (35) +instructions/step 326 (145) 133 (81) +loads/step 95 (39) 38 (16) +IPC 2.87 (2.60) 1.98 (2.33) +L1d misses/step 1.4 0.34 +FP-RF dispatch stall 40 cyc (35 %) 36 cyc (54 %) (2D: 20 cyc, 57 %) +flops/step (FMA = 2) n/a 196 (84) +``` + +Static K = 4 lane body (final build): 3D 122 instructions per lane-step +(110 GPR, 107 FP ALU, 74 loads incl. 8 constant operands, 59 FMA, 24 shuffles, +24 vector-int, 20 cvt per group of four; 2 stack operands per lane-step); +2D 78 per lane-step. The remaining cost is the dependent chain per step +(~90 cycles in 3D: two samples of cvt -> extract -> imul/add -> load -> three +lerp levels, plus midpoint and update) against the reorder window and the FP +register file; four lanes recover ~1.35 lane-steps of overlap. + +Thread balance (profile build, `[bioimage flow trace]` per-thread lines): +2T physical on the 3D fixture 0.9 % time / 1.1 % steps imbalance; stripes 1/4 +0.1 %. Static contiguous chunks remain adequate. Phase split at 1T on the 3D +fixture: pack 16.6 ms (2.6 %), tracing 96.5 %, init 1.7 %, scatter 0.8 %, +mask zero 1.0 %. + +Non-kernel: the parallel, GIL-free finiteness check cuts the `n_iter=0` call +from ~40 ms to a few ms on the 3D fixture (was 3.8 % at 1T, 7 % at 2T of the +previous kernel's runtime, growing with thread count). + +### Codegen trap: translation-unit growth disables inlining + +While the FMA TU held many template instantiations (K and layout variants for +the sweep), GCC 14 stopped inlining `step_packed`, `sample_packed` and even the +header's `sample_flow<3>`: the scalar-FMA fallback slowed from 1.05 s to 1.59 s +and the packed kernel from 0.65 s to 1.24 s, with out-of-line calls in the hot +loop. An out-of-line `sample_flow<3>` compiled with AVX is also exactly the +COMDAT/ODR hazard the dispatch design guards against. `BIOIMAGE_FORCE_INLINE` +on the per-step helpers removes the dependence on the unit-growth budget in both +translation units. + +### Rejected in this round + +- **Vector-integer address math** (`pmulld` + horizontal sums for the flat + offset, `pminsd`-clamped upper indices): 0.720 s vs 0.645 s for the + integer-register version at K=3 in 3D, and a 2D regression (0.113 s vs the + scalar 0.110 s). Vector-integer ops occupy the FP register file, which was + the bottleneck. +- **16-byte voxels (zero fourth lane, aligned loads, no lane clean-up)**: + −2 % in 3D for +33 % transient memory; split-line loads (3.3 per step, ~20 % + of corner loads) are not a measurable cost. Not worth the memory. +- **Weight-tree interpolation in 3D**: within noise (−1 %); the 3D loop is + bound by the reorder window (~150 micro-ops per lane-step against a + ~90-cycle chain), so the extra weight micro-ops offset the shorter chain. + Kept for 2D only, where it measured −12 %. +- **K = 6 lanes**: parity with K = 4 in both dimensions; K = 2 and K = 1 lose + (3D: 0.652 s / 1.19 s vs 0.645 s). +- **Padding away the boundary clamps**: only ~1 % (integer clamps were not on + the critical path), kept because it simplifies the sampler. + ## Correctness - The current suite reports 1057 passed / 8 skipped. The 17 @@ -440,11 +637,15 @@ Kept so these avenues are not blindly re-attempted. ## Optimization status and future validation -The 2026-07-13 Zen 3 counter run completed the final `perf` gate and did not -support another kernel implementation experiment. Flow-density optimization is -concluded for now at `e426920`; reopen it only for a demonstrated regression, a -new representative workload, or hardware-counter evidence of a substantial -bottleneck not covered above. +The 2026-07-13 Zen 3 counter run had closed the round at `e426920`. The +2026-09-02 round reopened it on new evidence (a per-step budget and FP +register-file stall counters, see "Packed-channel FMA tracer") and landed the +packed tracer. What remains is bounded by the per-step dependent chain against +the out-of-order window: further gains would need fewer FP-domain micro-ops per +step (e.g. cheaper lane extraction for the integer address path) or a shorter +chain, both of which are now measurable with `perf_kernel.py --steps`. +Reopen on a demonstrated regression, a new representative workload, or +counter evidence of a bottleneck not covered above. Cross-architecture validation remains useful but is not an active optimization target: confirm the FMA dispatch on macOS x86-64 and Windows x86-64, and confirm @@ -478,3 +679,29 @@ python development/flow/check_flow_density.py --dim both --repeats 3 `--restrict-to-mask` / `--no-restrict-to-mask`, and `--threads` to override defaults. The PASS gate is `rel_diff = mean(|ours-ref|)/mean(ref) ≤ 0.15` (`--rel-tol` to override). + +Kernel-level tooling added in the 2026-09-02 round (all in `development/flow/`, +none required by the test suite): + +```bash +# Kernel-only timing + hardware counters (perf stat attached around the timed +# calls only; five-event groups). --steps converts totals into per-step budgets +# using the executed-step count printed by a BIOIMAGE_PROFILE=ON build. +taskset -c 36 python development/flow/perf_kernel.py --case fixture3d --cpu 36 \ + --perf-group core,stall,fp --steps 33959237 + +# Paired A/B of two prebuilt _core .so files (fresh pinned subprocess per run, +# ABBA, noise estimate, density parity). Calibrate with --a X --b X first. +python development/flow/paired_bench.py --a base.so --b cand.so --cases all \ + --threads 1 --cpu 36 --repeats 4 --inner 3 + +# Randomized differential check (400 cases) between two builds or between the +# FMA path and BIOIMAGE_CPP_FLOW_FORCE_SCALAR=1. +python development/flow/differential_check.py --a base.so --b cand.so +python development/flow/differential_check.py --a cand.so --b cand.so \ + --env-b BIOIMAGE_CPP_FLOW_FORCE_SCALAR=1 +``` + +A `BIOIMAGE_PROFILE=ON` build prints, per call, a `[bioimage flow trace]` block +(executed steps, exit reasons, steps-per-particle histogram, per-thread time and +imbalance) in addition to the phase profile. diff --git a/development/flow/_flow_cases.py b/development/flow/_flow_cases.py new file mode 100644 index 0000000..6e5b11a --- /dev/null +++ b/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 diff --git a/development/flow/differential_check.py b/development/flow/differential_check.py new file mode 100644 index 0000000..8f1a32d --- /dev/null +++ b/development/flow/differential_check.py @@ -0,0 +1,198 @@ +"""Randomized differential check of two prebuilt ``_core`` modules. + +Both sides regenerate the same deterministic list of small cases (2D/3D, +length-1 axes, tiny and medium grids, several flow and mask families, Euler +and RK2, tol/dt/n_iter variants, mask restriction on/off, 1 and 4 threads), +compute densities and save them; the parent compares them case by case. + +Gates (accuracy bar, not bitwise identity): the density sum (particle count) +must be equal in every case, and the fraction of differing voxels must stay +below ``--max-diff-frac`` (default 1 %; the generator includes chaotic +scale-50 random flows where FMA-vs-non-FMA contraction moved 0.6 % of the +voxels by up to two counts). The report also states how many cases were +bitwise identical and the worst max |delta|. + +Examples:: + + python development/flow/differential_check.py --a base.so --b cand.so --n-cases 400 + python development/flow/differential_check.py --a base.so --b base.so \ + --env-b BIOIMAGE_CPP_FLOW_FORCE_SCALAR=1 # FMA vs scalar path +""" + +from __future__ import annotations + +import argparse +import json +import os +import subprocess +import sys +import tempfile + +import numpy as np + +from _flow_cases import pin_cpus, preload_core + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument("--a") + parser.add_argument("--b") + parser.add_argument("--env-a", action="append", default=[]) + parser.add_argument("--env-b", action="append", default=[]) + parser.add_argument("--n-cases", type=int, default=400) + parser.add_argument("--seed", type=int, default=0) + parser.add_argument("--max-diff-frac", type=float, default=1e-2, + help="max fraction of differing voxels per case (default 1%%)") + parser.add_argument("--cpu", default=None) + parser.add_argument("--json", default=None) + parser.add_argument("--worker", action="store_true") + parser.add_argument("--so", default=None) + parser.add_argument("--out", default=None) + return parser.parse_args() + + +SHAPES_2D = [(1, 7), (5, 1), (1, 1), (2, 2), (3, 4), (17, 9), (64, 64), (40, 96)] +SHAPES_3D = [(1, 4, 6), (3, 1, 5), (2, 3, 1), (1, 1, 3), (2, 2, 2), (3, 4, 5), (12, 40, 40), (8, 32, 32)] +FLOW_KINDS = ("zero", "int", "half", "normal0.05", "normal0.5", "normal5", "normal50", "swirl") +MASK_KINDS = ("ones", "p0.3", "sparse", "zeros") + + +def make_case(index: int, seed: int) -> dict: + rng = np.random.default_rng(seed * 100003 + index) + ndim = int(rng.integers(2, 4)) + shape = tuple(SHAPES_2D[int(rng.integers(len(SHAPES_2D)))] if ndim == 2 + else SHAPES_3D[int(rng.integers(len(SHAPES_3D)))]) + flow_kind = FLOW_KINDS[int(rng.integers(len(FLOW_KINDS)))] + mask_kind = MASK_KINDS[int(rng.integers(len(MASK_KINDS)))] + return { + "index": index, + "ndim": ndim, + "shape": shape, + "flow": flow_kind, + "mask": mask_kind, + "method": ("euler", "rk2")[int(rng.integers(2))], + "tol": (0.0, 0.005, 0.3)[int(rng.integers(3))], + "dt": (0.0, 0.05, 0.2, 1.0)[int(rng.integers(4))], + "n_iter": (0, 1, 7, 50)[int(rng.integers(4))], + "restrict_to_mask": bool(rng.integers(2)), + "threads": (1, 4)[int(rng.integers(2))], + } + + +def build_arrays(case: dict, seed: int) -> tuple[np.ndarray, np.ndarray]: + rng = np.random.default_rng(seed * 7919 + case["index"] + 1) + ndim, shape = case["ndim"], tuple(case["shape"]) + full = (ndim,) + shape + kind = case["flow"] + if kind == "zero": + flow = np.zeros(full, np.float32) + elif kind == "int": + flow = rng.integers(-3, 4, size=full).astype(np.float32) + elif kind == "half": + flow = rng.integers(-3, 4, size=full).astype(np.float32) + 0.5 + elif kind.startswith("normal"): + flow = rng.normal(scale=float(kind[len("normal"):]), size=full).astype(np.float32) + else: # swirl toward the center with rotation + grids = np.indices(shape, dtype=np.float32) + centers = [(s - 1) / 2.0 for s in shape] + d = [g - c for g, c in zip(grids, centers)] + r = np.sqrt(sum(x * x for x in d)) + 1e-3 + flow = np.stack([-x / r for x in d]).astype(np.float32) + flow[-1] += 0.3 * d[-2] / r + flow[-2] -= 0.3 * d[-1] / r + mk = case["mask"] + if mk == "ones": + mask = np.ones(shape, bool) + elif mk == "p0.3": + mask = rng.random(shape) > 0.3 + elif mk == "sparse": + mask = rng.random(shape) > 0.9 + else: + mask = np.zeros(shape, bool) + return np.ascontiguousarray(flow), mask + + +def worker(args) -> int: + pin_cpus(args.cpu) + preload_core(args.so) + import bioimage_cpp as bic + + densities = {} + meta = [] + for index in range(args.n_cases): + case = make_case(index, args.seed) + flow, mask = build_arrays(case, args.seed) + density = bic.flow.compute_flow_density( + flow, mask, n_iter=case["n_iter"], dt=case["dt"], tol=case["tol"], method=case["method"], + restrict_to_mask=case["restrict_to_mask"], sigma=None, number_of_threads=case["threads"], + ) + densities[f"c{index}"] = density + meta.append({**case, "sum": float(density.sum())}) + np.savez(args.out, **densities) + with open(args.out + ".json", "w") as handle: + json.dump({"so": sys.modules["bioimage_cpp._core"].__file__, "cases": meta, + "backend": getattr(sys.modules["bioimage_cpp._core"], "_flow_trace_backend", lambda: "unknown")()}, handle) + return 0 + + +def _spawn(args, so: str, env_extra: list[str], out: str) -> None: + cmd = [sys.executable, os.path.abspath(__file__), "--worker", "--so", so, "--n-cases", str(args.n_cases), + "--seed", str(args.seed), "--out", out] + if args.cpu: + cmd += ["--cpu", args.cpu] + env = dict(os.environ) + for item in env_extra: + key, _, value = item.partition("=") + env[key] = value + res = subprocess.run(cmd, capture_output=True, text=True, env=env, check=False) + if res.returncode != 0: + raise SystemExit(f"worker failed for {so}:\n{res.stderr}") + + +def main() -> int: + args = parse_args() + if args.worker: + return worker(args) + if not args.a or not args.b: + raise SystemExit("--a and --b are required") + with tempfile.TemporaryDirectory(prefix="flowdiff_", dir=os.environ.get("SCRATCH")) as tmp: + out_a, out_b = os.path.join(tmp, "a.npz"), os.path.join(tmp, "b.npz") + _spawn(args, args.a, args.env_a, out_a) + _spawn(args, args.b, args.env_b, out_b) + da, db = np.load(out_a), np.load(out_b) + with open(out_a + ".json") as fa, open(out_b + ".json") as fb: + ma, mb = json.load(fa), json.load(fb) + print(f"A = {ma['so']} backend={ma['backend']} {args.env_a or ''}") + print(f"B = {mb['so']} backend={mb['backend']} {args.env_b or ''}") + n_identical = 0 + failures = [] + worst_frac, worst_delta = 0.0, 0.0 + rows = [] + for ca, cb in zip(ma["cases"], mb["cases"], strict=True): + a, b = da[f"c{ca['index']}"], db[f"c{cb['index']}"] + n_diff = int(np.count_nonzero(a != b)) + frac = n_diff / max(a.size, 1) + delta = float(np.abs(a - b).max()) if a.size else 0.0 + identical = n_diff == 0 + n_identical += identical + worst_frac, worst_delta = max(worst_frac, frac), max(worst_delta, delta) + sums_equal = ca["sum"] == cb["sum"] + ok = sums_equal and frac <= args.max_diff_frac + rows.append({**ca, "n_diff": n_diff, "diff_frac": frac, "max_delta": delta, "sums_equal": sums_equal, "ok": ok}) + if not ok: + failures.append(rows[-1]) + n = len(rows) + print(f"cases={n} bitwise_identical={n_identical} ({100.0 * n_identical / n:.1f}%) " + f"worst_diff_frac={worst_frac:.5f} worst_max_delta={worst_delta:g} failures={len(failures)}") + for row in failures[:10]: + print(" FAIL", json.dumps({k: row[k] for k in ("index", "ndim", "shape", "flow", "mask", "method", "tol", "dt", + "n_iter", "restrict_to_mask", "threads", "n_diff", + "diff_frac", "max_delta", "sums_equal")})) + if args.json: + with open(args.json, "w") as handle: + json.dump({"a": ma, "b": mb, "rows": rows, "n_identical": n_identical}, handle, indent=1) + return 1 if failures else 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/development/flow/paired_bench.py b/development/flow/paired_bench.py new file mode 100644 index 0000000..9e21528 --- /dev/null +++ b/development/flow/paired_bench.py @@ -0,0 +1,189 @@ +"""Paired A/B benchmark of two prebuilt ``_core`` extension modules. + +Every measurement runs in a fresh, CPU-pinned subprocess that loads ``_core`` +from the requested ``.so`` (see ``_flow_cases.preload_core``), builds the case, +makes one warm call and ``--inner`` timed calls, and reports its min/median. +Per repeat the order is A B B A, so slow drift cancels. The verdict compares +the best-of-all-subprocess minimum of A and B against the per-side noise +(spread of the per-subprocess minima). + +Calibrate first with ``--a X.so --b X.so``: every verdict must read ``noise``. + +Example:: + + python development/flow/paired_bench.py --a base.so --b cand.so \ + --cases fixture3d,fixture2d --threads 1 --cpu 36 --repeats 4 --inner 3 +""" + +from __future__ import annotations + +import argparse +import hashlib +import json +import os +import subprocess +import sys +import time +from statistics import median + +import numpy as np + +from _flow_cases import ( + ADVERSARIAL_CASES, + CASES, + FIXTURE_CASES, + add_param_args, + build_case, + make_runner, + params_from_args, + pin_cpus, + preload_core, +) + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument("--a", help="baseline _core .so") + parser.add_argument("--b", help="candidate _core .so") + parser.add_argument("--env-a", action="append", default=[], help="VAR=VALUE for side A") + parser.add_argument("--env-b", action="append", default=[], help="VAR=VALUE for side B") + parser.add_argument("--cases", default="fixture3d,fixture2d", + help="comma list, or 'fixtures', 'adversarial', 'all'") + parser.add_argument("--threads", type=int, default=1) + parser.add_argument("--cpu", default=None, help="CPUs for the workers, e.g. 36 or 36,37") + parser.add_argument("--repeats", type=int, default=4, help="ABBA rounds per case") + parser.add_argument("--inner", type=int, default=3, help="timed calls per subprocess") + parser.add_argument("--wrapper", action="store_true", help="time the Python wrapper, not _core") + parser.add_argument("--json", default=None) + parser.add_argument("--noise-margin", type=float, default=0.5, help="extra %% added to 2*noise") + add_param_args(parser) + # worker mode + parser.add_argument("--worker", action="store_true") + parser.add_argument("--so", default=None) + parser.add_argument("--case", default=None) + parser.add_argument("--params", default="{}") + return parser.parse_args() + + +def worker(args) -> int: + pin_cpus(args.cpu) + preload_core(args.so) + import bioimage_cpp as bic + + params = json.loads(args.params) + flow, mask = build_case(args.case) + ndim = flow.shape[0] + bare = not args.wrapper + run = make_runner(bare, ndim) + mask_arg = np.ascontiguousarray(mask, dtype=np.uint8) if bare else mask + import gc + + gc.collect() + gc.disable() + run(flow, mask_arg, args.threads, **params) # warm + times = [] + result = None + for _ in range(args.inner): + start = time.perf_counter() + result = run(flow, mask_arg, args.threads, **params) + times.append(time.perf_counter() - start) + assert result is not None + print(json.dumps({ + "case": args.case, + "so": sys.modules["bioimage_cpp._core"].__file__, + "times": times, + "min": min(times), + "median": median(times), + "particles": float(result.sum()), + "sha256": hashlib.sha256(result.tobytes()).hexdigest()[:16], + "backend": getattr(sys.modules["bioimage_cpp._core"], "_flow_trace_backend", lambda: "unknown")(), + })) + return 0 + + +def _spawn(args, so: str, env_extra: list[str], case: str, params: dict) -> dict: + cmd = [sys.executable, os.path.abspath(__file__), "--worker", "--so", so, "--case", case, + "--threads", str(args.threads), "--inner", str(args.inner), "--params", json.dumps(params)] + if args.cpu: + cmd = ["taskset", "-c", args.cpu] + cmd + ["--cpu", args.cpu] + if args.wrapper: + cmd.append("--wrapper") + env = dict(os.environ) + for item in env_extra: + key, _, value = item.partition("=") + env[key] = value + out = subprocess.run(cmd, capture_output=True, text=True, env=env, check=False) + if out.returncode != 0: + raise SystemExit(f"worker failed for {case} ({so}):\n{out.stderr}") + line = [ln for ln in out.stdout.splitlines() if ln.startswith("{")][-1] + return json.loads(line) + + +def _resolve_cases(spec: str) -> list[str]: + if spec == "fixtures": + return list(FIXTURE_CASES) + if spec == "adversarial": + return list(ADVERSARIAL_CASES) + if spec == "all": + return list(FIXTURE_CASES) + list(ADVERSARIAL_CASES) + names = [c for c in spec.split(",") if c] + for name in names: + if name not in CASES: + raise SystemExit(f"unknown case {name!r}") + return names + + +def main() -> int: + args = parse_args() + if args.worker: + return worker(args) + if not args.a or not args.b: + raise SystemExit("--a and --b are required") + params = params_from_args(args) + cases = _resolve_cases(args.cases) + sides = {"A": (args.a, args.env_a), "B": (args.b, args.env_b)} + print(f"A = {args.a} {args.env_a or ''}\nB = {args.b} {args.env_b or ''}") + print(f"threads={args.threads} cpu={args.cpu} repeats={args.repeats} (ABBA) inner={args.inner} " + f"{'wrapper' if args.wrapper else 'bare'} params={params}") + header = (f"{'case':12s} {'thr':>3s} {'A_min':>9s} {'B_min':>9s} {'dMin%':>7s} " + f"{'A_med':>9s} {'B_med':>9s} {'dMed%':>7s} {'nzA%':>5s} {'nzB%':>5s} par verdict") + print(header) + print("-" * len(header)) + report = [] + for case in cases: + runs: dict[str, list[dict]] = {"A": [], "B": []} + for _ in range(args.repeats): + for side in ("A", "B", "B", "A"): + so, env_extra = sides[side] + runs[side].append(_spawn(args, so, env_extra, case, params)) + mins = {s: [r["min"] for r in runs[s]] for s in runs} + meds = {s: [r["median"] for r in runs[s]] for s in runs} + a_min, b_min = min(mins["A"]), min(mins["B"]) + a_med, b_med = median(meds["A"]), median(meds["B"]) + noise = {s: 100.0 * (max(mins[s]) - min(mins[s])) / min(mins[s]) for s in mins} + d_min = 100.0 * (b_min - a_min) / a_min + d_med = 100.0 * (b_med - a_med) / a_med + shas = {s: {r["sha256"] for r in runs[s]} for s in runs} + parity = "=" if shas["A"] == shas["B"] and len(shas["A"]) == 1 else ("~" if len(shas["A"]) == 1 and len(shas["B"]) == 1 else "!") + threshold = 2.0 * max(noise.values()) + args.noise_margin + verdict = "significant" if abs(d_min) > threshold else "noise" + row = {"case": case, "threads": args.threads, "a_min": a_min, "b_min": b_min, "d_min_pct": d_min, + "a_med": a_med, "b_med": b_med, "d_med_pct": d_med, "noise_a_pct": noise["A"], + "noise_b_pct": noise["B"], "parity": parity, "verdict": verdict, + "particles_a": runs["A"][0]["particles"], "particles_b": runs["B"][0]["particles"], + "backend_a": runs["A"][0]["backend"], "backend_b": runs["B"][0]["backend"]} + report.append(row) + print(f"{case:12s} {args.threads:3d} {a_min:9.4f} {b_min:9.4f} {d_min:+7.2f} " + f"{a_med:9.4f} {b_med:9.4f} {d_med:+7.2f} {noise['A']:5.2f} {noise['B']:5.2f} {parity} {verdict}", + flush=True) + print("parity: '=' identical densities, '~' each side deterministic but different, '!' nondeterministic") + if args.json: + with open(args.json, "w") as handle: + json.dump({"a": args.a, "b": args.b, "threads": args.threads, "cpu": args.cpu, + "repeats": args.repeats, "inner": args.inner, "params": params, "rows": report}, + handle, indent=2) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/development/flow/perf_kernel.py b/development/flow/perf_kernel.py new file mode 100644 index 0000000..8274707 --- /dev/null +++ b/development/flow/perf_kernel.py @@ -0,0 +1,311 @@ +"""Kernel-only timing and hardware-counter harness for compute_flow_density. + +Loads one case once, pre-converts the inputs, warms the kernel, then times N +calls. With ``--perf-group`` it attaches ``perf stat -p `` around only the +timed calls (one attach per event group, five events per group because +``nmi_watchdog=1`` leaves five programmable counters on Zen 3), so process +start-up, HDF5 loading and NumPy conversions are excluded from the counts. + +Examples:: + + taskset -c 36 python development/flow/perf_kernel.py --case fixture3d --perf-group all + python development/flow/perf_kernel.py --case fixture3d --overhead + python development/flow/perf_kernel.py --so /path/to/_core.so --case fixture2d --steps 4200000 + +``--steps S`` (the executed particle-step count printed by a +``BIOIMAGE_PROFILE=ON`` build) converts totals into per-step budgets. +""" + +from __future__ import annotations + +import argparse +import gc +import hashlib +import json +import os +import signal +import subprocess +import sys +import time +from statistics import median + +import numpy as np + +from _flow_cases import ( + CASES, + DEFAULTS, + add_param_args, + build_case, + make_runner, + params_from_args, + pin_cpus, + preload_core, +) + +PERF_GROUPS: dict[str, list[str]] = { + "core": ["cycles", "instructions", "branches", "branch-misses", "ls_dc_accesses"], + "l1l2": [ + "cycles", "L1-dcache-load-misses", "l2_cache_misses_from_dc_misses", + "ls_dmnd_fills_from_sys.lcl_l2", "ls_dmnd_fills_from_sys.int_cache", + ], + "far": [ + "cycles", "ls_dmnd_fills_from_sys.ext_cache_local", + "ls_dmnd_fills_from_sys.mem_io_local", "ls_dmnd_fills_from_sys.mem_io_remote", + "ls_pref_instr_disp", + ], + "fp": [ + "cycles", "fp_ret_sse_avx_ops.all", "fp_ret_sse_avx_ops.mac_flops", + "fp_ret_sse_avx_ops.add_sub_flops", "fp_ret_sse_avx_ops.mult_flops", + ], + "tlb": [ + "cycles", "fp_disp_faults.xmm_fill_fault", "ls_l1_d_tlb_miss.all", + "ls_l1_d_tlb_miss.tlb_reload_2m_l2_hit", "ls_l1_d_tlb_miss.tlb_reload_4k_l2_hit", + ], + "stall": [ + "cycles", "stalled-cycles-backend", "stalled-cycles-frontend", + "de_dis_dispatch_token_stalls1.fp_reg_file_rsrc_stall", + "de_dis_dispatch_token_stalls1.load_queue_rsrc_stall", + ], + "misal": ["cycles", "ls_misal_loads.ma64", "ls_misal_loads.ma4k", "ls_dc_accesses", "instructions"], +} + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument("--case", default="fixture3d", choices=sorted(CASES)) + parser.add_argument("--threads", type=int, default=1) + parser.add_argument("--warm", type=int, default=1) + parser.add_argument("--n-timed", type=int, default=8) + parser.add_argument("--bare", action="store_true", default=True, help="call _core directly (default)") + parser.add_argument("--wrapper", action="store_true", help="time the public Python wrapper instead") + parser.add_argument("--so", default=None, help="load bioimage_cpp._core from this path") + parser.add_argument("--cpu", default=None, help="pin to these CPUs, e.g. 36 or 36,37") + parser.add_argument("--perf-group", default=None, help="comma list of groups or 'all'") + parser.add_argument("--perf-events", default=None, help="explicit comma list of events (one group)") + parser.add_argument("--perf-out", default=None, help="directory for perf CSV files") + parser.add_argument("--perf-cpu", default=None, help="CPU to pin the perf process to") + parser.add_argument("--steps", type=float, default=None, help="executed particle-steps per call") + parser.add_argument("--json", default=None, help="write results to this JSON file") + parser.add_argument("--overhead", action="store_true", help="measure wrapper/binding overheads") + parser.add_argument("--attach-wait", type=float, default=0.0, help="print pid and sleep before timing") + parser.add_argument("--timeout", type=float, default=60.0) + add_param_args(parser) + return parser.parse_args() + + +def sha256(array: np.ndarray) -> str: + return hashlib.sha256(np.ascontiguousarray(array).tobytes()).hexdigest()[:16] + + +def _time_calls(run, flow, mask, threads, params, n: int) -> tuple[list[float], np.ndarray]: + times = [] + result = None + for _ in range(n): + start = time.perf_counter() + result = run(flow, mask, threads, **params) + times.append(time.perf_counter() - start) + return times, result + + +def _parse_perf_csv(path: str) -> dict[str, float | None]: + counts: dict[str, float | None] = {} + with open(path) as handle: + for line in handle: + line = line.strip() + if not line or line.startswith("#"): + continue + fields = line.split(",") + value, event = fields[0], fields[2] if len(fields) > 2 else fields[-1] + event = event.removesuffix(":u") + try: + counts[event] = float(value) + except ValueError: + counts[event] = None # / + return counts + + +# The fp_ret_sse_avx_ops.* umasks cannot be scheduled as one hardware group on +# this PMU (every member reports ), so that set is attached +# ungrouped and multiplexed instead. +UNGROUPED = {"fp"} + + +def _perf_attach(events: list[str], out_csv: str, perf_cpu: str | None, grouped: bool = True) -> subprocess.Popen: + if grouped: + group = "{" + ",".join(f"{e}:u" for e in events) + "}" + else: + group = ",".join(f"{e}:u" for e in events) + cmd = ["perf", "stat", "-x", ",", "-o", out_csv, "-p", str(os.getpid()), "-e", group] + proc = subprocess.Popen(cmd, stdout=subprocess.DEVNULL, stderr=subprocess.PIPE) + if perf_cpu: + try: + os.sched_setaffinity(proc.pid, {int(c) for c in perf_cpu.split(",")}) + except OSError: + pass + time.sleep(1.0) # let perf open the counters before the timed calls + return proc + + +def _perf_detach(proc: subprocess.Popen) -> str: + proc.send_signal(signal.SIGINT) + _, err = proc.communicate(timeout=30) + return err.decode(errors="replace") + + +def _derived(counts: dict[str, float | None], n_calls: int, steps: float | None) -> dict[str, float]: + out: dict[str, float] = {} + c = counts.get("cycles") + i = counts.get("instructions") + if c and i: + out["ipc"] = i / c + loads = counts.get("ls_dc_accesses") + if loads and i: + out["loads_per_instr"] = loads / i + l1m = counts.get("L1-dcache-load-misses") + if loads and l1m is not None: + out["l1d_miss_pct"] = 100.0 * l1m / loads + l2m = counts.get("l2_cache_misses_from_dc_misses") + if l1m and l2m is not None: + out["l2_miss_pct_of_l1_misses"] = 100.0 * l2m / l1m + br = counts.get("branches") + brm = counts.get("branch-misses") + if br and brm is not None: + out["branch_miss_pct"] = 100.0 * brm / br + if steps: + denom = n_calls * steps + for key, value in counts.items(): + if value is not None: + out[f"{key}/step"] = value / denom + return out + + +def run_overhead(args, flow, mask, mask_u8, ndim, params) -> dict: + import bioimage_cpp as bic + + def best(fn, n=10): + ts = [] + for _ in range(n): + t = time.perf_counter() + fn() + ts.append(time.perf_counter() - t) + return min(ts) + + bare = make_runner(True, ndim) + wrapper = make_runner(False, ndim) + threads = args.threads + n_iter0 = {**params, "n_iter": 0} + results = { + "wrapper_s": best(lambda: wrapper(flow, mask, threads, **params), 3), + "bare_s": best(lambda: bare(flow, mask_u8, threads, **params), 3), + "bare_n_iter0_s": best(lambda: bare(flow, mask_u8, threads, **n_iter0)), + "np_isfinite_all_s": best(lambda: np.isfinite(flow).all()), + "mask_astype_u8_s": best(lambda: mask.astype(np.uint8)), + "flow_ascontiguous_noop_s": best(lambda: np.ascontiguousarray(flow, dtype=np.float32)), + "flow_bytes": int(flow.nbytes), + } + _ = bic + return results + + +def main() -> int: + args = parse_args() + cpus = pin_cpus(args.cpu) + preload_core(args.so) + import bioimage_cpp as bic + + flow, mask = build_case(args.case) + ndim = flow.shape[0] + mask_u8 = np.ascontiguousarray(mask, dtype=np.uint8) + params = params_from_args(args) + bare = not args.wrapper + run = make_runner(bare, ndim) + mask_arg = mask_u8 if bare else mask + + info = { + "case": args.case, + "shape": list(flow.shape[1:]), + "particles_in_mask": int(mask.sum()), + "threads": args.threads, + "cpus": cpus, + "bare": bare, + "params": {**DEFAULTS, **params}, + "so": sys.modules["bioimage_cpp._core"].__file__, + "backend": getattr(sys.modules["bioimage_cpp._core"], "_flow_trace_backend", lambda: "unknown")(), + "pid": os.getpid(), + } + print(json.dumps(info), flush=True) + + if args.overhead: + results = run_overhead(args, flow, mask, mask_u8, ndim, params) + info["overhead"] = results + for key, value in results.items(): + print(f" {key:28s} {value:.6f}" if isinstance(value, float) else f" {key:28s} {value}") + if args.json: + with open(args.json, "w") as handle: + json.dump(info, handle, indent=2) + return 0 + + gc.collect() + gc.disable() + if args.warm: + _time_calls(run, flow, mask_arg, args.threads, params, args.warm) + if args.attach_wait > 0: + print(f"pid {os.getpid()} sleeping {args.attach_wait}s for manual attach", flush=True) + time.sleep(args.attach_wait) + + groups: list[tuple[str, list[str]]] = [] + if args.perf_events: + groups.append(("custom", args.perf_events.split(","))) + if args.perf_group: + names = list(PERF_GROUPS) if args.perf_group == "all" else args.perf_group.split(",") + if args.perf_group == "all": + names = [n for n in names if n != "misal"] + groups.extend((n, PERF_GROUPS[n]) for n in names) + + all_times: list[float] = [] + result = None + counters: dict[str, dict] = {} + if not groups: + all_times, result = _time_calls(run, flow, mask_arg, args.threads, params, args.n_timed) + else: + out_dir = args.perf_out or "." + os.makedirs(out_dir, exist_ok=True) + for name, events in groups: + csv_path = os.path.join(out_dir, f"{args.case}_{args.threads}T_{name}.csv") + proc = _perf_attach(events, csv_path, args.perf_cpu, grouped=name not in UNGROUPED) + times, result = _time_calls(run, flow, mask_arg, args.threads, params, args.n_timed) + err = _perf_detach(proc) + all_times.extend(times) + counts = _parse_perf_csv(csv_path) + counters[name] = {"counts": counts, "derived": _derived(counts, args.n_timed, args.steps), + "wall_min": min(times), "wall_median": median(times)} + if err.strip() and "Warning" not in err: + counters[name]["perf_stderr"] = err.strip()[:500] + gc.enable() + + assert result is not None + info.update({ + "n_timed": args.n_timed, + "wall_min": min(all_times), + "wall_median": median(all_times), + "particles_out": float(result.sum()), + "density_sha256": sha256(result), + "counters": counters, + }) + print(f"{args.case} threads={args.threads}: min={min(all_times):.4f}s median={median(all_times):.4f}s " + f"particles={info['particles_out']:.0f} sha={info['density_sha256']}") + for name, block in counters.items(): + print(f" [{name}] wall_min={block['wall_min']:.4f}s") + for event, value in block["counts"].items(): + per_call = "" if value is None else f" ({value / args.n_timed:,.0f}/call)" + print(f" {event:52s} {'' if value is None else f'{value:,.0f}'}{per_call}") + for key, value in block["derived"].items(): + print(f" -> {key:48s} {value:,.4f}") + if args.json: + with open(args.json, "w") as handle: + json.dump(info, handle, indent=2) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/include/bioimage_cpp/detail/finite.hxx b/include/bioimage_cpp/detail/finite.hxx new file mode 100644 index 0000000..820e9b9 --- /dev/null +++ b/include/bioimage_cpp/detail/finite.hxx @@ -0,0 +1,59 @@ +#pragma once + +#include "bioimage_cpp/detail/threading.hxx" + +#include +#include +#include +#include +#include + +namespace bioimage_cpp::detail { + +// True iff every float in [begin, end) is finite. The inner loop is +// branch-free (an OR-reduction over the exponent-all-ones test on the bit +// pattern) so the compiler vectorizes it on the SSE2 baseline; the block +// granularity gives an early exit once a non-finite value has been seen. +inline bool all_finite_range( + const float *data, const std::size_t begin, const std::size_t end +) noexcept { + constexpr std::size_t block = 4096; + for (std::size_t start = begin; start < end; start += block) { + const std::size_t stop = std::min(end, start + block); + std::uint32_t any_special = 0; + for (std::size_t i = start; i < stop; ++i) { + std::uint32_t bits; + std::memcpy(&bits, data + i, sizeof(bits)); + any_special |= static_cast((bits & 0x7f800000u) == 0x7f800000u); + } + if (any_special != 0) { + return false; + } + } + return true; +} + +// Parallel finiteness check over `n` floats. Worker threads are only used +// when every worker gets at least ~1M elements, so small arrays are checked +// inline without spawning threads. +inline bool all_finite( + const float *data, const std::size_t n, const std::size_t number_of_threads +) { + constexpr std::size_t min_elements_per_thread = std::size_t(1) << 20; + const std::size_t max_useful = std::max(1, n / min_elements_per_thread); + const std::size_t n_threads = + normalize_thread_count(std::min(number_of_threads, max_useful), n); + if (n_threads <= 1) { + return all_finite_range(data, 0, n); + } + std::vector chunk_ok(n_threads, 1); + parallel_for_chunks( + n_threads, n, + [&](const std::size_t thread_id, const std::size_t begin, const std::size_t end) { + chunk_ok[thread_id] = all_finite_range(data, begin, end) ? 1 : 0; + } + ); + return std::all_of(chunk_ok.begin(), chunk_ok.end(), [](std::uint8_t ok) { return ok != 0; }); +} + +} // namespace bioimage_cpp::detail diff --git a/include/bioimage_cpp/detail/force_inline.hxx b/include/bioimage_cpp/detail/force_inline.hxx new file mode 100644 index 0000000..e080dcc --- /dev/null +++ b/include/bioimage_cpp/detail/force_inline.hxx @@ -0,0 +1,16 @@ +#pragma once + +// BIOIMAGE_FORCE_INLINE: request inlining regardless of the compiler's +// per-translation-unit growth budget. Use it only for small helpers that sit +// inside a hot loop and whose being called out-of-line would change codegen +// materially (e.g. a per-sample interpolation kernel). It also removes a +// link-time hazard for ISA-specialized translation units: a helper that is not +// inlined becomes a weak COMDAT symbol, and the linker may pick the copy +// compiled with a different instruction set. +#if defined(_MSC_VER) && !defined(__clang__) +#define BIOIMAGE_FORCE_INLINE __forceinline +#elif defined(__GNUC__) || defined(__clang__) +#define BIOIMAGE_FORCE_INLINE inline __attribute__((always_inline)) +#else +#define BIOIMAGE_FORCE_INLINE inline +#endif diff --git a/include/bioimage_cpp/flow/flow_density.hxx b/include/bioimage_cpp/flow/flow_density.hxx index 3820aed..f743dea 100644 --- a/include/bioimage_cpp/flow/flow_density.hxx +++ b/include/bioimage_cpp/flow/flow_density.hxx @@ -1,6 +1,7 @@ #pragma once #include "bioimage_cpp/array_view.hxx" +#include "bioimage_cpp/detail/force_inline.hxx" #include "bioimage_cpp/detail/profile.hxx" #include "bioimage_cpp/detail/threading.hxx" @@ -8,8 +9,15 @@ #include #include #include +#include +#include #include +#ifdef BIOIMAGE_PROFILE +#include +#include +#endif + #if defined(BIOIMAGE_FLOW_FMA_DISPATCH) && defined(_MSC_VER) #include #include @@ -39,6 +47,125 @@ GridLayout make_grid_layout( return layout; } +// Executed-step accounting for profile builds. `record(steps, exit)` is called +// once per traced particle with the number of integration steps whose samples +// were computed and why tracing stopped. Outside BIOIMAGE_PROFILE builds the +// Null variant is used and every call site is guarded by `if constexpr`, so +// the production kernels carry no counters. +struct NullTraceStats { + static constexpr bool enabled = false; + enum Exit : int { Converged = 0, LeftMask = 1, MaxIter = 2 }; + explicit NullTraceStats(std::size_t = 0) noexcept {} + void record(std::size_t, int) noexcept {} + void merge(const NullTraceStats &) noexcept {} + static void report(const std::vector &, std::size_t, bool) noexcept {} +}; + +#ifdef BIOIMAGE_PROFILE +struct TraceStats { + static constexpr bool enabled = true; + enum Exit : int { Converged = 0, LeftMask = 1, MaxIter = 2 }; + + std::uint64_t steps = 0; + std::uint64_t particles = 0; + std::array exits{}; + std::vector histogram; // histogram[s] = particles that executed s steps + double seconds = 0.0; + + explicit TraceStats(std::size_t n_iter = 0) : histogram(n_iter + 1, 0) {} + + void record(std::size_t n_steps, int reason) noexcept { + steps += n_steps; + ++particles; + ++exits[static_cast(reason)]; + if (n_steps < histogram.size()) { + ++histogram[n_steps]; + } + } + + void merge(const TraceStats &other) { + steps += other.steps; + particles += other.particles; + for (std::size_t i = 0; i < exits.size(); ++i) { + exits[i] += other.exits[i]; + } + if (histogram.size() < other.histogram.size()) { + histogram.resize(other.histogram.size(), 0); + } + for (std::size_t i = 0; i < other.histogram.size(); ++i) { + histogram[i] += other.histogram[i]; + } + seconds += other.seconds; + } + + // Print the merged totals plus one line per worker (load balance). + static void report( + const std::vector &per_thread, const std::size_t n_iter, const bool rk2 + ) { + TraceStats total(n_iter); + for (const auto &t : per_thread) { + total.merge(t); + } + const double particles = static_cast(std::max(total.particles, 1)); + std::fprintf(stderr, "[bioimage flow trace]\n"); + std::fprintf( + stderr, " particles %llu n_iter %zu %s\n", + static_cast(total.particles), n_iter, rk2 ? "rk2" : "euler" + ); + std::fprintf( + stderr, " steps %llu mean_steps_per_particle %.3f samples_per_step %d\n", + static_cast(total.steps), total.steps / particles, rk2 ? 2 : 1 + ); + std::fprintf( + stderr, " exits converged %llu (%.1f%%) left_mask %llu (%.1f%%) hit_n_iter %llu (%.1f%%)\n", + static_cast(total.exits[0]), 100.0 * total.exits[0] / particles, + static_cast(total.exits[1]), 100.0 * total.exits[1] / particles, + static_cast(total.exits[2]), 100.0 * total.exits[2] / particles + ); + std::fprintf(stderr, " histogram(steps:particles)"); + std::size_t lo = 0; + for (std::size_t width = 1; lo < total.histogram.size(); width *= 2) { + const std::size_t hi = std::min(total.histogram.size(), lo + width); + std::uint64_t count = 0; + for (std::size_t s = lo; s < hi; ++s) { + count += total.histogram[s]; + } + if (hi - lo == 1) { + std::fprintf(stderr, " %zu:%llu", lo, static_cast(count)); + } else { + std::fprintf(stderr, " %zu-%zu:%llu", lo, hi - 1, static_cast(count)); + } + lo = hi; + } + std::fprintf(stderr, "\n"); + double max_seconds = 0.0, sum_seconds = 0.0; + std::uint64_t max_steps = 0; + for (std::size_t t = 0; t < per_thread.size(); ++t) { + std::fprintf( + stderr, " thread %zu %.4f s steps %llu particles %llu\n", t, per_thread[t].seconds, + static_cast(per_thread[t].steps), + static_cast(per_thread[t].particles) + ); + max_seconds = std::max(max_seconds, per_thread[t].seconds); + sum_seconds += per_thread[t].seconds; + max_steps = std::max(max_steps, per_thread[t].steps); + } + if (per_thread.size() > 1) { + const double mean_seconds = sum_seconds / per_thread.size(); + const double mean_steps = static_cast(total.steps) / per_thread.size(); + std::fprintf( + stderr, " imbalance time max/mean-1 %.1f%% steps max/mean-1 %.1f%%\n", + 100.0 * (max_seconds / std::max(mean_seconds, 1e-12) - 1.0), + 100.0 * (max_steps / std::max(mean_steps, 1e-12) - 1.0) + ); + } + } +}; +using ActiveTraceStats = TraceStats; +#else +using ActiveTraceStats = NullTraceStats; +#endif + // Linearly interpolate all D flow channels at `position` and write the result // to `out`. This is explicit bilinear (D==2) / trilinear (D==3) sampling rather // than a generic 2^D corner table: the lower/upper index and fractional weight @@ -52,7 +179,7 @@ GridLayout make_grid_layout( // software routine while the cast is a single instruction. At an upper boundary // the lower and upper index coincide, matching nearest-boundary behavior. template -inline void sample_flow( +BIOIMAGE_FORCE_INLINE void sample_flow( const std::array &channels, const std::array &position, const GridLayout &grid, @@ -115,7 +242,7 @@ inline void sample_flow( } template -inline std::ptrdiff_t round_to_flat_index( +BIOIMAGE_FORCE_INLINE std::ptrdiff_t round_to_flat_index( const std::array &position, const GridLayout &grid ) { @@ -139,7 +266,7 @@ inline std::ptrdiff_t round_to_flat_index( } template -inline bool position_is_in_mask( +BIOIMAGE_FORCE_INLINE bool position_is_in_mask( const std::array &position, const GridLayout &grid, const std::uint8_t *mask @@ -165,23 +292,25 @@ namespace detail { // are compile-time parameters so the per-step branches on them fold away and // the sampler/RK2/convergence/mask code inlines into one specialized loop. // -// CodegenVariant must match the instantiating trace_all (see there): the -// compiler is free to emit this function out-of-line (observed with GCC 14 -// under mild size pressure), and without the tag the FMA translation unit -// would emit AVX code under the same weak symbol name as the portable -// instantiation, letting COMDAT selection ship AVX code to the portable -// fallback path (SIGILL on pre-AVX CPUs) or silently discard the FMA kernel. +// CodegenVariant must match the instantiating trace_all (see there): without +// the tag the FMA translation unit would emit AVX code under the same weak +// symbol name as the portable instantiation, letting COMDAT selection ship +// AVX code to the portable fallback path (SIGILL on pre-AVX CPUs) or silently +// discard the FMA kernel. The per-step helpers are additionally force-inlined +// because GCC 14 was observed to stop inlining them once the FMA translation +// unit grew (which both slowed the loop ~1.5x and recreated the hazard). template < std::size_t D, bool UseRK2, bool CheckConvergence, bool RestrictToMask, bool CodegenVariant = false> -inline void trace_particle( +BIOIMAGE_FORCE_INLINE void trace_particle( std::array &position, const std::array &channels, const GridLayout &grid, const std::uint8_t *mask, const std::size_t n_iter, const float dt, - const float tol + const float tol, + ActiveTraceStats &stats ) { const auto clip = [&grid](std::array &p) { for (std::size_t axis = 0; axis < D; ++axis) { @@ -192,6 +321,8 @@ inline void trace_particle( } } }; + [[maybe_unused]] std::size_t executed_steps = 0; + [[maybe_unused]] int exit_reason = ActiveTraceStats::MaxIter; // When restricting to the mask, only in-mask (hence in-bounds) endpoints are // ever committed and the seed is an in-bounds integer voxel, so `position` @@ -203,6 +334,9 @@ inline void trace_particle( } for (std::size_t iter = 0; iter < n_iter; ++iter) { + if constexpr (ActiveTraceStats::enabled) { + ++executed_steps; + } std::array step{}; sample_flow(channels, position, grid, step); @@ -224,6 +358,9 @@ inline void trace_particle( } } if (max_step < tol) { + if constexpr (ActiveTraceStats::enabled) { + exit_reason = ActiveTraceStats::Converged; + } break; } } @@ -238,6 +375,9 @@ inline void trace_particle( // at its last in-mask position (only the endpoint is mask-tested, // not the RK2 midpoint). if (!position_is_in_mask(proposed, grid, mask)) { + if constexpr (ActiveTraceStats::enabled) { + exit_reason = ActiveTraceStats::LeftMask; + } break; } } else { @@ -245,6 +385,9 @@ inline void trace_particle( } position = proposed; } + if constexpr (ActiveTraceStats::enabled) { + stats.record(executed_steps, exit_reason); + } } // Trace K consecutive particles in lockstep. The trajectories are independent, @@ -257,16 +400,22 @@ inline void trace_particle( template < std::size_t D, std::size_t K, bool UseRK2, bool CheckConvergence, bool RestrictToMask, bool CodegenVariant = false> -inline void trace_particle_block( +BIOIMAGE_FORCE_INLINE void trace_particle_block( std::array *positions, const std::array &channels, const GridLayout &grid, const std::uint8_t *mask, const std::size_t n_iter, const float dt, - const float tol + const float tol, + ActiveTraceStats &stats ) { static_assert(K >= 2, "use trace_particle for single trajectories"); + [[maybe_unused]] std::array lane_steps{}; + [[maybe_unused]] std::array lane_exit{}; + if constexpr (ActiveTraceStats::enabled) { + lane_exit.fill(ActiveTraceStats::MaxIter); + } const auto clip = [&grid](std::array &p) { for (std::size_t axis = 0; axis < D; ++axis) { @@ -296,6 +445,9 @@ inline void trace_particle_block( if (!alive[k]) { continue; } + if constexpr (ActiveTraceStats::enabled) { + ++lane_steps[k]; + } std::array step{}; sample_flow(channels, pos[k], grid, step); @@ -318,6 +470,9 @@ inline void trace_particle_block( } if (max_step < tol) { alive[k] = false; + if constexpr (ActiveTraceStats::enabled) { + lane_exit[k] = ActiveTraceStats::Converged; + } continue; } } @@ -330,6 +485,9 @@ inline void trace_particle_block( if constexpr (RestrictToMask) { if (!position_is_in_mask(proposed, grid, mask)) { alive[k] = false; + if constexpr (ActiveTraceStats::enabled) { + lane_exit[k] = ActiveTraceStats::LeftMask; + } continue; } } else { @@ -350,6 +508,11 @@ inline void trace_particle_block( for (std::size_t k = 0; k < K; ++k) { positions[k] = pos[k]; } + if constexpr (ActiveTraceStats::enabled) { + for (std::size_t k = 0; k < K; ++k) { + stats.record(lane_steps[k], lane_exit[k]); + } + } } // CodegenVariant gives separately compiled ISA variants a distinct linker @@ -372,10 +535,21 @@ void trace_all( // across all integration steps in one fan-out. Trajectories are independent // until the sequential scatter, so no global alive state, per-step barrier, // or per-step alive scan is needed. + ActiveTraceStats null_stats{}; + std::vector per_thread_stats; + if constexpr (ActiveTraceStats::enabled) { + per_thread_stats.assign(n_threads, ActiveTraceStats(n_iter)); + } ::bioimage_cpp::detail::parallel_for_chunks( n_threads, positions.size(), - [&](const std::size_t, const std::size_t begin, const std::size_t end) { + [&](const std::size_t thread_id, const std::size_t begin, const std::size_t end) { + ActiveTraceStats *stats = &null_stats; + (void)thread_id; +#ifdef BIOIMAGE_PROFILE + stats = &per_thread_stats[thread_id]; + const auto chunk_start = std::chrono::steady_clock::now(); +#endif // Lockstep interleaving only pays for RK2: its two dependent // samples per step leave latency bubbles that other lanes fill. // The shorter Euler chain measured ~10% slower when interleaved @@ -389,18 +563,34 @@ void trace_all( trace_particle_block< D, K, UseRK2, CheckConvergence, RestrictToMask, CodegenVariant>( - &positions[i], channels, grid, mask, n_iter, dt, tol + &positions[i], channels, grid, mask, n_iter, dt, tol, *stats ); } } for (; i < end; ++i) { trace_particle< D, UseRK2, CheckConvergence, RestrictToMask, CodegenVariant>( - positions[i], channels, grid, mask, n_iter, dt, tol + positions[i], channels, grid, mask, n_iter, dt, tol, *stats ); } +#ifdef BIOIMAGE_PROFILE + stats->seconds += std::chrono::duration( + std::chrono::steady_clock::now() - chunk_start + ).count(); +#endif } ); + if constexpr (ActiveTraceStats::enabled) { + ActiveTraceStats::report(per_thread_stats, n_iter, UseRK2); + } +} + +// Runtime opt-out of the FMA-specialized tracer (mirrors +// BIOIMAGE_CPP_FILTERS_FORCE_SCALAR). Read on every call, not cached, so tests +// can toggle it in-process. +inline bool force_scalar_requested() noexcept { + const char *value = std::getenv("BIOIMAGE_CPP_FLOW_FORCE_SCALAR"); + return value != nullptr && value[0] != '\0' && std::strcmp(value, "0") != 0; } #if defined(BIOIMAGE_FLOW_FMA_DISPATCH) @@ -470,7 +660,7 @@ bool try_trace_all_fma( const float dt, const float tol ) { - if (!runtime_fma_supported()) { + if (force_scalar_requested() || !runtime_fma_supported()) { return false; } if constexpr (D == 2) { @@ -489,6 +679,18 @@ bool try_trace_all_fma( } // namespace detail +// Name of the tracer the default RK2/convergence/mask path will use on this +// machine: "fma" (runtime-dispatched FMA translation unit) or "scalar" +// (portable kernel). Exposed to Python as `_core._flow_trace_backend`. +inline const char *trace_backend() noexcept { +#if defined(BIOIMAGE_FLOW_FMA_DISPATCH) + if (!detail::force_scalar_requested() && detail::runtime_fma_supported()) { + return "fma"; + } +#endif + return "scalar"; +} + // Preconditions (validated in the binding layer): // * flow.ndim() == D + 1, flow.shape[0] == D, flow.shape[1..] == fg_mask.shape // * fg_mask.ndim() == D and density.shape == fg_mask.shape diff --git a/src/bindings/flow.cxx b/src/bindings/flow.cxx index 21d3b89..8108600 100644 --- a/src/bindings/flow.cxx +++ b/src/bindings/flow.cxx @@ -2,6 +2,7 @@ #include "ndarray.hxx" #include "bioimage_cpp/array_view.hxx" +#include "bioimage_cpp/detail/finite.hxx" #include "bioimage_cpp/detail/grid.hxx" #include "bioimage_cpp/flow/flow_density.hxx" @@ -91,16 +92,12 @@ DensityArray compute_flow_density_t( throw std::invalid_argument("number_of_threads must be >= 1"); } // Single authoritative finiteness check over the (contiguous) flow buffer; - // the Python wrapper deliberately does not repeat this scan. + // the Python wrapper deliberately does not repeat this scan. It runs below, + // inside the GIL release, in parallel with the caller's thread count. std::size_t flow_size = D; for (std::size_t axis = 0; axis < D; ++axis) { flow_size *= flow.shape(axis + 1); } - for (std::size_t index = 0; index < flow_size; ++index) { - if (!std::isfinite(flow.data()[index])) { - throw std::invalid_argument("flow must contain only finite values"); - } - } std::vector out_shape(D); std::vector view_shape(D); @@ -119,6 +116,13 @@ DensityArray compute_flow_density_t( { nb::gil_scoped_release release; + // `flow` (the ndarray argument) keeps the buffer alive for the whole + // call, so reading it without the GIL is safe. + if (!bioimage_cpp::detail::all_finite( + flow.data(), flow_size, static_cast(number_of_threads) + )) { + throw std::invalid_argument("flow must contain only finite values"); + } if constexpr (D == 2) { flow::compute_flow_density_2d( flow_view, @@ -151,6 +155,12 @@ DensityArray compute_flow_density_t( } // namespace void bind_flow(nb::module_ &m) { + m.def( + "_flow_trace_backend", + &flow::trace_backend, + "Return the internal flow tracing backend selected for the default path " + "('fma' or 'scalar')." + ); m.def( "_compute_flow_density_2d_float32", &compute_flow_density_t<2>, diff --git a/src/bioimage_cpp/flow/_flow.py b/src/bioimage_cpp/flow/_flow.py index 7cd9989..d875401 100644 --- a/src/bioimage_cpp/flow/_flow.py +++ b/src/bioimage_cpp/flow/_flow.py @@ -89,9 +89,19 @@ def compute_flow_density( Optional physical spacing. For 3D data and scalar ``sigma``, smoothing uses ``sigma / spacing`` per axis, matching the reference convention. number_of_threads: - Number of threads used for the particle-tracing iteration. The final - density scatter and the (optional) Gaussian smoothing are not - parallelized here. Results are deterministic regardless of the value. + Number of threads used for the particle-tracing iteration and for + validating the flow field. The final density scatter and the + (optional) Gaussian smoothing are not parallelized here. Results are + deterministic regardless of the value. + + Notes + ----- + On x86-64 CPUs with AVX and FMA, the default settings (``method="rk2"``, + ``tol > 0``, ``restrict_to_mask=True``) use a vectorized tracer that first + repacks the flow field into channel-last storage. This allocates a + temporary buffer of roughly the size of ``flow`` for the duration of the + call. Set the environment variable ``BIOIMAGE_CPP_FLOW_FORCE_SCALAR=1`` to + use the portable scalar tracer instead. Returns ------- diff --git a/src/cpp/flow/flow_density_fma.cxx b/src/cpp/flow/flow_density_fma.cxx index f321d35..25e5779 100644 --- a/src/cpp/flow/flow_density_fma.cxx +++ b/src/cpp/flow/flow_density_fma.cxx @@ -1,10 +1,504 @@ #include "bioimage_cpp/flow/flow_density.hxx" +#include "bioimage_cpp/detail/force_inline.hxx" +#include "bioimage_cpp/detail/profile.hxx" +#include "bioimage_cpp/detail/threading.hxx" + #if !defined(BIOIMAGE_FLOW_FMA_DISPATCH) #error "flow_density_fma.cxx must only be built with BIOIMAGE_FLOW_FMA_DISPATCH" #endif +#if !defined(__AVX__) || !(defined(__FMA__) || defined(_MSC_VER)) +#error "flow_density_fma.cxx must be compiled with AVX and FMA enabled (-mavx -mfma or /arch:AVX2)" +#endif + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +// Runtime-dispatched tracer for the default path (RK2, convergence check, mask +// restriction), selected by try_trace_all_fma when the CPU supports AVX+FMA. +// +// Packed-channel kernel. The flow field is repacked once into channel-last +// storage (D floats per voxel, no padding between channels) with the grid +// padded by one duplicated voxel along every axis. Every interpolation corner +// is then one 16-byte load holding all channels, the per-channel scalar lerps +// become packed operations, and the duplicated border makes the upper corner +// index always valid (the scalar kernel's nearest-boundary behaviour without a +// clamp). The particle position lives in one 4-lane register, so the RK2 +// midpoint, clipping, the convergence test and the mask lookup are packed too, +// while integer address arithmetic runs in general-purpose registers to keep it +// out of the FP register file. The per-lane arithmetic is the same IEEE +// single-precision operation sequence as the scalar kernel (fused lerp, fused +// midpoint, unfused position update). +// +// The scalar-FMA instantiation of the header kernel is kept as fallback when +// the packed buffer cannot be allocated or its offsets exceed int32. namespace bioimage_cpp::flow::detail { +namespace { + +constexpr std::ptrdiff_t kPackedTailSlack = 4; // floats past the last voxel + +template +struct PackedGrid { + const float *data = nullptr; + __m128 upper_f; // shape-1 per axis (float), 0 in unused lanes + __m128 lane_mask; // all-ones in the D used lanes + std::array pstride{}; // float-element strides of the padded packed grid + std::array mstride{}; // element strides of the mask grid +}; + +// Repack channel-first flow into channel-last storage on a grid padded by one +// voxel per axis (padding replicates the last voxel/row/plane). Returns false, +// and the caller falls back to the scalar kernel, if offsets would not fit +// int32 lanes or the buffer cannot be allocated. +template +bool make_packed_grid( + const std::array &channels, + const GridLayout &grid, + const std::size_t n_threads, + std::unique_ptr &storage, + PackedGrid &out +) { + static_assert(D == 2 || D == 3, "packed tracer supports 2D and 3D"); + constexpr std::ptrdiff_t int_max = std::numeric_limits::max(); + constexpr std::ptrdiff_t channels_per_voxel = static_cast(D); + + std::array padded{}; + std::ptrdiff_t n_padded = 1; + for (std::size_t axis = 0; axis < D; ++axis) { + padded[axis] = grid.shape[axis] + 1; + n_padded *= padded[axis]; + } + const std::ptrdiff_t total = channels_per_voxel * n_padded + kPackedTailSlack; + std::array pstride{}; + pstride[D - 1] = channels_per_voxel; + for (std::size_t axis = D - 1; axis > 0; --axis) { + pstride[axis - 1] = pstride[axis] * padded[axis]; + } + if (total > int_max) { + return false; + } + for (std::size_t axis = 0; axis < D; ++axis) { + if (grid.strides[axis] > int_max) { + return false; + } + } + try { + storage.reset(new float[static_cast(total)]); + } catch (const std::bad_alloc &) { + return false; + } + float *packed = storage.get(); + for (std::ptrdiff_t i = total - kPackedTailSlack; i < total; ++i) { + packed[i] = 0.0f; + } + + // Copy one source row (clamped y/z) into a padded row and duplicate its + // last voxel. Rows are the unit of work; planes/rows beyond the source + // extent replicate the last source plane/row. + const std::ptrdiff_t width = grid.shape[D - 1]; + const auto fill_row = [&](float *o, const std::ptrdiff_t src_offset) { + if constexpr (D == 3) { + const float *c0 = channels[0] + src_offset; + const float *c1 = channels[1] + src_offset; + const float *c2 = channels[2] + src_offset; + for (std::ptrdiff_t x = 0; x < width; ++x) { + o[3 * x] = c0[x]; + o[3 * x + 1] = c1[x]; + o[3 * x + 2] = c2[x]; + } + o[3 * width] = c0[width - 1]; + o[3 * width + 1] = c1[width - 1]; + o[3 * width + 2] = c2[width - 1]; + } else { + const float *c0 = channels[0] + src_offset; + const float *c1 = channels[1] + src_offset; + for (std::ptrdiff_t x = 0; x < width; ++x) { + o[2 * x] = c0[x]; + o[2 * x + 1] = c1[x]; + } + o[2 * width] = c0[width - 1]; + o[2 * width + 1] = c1[width - 1]; + } + }; + + if constexpr (D == 3) { + const std::ptrdiff_t n_planes = padded[0]; + const auto pack_threads = ::bioimage_cpp::detail::normalize_thread_count( + n_threads, static_cast(n_planes) + ); + ::bioimage_cpp::detail::parallel_for_chunks( + pack_threads, static_cast(n_planes), + [&](const std::size_t, const std::size_t begin, const std::size_t end) { + for (std::size_t zp = begin; zp < end; ++zp) { + const std::ptrdiff_t z = std::min(static_cast(zp), grid.shape[0] - 1); + for (std::ptrdiff_t yp = 0; yp < padded[1]; ++yp) { + const std::ptrdiff_t y = std::min(yp, grid.shape[1] - 1); + fill_row( + packed + static_cast(zp) * pstride[0] + yp * pstride[1], + z * grid.strides[0] + y * grid.strides[1] + ); + } + } + } + ); + } else { + const std::ptrdiff_t n_rows = padded[0]; + const auto pack_threads = ::bioimage_cpp::detail::normalize_thread_count( + n_threads, static_cast(n_rows) + ); + ::bioimage_cpp::detail::parallel_for_chunks( + pack_threads, static_cast(n_rows), + [&](const std::size_t, const std::size_t begin, const std::size_t end) { + for (std::size_t yp = begin; yp < end; ++yp) { + const std::ptrdiff_t y = std::min(static_cast(yp), grid.shape[0] - 1); + fill_row(packed + static_cast(yp) * pstride[0], y * grid.strides[0]); + } + } + ); + } + + out.data = packed; + for (std::size_t axis = 0; axis < D; ++axis) { + out.pstride[axis] = static_cast(pstride[axis]); + out.mstride[axis] = static_cast(grid.strides[axis]); + } + if constexpr (D == 3) { + out.upper_f = _mm_setr_ps(grid.upper[0], grid.upper[1], grid.upper[2], 0.0f); + out.lane_mask = _mm_castsi128_ps(_mm_setr_epi32(-1, -1, -1, 0)); + } else { + out.upper_f = _mm_setr_ps(grid.upper[0], grid.upper[1], 0.0f, 0.0f); + out.lane_mask = _mm_castsi128_ps(_mm_setr_epi32(-1, -1, 0, 0)); + } + return true; +} + +// lo + f * (hi - lo) with one rounding, the same contraction the scalar +// kernel's lerp compiles to under -mfma. +BIOIMAGE_FORCE_INLINE __m128 lerp(const __m128 lo, const __m128 hi, const __m128 f) { + return _mm_fmadd_ps(f, _mm_sub_ps(hi, lo), lo); +} + +// Lower corner (by truncation: position lanes are clipped and nonnegative, so +// this equals floor, as in the scalar kernel), fractional part, and the flat +// float offset of the lower corner in the padded packed grid. The integer +// work runs in general-purpose registers. +struct Stencil { + __m128 frac; + int base; +}; + +template +BIOIMAGE_FORCE_INLINE Stencil make_stencil(const PackedGrid &g, const __m128 p) { + const __m128i i0 = _mm_cvttps_epi32(p); + const __m128 frac = _mm_sub_ps(p, _mm_cvtepi32_ps(i0)); + int base = _mm_cvtsi128_si32(i0) * g.pstride[0] + _mm_extract_epi32(i0, 1) * g.pstride[1]; + if constexpr (D == 3) { + base += _mm_extract_epi32(i0, 2) * g.pstride[2]; + } + return Stencil{frac, base}; +} + +// All channels of the flow at packed position p = [z, y, x, 0] (3D) or +// [y, x, 0, 0] (2D). Lanes beyond D hold finite junk (the next voxel's first +// channel or tail slack); the caller cleans them where it matters. +// 3D: nested trilinear lerps, the same association order as the scalar kernel. +// (A weighted-sum tree with a shorter dependent chain was measured within noise +// here because the 3D loop is bound by the reorder window, not by FP register +// pressure; see the 2D sampler.) +BIOIMAGE_FORCE_INLINE __m128 sample_packed(const PackedGrid<3> &g, const __m128 p) { + const Stencil s = make_stencil<3>(g, p); + constexpr int dx = 3; + const int dy = g.pstride[1]; + const int dz = g.pstride[0]; + const float *p0 = g.data + s.base; + const float *p1 = p0 + dz; + // Each load is [c0, c1, c2, junk]: junk is the next voxel's first channel + // or the tail slack, always finite. + const __m128 c000 = _mm_loadu_ps(p0); + const __m128 c001 = _mm_loadu_ps(p0 + dx); + const __m128 c010 = _mm_loadu_ps(p0 + dy); + const __m128 c011 = _mm_loadu_ps(p0 + dy + dx); + const __m128 c100 = _mm_loadu_ps(p1); + const __m128 c101 = _mm_loadu_ps(p1 + dx); + const __m128 c110 = _mm_loadu_ps(p1 + dy); + const __m128 c111 = _mm_loadu_ps(p1 + dy + dx); + const __m128 fz = _mm_shuffle_ps(s.frac, s.frac, _MM_SHUFFLE(0, 0, 0, 0)); + const __m128 fy = _mm_shuffle_ps(s.frac, s.frac, _MM_SHUFFLE(1, 1, 1, 1)); + const __m128 fx = _mm_shuffle_ps(s.frac, s.frac, _MM_SHUFFLE(2, 2, 2, 2)); + const __m128 c00 = lerp(c000, c001, fx); + const __m128 c01 = lerp(c010, c011, fx); + const __m128 c10 = lerp(c100, c101, fx); + const __m128 c11 = lerp(c110, c111, fx); + const __m128 c0 = lerp(c00, c01, fy); + const __m128 c1 = lerp(c10, c11, fy); + return lerp(c0, c1, fz); +} + +// 2D: bilinear interpolation as a weighted sum of the four corners. The four +// weights are formed from the fractional coordinates off the load path, so the +// dependent chain after the loads is mul -> fma -> add (11 cycles) instead of +// two nested lerp levels (16 cycles). The 2D loop is bound by FP register-file +// occupancy, where this measured 12 % faster. The association order differs +// from the scalar kernel's nested lerps; the value is the same bilinear +// interpolation. +BIOIMAGE_FORCE_INLINE __m128 sample_packed(const PackedGrid<2> &g, const __m128 p) { + const Stencil s = make_stencil<2>(g, p); + const float *r0p = g.data + s.base; + // Row loads: [c0(x0), c1(x0), c0(x1), c1(x1)]; x1 = x0 + 1 is valid thanks + // to the padding column. + const __m128 r0 = _mm_loadu_ps(r0p); + const __m128 r1 = _mm_loadu_ps(r0p + g.pstride[0]); + const __m128 f = s.frac; // [fy, fx, 0, 0] + const __m128 omf = _mm_sub_ps(_mm_set1_ps(1.0f), f); // [gy, gx, 1, 1] + const __m128 vx = _mm_shuffle_ps(omf, f, _MM_SHUFFLE(1, 1, 1, 1)); // [gx, gx, fx, fx] + const __m128 w0 = _mm_mul_ps(vx, _mm_shuffle_ps(omf, omf, _MM_SHUFFLE(0, 0, 0, 0))); // row y0 + const __m128 w1 = _mm_mul_ps(vx, _mm_shuffle_ps(f, f, _MM_SHUFFLE(0, 0, 0, 0))); // row y1 + const __m128 t = _mm_fmadd_ps(r1, w1, _mm_mul_ps(r0, w0)); // [x0 terms | x1 terms] + return _mm_add_ps(t, _mm_movehl_ps(t, t)); // lanes 0,1 = the sum +} + +struct StepConstants { + __m128 zero; + __m128 half; + __m128 half_dt; + __m128 dt; + __m128 tol; + __m128 sign_bit; +}; + +inline StepConstants make_step_constants(const float dt, const float tol) { + return StepConstants{ + _mm_setzero_ps(), + _mm_set1_ps(0.5f), + _mm_set1_ps(0.5f * dt), + _mm_set1_ps(dt), + _mm_set1_ps(tol), + _mm_set1_ps(-0.0f), + }; +} + +constexpr int kStillAlive = -1; + +// One RK2 step of a single particle with convergence and mask termination. +// Returns kStillAlive, or the ActiveTraceStats exit reason if tracing stops. +template +BIOIMAGE_FORCE_INLINE int step_packed( + __m128 &pos, const PackedGrid &g, const StepConstants &c, const std::uint8_t *mask +) { + // The junk lane of the first sample is neutralised by the midpoint clip + // (upper bound 0 in unused lanes); the second sample is masked explicitly + // because it feeds the convergence test, the bounds test and the position. + __m128 step = sample_packed(g, pos); + __m128 mid = _mm_fmadd_ps(c.half_dt, step, pos); + mid = _mm_min_ps(_mm_max_ps(mid, c.zero), g.upper_f); + step = _mm_and_ps(sample_packed(g, mid), g.lane_mask); + + // Scalar kernel: converged iff max_axis |dt*step| < tol, i.e. no lane has + // |dt*step| >= tol. + const __m128 disp = _mm_mul_ps(c.dt, step); + const __m128 abs_disp = _mm_andnot_ps(c.sign_bit, disp); + if (_mm_movemask_ps(_mm_cmpge_ps(abs_disp, c.tol)) == 0) { + return ActiveTraceStats::Converged; + } + + const __m128 proposed = _mm_add_ps(pos, disp); + const __m128 outside = + _mm_or_ps(_mm_cmplt_ps(proposed, c.zero), _mm_cmpgt_ps(proposed, g.upper_f)); + if (_mm_movemask_ps(outside) != 0) { + return ActiveTraceStats::LeftMask; + } + // Round half up to the mask voxel (proposed is in-bounds and nonnegative + // here, so the truncating conversion equals floor(x + 0.5)). + const __m128i rounded = _mm_cvttps_epi32(_mm_add_ps(proposed, c.half)); + int flat = _mm_cvtsi128_si32(rounded) * g.mstride[0] + _mm_extract_epi32(rounded, 1) * g.mstride[1]; + if constexpr (D == 3) { + flat += _mm_extract_epi32(rounded, 2) * g.mstride[2]; + } + if (mask[flat] == 0) { + return ActiveTraceStats::LeftMask; + } + pos = proposed; + return kStillAlive; +} + +template +BIOIMAGE_FORCE_INLINE __m128 load_position(const std::array &p) { + if constexpr (D == 3) { + return _mm_setr_ps(p[0], p[1], p[2], 0.0f); + } else { + return _mm_setr_ps(p[0], p[1], 0.0f, 0.0f); + } +} + +template +BIOIMAGE_FORCE_INLINE void store_position(std::array &p, const __m128 v) { + alignas(16) float tmp[4]; + _mm_store_ps(tmp, v); + for (std::size_t axis = 0; axis < D; ++axis) { + p[axis] = tmp[axis]; + } +} + +template +BIOIMAGE_FORCE_INLINE void for_each_lane(F &&f, std::index_sequence) { + (f(std::integral_constant{}), ...); +} + +// Trace K consecutive particles in lockstep (see trace_particle_block in the +// header). Lanes are addressed with compile-time indices so their state stays +// in registers. +template +void trace_block_packed( + std::array *positions, + const PackedGrid &g, + const StepConstants &c, + const std::uint8_t *mask, + const std::size_t n_iter, + ActiveTraceStats &stats +) { + __m128 pos[K]; + bool alive[K]; + [[maybe_unused]] std::size_t lane_steps[K] = {}; + [[maybe_unused]] int lane_exit[K] = {}; + for_each_lane([&](auto k) { + pos[k] = load_position(positions[k]); + alive[k] = true; + if constexpr (ActiveTraceStats::enabled) { + lane_exit[k] = ActiveTraceStats::MaxIter; + } + }, std::make_index_sequence{}); + + for (std::size_t iter = 0; iter < n_iter; ++iter) { + bool any_alive = false; + for_each_lane([&](auto k) { + if (!alive[k]) { + return; + } + if constexpr (ActiveTraceStats::enabled) { + ++lane_steps[k]; + } + const int exit_reason = step_packed(pos[k], g, c, mask); + if (exit_reason != kStillAlive) { + alive[k] = false; + if constexpr (ActiveTraceStats::enabled) { + lane_exit[k] = exit_reason; + } + return; + } + any_alive = true; + }, std::make_index_sequence{}); + if (!any_alive) { + break; + } + } + + for_each_lane([&](auto k) { + store_position(positions[k], pos[k]); + if constexpr (ActiveTraceStats::enabled) { + stats.record(lane_steps[k], lane_exit[k]); + } + }, std::make_index_sequence{}); +} + +template +void trace_all_packed( + std::vector> &positions, + const PackedGrid &g, + const std::uint8_t *mask, + const std::size_t n_threads, + const std::size_t n_iter, + const float dt, + const float tol +) { + const StepConstants c = make_step_constants(dt, tol); + ActiveTraceStats null_stats{}; + std::vector per_thread_stats; + if constexpr (ActiveTraceStats::enabled) { + per_thread_stats.assign(n_threads, ActiveTraceStats(n_iter)); + } + ::bioimage_cpp::detail::parallel_for_chunks( + n_threads, positions.size(), + [&](const std::size_t thread_id, const std::size_t begin, const std::size_t end) { + ActiveTraceStats *stats = &null_stats; + (void)thread_id; +#ifdef BIOIMAGE_PROFILE + stats = &per_thread_stats[thread_id]; + const auto chunk_start = std::chrono::steady_clock::now(); +#endif + std::size_t i = begin; + if constexpr (K > 1) { + for (; i + K <= end; i += K) { + trace_block_packed(&positions[i], g, c, mask, n_iter, *stats); + } + } + for (; i < end; ++i) { + trace_block_packed(&positions[i], g, c, mask, n_iter, *stats); + } +#ifdef BIOIMAGE_PROFILE + stats->seconds += std::chrono::duration( + std::chrono::steady_clock::now() - chunk_start + ).count(); +#endif + } + ); + if constexpr (ActiveTraceStats::enabled) { + ActiveTraceStats::report(per_thread_stats, n_iter, true); + } +} + +// Lanes traced in lockstep per group. K = 4 measured 1-6 % faster than K = 3 +// on the registered fixtures and every adversarial case (2026-09-02, Zen 3); +// K = 6 was parity, K = 2 and K = 1 lose. +constexpr std::size_t kPackedLanes = 4; + +template +void trace_all_fma( + std::vector> &positions, + const std::array &channels, + const GridLayout &grid, + const std::uint8_t *mask, + const std::size_t n_threads, + const std::size_t n_iter, + const float dt, + const float tol +) { + BIOIMAGE_PROFILE_INIT(profiler); + { + std::unique_ptr storage; + PackedGrid packed{}; + bool ready = false; + { + BIOIMAGE_PROFILE_SCOPE(profiler, "fma_pack"); + ready = make_packed_grid(channels, grid, n_threads, storage, packed); + } + if (ready) { + BIOIMAGE_PROFILE_SCOPE(profiler, "fma_trace_packed"); + trace_all_packed(positions, packed, mask, n_threads, n_iter, dt, tol); + BIOIMAGE_PROFILE_REPORT_NAMED(profiler, "[bioimage profile: flow fma]"); + return; + } + } + // Fallback: packed offsets would not fit int32 or the buffer could not be + // allocated. Same header kernel as the portable path, compiled with FMA. + { + BIOIMAGE_PROFILE_SCOPE(profiler, "fma_trace_scalar"); + trace_all( + positions, channels, grid, mask, n_threads, n_iter, dt, tol + ); + } + BIOIMAGE_PROFILE_REPORT_NAMED(profiler, "[bioimage profile: flow fma]"); +} + +} // namespace void trace_all_fma_2d( std::vector> &positions, @@ -16,9 +510,7 @@ void trace_all_fma_2d( const float dt, const float tol ) { - trace_all<2, true, true, true, true>( - positions, channels, grid, mask, n_threads, n_iter, dt, tol - ); + trace_all_fma<2>(positions, channels, grid, mask, n_threads, n_iter, dt, tol); } void trace_all_fma_3d( @@ -31,9 +523,7 @@ void trace_all_fma_3d( const float dt, const float tol ) { - trace_all<3, true, true, true, true>( - positions, channels, grid, mask, n_threads, n_iter, dt, tol - ); + trace_all_fma<3>(positions, channels, grid, mask, n_threads, n_iter, dt, tol); } } // namespace bioimage_cpp::flow::detail diff --git a/tests/test_flow.py b/tests/test_flow.py index afc68d5..384a689 100644 --- a/tests/test_flow.py +++ b/tests/test_flow.py @@ -235,6 +235,72 @@ def test_flow_rejects_non_finite_values(): bic.flow.compute_flow_density(flow, np.ones((3, 3), bool)) +@pytest.mark.parametrize("value", [np.nan, np.inf, -np.inf]) +@pytest.mark.parametrize("shape", [(3, 3), (2, 3, 4), (1, 1), (1, 1, 1)]) +@pytest.mark.parametrize("threads", [1, 4]) +def test_flow_rejects_non_finite_values_anywhere(value, shape, threads): + ndim = len(shape) + flow = np.zeros((ndim,) + shape, dtype=np.float32) + # last element of the last channel: the position a chunked scan is most + # likely to miss + flow[(ndim - 1,) + tuple(s - 1 for s in shape)] = value + with pytest.raises(ValueError, match="finite"): + bic.flow.compute_flow_density( + flow, np.ones(shape, bool), number_of_threads=threads + ) + flow[...] = 0.0 + flow[(0,) + (0,) * ndim] = value + with pytest.raises(ValueError, match="finite"): + bic.flow.compute_flow_density( + flow, np.ones(shape, bool), number_of_threads=threads + ) + + +def test_large_finite_flow_is_accepted_multithreaded(): + # Large enough that the parallel finiteness scan actually fans out. + shape = (8, 384, 384) + flow = np.full((3,) + shape, 0.25, dtype=np.float32) + density = bic.flow.compute_flow_density( + flow, np.ones(shape, bool), n_iter=1, number_of_threads=4 + ) + assert density.sum() == np.prod(shape) + + +def test_forced_scalar_matches_fma_backend(monkeypatch): + monkeypatch.delenv("BIOIMAGE_CPP_FLOW_FORCE_SCALAR", raising=False) + if bic._core._flow_trace_backend() != "fma": + pytest.skip("FMA flow backend is not available on this CPU/build") + + rng = np.random.default_rng(2026) + cases = [] + for shape in [(48, 64), (6, 24, 32)]: + ndim = len(shape) + flow = rng.normal(scale=0.7, size=(ndim,) + shape).astype(np.float32) + mask = rng.random(shape) > 0.3 + cases.append((flow, mask)) + + def run_all(): + return [ + bic.flow.compute_flow_density(flow, mask, number_of_threads=threads) + for flow, mask in cases + for threads in (1, 4) + ] + + automatic = run_all() + monkeypatch.setenv("BIOIMAGE_CPP_FLOW_FORCE_SCALAR", "1") + assert bic._core._flow_trace_backend() == "scalar" + scalar = run_all() + monkeypatch.setenv("BIOIMAGE_CPP_FLOW_FORCE_SCALAR", "0") + assert bic._core._flow_trace_backend() == "fma" + + # The two backends may differ in floating-point contraction, so the bar is + # particle conservation plus near-identical densities, not bitwise equality. + for got, expected in zip(automatic, scalar, strict=True): + assert got.sum() == expected.sum() + differing = np.count_nonzero(got != expected) + assert differing <= max(1, got.size // 1000) + + def test_rk2_runs(): rng = np.random.default_rng(0) flow = rng.normal(scale=0.3, size=(2, 12, 12)).astype(np.float32) From 126df39bff499be1032a61d4b00e7964caaa8578 Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Thu, 3 Sep 2026 09:08:22 +0200 Subject: [PATCH 2/2] Reproduce on laptop CPU, add prompts for MAC and WIN optimization --- development/flow/PERFORMANCE_NOTES.md | 41 +++++++++- development/flow/PORTING_CAMPAIGNS.md | 113 ++++++++++++++++++++++++++ 2 files changed, 153 insertions(+), 1 deletion(-) create mode 100644 development/flow/PORTING_CAMPAIGNS.md diff --git a/development/flow/PERFORMANCE_NOTES.md b/development/flow/PERFORMANCE_NOTES.md index 4416608..234edc7 100644 --- a/development/flow/PERFORMANCE_NOTES.md +++ b/development/flow/PERFORMANCE_NOTES.md @@ -334,7 +334,8 @@ All A/B numbers come from `development/flow/paired_bench.py` (fresh pinned subprocess per measurement, A B B A per repeat, prebuilt `.so` files swapped through `sys.modules`, never rebuilt between sides) and from the kernel-only counter harness `development/flow/perf_kernel.py`. This machine cannot measure -beyond two physical cores; 8T-class numbers remain to be taken elsewhere. +beyond two physical cores; the 4T/8T numbers come from the Tiger Lake +reproduction below. ### What changed @@ -476,6 +477,44 @@ Non-kernel: the parallel, GIL-free finiteness check cuts the `n_iter=0` call from ~40 ms to a few ms on the 3D fixture (was 3.8 % at 1T, 7 % at 2T of the previous kernel's runtime, growing with thread count). +### Reproduction on Tiger Lake (2026-09-02) + +Same harness, same flags, same pair of commits (`492ff9b` vs `54fd4a9`) on the +i7-1185G7 laptop (governor `powersave`, 4 cores / 8 threads). Calibration +(base vs base) read `noise` on both fixtures. Best-of-subprocess minima, all +rows `significant` unless noted: + +``` +case threads previous packed change Zen 3 parity +fixture3d 1 1.5545 s 1.0342 s -33.5 % -39.4 identical +fixture2d 1 0.1664 s 0.1140 s -31.5 % -37.3 identical +fixture3d 4 phys 0.4728 s 0.3277 s -30.7 % -41.6 identical +fixture3d 8 0.3588 s 0.2610 s -27.2 % n/a identical +stripes 1/4 1 0.7611 s 0.6532 s -14.2 % -17.1 identical +stripes 3/4 1 1.9156 s 1.1343 s -40.8 % -40.0 identical +random s=1 1 2.1683 s 1.3579 s -37.4 % -42.1 identical +random s=10 1 2.2315 s 1.5259 s -31.6 % -44.1 identical +random s=5 1 1.5852 s 1.0519 s -33.6 % -40.7 identical +random s=20 1 2.0508 s 1.6521 s -19.4 % -40.4 identical +random s=40 1 1.3517 s 1.1559 s -14.5 % -33.5 identical +2D random 10 1 0.2586 s 0.1712 s -33.8 % -36.1 differs (2D tree) +``` + +The gain carries over but is about 5 points smaller than on Zen 3 for the +fixtures and structured cases, and clearly smaller for the high-scale random +flows (s=20: -19 % vs -40 %; s=40: -15 % vs -34 %), so those two rows should +not be quoted as general numbers. The 2D fixture at 4T/8T finishes in under +45 ms; the per-side spread there is 30-46 %, so the harness reports `noise` +although the minima dropped 27-30 %. + +Correctness on this host: `differential_check.py` is 400/400 bitwise +identical both previous-vs-packed and packed-vs-`FORCE_SCALAR` (Zen 3 saw +399/400 for the latter); the `check_flow_density.py` gate is unchanged; the +full suite is 1466 passed. Note that `BIOIMAGE_CPP_FLOW_FORCE_SCALAR=1` +bypasses the whole FMA translation unit and runs the portable SSE2 kernel +(1.93 s on the 3D fixture at 1T), which is slower than the previous kernel's +scalar-FMA default (1.55 s); this is expected, not a regression. + ### Codegen trap: translation-unit growth disables inlining While the FMA TU held many template instantiations (K and layout variants for diff --git a/development/flow/PORTING_CAMPAIGNS.md b/development/flow/PORTING_CAMPAIGNS.md new file mode 100644 index 0000000..aef9798 --- /dev/null +++ b/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 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. +```