def evaluate(self, meanfunc, covfunc, likfunc, x, y, nargout=1):
        if not isinstance(covfunc, cov.FITCOfKernel):
            raise Exception('Only covFITC supported.')  # check cov
        tol = 1e-4; max_sweep = 10; min_sweep = 2       # tolerance to stop EP iterations
        inffunc = EP()

        diagK,Kuu,Ku = covfunc.getCovMatrix(x=x, mode='train')  # evaluate the covariance matrix
        m = meanfunc.getMean(x)                         # evaluate the mean vector

        if likfunc.hyp:                                 # hard coded inducing inputs noise
            sn2  = np.exp(2.*likfunc.hyp[-1])
            snu2 = 1.e-6*sn2                            # similar to infFITC
        else:
            snu2 = 1.e-6

        n, D = x.shape; nu = Kuu.shape[0]
        rot180   = lambda A: np.rot90(np.rot90(A))      # little helper functions
        #chol_inv = lambda A: np.linalg.solve( rot180( np.linalg.cholesky(rot180(A)) ),np.eye(nu)) # chol(inv(A))
        chol_inv = lambda A: rot180( np.linalg.solve(jitchol(rot180(A)), np.eye(nu)) ) # chol(inv(A))

        R0 = chol_inv(Kuu+snu2*np.eye(nu))              # initial R, used for refresh O(nu^3)
        V  = np.dot(R0,Ku); d0 = diagK - np.array([(V*V).sum(axis=0)]).T # initial d, needed for refresh O(n*nu^2)

        # A note on naming: variables are given short but descriptive names in
        # accordance with Rasmussen & Williams "GPs for Machine Learning" (2006): mu
        # and s2 are mean and variance, nu and tau are natural parameters. A leading t
        # means tilde, a subscript _ni means "not i" (for cavity parameters), or _n
        # for a vector of cavity parameters.

        # marginal likelihood for ttau = tnu = zeros(n,1); equals n*log(2) for likCum*
        nlZ0 = -1.* likfunc.evaluate(y, m, np.reshape(diagK,(diagK.shape[0],1)), inffunc).sum()
        if self.last_ttau is None:                      # find starting point for tilde parameters
            ttau  = np.zeros((n,1))                     # initialize to zero if we have no better guess
            tnu   = np.zeros((n,1))
            [d,P,R,nn,gg] = self._epfitcRefresh(d0,Ku,R0,V,ttau,tnu)   # compute initial repres.
            nlZ = nlZ0
        else:
            ttau = self.last_ttau                       # try the tilde values from previous call
            tnu  = self.last_tnu
            [d,P,R,nn,gg] = self._epfitcRefresh(d0,Ku,R0,V,ttau,tnu) # compute initial repres.
            nlZ = self._epfitcZ(d,P,R,nn,gg,ttau,tnu,d0,R0,Ku,y,likfunc,m,inffunc)[0]
            if nlZ > nlZ0:                              # if zero is better ..
                ttau = np.zeros((n,1))                  # .. then initialize with zero instead
                tnu  = np.zeros((n,1))
                [d,P,R,nn,gg] = self._epfitcRefresh(d0,Ku,R0,V,ttau,tnu) # initial repres.
                nlZ = nlZ0

        nlZ_old = np.inf; sweep = 0                     # converged, max. sweeps or min. sweeps?
        while (np.abs(nlZ-nlZ_old) > tol and sweep < max_sweep) or (sweep < min_sweep):
            nlZ_old = nlZ
            sweep += 1
            rperm = list(range(n))                            # randperm(n)
            for ii in rperm:                            # iterate EP updates (in random order) over examples
                p_i = np.reshape(P[:,ii],(P.shape[0],1))
                t = np.dot(R,np.dot(R0,p_i))            # temporary variables
                sigma_i = d[ii] + np.dot(t.T,t); mu_i = nn[ii] + np.dot(p_i.T,gg) # post moments O(nu^2)
                tau_ni = old_div(1,sigma_i) - ttau[ii]                   #  first find the cavity distribution ..
                nu_ni  = old_div(mu_i,sigma_i) + m[ii]*tau_ni - tnu[ii]  # .. params tau_ni and nu_ni
                # compute the desired derivatives of the indivdual log partition function
                vargout = likfunc.evaluate(y[ii], old_div(nu_ni,tau_ni), old_div(1,tau_ni), inffunc, None, 3)
                lZ = vargout[0]; dlZ = vargout[1]; d2lZ = vargout[2]
                ttau_i = old_div(-d2lZ,(1.+old_div(d2lZ,tau_ni)))
                ttau_i = max(ttau_i,0)                  # enforce positivity i.e. lower bound ttau by zero
                tnu_i  = old_div(( dlZ + (m[ii]-old_div(nu_ni,tau_ni))*d2lZ ),(1.+old_div(d2lZ,tau_ni)))
                [d,P[:,ii],R,nn,gg,ttau,tnu] = self._epfitcUpdate(d,P[:,ii],R,nn,gg,ttau,tnu,ii,ttau_i,tnu_i,m,d0,Ku,R0)# update representation

            # recompute since repeated rank-one updates can destroy numerical precision
            [d,P,R,nn,gg] = self._epfitcRefresh(d0,Ku,R0,V,ttau,tnu)
            [nlZ,nu_n,tau_n] = self._epfitcZ(d,P,R,nn,gg,ttau,tnu,d0,R0,Ku,y,likfunc,m,inffunc)
        if sweep == max_sweep:
            pass
            # print '[warning] maximum number of sweeps reached in function infEP'

        self.last_ttau = ttau
        self.last_tnu = tnu       # remember for next call
        post = postStruct()
        post.sW = np.sqrt(ttau)   # unused for FITC_EP prediction with gp.m
        dd = old_div(1,(d0+old_div(1,ttau)))
        alpha = tnu/ttau*dd
        RV = np.dot(R,V)
        R0tV = np.dot(R0.T,V)
        alpha = alpha - np.dot(RV.T,np.dot(RV,alpha))*dd      # long alpha vector for ordinary infEP
        post.alpha = np.dot(R0tV,alpha)                       # alpha = R0'*V*inv(Kt+diag(1./ttau))*(tnu./ttau)
        B = R0tV*np.tile(dd.T,(nu,1)); L = np.dot(B,R0tV.T); B = np.dot(B,RV.T)
        post.L = np.dot(B,B.T) - L                            # L = -R0'*V*inv(Kt+diag(1./ttau))*V'*R0

        if nargout>2:                                         # do we want derivatives?
            dnlZ = dnlZStruct(meanfunc, covfunc, likfunc)     # allocate space for derivatives
            RVdd = RV*np.tile(dd.T,(nu,1))
            for ii in range(len(covfunc.hyp)):
                ddiagK,dKuu,dKu = covfunc.getDerMatrix(x=x, mode='train', der=ii)
                dA = 2*dKu.T - np.dot(R0tV.T,dKuu)            # dQ = dA*R0tV
                w = np.atleast_2d((dA*R0tV.T).sum(axis=1)).T  # w = diag(dQ)
                v = ddiagK - w                                # v = diag(dK)-diag(dQ)
                z = np.dot(dd.T,(v+w)) - np.dot(np.atleast_2d((RVdd*RVdd).sum(axis=0)), v) \
                       - (np.dot(RVdd,dA).T * np.dot(R0tV,RVdd.T)).sum()
                dnlZ.cov[ii] = old_div((z - np.dot(np.dot(alpha.T,dA),np.dot(R0tV,alpha)) - np.dot(alpha.T,(alpha*v))),2.)
                dnlZ.cov[ii] = dnlZ.cov[ii][0,0]
            for ii in range(len(likfunc.hyp)):                # likelihood hypers
                dlik = likfunc.evaluate(y, old_div(nu_n,tau_n)+m, old_div(1,tau_n), inffunc, ii, 1)
                dnlZ.lik[ii] = -dlik.sum()
                if ii == len(likfunc.hyp)-1:
                    # since snu2 is a fixed fraction of sn2, there is a covariance-like term
                    # in the derivative as well
                    v = np.atleast_2d((R0tV*R0tV).sum(axis=0)).T
                    z = (np.dot(RVdd,R0tV.T)**2).sum()  - np.dot(np.atleast_2d((RVdd*RVdd).sum(axis=0)),v)
                    z = z + np.dot(post.alpha.T,post.alpha) - np.dot(alpha.T,(v*alpha))
                    dnlZ.lik[ii] += snu2*z
                    dnlZ.lik[ii] = dnlZ.lik[ii][0,0]
            [junk,dlZ] = likfunc.evaluate(y, old_div(nu_n,tau_n), old_div(1,tau_n), inffunc, None, 2) # mean hyps
            for ii in range(len(meanfunc.hyp)):
                dm = meanfunc.getDerMatrix(x, ii)
                dnlZ.mean[ii] = -np.dot(dlZ.T,dm)
                dnlZ.mean[ii] = dnlZ.mean[ii][0,0]

            return post, nlZ[0,0], dnlZ
        else:
            return post, nlZ[0,0]
dataset/ETHPy150Open marionmari/pyGPs/pyGPs/Core/inf.py/FITC_EP.evaluate/swapped operands of `-`