def evaluate(self, meanfunc, covfunc, likfunc, x, y, nargout=1):
        tol = 1e-6                           # tolerance for when to stop the Newton iterations
        smax = 2; Nline = 20; thr = 1e-4     # line search parameters
        maxit = 20                           # max number of Newton steps in f
        inffunc = self
        K = covfunc.getCovMatrix(x=x, mode='train')    # evaluate the covariance matrix
        m = meanfunc.getMean(x)              # evaluate the mean vector
        n, D = x.shape
        Psi_old = np.inf                     # make sure while loop starts by the largest old objective val
        if self.last_alpha is None:          # find a good starting point for alpha and f
            alpha = np.zeros((n,1))
            f = np.dot(K,alpha) + m          # start at mean if sizes not match
            vargout = likfunc.evaluate(y, f, None, inffunc, None, 3)
            lp = vargout[0]; dlp = vargout[1]; d2lp = vargout[2]
            W= -d2lp
            Psi_new = -lp.sum()
        else:
            alpha = self.last_alpha
            f = np.dot(K,alpha) + m                       # try last one
            vargout = likfunc.evaluate(y, f, None, inffunc, None, 3)
            lp = vargout[0]; dlp = vargout[1]; d2lp = vargout[2]
            W= -d2lp
            Psi_new = lp.sum() - old_div(np.dot(alpha.T,(f-m)),2.) # objective for last alpha
            vargout = - likfunc.evaluate(y, m, None, inffunc, None, 1)
            Psi_def =  vargout[0]                         # objective for default init f==m
            if Psi_def < Psi_new:                         # if default is better, we use it
                alpha = np.zeros((n,1))
                f = np.dot(K,alpha) + m
                vargout = likfunc.evaluate(y, f, None, inffunc, None, 3)
                lp = vargout[0]; dlp = vargout[1]; d2lp = vargout[2]
                W=-d2lp; Psi_new = -lp.sum()
        isWneg = np.any(W<0)       # flag indicating whether we found negative values of W
        it = 0                     # this happens for the Student's t likelihood

        while (Psi_old - Psi_new > tol) and it<maxit:      # begin Newton
            Psi_old = Psi_new; it += 1
            if isWneg:                   # stabilise the Newton direction in case W has negative values
                W = np.maximum(W,0)      # stabilise the Hessian to guarantee postive definiteness
                tol = 1e-10              # increase accuracy to also get the derivatives right
            #sW = np.sqrt(W); L = np.linalg.cholesky(np.eye(n) + np.dot(sW,sW.T)*K).T
            sW = np.sqrt(W); L = jitchol(np.eye(n) + np.dot(sW,sW.T)*K).T
            b = W*(f-m) + dlp;
            dalpha = b - sW*solve_chol(L,sW*np.dot(K,b)) - alpha
            vargout = brentmin(0,smax,Nline,thr,self._Psi_line,4,dalpha,alpha,K,m,likfunc,y,inffunc)
            s = vargout[0]
            Psi_new = vargout[1]
            Nfun = vargout[2]
            alpha = vargout[3]
            f = vargout[4]
            dlp = vargout[5]
            W = vargout[6]
            isWneg = np.any(W<0)
        self.last_alpha = alpha                 # remember for next call
        vargout = likfunc.evaluate(y,f,None,inffunc,None,4)
        lp = vargout[0]; dlp = vargout[1]; d2lp = vargout[2]; d3lp = vargout[3]
        W = -d2lp; isWneg = np.any(W<0)
        post = postStruct()
        post.alpha = alpha                      # return the posterior parameters
        post.sW = np.sqrt(np.abs(W))*np.sign(W) # preserve sign in case of negative
        if isWneg:
            [ldA,iA,post.L] = self._logdetA(K,W,3)
            nlZ = old_div(np.dot(alpha.T,(f-m)),2.) - lp.sum() + old_div(ldA,2.)
            nlZ = nlZ[0]
        else:
            sW = post.sW
            #post.L = np.linalg.cholesky(np.eye(n)+np.dot(sW,sW.T)*K).T
            post.L = jitchol(np.eye(n)+np.dot(sW,sW.T)*K).T
            nlZ = old_div(np.dot(alpha.T,(f-m)),2.) + (np.log(np.diag(post.L))-np.reshape(lp,(lp.shape[0],))).sum()
            nlZ = nlZ[0]
        if nargout>2:                                           # do we want derivatives?
            dnlZ = dnlZStruct(meanfunc, covfunc, likfunc)       # allocate space for derivatives
            if isWneg:                                          # switch between Cholesky and LU decomposition mode
                Z = -post.L                                     # inv(K+inv(W))
                g = old_div(np.atleast_2d((iA*K).sum(axis=1)).T,2)      # deriv. of ln|B| wrt W; g = diag(inv(inv(K)+diag(W)))/2
            else:
                Z = np.tile(sW,(1,n))*solve_chol(post.L,np.diag(np.reshape(sW,(sW.shape[0],)))) #sW*inv(B)*sW=inv(K+inv(W))
                C = np.linalg.solve(post.L.T,np.tile(sW,(1,n))*K)            # deriv. of ln|B| wrt W
                g = old_div(np.atleast_2d((np.diag(K)-(C**2).sum(axis=0).T)).T,2.)   # g = diag(inv(inv(K)+W))/2
            dfhat = g* d3lp                 # deriv. of nlZ wrt. fhat: dfhat=diag(inv(inv(K)+W)).*d3lp/2
            for ii in range(len(covfunc.hyp)):                               # covariance hypers
                dK = covfunc.getDerMatrix(x=x, mode='train', der=ii)
                dnlZ.cov[ii] = old_div((Z*dK).sum(),2.) - old_div(np.dot(alpha.T,np.dot(dK,alpha)),2.)   # explicit part
                b = np.dot(dK,dlp)                                           # b-K*(Z*b) = inv(eye(n)+K*diag(W))*b
                dnlZ.cov[ii] -= np.dot(dfhat.T,b-np.dot(K,np.dot(Z,b)))      # implicit part
                dnlZ.cov[ii] = dnlZ.cov[ii][0,0]
            for ii in range(len(likfunc.hyp)):                               # likelihood hypers
                [lp_dhyp,dlp_dhyp,d2lp_dhyp] = likfunc.evaluate(y,f,None,inffunc,ii,3)
                dnlZ.lik[ii] = -np.dot(g.T,d2lp_dhyp) - lp_dhyp.sum()        # explicit part
                b = np.dot(K,dlp_dhyp)                                       # b-K*(Z*b) = inv(eye(n)+K*diag(W))*b
                dnlZ.lik[ii] -= np.dot(dfhat.T,b-np.dot(K,np.dot(Z,b)))      # implicit part
                dnlZ.lik[ii] = dnlZ.lik[ii][0,0]
            for ii in range(len(meanfunc.hyp)):                              # mean hypers
                dm = meanfunc.getDerMatrix(x, ii)
                dnlZ.mean[ii] = -np.dot(alpha.T,dm)                          # explicit part
                dnlZ.mean[ii] -= np.dot(dfhat.T,dm-np.dot(K,np.dot(Z,dm)))   # implicit part
                dnlZ.mean[ii] = dnlZ.mean[ii][0,0]
            return post,nlZ[0],dnlZ
        else:
            return post, nlZ[0]
dataset/ETHPy150Open marionmari/pyGPs/pyGPs/Core/inf.py/Laplace.evaluate/swapped operands of `-`