Returns the geometric attenuation term defined in equation 3
555461
def _get_geometric_attenuation_term(self, C, mag, rrup):
        
        return (C["c5"] + C["c6"] * mag) * np.log(np.sqrt((rrup ** 2.) +
                                                          (C["c7"] ** 2.)))
114947
def highpass(cutoff):
  
  R = thub(exp(cutoff - pi), 2)
  return (1 - R) / (1 + R * z ** -1)
1310427
def pgv(self, centimeters=False):
        """"""
        if self._pgv is None:
            # Compute transfer function -- only over non-zero frequencies
            mask = ~np.isclose(self.angular_freqs, 0)
            tf = np.zeros_like(mask, dtype=np.complex)
            tf[mask] = 1 / (self.angular_freqs[mask] * 1j)
            if centimeters:
                self._pgv = GRAVITY * 100 * self.calc_peak(tf)
            else:
                self._pgv = GRAVITY * self.calc_peak(tf)
        return self._pgv
1354018
def _jmomentsurfaceIntegrand(vz,vR,vT,R,z,df,sigmaR1,gamma,sigmaz1,n,m,o): #pragma: no cover because this is too slow; a warning is shown
    """"""
    return df(R,vR*sigmaR1,vT*sigmaR1*gamma,z,vz*sigmaz1,use_physical=False,
              func= (lambda x,y,z: x**n*y**m*z**o))
1321471
def geometricBar(weights, alldistribT):
    """"""
    assert(len(weights) == alldistribT.shape[1])
    return np.exp(np.dot(np.log(alldistribT), weights.T))
862422
def lowpass(cutoff):
  
  R = thub(exp(cutoff - pi), 2)
  G = (R + 1) / 2
  return G * (1 + z ** -1) / (1 + R * z ** -1)
934869
def hessian(self, x, y, Rs, theta_Rs, r_core, center_x=0, center_y=0):

        

        if Rs < 0.0001:
            Rs = 0.0001
        x_ = x - center_x
        y_ = y - center_y
        R = np.sqrt(x_ ** 2 + y_ ** 2)

        rho0 = self._alpha2rho0(theta_Rs=theta_Rs, Rs=Rs, r_core=r_core)

        kappa = self.density_2d(x_, y_, Rs, rho0, r_core)

        gamma1, gamma2 = self.cBurkGamma(R, Rs, rho0, r_core, x_, y_)
        f_xx = kappa + gamma1
        f_yy = kappa - gamma1
        f_xy = gamma2
        return f_xx, f_yy, f_xy
1128890
def _AtildeTerm(N, h, m, k, dW, Km0, Pm0):
    """"""
    M = m*(m-1)//2
    Xk = np.random.normal(0.0, 1.0, (N, m, 1))
    Yk = np.random.normal(0.0, 1.0, (N, m, 1))
    factor1 = np.dot(Km0, Pm0 - np.eye(m**2))
    factor1 = broadcast_to(factor1, (N, M, m**2))
    factor2 = _kp(Yk + np.sqrt(2.0/h)*dW, Xk)
    return _dot(factor1, factor2)/k
1346151
def equation(self):
        
        mat = np.empty((self.nunknowns, self.model.neq))
        rhs = np.zeros(self.nunknowns)  # Needs to be initialized to zero
        for icp in range(self.ncp):
            istart = icp * self.nlayers
            ieq = 0
            for e in self.model.elementlist:
                if e.nunknowns > 0:
                    mat[istart:istart + self.nlayers, ieq:ieq + e.nunknowns] = \
                        e.potinflayers(self.xcin[icp], self.ycin[icp], self.layers, aq=self.aqin) / self.aqin.Tcol - \
                        e.potinflayers(self.xcout[icp], self.ycout[icp], self.layers, aq=self.aqout) / self.aqout.Tcol
                    ieq += e.nunknowns
                else:
                    rhs[istart:istart + self.nlayers] -= \
                        e.potentiallayers(self.xcin[icp], self.ycin[icp], self.layers, aq=self.aqin) / self.aqin.T - \
                        e.potentiallayers(self.xcout[icp], self.ycout[icp], self.layers, aq=self.aqout) / self.aqout.T
        return mat, rhs
156632
def _get_msge_with_gradient(data, delta, xvschema, skipstep, p):
    
    t, m, l = data.shape

    n = (l - p) * t
    underdetermined = n < m * p

    if underdetermined:
        return _msge_with_gradient_underdetermined(data, delta, xvschema,
                                                   skipstep, p)
    else:
        return _msge_with_gradient_overdetermined(data, delta, xvschema,
                                                  skipstep, p)
159423
def get_magnitude_term(self, C, rup):
        
        b0, stress_drop = self._get_sof_terms(C, rup.rake)
        if rup.mag <= C["m1"]:
            return b0
        else:
            # Calculate moment (equation 5)
            m_0 = 10.0 ** (1.5 * rup.mag + 16.05)
            # Get stress-drop scaling (equation 6)
            if rup.mag > C["m2"]:
                stress_drop += (C["b2"] * (C["m2"] - self.CONSTANTS["mstar"]) +
                                (C["b3"] * (rup.mag - C["m2"])))
            else:
                stress_drop += (C["b2"] * (rup.mag - self.CONSTANTS["mstar"]))
            stress_drop = np.exp(stress_drop)
            # Get corner frequency (equation 4)
            f0 = 4.9 * 1.0E6 * 3.2 * ((stress_drop / m_0) ** (1. / 3.))
            return 1. / f0
1146634
def _get_distance_term(self, C, rrup, backarc):
        
        # Geometric attenuation function
        distance_scale = -np.log10(np.sqrt(rrup ** 2 + 3600.0))
        # Anelastic attenuation in the backarc
        distance_scale[backarc] += (C["c2"] * rrup[backarc])
        # Anelastic Attenuation in the forearc
        idx = np.logical_not(backarc)
        distance_scale[idx] += (C["c1"] * rrup[idx])
        return distance_scale
1273455
def _ncc_c_3dim(x, y):
    
    den = norm(x, axis=1)[:, None] * norm(y, axis=1)
    den[den == 0] = np.Inf
    x_len = x.shape[-1]
    fft_size = 1 << (2*x_len-1).bit_length()
    cc = ifft(fft(x, fft_size) * np.conj(fft(y, fft_size))[:, None])
    cc = np.concatenate((cc[:,:,-(x_len-1):], cc[:,:,:x_len]), axis=2)
    return np.real(cc) / den.T[:, :, None]
715444
def ndtr(a):
    
    sqrth = math.sqrt(2) / 2
    x = float(a) * sqrth
    z = abs(x)
    if z < sqrth:
        y = 0.5 + 0.5 * math.erf(x)
    else:
        y = 0.5 * math.erfc(z)
        if x > 0:
            y = 1 - y
    return y
391709
def _compute_attenuation(self, rup, dists, imt, C):
        

        vec = np.ones(len(dists.rrup))

        a1 = (np.log10(np.sqrt(dists.rrup ** 2.0 + C['c11'] ** 2.0)),
              np.log10(70. * vec))

        a = np.column_stack([a1[0], a1[1]])

        b3 = (np.log10(np.sqrt(dists.rrup ** 2.0 + C['c11'] ** 2.0) /
                      (70. * vec)),
              np.log10((140. / 70.) * vec))

        b2 = np.column_stack([b3[0], b3[1]])
        b1 = ([np.min(b2, axis=1), 0. * vec])
        b = np.column_stack([b1[0], b1[1]])

        c1 = (np.log10(np.sqrt(dists.rrup ** 2.0 + C['c11'] ** 2.0) /
              (140.) * vec), 0. * vec)
        c = np.column_stack([c1[0], c1[1]])

        return (((C['c4'] + C['c5'] * rup.mag) * np.min(a, axis=1)) +
                ((C['c6'] + C['c7'] * rup.mag) * np.max(b, axis=1)) +
                ((C['c8'] + C['c9'] * rup.mag) * np.max(c, axis=1)))
823255
def equation(self):
        
        mat = np.empty((self.nunknowns, self.model.neq))
        rhs = np.zeros(self.nunknowns)  # Needs to be initialized to zero
        for icp in range(self.ncp):
            istart = icp * self.nlayers
            ieq = 0
            for e in self.model.elementlist:
                if e.nunknowns > 0:
                    qx, qy = e.disvecinflayers(self.xcout[icp], self.ycout[icp], self.layers)
                    mat[istart:istart + self.nlayers, ieq:ieq + e.nunknowns] = \
                        qx * self.cosnorm[icp] + qy * self.sinnorm[icp]
                    ieq += e.nunknowns
                else:
                    qx, qy = e.disveclayers(self.xcout[icp], self.ycout[icp], self.layers)
                    rhs[istart:istart + self.nlayers] -= qx * self.cosnorm[icp] + qy * self.sinnorm[icp]
        return mat, rhs
815051
def get_distance_term(self, C, rup, rjb, imt):
        
        h = self._get_h(C, rup.hypo_depth)
        rval = np.sqrt(rjb ** 2. + h ** 2.)
        c3 = self.get_distance_coefficients(C, imt)

        f_r = (C["c1"] + C["c2"] * (rup.mag - self.CONSTANTS["Mref"])) *\
            np.log(rval / self.CONSTANTS["Rref"]) +\
            c3 * (rval - self.CONSTANTS["Rref"])
        return f_r