Performance fixes from the 2026-08 audit (float32 windows, batched irfft, rfft whitening) — timestamps stay float64

[mdenolle]

Summary

Performance fixes from the 2026-08 profiling audit of the S1 cross-correlation pipeline, in one branch as requested. The headline change reverts the waveform window matrix to float32 while keeping timestamps float64 — please read the context below.

Context: issue #352 / PR #353 (please review, @niyiyu; cc @koepflma (Manuela Koepfli))

Issue #352 (reported by Manuela Koepfli) showed irregular CCF time attributes: 256/384 s increments instead of the configured 300 s step. The root cause was that the timestamp array dataS_t was float32 — a float32 unix epoch at current dates only resolves ~256 s. PR #353 (Yiyu Ni) correctly fixed this by moving window cutting to float64.

That fix is preserved and is now locked in by a regression test: dataS_t stays float64, and tests/test_noise_module.py::test_cut_trace_make_stat_timestamps_are_regular builds timestamps at a 2025 epoch with step=300 and asserts np.diff(dataS_t) is exactly 300 — so the #352 fix can never silently regress. The same commit also found and fixed one remaining float32 timestamp array in correlate()'s substack_len != cc_len branch (t_corr was still np.float32).

However, PR #353 also flipped the waveform window matrix dataS to float64, which was not needed for #352 and doubles the memory and FFT cost of everything downstream. This PR takes dataS back to float32. A numerical check at campaign geometry (fs=40 Hz, cc_len=1800 s, 188 windows, RMA whitening, ±32 s lags) shows the float32 waveform path recovers CCFs with correlation r=0.9999996 against float64, residuals at the 1e-4-of-peak level, and identical dv/v from a stretching grid search.

Changes

1. dataS back to float32 (cut_trace_make_stat); dataS_t stays float64, with the regression test above.
2. Batched hermitian irfft in correlate(): the per-window loop of nwin individual scipy.fftpack.ifft(Nfft) complex calls is replaced with one scipy.fft.irfft(axis=1, workers=-1) over the half spectra (the manually mirrored spectrum was hermitian with a zero Nyquist bin, which is exactly irfft's contract). The frequency-domain mean removal now accumulates in float64. This also removes the GIL bottleneck that limited correlate.py's ThreadPoolExecutor scaling (audit: 0.64 s → 0.04 s per pair-day; ~1.7x → near-linear on 4 threads).
3. Dead allocations removed: the unused arr_out in whiten() (433 MB transient in the 2D branch at campaign geometry — allocated, filled, and discarded) and the dead np.zeros(nwin * Nfft2) in correlate() (54 MB, immediately overwritten).
4. whiten_2D uses rfft instead of a full complex fftn: the legacy code zeroed everything outside the whitened positive-frequency band anyway, so only that band needs to be computed and filled. Output shape and values are unchanged. moving_ave_2D now uses scipy.ndimage.uniform_filter1d, verified numerically identical (max diff ~1e-14) to the padded convolve2d for both odd and even window lengths.
5. detrend vectorized: the per-row loop becomes two matrix products (same QR-based operator, in-place update so dtype is preserved; bitwise-identical in the benchmark). The redundant mad() != 0 gate in cut_trace_make_stat is dropped — the std check covers it and mad() is a full O(n log n) median over the day-long trace.
6. Autocorrelation detection: correlate() now takes an is_autocorr flag derived from iiS == iiR at the call site instead of comparing the two full amplitude spectra (~25 ms per pair). The old comparison is kept as a fallback when the flag isn't passed (external callers), but note it was only ever reliable for the XCORR method — for DECONV/COHERENCY the smoothed source spectrum never compares equal, so autocorrelations were misclassified as cross-correlations there; the explicit flag is semantically correct for all methods.

A/B numbers (full-pipeline, seeded synthetic day, RMA/RMA, XCORR)

From tests/test_ccf_regression.py, which runs the complete legacy float64 path (per-window fftpack.ifft, fftn whitening, convolve2d smoothing — verbatim copies of the pre-change code) against the new path:

Metric	Value
per-window CCF correlation (22 windows)	min 0.999999881
stacked CCF correlation	1.000000000 (9 decimals)
max stacked residual / peak	2.5e-8
stretching grid search (±5%, legacy vs new stack)	dv/v = -0.002% (grid resolution 0.004%)
dv/v of a 0.5%-stretched target vs either stack	identical to machine precision (diff = 0.0)
timestamps	identical between paths, exactly regular (Δt = step)

Micro-benchmarks at campaign geometry (188×72000, Apple Silicon):

• correlate() inverse FFTs: 0.39 s → 0.11 s
• whiten_2D: 0.74 s → 0.10 s
• moving_ave_2D (188×5000, N=10): 31 ms → 3 ms
• detrend: 161 ms → 94 ms (identical output)

Tests

• New: tests/test_ccf_regression.py (A/B full-pipeline guard, substack and no-substack branches) and tests/test_noise_module.py::test_cut_trace_make_stat_timestamps_are_regular (Irregular time steps in CCF #352 guard).
• All existing tests pass locally (pytest tests/. integration_tests/.).

@niyiyu and @koepflma — since the precision change touches the ground your #353 fix covered, requesting your review in particular: timestamps remain float64 everywhere (now including the substack branch), and the float32 change is confined to the waveform windows with the A/B numbers above.

🤖 Generated with Claude Code

[Copilot]

Pull request overview

This PR applies performance-focused updates to NoisePy’s S1 cross-correlation pipeline while explicitly preserving the float64 timestamp fix from issue #352, adding regression coverage to prevent future dtype regressions.

Changes:

• Restore waveform window matrix (dataS) to float32 while keeping window timestamps (dataS_t) float64, with a new regression test for regular timestamp stepping.
• Speed up correlation/whitening by batching hermitian inverse FFT via scipy.fft.irfft, switching whiten_2D to rfft, and replacing convolve2d smoothing with uniform_filter1d.
• Add full A/B pipeline regression tests comparing legacy vs new implementations (including both substack and no-substack paths).

Reviewed changes

Copilot reviewed 4 out of 4 changed files in this pull request and generated 2 comments.

File	Description
tests/test_noise_module.py	Adds regression test ensuring window timestamps remain float64 and regular while waveform windows are float32.
tests/test_ccf_regression.py	Adds full-pipeline A/B regression coverage against legacy implementations to guard numerical equivalence.
src/noisepy/seis/noise_module.py	Implements float32 windowing, batched irfft path in correlate(), rfft-based whitening, uniform_filter1d smoothing, and vectorized detrend.
src/noisepy/seis/correlate.py	Threads an is_autocorr flag from call site into noise_module.correlate() to avoid expensive spectrum comparisons.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[codecov]

Codecov Report

❌ Patch coverage is 87.80488% with 5 lines in your changes missing coverage. Please review.
✅ Project coverage is 84.67%. Comparing base (8914125) to head (1b31aaa).

Files with missing lines	Patch %	Lines
src/noisepy/seis/noise_module.py	87.17%	3 Missing and 2 partials ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #356      +/-   ##
==========================================
- Coverage   84.84%   84.67%   -0.18%     
==========================================
  Files           9        9              
  Lines        1986     1970      -16     
  Branches      298      296       -2     
==========================================
- Hits         1685     1668      -17     
- Misses        191      192       +1     
  Partials      110      110              

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[niyiyu]

Agreed that dataS is not necessary to be np.float64. Merging now.