Model function Gaussian PSF model.
160723
def evaluate(self, x, y, flux, x_0, y_0, sigma):
        """"""

        return (flux / 4 *
                ((self._erf((x - x_0 + 0.5) / (np.sqrt(2) * sigma)) -
                  self._erf((x - x_0 - 0.5) / (np.sqrt(2) * sigma))) *
                 (self._erf((y - y_0 + 0.5) / (np.sqrt(2) * sigma)) -
                  self._erf((y - y_0 - 0.5) / (np.sqrt(2) * sigma)))))
1360060
def _compute_state_estimate(self):
        
        self.x.fill(0)
        for f, mu in zip(self.filters, self.mu):
            self.x += f.x * mu

        self.P.fill(0)
        for f, mu in zip(self.filters, self.mu):
            y = f.x - self.x
            self.P += mu * (outer(y, y) + f.P)
311443
def make_psf_kernel(psf, npix, cdelt, xpix, ypix, psf_scale_fn=None, normalize=False):
    

    egy = psf.energies
    x = make_pixel_distance(npix, xpix, ypix)
    x *= cdelt

    k = np.zeros((len(egy), npix, npix))
    for i in range(len(egy)):
        k[i] = psf.eval(i, x, scale_fn=psf_scale_fn)

    if normalize:
        k /= (np.sum(k, axis=0)[np.newaxis, ...] * np.radians(cdelt) ** 2)

    return k
255673
def spherical_variogram_model(m, d):
    """"""
    psill = float(m[0])
    range_ = float(m[1])
    nugget = float(m[2])
    return np.piecewise(d, [d <= range_, d > range_],
                        [lambda x: psill * ((3.*x)/(2.*range_) - (x**3.)/(2.*range_**3.)) + nugget, psill + nugget])
1180377
def function(self,p):
        """"""
        return 0.5 + 0.5*np.sin(p.frequency*2*pi*self.pattern_y + p.phase)
532389
def fgmc(log_fg_ratios, mu_log_vt, sigma_log_vt, Rf, maxfg):
    '''
    
    '''

    Lb = np.random.uniform(0., maxfg, len(Rf))
    pquit = 0

    while pquit < 0.1:
        # quit when the posterior on Lf is very close to its prior

        nsamp = len(Lb)
        Rf_sel = np.random.choice(Rf, nsamp)
        vt = np.random.lognormal(mu_log_vt, sigma_log_vt, len(Rf_sel))

        Lf = Rf_sel * vt

        log_Lf, log_Lb = log(Lf), log(Lb)

        plR = 0
        for lfr in log_fg_ratios:
            plR += np.logaddexp(lfr + log_Lf, log_Lb)

        plR -= (Lf + Lb)
        plRn = plR - max(plR)

        idx = np.exp(plRn) > np.random.random(len(plRn))

        pquit = ss.stats.ks_2samp(Lb, Lb[idx])[1]

        Lb = Lb[idx]

    return Rf_sel[idx], Lf[idx], Lb
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)
564655
def sparse_GP_regression_2D(num_samples=400, num_inducing=50, max_iters=100, optimize=True, plot=True, nan=False):
    """"""
    np.random.seed(1234)
    X = np.random.uniform(-3., 3., (num_samples, 2))
    Y = np.sin(X[:, 0:1]) * np.sin(X[:, 1:2]) + np.random.randn(num_samples, 1) * 0.05
    if nan:
        inan = np.random.binomial(1,.2,size=Y.shape)
        Y[inan] = np.nan

    # construct kernel
    rbf = GPy.kern.RBF(2)

    # create simple GP Model
    m = GPy.models.SparseGPRegression(X, Y, kernel=rbf, num_inducing=num_inducing)

    # contrain all parameters to be positive (but not inducing inputs)
    m['.*len'] = 2.

    m.checkgrad()

    # optimize
    if optimize:
        m.optimize('tnc', messages=1, max_iters=max_iters)

    # plot
    if plot:
        m.plot()

    print(m)
    return m
880909
function Gaussian(parms) {
      if (parms == null) {
        parms = {};
      }
      if (parms.pi !== void 0) {
        this.pi = parms["pi"];
        this.tau = parms["tau"];
      } else if (parms.mu !== void 0) {
        this.pi = pow(parms["sigma"], -2);
        this.tau = this.pi * parms["mu"];
      } else {
        this.pi = 0;
        this.tau = 0;
      }
      if (isNaN(this.pi) || isNaN(this.tau)) {
        throw new Error("Gaussian parms can not be NaN");
      }
    }
811708
def _create_model(self, X, Y):
        

        # --- define kernel
        self.input_dim = X.shape[1]
        if self.kernel is None:
            kern = GPy.kern.RBF(self.input_dim, variance=1.)
        else:
            kern = self.kernel
            self.kernel = None

        # --- define model
        noise_var = Y.var()*0.01 if self.noise_var is None else self.noise_var
        self.model = GPy.models.GPRegression(X, Y, kernel=kern, noise_var=noise_var)

        # --- Define prior on the hyper-parameters for the kernel (for integrated acquisitions)
        self.model.kern.set_prior(GPy.priors.Gamma.from_EV(2.,4.))
        self.model.likelihood.variance.set_prior(GPy.priors.Gamma.from_EV(2.,4.))

        # --- Restrict variance if exact evaluations of the objective
        if self.exact_feval:
            self.model.Gaussian_noise.constrain_fixed(1e-6, warning=False)
        else:
            self.model.Gaussian_noise.constrain_positive(warning=False)
1166655
def psf_class(self):
        
        if self._psf_type == 'GAUSSIAN':
            psf_type = "GAUSSIAN"
            fwhm = self._seeing
            kwargs_psf = {'psf_type': psf_type, 'fwhm': fwhm}
        elif self._psf_type == 'PIXEL':
            if self._psf_model is not None:
                kwargs_psf = {'psf_type': "PIXEL", 'kernel_point_source': self._psf_model}
            else:
                raise ValueError("You need to create the class instance with a psf_model!")
        else:
            raise ValueError("psf_type %s not supported!" % self._psf_type)
        psf_class = PSF(kwargs_psf)
        return psf_class
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)
120950
def cut_psf(psf_data, psf_size):
    
    kernel = image_util.cut_edges(psf_data, psf_size)
    kernel = kernel_norm(kernel)
    return kernel
213670
def _calc_inst_pmf(self):
        """"""
        # Easy vars
        t = self.t_
        epsilon = self.epsilon
        alpha = self.alpha
        preds = self._preds_avg_in_strata
        weights = self.strata.weights_[:,np.newaxis]
        p1 = self._BB_model.theta_[:,np.newaxis]
        p0 = 1 - p1
        if t==0:
            F = self._F_guess[self.opt_class]
        else:
            F = self._estimate[t - 1, self.opt_class]
            # Fill in non-finite estimates with the initial guess
            nonfinite = ~np.isfinite(F)
            F[nonfinite] = self._F_guess[self.opt_class][nonfinite]

        # Calculate optimal instrumental pmf
        sqrt_arg = np.sum(preds * (alpha**2 * F**2 * p0 + (1 - F)**2 * p1) + \
                          (1 - preds) * (1 - alpha)**2 * F**2 * p1, \
                          axis=1, keepdims=True) #: sum is over classifiers
        inst_pmf = weights * np.sqrt(sqrt_arg)
        # Normalize
        inst_pmf /= np.sum(inst_pmf)
        # Epsilon-greedy: (1 - epsilon) q + epsilon * p
        inst_pmf *= (1 - epsilon)
        inst_pmf += epsilon * weights

        if self.record_inst_hist:
            self._inst_pmf[:,t] = inst_pmf.ravel()
        else:
            self._inst_pmf = inst_pmf.ravel()
383521
def obfn_dfd(self):
        
        """
        XF = sl.rfftn(self.obfn_fvar(), mp_Nv, mp_axisN)
        DX = np.moveaxis(sl.irfftn(sl.inner(mp_Df, XF, mp_axisM),
                                   mp_Nv, mp_axisN), mp_axisM,
                         self.cri.axisM)
        return np.sum((self.W*(DX-self.S))**2)/2.0
470663
def gauss_box_model_deriv(x, amplitude=1.0, mean=0.0, stddev=1.0, hpix=0.5):
    """"""
    z = (x - mean) / stddev
    z2 = z + hpix / stddev
    z1 = z - hpix / stddev

    da = norm.cdf(z2) - norm.cdf(z1)

    fp2 = norm_pdf_t(z2)
    fp1 = norm_pdf_t(z1)

    dl = -amplitude / stddev * (fp2 - fp1)
    ds = -amplitude / stddev * (fp2 * z2 - fp1 * z1)
    dd = amplitude / stddev * (fp2 + fp1)

    return da, dl, ds, dd
389129
def stSpectralRollOff(X, c, fs):
    """"""
    totalEnergy = numpy.sum(X ** 2)
    fftLength = len(X)
    Thres = c*totalEnergy
    # Ffind the spectral rolloff as the frequency position 
    # where the respective spectral energy is equal to c*totalEnergy
    CumSum = numpy.cumsum(X ** 2) + eps
    [a, ] = numpy.nonzero(CumSum > Thres)
    if len(a) > 0:
        mC = numpy.float64(a[0]) / (float(fftLength))
    else:
        mC = 0.0
    return (mC)