%2) Detect and remove outlier spectra based on residuals determined by RQ-GPR regression
%3) Detect and remove outlier spectra based on anomalies
%*****************************************************************************************************

close all
clear

for k=1:10
tic
% Definitions spectra and glucose reading arrays to operate on ********************
load('spect_d1s0119c_acdc');
xo=xacdc0119;
yo=y0119;
idacq = "d1" %select algorithm rq4,3,2,1
idses = "s"+"0122c_acdc" %select session(s)
idr0=40*[1:50]; %these test spectra to remove before bad spectra removal for all sessions
%idr0=15*[1:10]; %these test spectra to remove before bad spectra removal fr single session
%***********************************************************************************************

[N M] = size(xo);
avgxo = mean(xo);
stdxo=std(xo);
%*********************************************************************************

%Detect outlier spectra via residuals and remove (2) ***********************
tic
k2=5*k;
%k2 =20; %select mg/dL threshold value for residuals
%perform calibration to determine residuals
[trainedModel, validationRMSE] = trainRegressionModelRQ(xo, yo);
partitionedModel = crossval(trainedModel.RegressionGP, 'KFold', 5);
zoqr = kfoldPredict(partitionedModel);
residual = yo-zoqr;
TMrq0=trainedModel; %save model

% Identify
[I J] = size(residual);
for i=1:I
    if abs(residual(i)) >= k2 oids2(i) = 1; else oids2(i) = 0; end
end
idr2=find(oids2);

% Remove
xo2=xo;
yo2=yo;
xo2b=[];
yo2b=[];
xo2t=[];
yo2t=[];

%get outliers and test spectra
idremove2=flip(idr2); %remove lines starting from the bottom up 
for j=1:length(idremove2)
  xo2b(j,:) = xo(idremove2(j),:);
  yo2b(j) = yo(idremove2(j));
end

idremove0=flip(idr0); %remove lines starting from the bottom to avoid messing up
for j=1:length(idremove0)
  xo2t(j,:) = xo(idremove0(j),:);
  yo2t(j) = yo(idremove0(j));
end

idremove=flip(union(idr2,idr0)); %remove lines starting from the bottom to avoid messing up
for j=1:J
   xo2(idremove(j),:) = [];
   yo2(idremove(j)) = [];
end
R2 = 100*length(idr2) / N; 
toc

% Save data and parameters
yo2t=yo2t';
yo2b=yo2b;
fndat = "dat_"+ idacq+idses;
save(fndat)
%*******************************************************************************************

%perform calibration on xo2,yo2 (2)
tic
[trainedModel, validationRMSE] = trainRegressionModelRQ(xo2, yo2);
partitionedModel = crossval(trainedModel.RegressionGP, 'KFold', 5);
zorq2 = kfoldPredict(partitionedModel);
figure
%plot(zorq2, yo2, 'bo')
mdl = fitlm(zorq2,yo2); %returns a linear regression model of the responses y, fit to the data matrix X.
anova(mdl,'summary'); %test fit by anova
plot(mdl)
title('mdlrq2')
MDLrq2 = mdl;
TMrq2 = trainedModel;
MARDrq2 = 100*mean(abs(yo2-zorq2)./zorq2);
SEPrq2 = validationRMSE;

fncal = "cal_"+ idacq+idses+k;
save(fncal);
%fnbad= "bad_"+ idacq+idses;
%save(fnbad,'idacq','idses','TMrq2','xo2b','xo3b','yo2b','yo3b','idr0','idr2','idr3','avgxo','avgxo2','avgxo3');
toc
end
%*******************************************************************************************


