def gp_posterior_moment_function(m, k, x, y, k_sparse=None, pseudoinputs=None, noise=None):

    # Prior
    # FIXME: We are ignoring the covariance of mu now..
    mu = m(x)[0]
    ## if np.ndim(mu) == 1:
    ##     mu = np.asmatrix(mu).T
    ## else:
    ##     mu = np.asmatrix(mu)
    
    K_noise = None
    
    if noise != None:
        if K_noise is None:
            K_noise = noise
        else:
            K_noise += noise
            
    if k_sparse != None:
        if K_noise is None:
            K_noise = k_sparse(x,x)[0]
        else:
            K_noise += k_sparse(x,x)[0]

    if pseudoinputs != None:
        p = pseudoinputs
        #print('in pseudostuff')
        #print(K_noise)
        #print(np.shape(K_noise))
        K_pp = k(p,p)[0]
        K_xp = k(x,p)[0]
        U = utils.chol(K_noise)

        # Compute Lambda
        Lambda = K_pp + np.dot(K_xp.T, utils.chol_solve(U, K_xp))
        U_lambda = utils.chol(Lambda)

        # Compute statistics for posterior predictions
        #print(np.shape(U_lambda))
        #print(np.shape(y))
        z = utils.chol_solve(U_lambda,
                       np.dot(K_xp.T,
                              utils.chol_solve(U,
                                         y - mu)))
        U = utils.chol(K_pp)

        # Now we can forget the location of the observations and
        # consider only the pseudoinputs when predicting.
        x = p

        
    else:
        K = K_noise
        if K is None:
            K = k(x,x)[0]
        else:
            try:
                K += k(x,x)[0]
            except:
                K = K + k(x,x)[0]

        # Compute posterior GP
        N = len(y)
        U = None
        z = None
        if N > 0:
            U = utils.chol(K)
            z = utils.chol_solve(U, y-mu)

    def get_moments(h, covariance=1, mean=True):

        K_xh = k(x, h)[0]
        if k_sparse != None:
            try:
                # This may not work, for instance, if either one is a
                # sparse matrix.
                K_xh += k_sparse(x, h)[0]
            except:
                K_xh = K_xh + k_sparse(x, h)[0]
        
        # NumPy has problems when mixing matrices and arrays.
        # Matrices may appear, for instance, when you sum an array and
        # a sparse matrix.  Make sure the result is either an array or
        # a sparse matrix (not dense matrix!), because matrix objects
        # cause lots of problems:
        #
        # array.dot(array) = array
        # matrix.dot(array) = matrix
        # sparse.dot(array) = array
        if not sp.issparse(K_xh):
            K_xh = np.asarray(K_xh)

        # Function for computing posterior moments
        if mean:
            # Mean vector
            # FIXME: Ignoring the covariance of prior mu
            m_h = m(h)[0]
            
            if z != None:
                m_h += K_xh.T.dot(z)
                
        else:
            m_h = None

        # Compute (co)variance matrix/vector
        if covariance:
            if covariance == 1:
                ## Compute variance vector
                
                k_h = k(h)[0]
                if k_sparse != None:
                    k_h += k_sparse(h)[0]
                if U != None:
                    if isinstance(K_xh, np.ndarray):
                        k_h -= np.einsum('i...,i...',
                                         K_xh,
                                         utils.chol_solve(U, K_xh))
                    else:
                        # TODO: This isn't very efficient way, but
                        # einsum doesn't work for sparse matrices..
                        # This may consume A LOT of memory for sparse
                        # matrices.
                        k_h -= np.asarray(K_xh.multiply(utils.chol_solve(U, K_xh))).sum(axis=0)
                if pseudoinputs != None:
                    if isinstance(K_xh, np.ndarray):
                        k_h += np.einsum('i...,i...',
                                         K_xh,
                                         utils.chol_solve(U_lambda, K_xh))
                    else:
                        # TODO: This isn't very efficient way, but
                        # einsum doesn't work for sparse matrices..
                        # This may consume A LOT of memory for sparse
                        # matrices.
                        k_h += np.asarray(K_xh.multiply(utils.chol_solve(U_lambda, K_xh))).sum(axis=0)
                # Ensure non-negative variances        
                k_h[0 < k_h] = 0
                
                return (m_h, k_h)
                    
            elif covariance == 2:
                ## Compute full covariance matrix
                
                K_hh = k(h,h)[0]
                if k_sparse != None:
                    K_hh += k_sparse(h)[0]
                if U != None:
                    K_hh -= K_xh.T.dot(utils.chol_solve(U,K_xh))
                    #K_hh -= np.dot(K_xh.T, utils.chol_solve(U,K_xh))
                if pseudoinputs != None:
                    K_hh += K_xh.T.dot(utils.chol_solve(U_lambda, K_xh))
                    #K_hh += np.dot(K_xh.T, utils.chol_solve(U_lambda, K_xh))
                return (m_h, K_hh)
        else:
            return (m_h, None)


    return get_moments


# Constant function using GP mean protocol
dataset/ETHPy150Open bayespy/bayespy/bayespy/inference/vmp/nodes/GaussianProcesses.py/gp_posterior_moment_function/swapped operands of `<`