Optimize fast_calc water-balance trio (NDWS/NDWL0/NDWL50) with an Rcpp single-pass kernel

[peetmate]

Problem

The fast_calc water-balance trio — NDWS / NDWL0 / NDWL50 — dominates the 04_indices runtime: ~820–835 s/GCM in the last bake, vs 12–42 s for every other index. On cglabs (CPU-only, no GPU) this is the bottleneck for any trio re-bake (e.g. after the hazards#19 AVAIL fix).

Root cause (CPU)

eabyep_calc (the EABYEP daily soil-water balance) is looped over ~30 days × 240 months per GCM. Each day runs ~10 full-raster terra ops (min/max/arith) plus an AVAIL deepcopy → on the order of hundreds of raster materializations + temp I/O per month. The per-op terra overhead, not the arithmetic, dominates.

Proposed optimization — Rcpp single-pass kernel

The water balance is a per-cell daily recursion (state = AVAIL carried day→day), so it's an ideal single-pass kernel:

• Pass daily inputs as cell×day matrices (pr, evap/ETMAX) + soil capacity/sat vectors + initial AVAIL vector.
• Loop cell×day inside C++, compute ERATIO (and AVAIL/logging) per cell-day, return the ERATIO matrix + final AVAIL.
• NDWS = rowSums(eratio < 0.5); NDWL from the logging term.
• Replaces ~300 sequential terra passes/month with one call.

Same pattern as the R/2.1 §3.4 trend speedup (~63×). Expected: large speedup, on cglabs, no GPU needed.

Constraints / lessons (from prior Rcpp work)

• Kernel must be an installed package (not sourceCpp-into-env) for future/parallel workers — repo already has a trendkernel package; add the water-balance kernel there or a sibling package.
• Preserve the hazards#19 fix semantics: deterministic prior-month AVAIL seed; the C++ recursion makes the day→day state explicit (drops the <<- global).
• NA handling must match terra (NA-aware min/max).
• Validate against the current R output on a real tile before any long run (synthetic probe + bitwise/round compare), per the validate-before-long-run discipline.

Relationship to the GPU issue

Complementary, not either/or:

• This (Rcpp): trio fast on cglabs immediately, stays in R, no new hardware.
• GPU/CuPy issue: bigger bias-correction/downscaling play + Python/xclim convergence, but needs SCIO/Kaggle/Alliance-MPL GPUs (cglabs has none).

Deliverable

• Rcpp water-balance kernel (installed package) + refactored fast_calc_NDWS/NDWL0/NDWL50 to call it.
• Benchmark vs current on a real GCM/month; confirm output equivalence.
• Decision input: does Rcpp close the gap enough that GPU is only needed for bias-correction/downscaling?