<commit_before># -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2013-2021 GEM Foundation, <NAME>
#
# OpenQuake 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.
#
# OpenQuake is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>.

"""
Module exports :class:`Lin2009`
"""
import numpy as np

from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
from openquake.hazardlib import const
from openquake.hazardlib.imt import PGA, SA


class Lin2009(GMPE):
    """
    Implements GMPE developed by <NAME> and published as "Ground-motion
    attenuation relationship and path-effect study using Taiwan Data set"
    (Ph.D. dissertation of National Central University, Taiwan).
    This class implements the equations for 'crustal events'.
    """

    #: Supported tectonic region type is active shallow crust.
    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST

    #: Supported intensity measure types are spectral acceleration,
    #: and peak ground acceleration, see Table 4.1 in pages 48-49.
    DEFINED_FOR_INTENSITY_MEASURE_TYPES = set([
        PGA,
        SA
    ])

    #: Supported intensity measure component is geometric mean
    #: of two horizontal components, see equation 4.1 page 46.
    DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.AVERAGE_HORIZONTAL

    #: Supported standard deviation types is total, see equation 4.1 page 46.
    DEFINED_FOR_STANDARD_DEVIATION_TYPES = set([
        const.StdDev.TOTAL
    ])

    #: Required site parameter is only Vs30 (used to distinguish rock).
    REQUIRES_SITES_PARAMETERS = {'vs30'}

    #: Required rupture parameters are magnitude and rake see
    #: equation 4.1 page 46.
    REQUIRES_RUPTURE_PARAMETERS = {'mag', 'rake'}

    #: Required distance measure is rupture distance, see equation 4.1
    #: page 46.
    REQUIRES_DISTANCES = {'rrup'}

    def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        """
        See :meth:`superclass method
        <.base.GroundShakingIntensityModel.get_mean_and_stddevs>`
        for spec of input and result values.
        """
        C = self.COEFFS[imt]
        mean = (
            self._get_magnitude_term(C, rup.mag) +
            self._get_distance_term(C, rup.mag, dists.rrup) +
            self._get_style_of_faulting_term(C, rup.rake) +
            self._get_site_response_term(C, sites.vs30))

        stddevs = self._get_stddevs(C, stddev_types, len(sites.vs30))
        return mean, stddevs

    def _get_magnitude_term(self, C, mag):
        """
        Returns the magnitude scaling term.
        """
        lny = C['C1'] + (C['C3'] * ((8.5 - mag) ** 2.))
        if mag > 6.3:
            return lny + (-C['H'] * C['C5']) * (mag - 6.3)
        else:
            return lny + C['C2'] * (mag - 6.3)

    def _get_distance_term(self, C, mag, rrup):
        """
        Returns the distance scaling term
        """
        return (C['C4'] + C['C5'] * (mag - 6.3)) *\
            np.log(np.sqrt(rrup ** 2. + np.exp(C['H']) ** 2.))

    def _get_style_of_faulting_term(self, C, rake):
        """
        Returns the style of faulting factor
        """
        f_n, f_r = self._get_fault_type_dummy_variables(rake)
        return C['C6'] * f_n + C['C7'] * f_r

    def _get_fault_type_dummy_variables(self, rake):
        """
        Defines the fault type dummy variables for normal faulting (f_n) and
        reverse faulting (f_r) from rake. Classification based on that
        found in the original fortran code of Lin (2009)
        """
        f_n, f_r = 0, 0
        if rake >= -120 and rake <= -60:
            # normal
            f_n = 1
        elif rake >= 30 and rake <= 150:
            # reverse
            f_r = 1
        return f_n, f_r

    def _get_site_response_term(self, C, vs30):
        """
        Returns the site amplification term
        """
        return C['C8'] * np.log(vs30 / 1130.0)

    def _get_stddevs(self, C, stddev_types, nsites):
        """
        Compute total standard deviation, see table 4.2, page 50.
        """
        stddevs = []
        for stddev_type in stddev_types:
            assert stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
            if stddev_type == const.StdDev.TOTAL:
                stddevs.append(C['sigma'] + np.zeros(nsites, dtype=float))
        return stddevs

    #: Coefficient table for rock sites, see table 3 page 227.
    COEFFS = CoeffsTable(sa_damping=5.0, table="""\
    IMT       C1      C2       C3       C4      C5       H       C6      C7       C8  sigma
    pga   1.0109  0.3822   0.0000  -1.1634  0.1722  1.5184  -0.1907  0.1322  -0.4741  0.5363
    0.01  1.0209  0.3822  -0.0003  -1.1633  0.1722  1.5184  -0.1922  0.1314  -0.4738  0.5360
    0.02  1.0416  0.3822   0.0017  -1.1668  0.1722  1.5184  -0.1942  0.1311  -0.4700  0.5349
    0.03  1.1961  0.3822   0.0038  -1.2028  0.1722  1.5184  -0.1990  0.1314  -0.4741  0.5419
    0.04  1.3834  0.3822   0.0087  -1.2499  0.1722  1.5184  -0.1959  0.1362  -0.4806  0.5506
    0.05  1.5612  0.3822   0.0153  -1.2957  0.1722  1.5184  -0.1922  0.1417  -0.4911  0.5607
    0.06  1.6907  0.3822   0.0210  -1.3218  0.1722  1.5184  -0.1984  0.1500  -0.4900  0.5718
    0.07  1.7673  0.3822   0.0261  -1.3336  0.1722  1.5184  -0.2011  0.1557  -0.4920  0.5800<commit_msg>Refactored lin_2009<commit_after># -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2013-2021 GEM Foundation, <NAME>
#
# OpenQuake 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.
#
# OpenQuake is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>.

"""
Module exports :class:`Lin2009`
"""
import numpy as np

from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
from openquake.hazardlib import const
from openquake.hazardlib.imt import PGA, SA


def _get_distance_term(C, mag, rrup):
    """
    Returns the distance scaling term
    """
    return (C['C4'] + C['C5'] * (mag - 6.3)) *\
        np.log(np.sqrt(rrup ** 2. + np.exp(C['H']) ** 2.))


def _get_fault_type_dummy_variables(rake):
    """
    Defines the fault type dummy variables for normal faulting (f_n) and
    reverse faulting (f_r) from rake. Classification based on that
    found in the original fortran code of Lin (2009)
    """
    f_n, f_r = 0, 0
    if rake >= -120 and rake <= -60:
        # normal
        f_n = 1
    elif rake >= 30 and rake <= 150:
        # reverse
        f_r = 1
    return f_n, f_r


def _get_magnitude_term(C, mag):
    """
    Returns the magnitude scaling term.
    """
    lny = C['C1'] + (C['C3'] * ((8.5 - mag) ** 2.))
    if mag > 6.3:
        return lny + (-C['H'] * C['C5']) * (mag - 6.3)
    else:
        return lny + C['C2'] * (mag - 6.3)


def _get_site_response_term(C, vs30):
    """
    Returns the site amplification term
    """
    return C['C8'] * np.log(vs30 / 1130.0)


def _get_stddevs(C, stddev_types, nsites):
    """
    Compute total standard deviation, see table 4.2, page 50.
    """
    stddevs = []
    for stddev_type in stddev_types:
        if stddev_type == const.StdDev.TOTAL:
            stddevs.append(C['sigma'] + np.zeros(nsites, dtype=float))
    return stddevs


def _get_style_of_faulting_term(C, rake):
    """
    Returns the style of faulting factor
    """
    f_n, f_r = _get_fault_type_dummy_variables(rake)
    return C['C6'] * f_n + C['C7'] * f_r


class Lin2009(GMPE):
    """
    Implements GMPE developed by <NAME> and published as "Ground-motion
    attenuation relationship and path-effect study using Taiwan Data set"
    (Ph.D. dissertation of National Central University, Taiwan).
    This class implements the equations for 'crustal events'.
    """

    #: Supported tectonic region type is active shallow crust.
    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST

    #: Supported intensity measure types are spectral acceleration,
    #: and peak ground acceleration, see Table 4.1 in pages 48-49.
    DEFINED_FOR_INTENSITY_MEASURE_TYPES = {PGA, SA}

    #: Supported intensity measure component is geometric mean
    #: of two horizontal components, see equation 4.1 page 46.
    DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.AVERAGE_HORIZONTAL

    #: Supported standard deviation types is total, see equation 4.1 page 46.
    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {const.StdDev.TOTAL}

    #: Required site parameter is only Vs30 (used to distinguish rock).
    REQUIRES_SITES_PARAMETERS = {'vs30'}

    #: Required rupture parameters are magnitude and rake see
    #: equation 4.1 page 46.
    REQUIRES_RUPTURE_PARAMETERS = {'mag', 'rake'}

    #: Required distance measure is rupture distance, see equation 4.1
    #: page 46.
    REQUIRES_DISTANCES = {'rrup'}

    def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        """
        See :meth:`superclass method
        <.base.GroundShakingIntensityModel.get_mean_and_stddevs>`
        for spec of input and result values.
        """
        C = self.COEFFS[imt]
        mean = (
            _get_magnitude_term(C, rup.mag) +
            _get_distance_term(C, rup.mag, dists.rrup) +
            _get_style_of_faulting_term(C, rup.rake) +
            _get_site_response_term(C, sites.vs30))

        stddevs = _get_stddevs(C, stddev_types, len(sites.vs30))
        return mean, stddevs

    #: Coefficient table for rock sites, see table 3 page 227.
    COEFFS = CoeffsTable(sa_damping=5.0, table="""\
    IMT       C1      C2       C3       C4      C5       H       C6      C7       C8  sigma
    pga   1.0109  0.3822   0.0000  -1.1634  0.1722  1.5184  -0.1907  0.1322  -0.4741  0.5363
    0.01  1.0209  0.3822  -0.0003  -1.1633  0.1722  1.5184  -0.1922  0.1314  -0.4738  0.5360
    0.02  1.0416  0.3822   0.0017  -1.1668  0.1722  1.5184  -0.1942  0.1311  -0.4700  0.5349
    0.03  1.1961  0.3822   0.0038  -1.2028  0.1722  1.5184  -0.1990  0.1314  -0.4741  0.5419
    0.04  1.3834  0.3822   0.0087  -1.2499  0.1722  1.5184  -0.1959  0.1362  -0.4806  0.5506
    0.05  1.5612  0.3822   0.0153  -1.2957  0.1722  1.5184  -0.1922  0.1417  -0.4911  0.5607
    0.06  1.6907  0.3822   0.0210  -1.3218  0.1722  1.5184  -0.1984  0.1500  -0.4900  0.5718
    0.07  1.7673  0.3822   0.0261  -1.3336  0.1722  1.5184  -0.2011  0.1557  -0.4920  0.5800
{'fraction_non_alphanumeric': 0.10724406659582005, 'fraction_numerical': 0.09582004959263195, 'mean_word_length': 2.4148775325068037, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3087874', 'n_tokens_mistral': 4699, 'n_tokens_neox': 4079, 'n_words': 1164}