'This function evaluates the Gaussian Process model at x.
Parameters
X : array_like
An array with shape (n_eval, n_features) giving the point(s) at
which the prediction(s) should be made.
eval_MSE : boolean, optional
A boolean specifying whether the Mean Squared Error should be
evaluated or not.
Default assumes evalMSE = False and evaluates only the BLUP (mean
prediction).
batch_size : integer, optional
An integer giving the maximum number of points that can be
evaluated simultaneously (depending on the available memory).
Default is None so that all given points are evaluated at the same
time.
Returns
y : array_like, shape (n_samples, ) or (n_samples, n_targets)
An array with shape (n_eval, ) if the Gaussian Process was trained
on an array of shape (n_samples, ) or an array with shape
(n_eval, n_targets) if the Gaussian Process was trained on an array
of shape (n_samples, n_targets) with the Best Linear Unbiased
Prediction at x.
MSE : array_like, optional (if eval_MSE == True)
An array with shape (n_eval, ) or (n_eval, n_targets) as with y,
with the Mean Squared Error at x.'

def predict(self, X, eval_MSE=False, batch_size=None):

    check_is_fitted(self, 'X')
    X = check_array(X)
    (n_eval, _) = X.shape
    (n_samples, n_features) = self.X.shape
    (n_samples_y, n_targets) = self.y.shape
    self._check_params(n_samples)
    if (X.shape[1] != n_features):
        raise ValueError(('The    number    of    features    in    X    (X.shape[1]    =    %d)    should    match    the    number    of    features    used    for    fit()    which    is    %d.' % (X.shape[1], n_features)))
    if (batch_size is None):
        X = ((X - self.X_mean) / self.X_std)
        y = np.zeros(n_eval)
        if eval_MSE:
            MSE = np.zeros(n_eval)
        dx = manhattan_distances(X, Y=self.X, sum_over_features=False)
        f = self.regr(X)
        r = self.corr(self.theta_, dx).reshape(n_eval, n_samples)
        y_ = (np.dot(f, self.beta) + np.dot(r, self.gamma))
        y = (self.y_mean + (self.y_std * y_)).reshape(n_eval, n_targets)
        if (self.y_ndim_ == 1):
            y = y.ravel()
        if eval_MSE:
            C = self.C
            if (C is None):
                if self.verbose:
                    print("This    GaussianProcess    used    'light'    storage    mode    at    instantiation.    Need    to    recompute    autocorrelation    matrix...")
                (reduced_likelihood_function_value, par) = self.reduced_likelihood_function()
                self.C = par['C']
                self.Ft = par['Ft']
                self.G = par['G']
            rt = linalg.solve_triangular(self.C, r.T, lower=True)
            if (self.beta0 is None):
                u = linalg.solve_triangular(self.G.T, (np.dot(self.Ft.T, rt) - f.T), lower=True)
            else:
                u = np.zeros((n_targets, n_eval))
            MSE = np.dot(self.sigma2.reshape(n_targets, 1), ((1.0 - (rt ** 2.0).sum(axis=0)) + (u ** 2.0).sum(axis=0))[np.newaxis, :])
            MSE = np.sqrt(((MSE ** 2.0).sum(axis=0) / n_targets))
            MSE[(MSE < 0.0)] = 0.0
            if (self.y_ndim_ == 1):
                MSE = MSE.ravel()
            return (y, MSE)
        else:
            return y
    else:
        if ((type(batch_size) is not int) or (batch_size <= 0)):
            raise Exception('batch_size    must    be    a    positive    integer')
        if eval_MSE:
            (y, MSE) = (np.zeros(n_eval), np.zeros(n_eval))
            for k in range(max(1, int((n_eval / batch_size)))):
                batch_from = (k * batch_size)
                batch_to = min([(((k + 1) * batch_size) + 1), (n_eval + 1)])
                (y[batch_from:batch_to], MSE[batch_from:batch_to]) = self.predict(X[batch_from:batch_to], eval_MSE=eval_MSE, batch_size=None)
            return (y, MSE)
        else:
            y = np.zeros(n_eval)
            for k in range(max(1, int((n_eval / batch_size)))):
                batch_from = (k * batch_size)
                batch_to = min([(((k + 1) * batch_size) + 1), (n_eval + 1)])
                y[batch_from:batch_to] = self.predict(X[batch_from:batch_to], eval_MSE=eval_MSE, batch_size=None)
            return y
