%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Detect and remove outlier spectra based on residuals
%Detect and remove outlier spectra based on anomalies
%*****************************************************************************************************

close all
clear;

% Definitions spectra and glucose reading arrays to operate on ********************
%idacq = "d1"; %select acquisition
%idses = "s0122c2"; %select session
%fnin = "spect_"+idacq+idses;
%load(fnin);
load('spe_d1r0c0sAlln'); load('sgl_vall');
xo=x0119n;
yo=y0119;
gid=
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 =25; %select mg/dL threshold value for residuals
%perform calibration to determine residuals
[trainedModel, validationRMSE] = trainRegressionModelE(xo, yo);
partitionedModel = crossval(trainedModel.RegressionGP, 'KFold', 5);
zoe0 = kfoldPredict(partitionedModel);
residualoe0= yo-zoe0;
TMe0=trainedModel; %save model
resper = 100*k2 / max(residualoe0);

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

% Remove
xo2=xo;
yo2=yo;
xo2b=[];
yo2b=[];
xtest2=[];
ytest2=[];

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

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

idremove=flip(union(idr2e,idr0)); %remove lines starting from the bottom to avoid messing up
for j=1:J
   xo2(idremove(j),:) = [];
   yo2(idremove(j)) = [];
end

R2 = length(idr2e) / N; 
toc
%*******************************************************************************************

% Detect outlier spectra via anomaly detecttion and remove (3) ***********
tic
k3=0.05; %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; 
%*******************************************************************************************

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

% Plot data
figure;plot(xo')
figure;plot(avgxo')
figure;plot(stdxo')

figure
plot(xo','r+')
hold on
plot(xo2','bo')
hold off
title("Red = Removed Spectra via Residuals")
avgxo2 = mean(xo2);
stdxo2=std(xo2);
figure;plot(avgxo2')
figure;plot(stdxo2')

%*******************************************************************************************
com23 = intersect(idr2e,idr3); cp23=100*length(com23)/length(union(idr2e,idr3));
percentidremove = 100*com23/23;
ytest2=ytest2';
yo2b=yo2b';
% Save data and parameters
%fndat = "dat_"+ idacq+idses;
%save(fndat)
%*******************************************************************************************
%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%Performs calibration on  xo2,yo2
%*****************************************************************************************************
%close all;
%clear;

%perform calibration on xo2,yo2 (2)
tic
[trainedModel, validationRMSE] = trainRegressionModelE(xo2, yo2);
partitionedModel = crossval(trainedModel.RegressionGP, 'KFold', 5);
zoe2 = kfoldPredict(partitionedModel);
figure
mdlcal = fitlm(zoe2,yo2); %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')
MDLe2 = mdlcal;
TMe2 = trainedModel;
MARDe2 = 100*mean(abs(yo2-zoe2)./zoe2);
SEPe2 = validationRMSE;
toc

[N M] = size(zoe2);
rn=[1:N];
rn=rn';
figure; plot(rn,yo2,'-ro',rn,zoe2,'-bo')
%fncal = "cal_"+ idacq+idses;
%save(fncal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 %Performs prediction using calibration data saved in trainedModel 
%********************************************************************************

%close all;
%clear;
%idtst = "s2022c"
%fntest="spect_"+idacq+idses;
%load(fntest); %load test spectra for testing
%xtest=x2022;
%ytest=y2022;
%xtest=xtest2;
%ytest=ytest2;
%xtest=xo2b;
%ytest=yo2b;
%xtest=xo4;
%ytest=yo4;
xtest=x2022n;
ytest=y2022;

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 = TMe2 %select calibration model
w(ip) = trainedModel.predictFcn(v(ip,:));
RDp(ip) = 100*(w(ip)-u(ip))/u(ip);
end
%u=u';
MARDp = 100*mean(abs(w-u)./u);
w=w';
[totalClark, percentClark] = clarke(u,w);

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')

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

%fnpred = "prd_" +idacq+idses+idtst+'x'; 
%save(fnpred);
save("prd_d1r0c0svAllonet0Group"+gid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

