SR_GPR.m

%%______________________________________________________________________________
%% This function is to implement GP training and prediction
%% This function is created by Kai Zhou 7/31/2019
%%______________________________________________________________________________
function [y_star, y_star_cov,b_opt] = SR_GPR(xL,yL,x_star)
global k1
k1 = @ (x1, x2, b1, fi1) fi1^2 * exp(-b1 * (x1 - x2)' * (x1 - x2));

%%Note, bounds need to be manually adjusted in order to aviod numerial
%%instability
lb = zeros(1,2)+1e-3;
ub = zeros(1,2)+2000;
ub(end) = 5;
%%_________

global x1
global y1
global numInp
x1 = xL;
y1 = yL;

options = optimoptions('particleswarm','MinNeighborsFraction',1);
options.SwarmSize = 100;
options.MaxIterations = 30000;
tic
[b_opt,fval,exitflag,output] = particleswarm(@Obj_PSO,2,lb,ub,options);
formatstring = 'particleswarm reached the value %f using %d function evaluations.\n';
fprintf(formatstring,fval,output.funccount)

fi1_opt = b_opt(end);
b1_opt = b_opt(1);


nTest = length(x_star(:,1));
numLow = length(xL(:,1));
Temp = zeros(length(x_star(:,1)),1)+1;
Temp1 = zeros(length(xL(:,1)),1)+1;
H = [Temp x_star];
HH = [Temp1 xL];


for i = 1 : nTest
    for j = 1 : numLow
        t(i, j) = k1(x_star(i, :)', xL(j, :)', b1_opt,fi1_opt);
    end
end


for i = 1 : numLow
    for j = 1 : numLow
        V(i,j) =  k1(xL(i, :)', xL(j, :)', b1_opt,fi1_opt);
    end
end


Temp3 = HH'*((V+1e-9*eye(size(V)))\HH);
Beta = ((Temp3 + 1e-9*eye(size(Temp3)))\HH')*((V+1e-9*eye(size(V)))\yL);
Beta = zeros(numInp+1,1);
y_star = H*Beta + t*inv(V+1e-8 * eye(size(V)))* (yL-HH*Beta);


for i1 = 1:nTest
    for i2 = 1:nTest
        C_1(i1,i2) =  k1(x_star(i1, :)', x_star(i2, :)', b1_opt,fi1_opt);
    end
end
y_star_cov = C_1 - t/(V + 1e-8 * eye(size(V)))*t';
end