<commit_before>"""
Module exports :class:`SadighEtAl1997`.
"""
from __future__ import division
from math import log, exp

from nhe.attrel.base import GMPE, CoeffsTable, SACoeffsTable
from nhe import const
from nhe.imt import PGA, SA


class SadighEtAl1997(GMPE):
    """
    Implements GMPE developed by <NAME>, <NAME>,
    <NAME>, and <NAME> (1997) and published as "Attenuation
    relationships for shallow crustal earthquakes based on California
    strong motion data", Seismological Research Letters, 68(1), 180-189.
    """
    #: Supported tectonic region type is only active shallow crust,
    #: since data consists of California earthquakes mainly.
    DEFINED_FOR_TECTONIC_REGION_TYPES = set([
        const.TRT.ACTIVE_SHALLOW_CRUST
    ])

    #: Supported intensity measure types are spectral acceleration,
    #: and peak ground acceleration, see page 180.
    DEFINED_FOR_INTENSITY_MEASURE_TYPES = set([
        PGA,
        SA
    ])

    #: Supported intensity measure component is the geometric mean
    #: of two horizontal components :attr:`~nhe.const.IMC.AVERAGE_HORIZONTAL`,
    #: see page 180.
    DEFINED_FOR_INTENSITY_MEASURE_COMPONENTS = set([
        const.IMC.AVERAGE_HORIZONTAL
    ])

    #: Supported standard deviation type is only total, see table 3.
    DEFINED_FOR_STANDARD_DEVIATION_TYPES = set([
        const.StdDev.TOTAL
    ])

    #: Required site parameter is only Vs30 (used to distinguish rock
    #: and deep soil).
    REQUIRES_SITE_PARAMETERS = set(('vs30', ))

    #: Required rupture parameters are magnitude and rake (eq. 1).
    REQUIRES_RUPTURE_PARAMETERS = set(('rake', 'mag'))

    #: Required distance measure is RRup (eq. 1).
    REQUIRES_DISTANCES = set(('rrup', ))

    #: If site vs30 is more than 750 m/s -- treat the soil as rock.
    #: See page 180.
    ROCK_VS30 = 750

    #: Magnitude value to separate coefficients table because of near field
    #: saturation effect is 6.5. See page 184.
    NEAR_FIELD_SATURATION_MAG = 6.5


    def get_mean_and_stddevs(self, ctx, imt, stddev_types, component_type):
        """
        See :meth:`superclass method
        <nhe.attrel.base.AttenuationRelationship.get_mean_and_stddevs>`
        for spec of input and result values.
        """
        is_rock = ctx.site_vs30 > self.ROCK_VS30
        mean = self._get_mean(ctx.rup_mag, ctx.rup_rake, ctx.dist_rrup,
                              is_rock, imt, component_type)
        stddevs = [self._get_stddev(ctx.rup_mag, is_rock, imt, stddev_type)
                   for stddev_type in stddev_types]
        return mean, stddevs

    def _get_mean(self, mag, rake, rrup, is_rock, imt, component_type):
        # TODO: document
        assert component_type in self.DEFINED_FOR_INTENSITY_MEASURE_COMPONENTS

        # GMPE differentiates strike-slip, reverse and normal ruptures,
        # but combines normal and strike-slip into one category. See page 180.
        is_reverse = (45 <= rake <= 135)

        if not is_rock:
            # deep soil sites, equation from table 4
            if mag <= self.NEAR_FIELD_SATURATION_MAG:
                c4 = self.COEFFS_SOIL_IMT_INDEPENDENT['c4lowmag']
                c5 = self.COEFFS_SOIL_IMT_INDEPENDENT['c5lowmag']
            else:
                c4 = self.COEFFS_SOIL_IMT_INDEPENDENT['c4himag']
                c5 = self.COEFFS_SOIL_IMT_INDEPENDENT['c5himag']
            c2 = self.COEFFS_SOIL_IMT_INDEPENDENT['c2']
            c3 = self.COEFFS_SOIL_IMT_INDEPENDENT['c3']
            C = self.COEFFS_SOIL[imt]
            if is_reverse:
                c1 = self.COEFFS_SOIL_IMT_INDEPENDENT['c1r']
                c6 = C['c6r']
            else:
                c1 = self.COEFFS_SOIL_IMT_INDEPENDENT['c1ss']
                c6 = C['c6ss']
            return (c1 + c2 * mag - c3 * log(rrup + c4 * exp(c5 * mag))
                    + c6 + C['c7'] * ((8.5 - mag) ** 2.5))

        # rock sites, equation from table 2
        if mag <= self.NEAR_FIELD_SATURATION_MAG:
            C = self.COEFFS_ROCK_LOWMAG[imt]
        else:
            C = self.COEFFS_ROCK_HIMAG[imt]
        mean = (
            C['c1'] + C['c2'] * mag + C['c3'] * ((8.5 - mag) ** 2.5)
            + C['c4'] * log(rrup + exp(C['c5'] + C['c6'] * mag))
            + C['c7'] * log(rrup + 2)
        )
        if is_reverse:
            # footnote in table 2 says that for reverse ruptures
            # the mean amplitude value should be multiplied by 1.2
            mean += 0.1823215567939546  # == log(1.2)
        return mean

    def _get_stddev(self, mag, is_rock, imt, stddev_type):
        # TODO: document
        assert stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES

        if is_rock:
            C = self.COEFFS_ROCK_STDDERR[imt]
            if mag > C['maxmag']:
                return C['maxsigma']
            else:
                return C['sigma0'] + C['magfactor'] * mag

        # footnote from table 4 says that stderr for magnitudes over 7
        # is equal to one of magnitude 7.
        if mag > 7:
            mag = 7

        C = self.COEFFS_SOIL[imt]
        return C['sigma0'] + C['magfactor'] * mag

    #: Coefficients tables for rock sites (table 2), for magnitude
    #: values of :attr:`NEAR_FIELD_SATURATION_MAG` and below. Damping
    #: for spectral acceleration here and in other SA-tables is 5%,
    #: see "introduction" section.
    #:
    #: Here we use lowest SA period of 0.075 s instead of 0.07 which
    #: is listed in table 2 for consistency with table 4.
    COEFFS_ROCK_LOWMAG = CoeffsTable("""\
        IMT        c1     c2    c3      c4      c5       c6      c7
        PGA       -0.624  1.0   0.000  -2.100   1.29649  0.250   0.0
    """).chain(SACoeffsTable("""\
    period damping c1     c2    c3      c4      c5       c6      c7
     0.075    5    0.110  1.0   0.006  -2.128   1.29649  0.250  -0.082
     0.10     5    0.275  1.0   0.006  -2.148   1.29649  0.250  -0.041
     0.20     5    0.153  1.0  -0.004  -2.080   1.29649  0.250   0.0
     0.30     5   -0.057  1.0  -0.017  -2.028   1.29649  0.250   0.0
     0.40     5   -0.298  1.0  -0.028  -1.990   1.29649  0.250   0.0<commit_msg>attrel/sadigh_1997: some refactoring and docs<commit_after>"""
Module exports :class:`SadighEtAl1997`.
"""
from __future__ import division
from math import log, exp

from nhe.attrel.base import GMPE, CoeffsTable, SACoeffsTable
from nhe import const
from nhe.imt import PGA, SA


class SadighEtAl1997(GMPE):
    """
    Implements GMPE developed by <NAME>, <NAME>,
    <NAME>, and <NAME> (1997) and published as "Attenuation
    relationships for shallow crustal earthquakes based on California
    strong motion data", Seismological Research Letters, 68(1), 180-189.
    """
    #: Supported tectonic region type is only active shallow crust,
    #: since data consists of California earthquakes mainly.
    DEFINED_FOR_TECTONIC_REGION_TYPES = set([
        const.TRT.ACTIVE_SHALLOW_CRUST
    ])

    #: Supported intensity measure types are spectral acceleration,
    #: and peak ground acceleration, see page 180.
    DEFINED_FOR_INTENSITY_MEASURE_TYPES = set([
        PGA,
        SA
    ])

    #: Supported intensity measure component is the geometric mean
    #: of two horizontal components :attr:`~nhe.const.IMC.AVERAGE_HORIZONTAL`,
    #: see page 180.
    DEFINED_FOR_INTENSITY_MEASURE_COMPONENTS = set([
        const.IMC.AVERAGE_HORIZONTAL
    ])

    #: Supported standard deviation type is only total, see table 3.
    DEFINED_FOR_STANDARD_DEVIATION_TYPES = set([
        const.StdDev.TOTAL
    ])

    #: Required site parameter is only Vs30 (used to distinguish rock
    #: and deep soil).
    REQUIRES_SITE_PARAMETERS = set(('vs30', ))

    #: Required rupture parameters are magnitude and rake (eq. 1).
    REQUIRES_RUPTURE_PARAMETERS = set(('rake', 'mag'))

    #: Required distance measure is RRup (eq. 1).
    REQUIRES_DISTANCES = set(('rrup', ))

    #: If site vs30 is more than 750 m/s -- treat the soil as rock.
    #: See page 180.
    ROCK_VS30 = 750

    #: Magnitude value to separate coefficients table because of near field
    #: saturation effect is 6.5. See page 184.
    NEAR_FIELD_SATURATION_MAG = 6.5

    def get_mean_and_stddevs(self, ctx, imt, stddev_types, component_type):
        """
        See :meth:`superclass method
        <nhe.attrel.base.AttenuationRelationship.get_mean_and_stddevs>`
        for spec of input and result values.
        """
        assert component_type in self.DEFINED_FOR_INTENSITY_MEASURE_COMPONENTS
        assert all(stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
                   for stddev_type in stddev_types)

        # GMPE differentiates strike-slip, reverse and normal ruptures,
        # but combines normal and strike-slip into one category. See page 180.
        is_reverse = (45 <= ctx.rup_rake <= 135)

        if ctx.site_vs30 > self.ROCK_VS30:
            # site is rock
            fmean = self._get_mean_rock
            fstddev = self._get_stddev_rock
        else:
            # site is deep soil
            fmean = self._get_mean_deep_soil
            fstddev = self._get_stddev_deep_soil

        mean = fmean(ctx.rup_mag, ctx.rup_rake, ctx.dist_rrup, is_reverse, imt)
        stddevs = [fstddev(ctx.rup_mag, imt) for stddev_type in stddev_types]
        return mean, stddevs

    def _get_mean_deep_soil(self, mag, rake, rrup, is_reverse, imt):
        """
        Calculate and return the mean intensity for deep soil sites.

        Implements an equation from table 4.
        """
        if mag <= self.NEAR_FIELD_SATURATION_MAG:
            c4 = self.COEFFS_SOIL_IMT_INDEPENDENT['c4lowmag']
            c5 = self.COEFFS_SOIL_IMT_INDEPENDENT['c5lowmag']
        else:
            c4 = self.COEFFS_SOIL_IMT_INDEPENDENT['c4himag']
            c5 = self.COEFFS_SOIL_IMT_INDEPENDENT['c5himag']
        c2 = self.COEFFS_SOIL_IMT_INDEPENDENT['c2']
        c3 = self.COEFFS_SOIL_IMT_INDEPENDENT['c3']
        C = self.COEFFS_SOIL[imt]
        if is_reverse:
            c1 = self.COEFFS_SOIL_IMT_INDEPENDENT['c1r']
            c6 = C['c6r']
        else:
            c1 = self.COEFFS_SOIL_IMT_INDEPENDENT['c1ss']
            c6 = C['c6ss']
        return (c1 + c2 * mag - c3 * log(rrup + c4 * exp(c5 * mag))
                + c6 + C['c7'] * ((8.5 - mag) ** 2.5))

    def _get_mean_rock(self, mag, rake, rrup, is_reverse, imt):
        """
        Calculate and return the mean intensity for rock sites.

        Implements an equation from table 2.
        """
        if mag <= self.NEAR_FIELD_SATURATION_MAG:
            C = self.COEFFS_ROCK_LOWMAG[imt]
        else:
            C = self.COEFFS_ROCK_HIMAG[imt]
        mean = (
            C['c1'] + C['c2'] * mag + C['c3'] * ((8.5 - mag) ** 2.5)
            + C['c4'] * log(rrup + exp(C['c5'] + C['c6'] * mag))
            + C['c7'] * log(rrup + 2)
        )
        if is_reverse:
            # footnote in table 2 says that for reverse ruptures
            # the mean amplitude value should be multiplied by 1.2
            mean += 0.1823215567939546  # == log(1.2)
        return mean

    def _get_stddev_rock(self, mag, imt):
        """
        Calculate and return total standard deviation for rock sites.

        Implements formulae from table 3.
        """
        C = self.COEFFS_ROCK_STDDERR[imt]
        if mag > C['maxmag']:
            return C['maxsigma']
        else:
            return C['sigma0'] + C['magfactor'] * mag

    def _get_stddev_deep_soil(self, mag, imt):
        """
        Calculate and return total standard deviation for deep soil sites.

        Implements formulae from the last column of table 4.
        """
        # footnote from table 4 says that stderr for magnitudes over 7
        # is equal to one of magnitude 7.
        if mag > 7:
            mag = 7
        C = self.COEFFS_SOIL[imt]
        return C['sigma0'] + C['magfactor'] * mag

    #: Coefficients tables for rock sites (table 2), for magnitude
    #: values of :attr:`NEAR_FIELD_SATURATION_MAG` and below. Damping
    #: for spectral acceleration here and in other SA-tables is 5%,
    #: see "introduction" section.
    #:
    #: Here we use lowest SA period of 0.075 s instead of 0.07 which
    #: is listed in table 2 for consistency with table 4.
    COEFFS_ROCK_LOWMAG = CoeffsTable("""\
        IMT        c1     c2    c3      c4      c5       c6      c7
        PGA       -0.624  1.0   0.000  -2.100   1.29649  0.250   0.0
    """).chain(SACoeffsTable("""\
    period damping c1     c2    c3      c4      c5       c6      c7
     0.075    5    0.110  1.0   0.006  -2.128   1.29649  0.250  -0.082
     0.10     5    0.275  1.0   0.006  -2.148   1.29649  0.250  -0.041
     0.20     5    0.153  1.0  -0.004  -2.080   1.29649  0.250   0.0
     0.30     5   -0.057  1.0  -0.017  -2.028   1.29649  0.250   0.0
     0.40     5   -0.298  1.0  -0.028  -1.990   1.29649  0.250   0.0
{'fraction_non_alphanumeric': 0.09739251901288219, 'fraction_numerical': 0.05416731336334006, 'mean_word_length': 2.1180740382288894, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7454872', 'n_tokens_mistral': 4906, 'n_tokens_neox': 4429, 'n_words': 1198}