function y_pred = get_ypred_prs(x,y,x_test,nvar,npoly)
% This function builds an PRS metamodel and predicts the function value at
% test points.

% ---Variable Descriptions---
% x         = normalized design(training points): in an 
%           [num_of_points by num_of_variables] matrix
% y         = vector of function values at each x
% x_test    = test points: in an 
%           [num_of_points by num_of_variables] matrix
% nvar      = the number of design variables
% npoly     = user chosen polynomial degree

[X, nbeta1] = create_X_matrix(x, nvar, npoly);
b = X\y;   % solve for the coefficients

% now make predictions
[Xt, nbeta1] = create_X_matrix(x_test, nvar, npoly);
y_pred = Xt*b;

% Create X matrix
function [X, nbeta1] = create_X_matrix(vrn, nvar, npoly)
n_term = 1;
n_beta(1) = 1;
for i=1:10
    %Number of terms in the response surface for nvar variable problem
    n_term = n_term*(nvar+i); 
    n_beta(i+1) = n_term/factorial(i);
end
nbeta1 = n_beta(npoly+1);

[ndata,n] = size(vrn);

%Create X
X(:,1) = ones(ndata,1);
X = [X vrn];
X1 = vrn;
nc1 = nvar;
n_loc = [1:nvar];
n_loc1 = 1;
for i = 2:npoly
    [nr,nc]=size(X1);
    clear X2;
    ctr=1;
    for k=1:nvar
        l_ctr = 0;
        for j=n_loc(k):nc
            X2(:,ctr) = vrn(:,k).*X1(:,j);
            ctr = ctr+1;
            l_ctr = l_ctr+1;
        end
        n_loc1(k+1)=l_ctr + n_loc1(k);
    end
    nc1 = nc;
    X = [X X2];
    X1 = X2;
    n_loc = n_loc1;
end