Root mean squared prop: See Adagrad paper for details.
494788
def rmsprop(grad, x, callback=None, num_iters=100,
            step_size=0.1, gamma=0.9, eps=10**-8):
    """"""
    avg_sq_grad = np.ones(len(x))
    for i in range(num_iters):
        g = grad(x, i)
        if callback: callback(x, i, g)
        avg_sq_grad = avg_sq_grad * gamma + g**2 * (1 - gamma)
        x = x - step_size * g/(np.sqrt(avg_sq_grad) + eps)
    return x
656512
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        C = self.COEFFS[imt]

        # clip rhypo at 10 (this is the minimum distance used in
        # deriving the equation), see page 22, this avoids singularity
        # in mean value equation
        rhypo = dists.rhypo.copy()
        rhypo[rhypo < 10] = 10

        # convert magnitude from Mblg to Mw
        mag = rup.mag * 0.98 - 0.39 if rup.mag <= 5.5 else \
              2.715 - 0.277 * rup.mag + 0.127 * rup.mag * rup.mag

        # functional form as explained in 'Youngs_fit_to_AB95lookup.doc'
        f1 = np.minimum(np.log(rhypo), np.log(70.))
        f2 = np.maximum(np.log(rhypo / 130.), 0)
        mean = (
            C['c1'] + C['c2'] * mag + C['c3'] * mag ** 2 +
            (C['c4'] + C['c5'] * mag) * f1 +
            (C['c6'] + C['c7'] * mag) * f2 +
            C['c8'] * rhypo
        )

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

        return mean, stddevs
466385
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]
        imean = (self._get_magnitude_scaling_term(C, rup.mag) +
                 self._get_distance_scaling_term(C, dists.rjb, rup.mag))
        # 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 imt.name in "SA PGA":
            mean = np.log((10.0 ** (imean - 2.0)) / g)
        else:
            mean = np.log(10 ** imean)

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

        return mean, stddevs
662688
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        # extracting dictionary of coefficients (for soil amplification)
        # specific to required intensity measure type
        C_SR = self.COEFFS_SOIL_RESPONSE[imt]

        # compute median PGA on rock (in g), needed to compute non-linear site
        # amplification
        C = self.COEFFS_AC10[PGA()]
        pga4nl = np.exp(
            self._compute_mean(C, rup.mag, dists.rjb, rup.rake)) * 1e-2 / g

        # compute full mean value by adding site amplification terms
        # (but avoiding recomputing mean on rock for PGA)
        if imt == PGA():
            mean = (np.log(pga4nl) +
                    self._get_site_amplification_linear(sites.vs30, C_SR) +
                    self._get_site_amplification_non_linear(sites.vs30, pga4nl,
                                                            C_SR))
        else:
            C = self.COEFFS_AC10[imt]
            mean = (self._compute_mean(C, rup.mag, dists.rjb, rup.rake) +
                    self._get_site_amplification_linear(sites.vs30, C_SR) +
                    self._get_site_amplification_non_linear(sites.vs30, pga4nl,
                                                            C_SR))

        # convert from cm/s**2 to g for SA (PGA is already computed in g)
        if imt.name == "SA":
            mean = np.log(np.exp(mean) * 1e-2 / g)

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

        return mean, stddevs
194842
def predict_withGradients(self, X):
        
        if X.ndim==1: X = X[None,:]
        ps = self.model.param_array.copy()
        means = []
        stds = []
        dmdxs = []
        dsdxs = []
        for s in self.hmc_samples:
            if self.model._fixes_ is None:
                self.model[:] = s
            else:
                self.model[self.model._fixes_] = s
            self.model._trigger_params_changed()
            m, v = self.model.predict(X)
            std = np.sqrt(np.clip(v, 1e-10, np.inf))
            dmdx, dvdx = self.model.predictive_gradients(X)
            dmdx = dmdx[:,:,0]
            dsdx = dvdx / (2*std)
            means.append(m)
            stds.append(std)
            dmdxs.append(dmdx)
            dsdxs.append(dsdx)
        self.model.param_array[:] = ps
        self.model._trigger_params_changed()
        return means, stds, dmdxs, dsdxs
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
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
747151
def transform(self, X):
        """"""
        if not self.is_fit:
            raise ValueError("The scaler has not been fit yet.")
        return (X-self.mean) / (self.std + 10e-7)
287606
def _pca(x, k=2):
    ""
    x = x-torch.mean(x,0)
    U,S,V = torch.svd(x.t())
    return torch.mm(x,U[:,:k])
354810
def nb_r_deriv(r, data_row):
    
    n = len(data_row)
    d = sum(digamma(data_row + r)) - n*digamma(r) + n*np.log(r/(r+np.mean(data_row)))
    return d
1113352
def gradients_X(self, dL_dK, X, X2, target):
        """"""
        self._K_computations(X, X2)
        arg = self._K_poly_arg
        if X2 is None:
            target += 2*self.weight_variance*self.degree*self.variance*(((X[None,:, :])) *(arg**(self.degree-1))[:, :, None]*dL_dK[:, :, None]).sum(1)
        else:
            target += self.weight_variance*self.degree*self.variance*(((X2[None,:, :])) *(arg**(self.degree-1))[:, :, None]*dL_dK[:, :, None]).sum(1)
1328626
def _hp_ars(x,params):
    """"""
    mO, sO2= params
    return -(x-mO)/sO2+1./x
284425
def _mdens(self,m):
        """"""
        return (self.a/m)**self.alpha/(1.+m/self.a)**(self.betaminusalpha)
395359
def coreBurkAlpha(self, R, Rs, rho0, r_core, ax_x, ax_y):
        
        x = R * Rs ** -1
        p = Rs * r_core ** -1

        gx = self._G(x, p)

        a = 2 * rho0 * Rs ** 2 * gx / x

        return a * ax_x / R, a * ax_y / R
1260282
def _compute_acq_withGradients(self, x):
        
        means, stds, dmdxs, dsdxs = self.model.predict_withGradients(x)
        fmins = self.model.get_fmin()
        f_acqu = None
        df_acqu = None
        for m, s, fmin, dmdx, dsdx in zip(means, stds, fmins, dmdxs, dsdxs):
            phi, Phi, u = get_quantiles(self.jitter, fmin, m, s)
            f = Phi
            df = -(phi/s)* (dmdx + dsdx * u)
            if f_acqu is None:
                f_acqu = f
                df_acqu = df
            else:
                f_acqu += f
                df_acqu += df
        return f_acqu/(len(means)), df_acqu/(len(means))
1168549
def regress_components( data, components, remove_mean=True ):
    
    timepoints = data.shape[0]
    betas = np.linalg.pinv(components).dot(data)
    if remove_mean:
        datahat = components.dot(betas)
    else:  # disregard the first layer of X, which is degree 0
        datahat = components[:, 1:].dot(betas[1:, ...])
    regressed_data = data - datahat
    return regressed_data
1134771
def smeftpar(scale, C, basis):
    """"""
    # start with a zero dict and update it with the input values
    MW = p['m_W']
    # MZ = p['m_Z']
    GF = p['GF']
    Mh = p['m_h']
    vb = sqrt(1 / sqrt(2) / GF)
    v = vb  # TODO
    _d = vMh2_to_m2Lambda(v=v, Mh2=Mh**2, C=C)
    m2 = _d['m2'].real
    Lambda = _d['Lambda'].real
    gsbar = sqrt(4 * pi * p['alpha_s'])
    gs = (1 - C['phiG'] * (v**2)) * gsbar
    gbar = 2 * MW / v
    g = gbar * (1 - C['phiW'] * (v**2))
    ebar = sqrt(4 * pi * p['alpha_e'])
    gp = get_gpbar(ebar, gbar, v, C)
    c = {}
    c['m2'] = m2
    c['Lambda'] = Lambda
    c['g'] = g
    c['gp'] = gp
    c['gs'] = gs
    K = ckmutil.ckm.ckm_tree(p['Vus'], p['Vub'], p['Vcb'], p['delta'])
    if basis == 'Warsaw':
        Mu = K.conj().T @ np.diag([p['m_u'], p['m_c'], p['m_t']])
        Md = np.diag([p['m_d'], p['m_s'], p['m_b']])
    elif basis == 'Warsaw up':
        Mu = np.diag([p['m_u'], p['m_c'], p['m_t']])
        Md = K @ np.diag([p['m_d'], p['m_s'], p['m_b']])
    else:
        raise ValueError("Basis '{}' not supported".format(basis))
    Me = np.diag([p['m_e'], p['m_mu'], p['m_tau']])
    c['Gd'] = Md / (v / sqrt(2)) + C['dphi'] * (v**2) / 2
    c['Gu'] = Mu / (v / sqrt(2)) + C['uphi'] * (v**2) / 2
    c['Ge'] = Me / (v / sqrt(2)) + C['ephi'] * (v**2) / 2
    return c