Source code for openquake.hmtk.seismicity.max_magnitude.kijko_sellevol_bayes

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4

#
# LICENSE
#
# Copyright (C) 2010-2023 GEM Foundation, G. Weatherill, M. Pagani, D. Monelli
#
# The Hazard Modeller's Toolkit (openquake.hmtk) is free software: you can
# redistribute it and/or modify it under the terms of the GNU Affero General
# Public License as published by the Free Software Foundation, either version 3
# of the License, or (at your option) any later version.
#
# You should have received a copy of the GNU Affero General Public License
# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>
#
# DISCLAIMER
#
# The software Hazard Modeller's Toolkit (openquake.hmtk) provided herein is
# released as a prototype implementation on behalf of scientists and engineers
# working within the GEM Foundation (Global Earthquake Model).
#
# It is distributed for the purpose of open collaboration and in the hope that
# it will be useful to the scientific, engineering, disaster risk and software
# design communities.
#
# The software is NOT distributed as part of GEM's OpenQuake suite
# (https://www.globalquakemodel.org/tools-products) and must be considered as a
# separate entity. The software provided herein is designed and implemented
# by scientific staff. It is not developed to the design standards, nor
# subject to same level of critical review by professional software developers,
# as GEM's OpenQuake software suite.
#
# Feedback and contribution to the software is welcome, and can be directed to
# the hazard scientific staff of the GEM Model Facility
# (hazard@globalquakemodel.org).
#
# The Hazard Modeller's Toolkit (openquake.hmtk) is therefore distributed WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# The GEM Foundation, and the authors of the software, assume no liability for
# use of the software.

"""
Module :mod:`openquake.hmtk.seismicity.max_magnitude.kijko_sellevol_bayes`
implements the Kijko & Sellevol (1989) method for estimating maximum magnitude
from observed seismicity with uncertain b-value
"""

import numpy as np
from math import fabs
from scipy.integrate import quadrature
from openquake.hmtk.seismicity.max_magnitude.base import (
    BaseMaximumMagnitude,
    MAX_MAGNITUDE_METHODS,
    _get_observed_mmax,
    _get_magnitude_vector_properties,
)


[docs]def check_config(config, data):
    """Check config file inputs

    :param dict config:
         Configuration settings for the function

    """
    essential_keys = ["input_mmin", "b-value", "sigma-b"]
    for key in essential_keys:
        if key not in config:
            raise ValueError(
                "For KijkoSellevolBayes the key %s needs to "
                "be set in the configuation" % key
            )
    if "tolerance" not in config.keys() or not config["tolerance"]:
        config["tolerance"] = 1e-5

    if not config.get("maximum_iterations", False):
        config["maximum_iterations"] = 1000

    if config["input_mmin"] < np.min(data["magnitude"]):
        config["input_mmin"] = np.min(data["magnitude"])

    if fabs(config["sigma-b"] < 1e-15):
        raise ValueError("Sigma-b must be greater than zero!")

    return config


[docs]@MAX_MAGNITUDE_METHODS.add(
    "get_mmax",
    **{
        "input_mmin": lambda cat: np.min(cat.data["magnitude"]),
        "input_mmax": lambda cat: cat.data["magnitude"][
            np.argmax(cat.data["magnitude"])
        ],
        "input_mmax_uncertainty": lambda cat: cat.get_observed_mmax_sigma(0.2),
        "b-value": float,
        "sigma-b": float,
        "maximum_iterations": 1000,
        "tolerance": 1e-5,
    },
)
class KijkoSellevolBayes(BaseMaximumMagnitude):
    """
    Class to implement Kijko & Sellevol Bayesian estimator of Mmax, with
    uncertain b-value
    """

[docs]    def get_mmax(self, catalogue, config):
        """Calculate maximum magnitude

        :returns: **mmax** Maximum magnitude and **mmax_sig** corresponding uncertainty
        :rtype: Float
        """
        # Check configuration file
        config = check_config(config, catalogue.data)
        # Negative b-values will return nan - this simply skips the integral
        if config["b-value"] <= 0.0:
            return np.nan, np.nan

        obsmax, obsmaxsig = _get_observed_mmax(catalogue.data, config)

        beta = config["b-value"] * np.log(10.0)
        sigbeta = config["sigma-b"] * np.log(10.0)

        neq, mmin = _get_magnitude_vector_properties(catalogue.data, config)

        pval = beta / (sigbeta**2.0)
        qval = (beta / sigbeta) ** 2.0

        mmax = np.copy(obsmax)
        d_t = np.inf
        iterator = 0
        while d_t > config["tolerance"]:
            rval = pval / (pval + mmax - mmin)
            ldelt = (1.0 / (1.0 - (rval**qval))) ** neq
            delta = (
                ldelt
                * quadrature(
                    self._ksb_intfunc, mmin, mmax, args=(neq, mmin, pval, qval)
                )[0]
            )

            tmmax = obsmax + delta
            d_t = np.abs(tmmax - mmax)
            mmax = np.copy(tmmax)
            iterator += 1
            if iterator > config["maximum_iterations"]:
                print(
                    "Kijko-Sellevol-Bayes estimator reached"
                    "maximum # of iterations"
                )
                d_t = -np.inf

        return mmax.item(), np.sqrt(obsmaxsig**2.0 + delta**2.0)

    def _ksb_intfunc(self, mval, neq, mmin, pval, qval):
        """
        Integral function inside Kijko-Sellevol-Bayes estimator
        (part of Eq. 10 in Kijko, 2004 - section 3.2)

        :param float mval:
            Magnitude
        :param float neq:
            Number of Earthquakes
        :param float mmin:
            Minimum Magnitude
        :param float pval:
            p-value (see Kijko, 2004 - section 3.2)
        :param float qval:
            q-value (see Kijki, 2004 - section 3.2)
        :returns:
            Output of function integrand
        """
        func1 = (1.0 - ((pval / (pval + mval - mmin)) ** qval)) ** neq
        return func1