function rul=GP_Battery(funcOrd,h0)
% Ex) rul=GP_Battery(1, 0.05);
clear global; global DegraUnit ...
TimeUnit time y thres signiLevel ns ny nt xTrain yTrain X dof 
%=== PROBLEM DEFINITION 1 (Required Variables) ==============
WorkName='Battery_GP';   % work results are saved by WorkName
DegraUnit='C/1 Capacity';                  % degradation unit
TimeUnit='Cycles';          % time unit (Cycles, Weeks, etc.)
time=[0:5:200]';    % time at both measurement and prediction
y=[1.00 0.99 0.99 0.94 0.95 0.94 0.91 0.91 0.87 0.86]';
%y=exp(-0.003.*time(1:10));
nT=1;        % num. of training set including the current one
nt=[10]';%[nT x 1]: no. of training data in each training set
thres=0.7;                       % threshold (critical value)
degraTrue=exp(-0.003.*time);     % true values of degradation
signiLevel=5;          % significance level for C.I. and P.I.
ns=5e3;                         % number of particles/samples
%=== PROBLEM DEFINITION 2 (Training Matrix) =================
y0=(y-min(y))/(max(y)-min(y));
time0=(time-min(time))/(max(time)-min(time));
xTrain=time0(1:length(y)); yTrain=y0;
%============================================================
% % % 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=time0(ny-1+k);
 [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)^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