def evaluate(self, meanfunc, covfunc, likfunc, x, y, nargout=1):
        tol = 1e-4; max_sweep = 10; min_sweep = 2 # tolerance to stop EP iterations
        n = x.shape[0]
        inffunc = self
        K = covfunc.getCovMatrix(x=x, mode='train') # evaluate the covariance matrix
        m = meanfunc.getMean(x)                   # evaluate the mean vector
        nlZ0 = -likfunc.evaluate(y, m, np.reshape(np.diag(K),(np.diag(K).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))
            Sigma = K                             # initialize Sigma and mu, the parameters of ..
            mu    = np.zeros((n,1))               # .. the Gaussian posterior approximation
            nlZ   = nlZ0
        else:
            ttau = self.last_ttau                 # try the tilde values from previous call
            tnu  = self.last_tnu
            Sigma, mu, nlZ, L = self._epComputeParams(K, y, ttau, tnu, likfunc, m, inffunc)
            if nlZ > nlZ0:                        # if zero is better ..
                ttau = np.zeros((n,1))            # .. then initialize with zero instead
                tnu  = np.zeros((n,1))
                Sigma = K                         # initialize Sigma and mu, the parameters of ..
                mu = np.zeros((n,1))              # .. the Gaussian posterior approximation
                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 = range(n)                     # randperm(n)
            for ii in rperm:                      # iterate EP updates (in random order) over examples
                tau_ni = old_div(1,Sigma[ii,ii]) - ttau[ii]#  first find the cavity distribution ..
                nu_ni  = old_div(mu[ii],Sigma[ii,ii]) + m[ii]*tau_ni - tnu[ii]    # .. params tau_ni and nu_ni
                # compute the desired derivatives of the indivdual log partition function
                lZ,dlZ,d2lZ = likfunc.evaluate(y[ii], old_div(nu_ni,tau_ni), old_div(1,tau_ni), inffunc, None, 3)
                ttau_old = copy(ttau[ii])         # then find the new tilde parameters, keep copy of old
                ttau[ii] = old_div(-d2lZ,(1.+old_div(d2lZ,tau_ni)))
                ttau[ii] = max(ttau[ii],0)        # enforce positivity i.e. lower bound ttau by zero
                tnu[ii]  = old_div(( dlZ + (m[ii]-old_div(nu_ni,tau_ni))*d2lZ ),(1.+old_div(d2lZ,tau_ni)))
                ds2 = ttau[ii] - ttau_old         # finally rank-1 update Sigma ..
                si  = np.reshape(Sigma[:,ii],(Sigma.shape[0],1))
                Sigma = Sigma - (1.+ds2*si[ii])*np.dot(si,si.T)/ds2   # takes 70# of total time
                mu = np.dot(Sigma,tnu)                                # .. and recompute mu
            # recompute since repeated rank-one updates can destroy numerical precision
            Sigma, mu, nlZ, L = self._epComputeParams(K, y, ttau, tnu, 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
        sW = np.sqrt(ttau); alpha = tnu-sW*solve_chol(L,sW*np.dot(K,tnu))
        post = postStruct()
        post.alpha = alpha                                  # return the posterior params
        post.sW    = sW
        post.L     = L
        if nargout>2:                                       # do we want derivatives?
            dnlZ = dnlZStruct(meanfunc, covfunc, likfunc)   # allocate space for derivatives
            ssi  = np.sqrt(ttau)
            V = np.linalg.solve(L.T,np.tile(ssi,(1,n))*K)
            Sigma = K - np.dot(V.T,V)
            mu = np.dot(Sigma,tnu)
            Dsigma = np.reshape(np.diag(Sigma),(np.diag(Sigma).shape[0],1))
            tau_n = old_div(1,Dsigma)-ttau                           # compute the log marginal likelihood
            nu_n  = old_div(mu,Dsigma)-tnu                           # vectors of cavity parameters
            F = np.dot(alpha,alpha.T) - np.tile(sW,(1,n))* \
                solve_chol(L,np.diag(np.reshape(sW,(sW.shape[0],))))   # covariance hypers
            for jj in range(len(covfunc.hyp)):
                dK = covfunc.getDerMatrix(x=x, mode='train', der=jj)
                dnlZ.cov[jj] = old_div(-(F*dK).sum(),2.)
            for ii in range(len(likfunc.hyp)):
                dlik = likfunc.evaluate(y, old_div(nu_n,tau_n), old_div(1,tau_n), inffunc, ii)
                dnlZ.lik[ii] = -dlik.sum()
            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], dnlZ
        else:
            return post, nlZ[0]
dataset/ETHPy150Open marionmari/pyGPs/pyGPs/Core/inf.py/EP.evaluate/swapped operands of `/`