Optimize hot kernels and L-BFGS minimizer; fix correctness issues

CMakeLists.txt

@@ -2,7 +2,14 @@ cmake_minimum_required(VERSION 3.14)
 if(POLICY CMP0135)
     cmake_policy(SET CMP0135 NEW)
 endif()
-project(my-fantastic-simulation-code LANGUAGES C CXX VERSION 0.1)
+project(qew LANGUAGES C CXX VERSION 0.1)
+
+# Bandwidth-bound HPC code: an accidental no-flags build (-O0) is never what
+# anyone wants, so default to Release when the user does not choose.
+if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
+    set(CMAKE_BUILD_TYPE Release CACHE STRING "Build type" FORCE)
+    message(STATUS "CMAKE_BUILD_TYPE not set, defaulting to Release")
+endif()
 
 # Kokkos 5.x requires C++20 or newer.
 set(CMAKE_CXX_STANDARD 20)

docs/fix_plan.md

@@ -0,0 +1,176 @@
+# Implementation plan for `optimization_notes.md`
+
+> **Status (June 2026): all phases 0–4 are implemented on this branch.**
+> Phase 1 was settled by the paper (arXiv:2410.21838 eq. (10):
+> `lambda_p/xi = (Gamma/U0)^{2/3}`, varied by scanning Gamma at fixed
+> disorder): all drivers now take the dimensionless ratio as
+> `pinning_length_over_xi` and share `Gamma = A (lambda_p/xi)^{3/2}`, with
+> per-driver mapping tests; line tensions are numerically unchanged.
+> Measured on the benchmark (Serial): Phase 2 took the end-to-end
+> preconditioned L-BFGS relaxation (n = 4096, arclength) from 5.73 s / 12893
+> objective evals to 4.77 s / 10424 at an identical 2470 iterations; Phase 3
+> took it to 3.55 s (38% total) and the hot kernels (n = 262144) from
+> grad 33.5 / dyn 52.4 ms to 26.5 / 34.4 ms, with the L-BFGS trajectory
+> bit-identical throughout. Phase 4's RK change removes one of two
+> reductions per sweep; the Langevin RNG chunking is serial-neutral and its
+> GPU benefit still needs the CUDA runner (see below). Phase 5 remains
+> won't-fix unless profiling says otherwise.
+
+Companion to [`optimization_notes.md`](optimization_notes.md); item numbers
+below refer to that document. Work is grouped into five phases that are
+independently land-able, ordered so that correctness fixes go first and each
+perf change is benchmarked in isolation. Every phase ends with the full
+`ctest` suite on a Serial **and** an OpenMP build (the threaded run doubles as
+the race check), plus a before/after run of `executables/benchmark` for the
+perf phases.
+
+Ground rules:
+
+- The gradient-consistency test and the golden-file tests are the correctness
+  gates; they must pass unchanged (same tolerances) in every phase.
+- Tests that assert `n_func_evals` counts will be updated *deliberately* in
+  Phase 2, in the same commit as the behaviour change, with the new expected
+  count derived by hand.
+- One commit (or small commit series) per phase, so a regression bisects to a
+  phase.
+
+## Phase 0 — Safety & build hygiene (items 1, 3, 4)
+
+Small, riskless, do first.
+
+1. **Odd-`n` guard** (`rosso_krauth.h`): in `rk_block()`, before launching
+   kernels, `throw std::invalid_argument` (host code) when `n % 2 != 0`,
+   citing the red-black neighbour conflict. A 3-colour sweep for odd `n` is
+   *not* worth the complexity until someone needs odd lines. Add a death/throw
+   test.
+2. **Integer ramp counter** (`qew_runs.h::depinning_ramp`): iterate
+   `for (int k = 0; k * f_step <= f_max; ++k)` with `f = k * f_step`. The
+   default `f_step = 1.0` produces bit-identical forces, so no test or HDF5
+   key changes.
+3. **CMake**: default `CMAKE_BUILD_TYPE=Release` when unset (with a STATUS
+   message); rename the project to `qew`. Leave `hello.cpp` /
+   `kokkos_smoketest.cpp` in place for now — they are referenced by tests and
+   removing them is cosmetic churn; fold into a later cleanup if desired.
+
+## Phase 1 — Decision point: line-tension convention (item 2)
+
+Blocked on information outside the repo: the `../Overleaf/qEW/qew_*.py`
+sources. **Needs Lars to confirm** which convention each Python driver uses.
+Then, in one commit:
+
+- Rename the field in whichever config struct holds the ratio to
+  `pinning_length_over_xi` (or convert at the call site so all three drivers
+  use one formula).
+- Add a tiny unit test per driver pinning the `pinning_length →
+  line_tension` mapping, so the convention can never silently diverge again.
+
+Nothing later depends on this; it can land any time.
+
+## Phase 2 — L-BFGS overhead (items 5–7)
+
+All inside `lbfgs.h`, in three commits, benchmarking `qew_static`-shaped runs
+(preconditioned, n = 4096) after each:
+
+1. **Carry the objective forward** (item 5). Hoist one
+   `phi = objective(h, ...)` above the iteration loop; on line-search
+   acceptance set `phi = phi_accepted`; delete the per-iteration `phi0`
+   recompute. The line-search-failure `continue` path restores `h = x0`, so
+   the carried value stays valid there too. Arithmetic is unchanged →
+   golden L-BFGS test passes bit-for-bit; only `n_func_evals` assertions in
+   `test_lbfgs.cpp` need recounting.
+2. **Fuse trial step into the objective** (item 6). Add a local helper
+   `trial_objective(x0, d, alpha, p, noise)` — the `objective()` reduce with
+   `h(i)` replaced by `x0(i) + alpha * d(i)` (and the same for the `ip`
+   neighbour, which is the identical FP expression, so results stay
+   bit-identical). Commit `h` with a single `parallel_for` only on
+   acceptance; the `deep_copy(x0, h)` save and `deep_copy(h, x0)` undo both
+   disappear. Keep `objective()` itself untouched (public API, used by
+   tests/benchmark).
+3. **Reduction fusion + history swap** (item 7).
+   - One `parallel_reduce` with three reducers computing `sy`, `ss`, `yy`
+     while also writing `s_tmp`/`y_tmp` (replaces one `parallel_for` + three
+     reduces).
+   - `std::swap(S[row], s_tmp)` / `std::swap(Y[row], y_tmp)` instead of
+     `deep_copy` (safe: the FFT plans are bound to `q`/`d` only).
+   - Fold `q -= a*y` into the subsequent dot where profiling on a GPU build
+     shows the sync count matters; otherwise skip — this is the only
+     speculative sub-item, and it slightly obfuscates the two-loop recursion.
+
+## Phase 3 — Hot-kernel work reduction (items 8, 9, 10)
+
+The order matters: land 10 (mechanical, bit-identical) first so the bigger
+rewrites diff cleanly against it.
+
+1. **Modulo removal** (item 10, `filtered_noise.h::sample`,
+   `qew_model.h`, `rosso_krauth.h`): replace `(i±1) % n` stencil wraps with
+   compare-and-adjust; in `sample()`, clamp the `gx == nx` rounding edge case
+   explicitly (`if (ix >= nx) { ix -= nx; gx -= nx; }` before computing
+   `fx`), then build the four neighbour indices with conditional wraps.
+   Same arithmetic → golden noise-interpolation test must pass bit-for-bit.
+2. **Value-only / dy-only sampling** (item 9). Add to `DeviceNoise`:
+   `sample_value(x, y)` (no `dw`, no derivative accumulators) and keep
+   `sample()` as-is. In `qew_model.h`, route calls through C++20
+   `requires`-detected helpers (`noise_value(noise, x, y)` falls back to
+   `noise.sample(x, y).v` for the analytic test stubs, so no test functor
+   has to change). Split `bspline_weights` into weights-only and
+   weights+derivatives variants. `objective()` switches to the value path.
+3. **Edge-pass gradient** (item 8). Introduce
+   `struct GradientWorkspace { View1D edge; }` (sized lazily to `n`).
+   `gradient(h, p, noise, grad, ws)` becomes two kernels:
+   - edge kernel: for edge `(i, i+1)` store
+     `edge(i) = 0.5 * dv_dy(midpoint)` — plus, for the arclength model, the
+     normalized slope `dh/sqrt(1+dh²)`, packed into a second view or a 2-D
+     workspace, since both adjacent sites reuse it;
+   - combine kernel: `grad(i) = dx * (line_force(from edge slopes or h) +
+     edge(im) + edge(i) - f)`.
+   Keep a back-compat `gradient(h, p, noise, grad)` overload that allocates a
+   workspace internally (used by tests); `lbfgs_minimize`, `fire_minimize`,
+   and `LangevinDynamics` own a persistent workspace. Floating-point sums
+   re-associate slightly (`0.5*(sl+sr)` becomes `0.5*sl + 0.5*sr`), so the
+   gradient-consistency test (tolerance-based) is the gate, and the golden
+   L-BFGS comparison tolerance may need a documented epsilon bump — verify
+   it stays within the existing tolerance first.
+   `rk_site_velocity` is intentionally untouched (single-site re-evaluation
+   cannot share edges).
+
+Acceptance: benchmark `grad` ms/eval roughly halves on OpenMP; `obj`
+improves measurably from item 9; no test tolerance loosened beyond what is
+justified in the commit message.
+
+## Phase 4 — Rosso–Krauth & dynamics (items 11, 12)
+
+1. **Drop the per-sweep mean reduction** (item 11): track
+   `cum_adv += total_adv` and compute `res.mean_h = mean_h0 + cum_adv / n`.
+   Bit-identical apart from summation order of the mean; the RK tests
+   (washboard `f_c`, runaway) are tolerance-based and unaffected.
+2. **Langevin RNG** (item 12): measure first. On the OpenMP build the pool
+   is per-thread and fine; the concern is GPU lock contention. Run the
+   benchmark's `dyn` timing on the CUDA CI runner; only if `dyn − grad`
+   overhead is significant, switch to a pre-generated normal-deviate view
+   filled by a dedicated kernel where each thread owns one pool state for the
+   whole row-chunk (avoids per-element lock churn), fused with the update
+   step. The equipartition test gates correctness; exact trajectories may
+   change (different draw order) — that test is statistical, so no fixture
+   updates expected.
+
+## Phase 5 — Optional cleanups (items 13–15)
+
+FIRE reduction fusion, host-side noise-setup loops, `to_host_vec` memcpy.
+Do only if profiling ever shows them; otherwise close as won't-fix and note
+that in `optimization_notes.md`.
+
+## Sequencing summary
+
+| Phase | Contents | Risk | Blocked on |
+|---|---|---|---|
+| 0 | odd-`n` guard, ramp counter, CMake | none | — |
+| 1 | line-tension convention | none | Python sources (Lars) |
+| 2 | L-BFGS: carry φ, fused trial, fused reduces, view swap | low | — |
+| 3 | modulo removal → sample variants → edge-pass gradient | medium | — |
+| 4 | RK mean tracking; Langevin RNG (measure-first) | low | GPU runner for 12 |
+| 5 | FIRE / setup / copy cleanups | — | profiling evidence |
+
+Phases 0, 2, 3, 4 are mutually independent except that Phase 3's edge-pass
+should rebase on Phase 2's line-search changes (both touch the `lbfgs.h`
+call sites). Suggested landing order: 0 → 2 → 3 → 4, with 1 whenever the
+Python check happens.

docs/optimization_notes.md

@@ -0,0 +1,159 @@
+# Optimization potential & known issues
+
+Technical notes from a code review of the qEW port (June 2026). Items are
+ordered by expected impact.
+
+> **Status:** items 1–12 are addressed on this branch — see `fix_plan.md`
+> for the phase breakdown and measured numbers. Item 12's GPU benefit
+> (RNG lock relief) still needs verification on the CUDA runner; the minor
+> items 13–15 are won't-fix unless profiling says otherwise. File references point at the current `main`-line
+sources. Nothing here changes physics; everything preserves the documented
+analytic consistency between `objective()` and `gradient()` and the golden-file
+contracts.
+
+The hot path, per the benchmark (`executables/benchmark`), is the bicubic noise
+gather inside `gradient()` / `objective()` — memory-bandwidth-bound (16
+coefficient loads + ~70 FLOPs per `DeviceNoise::sample()`). Most of the wins
+below attack either (a) the number of `sample()` calls, or (b) per-iteration
+kernel-launch / device-sync overhead in the minimizers.
+
+## Bugs / correctness hazards
+
+1. **Red–black sweep races for odd `n`** — `src/rosso_krauth.h:147`
+   (`if ((i & 1) != color) return;`). On a periodic ring with odd `n`, sites
+   `0` and `n-1` have the same parity but are nearest neighbours, so the same
+   colour pass updates both while each reads the other's `h`. All current
+   configs use even `nx`, so this is latent — but nothing enforces it.
+   *Fix:* assert `n % 2 == 0` in `rk_block()` (cheap), or fall back to a
+   3-colour sweep for odd `n`.
+
+2. **Inconsistent pinning-length → line-tension map across run drivers** —
+   `src/qew_runs.h`: `static_sweep()` and `depinning_ramp()` used
+   `line_tension = pinning_length^1.5 * amplitude`, while `dynamic_anneal()`
+   used `(pinning_length / xi)^1.5 * amplitude` — the same struct-field name
+   meant two different physical quantities. **Resolved:** the paper
+   (arXiv:2410.21838, eq. (10)) controls the model through the dimensionless
+   Larkin ratio `lambda_p/xi = (Gamma/U0)^{2/3}`, so all drivers now take a
+   `pinning_length_over_xi` field and share `Gamma = A * (lambda_p/xi)^{3/2}`;
+   each driver's mapping is pinned by a unit test. HDF5 key strings are
+   unchanged and their values carry the ratio (already true for static and
+   depinning output; dynamic output previously recorded the physical length).
+
+3. **Floating-point loop counter in the force ramp** —
+   `src/qew_runs.h:214`: `for (real_t f = 0; f <= c.f_max; f += c.f_step)`.
+   Exact for the default `f_step = 1.0`, but for a non-representable step
+   (e.g. `0.1`) accumulated rounding can drop or shift the final step and
+   perturb the HDF5 key strings (`driving_force=%.4g`). *Fix:* integer index,
+   `f = k * f_step`.
+
+4. **No default `CMAKE_BUILD_TYPE`** — root `CMakeLists.txt` never sets a
+   fallback, so configuring without `-DCMAKE_BUILD_TYPE=Release` silently
+   produces an unoptimized (`-O0`) build of a bandwidth-bound HPC code. *Fix:*
+   default to `Release` when unset (warn). Cosmetic: the project is still
+   named `my-fantastic-simulation-code` (line 5), and `hello.cpp` /
+   `kokkos_smoketest.cpp` scaffolding is still compiled into the production
+   `lib` target.
+
+## Optimization potential
+
+### Minimizer (L-BFGS, `src/lbfgs.h`)
+
+5. **Redundant objective evaluation every iteration** — line 200:
+   `phi0 = objective(h, ...)` is recomputed at the start of each line search,
+   but `h` is exactly the point accepted by the previous iteration's line
+   search, whose `phi` is already known. Carrying the accepted value forward
+   removes one full objective (one `sample()` per site + reduction) per
+   iteration — when the first trial step is accepted (the common,
+   preconditioned case) this *halves* the objective-evaluation cost.
+
+6. **Fuse the trial step into the objective reduce** — each Armijo trial does
+   a `parallel_for` (write `h = x0 + a*d`, line 212) plus a `parallel_reduce`
+   (objective). Evaluating the objective directly at `x0(i) + a*d(i)` inside
+   the reduce, committing `h` only on acceptance, saves one kernel launch and
+   `n` stores per trial and removes the `deep_copy(h, x0)` undo path
+   (line 224) and the `deep_copy(x0, h)` save (line 199).
+
+7. **Fewer reductions / device syncs per iteration** —
+   - `sy`, `ss`, `yy` (lines 245–247) are three separate `parallel_reduce`s
+     over the same two vectors; fuse into one reduce with a 3-component
+     reducer. The `lbfgs_sy` `parallel_for` (line 240) can be folded in too.
+   - The two-loop recursion issues `2m` blocking dot products per iteration,
+     each a device→host sync; for the production `n = 4096` a GPU run is
+     latency- rather than bandwidth-limited here. The `axpy` can be fused
+     with the next loop's dot (`q -= a*y` and `s·q` in one kernel).
+   - History update: `deep_copy(S[row], s_tmp)` / `deep_copy(Y[row], y_tmp)`
+     (lines 256–257) copy two full vectors per iteration; `std::swap`-ing the
+     `View`s into the slot is free. (The FFT plans are bound to `q`/`d`, not
+     `s_tmp`/`y_tmp`, so swapping is safe.)
+
+### Force kernel (`src/qew_model.h`)
+
+8. **Each noise edge is sampled twice in `gradient()`** — site `i` samples the
+   midpoints of edges `(i-1,i)` and `(i,i+1)` (lines 105–106), so every edge
+   is evaluated by both adjacent sites: `2n` bicubic gathers for `n` distinct
+   edges. A two-pass scheme — edge kernel writing `dv_dy` (and the arclength
+   slope term) per edge into a scratch view, then a cheap combine kernel —
+   does `n` gathers + `~3n` doubles of scratch traffic instead of `2n`
+   gathers (16 loads each). Expect close to 2× on the dominant kernel.
+   (`objective()` already samples each edge exactly once; `rk_site_velocity`
+   re-evaluates at varying `h` and cannot share.)
+
+9. **`sample()` computes derivatives the caller discards** —
+   `objective()` uses only `s.v`, yet `DeviceNoise::sample()`
+   (`src/filtered_noise.h:67`) always accumulates `dv_dgx` *and* `dv_dgy` and
+   `bspline_weights()` always fills `dw[4]`. `gradient()` never uses `dv_dx`
+   anywhere in the codebase. Splitting into value-only / value+dy variants
+   (template tag or separate methods) trims ~⅓–½ of the FLOPs in the 4×4
+   inner loop. Gain is modest where memory-bound, but the Armijo line search
+   is objective-only and benefits directly.
+
+10. **Integer modulo in the inner loops** — each `sample()` does 8 `%`
+    operations for the neighbour wrap (`src/filtered_noise.h:79–80`), plus the
+    stencils' `(i±1) % n`; integer division is expensive on GPUs. After the
+    positive wrap, `ix ∈ [0, nx]` (the `gx == nx` rounding edge case is
+    currently absorbed by `ix % nx`), so conditional add/subtract wraps are
+    sufficient — keep an explicit clamp for that edge case. Power-of-two
+    grids could use bit masks.
+
+### Rosso–Krauth (`src/rosso_krauth.h`)
+
+11. **Redundant mean-`h` reduction per sweep** — `rk_block()` runs two
+    reductions per sweep (lines 161–167): the advance sum and a fresh mean of
+    `h`. Since every move is forward, `mean_h = mean_h0 + Σ_sweeps
+    total_adv / n`; tracking the cumulative advance eliminates the second
+    reduction. Sweeps run up to `max_sweeps = 2·10⁵`, so this halves the
+    per-sweep reduction/sync count.
+
+### Langevin dynamics (`src/dynamics.h`)
+
+12. **RNG pool acquisition per site per step** — `pool.get_state()` /
+    `free_state()` inside the `parallel_for` (lines 48–50) acquires a locked
+    state per site per step. With the benchmark's `nx = 262144` on a GPU the
+    default-sized `Random_XorShift64_Pool` serializes on lock contention.
+    Alternatives: size the pool to the launch width, or pre-fill a normal
+    deviate view per step with `Kokkos::fill_random` and fuse the update.
+    Also, the `gradient()` call inside `step()` inherits item 8's 2× edge
+    duplication.
+
+### Minor / one-time
+
+13. **FIRE** (`src/fire.h`) — `nv2`/`nf2` are two separate reductions
+    (lines 103–108) and the v-step (line 79) could fuse with the power
+    reduction (line 85). Low priority: FIRE is the fallback minimizer.
+14. **Noise setup** (`src/filtered_noise.cpp`) — phase generation, mean and
+    variance are serial host loops over `nx·ny`; `build_coefficients()` stages
+    through three complex 2-D buffers + host mirrors. One-time cost
+    (`setup ≈` seconds in the benchmark), only worth touching if many seeds
+    are ever swept in one run.
+15. **`to_host_vec`** (`src/qew_runs.h:28`) — extra element-wise copy from the
+    mirror into a `std::vector`; could `memcpy` from the mirror when layouts
+    match. Negligible at current snapshot rates.
+
+## Suggested order of attack
+
+Items 5–7 (L-BFGS overhead) and 8 (edge-duplication) are independent,
+low-risk, and cover the production `qew_static` runtime; the
+gradient-consistency and golden-file tests pin their correctness. Item 11 is
+a five-line change covering `qew_depinning`. Item 1 should be fixed (or
+asserted away) before anyone runs an odd-`n` line. Items 2 and 4 are
+correctness-of-use issues rather than speed, and are cheap insurance.