feat: new class structure for a values

seismostats/__init__.py

@@ -2,9 +2,9 @@
 
 # flake8: noqa
 
-from seismostats.analysis.bvalue import estimate_b
 # analysis
-from seismostats.analysis.estimate_a import estimate_a
+from seismostats.analysis.avalue import estimate_a
+from seismostats.analysis.bvalue import estimate_b
 # seismicity
 from seismostats.catalogs.catalog import Catalog, ForecastCatalog
 from seismostats.catalogs.rategrid import ForecastGRRateGrid, GRRateGrid

seismostats/analysis/__init__.py

@@ -1,4 +1,2 @@
 # flake8: noqa
-from seismostats.analysis.estimate_a import (estimate_a_classic,
-                                             estimate_a_positive)
 from seismostats.analysis.estimate_mc import mc_ks, mc_max_curvature

seismostats/analysis/avalue/__init__.py

@@ -0,0 +1,29 @@
+# flake8: noqa
+
+import numpy as np
+
+from seismostats.analysis.avalue.base import AValueEstimator
+from seismostats.analysis.avalue.classic import ClassicAValueEstimator
+from seismostats.analysis.avalue.more_positive import \
+    AMorePositiveAValueEstimator
+from seismostats.analysis.avalue.positive import APositiveAValueEstimator
+
+
+def estimate_a(
+    magnitudes: np.ndarray,
+    mc: float,
+    delta_m: float,
+    scaling_factor: float | None = None,
+    m_ref: float | None = None,
+    b_value: float | None = None,
+    method: AValueEstimator = ClassicAValueEstimator,
+    *args,
+    **kwargs
+) -> float:
+
+    estimator = method()
+    estimator.calculate(magnitudes, mc=mc, delta_m=delta_m,
+                        scaling_factor=scaling_factor, m_ref=m_ref,
+                        b_value=b_value, *args, **kwargs)
+
+    return estimator.a_value

seismostats/analysis/avalue/base.py

@@ -0,0 +1,132 @@
+import warnings
+from abc import ABC, abstractmethod
+
+import numpy as np
+from typing_extensions import Self
+
+from seismostats.utils._config import get_option
+
+# from seismostats.utils.binning import binning_test
+
+
+class AValueEstimator(ABC):
+
+    def __init__(self) -> Self:
+        self.magnitudes: np.ndarray | None = None
+        self.mc: float | None = None
+        self.delta_m: float | None = None
+        self.scaling_factor: float | None = None
+        self.m_ref: float | None = None
+        self.b_value: float | None = None
+
+        self.__a_value: float | None = None
+
+    def calculate(self,
+                  magnitudes: np.ndarray,
+                  mc: float,
+                  delta_m: float,
+                  *args,
+                  scaling_factor: float | None = None,
+                  m_ref: float | None = None,
+                  b_value: float | None = None,
+                  **kwargs) -> float:
+        '''
+        Return the a-value estimate.
+
+        Args:
+            magnitudes:     Array of magnitudes
+            mc:             Completeness magnitude
+            delta_m:        Discretization of magnitudes.
+            scaling_factor: Scaling factor
+                            If given, this is used to normalize the number of
+                            observed events. For example: Volume or area of the
+                            region considered or length of the time interval,
+                            given in the unit of interest.
+            m_ref:          Reference magnitude for which the a-value
+                            is estimated.
+            b_value:        B-value of the Gutenberg-Richter law. Only relevant
+                            when m_ref is not None.
+
+        Returns:
+            a: a-value of the Gutenberg-Richter law.
+        '''
+
+        self.magnitudes = np.array(magnitudes)
+        self.mc = mc
+        self.delta_m = delta_m
+        self.scaling_factor = scaling_factor
+        self.m_ref = m_ref
+        self.b_value = b_value
+
+        self._sanity_checks()
+        self._filtering()
+
+        self.__a_value = self._estimate(*args, **kwargs)
+        self.__a_value = self._reference_scaling(self.__a_value)
+
+        return self.__a_value
+
+    @abstractmethod
+    def _estimate(self, *args, **kwargs) -> float:
+        '''
+        Specific implementation of the a-value estimator.
+        '''
+        pass
+
+    def _filtering(self) -> np.ndarray:
+        '''
+        Filter out magnitudes below the completeness magnitude.
+        '''
+        self.idx = self.magnitudes >= self.mc - self.delta_m / 2
+        self.magnitudes = self.magnitudes[self.idx]
+
+    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.'
+                )
+
+    def _reference_scaling(self, a: float) -> float:
+        '''
+        Args:
+            a: a-value
+
+        Returns:
+            a: scaled a-value
+        '''
+        # scale to reference magnitude
+        if self.m_ref is not None:
+            if self.b_value is None:
+                raise ValueError(
+                    "b_value must be provided if m_ref is given")
+            a = a - self.b_value * (self.m_ref - self.mc)
+
+        # scale to reference time-interal or volume of interest
+        if self.scaling_factor is not None:
+            a = a - np.log10(self.scaling_factor)
+        return a
+
+    @property
+    def a_value(self) -> float:
+        '''
+        Returns the a value of the Gutenberg-Richter law.
+        '''
+        if self.__a_value is None:
+            raise AttributeError('Please calculate the a value first.')
+        return self.__a_value

seismostats/analysis/avalue/classic.py

@@ -0,0 +1,21 @@
+import numpy as np
+
+from seismostats.analysis.avalue.base import AValueEstimator
+
+
+class ClassicAValueEstimator(AValueEstimator):
+    '''
+    Return the a-value of the Gutenberg-Richter (GR) law.
+
+    .. math::
+        N(m) = 10 ^ {a - b \\cdot (m - m_{ref})}
+
+    where N(m) is the number of events with magnitude larger or equal to m
+    that occurred in the timeframe of the catalog.
+    '''
+
+    def __init__(self):
+        super().__init__()
+
+    def _estimate(self) -> float:
+        return np.log10(len(self.magnitudes))

seismostats/analysis/avalue/more_positive.py

@@ -0,0 +1,104 @@
+import warnings
+
+import numpy as np
+
+from seismostats.analysis.avalue.base import AValueEstimator
+from seismostats.analysis.bvalue.utils import find_next_larger
+from seismostats.utils._config import get_option
+
+
+class AMorePositiveAValueEstimator(AValueEstimator):
+    '''
+    Return the a-value of the Gutenberg-Richter (GR) law using only the
+    earthquakes with magnitude m_i >= m_i-1 + dmc.
+
+    Source:
+        van der Elst and Page 2023 (JGR: Solid Earth, Vol 128, Issue 10).
+    '''
+
+    def __init__(self):
+        super().__init__()
+
+    def calculate(self,
+                  magnitudes: np.ndarray,
+                  mc: float,
+                  delta_m: float,
+                  times: np.ndarray,
+                  b_value: float,
+                  scaling_factor: float | None = None,
+                  m_ref: float | None = None,
+                  dmc: float | None = None,
+                  ) -> float:
+        '''
+        Args:
+            magnitudes:     Array of magnitudes
+            mc:             Completeness magnitude
+            delta_m:        Discretization of magnitudes
+            times:          Vector of times of the events, in any format
+                            (datetime, float, etc.).
+            b_value:        B-value of the Gutenberg-Richter law.
+            scaling_factor: Scaling factor
+                            If given, this is used to normalize the number of
+                            observed events. For example: Volume or area of the
+                            region considered or length of the time interval,
+                            given in the unit of interest.
+            m_ref:          Reference magnitude for which the a-value
+                            is estimated.
+            dmc:            Minimum magnitude difference between consecutive
+                            events. If None, the default value is delta_m.
+
+        Returns:
+            a_pos: a-value of the Gutenberg-Richter distribution
+        '''
+
+        if not b_value:
+            # if using estimate_a function, b_value can be None even though its
+            # not passed to the function
+            raise ValueError("b_value must be given")
+
+        return super().calculate(magnitudes,
+                                 mc=mc,
+                                 delta_m=delta_m,
+                                 times=times,
+                                 scaling_factor=scaling_factor,
+                                 m_ref=m_ref,
+                                 b_value=b_value,
+                                 dmc=dmc,
+                                 )
+
+    def _estimate(self, times: np.ndarray, dmc: float | None = None) -> float:
+        if dmc is None:
+            dmc = self.delta_m
+        elif dmc < 0:
+            raise ValueError("dmc must be larger or equal to 0")
+        elif dmc < self.delta_m and get_option("warnings") is True:
+            warnings.warn("dmc is smaller than delta_m, not recommended")
+
+        times = np.array(times)
+        times = times[self.idx]
+
+        # order the magnitudes and times
+        srt = np.argsort(times)
+        self.magnitudes = self.magnitudes[srt]
+        times = times[srt]
+
+        # differences
+        idx_next_larger = find_next_larger(self.magnitudes, self.delta_m, dmc)
+        time_diffs = times[idx_next_larger] - times
+
+        # deal with events which do not have a next larger event
+        idx_no_next = idx_next_larger == 0
+        time_diffs[idx_no_next] = times[-1] - times[idx_no_next]
+
+        # estimate the number of events within the time interval
+        total_time = times[-1] - times[0]
+
+        # scale the time
+        tau = time_diffs * 10**(-self.b_value
+                                * (self.magnitudes + dmc - self.mc))
+
+        time_factor = sum(tau / total_time)
+        n_more_pos = sum(~idx_no_next) / time_factor
+
+        # estimate a-value
+        return np.log10(n_more_pos)