resolve merge conflict with main

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
 ```

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
 ```

docs/source/reference/utils.md

@@ -23,7 +23,7 @@
     :toctree: api/
 
     utils.simulate_magnitudes
-    utils.simulated_magnitudes_binned
+    utils.simulate_magnitudes_binned
 
 ```
 

docs/source/user/getting_started.md

@@ -8,28 +8,28 @@ 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
 ```
 
 ## Using SeismoStats in another code
 
 ### 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
 ```

seismostats/analysis/__init__.py

@@ -5,5 +5,11 @@
                                                 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,
+                                                beta_to_b_value,
+                                                b_value_to_beta
+                                                )
 from seismostats.analysis.estimate_mc import mc_ks, mc_max_curvature

seismostats/analysis/estimate_beta.py

@@ -253,6 +253,7 @@ def differences(magnitudes: np.ndarray) -> np.ndarray:
 def estimate_b_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,
@@ -270,6 +271,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.
+        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
                 exponential distribution [p(M) = exp(-beta*(M-mc))].
@@ -286,13 +289,95 @@ def estimate_b_positive(
         n:      number of events used for the estimation
     """
 
+    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
+    # previous one. delta_m is added to avoid numerical errors
+    mag_diffs = abs(mag_diffs[mag_diffs > dmc - delta_m / 2])
+
+    out = estimate_b_tinti(
+        mag_diffs,
+        mc=dmc,
+        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 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,
+) -> 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
+    """
+
+    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.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]
+            if mag_diff_loop > dmc - delta_m / 2:
+                mag_diffs[ii] = mag_diff_loop
+                break
+
     # 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,
-        mc=delta_m,
+        mc=dmc,
         delta_m=delta_m,
         b_parameter=b_parameter,
         return_std=return_std,
@@ -307,6 +392,61 @@ def estimate_b_positive(
         return out
 
 
+def make_more_incomplete(
+    magnitudes: np.ndarray,
+    times: np.array,
+    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
+    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).
+        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 kept are
+                returned
+
+    Returns:
+        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)
+    magnitudes = magnitudes[idx_sort]
+    times = times[idx_sort]
+
+    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_loop = magnitudes[idx_close] > magnitudes[ii]
+
+        # if there are any, remove the current event
+        if sum(idx_loop) > 0:
+            idx[ii] = False
+
+    magnitudes = magnitudes[idx]
+    times = times[idx]
+
+    if return_idx is True:
+        return magnitudes, times, idx
+
+    return magnitudes, times
+
+
 def estimate_b_laplace(
     magnitudes: np.ndarray,
     delta_m: float = 0,
@@ -344,14 +484,22 @@ 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,
         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 estimate_b_weichert(
     magnitudes: np.ndarray,

seismostats/analysis/estimate_mc.py

@@ -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
     """