function rul=GP_Crack(funcOrd,h0)
% Ex) rul=GP_Crack(1, 0.2);
clear global; global DegraUnit ...
TimeUnit time y thres signiLevel ns ny nt xTrain yTrain X dof 
%=== PROBLEM DEFINITION 1 (Required Variables) ==============
WorkName='Crack_GP';   % work results are saved by WorkName
DegraUnit='Crack size (m)';                % degradation unit
TimeUnit='Cycles';          % time unit (Cycles, Weeks, etc.)
time=[0:100:2500, 0:100:1600, 1700:100:4000]';
y=[0.01 0.0104 0.0107 0.0111 0.0116 0.0120 0.0125  ...
 0.0131 0.0137 0.0143 0.0150 0.0158 0.0166 0.0176   ...
 0.0186 0.0198 0.0211 0.0226 0.0243 0.0263 0.0287   ...
 0.0314 0.0347 0.0387 0.0437 0.0501 ...
 0.0100 0.0103 0.0106 0.0109 0.0112 0.0116 0.012 0.0124 ...
 0.0128 0.0133 0.0138 0.0143 0.0149 0.0155 0.0162 0.0169]';
nT=2;        % num. of training set including the current one
nt=[26 16]';      % no. of training data in each training set
thres=0.05;                      % threshold (critical value)
mTrue=3.8; cTrue=1.5e-10; aTrue=0.01; dsig=75; t=0:100:3500;
degraTrue=(t'.*cTrue.*(1-mTrue./2).*(dsig*sqrt(pi)).^mTrue...
 +aTrue.^(1-mTrue./2)).^(2./(2-mTrue));
signiLevel=5;          % significance level for C.I. and P.I.
ns=5e3;                         % number of particles/samples
%=== PROBLEM DEFINITION 2 (Training Matrix) =================
nx=3;                    % num. of previous data as the input
y0=(y-min(y))/(max(y)-min(y));
nm=size(y,2); a=0; b=0;
for l=1:nT;
 for k=1:nt(l)-nx;
  xTrain(k+a,1:nx*nm)=reshape(y0(k+b:(k-1)+nx+b,:),1,nx*nm);
  yTrain(k+a,:)=y0(k+nx+b,:);
 end; a=size(xTrain,1); b=sum(nt(1:l));
end
%============================================================
% % % 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 ======
input=y0(end-nx:end,:);
for k=1:length(time(ny:end));
 xNew=reshape(input(k:(k-1)+nx,:),1,nx*nm);
 [zMean(k,:),zSig(k,:)]= ...
                      GPSIM(xNew,funcOrd,hH,R,thetaH,sigmaH);
 if k>1 input(k+nx,:)=zMean(k,:); end
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)^1.9);
 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