<commit_before># coding: utf-8
# nhlib: A New Hazard Library
# Copyright (C) 2012 GEM Foundation
#
# This program 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.
#
# This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
"""
Module exports :class:`CauzziFaccioli2008`.
"""
from __future__ import division

import numpy as np
# standard acceleration of gravity in m/s**2
from scipy.constants import g

from nhlib.gsim.base import CoeffsTable, GMPE
from nhlib import const
from nhlib.imt import PGA, PGV, SA


class CauzziFaccioli2008(GMPE):
    """
    Implements GMPE developed by <NAME> and <NAME> and published
    as "Broadband (0.05 to 20s) prediction of displacement response spectra
    based on worldwide digital records" (Journal of Seismology, 2008, volume
    12, pages 453-475).
    This class implements the prediction equations for horizontal peak ground
    acceleration, peak ground velocity and 5%-damped spectral acceleration
    (equation 2, page 462, plus faulting style term as given in equation 5,
    page 465). Hypocentral distances are clipped at 15 km (Faccioli's personal
    communication).
    Spectral acceleration (SA) values are obtained from displacement response
    spectrum  (DSR) values (as provided by the original equations) using the
    following formula ::

        SA = DSR * 2 * π / (T ** 2)

    """

    #: Supported tectonic region type is active shallow crust,
    #: see end of 'Introduction', page 454.
    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST

    #: Supported intensity measure types are spectral acceleration, peak
    #: ground acceleration and peak ground velocity.
    #: The original paper provides coefficients for PGA and PGV, while SA
    #: is obtained from displacement response spectrum values.
    DEFINED_FOR_INTENSITY_MEASURE_TYPES = set([
        PGA,
        PGV,
        SA])

    #: Supported intensity measure component is the geometric mean of two
    #: horizontal components :attr:`~nhlib.const.IMC.AVERAGE_HORIZONTAL`,
    #: see paragraph 'On functional forms', page 462.
    DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.AVERAGE_HORIZONTAL

    #: Supported standard deviation type is only total, see paragraph 'On
    #: functional forms', page 462.
    DEFINED_FOR_STANDARD_DEVIATION_TYPES = set([
        const.StdDev.TOTAL
    ])

    #: Required site parameter is only Vs30 (used to distinguish rock
    #: and deep soils), see paragraph 'On functional forms', page 463.
    REQUIRES_SITES_PARAMETERS = set(('vs30', ))

    #: Required rupture parameters are magnitude and rake, see paragraph 'On
    #: functional forms', page 463
    REQUIRES_RUPTURE_PARAMETERS = set(('rake', 'mag'))

    #: Required distance measure is Rhypo, see paragraph 'Distance', page 456.
    REQUIRES_DISTANCES = set(('rhypo', ))

    def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        """
        See :meth:`superclass method
        <nhlib.gsim.base.GroundShakingIntensityModel.get_mean_and_stddevs>`
        for spec of input and result values.
        """
        # extract dictionaries of coefficients specific to required
        # intensity measure type
        C = self.COEFFS[imt]

        mean = self._compute_mean(C, rup.mag, dists.rhypo, sites.vs30,
                                  rup.rake, imt)

        # convert from cm/s**2 to g for SA and PGA (PGV is already in cm/s)
        # and also convert from base 10 to base e.
        if isinstance(imt, PGA):
            mean = np.log((10 ** mean) / g)
        elif isinstance(imt, SA):
            mean = np.log((10 ** mean) * ((2 * np.pi / imt.period) ** 2) *
                          1e-2 / g)
        else:
            mean = np.log(10 ** mean)

        stddevs = self._get_stddevs(C, stddev_types, sites.vs30.shape[0])

        return mean, stddevs

    def _compute_term_1_2(self, C, mag):
        """
        Compute terms 1 and 2 in equation 2 page 462.
        """
        return C['a1'] + C['a2'] * mag

    def _compute_term_3(self, C, rhypo):
        """
        Compute term 3 in equation 2 page 462.
        Distances are clipped at 15 km (as per Ezio Faccioli's personal
        communication.)
        """
        d = rhypo
        d[d <= 15.0] = 15.0

        return C['a3'] * np.log10(d)

    def _compute_site_term(self, C, vs30):
        """
        Compute site term as a function of vs30: 4th, 5th and 6th terms in
        equation 2 page 462.
        """
        # for rock values the site term is zero
        site_term = np.zeros_like(vs30)

         # hard soil
        site_term[(vs30 >= 360) & (vs30 < 800)] = C['aB']

        # medium soil
        site_term[(vs30 >= 180) & (vs30 < 360)] = C['aC']

        # soft soil
        site_term[vs30 < 180] = C['aD']

        return site_term

    def _compute_faulting_style_term(self, C, rake):
        """
        Compute faulting style term as a function of rake angle value as given
        in equation 5 page 465.
        """
        if rake > -120.0 and rake <= -60.0:
            return C['aN']
        elif rake > 30.0 and rake <= 150.0:
            return C['aR']
        else:
            return C['aS']

    def _compute_mean(self, C, mag, rhypo, vs30, rake, imt):
        """
        Compute mean value for PGV, PGA and Displacement responce spectrum,
        as given in equation 2, page 462 with the addition of the faulting
        style term as given in equation 5, page 465.
        """
        mean = (self._compute_term_1_2(C, mag) +
                self._compute_term_3(C, rhypo) +
                self._compute_site_term(C, vs30) +
                self._compute_faulting_style_term(C, rake))

        return mean

    def _get_stddevs(self, C, stddev_types, num_sites):
        """
        Return total standard deviation.
        """
        stddevs = []
        for stddev_type in stddev_types:
            assert stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
            stddevs.append(np.log(10 ** C['sigmaTM']) + np.zeros(num_sites))

        return stddevs

    #: Coefficient table constructed from the electronic suplements of the
    #: original paper.
    COEFFS = CoeffsTable(sa_damping=5, table="""\
    IMT       a1          a2         aN          aR          aS          a3          aB         aC         aD         sigmaTM
    pgv       -2.0500000  0.7710000  -0.0300000  0.0100000   0.0100000   -1.4190000  0.1600000  0.3600000  0.4800000  0.3090000
    pga       -1.2960000  0.5560000  -0.0600000  0.0940000   -0.0130000  -1.5820000  0.2200000  0.3040000  0.3320000  0.3410000
    0.05      -2.8853436  0.5241770  -0.1002050  0.1323607   -0.0104223  -1.7131617  0.1590351  0.1913741  0.1662669  0.3569631
    0.10      -1.9083649  0.4879421  -0.0893446  0.1273364   -0.0141940  -1.7141846  0.2191124  0.2180647  0.1599583  0.3698820
    0.15      -1.9913830  0.5218353  -0.0733947  0.1189461   -0.0178601  -1.6011819  0.2516196  0.3149463  0.2929666  0.3701501
    0.20      -2.1855097  0.5575823  -0.0410869  0.1014822   -0.0242622  -1.5157781  0.2795653  0.4209394  0.3950845  0.3735782
    0.25      -2.2986747  0.5918757  -0.0505207  0.0853928   -0.0144371  -1.4772492  0.2468365  0.4352964  0.4414190  0.3658115
    0.30      -2.3197219  0.6141047  -0.0423759  0.0653264   -0.0096095  -1.4608963  0.2027381  0.4325227  0.4545156  0.3574604
    0.35      -2.3532043  0.6301540  -0.0290204  0.0569969   -0.0117476  -1.4526222  0.1948161  0.4572770  0.5023281  0.3576558
    0.40      -2.4625603  0.6520124  -0.0248118  0.0514677   -0.0113052  -1.4182980  0.1736606  0.4457811  0.5335752  0.3571774
    0.45      -2.4967480  0.6641638  -0.0098857  0.0366742   -0.0112005  -1.4072543  0.1769518  0.4470774  0.5839604  0.3542876
    0.50      -2.5409604  0.6736385  -0.0050996  0.0283795   -0.0095620  -1.3894703  0.1855959  0.4482634  0.6360733  0.3585923<commit_msg>renamed sigmaTM to sigma<commit_after># coding: utf-8
# nhlib: A New Hazard Library
# Copyright (C) 2012 GEM Foundation
#
# This program 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.
#
# This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
"""
Module exports :class:`CauzziFaccioli2008`.
"""
from __future__ import division

import numpy as np
# standard acceleration of gravity in m/s**2
from scipy.constants import g

from nhlib.gsim.base import CoeffsTable, GMPE
from nhlib import const
from nhlib.imt import PGA, PGV, SA


class CauzziFaccioli2008(GMPE):
    """
    Implements GMPE developed by <NAME> and <NAME> and published
    as "Broadband (0.05 to 20s) prediction of displacement response spectra
    based on worldwide digital records" (Journal of Seismology, 2008, volume
    12, pages 453-475).
    This class implements the prediction equations for horizontal peak ground
    acceleration, peak ground velocity and 5%-damped spectral acceleration
    (equation 2, page 462, plus faulting style term as given in equation 5,
    page 465). Hypocentral distances are clipped at 15 km (Faccioli's personal
    communication).
    Spectral acceleration (SA) values are obtained from displacement response
    spectrum  (DSR) values (as provided by the original equations) using the
    following formula ::

        SA = DSR * 2 * π / (T ** 2)

    """

    #: Supported tectonic region type is active shallow crust,
    #: see end of 'Introduction', page 454.
    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.ACTIVE_SHALLOW_CRUST

    #: Supported intensity measure types are spectral acceleration, peak
    #: ground acceleration and peak ground velocity.
    #: The original paper provides coefficients for PGA and PGV, while SA
    #: is obtained from displacement response spectrum values.
    DEFINED_FOR_INTENSITY_MEASURE_TYPES = set([
        PGA,
        PGV,
        SA])

    #: Supported intensity measure component is the geometric mean of two
    #: horizontal components :attr:`~nhlib.const.IMC.AVERAGE_HORIZONTAL`,
    #: see paragraph 'On functional forms', page 462.
    DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.AVERAGE_HORIZONTAL

    #: Supported standard deviation type is only total, see paragraph 'On
    #: functional forms', page 462.
    DEFINED_FOR_STANDARD_DEVIATION_TYPES = set([
        const.StdDev.TOTAL
    ])

    #: Required site parameter is only Vs30 (used to distinguish rock
    #: and deep soils), see paragraph 'On functional forms', page 463.
    REQUIRES_SITES_PARAMETERS = set(('vs30', ))

    #: Required rupture parameters are magnitude and rake, see paragraph 'On
    #: functional forms', page 463
    REQUIRES_RUPTURE_PARAMETERS = set(('rake', 'mag'))

    #: Required distance measure is Rhypo, see paragraph 'Distance', page 456.
    REQUIRES_DISTANCES = set(('rhypo', ))

    def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        """
        See :meth:`superclass method
        <nhlib.gsim.base.GroundShakingIntensityModel.get_mean_and_stddevs>`
        for spec of input and result values.
        """
        # extract dictionaries of coefficients specific to required
        # intensity measure type
        C = self.COEFFS[imt]

        mean = self._compute_mean(C, rup.mag, dists.rhypo, sites.vs30,
                                  rup.rake, imt)

        # convert from cm/s**2 to g for SA and PGA (PGV is already in cm/s)
        # and also convert from base 10 to base e.
        if isinstance(imt, PGA):
            mean = np.log((10 ** mean) / g)
        elif isinstance(imt, SA):
            mean = np.log((10 ** mean) * ((2 * np.pi / imt.period) ** 2) *
                          1e-2 / g)
        else:
            mean = np.log(10 ** mean)

        stddevs = self._get_stddevs(C, stddev_types, sites.vs30.shape[0])

        return mean, stddevs

    def _compute_term_1_2(self, C, mag):
        """
        Compute terms 1 and 2 in equation 2 page 462.
        """
        return C['a1'] + C['a2'] * mag

    def _compute_term_3(self, C, rhypo):
        """
        Compute term 3 in equation 2 page 462.
        Distances are clipped at 15 km (as per Ezio Faccioli's personal
        communication.)
        """
        d = rhypo
        d[d <= 15.0] = 15.0

        return C['a3'] * np.log10(d)

    def _compute_site_term(self, C, vs30):
        """
        Compute site term as a function of vs30: 4th, 5th and 6th terms in
        equation 2 page 462.
        """
        # for rock values the site term is zero
        site_term = np.zeros_like(vs30)

         # hard soil
        site_term[(vs30 >= 360) & (vs30 < 800)] = C['aB']

        # medium soil
        site_term[(vs30 >= 180) & (vs30 < 360)] = C['aC']

        # soft soil
        site_term[vs30 < 180] = C['aD']

        return site_term

    def _compute_faulting_style_term(self, C, rake):
        """
        Compute faulting style term as a function of rake angle value as given
        in equation 5 page 465.
        """
        if rake > -120.0 and rake <= -60.0:
            return C['aN']
        elif rake > 30.0 and rake <= 150.0:
            return C['aR']
        else:
            return C['aS']

    def _compute_mean(self, C, mag, rhypo, vs30, rake, imt):
        """
        Compute mean value for PGV, PGA and Displacement responce spectrum,
        as given in equation 2, page 462 with the addition of the faulting
        style term as given in equation 5, page 465.
        """
        mean = (self._compute_term_1_2(C, mag) +
                self._compute_term_3(C, rhypo) +
                self._compute_site_term(C, vs30) +
                self._compute_faulting_style_term(C, rake))

        return mean

    def _get_stddevs(self, C, stddev_types, num_sites):
        """
        Return total standard deviation.
        """
        stddevs = []
        for stddev_type in stddev_types:
            assert stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
            stddevs.append(np.log(10 ** C['sigma']) + np.zeros(num_sites))

        return stddevs

    #: Coefficient table constructed from the electronic suplements of the
    #: original paper.
    COEFFS = CoeffsTable(sa_damping=5, table="""\
    IMT       a1          a2         aN          aR          aS          a3          aB         aC         aD         sigma
    pgv       -2.0500000  0.7710000  -0.0300000  0.0100000   0.0100000   -1.4190000  0.1600000  0.3600000  0.4800000  0.3090000
    pga       -1.2960000  0.5560000  -0.0600000  0.0940000   -0.0130000  -1.5820000  0.2200000  0.3040000  0.3320000  0.3410000
    0.05      -2.8853436  0.5241770  -0.1002050  0.1323607   -0.0104223  -1.7131617  0.1590351  0.1913741  0.1662669  0.3569631
    0.10      -1.9083649  0.4879421  -0.0893446  0.1273364   -0.0141940  -1.7141846  0.2191124  0.2180647  0.1599583  0.3698820
    0.15      -1.9913830  0.5218353  -0.0733947  0.1189461   -0.0178601  -1.6011819  0.2516196  0.3149463  0.2929666  0.3701501
    0.20      -2.1855097  0.5575823  -0.0410869  0.1014822   -0.0242622  -1.5157781  0.2795653  0.4209394  0.3950845  0.3735782
    0.25      -2.2986747  0.5918757  -0.0505207  0.0853928   -0.0144371  -1.4772492  0.2468365  0.4352964  0.4414190  0.3658115
    0.30      -2.3197219  0.6141047  -0.0423759  0.0653264   -0.0096095  -1.4608963  0.2027381  0.4325227  0.4545156  0.3574604
    0.35      -2.3532043  0.6301540  -0.0290204  0.0569969   -0.0117476  -1.4526222  0.1948161  0.4572770  0.5023281  0.3576558
    0.40      -2.4625603  0.6520124  -0.0248118  0.0514677   -0.0113052  -1.4182980  0.1736606  0.4457811  0.5335752  0.3571774
    0.45      -2.4967480  0.6641638  -0.0098857  0.0366742   -0.0112005  -1.4072543  0.1769518  0.4470774  0.5839604  0.3542876
    0.50      -2.5409604  0.6736385  -0.0050996  0.0283795   -0.0095620  -1.3894703  0.1855959  0.4482634  0.6360733  0.3585923
{'fraction_non_alphanumeric': 0.09201688165012185, 'fraction_numerical': 0.14016524995541818, 'mean_word_length': 2.315727236893969, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4104034', 'n_tokens_mistral': 7140, 'n_tokens_neox': 5663, 'n_words': 1730}