%% [GP]: MATLAB Code for Gaussian Process
 
 
function rul=GP(funcOrd,h0)
clear global; global DegraUnit ...
TimeUnit time y thres signiLevel ns ny nt xTrain yTrain X dof 
%=== PROBLEM DEFINITION 1 (Required Variables) ==============
WorkName=' ';            % work results are saved by WorkName
DegraUnit=' ';                             % degradation unit
TimeUnit=' ';               % time unit (Cycles, Weeks, etc.)
time=[ ]';     %[cv]: time at both measurement and prediction
y=[ ]';                              %[ny x 1]: measured data
nT= ;        % num. of training set including the current one
nt=[ ]';%[nT x 1]: num. of training data in each training set
thres= ;                         % threshold (critical value)
degraTrue=[ ]';             %[cv]: true values of degradation
signiLevel= ;          % significance level for C.I. and P.I.
ns= ;                           % number of particles/samples
%=== PROBLEM DEFINITION 2 (Training Matrix) =================
xTrain=[ ]; yTrain=[ ];
%============================================================
% % % PROGNOSIS using GP
ny=length(y);
X=GLOBAL(xTrain,funcOrd);         %% Determine Hyperparameter
dof=size(yTrain,1)-size(X,2);
hMax=2; hBound=hMax*ones(1,length(h0));
hH=fmincon(@FUNC,h0,[],[],[],[],-hBound,hBound);
[R,thetaH,sigmaH]=RTHESIG(hH);         %% R, Theta, and Sigma
%=== PROBLEM DEFINITION 2 (Input Matrix for Prediction ======
for k=1:length(time(ny:end));       %% Degradation Prediction
 xNew=[ ];
 [zMean(k,:),zSig(k,:)]= ...
                      GPSIM(xNew,funcOrd,hH,R,thetaH,sigmaH); 
end;  
%============================================================
tDist=trnd(dof,size(zMean,1),ns);   %% Prediction Uncertainty
zHat0=repmat(zMean,1,ns)+tDist.*repmat(zSig,1,ns);
zHat=zHat0*(max(y)-min(y))+min(y);
degraPredi=zHat;
% % % POST-PROCESSING
rul=POST([],degraPredi,degraTrue);       %% RUL & Result Disp
Name=[WorkName ' at ' num2str(time(ny)) '.mat']; save(Name);
end
% % % GLOBAL FUNCTION
function X=GLOBAL(x0,funcOrd)
 nx=size(x0,1);
 if funcOrd==0; X=ones(nx,1);
 elseif funcOrd==1; X=[ones(nx,1) x0];
 elseif funcOrd==2; X=[ones(nx,1) x0 x0.^2];
 end
end
% % % OBJECTIVE FUNCTION TO FIND H
function objec=FUNC(h)
 global dof
 [R,~,sigma]=RTHESIG(h);
 objec=log(sigma^(2*dof)*det(R));
end
% % % R, THETA, and SIGMA
function [R,theta,sigma]=RTHESIG(h)
 global xTrain yTrain X dof
 R=CORREL(xTrain,h); Rinv=R^-1;
 theta=(X'*Rinv*X)\(X'*Rinv*yTrain);
 sigma=sqrt(1/dof*((yTrain-X*theta)'*Rinv*(yTrain-X*theta)));
end
% % % CORRELATION FUNCTION
function R=CORREL(x0,h)
 global xTrain
 for k=1:size(xTrain,1); for l=1:size(x0,1);
  R(k,l)=exp(-(norm(xTrain(k,:)-x0(l,:))/h)^2);
 end; end;
end
% % % GP SIMULATION at NEW INPUTS
function [zMean,zSig]=GPSIM(x,funcOrd,h,R,theta,sigma)
 global yTrain X
 Rinv=R^-1;
 xi=GLOBAL(x,funcOrd);
 r=CORREL(x,h);
 zMean=xi*theta+r'*Rinv*(yTrain-X*theta);
 w=Rinv*r+Rinv*X*((X'*Rinv*X)\(xi'-X'*Rinv*r));
 zSig=sigma*sqrt(w'*R*w-2*w'*r+1);   
end