Yield imls/IMT and poes/IMT containg mean and stddev for all sites
199538
def extract_mean_std_curves(dstore, what):
    
    getter = getters.PmapGetter(dstore)
    arr = getter.get_mean().array
    for imt in getter.imtls:
        yield 'imls/' + imt, getter.imtls[imt]
        yield 'poes/' + imt, arr[:, getter.imtls(imt)]
883791
def _get_stddevs(self, sites, rup, C, stddev_types, ln_y_ref, exp1, exp2):
        
        ret = []
        for stddev_type in stddev_types:
            assert stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
            if stddev_type == const.StdDev.TOTAL:
                # eq. 13
                ret.append(0.65 * np.ones_like(sites.vs30))
        return ret
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
845819
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)

        C = self.COEFFS[imt]
        mean = self._compute_mean(C, rup.mag, dists.rrup)
        if isinstance(imt, SA) or isinstance(imt, PGA):
            # Convert from m/s**2 to g
            mean = mean - np.log(g)
        elif isinstance(imt, PGV):  # Convert from m/s to cm/s
            mean = mean + np.log(100.0)
        stddevs = self._get_stddevs(C, stddev_types, rup.mag,
                                    dists.rrup.shape[0])

        return mean, stddevs
1106855
def get_stddevs(self, C, stddev_shape, stddev_types, sites):
        
        stddevs = []
        tau = C["tau_event"]
        sigma_s = np.zeros(sites.vs30measured.shape, dtype=float)
        sigma_s[sites.vs30measured] += C["sigma_s_obs"]
        sigma_s[np.logical_not(sites.vs30measured)] += C["sigma_s_inf"]
        phi = np.sqrt(C["phi0"] ** 2.0 + sigma_s ** 2.)
        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(tau ** 2. + phi ** 2.) +
                               np.zeros(stddev_shape))
            elif stddev_type == const.StdDev.INTRA_EVENT:
                stddevs.append(phi + np.zeros(stddev_shape))
            elif stddev_type == const.StdDev.INTER_EVENT:
                stddevs.append(tau + np.zeros(stddev_shape))
        return stddevs
1021093
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        # get the necessary set of coefficients
        C = self.COEFFS[imt]
        # compute median sa on rock (vs30=1180m/s). Used for site response
        # term calculation
        sa1180 = np.exp(self._get_sa_at_1180(C, imt, sites, rup, dists))

        # get the mean value
        mean = (self._get_basic_term(C, rup, dists) +
                self._get_faulting_style_term(C, rup) +
                self._get_site_response_term(C, imt, sites.vs30, sa1180) +
                self._get_hanging_wall_term(C, dists, rup) +
                self._get_top_of_rupture_depth_term(C, imt, rup) +
                self._get_soil_depth_term(C, sites.z1pt0 / METRES_PER_KM,
                                          sites.vs30)
                )
        mean += self._get_regional_term(C, imt, sites.vs30, dists.rrup)
        # get standard deviations
        stddevs = self._get_stddevs(C, imt, rup, sites, stddev_types, sa1180,
                                    dists)
        return mean, stddevs
669975
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)

        mean = np.zeros_like(sites.vs30)
        stddevs = [np.zeros_like(sites.vs30) for _ in stddev_types]

        idx_rock = sites.vs30 >= self.ROCK_VS30
        idx_soil = sites.vs30 < self.ROCK_VS30

        if idx_rock.any():
            C = self.COEFFS_ROCK[imt]
            self._compute_mean(C, self.CONSTS['A1_rock'],
                               self.CONSTS['A2_rock'], self.CONSTS['A3_rock'],
                               self.CONSTS['A4_rock'], self.CONSTS['A5_rock'],
                               self.CONSTS['A6_rock'], rup.mag, rup.hypo_depth,
                               dists.rrup, mean, idx_rock)
            self._compute_std(C, rup.mag, stddevs, idx_rock)

            if imt == SA(period=4.0, damping=5.0):
                mean = mean / 0.399

        if idx_soil.any():
            C = self.COEFFS_SOIL[imt]
            self._compute_mean(C, self.CONSTS['A1_soil'],
                               self.CONSTS['A2_soil'], self.CONSTS['A3_soil'],
                               self.CONSTS['A4_soil'], self.CONSTS['A5_soil'],
                               self.CONSTS['A6_soil'], rup.mag, rup.hypo_depth,
                               dists.rrup, mean, idx_soil)
            self._compute_std(C, rup.mag, stddevs, idx_soil)

        return mean, stddevs
399810
def get_mean_and_stddevs(self, sctx, rctx, dctx, imt, stddev_types):
        
        return self.kwargs[str(imt)].get_mean_and_stddevs(
            sctx, rctx, dctx, imt, stddev_types)
414232
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        C = self.COEFFS[imt]

        mag = rup.mag - 6
        d = np.sqrt(dists.rjb ** 2 + C['c7'] ** 2)
        mean = np.zeros_like(d)

        mean += C['c1'] + C['c2'] * mag + C['c3'] * mag ** 2 + C['c6']

        idx = d <= 100.
        mean[idx] = mean[idx] + C['c5'] * np.log10(d[idx])

        idx = d > 100.
        mean[idx] = (mean[idx] + C['c5'] * np.log10(100.) -
                     np.log10(d[idx] / 100.) + C['c4'] * (d[idx] - 100.))

        # convert from log10 to ln and from cm/s**2 to g
        mean = np.log((10.0 ** (mean - 2.0)) / g)

        stddevs = self._get_stddevs(C, stddev_types,  dists.rjb.shape[0])

        return mean, stddevs
382329
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        

        # assign constant
        log10e = np.log10(np.e)

        # Distance term
        R = np.sqrt(dists.rjb ** 2 + 11.29 ** 2)
        # Magnitude term
        M = rup.mag - 6

        # Site term only distinguishes between lava and ash;
        # since ash sites have Vs30 in the range 60-200m/s,
        # we use this upper value as class separator
        S = np.zeros(R.shape)
        S[sites.vs30 <= 200] = 1

        # Mean ground motion (natural log)
        # call super
        mean, stddevs = super().get_mean_and_stddevs(sites, rup, dists,
                                                     imt, stddev_types)

        if rup.mag > 7. and rup.mag <= 7.7:
            mean = (0.171 * (1 - M)) / log10e + mean

        elif rup.mag > 7.7:
            mean = (0.1512 + 0.387 * (1 - M)) / log10e + mean

        # define natural log of SA 0.3 sec and 0.2 sec
        if isinstance(imt, SA):
            if imt.period == 0.3:
                mean = np.log(2.2) + mean

            if imt.period == 0.2:
                mean = np.log(2.5) + mean

        return mean, stddevs
874304
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        R = (dists.rrup)
        M = rup.mag
        # get constants
        Ssr = self.get_Ssr_term(sites.vs30)
        Shr = self.get_Shr_term(sites.vs30)
        rake = rup.rake
        F = self.get_fault_term(rake)

        # compute mean
        mean = -3.512 + (0.904 * M) - (1.328 * np.log(np.sqrt(R**2
               + (0.149 * np.exp(0.647 * M))**2))) \
               + (1.125 - 0.112 * np.log(R) - 0.0957 * M) * F \
               + (0.440 - 0.171 * np.log(R)) * Ssr \
               + (0.405 - 0.222 * np.log(R)) * Shr
        stddevs = self.get_stddevs(mean, stddev_types)
        return mean, stddevs
714720
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)

        C = self.COEFFS[imt]
        mean = (self._get_magnitude_scaling(C, rup.mag) +
                self._get_distance_scaling(C, rup.mag, dists.rhypo))
        if imt.name in "SA PGA":
            mean = np.log(np.exp(mean) / (100.0 * g))
        stddevs = self._compute_std(C, stddev_types, len(dists.rhypo))
        return mean, stddevs
432881
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]
        # 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)

        mean = self._get_mean(sites, C, ln_y_ref, exp1, exp2, v1)
        mean += convert_to_LHC(imt)

        stddevs = self._get_stddevs(sites, rup, C, stddev_types,
                                    ln_y_ref, exp1, exp2)
        return mean, 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
93906
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        # Check that the requested standard deviation type is available
        assert all(stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
                   for stddev_type in stddev_types)
        #
        # Set parameters
        magn = rup.mag
        epi = dists.repi
        theta = dists.azimuth
        #
        # Convert Mw into Ms
        if magn < 6.58:
            mag = (magn - 0.59) / 0.86
        else:
            mag = (magn + 2.42) / 1.28
        #
        # Set coefficients
        coeff = self.COEFFS[imt]
        a1ca, a1cb, a1cc, a1cd, a1ce, a2ca, a2cb, a2cc, a2cd, a2ce = \
            gc(coeff, mag)
        #
        # Get correction coefficients. Here for each site we find the
        # the geometry of the ellipses
        ras = []
        for epi, theta in zip(dists.repi, dists.azimuth):
            res = get_ras(epi, theta, mag, coeff)
            ras.append(res)
        ras = np.array(ras)
        rbs = rbf(ras, coeff, mag)
        #
        # Compute values of ground motion for the two cases. The value of
        # 225 is hardcoded under the assumption that the hypocentral depth
        # corresponds to 15 km (i.e. 15**2)
        mean1 = (a1ca + a1cb * mag +
                 a1cc * np.log((ras**2+225)**0.5 +
                               a1cd * np.exp(a1ce * mag)))
        mean2 = (a2ca + a2cb * mag +
                 a2cc * np.log((rbs**2+225)**0.5 +
                               a2cd * np.exp(a2ce * mag)))
        #
        # Get distances
        x = (mean1 * np.sin(np.radians(dists.azimuth)))**2
        y = (mean2 * np.cos(np.radians(dists.azimuth)))**2
        mean = mean1 * mean2 / np.sqrt(x+y)
        if imt.name == "PGA":
            mean = np.exp(mean)/g/100
        elif imt.name == "PGV":
            mean = np.exp(mean)
        else:
            raise ValueError('Unsupported IMT')
        #
        # Get the standard deviation
        stddevs = self._compute_std(coeff, stddev_types, len(dists.repi))
        #
        # Return results
        return np.log(mean), stddevs
1126258
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_term(C, rup, dists.rjb, imt) +
                self.get_site_amplification(C, sites))
        # GMPE originally in cm/s/s - convert to g
        if imt.name in "PGA SA":
            mean -= np.log(100.0 * g)
        stddevs = self.get_stddevs(C, dists.rjb.shape, stddev_types, sites)
        if self.sigma_mu_epsilon:
            # Apply the epistemic uncertainty factor (sigma_mu) multiplied by
            # the number of standard deviations
            sigma_mu = self.get_sigma_mu_adjustment(C, imt, rup, dists)
            # Cap sigma_mu at 0.5 ln units
            sigma_mu[sigma_mu > 0.5] = 0.5
            # Sigma mu should not be less than the standard deviation of the
            # fault-to-fault variability
            sigma_mu[sigma_mu < C["tau_fault"]] = C["tau_fault"]
            mean += (self.sigma_mu_epsilon * sigma_mu)
        return mean, stddevs
1063321
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        C = self.COEFFS[imt]
        C_SIG = self.SIGMA_COEFFS[imt]

        mean = (self.get_magnitude_scaling_term(C, rup.mag) +
                self.get_distance_scaling_term(C, dists.rhypo))

        std_devs = self.get_stddevs(C_SIG, stddev_types, len(dists.rhypo))

        #: Mean ground motions initially returned in cm/s/s (for PGA, SA)
        #: and cm/s for PGV
        if not imt.name == "PGV":
            # Convert mean from log(cm/s/s) to g
            mean = np.log(np.exp(mean) / (100. * g))

        return mean, std_devs