Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Energy ratio adjustments #3243

Merged
merged 2 commits into from Dec 20, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
66 changes: 46 additions & 20 deletions obspy/signal/tests/test_trigger.py
Expand Up @@ -4,11 +4,13 @@
"""
import gzip
import os
import re
import unittest
import warnings
from ctypes import ArgumentError

import numpy as np
import pytest
from numpy.testing import assert_array_almost_equal, assert_array_equal

from obspy import Stream, UTCDateTime, read
Expand Down Expand Up @@ -598,8 +600,8 @@ def test_classic_sta_lta_c_python(self):
class EnergyRatioTestCase(unittest.TestCase):

def test_all_zero(self):
a = np.zeros(100)
for nsta in range(len(a)):
a = np.zeros(10)
for nsta in range(1, len(a) // 2):
with self.subTest(nsta=nsta):
er = energy_ratio(a, nsta=nsta)
assert_array_equal(er, 0)
Expand All @@ -611,15 +613,8 @@ def test_arange(self):
er_expected = [0., 0., 0., 10., 5.5, 3.793103, 2.98, 2.519481, 0., 0.]
assert_array_almost_equal(er, er_expected)

def test_large_nsta(self):
a = np.arange(100)
for nsta in range(len(a) // 2 + 1, len(a)):
with self.subTest(nsta=nsta):
er = energy_ratio(a, nsta=nsta)
assert_array_equal(er, 0)

def test_all_ones(self):
a = np.ones(100, dtype=np.float32)
a = np.ones(10, dtype=np.float32)
# Forward and backward entries are symmetric -> expecting output '1'
# Fill nsta on both sides with zero to return same length
for nsta in range(1, len(a) // 2 + 1):
Expand All @@ -629,12 +624,31 @@ def test_all_ones(self):
er_exp[nsta: len(a) - nsta + 1] = 1
assert_array_equal(er, er_exp)

def test_nsta_too_large(self):
a = np.empty(10)
nsta = 6
for nsta in (6, 10, 20):
expected_msg = re.escape(
f'nsta ({nsta}) must not be larger than half the length of '
f'the data (10 samples).')
with pytest.raises(ValueError, match=expected_msg):
energy_ratio(a, nsta)

def test_nsta_zero_or_less(self):
a = np.empty(10)
nsta = 6
for nsta in (0, -1, -10):
expected_msg = re.escape(
f'nsta ({nsta}) must not be equal to or less than zero.')
with pytest.raises(ValueError, match=expected_msg):
energy_ratio(a, nsta)


class ModifiedEnergyRatioTestCase(unittest.TestCase):

def test_all_zero(self):
a = np.zeros(100)
for nsta in range(len(a)):
a = np.zeros(10)
for nsta in range(1, len(a) // 2):
with self.subTest(nsta=nsta):
er = modified_energy_ratio(a, nsta=nsta)
assert_array_equal(er, 0)
Expand All @@ -647,15 +661,8 @@ def test_arange(self):
5485.637866, 0., 0.]
assert_array_almost_equal(er, er_expected)

def test_large_nsta(self):
a = np.arange(100)
for nsta in range(len(a) // 2 + 1, len(a)):
with self.subTest(nsta=nsta):
er = modified_energy_ratio(a, nsta=nsta)
assert_array_equal(er, 0)

def test_all_ones(self):
a = np.ones(100, dtype=np.float32)
a = np.ones(10, dtype=np.float32)
# Forward and backward entries are symmetric -> expecting output '1'
# Fill nsta on both sides with zero to return same length
for nsta in range(1, len(a) // 2 + 1):
Expand All @@ -665,6 +672,25 @@ def test_all_ones(self):
er_exp[nsta: len(a) - nsta + 1] = 1
assert_array_equal(er, er_exp)

def test_nsta_too_large(self):
a = np.empty(10)
nsta = 6
for nsta in (6, 10, 20):
expected_msg = re.escape(
f'nsta ({nsta}) must not be larger than half the length of '
f'the data (10 samples).')
with pytest.raises(ValueError, match=expected_msg):
energy_ratio(a, nsta)

def test_nsta_zero_or_less(self):
a = np.empty(10)
nsta = 6
for nsta in (0, -1, -10):
expected_msg = re.escape(
f'nsta ({nsta}) must not be equal to or less than zero.')
with pytest.raises(ValueError, match=expected_msg):
energy_ratio(a, nsta)


def suite():
return unittest.makeSuite(TriggerTestCase, 'test')
Expand Down
12 changes: 9 additions & 3 deletions obspy/signal/trigger.py
Expand Up @@ -305,10 +305,16 @@ def energy_ratio(a, nsta):

.. seealso:: [Han2009]_
"""
if nsta > len(a) // 2 or nsta == 0:
if nsta > len(a) // 2:
# Half forward, half backward -> empty medium
# If nsta is zero, the sum is undefined -> return zero
return np.zeros(len(a), dtype=np.float64)
msg = (
f'nsta ({nsta}) must not be larger than half the length of the '
f'data ({len(a)} samples).')
raise ValueError(msg)
if nsta <= 0:
# If nsta is zero, the sum is undefined
msg = f'nsta ({nsta}) must not be equal to or less than zero.'
raise ValueError(msg)
sig_power = np.r_[0, np.cumsum(a ** 2, dtype=np.float64)]
energy_diff = sig_power[nsta:] - sig_power[:len(sig_power) - nsta]
er = np.zeros(len(a), dtype=np.float64)
Expand Down