From 42a7a755d12315236b6432a20f8d59c16eb36ac0 Mon Sep 17 00:00:00 2001
From: Aron Mirwald <aron.mirwald@sed.ethz.ch>
Date: Thu, 6 Jun 2024 13:27:20 +0200
Subject: [PATCH 01/15] include b-more-pos and b-more-incomplete. tests missing

---
 seismostats/analysis/estimate_beta.py | 166 +++++++++++++++++++++++++-
 1 file changed, 165 insertions(+), 1 deletion(-)

diff --git a/seismostats/analysis/estimate_beta.py b/seismostats/analysis/estimate_beta.py
index 86ca52d..c8036a9 100644
--- a/seismostats/analysis/estimate_beta.py
+++ b/seismostats/analysis/estimate_beta.py
@@ -6,6 +6,7 @@
 import numpy as np
 import pandas as pd
 from scipy.optimize import minimize
+import datetime as dt
 
 from seismostats.utils._config import get_option
 
@@ -307,6 +308,161 @@ def estimate_b_positive(
         return out
 
 
+def estimate_b_more_positive(
+    magnitudes: np.ndarray,
+    delta_m: float = 0,
+    b_parameter: str = "b_value",
+    return_std: bool = False,
+    return_n: bool = False,
+) -> float | tuple[float, float] | tuple[float, float, float]:
+    """Return the b-value estimate calculated using the
+    next positive differences (this means that almost every magnitude has a
+    difference, as opposed to the b-positive method which results in half the
+    data).
+
+    Source:
+        E. Lippiello and G. Petrillo. Journal of Geophysical Research: Solid
+        Earth, 129(2):e2023JB027849, 2024.
+
+    Args:
+        magnitudes: array of magnitudes, sorted in time (first
+                    entry is the earliest earthquake).
+                    To achieve the effect
+                    of reduced STAI, the magnitudes must be ordered in time.
+        delta_m:    discretization of magnitudes. default is no discretization.
+        b_parameter:either 'b-value', then the corresponding value  of the
+                Gutenberg-Richter law is returned, otherwise 'beta' from the
+                exponential distribution [p(M) = exp(-beta*(M-mc))].
+        return_std: if True the standard deviation of beta/b-value (see above)
+                is returned.
+        return_n:   if True the number of events used for the estimation is
+                returned.
+
+    Returns:
+        b:      maximum likelihood beta or b-value, depending on value of
+                input variable 'b_parameter'. Note that the difference is just a
+                factor [b_value = beta * log10(e)]
+        std:    Shi and Bolt estimate of the beta/b-value estimate
+        n:      number of events used for the estimation
+    """
+    mag_diffs = np.zeros(len(magnitudes) - 1)
+    for ii in range(len(magnitudes) - 1):
+        for jj in range(ii + 1, len(magnitudes)):
+            mag_diff_loop = magnitudes[jj] - magnitudes[ii]
+            # print(mag_diff_loop, "diff loop")
+            if mag_diff_loop > 0:
+                mag_diffs[ii] = mag_diff_loop
+                # print("take the value")
+                break
+
+    # print(mag_diffs)
+
+    # only take the values where the next earthquake is larger
+    mag_diffs = abs(mag_diffs[mag_diffs > 0])
+
+    out = estimate_b_tinti(
+        mag_diffs,
+        mc=delta_m,
+        delta_m=delta_m,
+        b_parameter=b_parameter,
+        return_std=return_std,
+    )
+
+    if return_n:
+        if type(out) is tuple:
+            return out + tuple([len(mag_diffs)])
+        else:
+            return out, len(mag_diffs)
+    else:
+        return out
+
+
+def make_more_incomplete(
+    magnitudes: np.ndarray,
+    times: dt.datetime,
+    delta_t: np.timedelta64 = np.timedelta64(60, "s"),
+) -> tuple[np.ndarray, np.ndarray]:
+    # filter out events where the previous event is larger and less than
+    # delta_t seconds away
+    incomplete = False
+    while incomplete is False:
+        incomplete = True
+        idx_del = []
+        for ii in range(1, len(magnitudes)):
+            if magnitudes[ii] < magnitudes[ii - 1]:
+                if times[ii] - times[ii - 1] < delta_t:
+                    idx_del.append(ii)
+                    incomplete = False
+
+        magnitudes = np.delete(magnitudes, idx_del)
+        times = np.delete(times, idx_del)
+
+        if len(magnitudes) < 2:
+            break
+
+    return magnitudes, times
+
+
+def estimate_b_more_incomplete(
+    magnitudes: np.ndarray,
+    times: dt.datetime,
+    delta_t: np.timedelta64 = np.timedelta64(60, "s"),
+    delta_m: float = 0,
+    b_parameter: str = "b_value",
+    return_std: bool = False,
+    return_n: bool = False,
+) -> float | tuple[float, float] | tuple[float, float, float]:
+    """Return the b-value estimate calculated using the b-more-incomplete
+    method proposed by Lippiello and Petrillo (2024). This method first filters
+    out events where the previous event is larger and less than delta_t seconds
+    away and then calculates the b-value using b-more-positive.
+
+    Source:
+        E. Lippiello and G. Petrillo. Journal of Geophysical Research: Solid
+        Earth, 129(2):e2023JB027849, 2024.
+
+    Args:
+        magnitudes: array of magnitudes, sorted in time (first
+                entry is the earliest earthquake).
+                To achieve the effect of reduced STAI, the magnitudes must be
+                ordered in time.
+        times:      array of datetime objects of occurrence of each earthquake
+        delta_t:    time window in seconds to filter out events. default is 60
+                seconds.
+        delta_m:    discretization of magnitudes. default is no discretization.
+        b_parameter:either 'b-value', then the corresponding value  of the
+                Gutenberg-Richter law is returned, otherwise 'beta' from the
+                exponential distribution [p(M) = exp(-beta*(M-mc))].
+        return_std: if True the standard deviation of beta/b-value (see above)
+                is returned.
+        return_n:   if True the number of events used for the estimation is
+                returned.
+
+    Returns:
+        b:      maximum likelihood beta or b-value, depending on value of
+                input variable 'b_parameter'. Note that the difference is just a
+                factor [b_value = beta * log10(e)]
+        std:    Shi and Bolt estimate of the beta/b-value estimate
+        n:      number of events used for the estimation
+    """
+
+    # filter out events where the previous event is larger and less than
+    # delta_t seconds away
+
+    magnitudes, times = make_more_incomplete(magnitudes, times, delta_t)
+
+    if len(magnitudes) < 2:
+        raise ValueError("not enough events to estimate b-value")
+
+    return estimate_b_more_positive(
+        magnitudes,
+        delta_m=delta_m,
+        b_parameter=b_parameter,
+        return_std=return_std,
+        return_n=return_n,
+    )
+
+
 def estimate_b_laplace(
     magnitudes: np.ndarray,
     delta_m: float = 0,
@@ -344,7 +500,7 @@ def estimate_b_laplace(
     mag_diffs = abs(mag_diffs)
     mag_diffs = mag_diffs[mag_diffs > 0]
 
-    return estimate_b_tinti(
+    out = estimate_b_tinti(
         mag_diffs,
         mc=delta_m,
         delta_m=delta_m,
@@ -352,6 +508,14 @@ def estimate_b_laplace(
         return_std=return_std,
     )
 
+    if return_n:
+        if type(out) is tuple:
+            return out + tuple([len(mag_diffs)])
+        else:
+            return out, len(mag_diffs)
+    else:
+        return out
+
 
 def estimate_b_weichert(
     magnitudes: np.ndarray,

From d9627ff5fb40b72cd0f04080f39a4ddab309576c Mon Sep 17 00:00:00 2001
From: Aron Mirwald <aron.mirwald@sed.ethz.ch>
Date: Thu, 6 Jun 2024 13:53:01 +0200
Subject: [PATCH 02/15] add docstring and testing functions

---
 seismostats/analysis/estimate_beta.py         | 23 ++++++--
 .../analysis/tests/test_estimate_beta.py      | 55 +++++++++++++++++++
 2 files changed, 74 insertions(+), 4 deletions(-)

diff --git a/seismostats/analysis/estimate_beta.py b/seismostats/analysis/estimate_beta.py
index c8036a9..29e10b1 100644
--- a/seismostats/analysis/estimate_beta.py
+++ b/seismostats/analysis/estimate_beta.py
@@ -382,8 +382,23 @@ def make_more_incomplete(
     times: dt.datetime,
     delta_t: np.timedelta64 = np.timedelta64(60, "s"),
 ) -> tuple[np.ndarray, np.ndarray]:
-    # filter out events where the previous event is larger and less than
-    # delta_t seconds away
+    """Return filtered magnitudes and times. Filter the magnitudes and times in
+    the following way: If an earthquake is smaller than the previous one and
+    less than ``delta_t`` away, the earthquake is removed.
+
+    Args:
+        magnitudes: array of magnitudes, sorted in time (first
+                entry is the earliest earthquake).
+        times:      array of datetime objects of occurrence of each earthquake
+        delta_t:    time window in seconds to filter out events. default is 60
+                seconds.
+
+    Returns:
+        magnitudes: filtered array of magnitudes
+        times:      filtered array of datetime objects
+
+    """
+
     incomplete = False
     while incomplete is False:
         incomplete = True
@@ -414,8 +429,8 @@ def estimate_b_more_incomplete(
 ) -> float | tuple[float, float] | tuple[float, float, float]:
     """Return the b-value estimate calculated using the b-more-incomplete
     method proposed by Lippiello and Petrillo (2024). This method first filters
-    out events where the previous event is larger and less than delta_t seconds
-    away and then calculates the b-value using b-more-positive.
+    out events where the previous event is larger and less than ``delta_t``
+    seconds away and then calculates the b-value using b-more-positive.
 
     Source:
         E. Lippiello and G. Petrillo. Journal of Geophysical Research: Solid
diff --git a/seismostats/analysis/tests/test_estimate_beta.py b/seismostats/analysis/tests/test_estimate_beta.py
index 2d9c3c7..178d090 100644
--- a/seismostats/analysis/tests/test_estimate_beta.py
+++ b/seismostats/analysis/tests/test_estimate_beta.py
@@ -3,6 +3,7 @@
 
 import numpy as np
 import pytest
+import datetime as dt
 
 # import functions to be tested
 from seismostats.analysis.estimate_beta import (
@@ -14,6 +15,8 @@
     estimate_b_utsu,
     estimate_b_weichert,
     shi_bolt_confidence,
+    estimate_b_more_positive,
+    make_more_incomplete,
 )
 from seismostats.utils.simulate_distributions import simulate_magnitudes_binned
 
@@ -194,6 +197,58 @@ def test_estimate_b_laplace(
     assert abs(b - b_estimate) / b <= precision
 
 
+@pytest.mark.parametrize(
+    "n,b,mc,delta_m,b_parameter,precision",
+    [
+        (100000, 1.2 * np.log(10), 3, 0, "beta", 0.01),
+        (100000, 1, 3, 0.1, "b_value", 0.02),
+    ],
+)
+def test_estimate_b_more_positive(
+    n: int,
+    b: float,
+    mc: float,
+    delta_m: float,
+    b_parameter: str,
+    precision: float,
+):
+    mags = simulate_magnitudes_binned(
+        n, b, mc, delta_m, b_parameter=b_parameter
+    )
+    b_estimate = estimate_b_more_positive(
+        mags, delta_m=delta_m, b_parameter=b_parameter
+    )
+    assert abs(b - b_estimate) / b <= precision
+
+
+def test_make_more_incomplete():
+    magnitudes = [1, 2, 20, 3, 4, 9, 3]
+    times = [
+        dt.datetime(2020, 1, 1),
+        dt.datetime(2020, 1, 2),
+        dt.datetime(2020, 1, 3),
+        dt.datetime(2020, 1, 4),
+        dt.datetime(2020, 1, 5),
+        dt.datetime(2020, 1, 6),
+        dt.datetime(2020, 1, 7),
+    ]
+
+    mags_inc, times_inc = make_more_incomplete(
+        magnitudes, times, delta_t=np.timedelta64(49, "h")
+    )
+
+    assert (mags_inc == [1, 2, 20, 9]).all()
+    assert (
+        times_inc
+        == [
+            dt.datetime(2020, 1, 1),
+            dt.datetime(2020, 1, 2),
+            dt.datetime(2020, 1, 3),
+            dt.datetime(2020, 1, 6),
+        ]
+    ).all()
+
+
 def _create_test_catalog_poisson(a_val_true: float, b_val_true: float):
     """
     creates a synthetic catalog with magnitudes

From 764922f0cb832536022851d33d57f930e4d041ad Mon Sep 17 00:00:00 2001
From: Aron Mirwald <aron.mirwald@sed.ethz.ch>
Date: Fri, 7 Jun 2024 10:36:39 +0200
Subject: [PATCH 03/15] added d_mc to b-positive

---
 seismostats/analysis/estimate_beta.py | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/seismostats/analysis/estimate_beta.py b/seismostats/analysis/estimate_beta.py
index 29e10b1..19223af 100644
--- a/seismostats/analysis/estimate_beta.py
+++ b/seismostats/analysis/estimate_beta.py
@@ -254,6 +254,7 @@ def differences(magnitudes: np.ndarray) -> np.ndarray:
 def estimate_b_positive(
     magnitudes: np.ndarray,
     delta_m: float = 0,
+    d_mc: float | None = None,
     b_parameter: str = "b_value",
     return_std: bool = False,
     return_n: bool = False,
@@ -271,6 +272,8 @@ def estimate_b_positive(
                     To achieve the effect
                     of reduced STAI, the magnitudes must be ordered in time.
         delta_m:    discretization of magnitudes. default is no discretization.
+        d_mc:       cutoff value for the differences (diffferences below this
+                value are not considered). If None, the cutoff is set to delta_m
         b_parameter:either 'b-value', then the corresponding value  of the
                 Gutenberg-Richter law is returned, otherwise 'beta' from the
                 exponential distribution [p(M) = exp(-beta*(M-mc))].
@@ -287,9 +290,13 @@ def estimate_b_positive(
         n:      number of events used for the estimation
     """
 
+    if d_mc is None:
+        d_mc = delta_m
+
     mag_diffs = np.diff(magnitudes)
-    # only take the values where the next earthquake is larger
-    mag_diffs = abs(mag_diffs[mag_diffs > 0])
+    # only take the values where the next earthquake is d_mc larger than the
+    # previous one. delta_m is added to avoid numerical errors
+    mag_diffs = abs(mag_diffs[mag_diffs > d_mc - delta_m / 2])
 
     out = estimate_b_tinti(
         mag_diffs,

From ce3b17256582682c9e19df8f13bb8480525604d3 Mon Sep 17 00:00:00 2001
From: Aron Mirwald <aron.mirwald@sed.ethz.ch>
Date: Fri, 7 Jun 2024 10:37:56 +0200
Subject: [PATCH 04/15] more changes to include d_mc correclty

---
 seismostats/analysis/estimate_beta.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/seismostats/analysis/estimate_beta.py b/seismostats/analysis/estimate_beta.py
index 19223af..93161b7 100644
--- a/seismostats/analysis/estimate_beta.py
+++ b/seismostats/analysis/estimate_beta.py
@@ -300,7 +300,7 @@ def estimate_b_positive(
 
     out = estimate_b_tinti(
         mag_diffs,
-        mc=delta_m,
+        mc=d_mc,
         delta_m=delta_m,
         b_parameter=b_parameter,
         return_std=return_std,

From 57cd97bd0fffd57cb67dc173296a14f85e56b650 Mon Sep 17 00:00:00 2001
From: Aron Mirwald <aron.mirwald@sed.ethz.ch>
Date: Fri, 7 Jun 2024 10:43:35 +0200
Subject: [PATCH 05/15] estimate b_more_incomplete is probably not needed,
 makes it too complicated

---
 seismostats/analysis/estimate_beta.py | 64 ++-------------------------
 1 file changed, 4 insertions(+), 60 deletions(-)

diff --git a/seismostats/analysis/estimate_beta.py b/seismostats/analysis/estimate_beta.py
index 93161b7..be7f6cb 100644
--- a/seismostats/analysis/estimate_beta.py
+++ b/seismostats/analysis/estimate_beta.py
@@ -393,6 +393,10 @@ def make_more_incomplete(
     the following way: If an earthquake is smaller than the previous one and
     less than ``delta_t`` away, the earthquake is removed.
 
+    Source:
+        E. Lippiello and G. Petrillo. Journal of Geophysical Research: Solid
+-       Earth, 129(2):e2023JB027849, 2024.
+
     Args:
         magnitudes: array of magnitudes, sorted in time (first
                 entry is the earliest earthquake).
@@ -425,66 +429,6 @@ def make_more_incomplete(
     return magnitudes, times
 
 
-def estimate_b_more_incomplete(
-    magnitudes: np.ndarray,
-    times: dt.datetime,
-    delta_t: np.timedelta64 = np.timedelta64(60, "s"),
-    delta_m: float = 0,
-    b_parameter: str = "b_value",
-    return_std: bool = False,
-    return_n: bool = False,
-) -> float | tuple[float, float] | tuple[float, float, float]:
-    """Return the b-value estimate calculated using the b-more-incomplete
-    method proposed by Lippiello and Petrillo (2024). This method first filters
-    out events where the previous event is larger and less than ``delta_t``
-    seconds away and then calculates the b-value using b-more-positive.
-
-    Source:
-        E. Lippiello and G. Petrillo. Journal of Geophysical Research: Solid
-        Earth, 129(2):e2023JB027849, 2024.
-
-    Args:
-        magnitudes: array of magnitudes, sorted in time (first
-                entry is the earliest earthquake).
-                To achieve the effect of reduced STAI, the magnitudes must be
-                ordered in time.
-        times:      array of datetime objects of occurrence of each earthquake
-        delta_t:    time window in seconds to filter out events. default is 60
-                seconds.
-        delta_m:    discretization of magnitudes. default is no discretization.
-        b_parameter:either 'b-value', then the corresponding value  of the
-                Gutenberg-Richter law is returned, otherwise 'beta' from the
-                exponential distribution [p(M) = exp(-beta*(M-mc))].
-        return_std: if True the standard deviation of beta/b-value (see above)
-                is returned.
-        return_n:   if True the number of events used for the estimation is
-                returned.
-
-    Returns:
-        b:      maximum likelihood beta or b-value, depending on value of
-                input variable 'b_parameter'. Note that the difference is just a
-                factor [b_value = beta * log10(e)]
-        std:    Shi and Bolt estimate of the beta/b-value estimate
-        n:      number of events used for the estimation
-    """
-
-    # filter out events where the previous event is larger and less than
-    # delta_t seconds away
-
-    magnitudes, times = make_more_incomplete(magnitudes, times, delta_t)
-
-    if len(magnitudes) < 2:
-        raise ValueError("not enough events to estimate b-value")
-
-    return estimate_b_more_positive(
-        magnitudes,
-        delta_m=delta_m,
-        b_parameter=b_parameter,
-        return_std=return_std,
-        return_n=return_n,
-    )
-
-
 def estimate_b_laplace(
     magnitudes: np.ndarray,
     delta_m: float = 0,

From 18e45d4f3fa63f05ce987e1e3e1304c1d706365e Mon Sep 17 00:00:00 2001
From: Aron Mirwald <aron.mirwald@sed.ethz.ch>
Date: Fri, 7 Jun 2024 10:51:06 +0200
Subject: [PATCH 06/15] rename d_mc to dmc and include proper testing

---
 seismostats/analysis/estimate_beta.py            | 12 ++++++------
 seismostats/analysis/tests/test_estimate_beta.py |  9 +++++----
 2 files changed, 11 insertions(+), 10 deletions(-)

diff --git a/seismostats/analysis/estimate_beta.py b/seismostats/analysis/estimate_beta.py
index be7f6cb..b0f3b80 100644
--- a/seismostats/analysis/estimate_beta.py
+++ b/seismostats/analysis/estimate_beta.py
@@ -254,7 +254,7 @@ def differences(magnitudes: np.ndarray) -> np.ndarray:
 def estimate_b_positive(
     magnitudes: np.ndarray,
     delta_m: float = 0,
-    d_mc: float | None = None,
+    dmc: float | None = None,
     b_parameter: str = "b_value",
     return_std: bool = False,
     return_n: bool = False,
@@ -272,7 +272,7 @@ def estimate_b_positive(
                     To achieve the effect
                     of reduced STAI, the magnitudes must be ordered in time.
         delta_m:    discretization of magnitudes. default is no discretization.
-        d_mc:       cutoff value for the differences (diffferences below this
+        dmc:       cutoff value for the differences (diffferences below this
                 value are not considered). If None, the cutoff is set to delta_m
         b_parameter:either 'b-value', then the corresponding value  of the
                 Gutenberg-Richter law is returned, otherwise 'beta' from the
@@ -290,17 +290,17 @@ def estimate_b_positive(
         n:      number of events used for the estimation
     """
 
-    if d_mc is None:
-        d_mc = delta_m
+    if dmc is None:
+        dmc = delta_m
 
     mag_diffs = np.diff(magnitudes)
     # only take the values where the next earthquake is d_mc larger than the
     # previous one. delta_m is added to avoid numerical errors
-    mag_diffs = abs(mag_diffs[mag_diffs > d_mc - delta_m / 2])
+    mag_diffs = abs(mag_diffs[mag_diffs > dmc - delta_m / 2])
 
     out = estimate_b_tinti(
         mag_diffs,
-        mc=d_mc,
+        mc=dmc,
         delta_m=delta_m,
         b_parameter=b_parameter,
         return_std=return_std,
diff --git a/seismostats/analysis/tests/test_estimate_beta.py b/seismostats/analysis/tests/test_estimate_beta.py
index 178d090..232dea4 100644
--- a/seismostats/analysis/tests/test_estimate_beta.py
+++ b/seismostats/analysis/tests/test_estimate_beta.py
@@ -150,10 +150,10 @@ def test_differences(magnitudes: np.ndarray, mag_diffs: np.ndarray):
 
 
 @pytest.mark.parametrize(
-    "n,b,mc,delta_m,b_parameter,precision",
+    "n,b,mc,delta_m,b_parameter,dmc,precision",
     [
-        (1000000, 1.2 * np.log(10), 3, 0, "beta", 0.005),
-        (1000000, np.log(10), 3, 0.1, "beta", 0.01),
+        (1000000, 1.2, 3, 0, "b_value", None, 0.005),
+        (1000000, np.log(10), 3, 0.1, "beta", 1, 0.01),
     ],
 )
 def test_estimate_b_positive(
@@ -162,13 +162,14 @@ def test_estimate_b_positive(
     mc: float,
     delta_m: float,
     b_parameter: str,
+    dmc: float,
     precision: float,
 ):
     mags = simulate_magnitudes_binned(
         n, b, mc, delta_m, b_parameter=b_parameter
     )
     b_estimate = estimate_b_positive(
-        mags, delta_m=delta_m, b_parameter=b_parameter
+        mags, delta_m=delta_m, dmc=dmc, b_parameter=b_parameter
     )
     assert abs(b - b_estimate) / b <= precision
 

From ec1f1e4de8e31353037903c1b87a73d7953ed0a9 Mon Sep 17 00:00:00 2001
From: Aron Mirwald <aron.mirwald@sed.ethz.ch>
Date: Fri, 7 Jun 2024 11:38:54 +0200
Subject: [PATCH 07/15] include dmc in b-more-pos

---
 seismostats/analysis/estimate_beta.py            | 9 +++++++--
 seismostats/analysis/tests/test_estimate_beta.py | 9 +++++----
 2 files changed, 12 insertions(+), 6 deletions(-)

diff --git a/seismostats/analysis/estimate_beta.py b/seismostats/analysis/estimate_beta.py
index b0f3b80..f729c04 100644
--- a/seismostats/analysis/estimate_beta.py
+++ b/seismostats/analysis/estimate_beta.py
@@ -318,6 +318,7 @@ def estimate_b_positive(
 def estimate_b_more_positive(
     magnitudes: np.ndarray,
     delta_m: float = 0,
+    dmc: float | None = None,
     b_parameter: str = "b_value",
     return_std: bool = False,
     return_n: bool = False,
@@ -352,12 +353,16 @@ def estimate_b_more_positive(
         std:    Shi and Bolt estimate of the beta/b-value estimate
         n:      number of events used for the estimation
     """
+
+    if dmc is None:
+        dmc = delta_m
+
     mag_diffs = np.zeros(len(magnitudes) - 1)
     for ii in range(len(magnitudes) - 1):
         for jj in range(ii + 1, len(magnitudes)):
             mag_diff_loop = magnitudes[jj] - magnitudes[ii]
             # print(mag_diff_loop, "diff loop")
-            if mag_diff_loop > 0:
+            if mag_diff_loop > dmc - delta_m / 2:
                 mag_diffs[ii] = mag_diff_loop
                 # print("take the value")
                 break
@@ -369,7 +374,7 @@ def estimate_b_more_positive(
 
     out = estimate_b_tinti(
         mag_diffs,
-        mc=delta_m,
+        mc=dmc,
         delta_m=delta_m,
         b_parameter=b_parameter,
         return_std=return_std,
diff --git a/seismostats/analysis/tests/test_estimate_beta.py b/seismostats/analysis/tests/test_estimate_beta.py
index 232dea4..592d085 100644
--- a/seismostats/analysis/tests/test_estimate_beta.py
+++ b/seismostats/analysis/tests/test_estimate_beta.py
@@ -199,10 +199,10 @@ def test_estimate_b_laplace(
 
 
 @pytest.mark.parametrize(
-    "n,b,mc,delta_m,b_parameter,precision",
+    "n,b,mc,delta_m,b_parameter,dmc,precision",
     [
-        (100000, 1.2 * np.log(10), 3, 0, "beta", 0.01),
-        (100000, 1, 3, 0.1, "b_value", 0.02),
+        (100000, 1.2 * np.log(10), 3, 0, "beta", None, 0.01),
+        (100000, 1, 3, 0.1, "b_value", 1, 0.02),
     ],
 )
 def test_estimate_b_more_positive(
@@ -211,13 +211,14 @@ def test_estimate_b_more_positive(
     mc: float,
     delta_m: float,
     b_parameter: str,
+    dmc: float,
     precision: float,
 ):
     mags = simulate_magnitudes_binned(
         n, b, mc, delta_m, b_parameter=b_parameter
     )
     b_estimate = estimate_b_more_positive(
-        mags, delta_m=delta_m, b_parameter=b_parameter
+        mags, delta_m=delta_m, dmc=dmc, b_parameter=b_parameter
     )
     assert abs(b - b_estimate) / b <= precision
 

From 0b485616073d2cc95d9568585da6569ca8b2c19a Mon Sep 17 00:00:00 2001
From: Aron Mirwald <aron.mirwald@sed.ethz.ch>
Date: Sat, 8 Jun 2024 16:33:01 +0200
Subject: [PATCH 08/15] make_mre incomplete: magnitude must be ordered in time

---
 seismostats/analysis/estimate_beta.py            | 5 +++++
 seismostats/analysis/tests/test_estimate_beta.py | 6 +++---
 2 files changed, 8 insertions(+), 3 deletions(-)

diff --git a/seismostats/analysis/estimate_beta.py b/seismostats/analysis/estimate_beta.py
index f729c04..0d5b895 100644
--- a/seismostats/analysis/estimate_beta.py
+++ b/seismostats/analysis/estimate_beta.py
@@ -415,6 +415,11 @@ def make_more_incomplete(
 
     """
 
+    # sort magnitudes in time
+    idx_sort = np.argsort(times)
+    magnitudes = magnitudes[idx_sort]
+    times = times[idx_sort]
+
     incomplete = False
     while incomplete is False:
         incomplete = True
diff --git a/seismostats/analysis/tests/test_estimate_beta.py b/seismostats/analysis/tests/test_estimate_beta.py
index 592d085..e941d94 100644
--- a/seismostats/analysis/tests/test_estimate_beta.py
+++ b/seismostats/analysis/tests/test_estimate_beta.py
@@ -224,8 +224,8 @@ def test_estimate_b_more_positive(
 
 
 def test_make_more_incomplete():
-    magnitudes = [1, 2, 20, 3, 4, 9, 3]
-    times = [
+    magnitudes = np.array([1, 2, 20, 3, 4, 9, 3])
+    times = np.array([
         dt.datetime(2020, 1, 1),
         dt.datetime(2020, 1, 2),
         dt.datetime(2020, 1, 3),
@@ -233,7 +233,7 @@ def test_make_more_incomplete():
         dt.datetime(2020, 1, 5),
         dt.datetime(2020, 1, 6),
         dt.datetime(2020, 1, 7),
-    ]
+    ])
 
     mags_inc, times_inc = make_more_incomplete(
         magnitudes, times, delta_t=np.timedelta64(49, "h")

From 21a4a520d882fcbf6c2a1207277c9cdc4e5f46f8 Mon Sep 17 00:00:00 2001
From: Aron Mirwald <aron.mirwald@sed.ethz.ch>
Date: Sat, 8 Jun 2024 16:51:24 +0200
Subject: [PATCH 09/15] added possibility to return delted index in make more
 incomplete

---
 seismostats/analysis/estimate_beta.py | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/seismostats/analysis/estimate_beta.py b/seismostats/analysis/estimate_beta.py
index 0d5b895..75f2681 100644
--- a/seismostats/analysis/estimate_beta.py
+++ b/seismostats/analysis/estimate_beta.py
@@ -393,6 +393,7 @@ def make_more_incomplete(
     magnitudes: np.ndarray,
     times: dt.datetime,
     delta_t: np.timedelta64 = np.timedelta64(60, "s"),
+    return_idx: bool = False,
 ) -> tuple[np.ndarray, np.ndarray]:
     """Return filtered magnitudes and times. Filter the magnitudes and times in
     the following way: If an earthquake is smaller than the previous one and
@@ -408,10 +409,13 @@ def make_more_incomplete(
         times:      array of datetime objects of occurrence of each earthquake
         delta_t:    time window in seconds to filter out events. default is 60
                 seconds.
+        return_idx: if True the indices of the events that were removed are
+                returned
 
     Returns:
         magnitudes: filtered array of magnitudes
         times:      filtered array of datetime objects
+        idx:        indices of the events that were removed
 
     """
 
@@ -436,6 +440,9 @@ def make_more_incomplete(
         if len(magnitudes) < 2:
             break
 
+    if return_idx is True:
+        return magnitudes, times, idx_del
+
     return magnitudes, times
 
 

From 827ee7f1cfb2fafd9a40e9abb264d59040f1e322 Mon Sep 17 00:00:00 2001
From: Aron Mirwald <aron.mirwald@sed.ethz.ch>
Date: Sat, 8 Jun 2024 17:06:38 +0200
Subject: [PATCH 10/15] corrected error in make_more_inc

---
 seismostats/analysis/estimate_beta.py | 29 +++++++++++++--------------
 1 file changed, 14 insertions(+), 15 deletions(-)

diff --git a/seismostats/analysis/estimate_beta.py b/seismostats/analysis/estimate_beta.py
index 75f2681..cf597ee 100644
--- a/seismostats/analysis/estimate_beta.py
+++ b/seismostats/analysis/estimate_beta.py
@@ -424,21 +424,20 @@ def make_more_incomplete(
     magnitudes = magnitudes[idx_sort]
     times = times[idx_sort]
 
-    incomplete = False
-    while incomplete is False:
-        incomplete = True
-        idx_del = []
-        for ii in range(1, len(magnitudes)):
-            if magnitudes[ii] < magnitudes[ii - 1]:
-                if times[ii] - times[ii - 1] < delta_t:
-                    idx_del.append(ii)
-                    incomplete = False
-
-        magnitudes = np.delete(magnitudes, idx_del)
-        times = np.delete(times, idx_del)
-
-        if len(magnitudes) < 2:
-            break
+    idx_del = []
+    for ii in range(1, len(magnitudes)):
+        # get all times that are closer than delta_t
+        idx_close = np.where(times[ii] - times[:ii] < delta_t)[0]
+
+        # check if these events are larger than the current event
+        idx_del_loop = np.where(magnitudes[idx_close] > magnitudes[ii])[0]
+
+        # if there are any, remove the current event
+        if len(idx_del_loop) > 0:
+            idx_del.append(ii)
+
+    magnitudes = np.delete(magnitudes, idx_del)
+    times = np.delete(times, idx_del)
 
     if return_idx is True:
         return magnitudes, times, idx_del

From a76caaf5e41c3bbda43359fd5f41e9ef1ff555ce Mon Sep 17 00:00:00 2001
From: Nicolas Schmid <nicolas@breiten.ch>
Date: Thu, 20 Jun 2024 16:49:42 +0200
Subject: [PATCH 11/15] fix: use almost_equal for tests so we don't get
 floating point errors

---
 seismostats/analysis/tests/test_estimate_mc.py | 18 ++++++++----------
 1 file changed, 8 insertions(+), 10 deletions(-)

diff --git a/seismostats/analysis/tests/test_estimate_mc.py b/seismostats/analysis/tests/test_estimate_mc.py
index 54faeb7..c94e3ca 100644
--- a/seismostats/analysis/tests/test_estimate_mc.py
+++ b/seismostats/analysis/tests/test_estimate_mc.py
@@ -1,13 +1,11 @@
-from numpy.testing import assert_allclose, assert_equal
-import numpy as np
 import pickle
+
+import numpy as np
 import pytest
+from numpy.testing import assert_allclose, assert_almost_equal, assert_equal
 
-from seismostats.analysis.estimate_mc import (
-    empirical_cdf,
-    mc_max_curvature,
-    mc_ks,
-)
+from seismostats.analysis.estimate_mc import (empirical_cdf, mc_ks,
+                                              mc_max_curvature)
 
 
 @pytest.fixture
@@ -148,7 +146,7 @@ def test_estimate_mc_ks(mags, mcs):
         mags, delta_m=0.1, mcs_test=mcs, p_pass=0.1
     )
     assert_equal(1.1, best_mc)
-    assert_equal(2.242124985031149, best_beta)
+    assert_almost_equal(2.242124985031149, best_beta)
     assert_equal([0.8, 0.9, 1.0, 1.1], mcs_tested)
     assert_allclose(
         [
@@ -190,7 +188,7 @@ def test_estimate_mc_ks(mags, mcs):
         mags, delta_m=0.1, p_pass=0.1
     )
     assert_equal(1.1, best_mc)
-    assert_equal(2.242124985031149, best_beta)
+    assert_almost_equal(2.242124985031149, best_beta)
     assert_equal([1.0, 1.1], mcs_tested)
 
     # test when b-positive is used
@@ -198,7 +196,7 @@ def test_estimate_mc_ks(mags, mcs):
         mags, delta_m=0.1, b_method="positive"
     )
     assert_equal(1.5, best_mc)
-    assert_equal(3.2542240043462796, best_beta)
+    assert_almost_equal(3.2542240043462796, best_beta)
     assert_equal(len(mcs_tested), 6)
     assert_equal(len(betas), 6)
     assert_equal(len(ks_ds), 6)

From 45d45fed06207ca8a43c8a51c9ca770b9b7aff0e Mon Sep 17 00:00:00 2001
From: Aron Mirwald <aron.mirwald@sed.ethz.ch>
Date: Thu, 20 Jun 2024 19:05:31 +0200
Subject: [PATCH 12/15] leilas review, added a warning if dmc is smaller than
 delta_m and error when smaller than 0

---
 seismostats/analysis/estimate_beta.py         | 37 ++++++++++---------
 .../analysis/tests/test_estimate_beta.py      |  8 +++-
 2 files changed, 27 insertions(+), 18 deletions(-)

diff --git a/seismostats/analysis/estimate_beta.py b/seismostats/analysis/estimate_beta.py
index cf597ee..7b063f7 100644
--- a/seismostats/analysis/estimate_beta.py
+++ b/seismostats/analysis/estimate_beta.py
@@ -6,7 +6,6 @@
 import numpy as np
 import pandas as pd
 from scipy.optimize import minimize
-import datetime as dt
 
 from seismostats.utils._config import get_option
 
@@ -292,6 +291,10 @@ def estimate_b_positive(
 
     if dmc is None:
         dmc = delta_m
+    elif dmc < 0:
+        raise ValueError("dmc must be larger or equal to 0")
+    elif dmc < delta_m:
+        warnings.warn("dmc is smaller than delta_m, not recommended")
 
     mag_diffs = np.diff(magnitudes)
     # only take the values where the next earthquake is d_mc larger than the
@@ -356,21 +359,21 @@ def estimate_b_more_positive(
 
     if dmc is None:
         dmc = delta_m
+    elif dmc < 0:
+        raise ValueError("dmc must be larger or equal to 0")
+    elif dmc < delta_m:
+        warnings.warn("dmc is smaller than delta_m, not recommended")
 
-    mag_diffs = np.zeros(len(magnitudes) - 1)
+    mag_diffs = - np.ones(len(magnitudes) - 1) * delta_m
     for ii in range(len(magnitudes) - 1):
         for jj in range(ii + 1, len(magnitudes)):
             mag_diff_loop = magnitudes[jj] - magnitudes[ii]
-            # print(mag_diff_loop, "diff loop")
             if mag_diff_loop > dmc - delta_m / 2:
                 mag_diffs[ii] = mag_diff_loop
-                # print("take the value")
                 break
 
-    # print(mag_diffs)
-
     # only take the values where the next earthquake is larger
-    mag_diffs = abs(mag_diffs[mag_diffs > 0])
+    mag_diffs = abs(mag_diffs[mag_diffs > - delta_m / 2])
 
     out = estimate_b_tinti(
         mag_diffs,
@@ -391,7 +394,7 @@ def estimate_b_more_positive(
 
 def make_more_incomplete(
     magnitudes: np.ndarray,
-    times: dt.datetime,
+    times: np.array,
     delta_t: np.timedelta64 = np.timedelta64(60, "s"),
     return_idx: bool = False,
 ) -> tuple[np.ndarray, np.ndarray]:
@@ -409,13 +412,13 @@ def make_more_incomplete(
         times:      array of datetime objects of occurrence of each earthquake
         delta_t:    time window in seconds to filter out events. default is 60
                 seconds.
-        return_idx: if True the indices of the events that were removed are
+        return_idx: if True the indices of the events that were kept are
                 returned
 
     Returns:
         magnitudes: filtered array of magnitudes
         times:      filtered array of datetime objects
-        idx:        indices of the events that were removed
+        idx:        indices of the events that were kept
 
     """
 
@@ -424,23 +427,23 @@ def make_more_incomplete(
     magnitudes = magnitudes[idx_sort]
     times = times[idx_sort]
 
-    idx_del = []
+    idx = np.full(len(magnitudes), True)
     for ii in range(1, len(magnitudes)):
         # get all times that are closer than delta_t
         idx_close = np.where(times[ii] - times[:ii] < delta_t)[0]
 
         # check if these events are larger than the current event
-        idx_del_loop = np.where(magnitudes[idx_close] > magnitudes[ii])[0]
+        idx_loop = magnitudes[idx_close] > magnitudes[ii]
 
         # if there are any, remove the current event
-        if len(idx_del_loop) > 0:
-            idx_del.append(ii)
+        if sum(idx_loop) > 0:
+            idx[ii] = False
 
-    magnitudes = np.delete(magnitudes, idx_del)
-    times = np.delete(times, idx_del)
+    magnitudes = magnitudes[idx]
+    times = times[idx]
 
     if return_idx is True:
-        return magnitudes, times, idx_del
+        return magnitudes, times, idx
 
     return magnitudes, times
 
diff --git a/seismostats/analysis/tests/test_estimate_beta.py b/seismostats/analysis/tests/test_estimate_beta.py
index e941d94..e35af24 100644
--- a/seismostats/analysis/tests/test_estimate_beta.py
+++ b/seismostats/analysis/tests/test_estimate_beta.py
@@ -202,7 +202,7 @@ def test_estimate_b_laplace(
     "n,b,mc,delta_m,b_parameter,dmc,precision",
     [
         (100000, 1.2 * np.log(10), 3, 0, "beta", None, 0.01),
-        (100000, 1, 3, 0.1, "b_value", 1, 0.02),
+        (100000, 1, 3, 0.1, "b_value", 1, 0.04),
     ],
 )
 def test_estimate_b_more_positive(
@@ -250,6 +250,12 @@ def test_make_more_incomplete():
         ]
     ).all()
 
+    mags_inc, times_inc, idx = make_more_incomplete(
+        magnitudes, times, delta_t=np.timedelta64(49, "h"), return_idx=True
+    )
+
+    assert (mags_inc == magnitudes[idx]).all()
+
 
 def _create_test_catalog_poisson(a_val_true: float, b_val_true: float):
     """

From f1276bedd10455ea769926da99f0c07ea54e2603 Mon Sep 17 00:00:00 2001
From: Aron Mirwald <aron.mirwald@sed.ethz.ch>
Date: Thu, 20 Jun 2024 19:09:40 +0200
Subject: [PATCH 13/15] add new estimation function to analysis init

---
 seismostats/analysis/__init__.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/seismostats/analysis/__init__.py b/seismostats/analysis/__init__.py
index c8d6e2e..63494e2 100644
--- a/seismostats/analysis/__init__.py
+++ b/seismostats/analysis/__init__.py
@@ -5,5 +5,8 @@
                                                 estimate_b_positive,
                                                 estimate_b_tinti,
                                                 estimate_b_utsu,
-                                                estimate_b_weichert)
+                                                estimate_b_weichert,
+                                                estimate_b_kijko_smit,
+                                                estimate_b_more_positive,
+                                                make_more_incomplete)
 from seismostats.analysis.estimate_mc import mc_ks, mc_max_curvature

From bd35a5d3066142565efd898169bdaf5304632bae Mon Sep 17 00:00:00 2001
From: Aron Mirwald <aron.mirwald@sed.ethz.ch>
Date: Fri, 21 Jun 2024 14:52:28 +0200
Subject: [PATCH 14/15] added to docs, few changes to redfuces number of error
 messages when make server

---
 docs/source/reference/analysis.md     | 11 +++++++++++
 docs/source/reference/utils.md        |  2 +-
 docs/source/user/getting_started.md   | 10 +++++-----
 seismostats/analysis/__init__.py      |  5 ++++-
 seismostats/analysis/estimate_beta.py |  5 ++---
 seismostats/analysis/estimate_mc.py   | 13 +++++++------
 6 files changed, 30 insertions(+), 16 deletions(-)

diff --git a/docs/source/reference/analysis.md b/docs/source/reference/analysis.md
index 32a18bb..70eebe8 100644
--- a/docs/source/reference/analysis.md
+++ b/docs/source/reference/analysis.md
@@ -14,6 +14,7 @@
     shi_bolt_confidence
     analysis.estimate_b_tinti
     analysis.estimate_b_positive
+    analysis.estimate_b_more_positive
     analysis.estimate_b_utsu
     analysis.estimate_b_laplace
     analysis.estimate_b_weichert
@@ -27,4 +28,14 @@
 
     analysis.mc_ks
     analysis.mc_max_curvature
+```
+
+## Other
+```{eval-rst}
+.. autosummary::
+    :toctree: api/
+
+    analysis.make_more_incomplete
+    analysis.beta_to_b_value
+    analysis.b_value_to_beta
 ```
\ No newline at end of file
diff --git a/docs/source/reference/utils.md b/docs/source/reference/utils.md
index 430d5ec..2a1b675 100644
--- a/docs/source/reference/utils.md
+++ b/docs/source/reference/utils.md
@@ -23,7 +23,7 @@
     :toctree: api/
 
     utils.simulate_magnitudes
-    utils.simulated_magnitudes_binned
+    utils.simulate_magnitudes_binned
 
 ```
 
diff --git a/docs/source/user/getting_started.md b/docs/source/user/getting_started.md
index ec9769e..608edc3 100644
--- a/docs/source/user/getting_started.md
+++ b/docs/source/user/getting_started.md
@@ -8,11 +8,11 @@ We didn't reinvent the wheel and rely on existing libraries and packages to perf
 
 ### GEOS
 The plotting of the seismicity requires [GEOS](https://libgeos.org/), a C/C++ library for computational geometry. If `GEOS` is not installed on your machine, you will need to get it, for example on a linux machine with 
-```terminal
+```console
 sudo apt-get libgeos-dev
 ```
 or on a mac with
-```terminal
+```console
 brew install geos
 ```
 
@@ -20,16 +20,16 @@ brew install geos
 
 ### Install from source
 This way of installing `SeismoStats` in another environement allows you to use the static version. 
-```terminal
+```console
 pip install git+https://github.com/swiss-seismological-service/SeismoStats.git
 ```
 
 If you want to install a specific branch:
-```terminal
+```console
 pip install git+https://github.com/swiss-seismological-service/SeismoStats.git@feature/branch
 ```
 
 To update your environment to the latest version of `SeismoStats`:
-```terminal
+```console
 pip install --force-reinstall git+https://github.com/swiss-seismological-service/SeismoStats.git
 ```
\ No newline at end of file
diff --git a/seismostats/analysis/__init__.py b/seismostats/analysis/__init__.py
index 63494e2..553b6a1 100644
--- a/seismostats/analysis/__init__.py
+++ b/seismostats/analysis/__init__.py
@@ -8,5 +8,8 @@
                                                 estimate_b_weichert,
                                                 estimate_b_kijko_smit,
                                                 estimate_b_more_positive,
-                                                make_more_incomplete)
+                                                make_more_incomplete,
+                                                beta_to_b_value,
+                                                b_value_to_beta
+                                                )
 from seismostats.analysis.estimate_mc import mc_ks, mc_max_curvature
diff --git a/seismostats/analysis/estimate_beta.py b/seismostats/analysis/estimate_beta.py
index 7b063f7..61583cc 100644
--- a/seismostats/analysis/estimate_beta.py
+++ b/seismostats/analysis/estimate_beta.py
@@ -404,7 +404,7 @@ def make_more_incomplete(
 
     Source:
         E. Lippiello and G. Petrillo. Journal of Geophysical Research: Solid
--       Earth, 129(2):e2023JB027849, 2024.
+        Earth, 129(2):e2023JB027849, 2024.
 
     Args:
         magnitudes: array of magnitudes, sorted in time (first
@@ -419,8 +419,7 @@ def make_more_incomplete(
         magnitudes: filtered array of magnitudes
         times:      filtered array of datetime objects
         idx:        indices of the events that were kept
-
-    """
+        """
 
     # sort magnitudes in time
     idx_sort = np.argsort(times)
diff --git a/seismostats/analysis/estimate_mc.py b/seismostats/analysis/estimate_mc.py
index 930414a..227f7e1 100644
--- a/seismostats/analysis/estimate_mc.py
+++ b/seismostats/analysis/estimate_mc.py
@@ -129,7 +129,7 @@ def ks_test_gr(
     )
 
     for ii in range(n):
-        simulated = simulated_all[n_sample * ii : n_sample * (ii + 1)]
+        simulated = simulated_all[n_sample * ii: n_sample * (ii + 1)]
         _, y_th = cdf_discrete_GR(simulated, mc=mc, delta_m=delta_m, beta=beta)
         _, y_emp = empirical_cdf(simulated)
 
@@ -303,12 +303,13 @@ def mc_max_curvature(
           catalogues: Estimating the magnitude of completeness and its
           uncertainty.
           Bulletin of the Seismological Society of America, 95(2), pp.684-698.
+
     Args:
-        sample:             Magnitudes to test
-        delta_m:            Magnitude bins (sample has to be rounded to bins
-                            beforehand)
-        correction_factor:  Correction factor for the maximum curvature method
-        (default 0.2 after Woessner & Wiemer 2005)
+        sample:     Magnitudes to test
+        delta_m:    Magnitude bins (sample has to be rounded to bins beforehand)
+            correction_factor:  Correction factor for the maximum curvature
+            method (default 0.2 after Woessner & Wiemer 2005)
+
     Returns:
         mc:                 estimated completeness magnitude
     """

From cb5e89c47ef34b2ede04a66d8281042f5f3d9c9c Mon Sep 17 00:00:00 2001
From: Aron Mirwald <aron.mirwald@sed.ethz.ch>
Date: Fri, 21 Jun 2024 15:02:27 +0200
Subject: [PATCH 15/15] include Catalog methods in docs

---
 docs/source/reference/catalog.md | 26 +++++++++++++++++++++++++-
 1 file changed, 25 insertions(+), 1 deletion(-)

diff --git a/docs/source/reference/catalog.md b/docs/source/reference/catalog.md
index bb44810..1780939 100644
--- a/docs/source/reference/catalog.md
+++ b/docs/source/reference/catalog.md
@@ -13,11 +13,35 @@
     Catalog
 ```
 
-## Transformation
+## Modify Catalog
 
 ```{eval-rst}
 .. autosummary::
     :toctree: api/
 
     Catalog.bin_magnitudes
+    Catalog.strip
+    Catalog.drop_ids
+    Catalog.drop_uncertainties
+```
+
+## Estimate from Catalog
+
+```{eval-rst}
+.. autosummary::
+    :toctree: api/
+
+    Catalog.estimate_b
+    Catalog.estimate_mc
+```
+
+## Transform from or to other format
+
+```{eval-rst}
+.. autosummary::
+    :toctree: api/
+
+    Catalog.to_quakeml
+    Catalog.from_quakeml
+    Catalog.from_dict
 ```
\ No newline at end of file