de-projection of the Sersic profile based on
        Prugniel & Simien (1997)
        :return:
172657
def density(self, x, y, n_sersic, r_eff, k_eff, center_x=0, center_y=0):
        
        raise ValueError("not implemented! Use a Multi-Gaussian-component decomposition.")
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
591775
protected void computePsiDataPsiGradient(GrayF32 image1, GrayF32 image2,
											 GrayF32 deriv1x, GrayF32 deriv1y,
											 GrayF32 deriv2x, GrayF32 deriv2y,
											 GrayF32 deriv2xx, GrayF32 deriv2yy, GrayF32 deriv2xy,
											 GrayF32 du, GrayF32 dv,
											 GrayF32 psiData, GrayF32 psiGradient ) {
		int N = image1.width * image1.height;

		for( int i = 0; i < N; i++ ) {

			float du_ = du.data[i];
			float dv_ = dv.data[i];

			// compute Psi-data
			float taylor2 = image2.data[i] + deriv2x.data[i]*du_ + deriv2y.data[i]*dv_;
			float v = taylor2 - image1.data[i];

			psiData.data[i] = (float)(1.0/(2.0*Math.sqrt(v*v + EPSILON*EPSILON)));

			// compute Psi-gradient
			float dIx = deriv2x.data[i] + deriv2xx.data[i]*du_ + deriv2xy.data[i]*dv_ - deriv1x.data[i];
			float dIy = deriv2y.data[i] + deriv2xy.data[i]*du_ + deriv2yy.data[i]*dv_ - deriv1y.data[i];
			float dI2 = dIx*dIx +  dIy*dIy;

			psiGradient.data[i] = (float)(1.0/(2.0*Math.sqrt(dI2 + EPSILON*EPSILON)));
		}
	}
1055341
def _analytical_fit_adjacent(self, ref_dists):
        
        dists = ref_dists**2
        rot_dists = self._rotate_rows(dists)
        b = dists - rot_dists + self._partial_b
        self._b = b
        return self._pinvA.dot(b)
926181
def _bss_source_crit(s_true, e_spat, e_interf, e_artif):
    
    # energy ratios
    s_filt = s_true + e_spat
    sdr = _safe_db(np.sum(s_filt**2), np.sum((e_interf + e_artif)**2))
    sir = _safe_db(np.sum(s_filt**2), np.sum(e_interf**2))
    sar = _safe_db(np.sum((s_filt + e_interf)**2), np.sum(e_artif**2))
    return (sdr, sir, sar)
770538
def _calc_covar_matrix(self, profile):
        
        corr = self._calc_corr(profile)
        std = self._calc_ln_std(profile)
        # Modify the standard deviation by the truncated norm scale
        std *= randnorm.scale

        var = std ** 2
        covar = corr * std[:-1] * std[1:]

        # Main diagonal is the variance
        mat = diags([covar, var, covar], [-1, 0, 1]).toarray()

        return mat
420853
def _laplace_fit(self,obj_type):
        

        # Get Mode and Inverse Hessian information
        y = self.fit(method='PML',printer=False)

        if y.ihessian is None:
            raise Exception("No Hessian information - Laplace approximation cannot be performed")
        else:

            self.latent_variables.estimation_method = 'Laplace'

            theta, Y, scores, states, states_var, X_names = self._categorize_model_output(self.latent_variables.get_z_values())

            # Change this in future
            try:
                latent_variables_store = self.latent_variables.copy()
            except:
                latent_variables_store = self.latent_variables

            return LaplaceResults(data_name=self.data_name,X_names=X_names,model_name=self.model_name,
                model_type=self.model_type, latent_variables=latent_variables_store,data=Y,index=self.index,
                multivariate_model=self.multivariate_model,objective_object=obj_type, 
                method='Laplace',ihessian=y.ihessian,signal=theta,scores=scores,
                z_hide=self._z_hide,max_lag=self.max_lag,states=states,states_var=states_var)
1274856
def function(self, x, y, amp, R_sersic, n_sersic, e1, e2, center_x=0, center_y=0):
        
        #if n_sersic < 0.2:
        #    n_sersic = 0.2
        #if R_sersic < 10.**(-6):
        #    R_sersic = 10.**(-6)
        R_sersic = np.maximum(0, R_sersic)
        phi_G, q = param_util.ellipticity2phi_q(e1, e2)
        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)
        if isinstance(R_, int) or isinstance(R_, float):
            R_ = max(self._smoothing, R_)
        else:
            R_[R_ < self._smoothing] = self._smoothing
        k, bn = self.k_bn(n_sersic, R_sersic)
        R_frac = R_/R_sersic
        R_frac = R_frac.astype(np.float32)
        if isinstance(R_, int) or isinstance(R_, float):
            if R_frac > 100:
                result = 0
            else:
                exponent = -bn*(R_frac**(1./n_sersic)-1.)
                result = amp * np.exp(exponent)
        else:
            R_frac_real = R_frac[R_frac <= 100]
            exponent = -bn*(R_frac_real**(1./n_sersic)-1.)
            result = np.zeros_like(R_)
            result[R_frac <= 100] = amp * np.exp(exponent)
        return np.nan_to_num(result)
1185273
def spectral_centroid_and_spread(data, fs):
    """"""
    data = np.mean(data, axis=1)

    nFFT = len(data) // 2
    X = FFT(data, nFFT)

    ind = (np.arange(1, len(X) + 1)) * (fs / (2.0 * len(X)))

    Xt = X.copy()
    Xt = Xt / Xt.max()
    NUM = np.sum(ind * Xt)
    DEN = np.sum(Xt) + EPSILON

    # Centroid:
    C = (NUM / DEN)

    # Spread:
    S = np.sqrt(np.sum(((ind - C) ** 2) * Xt) / DEN)

    # Normalize:
    C = C / (fs / 2.0)
    S = S / (fs / 2.0)

    return (C, S)
644191
def least_squares(Cui, X, Y, regularization, num_threads=0):
    
    users, n_factors = X.shape
    YtY = Y.T.dot(Y)

    for u in range(users):
        X[u] = user_factor(Y, YtY, Cui, u, regularization, n_factors)
213506
def _calc_distortion(self):
        
        m = self._X.shape[0]
        self.distortion = 1/m * sum(
            linalg.norm(self._X[i, :] - self.centroids[self.clusters[i]])**2 for i in range(m)
        )
        return self.distortion
434337
function hingeLoss(p, weights) {
  const [label, features] = p;
  const grad = [];
  let dotProd = 0;

  for (let i = 0; i < features.length; i++)
    dotProd += (features[i] || 0) * (weights[i] || 0);

  if (label * dotProd < 1)
    for (let i = 0; i < features.length; i++)
      grad[i] = -label * (features[i] || 0);
  else
    for (let i = 0; i < features.length; i++)
      grad[i] = 0;

  return grad;
}
1049933
def pcolormesh(x, y, z, ax, infer_intervals=None, **kwargs):
    

    # decide on a default for infer_intervals (GH781)
    x = np.asarray(x)
    if infer_intervals is None:
        if hasattr(ax, 'projection'):
            if len(x.shape) == 1:
                infer_intervals = True
            else:
                infer_intervals = False
        else:
            infer_intervals = True

    if (infer_intervals and
            ((np.shape(x)[0] == np.shape(z)[1]) or
             ((x.ndim > 1) and (np.shape(x)[1] == np.shape(z)[1])))):
        if len(x.shape) == 1:
            x = _infer_interval_breaks(x, check_monotonic=True)
        else:
            # we have to infer the intervals on both axes
            x = _infer_interval_breaks(x, axis=1)
            x = _infer_interval_breaks(x, axis=0)

    if (infer_intervals and
            (np.shape(y)[0] == np.shape(z)[0])):
        if len(y.shape) == 1:
            y = _infer_interval_breaks(y, check_monotonic=True)
        else:
            # we have to infer the intervals on both axes
            y = _infer_interval_breaks(y, axis=1)
            y = _infer_interval_breaks(y, axis=0)

    primitive = ax.pcolormesh(x, y, z, **kwargs)

    # by default, pcolormesh picks "round" values for bounds
    # this results in ugly looking plots with lots of surrounding whitespace
    if not hasattr(ax, 'projection') and x.ndim == 1 and y.ndim == 1:
        # not a cartopy geoaxis
        ax.set_xlim(x[0], x[-1])
        ax.set_ylim(y[0], y[-1])

    return primitive
1113646
def balance_ss_model(F,L,Qc,H,Pinf,P0,dF=None,dQc=None,dPinf=None,dP0=None):
    

    bF,T,T_inv = balance_matrix(F)

    bL = np.dot( T_inv, L)
    bQc = Qc # not affected

    bH = np.dot(H, T)

    bPinf = np.dot(T_inv, np.dot(Pinf, T_inv.T))

    #import pdb; pdb.set_trace()
#    LL,islower = linalg.cho_factor(Pinf)
#    inds = np.triu_indices(Pinf.shape[0],k=1)
#    LL[inds] = 0.0
#    bLL = np.dot(T_inv, LL)
#    bPinf = np.dot( bLL, bLL.T)

    bP0 = np.dot(T_inv, np.dot(P0, T_inv.T))

    if dF is not None:
        bdF = np.zeros(dF.shape)
        for i in range(dF.shape[2]):
            bdF[:,:,i] = np.dot( T_inv, np.dot( dF[:,:,i], T))

    else:
        bdF = None

    if dPinf is not None:
        bdPinf = np.zeros(dPinf.shape)
        for i in range(dPinf.shape[2]):
            bdPinf[:,:,i] = np.dot( T_inv, np.dot( dPinf[:,:,i], T_inv.T))

#            LL,islower = linalg.cho_factor(dPinf[:,:,i])
#            inds = np.triu_indices(dPinf[:,:,i].shape[0],k=1)
#            LL[inds] = 0.0
#            bLL = np.dot(T_inv, LL)
#            bdPinf[:,:,i] = np.dot( bLL, bLL.T)


    else:
        bdPinf = None

    if dP0 is not None:
        bdP0 = np.zeros(dP0.shape)
        for i in range(dP0.shape[2]):
            bdP0[:,:,i] = np.dot( T_inv, np.dot( dP0[:,:,i], T_inv.T))
    else:
        bdP0 = None


    bdQc = dQc # not affected

    # (F,L,Qc,H,Pinf,P0,dF,dQc,dPinf,dP0)

    return bF, bL, bQc, bH, bPinf, bP0, bdF, bdQc, bdPinf, bdP0
646205
def _vmomentsurfacemass(self,R,n,m,romberg=False,nsigma=None,
                           relative=False,phi=0.,deriv=None):
        """"""
        #odd moments of vR are zero
        if isinstance(n,int) and n%2 == 1:
            return 0.
        if nsigma == None:
            nsigma= _NSIGMA
        logSigmaR= self.targetSurfacemass(R,log=True,use_physical=False)
        sigmaR2= self.targetSigma2(R,use_physical=False)
        sigmaR1= sc.sqrt(sigmaR2)
        logsigmaR2= sc.log(sigmaR2)
        if relative:
            norm= 1.
        else:
            norm= sc.exp(logSigmaR+logsigmaR2*(n+m)/2.)/self._gamma**m
        #Use the asymmetric drift equation to estimate va
        va= sigmaR2/2./R**self._beta*(1./self._gamma**2.-1.
                                      -R*self._surfaceSigmaProfile.surfacemassDerivative(R,log=True)
                                      -R*self._surfaceSigmaProfile.sigma2Derivative(R,log=True))
        if math.fabs(va) > sigmaR1: va = 0. #To avoid craziness near the center
        if deriv is None:
            if romberg:
                return sc.real(bovy_dblquad(_vmomentsurfaceIntegrand,
                                            self._gamma*(R**self._beta-va)/sigmaR1-nsigma,
                                            self._gamma*(R**self._beta-va)/sigmaR1+nsigma,
                                            lambda x: -nsigma, lambda x: nsigma,
                                            [R,self,logSigmaR,logsigmaR2,sigmaR1,
                                             self._gamma,n,m],
                                            tol=10.**-8)/sc.pi*norm/2.)
            else:
                return integrate.dblquad(_vmomentsurfaceIntegrand,
                                         self._gamma*(R**self._beta-va)/sigmaR1-nsigma,
                                         self._gamma*(R**self._beta-va)/sigmaR1+nsigma,
                                         lambda x: -nsigma, lambda x: nsigma,
                                         (R,self,logSigmaR,logsigmaR2,sigmaR1,
                                          self._gamma,n,m),
                                         epsrel=_EPSREL)[0]/sc.pi*norm/2.
        else:
            if romberg:
                return sc.real(bovy_dblquad(_vmomentderivsurfaceIntegrand,
                                            self._gamma*(R**self._beta-va)/sigmaR1-nsigma,
                                            self._gamma*(R**self._beta-va)/sigmaR1+nsigma,
                                            lambda x: -nsigma, lambda x: nsigma,
                                            [R,self,logSigmaR,logsigmaR2,sigmaR1,
                                             self._gamma,n,m,deriv],
                                            tol=10.**-8)/sc.pi*norm/2.)
            else:
                return integrate.dblquad(_vmomentderivsurfaceIntegrand,
                                         self._gamma*(R**self._beta-va)/sigmaR1-nsigma,
                                         self._gamma*(R**self._beta-va)/sigmaR1+nsigma,
                                         lambda x: -nsigma, lambda x: nsigma,
                                         (R,self,logSigmaR,logsigmaR2,sigmaR1,
                                          self._gamma,n,m,deriv),
                                         epsrel=_EPSREL)[0]/sc.pi*norm/2.
12313
def generate_gaussian_profile(seeing_fwhm):
    """"""
    FWHM_G = 2 * math.sqrt(2 * math.log(2))
    sigma = seeing_fwhm / FWHM_G
    amplitude = 1.0 / (2 * math.pi * sigma * sigma)
    seeing_model = Gaussian2D(amplitude=amplitude,
                              x_mean=0.0,
                              y_mean=0.0,
                              x_stddev=sigma,
                              y_stddev=sigma)
    return seeing_model
359292
def profile_poly2o(data, mask):
    """"""
    # lmfit
    params = lmfit.Parameters()
    params.add(name="mx", value=0)
    params.add(name="my", value=0)
    params.add(name="mxy", value=0)
    params.add(name="ax", value=0)
    params.add(name="ay", value=0)
    params.add(name="off", value=np.average(data[mask]))
    fr = lmfit.minimize(poly2o_residual, params, args=(data, mask))
    bg = poly2o_model(fr.params, data.shape)
    return bg