From 49934b353788a10d35bb08ef10f6290849af2fad Mon Sep 17 00:00:00 2001 From: "ucgmsim[bot]" <335857963+ucgmsim[bot]@users.noreply.github.com> Date: Fri, 2 Oct 2026 02:18:41 +1300 Subject: [PATCH 1/2] Validate tp and fix noise-duration unit check in calculate_snr MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Closes #236 ## What was wrong `IM/snr_calculation.py:calculate_snr` had two related bugs: 1. The "noise duration must be at least 1s" check compared `noise_duration` (a sample count, equal to `tp`) directly against `1`, instead of `noise_duration * dt` against `1.0` seconds. A noise window of e.g. 0.5s (`tp=50`, `dt=0.01`) was silently accepted; the check only ever fired when `tp <= 0`. 2. `tp` was never validated against the waveform length `nt`. If `tp >= nt`, the signal slice became empty and the call failed deep inside `ims.fourier_amplitude_spectra` with an unrelated `OverflowError` (`cannot convert float infinity to integer`) instead of a clear error at the point of the bad input. Also fixed the accompanying docstrings: `SNRResult.signal_duration` / `noise_duration` were documented as seconds but are actually raw sample counts (now typed and documented as such, per the issue's suggestion — the SNR calculation only uses their ratio, so the unit choice doesn't affect results), and the `tp` parameter was documented as `float` but is used as an integer index (now documented as `int`). ## Fix - Added an upfront `if not (0 < tp < nt): raise ValueError(...)` check in `calculate_snr`. - Changed the noise-duration check to `noise_duration * dt < 1.0`, matching the documented "less than 1s" rule. - Updated `SNRResult` and `calculate_snr` docstrings accordingly. ## Test Added `tests/test_ims.py::test_snr_noise_duration_in_seconds` (a 1000-sample, `dt=0.01` waveform with `tp=50`, i.e. 0.5s of noise, now raises `ValueError` instead of being silently accepted) and `tests/test_ims.py::test_snr_tp_validated_against_record_length` (parametrized over `tp` of `0`, `-1`, `1000`, `1200` against a 1000-sample waveform, all of which now raise a clear `ValueError` instead of `tp >= nt` crashing with `OverflowError` downstream). Both tests failed before the fix for the reasons described in the issue, and pass after it. ## For the reviewer - `tp == nt - 1` (e.g. `tp=999` for `nt=1000`) is intentionally still accepted: it satisfies `0 < tp < nt` and produces a very short (2-sample) signal window, per the issue's suggested fix of validating `0 < tp < nt` rather than requiring a minimum signal length. Worth double-checking this matches the intended behaviour. - `signal_duration`/`noise_duration` are kept as sample counts (not converted to seconds) to avoid changing the numeric values returned by `calculate_snr` beyond what the issue describes; only the type annotation and docstrings were corrected. _Cost $0.48, 27 turns._ 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude --- IM/snr_calculation.py | 24 ++++++++++++++---------- tests/test_ims.py | 22 ++++++++++++++++++++++ 2 files changed, 36 insertions(+), 10 deletions(-) diff --git a/IM/snr_calculation.py b/IM/snr_calculation.py index 1db331b4..bc460a85 100644 --- a/IM/snr_calculation.py +++ b/IM/snr_calculation.py @@ -31,13 +31,13 @@ class SNRResult(NamedTuple): for each component (000, 090, ver). The index represents frequencies and columns represent the different components.""" - signal_duration: float - """Duration of the signal portion in seconds, calculated as the number of - samples in the signal multiplied by the sampling interval (dt).""" + signal_duration: int + """Duration of the signal portion in samples (i.e. the number of samples + in the signal portion of the waveform).""" - noise_duration: float - """Duration of the noise portion in seconds, calculated as the number of - samples in the noise multiplied by the sampling interval (dt).""" + noise_duration: int + """Duration of the noise portion in samples (i.e. the number of samples + in the noise portion of the waveform).""" def calculate_snr( @@ -56,7 +56,7 @@ def calculate_snr( Waveform data as a NumPy array. dt : float The sampling rate of the waveform - tp : float + tp : int The index of the p-arrival frequencies : np.ndarray, optional The frequency vector to use for the SNR calculation, @@ -73,11 +73,15 @@ def calculate_snr( Raises ------ ValueError - If the noise duration is less than 1s and so SNR can't be computed. + If `tp` is not a valid index into the waveform, or if the noise + duration is less than 1s and so SNR can't be computed. """ + (_, _, nt) = waveform.shape + if not (0 < tp < nt): + raise ValueError(f"tp ({tp}) must be between 0 and nt ({nt}), exclusive") + # This extra time is to ensure that when a taper is applied, the signal part of the waveform # is not affected by the tapering. The tapering is applied to the signal and noise separately. - (_, _, nt) = waveform.shape tp_extra = (nt - tp) / 19 # Round up the tp_extra to the nearest highest integer tp_extra = int(np.ceil(tp_extra)) @@ -91,7 +95,7 @@ def calculate_snr( ) # Ensure the noise is not shorter than 1s, if not then skip the calculation - if noise_duration < 1: + if noise_duration * dt < 1.0: raise ValueError("Noise duration is less than 1s") # Apply Taper diff --git a/tests/test_ims.py b/tests/test_ims.py index c3a11b46..f0289519 100644 --- a/tests/test_ims.py +++ b/tests/test_ims.py @@ -284,6 +284,28 @@ def test_snr_benchmark() -> None: ) +def test_snr_noise_duration_in_seconds() -> None: + """A noise window under 1s (in seconds, not samples) must raise.""" + rng = np.random.default_rng(0) + waveform = rng.normal(size=(3, 1, 1000)) + dt = 0.01 + + # tp = 50 samples * dt = 0.01s is only 0.5s of noise, well under 1s. + with pytest.raises(ValueError): + snr_calculation.calculate_snr(waveform, dt, tp=50) + + +@pytest.mark.parametrize("tp", [0, -1, 1000, 1200]) +def test_snr_tp_validated_against_record_length(tp: int) -> None: + """tp must be validated against the number of samples in the waveform.""" + rng = np.random.default_rng(0) + waveform = rng.normal(size=(3, 1, 1000)) + dt = 0.01 + + with pytest.raises(ValueError): + snr_calculation.calculate_snr(waveform, dt, tp=tp) + + def test_all_ims_benchmark() -> None: """Compare benchmark IM calculation against current implementation.""" # Load the DataFrame From ce1419a06fa80899b0b85d376f43b3e0fe4a4830 Mon Sep 17 00:00:00 2001 From: joelridden Date: Mon, 5 Oct 2026 12:42:58 +1300 Subject: [PATCH 2/2] clean up --- IM/snr_calculation.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/IM/snr_calculation.py b/IM/snr_calculation.py index bc460a85..ff06320c 100644 --- a/IM/snr_calculation.py +++ b/IM/snr_calculation.py @@ -31,13 +31,13 @@ class SNRResult(NamedTuple): for each component (000, 090, ver). The index represents frequencies and columns represent the different components.""" - signal_duration: int - """Duration of the signal portion in samples (i.e. the number of samples - in the signal portion of the waveform).""" + signal_duration: float + """Duration of the signal portion in seconds, calculated as the number of + samples in the signal multiplied by the sampling interval (dt).""" - noise_duration: int - """Duration of the noise portion in samples (i.e. the number of samples - in the noise portion of the waveform).""" + noise_duration: float + """Duration of the noise portion in seconds, calculated as the number of + samples in the noise multiplied by the sampling interval (dt).""" def calculate_snr( @@ -87,15 +87,16 @@ def calculate_snr( tp_extra = int(np.ceil(tp_extra)) # Calculate signal and noise areas - signal_duration, noise_duration = (nt - max(tp - tp_extra, 0), tp) + signal_npts, noise_npts = (nt - max(tp - tp_extra, 0), tp) + signal_duration, noise_duration = (signal_npts * dt, noise_npts * dt) signal_acc, noise_acc = ( - waveform[:, :, nt - signal_duration :], - waveform[:, :, :noise_duration], + waveform[:, :, nt - signal_npts :], + waveform[:, :, :noise_npts], ) # Ensure the noise is not shorter than 1s, if not then skip the calculation - if noise_duration * dt < 1.0: + if noise_duration < 1.0: raise ValueError("Noise duration is less than 1s") # Apply Taper