diff --git a/seismostats/catalogs/catalog.py b/seismostats/catalogs/catalog.py index 0d6e4a1..ce6031e 100644 --- a/seismostats/catalogs/catalog.py +++ b/seismostats/catalogs/catalog.py @@ -328,6 +328,7 @@ def estimate_mc( self.magnitude.max(), delta_m) + # TODO change once we have a global estimate_mc mc_est = mc_ks(self.magnitude, mcs_test, delta_m, @@ -348,7 +349,7 @@ def estimate_b( weights: list | None = None, b_parameter: str = "b_value", return_std: bool = False, - method: str = "tinti", + method: str = "classic", return_n: bool = False, ) -> float | tuple[float, float] | tuple[float, float, float]: """ @@ -394,14 +395,17 @@ def estimate_b( if delta_m is None: delta_m = self.delta_m + # filter magnitudes above mc without changing the original dataframe + df = self[self.magnitude >= mc] + if method == "positive": # dataframe needs 'time' column to be sorted - if 'time' not in self.columns: + if 'time' not in df.columns: raise ValueError('"time" column needs to be set in order to \ use the b-positive method') - mags = self.sort_values("time").magnitude + mags = df.sort_values("time").magnitude else: - mags = self.magnitude + mags = df.magnitude b_estimate = estimate_b(mags, mc, diff --git a/seismostats/catalogs/tests/test_catalog.py b/seismostats/catalogs/tests/test_catalog.py index 63d79c1..53888b5 100644 --- a/seismostats/catalogs/tests/test_catalog.py +++ b/seismostats/catalogs/tests/test_catalog.py @@ -9,6 +9,7 @@ from seismostats.catalogs.catalog import (REQUIRED_COLS_CATALOG, Catalog, ForecastCatalog) from seismostats.utils.binning import bin_to_precision +from seismostats.analysis.estimate_beta import estimate_b RAW_DATA = {'name': ['Object 1', 'Object 2', 'Object 3'], 'magnitude': [10.0, 12.5, 8.2], @@ -84,8 +85,6 @@ def test_forecast_catalog_strip(): [ (np.array([0.235, -0.235, 4.499, 4.5, 6, 0.1, 1.6]), 0.1), - (np.array([0.235, -0.235, 4.499, 4.5, 6, 0.1, 1.6]), - None), (np.array([0.235, -0.235, 4.499, 5.5, 6, 0.1, 1.6]), 0.2), ([0.235, -0.235, 4.499, 5.5, 6, 0.1, 1.6], @@ -98,43 +97,49 @@ def test_catalog_bin(mag_values: np.ndarray, delta_m: float): assert (catalog.bin_magnitudes( delta_m)['magnitude'].tolist() == bin_to_precision(mag_values, delta_m)).all() + assert (catalog.bin_magnitudes()['magnitude'].tolist() + == bin_to_precision(mag_values, 0.1)).all() return_value = catalog.bin_magnitudes(delta_m, inplace=True) assert (catalog['magnitude'].tolist() == bin_to_precision(mag_values, delta_m)).all() assert return_value is None - assert catalog.delta_m == delta_m -@pytest.mark.parametrize( - "mag_values, delta_m", - [ - (np.array([0.235, -0.235, 4.499, 4.5, 6, 0.1, 1.6]), - 0), - ] -) -def test_catalog_bin_none(mag_values: np.ndarray, delta_m: float): - catalog = Catalog({'magnitude': mag_values}) - - with pytest.raises(ValueError): - catalog.bin_magnitudes(delta_m=delta_m) - - def test_catalog_estimate_mc(): + # TODO once global_mc method is implemented catalog = Catalog({'magnitude': [0.235, -0.235, 4.499, 4.5, 6, 0.1, 1.6]}) with pytest.raises(ValueError): catalog.estimate_mc() - -def test_catalog_estimate_b(): - catalog = Catalog({'magnitude': [0.235, -0.235, 4.499, 4.5, 6, 0.1, 1.6]}) +@pytest.mark.parametrize( + "mag_values, delta_m, mc", + [ + (np.array([0.0, 0.235, 0.238, 4.499, 4.5, 6, 0.1, 1.6]), + 0.001, 0.0) + ] +) +def test_catalog_estimate_b(mag_values, delta_m, mc): + catalog = Catalog({'magnitude': mag_values}) with pytest.raises(ValueError): catalog.estimate_b(mc=None, delta_m=None) catalog.estimate_b(mc=1.0, delta_m=None) catalog.estimate_b(mc=None, delta_m=0.1) + catalog.estimate_b(mc=1.0, delta_m=0.1, method='positive') + + b_value = estimate_b(catalog['magnitude'], mc=mc, delta_m=delta_m, method="classic") + return_value = catalog.estimate_b(mc=mc, delta_m=delta_m, method="classic") + assert catalog.b_value == b_value + assert return_value == b_value + + catalog.mc = mc + catalog.delta_m = delta_m + return_value = catalog.estimate_b(method="classic") + assert catalog.b_value == b_value + assert return_value == b_value def test_to_quakeml(): diff --git a/seismostats/utils/binning.py b/seismostats/utils/binning.py index cd6b62a..f92c9de 100644 --- a/seismostats/utils/binning.py +++ b/seismostats/utils/binning.py @@ -57,8 +57,7 @@ def bin_to_precision(x: np.ndarray | list, delta_x: float = 0.1) -> np.ndarray: if isinstance(x, list): x = np.array(x) - if delta_x is None: - delta_x = 0.1 + d = decimal.Decimal(str(delta_x)) decimal_places = abs(d.as_tuple().exponent) return np.round(normal_round_to_int(x / delta_x) * delta_x, decimal_places) diff --git a/seismostats/utils/tests/test_binning.py b/seismostats/utils/tests/test_binning.py index 326db99..b11b933 100644 --- a/seismostats/utils/tests/test_binning.py +++ b/seismostats/utils/tests/test_binning.py @@ -30,9 +30,6 @@ def test_normal_round(x: float, n: int, rounded_value: float): (np.array([0.235, -0.235, 4.499, 4.5, 6, 0.1, 1.6]), 0.1, np.array([0.2, -0.2, 4.5, 4.5, 6, 0.1, 1.6])), - (np.array([0.235, -0.235, 4.499, 4.5, 6, 0.1, 1.6]), - None, - np.array([0.2, -0.2, 4.5, 4.5, 6, 0.1, 1.6])), (np.array([0.235, -0.235, 4.499, 5.5, 6, 0.1, 1.6]), 0.2, np.array([0.2, -0.2, 4.4, 5.6, 6, 0.2, 1.6])), @@ -46,6 +43,10 @@ def test_bin_to_precision(x: np.ndarray, delta_x: float, y = bin_to_precision(x, delta_x) assert (y == rounded_value).all() + y = bin_to_precision(x) + z = bin_to_precision(x, 0.1) + assert (y == z).all() + def test_bin_to_precision_none(): with pytest.raises(ValueError):