Integrated GP-Lower Confidence Bound and its derivative
330630
def _compute_acq_withGradients(self, x):
        
        means, stds, dmdxs, dsdxs = self.model.predict_withGradients(x)
        f_acqu = None
        df_acqu = None
        for m, s, dmdx, dsdx in zip(means, stds, dmdxs, dsdxs):  
            f = -m + self.exploration_weight * s        
            df = -dmdx + self.exploration_weight * dsdx
            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))
469883
def leastsq_generalized(x, y, x2, y2, res_std=None, res2_std=None, fit_offset=False, fit_offset2=False, perc=None):
    
    if perc is not None:
        if not fit_offset and isinstance(perc, (list, tuple)): perc = perc[1]
        weights = csr_matrix(get_weight(x, y, perc)).astype(bool)
        x, y = weights.multiply(x).tocsr(), weights.multiply(y).tocsr()

    n_obs, n_var = x.shape
    offset, offset_ss = np.zeros(n_var, dtype="float32"), np.zeros(n_var, dtype="float32")
    gamma = np.ones(n_var, dtype="float32")

    if (res_std is None) or (res2_std is None): res_std, res2_std = np.ones(n_var), np.ones(n_var)
    ones, zeros = np.ones(n_obs), np.zeros(n_obs)

    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        x, y = np.vstack((make_dense(x)/res_std, x2/res2_std)), np.vstack((make_dense(y)/res_std, y2/res2_std))

    if fit_offset and fit_offset2:
        for i in range(n_var):
            A = np.c_[np.vstack((np.c_[ones/res_std[i], zeros], np.c_[zeros, ones/res2_std[i]])), x[:, i]]
            offset[i], offset_ss[i], gamma[i] = np.linalg.pinv(A.T.dot(A)).dot(A.T.dot(y[:, i]))
    elif fit_offset:
        for i in range(n_var):
            A = np.c_[np.hstack((ones/res_std[i], zeros)), x[:, i]]
            offset[i], gamma[i] = np.linalg.pinv(A.T.dot(A)).dot(A.T.dot(y[:, i]))
    elif fit_offset2:
        for i in range(n_var):
            A = np.c_[np.hstack((zeros, ones/res2_std[i])), x[:, i]]
            offset_ss[i], gamma[i] = np.linalg.pinv(A.T.dot(A)).dot(A.T.dot(y[:, i]))
    else:
        for i in range(n_var):
            A = np.c_[x[:, i]]
            gamma[i] = np.linalg.pinv(A.T.dot(A)).dot(A.T.dot(y[:, i]))

    offset[np.isnan(offset)] = 0
    offset_ss[np.isnan(offset_ss)] = 0
    gamma[np.isnan(gamma)] = 0

    return offset, offset_ss, gamma
382572
function ecp_proj(R, N, p, l) {
	var x = N[0], y = N[1], z = N[2];
	var xr = R[0], yr = R[1], zr = R[2];
	var t0 = _ecp_t0, t1 = _ecp_t1, t2 = _ecp_t2;
	fp_mul(t2, t0, fp_sqr(t1, fp_inv(t0, z, p, l), p, l), p, l);
	fp_mul(xr, x, t1, p, l);
	fp_mul(yr, y, t2, p, l);
	mpn_zero(zr, l), zr[0] = 1;
	return R;
}
337141
def Kdiag(self, X, target):
        """"""
        ly=1/self.lengthscaleY
        lu=np.sqrt(3)/self.lengthscaleU
        #ly=self.lengthscaleY
        #lu=self.lengthscaleU
        
        k1 = (2*lu+ly)/(lu+ly)**2
        k2 = (ly-2*lu + 2*lu-ly ) / (ly-lu)**2 
        k3 = 1/(lu+ly) + (lu)/(lu+ly)**2 

        np.add(self.varianceU*self.varianceY*(k1+k2+k3), target, target)
324454
def black_box_variational_inference(logprob, D, num_samples):
    """"""

    def unpack_params(params):
        # Variational dist is a diagonal Gaussian.
        mean, log_std = params[:D], params[D:]
        return mean, log_std

    def gaussian_entropy(log_std):
        return 0.5 * D * (1.0 + np.log(2*np.pi)) + np.sum(log_std)

    rs = npr.RandomState(0)
    def variational_objective(params, t):
        """Provides a stochastic estimate of the variational lower bound."""
        mean, log_std = unpack_params(params)
        samples = rs.randn(num_samples, D) * np.exp(log_std) + mean
        lower_bound = gaussian_entropy(log_std) + np.mean(logprob(samples, t))
        return -lower_bound

    gradient = grad(variational_objective)

    return variational_objective, gradient, unpack_params
1082636
def get_variable_errors(X, L, LX, Z, U, step_g, e_rel, e_abs=0):
    
    n = X.size
    p = Z.size
    e_pri2 = np.sqrt(p)*e_abs/L.spectral_norm + e_rel*np.max([l2(LX), l2(Z)])
    if step_g is not None:
        e_dual2 = np.sqrt(n)*e_abs/L.spectral_norm + e_rel*l2(L.T.dot(U)/step_g)
    else:
        e_dual2 = np.sqrt(n)*e_abs/L.spectral_norm + e_rel*l2(L.T.dot(U))
    return e_pri2, e_dual2
602121
def _compute_mle(self):
        
        min_y = np.min(self._interp.y)
        if self._interp.y[0] == min_y:
            self._mle = self._interp.x[0]
        elif self._interp.y[-1] == min_y:
            self._mle = self._interp.x[-1]
        else:
            argmin_y = np.argmin(self._interp.y)
            ix0 = max(argmin_y - 4, 0)
            ix1 = min(argmin_y + 4, len(self._interp.x) - 1)

            while np.sign(self._interp.derivative(self._interp.x[ix0])) == \
                    np.sign(self._interp.derivative(self._interp.x[ix1])):
                ix0 += 1

            self._mle = scipy.optimize.brentq(self._interp.derivative,
                                              self._interp.x[ix0],
                                              self._interp.x[ix1],
                                              xtol=1e-10 *
                                              np.median(self._interp.x))
463897
def gradF_t(X, Y, S, M_E, E, m0, rho):
    ''' 
    '''
    n, r = X.shape
    m, r = Y.shape

    XS = np.dot(X, S)
    YS = np.dot(Y, S.T)
    XSY = np.dot(XS, Y.T)

    Qx = np.dot(np.dot(X.T, ((M_E - XSY) * E)), YS) / n
    Qy = np.dot(np.dot(Y.T, ((M_E - XSY) * E).T), XS) / m

    W = np.dot((XSY - M_E) * E, YS) + np.dot(X, Qx) + rho * Gp(X, m0, r)
    Z = np.dot(((XSY - M_E) * E).T, XS) + np.dot(Y, Qy) + rho * Gp(Y, m0, r)

    return W, Z
990102
def root_mean_squared_error(pred:Tensor, targ:Tensor)->Rank0Tensor:
    ""
    pred,targ = flatten_check(pred,targ)
    return torch.sqrt(F.mse_loss(pred, targ))
1368338
def NegLnLike(x, time, flux, errors, kernel):
    '''
    

    '''

    gp = GP(kernel, x, white=True)
    gp.compute(time, errors)
    if OLDGEORGE:
        nll = -gp.lnlikelihood(flux)
        # NOTE: There was a bug on this next line! Used to be
        #
        #    ngr = -gp.grad_lnlikelihood(flux) / gp.kernel.pars
        #
        # But I think we want
        #
        # dlogL/dx =     dlogL/dlogx^2       * dlogx^2/dx^2 * dx^2/dx
        #          = gp.grad_lnlikelihood()  *     1/x^2    *   2x
        #          = 2 * gp.grad_lnlikelihood() / x
        #          = 2 * gp.grad_lnlikelihood() / np.sqrt(x^2)
        #          = 2 * gp.grad_lnlikelihood() / np.sqrt(gp.kernel.pars)
        #
        # (with a negative sign out front for the negative gradient).
        # So we probably weren't optimizing the GP correctly! This affects
        # all campaigns through C13. It's not a *huge* deal, since the sign
        # of the gradient was correct and the model isn't that sensitive to
        # the value of the hyperparameters, but it may have contributed to
        # the poor performance on super variable stars. In most cases it means
        # the solver takes longer to converge and isn't as good at finding
        # the minimum.
        ngr = -2 * gp.grad_lnlikelihood(flux) / np.sqrt(gp.kernel.pars)
    else:
        nll = -gp.log_likelihood(flux)
        ngr = -2 * gp.grad_log_likelihood(flux) / \
            np.sqrt(np.exp(gp.get_parameter_vector()))

    return nll, ngr
720558
def optimise_gradient_descent(x, a, c, tolerance=0.001):
    
    x_new = x
    x_old = x-1 # dummy value
    while np.abs(x_new - x_old).sum() > tolerance:
        x_old = x_new
        step_size = golden_section_search(lambda step: f(grad_desc_update(x_old, a, c, step), a, c), -1.0, 1.0)
        x_new = grad_desc_update(x_old, a, c, step_size)
    return x_new
426018
def inverse_prediction(self, s_g):
        
        return self.infer(self.conf.s_dims, self.conf.m_dims, s_g)
1247906
def xstep(self):
        
        """

        ngsit = 0
        gsrrs = np.inf
        while gsrrs > self.opt['GSTol'] and ngsit < self.opt['MaxGSIter']:
            self.X = self.GaussSeidelStep(self.S, self.X,
                                          self.cnst_AT(self.Y-self.U),
                                          self.rho, self.lcw, self.Wdf2)
            gsrrs = sl.rrs(
                self.rho*self.cnst_AT(self.cnst_A(self.X)) +
                self.Wdf2*self.X, self.Wdf2*self.S +
                self.rho*self.cnst_AT(self.Y - self.U))
            ngsit += 1

        self.xs = (ngsit, gsrrs)
314745
def fit(self,Params0=None,grad_threshold=1e-2):
        

        # GPVD initialization
        lik = limix.CLikNormalNULL()
        # Initial Params
        if Params0==None:
            n_params = self.C1.getNumberParams()
            n_params+= self.C2.getNumberParams()
            Params0 = SP.rand(n_params)
        # MultiGP Framework
        covar = []
        gp    = []
        mean  = []
        for i in range(self.N):
            covar.append(limix.CLinCombCF())
            covar[i].addCovariance(self.C1)
            covar[i].addCovariance(self.C2)
            coeff = SP.array([self.eigen[i],1])
            covar[i].setCoeff(coeff)
            mean.append(limix.CLinearMean(self.Yt[:,i],SP.eye(self.P)))
        gpVD = limix.CGPvarDecomp(covar[0],lik,mean[0],SP.ones(self.N),self.P,self.Yt,Params0)
        for i in range(self.N):
            gp.append(limix.CGPbase(covar[i],lik,mean[i]))
            gp[i].setY(self.Yt[:,i])
            gpVD.addGP(gp[i])
                
        # Optimization
        gpVD.initGPs()
        gpopt = limix.CGPopt(gpVD)
        LML0=-1.0*gpVD.LML()
        start_time = time.time()
        conv = gpopt.opt()
        time_train = time.time() - start_time
        LML=-1.0*gpVD.LML()
        LMLgrad = SP.linalg.norm(gpVD.LMLgrad()['covar'])
        Params = gpVD.getParams()['covar']
    
        # Check whether limix::CVarianceDecomposition.train() has converged
        if conv!=True or LMLgrad>grad_threshold or Params.max()>10:
            print('limix::CVarianceDecomposition::train has not converged')
            res=None
        else:
            res = {
                'Params0':          Params0,
                'Params':           Params,
                'LML':              SP.array([LML]),
                'LML0':             SP.array([LML0]),
                'LMLgrad':          SP.array([LMLgrad]),
                'time_train':       SP.array([time_train]),
                }
        return res
        pass
72154
def confidence_interval_continuous(
        point_estimate, stddev, sample_size, confidence=.95, **kwargs
):
    """"""
    alpha = ppf((confidence + 1) / 2, sample_size - 1)

    margin = stddev / sqrt(sample_size)
    return (point_estimate - alpha * margin, point_estimate + alpha * margin)
325343
def confidence_interval_min_width(mu, post, alpha=0.9):
    '''
    
    '''
    if not 0 < alpha < 1:
        raise ValueError("Confidence level must be in (0,1).")

    # choose a step size for the sliding confidence window
    alpha_step = 0.01

    # initialize the lower and upper limits
    mu_low = numpy.min(mu)
    mu_high = numpy.max(mu)

    # find the smallest window (by delta-mu) stepping by dalpha
    for ai in numpy.arange(0, 1 - alpha, alpha_step):
        ml = compute_lower_limit(mu, post, 1 - ai)
        mh = compute_upper_limit(mu, post, alpha + ai)
        if mh - ml < mu_high - mu_low:
            mu_low = ml
            mu_high = mh

    return mu_low, mu_high
1283546
def trainGP(self,fast=False,scales0=None,fixed0=None,lambd=None):
        
        assert self.n_terms>0, 'CVarianceDecomposition:: No variance component terms'

        if not self.init:		self.initGP(fast=fast)

        # set lambda
        if lambd!=None:		self.gp.setLambda(lambd)

        # set scales0
        if scales0!=None:
            self.setScales(scales0)
        # init gp params
        self.vd.initGPparams()
        # set fixed0
        if fixed0!=None:
            params = self.gp.getParams()
            params['dataTerm'] = fixed0
            self.gp.setParams(params)

        # LIMIX CVARIANCEDECOMPOSITION TRAINING
        conv =self.vd.trainGP()
        
        self.cache['Sigma']   = None
        self.cache['Hessian'] = None
            
        return conv