
%1) Detect and remove outlier spectra based on anomalies using One Cluster SVM
%*****************************************************************************************************
close all
clear;

% Definitions spectra and glucose reading arrays to operate on ********************
load('spe_d1r0c0s0122bDCn2'); 
load('sgl_y2b2');

%Select one group for processing at a time
%Group=1; xind=x1dc06n; yind=y206; xo=x1dc10141720n; yo=y210141720; 
%Group=2; xind=x1dc10n; yind=y210; xo=x1dc06141720n; yo=y206141720;
%Group=3; xind=x1dc14n; yind=y214; xo=x1dc06101720n; yo=y206101720; 
%Group=4; xind=x1dc17n; yind=y217; xo=x1dc06101420n; yo=y206101420; 
Group=5; xind=x1dc17n; yind=y220; xo=x1dc06101417n; yo=y206101417; 

k3 = 0.1; %select possible fraction of outliers
k2 = 0.2; %select % mg/dL threshold value for residuals

gid=Group;
idtest =[];

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

%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');
%*******************************************************************************************

