Cross-station CCF equivalence + three-archive demo notebook - #4
Open
mdenolle wants to merge 8 commits into
Open
Conversation
Extends the noisepy equivalence evidence from single-station to the cross-station case: CI.PASC (SCEDC) x BK.PKD (NCEDC) x II.PFO (EarthScope), 2022-01-02, three interstation pairs (~150/230/380 km), through real noisepy compute_fft/correlate at a 20 sps target so the Fourier-resample and sub-sample-alignment branches run inside the chain. All three pairs bit-identical (max abs diff 0.0). Two fidelity bugs found and fixed on the way: - segment_interpolate_np was 1 ulp off on ~30% of samples (2e-5 of CCF peak): numba's type unification computes (1 - nfric) in float64 while nfric*sig1[ii] stays float32; the port now encodes those semantics exactly and a reference-loop test locks them down. - preprocess_raw_np force-cast its output to float32; noisepy leaves float64 after resample when the fric branch does not fire. The adapter now keeps the chain's natural dtype (caught by the new sr=20 chain equivalence case). Also fixes the single-station harness, silently broken since the sample-alignment guard landed: both harnesses now snap the day window to the data's sample grid (CI.PASC starts 0.218 samples off midnight) and feed the same snapped window to both paths. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
A week of BHZ noise from CI.PASC (SCEDC S3), BK.PKD (NCEDC S3) and II.PFO (EarthScope S3), obspy-free from fetch through preprocessing, correlated with NoisePy's own compute_fft/correlate. Shows raw hour-long waveforms per archive, single-day vs 7-day stacked CCFs with 3/1.5 km/s moveouts, and a folded record section vs interstation distance. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
There was a problem hiding this comment.
Pull request overview
Extends the NoisePy numerical-equivalence validation from single-station to cross-station cross-correlation across three S3 archives, and tightens the numpy adapter to reproduce NoisePy/numba mixed-precision behavior and dtype propagation.
Changes:
- Add a new cross-station equivalence harness (
run_xcorr_eval.py) + an integration test wrapper to validate bit-level agreement through real NoisePy FFT/correlation. - Fix two fidelity mismatches in the adapter:
segment_interpolate_npmixed-precision semantics andpreprocess_raw_npoutput dtype preservation (float32 vs float64 depending on branch). - Update the single-station harness to snap day windows to the channel sample grid, and commit new benchmark result JSONs.
Reviewed changes
Copilot reviewed 8 out of 10 changed files in this pull request and generated 4 comments.
Show a summary per file
| File | Description |
|---|---|
| tests/precision/test_xcorr_equivalence.py | Adds integration test invoking the new cross-station harness. |
| tests/precision/test_preprocess_equivalence.py | Adds targeted tests for segment_interpolate_np semantics and exercises the resample branch in chain-equivalence testing. |
| seisfetch/contrib/noisepy_adapter.py | Updates adapter semantics for mixed precision in segment_interpolate_np and preserves natural dtype in preprocess_raw_np. |
| notebooks/README.md | Adds an entry for the new three-archive cross-correlation notebook. |
| benchmarks/results/xcorr_three_archives_2022-01-02.json | Stores cross-station harness results for the three interstation pairs. |
| benchmarks/results/ccf_eval_ci_pasc_2022-01-02.json | Stores updated single-station harness results after window snapping fix. |
| benchmarks/noisepy_eval/run_xcorr_eval.py | New cross-station harness (three archives, three station pairs) using NoisePy compute_fft/correlate. |
| benchmarks/noisepy_eval/run_ccf_eval.py | Adds snap_window() and uses it to align windows across paths for fair comparison. |
| .gitignore | Ignores local cache directories for the evaluation harnesses. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
+227
to
+230
| nf32 = np.float32(nfric) | ||
| sig2[1:-1] = (1.0 - np.float64(nf32)) * sig1[2:].astype(np.float64) + ( | ||
| nf32 * sig1[1:-1] | ||
| ).astype(np.float64) |
Contributor
Author
There was a problem hiding this comment.
Addressed in a258188, with a twist: the suggested plain expression is a trap under numpy 1.x. Value-based casting means a float64 scalar does not upcast a float32 array, so np.float64(...) * float32_array silently runs in float32 there (numpy 2 promotes, so the difference is invisible in this repo's test env — but noisepy pins numpy<2). The now-strict harness caught the 2e-5 CCF divergence immediately. Fixed with np.multiply(..., dtype=np.float64), which forces float64 on every numpy version and drops the intermediate astype copy you flagged. Verified bit-exact vs numba on numpy 1.26.4 and 2.5.1.
Comment on lines
+56
to
+59
| dt_ns = round(1e9 / seg.sampling_rate) | ||
| day0_ns = int(day0.timestamp() * 1e9) | ||
| start_ns = seg.starttime_ns + round((day0_ns - seg.starttime_ns) / dt_ns) * dt_ns | ||
| start = datetime.fromtimestamp(start_ns / 1e9, tz=timezone.utc) |
Contributor
Author
There was a problem hiding this comment.
Fixed in a258188: to_epoch_ns() computes exact integer nanoseconds via timedelta floor-division, used by both snap_window and path_b (which had the same timestamp()*1e9 issue).
Comment on lines
+100
to
+107
| cache.mkdir(parents=True, exist_ok=True) | ||
| f = cache / f"{net}.{sta}.{day}.ms" | ||
| if f.exists(): | ||
| return f.read_bytes() | ||
| kwargs = {} if dc == "earthscope" else {"location": loc, "channel": cha} | ||
| raw = S3OpenClient(datacenter=dc).get_raw(net, sta, day, **kwargs) | ||
| f.write_bytes(raw) | ||
| return raw |
Contributor
Author
There was a problem hiding this comment.
Fixed in a258188: cache keys now encode datacenter and location/channel (scedc.CI.PASC.00.BHZ.<day>.ms, earthscope.II.PFO.all.<day>.ms); the notebook uses the same scheme.
Comment on lines
+140
to
+154
| max_abs = float(np.abs(a - b).max()) | ||
| peak = float(np.abs(a).max()) | ||
| corr = float(np.corrcoef(a, b)[0, 1]) | ||
| pair_ok = corr > 0.99999 | ||
| ok &= pair_ok | ||
| results[f"{s}-{r}"] = { | ||
| "max_abs_diff": max_abs, | ||
| "max_abs_rel_to_peak": max_abs / peak if peak else 0.0, | ||
| "waveform_corr": corr, | ||
| "pass": bool(pair_ok), | ||
| } | ||
| print( | ||
| f"{s}-{r}: corr={corr:.9f} max|diff|/peak={max_abs / peak:.2e} " | ||
| f"-> {'PASS' if pair_ok else 'FAIL'}" | ||
| ) |
Contributor
Author
There was a problem hiding this comment.
Fixed in a258188: all-zero CCFs now raise explicitly, and the pass criterion is max_abs_diff == 0.0 (bit-identity), matching the precision-bank policy that any future drift must be seen and justified. Making it strict paid off immediately — it caught a numpy-1.x casting regression while addressing the comment above.
- snap_window/path_b: exact integer epoch-ns via timedelta division; dt.timestamp()*1e9 exceeds float64's 2**53 integer range and shifted the rounding target by hundreds of ns. - run_xcorr_eval: cache keys now encode datacenter + location/channel so a rerun with different parameters can never silently reuse the wrong bytes (notebook 06 uses the same scheme, re-executed). - Harness pass criterion is now bit-identity (max_abs_diff == 0.0), with an explicit error on all-zero CCFs instead of a ZeroDivisionError. - segment_interpolate_np: dropped the redundant temp copies, but NOT the way Copilot suggested — under numpy 1.x value-based casting a float64 scalar does not upcast a float32 array, so the plain expression silently ran the first product in float32 (2e-5 CCF divergence, caught by the now-strict harness; invisible under numpy 2). np.multiply with dtype=np.float64 forces float64 on every numpy version without an intermediate astype copy. Verified bit-exact vs numba on numpy 1.26.4 and vs the reference loop on numpy 2.5.1. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
January 2022 (93 station-days, 1.1 GB across the three archives). The comparison panel now shows 1 day / 1 week / 1 month: interstation Rayleigh-wave packets emerge cleanly in the month stack on all three pairs, stepping out with distance in the folded record section. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
January-February 2022 (177 station-days, 2 GB across the three archives). The comparison panel now spans 1 day / 1 week / 1 month / 2 months: the interstation Rayleigh-wave packets sharpen and the coda stabilizes between the one- and two-month stacks on all three pairs. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
CI.RPV (Rancho Palos Verdes, SCEDC S3, blank location code) added as a fourth station: the 52-km path across the LA basin produces a dominant symmetric Rayleigh-wave packet visible from a single day of data, and the four-pair record section now shows moveout from 52 to 455 km. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Notebook 06 now demonstrates the full obspy-free path Marine asked for: StationXML fetched from each archive's FDSN station service with seisfetch's raw-HTTP client, parsed and removed with seisfetch.contrib.response (remove_response_np, cross-checked at 7.6e-16 of peak vs obspy on this data), then every station-day converted to ground velocity before the NoisePy correlation chain. A new spectra figure shows all four instruments collapsing onto the same microseism spectrum in physical units. The prose now states explicitly that the bit-identity harness runs at rm_resp=NO and where the response port's own equivalence evidence lives. Fixes a real dirty-metadata bug found on the way: EarthScope II StationXML carries 4-digit fractional seconds, which fromisoformat rejects before Python 3.11 — _parse_iso_utc now normalizes fractions (regression test added). benchmarks/RESULTS.md gains a 'NoisePy equivalence' section rendering the bit-identity JSONs as tables and embedding the two notebook validation figures (HTML inliner learned PNGs). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
What
Extends the noisepy equivalence evidence from single-station to cross-station across three archives, and adds a demo notebook.
benchmarks/noisepy_eval/run_xcorr_eval.py: CI.PASC (SCEDC S3) × BK.PKD (NCEDC S3) × II.PFO (EarthScope S3), 2022-01-02, three interstation pairs (172/291/455 km) through real noisepycompute_fft/correlateat a 20 sps target — so the Fourier-resample and sub-sample-alignment branches run inside the chain, which the single-station eval never exercised. Result: bit-identical (max abs diff 0.0) on all three pairs (benchmarks/results/xcorr_three_archives_2022-01-02.json). Integration test wrapper added.Two fidelity bugs found and fixed
segment_interpolate_npwas 1 ulp off on ~30% of samples (≈2e-5 of CCF peak after stacking). Root cause: numba's type unification in noisepy's original computes(1 - nfric)in float64 whilenfric*sig1[ii]stays float32. The port now encodes those exact mixed-precision semantics; a reference-loop test locks them down, plus a direct test vs actual numba when noisepy is installed.preprocess_raw_npforce-cast its output to float32; noisepy leaves float64 after resample whenever the fric branch doesn't fire. The adapter now keeps the chain's natural dtype (caught by the newsr=20chain-equivalence case).Also fixes the single-station harness, silently broken since the sample-alignment guard landed: archive day files start sub-sample off midnight (PASC: 0.218 samples), so both harnesses now snap the day window to the data's sample grid and feed the same snapped window to both paths. Re-run: still bit-identical, same dv/v cells.
Notebook
notebooks/06_cross_correlation_three_archives.ipynb(executed, committed): a week of BHZ noise, one station per archive, obspy-free from fetch through preprocessing; raw waveforms per archive, single-day vs 7-day stacked CCFs with 3/1.5 km/s moveouts, folded record section vs distance.Tests
pytest -m "not integration"), lint clean🤖 Generated with Claude Code