include testing function for estimate_a

seismostats/__init__.py

@@ -4,6 +4,7 @@
 
 # analysis
 from seismostats.analysis.estimate_beta import estimate_b, shi_bolt_confidence
+from seismostats.analysis.estimate_a import estimate_a
 
 # plots
 from seismostats.plots.basics import (
@@ -23,4 +24,4 @@
 # utils
 from seismostats.utils.binning import bin_to_precision
 from seismostats.utils.filtering import cat_intersect_polygon
-from seismostats.utils.simulate_distributions import simulate_magnitudes
+from seismostats.utils.simulate_distributions import simulate_magnitudes, simulate_magnitudes_binned

seismostats/analysis/estimate_a.py

@@ -3,6 +3,7 @@
 
 import numpy as np
 import warnings
+from seismostats.utils._config import get_option
 
 
 def estimate_a(magnitudes: np.ndarray,
@@ -16,7 +17,7 @@ def estimate_a(magnitudes: np.ndarray,
     N = 10 ** (a - b * (m_ref - mc)) (1)
 
     where N is the number of events with magnitude greater than m_ref, which
-    occurred in the time interval T. T should be given as a float- to be
+    occurred in the time interval T. T should be given as a float- to begit 
     precise, it should be the time interval scaled to the time-unit of interest.
     E.g., if the number of events per year are of interest, T should be the
     number of years in which the events occurred.
@@ -41,9 +42,10 @@ def estimate_a(magnitudes: np.ndarray,
     if mc is None:
         mc = magnitudes.min()
     elif magnitudes.min() < mc:
-        warnings.warn(
-            "Completeness magnitude is higher than the lowest magnitude."
-            "Cutting the magnitudes to the completeness magnitude.")
+        if get_option("warnings") is True:
+            warnings.warn(
+                "Completeness magnitude is higher than the lowest magnitude."
+                "Cutting the magnitudes to the completeness magnitude.")
         magnitudes = magnitudes[magnitudes >= mc]
 
     a = np.log10(len(magnitudes))
@@ -52,7 +54,7 @@ def estimate_a(magnitudes: np.ndarray,
     if m_ref is not None:
         if b_value is None:
             raise ValueError(
-                "b_value must be provided if m0 is given")
+                "b_value must be provided if m_ref is given")
         a = a - b_value * (m_ref - mc)
 
     # scale to reference time-interval

seismostats/analysis/tests/test_estimate_a.py

@@ -0,0 +1,28 @@
+import numpy as np
+import warnings
+from seismostats.analysis.estimate_a import estimate_a
+
+
+def test_estimate_a():
+    mags = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
+    a = estimate_a(mags)
+    assert a == 1.0
+
+    # reference magnitude is given and b-value given
+    a = estimate_a(mags, mc=1, m_ref=0, b_value=1)
+    assert a == 2.0
+
+    # reference magnitude but no b-value
+    try:
+        a = estimate_a(mags, mc=1, m_ref=0)
+    except ValueError as e:
+        assert str(e) == "b_value must be provided if m_ref is given"
+
+    # reference time is given
+    a = estimate_a(mags, T=10)
+    assert a == 0.0
+
+    # magnitudes not cut at mc
+    with warnings.catch_warnings(record=True) as w:
+        estimate_a(mags, mc=2)
+        assert w[-1].category == UserWarning

seismostats/analysis/tests/test_estimate_beta.py

@@ -67,7 +67,7 @@ def test_estimate_b(
         )
     assert str(excinfo.value) == "magnitudes are not binned correctly"
 
-    # test that magnitudes smaller than mc lkead to error
+    # test that magnitudes smaller than mc lead to error
     with pytest.raises(AssertionError) as excinfo:
         estimate_b(
             mags,