%% [Metric]: MATLAB code for prognostics metrics
 
%=== Required Information ===================================
rul=[ ]';                                %[k x 1] or [k x ns]
cycle=[ ]';                        %[k x 1] at rul prediction
eolTrue=[ ];                                       % true EOL
t_s=[ ];                                     % starting cycle
t_e=[ ];                                       % ending cycle
alpha=[ ];
lambda=[ ];
%============================================================
%%% PH: Prognostic Horizon
rulTrue=eolTrue-cycle;
alphaZone=eolTrue*alpha;
loca=find( rulTrue-alphaZone<rul & rul<rulTrue+alphaZone );
a=find( [0 diff(loca')]>1,1,'last' ); if isempty(a); a=1; end
ph=eolTrue-cycle(loca(a));
disp(['PH= ' num2str(ph)]);
%%% ALA: Alpha-Lambda Accuracy
t_lam=t_s+(eolTrue-t_s)*lambda;
[~, idx_lam]=min(abs(cycle-t_lam)); t_lam=cycle(idx_lam);
a=find( rul(idx_lam)>rulTrue(idx_lam)*(1-alpha) ...
      & rul(idx_lam)<rulTrue(idx_lam)*(1+alpha) );
if isempty(a); 
    disp('ALA: False'); 
else disp ('ALA: True');
end
%%% RA & CRA: Relative & Cumulative Relative Accuracy 
ra=1-abs(rulTrue(idx_lam)-rul(idx_lam))/rulTrue(idx_lam);
si=find(cycle>=t_s,1); ei=find(cycle<=t_e,1,'last');
cra=mean(1-abs(rulTrue(si:ei)-rul(si:ei))./rulTrue(si:ei));
disp(['RA= ' num2str(ra) ', CRA= ' num2str(cra)])
%%% Cvg: Convergence
errCurve=abs( rulTrue(si:ei)-rul(si:ei) )./rulTrue(si:ei);
cycle1=[cycle(si+1:ei); eolTrue];
numer=sum( (cycle1 - cycle(si:ei)).*errCurve );
xc=0.5*sum( (cycle1.^2 - cycle(si:ei).^2).*errCurve )/numer;
yc=0.5*sum( (cycle1 - cycle(si:ei)).*errCurve.^2 )/numer;
cvg=sqrt( (xc-t_s)^2 + yc^2 );
disp(['Cvg= ' num2str(cvg)])