def predict(self, X, eval_MSE=True):

        X = array2d(X)
        nlevel = self.nlevel
        n_eval, n_features_X = X.shape
        
        # Calculate kriging mean and variance at level 0
        mu = np.zeros((n_eval, nlevel))

        f = self.regr(X)
        f0 = self.regr(X)
        dx = l1_cross_distances(X, Y=self.X[0])
        
        # Get regression function and correlation
        F = self.F[0]
        C = self.C[0]
        
        beta = self.beta[0]
        Ft = solve_triangular(C, F, lower=True)
        yt = solve_triangular(C, self.y[0], lower=True)
        r_ = self.corr(self.theta[0], dx).reshape(n_eval, self.n_samples[0])
        gamma = solve_triangular(C.T, yt - np.dot(Ft,beta), lower=False)

        # Scaled predictor
        mu[:,0]= (np.dot(f, beta) + np.dot(r_,gamma)).ravel()

        if eval_MSE:
            self.sigma2_rho = nlevel*[None]
            MSE = np.zeros((n_eval,nlevel))
            r_t = solve_triangular(C, r_.T, lower=True)
            G = self.G[0]
            
            u_ = solve_triangular(G.T, f.T - np.dot(Ft.T, r_t), lower=True)
            MSE[:,0] = self.sigma2[0] * (1 \
                        - (r_t**2).sum(axis=0) + (u_**2).sum(axis=0))
    
        # Calculate recursively kriging mean and variance at level i
        for i in range(1,nlevel):
            C = self.C[i]
            F = self.F[i]
            g = self.rho_regr(X)
            dx = l1_cross_distances(X, Y=self.X[i])
            r_ = self.corr(self.theta[i], dx).reshape(n_eval, self.n_samples[i])   
            f = np.vstack((g.T*mu[:,i-1], f0.T))

            Ft = solve_triangular(C, F, lower=True)
            yt = solve_triangular(C, self.y[i], lower=True)
            r_t = solve_triangular(C,r_.T, lower=True)
            G = self.G[i]
            beta = self.beta[i]

            # scaled predictor
            mu[:,i] = (np.dot(f.T, beta) \
                       + np.dot(r_t.T, yt - np.dot(Ft,beta))).ravel()
            
            if eval_MSE:
                Q_ = (np.dot((yt-np.dot(Ft,beta)).T, yt-np.dot(Ft,beta)))[0,0]
                u_ = solve_triangular(G.T, f - np.dot(Ft.T, r_t), lower=True)
                sigma2_rho = np.dot(g, \
                    self.sigma2[i]*linalg.inv(np.dot(G.T,G))[:self.q[i],:self.q[i]] \
                        + np.dot(beta[:self.q[i]], beta[:self.q[i]].T))
                sigma2_rho = (sigma2_rho * g).sum(axis=1)
                
                MSE[:,i] = sigma2_rho * MSE[:,i-1] \
                        + Q_/(2*(self.n_samples[i]-self.p[i]-self.q[i])) \
                        * (1 - (r_t**2).sum(axis=0)) \
                        + self.sigma2[i] * (u_**2).sum(axis=0)                       
        
        # scaled predictor                  
        for i in range(nlevel):# Predictor
            mu[:,i] = self.y_mean + self.y_std * mu[:,i]
            if eval_MSE:
                MSE[:,i] = self.y_std**2 * MSE[:,i]
            
        if eval_MSE:
            return mu[:,-1].reshape((n_eval,1)), MSE[:,-1].reshape((n_eval,1))
        else:
            return mu[:,-1].reshape((n_eval,1))
This function performs the predictions of the kriging model on X.
dataset/ETHPy150Open OpenMDAO/OpenMDAO-Framework/openmdao.lib/src/openmdao/lib/surrogatemodels/multifi_cokriging_surrogate.py/MultiFiCoKriging.predict