merge main and update documentation of binning_test since tolerance t…

…hreshold changed

seismostats/analysis/avalue/base.py

@@ -5,8 +5,7 @@
 from typing_extensions import Self
 
 from seismostats.utils._config import get_option
-
-# from seismostats.utils.binning import binning_test
+from seismostats.utils.binning import binning_test
 
 
 class AValueEstimator(ABC):
@@ -18,6 +17,7 @@ def __init__(self) -> Self:
         self.scaling_factor: float | None = None
         self.m_ref: float | None = None
         self.b_value: float | None = None
+        self.idx: np.ndarray | None = None
 
         self.__a_value: float | None = None
 
@@ -56,8 +56,12 @@ def calculate(self,
         self.m_ref = m_ref
         self.b_value = b_value
 
-        self._filter_magnitudes()
         self._sanity_checks()
+        self._filter_magnitudes()
+
+        if len(self.magnitudes) == 0:
+            self.__a_value = np.nan
+            return self.__a_value
 
         self.__a_value = self._estimate()
         self.__a_value = self._reference_scaling(self.__a_value)
@@ -75,32 +79,30 @@ def _filter_magnitudes(self) -> np.ndarray:
         '''
         Filter out magnitudes below the completeness magnitude.
         '''
-        idx = self.magnitudes >= self.mc - self.delta_m / 2
-        self.magnitudes = self.magnitudes[idx]
+        self.idx = (self.magnitudes >= self.mc - self.delta_m / 2).nonzero()[0]
+        self.magnitudes = self.magnitudes[self.idx]
 
-        return idx
+        if len(self.magnitudes) == 0:
+            if get_option('warnings') is True:
+                warnings.warn('No magnitudes above the completeness magnitude.')
 
     def _sanity_checks(self):
         '''
         Perform sanity checks on the input data.
         '''
-        # TODO: test that the magnitudes are binned correctly
-        # if self.delta_m == 0:
-        #     tolerance = 1e-08
-        # else:
-        #     tolerance = max(self.delta_m / 100, 1e-08)
-        # assert (
-        #     binning_test(self.magnitudes, self.delta_m, tolerance)
-        # )
-        # 'Magnitudes are not binned correctly.'
-
-        # test if lowest magnitude is much larger than mc
-        if get_option('warnings') is True:
-            if np.min(self.magnitudes) - self.mc > self.delta_m / 2:
-                warnings.warn(
-                    'No magnitudes in the lowest magnitude bin are present. '
-                    'Check if mc is chosen correctly.'
-                )
+        # test magnitude binnning
+        if len(self.magnitudes) > 0:
+            if not binning_test(self.magnitudes, self.delta_m,
+                                check_larger_binning=False):
+                raise ValueError('Magnitudes are not binned correctly.')
+
+            # give warnings
+            if get_option('warnings') is True:
+                if np.min(self.magnitudes) - self.mc > self.delta_m / 2:
+                    warnings.warn(
+                        'No magnitudes in the lowest magnitude bin are present.'
+                        'Check if mc is chosen correctly.'
+                    )
 
     def _reference_scaling(self, a: float) -> float:
         '''

seismostats/analysis/avalue/more_positive.py

@@ -73,23 +73,21 @@ def calculate(self,
                                  b_value=b_value,
                                  )
 
-    def _filter_magnitudes(self):
+    def _filter_magnitudes(self) -> np.ndarray:
         '''
         Filter out magnitudes below the completeness magnitude.
         '''
-        idx = super()._filter_magnitudes()
-
-        self.times = self.times[idx]
-
-        return idx
+        super()._filter_magnitudes()
+        self.times = self.times[self.idx]
 
     def _estimate(self) -> float:
         # order the magnitudes and times
         srt = np.argsort(self.times)
         self.magnitudes = self.magnitudes[srt]
         self.times = self.times[srt]
+        self.idx = self.idx[srt]
 
-        # differences
+        # find next larger event (if it exists)
         idx_next_larger = find_next_larger(
             self.magnitudes, self.delta_m, self.dmc)
         time_diffs = self.times[idx_next_larger] - self.times
@@ -108,5 +106,11 @@ def _estimate(self) -> float:
         time_factor = sum(tau / total_time)
         n_more_pos = sum(~idx_no_next) / time_factor
 
+        # make sure that all attributes are consistent
+        idx_next_larger = idx_next_larger[~idx_no_next]
+        self.magnitudes = self.magnitudes[idx_next_larger]
+        self.times = self.times[idx_next_larger]
+        self.idx = self.idx[idx_next_larger]
+
         # estimate a-value
         return np.log10(n_more_pos)

seismostats/analysis/avalue/positive.py

@@ -76,33 +76,32 @@ def _filter_magnitudes(self) -> np.ndarray:
         '''
         Filter out magnitudes below the completeness magnitude.
         '''
-        idx = super()._filter_magnitudes()
-
-        self.times = self.times[idx]
-
-        return idx
+        super()._filter_magnitudes()
+        self.times = self.times[self.idx]
 
     def _estimate(self) -> float:
-
         # order the magnitudes and times
-        idx = np.argsort(self.times)
-        self.magnitudes = self.magnitudes[idx]
-        self.times = self.times[idx]
+        srt = np.argsort(self.times)
+        self.magnitudes = self.magnitudes[srt]
+        self.times = self.times[srt]
+        self.idx = self.idx[srt]
 
         # differences
         mag_diffs = np.diff(self.magnitudes)
         time_diffs = np.diff(self.times)
 
         # only consider events with magnitude difference >= dmc
-        idx = mag_diffs > self.dmc - self.delta_m / 2
-        mag_diffs = mag_diffs[idx]
-        time_diffs = time_diffs[idx]
+        is_larger = mag_diffs >= self.dmc - self.delta_m / 2
+        time_diffs = time_diffs[is_larger]
+        self.magnitudes = self.magnitudes[1:][is_larger]
+        self.idx = self.idx[1:][is_larger]
 
         # estimate the number of events within the time interval
         total_time = self.times[-1] - self.times[0]
+        self.times = self.times[1:][is_larger]
 
         time_factor = sum(time_diffs / total_time)
-        n_pos = sum(idx) / time_factor
+        n_pos = sum(is_larger) / time_factor
 
         # estimate a-value
         return np.log10(n_pos)

seismostats/analysis/avalue/tests/test_avalue.py

@@ -15,7 +15,7 @@ def test_estimate_a():
     mags = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
 
     # arguments are passed correctly
-    a = estimate_a(mags, mc=1, delta_m=10.0)
+    a = estimate_a(mags, mc=1, delta_m=1)
     assert a == 1.0
 
     mags = np.array([0.9, 0.9, 0.9, 0.9, 10.9])
@@ -37,7 +37,7 @@ def test_estimate_a():
 def test_estimate_a_classic():
     mags = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
     estimator = ClassicAValueEstimator()
-    a_estimate = estimator.calculate(mags, mc=1, delta_m=10.0)
+    a_estimate = estimator.calculate(mags, mc=1, delta_m=1)
     assert a_estimate == 1.0
 
     # reference magnitude is given and b-value given
@@ -47,11 +47,10 @@ def test_estimate_a_classic():
 
     # reference magnitude but no b-value
     with pytest.raises(ValueError):
-        estimator.calculate(mags, mc=1, delta_m=0.0,
-                            m_ref=0)
+        estimator.calculate(mags, mc=1, delta_m=0.0001, m_ref=0)
 
     # reference time is given
-    a_estimate = estimator.calculate(mags, mc=1, delta_m=0.0,
+    a_estimate = estimator.calculate(mags, mc=1, delta_m=0.0001,
                                      scaling_factor=10)
     assert a_estimate == 0.0
 
@@ -91,6 +90,28 @@ def test_estimate_a_positive():
     with pytest.warns(UserWarning):
         estimator.calculate(mags, delta_m=1, mc=-1, times=times)
 
+    # test that index works correctly
+    mags = np.array([0, 0.9, -1, 0.2, 0.5])
+    times = np.array([datetime(2000, 1, 1),
+                      datetime(2000, 1, 2),
+                      datetime(2000, 1, 3),
+                      datetime(2000, 1, 4),
+                      datetime(2000, 1, 5)])
+    estimator.calculate(mags, mc=-1, delta_m=0.1, times=times)
+    assert (mags[estimator.idx] == np.array([0.9, 0.2, 0.5])).all()
+    assert (mags[estimator.idx] == estimator.magnitudes).all()
+
+    # test that time array is correctly used
+    mags = np.array([0, 0.9, 0.5, 0.2, -1])
+    times = np.array([datetime(2000, 1, 1),
+                      datetime(2000, 1, 2),
+                      datetime(2000, 1, 5),
+                      datetime(2000, 1, 4),
+                      datetime(2000, 1, 3)])
+    estimator.calculate(mags, mc=-1, delta_m=0.1, times=times)
+    assert (mags[estimator.idx] == np.array([0.9, 0.2, 0.5])).all()
+    assert (mags[estimator.idx] == estimator.magnitudes).all()
+
 
 @pytest.mark.filterwarnings("ignore")
 def test_estimate_a_more_positive():
@@ -121,3 +142,29 @@ def test_estimate_a_more_positive():
     a = estimator.calculate(
         mags, 0, 0.1, times, b_value=1, dmc=0.1)
     assert_almost_equal(10**a, 16.0)
+
+    # test that index works correctly
+    mags = np.array([0, 0.9, -1, 0.2, 0.5, 1])
+    times = np.array([datetime(2000, 1, 1),
+                      datetime(2000, 1, 2),
+                      datetime(2000, 1, 3),
+                      datetime(2000, 1, 4),
+                      datetime(2000, 1, 5),
+                      datetime(2000, 1, 6)])
+    estimator.calculate(mags, mc=-1, delta_m=0.1, times=times, b_value=1)
+    assert (mags[estimator.idx] == np.array([0.9, 1, 0.2, 0.5, 1])).all()
+    assert (estimator.idx == np.array([1, 5, 3, 4, 5])).all()
+    assert (mags[estimator.idx] == estimator.magnitudes).all()
+    assert (estimator.times == times[estimator.idx]).all()
+
+    # test that time array is correctly used
+    mags = np.array([0, 0.9, 0.5, 0.2, -1])
+    times = np.array([datetime(2000, 1, 1),
+                      datetime(2000, 1, 2),
+                      datetime(2000, 1, 5),
+                      datetime(2000, 1, 4),
+                      datetime(2000, 1, 3)])
+    estimator.calculate(mags, mc=-1, delta_m=0.1, times=times, b_value=1)
+    assert (mags[estimator.idx] == np.array([0.9, 0.2, 0.5])).all()
+    assert (mags[estimator.idx] == estimator.magnitudes).all()
+    assert (estimator.times == times[estimator.idx]).all()

seismostats/analysis/avalue/tests/test_base.py

@@ -8,24 +8,43 @@
 def test_estimate_a_warnings():
     mags = simulate_magnitudes_binned(n=100, b=1, mc=0, delta_m=0.1)
 
-    # TODO: test that uncorrect binninng leads to error
-    # with pytest.raises(AssertionError):
-    #     estimator = ClassicAValueEstimator()
-    #     estimator.calculate(mags, mc=0, delta_m=0.2)
+    # test that uncorrect binninng leads to error
+    with pytest.raises(ValueError):
+        estimator = ClassicAValueEstimator()
+        estimator.calculate(mags, mc=0, delta_m=0.2)
 
     # test that warning is raised if smallest magnitude is much larger than mc
     with pytest.warns(UserWarning):
         estimator = ClassicAValueEstimator()
         estimator.calculate(mags, mc=-1, delta_m=0.1)
 
+    # Magnitudes contain NaN values
+    mags1 = np.array([np.nan, 1, 2, 3, 4])
+    a1 = estimator.calculate(mags1, mc=1, delta_m=1)
+    mags2 = np.array([1, 2, 3, 4])
+    a2 = estimator.calculate(mags2, mc=1, delta_m=1)
+    assert (a1 == a2)
+
+    # No magnitudes above completeness magnitude
+    mags = np.array([0, 0.9, 0.1, 0.2, 0.5])
+    with pytest.warns(UserWarning):
+        estimator.calculate(mags, mc=1, delta_m=0.1)
+    assert (np.isnan(estimator.a_value))
+
 
 def test_by_reference():
-    mags = simulate_magnitudes_binned(n=100, b=1, mc=0, delta_m=0.1)
     estimator = ClassicAValueEstimator()
+
+    mags = simulate_magnitudes_binned(n=100, b=1, mc=0, delta_m=0.1)
     estimator.calculate(mags, mc=0, delta_m=0.1)
     estimator.magnitudes.sort()
     assert not np.array_equal(mags, estimator.magnitudes)
 
+    # test index is working
+    mags = np.array([0, 0.9, -1, 0.2, 0.5])
+    estimator.calculate(mags, mc=0.1, delta_m=0.1)
+    assert (estimator.magnitudes == mags[estimator.idx]).all()
+
 
 def test_reference_scaling():
     mags = simulate_magnitudes_binned(n=100, b=1, mc=0, delta_m=0.1)