Returns the intra-event standard deviation at the site, as defined in
        equation 15, page 147
492941
def _compute_intra_event_std(self, C, vs30, pga1100, sigma_pga):
        
        # Get intra-event standard deviation at the base of the site profile
        sig_lnyb = np.sqrt(C['s_lny'] ** 2. - C['s_lnAF'] ** 2.)
        sig_lnab = np.sqrt(sigma_pga ** 2. - C['s_lnAF'] ** 2.)
        # Get linearised relationship between f_site and ln PGA
        alpha = self._compute_intra_event_alpha(C, vs30, pga1100)

        return np.sqrt(
            (sig_lnyb ** 2.) +
            (C['s_lnAF'] ** 2.) +
            ((alpha ** 2.) * (sig_lnab ** 2.)) +
            (2.0 * alpha * C['rho'] * sig_lnyb * sig_lnab))
248038
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        # extracting dictionary of coefficients specific to required
        # intensity measure type.
        C = self.COEFFS[imt]
        if isinstance(imt, PGA):
            imt_per = 0.0
        else:
            imt_per = imt.period
        # Fix site parameters for consistent dS2S application.
        sites.vs30 = np.array([250])
        sites.z1pt0 = np.array([330])
        # intensity on a reference soil is used for both mean
        # and stddev calculations.
        ln_y_ref = self._get_ln_y_ref(rup, dists, C)
        # exp1 and exp2 are parts of eq. 7
        exp1 = np.exp(C['phi3'] * (sites.vs30.clip(-np.inf, 1130) - 360))
        exp2 = np.exp(C['phi3'] * (1130 - 360))
        # v1 is the period dependent site term. The Vs30 above which, the
        # amplification is constant
        v1 = self._get_v1(imt)
        # Get log-mean from regular unadjusted model
        b13a_mean = self._get_mean(sites, C, ln_y_ref, exp1, exp2, v1)
        # Adjust mean and standard deviation
        mean = b13a_mean + self._get_dL2L(imt_per) + self._get_dS2S(imt_per)
        mean += convert_to_LHC(imt)
        stddevs = self._get_adjusted_stddevs(sites, rup, C, stddev_types,
                                             ln_y_ref, exp1, exp2, imt_per)

        return mean, stddevs
109913
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        # compute median PGA on rock, needed to compute non-linear site
        # amplification
        C_pga = self.COEFFS[PGA()]
        median_pga = np.exp(
            self._compute_mean(C_pga, rup.mag, dists, rup.rake)
        )

        # compute full mean value by adding nonlinear site amplification terms
        C = self.COEFFS[imt]
        mean = (self._compute_mean(C, rup.mag, dists, rup.rake) +
                self._compute_non_linear_term(C, median_pga, sites))

        stddevs = self._get_stddevs(C, stddev_types, num_sites=sites.vs30.size)

        return mean + self.adjustment_factor, stddevs
1085939
def _get_stddevs(self, C, stddev_types, stddev_shape):
        
        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["sigtot"] + np.zeros(stddev_shape))
            elif stddev_type == const.StdDev.INTRA_EVENT:
                stddevs.append(C['sig2'] + np.zeros(stddev_shape))
            elif stddev_type == const.StdDev.INTER_EVENT:
                stddevs.append(C['sig1'] + np.zeros(stddev_shape))
        return stddevs
691473
def _get_stddevs(self, C, stddev_types, num_sites):
        
        # standard deviation is converted from log10 to ln
        std_total = np.log(10 ** C['sigma'])
        stddevs = []
        for _ in stddev_types:
            stddevs.append(np.zeros(num_sites) + std_total)

        return stddevs
245561
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        assert all(stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
                   for stddev_type in stddev_types)

        # extract dictionaries of coefficients specific to required
        # intensity measure type
        C = self.COEFFS[imt]

        # Deltas for Tectonic Region Type and rake angles
        delta_R, delta_S, delta_V, delta_I = self._get_deltas(rup.rake)

        mean = self._compute_mean(C, rup.mag, dists.rrup, rup.hypo_depth,
                                  delta_R, delta_S, delta_V, delta_I,
                                  sites.vs30)

        stddevs = self._get_stddevs(C, stddev_types, sites.vs30.size)

        return mean, stddevs
1367907
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        # extract dictionaries of coefficients specific to required
        # intensity measure type
        C = self.COEFFS[imt]

        mean = self._compute_mean(C, rup, dists, sites, imt)

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

        return mean, stddevs
293024
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        # extracting dictionary of coefficients specific to required
        # intensity measure type.

        C = self.COEFFS[imt]
        mean = (self._get_magnitude_scaling(C, rup.mag) +
                self._get_distance_scaling(C, dists, rup.mag) +
                self._get_site_term(C, sites.vs30))

        # Mean is returned in terms of m/s^2. Need to convert to g
        mean -= np.log(g)
        stddevs = self.get_stddevs(C, sites.vs30.shape, stddev_types)
        return mean + self.adjustment_factor, stddevs
248809
def _get_stddevs(self, C, C_PGA, pga1100, rup, sites, stddev_types):
        
        std_intra = self._compute_intra_event_std(C, C_PGA, pga1100, rup.mag,
                                                  sites.vs30,
                                                  sites.vs30measured)
        std_inter = self._compute_inter_event_std(C, C_PGA, pga1100, rup.mag,
                                                  sites.vs30)
        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(np.sqrt(std_intra ** 2 + std_inter ** 2))
            elif stddev_type == const.StdDev.INTRA_EVENT:
                stddevs.append(std_intra)
            elif stddev_type == const.StdDev.INTER_EVENT:
                stddevs.append(std_inter)
        return stddevs
469422
def _get_stddevs(self, coeffs, stddev_types, num_sites):
        
        stddevs = []
        for stddev_type in stddev_types:
            assert stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
            stddevs.append(coeffs['sigma'] + np.zeros(num_sites))
        return np.array(stddevs)
446669
def _get_stddevs(self, C, stddev_types, rup, imt, num_sites):
        
        stddevs = []
        for stddev_type in stddev_types:
            sigma_mean = self._compute_standard_dev(rup, imt, C)
            sigma_tot = np.sqrt((sigma_mean ** 2) + (C['SigmaReg'] ** 2))
            sigma_tot = np.log10(np.exp(sigma_tot))
            stddevs.append(sigma_tot + np.zeros(num_sites))
        return stddevs
105398
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        # extracting dictionary of coefficients specific to required
        # intensity measure type.
        C = self.COEFFS_SSLAB[imt]

        # cap magnitude values at 8.0, see page 1709
        mag = rup.mag
        if mag >= 8.0:
            mag = 8.0

        # compute PGA on rock (needed for site amplification calculation)
        G = 10 ** (0.301 - 0.01 * mag)
        pga_rock = self._compute_mean(self.COEFFS_SSLAB[PGA()], G, mag,
                                      rup.hypo_depth, dists.rrup, sites.vs30,
                                      # by passing pga_rock > 500 the soil
                                      # amplification is 0
                                      np.zeros_like(sites.vs30) + 600,
                                      PGA())
        pga_rock = 10 ** (pga_rock)

        # compute actual mean and convert from log10 to ln and units from
        # cm/s**2 to g
        mean = self._compute_mean(C, G, mag, rup.hypo_depth, dists.rrup,
                                  sites.vs30, pga_rock, imt)
        mean = np.log((10 ** mean) * 1e-2 / g)

        if imt.period == 4.0:
            mean /= 0.550

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

        return mean, stddevs
1028093
def _get_intra_event_phi(self, C, mag, rjb, vs30, num_sites):
        
        base_vals = np.zeros(num_sites)
        # Magnitude Dependent phi (Equation 17)
        if mag <= 4.5:
            base_vals += C["f1"]
        elif mag >= 5.5:
            base_vals += C["f2"]
        else:
            base_vals += (C["f1"] + (C["f2"] - C["f1"]) * (mag - 4.5))
        # Distance dependent phi (Equation 16)
        idx1 = rjb > C["R2"]
        base_vals[idx1] += C["DfR"]
        idx2 = np.logical_and(rjb > C["R1"], rjb <= C["R2"])
        base_vals[idx2] += (C["DfR"] * (np.log(rjb[idx2] / C["R1"]) /
                                        np.log(C["R2"] / C["R1"])))
        # Site-dependent phi (Equation 15)
        idx1 = vs30 <= self.CONSTS["v1"]
        base_vals[idx1] -= C["DfV"]
        idx2 = np.logical_and(vs30 >= self.CONSTS["v1"],
                              vs30 <= self.CONSTS["v2"])
        base_vals[idx2] -= (
            C["DfV"] * (np.log(self.CONSTS["v2"] / vs30[idx2]) /
                        np.log(self.CONSTS["v2"] / self.CONSTS["v1"])))
        return base_vals
485759
def _get_stddevs(self, C, stddev_types, num_sites, mag, c1_rrup,
                     log_phi_ss, mean_phi_ss):
        
        phi_ss = _compute_phi_ss(C, mag, c1_rrup, log_phi_ss, mean_phi_ss)

        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(np.sqrt(
                               C['tau'] * C['tau'] +
                               phi_ss * phi_ss) +
                               np.zeros(num_sites))

            elif stddev_type == const.StdDev.INTRA_EVENT:
                stddevs.append(phi_ss + np.zeros(num_sites))

            elif stddev_type == const.StdDev.INTER_EVENT:
                stddevs.append(C['tau'] + np.zeros(num_sites))
        return stddevs
372682
def _get_mean(self, sites, rup, dists, imt, stddev_types):
        
        # Zhao et al. 2006 - Vs30 + Rrup
        mean_zh06, stds1 = super().get_mean_and_stddevs(sites, rup, dists, imt,
                                                        stddev_types)
        #
        # Atkinson and Macias (2009) - Rrup
        gmpe = AtkinsonMacias2009()
        mean_am09, stds2 = gmpe.get_mean_and_stddevs(sites, rup, dists, imt,
                                                     stddev_types)
        #
        # Abrahamson et al. (2015) - Rrup + vs30 + backarc
        gmpe = AbrahamsonEtAl2015SInter()
        mean_ab15, stds3 = gmpe.get_mean_and_stddevs(sites, rup, dists, imt,
                                                     stddev_types)
        #
        # Ghofrani and Atkinson (2014) - Rrup + vs30
        gmpe = GhofraniAtkinson2014()
        mean_ga14, stds4 = gmpe.get_mean_and_stddevs(sites, rup, dists, imt,
                                                     stddev_types)
        # Computing adjusted mean and stds
        cff = self.SITE_COEFFS[imt]
        mean_adj = (np.log(np.exp(mean_zh06)*cff['mf'])*0.1 +
                    mean_am09*0.5 + mean_ab15*0.2 +
                    np.log(np.exp(mean_ga14)*cff['mf'])*0.2)
        return mean_adj
56231
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        assert all(stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
                   for stddev_type in stddev_types)

        # Compute SA with primed coeffs and PGA with both unprimed and
        # primed coeffs
        C = self.COEFFS_PRIMED[imt]
        C_PGA = self.COEFFS_PRIMED[PGA()]
        C_PGA_unprimed = self.COEFFS_UNPRIMED[PGA()]

        # Get S term to determine if consider site term is applied
        S = self._get_site_class(sites)
        # Abrahamson and Silva (1997) hanging wall term. This is not used
        # in the latest version of GMPE but is defined in functional form in
        # the paper so we keep it here as a placeholder
        f4HW = self._compute_f4(C, rup.mag, dists.rrup)

        # Flags for rake angles
        CN, CR = self._get_fault_mechanism_flags(rup.rake)

        # Get volcanic path distance which Rvol=0 for current implementation
        # of McVerry2006Asc, but kept here as placeholder for future use
        rvol = self._get_volcanic_path_distance(dists.rrup)

        # Get delta_C and delta_D terms for site class
        delta_C, delta_D = self._get_deltas(sites)
        # Compute lnPGA_ABCD primed
        lnPGAp_ABCD = self._compute_mean(C_PGA, S, rup.mag, dists.rrup, rvol,
                                         rup.hypo_depth, CN, CR, f4HW,
                                         delta_C, delta_D)

        # Compute lnPGA_ABCD unprimed
        lnPGA_ABCD = self._compute_mean(C_PGA_unprimed, S, rup.mag, dists.rrup,
                                        rvol, rup.hypo_depth, CN, CR, f4HW,
                                        delta_C, delta_D)

        # Compute lnSA_ABCD
        lnSAp_ABCD = self._compute_mean(C, S, rup.mag, dists.rrup, rvol,
                                        rup.hypo_depth, CN, CR, f4HW,
                                        delta_C, delta_D)

        # Stage 3: Equation 6 SA_ABCD(T). This is lnSA_ABCD
        # need to calculate final lnSA_ABCD from non-log values but return log
        mean = np.log(np.exp(lnSAp_ABCD) *
                      (np.exp(lnPGA_ABCD) / np.exp(lnPGAp_ABCD)))

        # Compute standard deviations
        C_STD = self.COEFFS_STD[imt]
        stddevs = self._get_stddevs(
            C_STD, rup.mag, stddev_types, sites
        )

        return mean, stddevs
487021
def _get_stddevs(self, C, stddev_types, num_sites):
        
        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