Partial derivative of site amplification term with respect to
        PGA on rock (equation 26), as described in the errata and not
        in the original paper.
193498
def _compute_partial_derivative_site_amp(self, C, pga1100, vs30):
        
        delta_amp = np.zeros_like(vs30)
        vlin = C['VLIN']
        c = self.CONSTS['c']
        b = C['b']
        n = self.CONSTS['n']

        idx = vs30 < vlin
        delta_amp[idx] = (- b * pga1100[idx] / (pga1100[idx] + c) +
                          b * pga1100[idx] / (pga1100[idx] + c *
                          ((vs30[idx] / vlin) ** n)))

        return delta_amp
680935
def _compute_imt1100(self, C, sites, rup, dists, get_pga_site=False):
        
        # Calculates simple site response term assuming all sites 1100 m/s
        fsite = (C['c10'] + (C['k2'] * C['n'])) * log(1100. / C['k1'])
        # Calculates the PGA on rock
        pga1100 = np.exp(self._compute_magnitude_term(C, rup.mag) +
                         self._compute_distance_term(C, rup, dists) +
                         self._compute_style_of_faulting_term(C, rup) +
                         self._compute_hanging_wall_term(C, rup, dists) +
                         self._compute_basin_response_term(C, sites.z2pt5) +
                         fsite)
        # If PGA at the site is needed then remove factor for rock and
        # re-calculate on correct site condition
        if get_pga_site:
            pga_site = np.exp(np.log(pga1100) - fsite)
            fsite = self._compute_shallow_site_response(C, sites, pga1100)
            pga_site = np.exp(np.log(pga_site) + fsite)
        else:
            pga_site = None
        return pga1100, pga_site
17947
def _get_site_scaling_term(self, C, vs30):
        
        site_amp = C["xi"] * np.log(1200.0) * np.ones(len(vs30))
        idx = vs30 < 1200.0
        site_amp[idx] = C["xi"] * np.log(vs30[idx])
        return site_amp
1073492
def _get_mean(self, vs30, mag, rrup, imt, scale_fac):
        
        C_HR, C_BC, C_SR, SC = self._extract_coeffs(imt)

        rrup = self._clip_distances(rrup)

        f0 = self._compute_f0_factor(rrup)
        f1 = self._compute_f1_factor(rrup)
        f2 = self._compute_f2_factor(rrup)

        pga_bc = self._get_pga_bc(
            f0, f1, f2, SC, mag, rrup, vs30, scale_fac
        )

        # compute mean values for hard-rock sites (vs30 >= 2000),
        # and non-hard-rock sites (vs30 < 2000) and add soil amplification
        # term
        mean = np.zeros_like(vs30)
        self._compute_mean(C_HR, f0, f1, f2, SC, mag, rrup,
                           vs30 >= 2000.0, mean, scale_fac)
        self._compute_mean(C_BC, f0, f1, f2, SC, mag, rrup,
                           vs30 < 2000.0, mean, scale_fac)
        self._compute_soil_amplification(C_SR, vs30, pga_bc, mean)

        # convert from base 10 to base e
        if imt == PGV():
            mean = np.log(10 ** mean)
        else:
            # convert from cm/s**2 to g
            mean = np.log((10 ** mean) * 1e-2 / g)

        return mean
818968
def _compute_nonlinear_soil_term(self, C, lnSA_AB, delta_C, delta_D):
        

        lnSA_CD = (
            # line 1 equation 4 without first term (lnSA_AB)
            C['c29'] * delta_C +
            # line 2 and 3
            (C['c30as'] * np.log(np.exp(lnSA_AB) + 0.03) + C['c43']) * delta_D
            )

        return lnSA_CD
1015420
def _compute_mean(self, C, rup, dists, sites, imt):
        
        # Convert rhypo to rrup
        rrup = rhypo_to_rrup(dists.rhypo, rup.mag)
        mean = (self._get_magnitude_scaling_term(C, rup.mag) +
                self._get_distance_scaling_term(C, rup.mag, rrup) +
                self._get_style_of_faulting_term(C, rup.rake) +
                self._get_site_amplification_term(C, sites.vs30))
        # convert from cm/s**2 to g for SA and from cm/s**2 to g for 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) * ((2 * np.pi / 0.01) ** 2) *
                          1e-2 / 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)

        return mean + self.adjustment_factor
946191
def _get_site_amplification_non_linear(self, vs30, pga4nl, C):
        
        # non linear slope
        bnl = self._compute_non_linear_slope(vs30, C)
        # compute the actual non-linear term
        return self._compute_non_linear_term(pga4nl, bnl)
936233
def _compute_mean(self, C, rup, dists, sites, imt):
        
        mean = (self._get_magnitude_scaling_term(C, rup.mag) +
                self._get_distance_scaling_term(C, rup.mag, dists.rrup) +
                self._get_site_amplification_term(C, sites.vs30))
        # convert from cm/s**2 to g for SA and from m/s**2 to g for PGA (PGV
        # is already in cm/s) and also convert from base 10 to base e.
        if imt.name == "PGA":
            mean = np.log((10 ** mean) * ((2 * np.pi / 0.01) ** 2) *
                          1e-2 / g)
        elif imt.name == "SA":
            mean = np.log((10 ** mean) * ((2 * np.pi / imt.period) ** 2) *
                          1e-2 / g)
        else:
            mean = np.log(10 ** mean)

        return mean
281789
def _compute_anelastic_attenuation_term(self, C, rrup, mag):
        
        r = (rrup**2. + (C['c5'] * np.exp(C['c6'] * mag +
                                          C['c7'] * (8.5 - mag)**2.5))**2.)**.5
        f3 = ((C['c4'] + C['c13'] * mag) * np.log(r) +
              (C['c8'] + C['c12'] * mag) * r)
        return f3
1209908
private void computeDir(double[] dir, double[] fg, QNInfo qn)
      throws SurpriseConvergence {
    System.arraycopy(fg, 0, dir, 0, fg.length);

    int mmm = qn.size();
    double[] as = new double[mmm];

    for (int i = mmm - 1; i >= 0; i--) {
      as[i] = qn.getRho(i) * ArrayMath.innerProduct(qn.getS(i), dir);
      plusAndConstMult(dir, qn.getY(i), -as[i], dir);
    }

    // multiply by hessian approximation
    qn.applyInitialHessian(dir);

    for (int i = 0; i < mmm; i++) {
      double b = qn.getRho(i) * ArrayMath.innerProduct(qn.getY(i), dir);
      plusAndConstMult(dir, qn.getS(i), as[i] - b, dir);
    }

    ArrayMath.multiplyInPlace(dir, -1);
  }
149199
def _compute_non_linear_term(self, C, pga_only, sites):
        
        Vref = self.CONSTS['Vref']
        Vcon = self.CONSTS['Vcon']
        c = self.CONSTS['c']
        n = self.CONSTS['n']
        lnS = np.zeros_like(sites.vs30)

        # equation (6a)
        idx = sites.vs30 < Vref
        lnS[idx] = (
            C['sb1'] * np.log(sites.vs30[idx] / Vref) +
            C['sb2'] * np.log(
                (pga_only[idx] + c * (sites.vs30[idx] / Vref) ** n) /
                ((pga_only[idx] + c) * (sites.vs30[idx] / Vref) ** n)
            )
        )

        # equation (6b)
        idx = sites.vs30 >= Vref
        new_sites = sites.vs30[idx]
        new_sites[new_sites > Vcon] = Vcon
        lnS[idx] = C['sb1'] * np.log(new_sites / Vref)

        return lnS
1290161
def _compute_site_response_term(self, C, sites, pga1000):
        
        vs_star = sites.vs30.copy()
        vs_star[vs_star > 1000.0] = 1000.
        arg = vs_star / C["vlin"]
        site_resp_term = C["theta12"] * np.log(arg)
        # Get linear scaling term
        idx = sites.vs30 >= C["vlin"]
        site_resp_term[idx] += (C["b"] * self.CONSTS["n"] * np.log(arg[idx]))
        # Get nonlinear scaling term
        idx = np.logical_not(idx)
        site_resp_term[idx] += (
            -C["b"] * np.log(pga1000[idx] + self.CONSTS["c"]) +
            C["b"] * np.log(pga1000[idx] + self.CONSTS["c"] *
                            (arg[idx] ** self.CONSTS["n"])))
        return site_resp_term
453836
def levenberg_marquardt_update(x, a, c, damping=0.001):
    
    hess = hessian(x, a)
    return x - np.linalg.inv(hess + damping*np.diag(hess)).dot(gradient(x, a, c))
1373392
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        # extract dictionaries of coefficients specific to required
        # intensity measure type and for PGA
        C = self.COEFFS[imt]
        C_PGA = self.COEFFS[PGA()]

        # compute median pga on rock (vs30=1100), needed for site response
        # term calculation
        pga1100 = np.exp(self._compute_imt1100(PGA(), sites, rup, dists))

        mean = (self._compute_base_term(C, rup, dists) +
                self._compute_faulting_style_term(C, rup) +
                self._compute_site_response_term(C, imt, sites, pga1100) +
                self._compute_hanging_wall_term(C, dists, rup) +
                self._compute_top_of_rupture_depth_term(C, rup) +
                self._compute_large_distance_term(C, dists, rup) +
                self._compute_soil_depth_term(C, imt, sites.z1pt0, sites.vs30))

        stddevs = self._get_stddevs(C, C_PGA, pga1100, rup, sites,
                                    stddev_types)

        return mean, stddevs
407111
def _compute_inter_event_std(self, C, C_PGA, pga1100, mag, vs30):
        
        tau_0 = self._compute_std_0(C['s3'], C['s4'], mag)
        tau_b_pga = self._compute_std_0(C_PGA['s3'], C_PGA['s4'], mag)
        delta_amp = self._compute_partial_derivative_site_amp(C, pga1100, vs30)

        std_inter = np.sqrt(tau_0 ** 2 + (delta_amp ** 2) * (tau_b_pga ** 2) +
                            2 * delta_amp * tau_0 * tau_b_pga * C['rho'])

        return std_inter
721648
def _get_nonlinear_site_term(self, C, vs30, pga_rock):
        
        v_s = np.copy(vs30)
        v_s[vs30 > 760.] = 760.
        # Nonlinear controlling parameter (equation 8)
        f_2 = C["f4"] * (np.exp(C["f5"] * (v_s - 360.)) -
                         np.exp(C["f5"] * 400.))
        fnl = self.CONSTS["f1"] + f_2 * np.log((pga_rock + self.CONSTS["f3"]) /
                                               self.CONSTS["f3"])
        return fnl
677333
def _compute_geometrical_spreading_term(self, C, rrup):
        
        f2 = np.ones_like(rrup)
        idx1 = np.nonzero(rrup <= 70.)
        idx2 = np.nonzero((rrup > 70.) & (rrup <= 130.))
        idx3 = np.nonzero(rrup > 130.)

        f2[idx1] = (C['c9'] * np.log(rrup[idx1] + 4.5))
        f2[idx2] = (C['c10'] * np.log(rrup[idx2]/70.) +
                    C['c9'] * np.log(rrup[idx2] + 4.5))
        f2[idx3] = (C['c11'] * np.log(rrup[idx3]/130.) +
                    C['c10'] * np.log(rrup[idx3]/70.) +
                    C['c9'] * np.log(rrup[idx3] + 4.5))
        return f2