Return C1 coeffs as function of Rrup as proposed by
    Rodriguez-Marek et al (2013)
    The C1 coeff are used to compute the single station sigma
14046
def _compute_C1_term(C, dists):
    

    c1_dists = np.zeros_like(dists)
    idx = dists < C['Rc11']
    c1_dists[idx] = C['phi_11']
    idx = (dists >= C['Rc11']) & (dists <= C['Rc21'])
    c1_dists[idx] = C['phi_11'] + (C['phi_21'] - C['phi_11']) * \
        ((dists[idx] - C['Rc11']) / (C['Rc21'] - C['Rc11']))
    idx = dists > C['Rc21']
    c1_dists[idx] = C['phi_21']
    return c1_dists
114757
def _alpha_eff(self, r_eff, n_sersic, k_eff):
        
        b = self.b_n(n_sersic)
        alpha_eff = n_sersic * r_eff * k_eff * b**(-2*n_sersic) * np.exp(b) * special.gamma(2*n_sersic)
        return -alpha_eff
1023921
def _get_distance_scaling_term(self, C, mag, rrup):
        
        return (C["r1"] + C["r2"] * mag) * np.log10(rrup + C["r3"])
1230136
def _sigma_ee_nonrel(self, gam, eps):
        
        s0 = 4 * r0 ** 2 * alpha / (15 * eps)
        x = 4 * eps / (gam ** 2 - 1)
        sigma_nonrel = s0 * self._F(x, gam)
        sigma_nonrel[np.where(eps >= 0.25 * (gam ** 2 - 1.0))] = 0.0
        sigma_nonrel[np.where(gam * np.ones_like(eps) < 1.0)] = 0.0
        return sigma_nonrel / mec2_unit
144133
def sigma_s2(self, r, R, r_ani, a, gamma, rho0_r0_gamma):
        
        beta = self._beta_ani(r, r_ani)
        return (1 - beta * R**2/r**2) * self.sigma_r2(r, a, gamma, rho0_r0_gamma, r_ani)
973720
def correlation(s, o):
    
    # s,o = filter_nan(s,o)
    if s.size == 0:
        corr = np.NaN
    else:
        corr = np.corrcoef(o, s)[0, 1]
    return corr
148616
def _compute_f2_factor(self, rrup):
        
        # f2 = max(log10(rrup/R2),0)
        f2 = np.log10(rrup / self.COEFFS_IMT_INDEPENDENT['R2'])
        f2[f2 < 0] = 0.0

        return f2
1104974
def _compute_term3(self, C, rrup):
        
        f3 = np.zeros_like(rrup)

        idx_between_70_130 = (rrup > 70) & (rrup <= 130)
        idx_greater_130 = rrup > 130

        f3[idx_between_70_130] = (
            C['c9'] * (np.log(rrup[idx_between_70_130]) - np.log(70))
        )

        f3[idx_greater_130] = (
            C['c9'] * (np.log(rrup[idx_greater_130]) - np.log(70)) +
            C['c10'] * (np.log(rrup[idx_greater_130]) - np.log(130))
        )

        return f3
315114
def sigma_r2(self, r, kwargs_profile, kwargs_anisotropy, kwargs_light):
        
        if self._mass_profile == 'power_law':
            if self._anisotropy_type == 'r_ani':
                if self._light_profile == 'Hernquist':
                    sigma_r = self.power_law_anisotropy(r, kwargs_profile, kwargs_anisotropy, kwargs_light)
                else:
                    raise ValueError('light profile %s not supported for Jeans solver' % self._light_profile)
            else:
                raise ValueError('anisotropy type %s not implemented in Jeans equation modelling' % self._anisotropy_type)
        else:
            raise ValueError('mass profile type %s not implemented in Jeans solver' % self._mass_profile)
        return sigma_r
300051
def _compute_distance_scaling(self, rup, dists, C):
        
        Mref = 4.5
        Rref = 1.0
        R = np.sqrt(dists.rjb ** 2 + C['h'] ** 2)
        return (C['c1'] + C['c2'] * (rup.mag - Mref)) * np.log(R / Rref) + \
            C['c3'] * (R - Rref)
1074451
def sigma_s2(self, r, R, kwargs_profile, kwargs_anisotropy, kwargs_light):
        
        beta = self.anisotropy.beta_r(r, kwargs_anisotropy)
        return (1 - beta * R**2/r**2) * self.sigma_r2(r, kwargs_profile, kwargs_anisotropy, kwargs_light)
893585
def covar_rescaling_factor_efficient(C):
    
    n = C.shape[0]
    P = sp.eye(n) - sp.ones((n,n))/float(n)
    CP = C - C.mean(0)[:, sp.newaxis]
    trPCP = sp.sum(P * CP)
    r = (n-1) / trPCP
    return r
302061
def crpss(self):
        
        crps_f = self.crps()
        crps_c = self.crps_climo()
        return 1.0 - float(crps_f) / float(crps_c)
518263
def _compute_stddevs(self, C, mag, rjb, imt, stddev_types):
        
        # aleatory uncertainty
        sigma_ale_m = np.interp(mag, [5.0, 5.5, 8.0],
                                [C['m50'], C['m55'], C['m80']])
        sigma_ale_rjb = np.interp(rjb, [5.0, 20.0], [C['r5'], C['r20']])
        sigma_ale = np.sqrt(sigma_ale_m ** 2 + sigma_ale_rjb ** 2)

        # epistemic uncertainty
        if imt.period < 1:
            sigma_epi = 0.36 + 0.07 * (mag - 6)
        else:
            sigma_epi = 0.34 + 0.06 * (mag - 6)

        sigma_total = np.sqrt(sigma_ale ** 2 + sigma_epi ** 2)

        stddevs = []
        for _ in stddev_types:
            stddevs.append(sigma_total)

        return stddevs
572393
def _get_distance_scaling(self, C, dists, mag):
        
        r_h = self._get_rh(C, dists)
        return (C["c1"] + C["c2"] * (mag - self.CONSTANTS["mref"])) *\
            np.log(r_h / self.CONSTANTS["rref"]) +\
            C["c3"] * (r_h - self.CONSTANTS["rref"])
605124
def _get_hanging_wall_coeffs_rrup(self, dists):
        
        fhngrrup = np.ones(len(dists.rrup))
        idx = dists.rrup > 0.0
        fhngrrup[idx] = (dists.rrup[idx] - dists.rjb[idx]) / dists.rrup[idx]
        return fhngrrup
94219
def function(self, x, y, amp, R_sersic, Re, n_sersic, gamma, e1, e2, center_x=0, center_y=0, alpha=3.):
        
        phi_G, q = param_util.ellipticity2phi_q(e1, e2)
        Rb = R_sersic
        x_shift = x - center_x
        y_shift = y - center_y

        cos_phi = np.cos(phi_G)
        sin_phi = np.sin(phi_G)

        xt1 = cos_phi*x_shift+sin_phi*y_shift
        xt2 = -sin_phi*x_shift+cos_phi*y_shift
        xt2difq2 = xt2/(q*q)
        R_ = np.sqrt(xt1*xt1+xt2*xt2difq2)
        #R_ = R_.astype(np.float32)
        if isinstance(R_, int) or isinstance(R_, float):
            R_ = max(self._smoothing, R_)
        else:
            R_[R_ < self._smoothing] = self._smoothing
        if isinstance(R_, int) or isinstance(R_, float):
            R = max(self._smoothing, R_)
        else:
            R=np.empty_like(R_)
            _R = R_[R_ > self._smoothing]  #in the SIS regime
            R[R_ <= self._smoothing] = self._smoothing
            R[R_ > self._smoothing] = _R

        k, bn = self.k_bn(n_sersic, Re)
        result = amp * (1 + (Rb / R) ** alpha) ** (gamma / alpha) * np.exp(-bn * (((R ** alpha + Rb ** alpha) / Re ** alpha) ** (1. / (alpha * n_sersic)) - 1.))
        return np.nan_to_num(result)