%*****************************************************************************************************
close all
clear;

load('spe_d1r0c0sAll64n'); load('sgl_vall'); 

%outlier thresholds k=0
%k3 = 0.05; %select possible fraction of anomalies
%k2 = 0.25; %select % mg/dL threshold value for residuals
%outlier thresholds k=1
k3 = 0.14; %select possible fraction of anomalies
k2 = 0.35; %select % mg/dL threshold value for residuals


%Select one group for processing at a time
gid=14; xind=x1pd2022n; yind=y2v2022; xo=x1pd0119n; yo=y2v0119; 
%gid=16; xind=x1dc2022n; yind=y2v2022; xo=x1dc0119n; yo=y2v0119; 


idtest =[];
gspectra = "DC";
gluco = "ONET"; %"onet";
%*********************************************************************************************

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1) Detect and remove outlier spectra based on anomalies using One Cluster SVM

%Fig1
figure;plot(xo')
%Fig2
figure;
avgxo = mean(xo);
plot(avgxo')
%Fig3
figure;
stdxo=std(xo);
plot(stdxo')
%*********************************************************************************************

idr0 = flip(idtest);
[Na Ma] = size(xo);
%*******************************************************************************************

% Detect outlier spectra via anomaly detecttion and remove (3) ***********
tic
%k3=0.2; %select possible fraction of outliers
[Na Ma] = size(xo);
[Mdl,tfa,sa] = ocsvm(xo, ContaminationFraction=k3, KernelScale="auto", StandardizeData=true);
%figure
%histogram(sa,Normalization="probability")
%xline(Mdl.ScoreThreshold,"k-", join(["Threshold =" Mdl.ScoreThreshold]))
%title("Histogram of Anomaly Scores for One-Class SVM")

rng("default") % For reproducibility
[T,loss] = tsne(xo,Standardize=false,Perplexity=10,Exaggeration=2);
%figure
%gscatter(T(:,1),T(:,2),tfa,"kr")
%title("t-distributed stochastic neighbor embedding")
%L = legend("Normal Points","Outliers");
%xlabel("Column T(1)")
%ylabel("Column T(2)")

[idx, scoresf] = fsulaplacian(xo);
%figure
%gscatter(xo(:,idx(1)),xo(:,idx(2)),tfa,"kr")
%title("Feature ranking by Laplacian scores")
%L = legend("Normal Points","Outliers");
%xlabel("Column idx(1)")
%ylabel("Column idx(2)")

% Identify
oids3 = tfa';
idr3= find(oids3);

% Remove
xo3=xo;
yo3=yo;
xo3b=[];
yo3b=[];
idremove3=flip(idr3); %remove lines starting from the bottom to avoid messing up

for j=1:length(idremove3)
  xo3b(j,:) = xo(idremove3(j),:);
  yo3b(j,:) = yo(idremove3(j),:);
end
for j=1:length(idremove3)
   xo3(idremove3(j),:) = [];
   yo3(idremove3(j),:) = [];
end
R3 = length(idr3) / Na; 
toc

%Fig4
figure
plot(xo','r+')
hold on
plot(xo3','b+')
title("Red = Removed Spectra via Anomaly")
%legend('Removed', 'Original');
hold off
%Fig5
avgxo3 = mean(xo3);
figure;plot(avgxo3')
%Fig6
stdxo3=std(xo3);
figure;plot(stdxo3')
%*******************************************************************************************

%********************************************************************************************
xcal1=xo; 
ycal1=yo; 
%xcal2=xo2
%ycal2=yo2;
xcal3=xo3;
ycal3=yo3;
%*******************************************************************************************

% Save data and parameters***********************************************************
%save('cal_model_res','TMe0');
%save('dat_d1r0c0sUnivPDS_EGPR64');
%*******************************************************************************************
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%2) Detect and remove outlier spectra based on residuals determined by Exp-GPR regression
%*****************************************************************************************************
xo=xcal3;
yo=ycal3;

%*********************************************************************************************

avgxo = mean(xo);
stdxo=std(xo);
idr0 = flip(idtest);
[N M] = size(xo);

xo2=xo;
yo2=yo;

%*********************************************************************************************
tic
%perform calibration to determine residuals
%k2 = 0.2; %select % mg/dL threshold value for residuals
[trainedModel, validationRMSE] = trainRegressionModelEGPR64(xo2, yo2);
partitionedModel = crossval(trainedModel.RegressionGP, 'KFold', 5);
zoe0 = kfoldPredict(partitionedModel);
residualoe0= yo2-zoe0;
TMe640=trainedModel; %save model
resper = 100*k2 / max(residualoe0);

% Identify
[I J] = size(residualoe0);
for i=1:I
  if  zoe0(i) - yo2(i)*(1+k2)  >= 0  oids2(i) = 1; else oids2(i) = 0; end
 end
idr2=find(oids2);

% Remove
idremove2=flip(idr2);
idremove2=idremove2';
for j=1:length(idremove2)
  xo2b(j,:) = xo(idremove2(j),:);
  yo2b(j) = yo(idremove2(j));
end
for j=1:length(idremove2)
   xo2(idremove2(j),:) = [];
   yo2(idremove2(j)) = [];
end
yo2b=yo2b';
toc
%***************************************************************************************

%xcal1=xo; 
%ycal1=yo; 
xcal2=xo2;
ycal2=yo2;
%xcal2=xo3;
%ycal2=yo3;
%*******************************************************************************************

% Save data and parameters
%save('cal_model_res','TMe0');
%save('dat_d1r0c0sPamDC');
%*******************************************************************************************
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%3) Perform calibration on  xcal,ycal using EGPR model
%*****************************************************************************************************

%xcal=xcal1;
%ycal=ycal1;
xcal=xcal2;
ycal=ycal2;
%xcal=xcal3;
%ycal=ycal3;
%xcal=xcal4;
%ycal=ycal4;

tic
[trainedModel, validationRMSE] = trainRegressionModelEGPR64(xcal, ycal);
partitionedModel = crossval(trainedModel.RegressionGP, 'KFold', 5);
zcal = kfoldPredict(partitionedModel);
figure
mdlcal = fitlm(ycal,zcal); %returns a linear regression model of the responses y, fit to the data matrix X.
anova(mdlcal,'summary'); %test fit by anova
plot(mdlcal)
title('mdle2 cal')

TMcal64 = trainedModel;
EMARDcal = 100*mean(abs(ycal-zcal)./zcal);
ERDcal=100*max((ycal-zcal)./zcal);
ESEPcal = validationRMSE;
toc

[N M] = size(zcal);
rn=[1:N];
rn=rn';
figure; plot(rn,ycal,'-ro',rn,zcal,'-bo');title('Validation'); xlabel('Record Number]'); ylabel('Concentration [mg/dL]');
legend('Reference', 'Predicted');

%save('cal_d1r0c0sPamDC');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%4) Find anomalies in new data sets given previously trained anomaly OCSVM.
%******************************************************************************************

tic
%Run data_anoalyocsvm first to get mdl, x, and anomalies
Xa = xo3; %this is selected in data_anomalyocsvm.mqw Ya=y03
Xn=xind; %select array other than x to test for novelties
Yn=yind; %select array other than x to test for novelties

threshold = Mdl.ScoreThreshold; %comes from data_anomalyocsvm.m
[Nn Mn] = size(Xa);

%object Mdl and returns the logical array tf, whose elements are true when an anomaly is
%detected in the corresponding row of x.n
rng("default") % For reproducibility 
[tfn,sn] = isanomaly(Mdl,Xn,ScoreThreshold= threshold); % A negative score value with large magnitude indicates a normal observation, and a large positive value indicates an anomaly.
figure
histogram(sn, Normalization="probability")
hold on 
histogram(sn,Normalization="probability")
xline(threshold,"k-",join(["Threshold =" threshold ]))
legend("Training data","Test data",Location="southeast")
title("Histograms of Anomaly Scores for One-Class SVM")
hold off

idn = find(tfn);
idn=idn';
NFn = 100*sum(tfn) / Nn

% Remove
%Xn=xo;
xdc2022b=[];
ydc2022b=[];
xdc2022n=Xn;
ydc2022n=Yn;
xo=Xn;
yo=Yn;
idremoven = flip(idn);
J=length(idremoven);
%get outliers and test spectra
idremoven=flip(idn); %remove lines starting from the bottom up 
for j=1:J
  xdc2022b(j,:) = xo(idremoven(j),:);
  ydc2022b(j) = yo(idremoven(j));
end

for j=1:J
   xdc2022n(idremoven(j),:) = [];
   ydc2022n(idremoven(j)) = [];
end
xprd1=xdc2022n;
yprd1=ydc2022n';

Rn = 100*length(idn) / Nn;
toc
%*******************************************************************************
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%5) Perform prediction using saved trained EGPR model 
%********************************************************************************

xtest=xprd1;
ytest=yprd1;
p=6;

v=[];
u=[];
w=[];

for ip=1:length(ytest)
v(ip,:) = xtest(ip,:); %select x spectra data to use for prediction
u(ip) = ytest(ip);; %select y glucose data to use for error calculation if available
%u(ip)=u(ip)';
trainedModel = TMcal64 %select calibration model
w(ip) = trainedModel.predictFcn(v(ip,:));
end

EMARDprd = 100*mean(abs(avp(w,p)./avp(u,p)-1));
ERDprd = 100*max((w-u)./u);

%figure 1
w=w';
[totalClark, percentClark] = clarke(u,w);

%figure 2
%mdlpred = fitlm(w,u); %returns a linear regression model of the responses y, fit to the data matrix X.
%anova(mdlpred,'summary'); %test fit by anova
%figure;plot(mdlpred)
%title('mdle2 pred')

%figure 3;
[N M] = size(w);
m=[1:N];
m=m';
figure; plot(m,u,'-ro',m,w,'-bo')

%************************************************************************
%Fig6
%figure;plot(m,uavg,'-ro',m,wavg,'-bo')
%title('Prediction'); xlabel('Record Number]'); ylabel('Concentration [mg/dL]');
%legend('Reference', 'Predicted');

%uavg = avp(u',p);
%wavg = avp(w,p); 
%[N2 M2] = size(uavg);
%m=[1:N2];

%Fig4
%[totalClark, percentClark] = clarke(uavg,wavg);

%Fig5
%mdlpred = fitlm(uavg,wavg); %returns a linear regression model of the responses y, fit to the data matrix X.
%anova(mdlpred,'summary'); %test fit by anova
%figure;plot(mdlpred);title('mdle2 pred')
%************************************************************************

save("prd_d1r0c0svAllonet1Group"+gid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


