%
%1) Removes spectra based on difference from the mean
%2) Removes spectra based on residuals determined by Exp-GPR regression
%3) Removes spectra based on anomalies
%4) Removes spectra based on x*b' by LR-PLS regression and difference from mean
%5) Removes spectra based on xtrans of Exp-GPR and difference from mean
%6) Removes spectra based on local outlier function
%*****************************************************************************************************

close all
clear;

% Definitions spectra and glucose reading arrays to operate on ********************
load('spect_d1s0119c');
xo=x0119;
yo=y0119;
idacq = "d1";%select algorithm rq4,3,2,1
idses = "s"+"0119c_CAL"; %select session(s)
idtest=[ ]; %45*[1:50]; %these test spectra to remove before bad spectra removal
%***********************************************************************************************

avgxo = mean(xo);
stdxo=std(xo);
idr0 = flip(idtest);

%Detect outlier spectra via dfifference from mean and remove (1) ***********************
tic
k1a = 0.4;% Define removal parameters
k1b = 15; % Define removal parameters
[N, M] = size(xo);
for n=1:N
    for m=1:M
        if abs(xo(n,m) - avgxo(m)) > abs(k1a*xo(n,m)); loc1(n,m) = 1; else loc1(n,m) = 0; end
    end
end

% Identify
loc1=loc1';
S1=sum(loc1);
S1=S1';
S1x=sum(S1);
for n=1:N
      if S1(n) >= k1b oids1(n) = 1; else oids1(n) = 0; end
end
idr1=find(oids1);

% Remove
xo1=xo;
yo1=yo;
xo1b=xo;
yo1b=yo;
idremove1=flip(idr1); %remove lines starting from the bottom to avoid messing up

for j=1:length(idremove1)
  xo1b(j,:) = xo(idremove1(j),:);
  yo1b(j) = yo(idremove1(j));
end
for j=1:length(idremove1)
   xo1(idremove1(j),:) = [];
   yo1(idremove1(j)) = [];
end

R1 = 100*length(idr1) / N; 
toc
%*********************************************************************************

%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
idr2=find(oids2);

% Remove
xo2=xo;
yo2=yo;
idremove2=flip(union(idr2,idr0)); %remove lines starting from the bottom to avoid messing up
%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

R2 = 100*length(idr2) / N; 
toc
%*******************************************************************************************

% Detect outlier spectra via anomaly detecttion and remove (3) ***********
tic
k3=0.08; %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;
idremove3=flip(union(idr3,idr0)); %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 = 100*length(idr3) / Na; 
toc
%*******************************************************************************************

%Detect outlier spectra via x*b' and remove (4) ********************
tic
k4a = 0.6;
k4b = 10;
xo4=xo;
yo4=yo;
% perform calibration on xo4,yo4 to detect bad spectra
[trainedModel, validationRMSE] = trainRegressionModelLR(xo4,yo4);
zolr4 = trainedModel.LinearModel.Fitted;
T = trainedModel.LinearModel.Coefficients(2:59,1);
TMlp0 = trainedModel;
b4 = table2array(T);
xb4=xo4.*b4'; %

% Identify
avgxb4 = mean(xb4);
[N, M] = size(xb4);
for n=1:N
    for m=1:M
        if abs(xb4(n,m) - avgxb4(m)) > abs(k4a*xb4(n,m)); loc4(n,m) = 1; else loc4(n,m) = 0; end
    end
end

loc4=loc4';
S4=sum(loc4);
S4=S4';
S4x=sum(S4);

for n=1:N
      if S4(n) >= k4b oids4(n) = 1; else oids4(n) = 0; end
end
idr4= find(oids4);

% Remove
xo4=xo;
yo4=yo;
idremove4=flip(union(idr4,idr0)); %remove lines starting from the bottom to avoid messing up

for j=1:length(idr4)
  xo4b(j,:) = xo(idremove4(j),:);
  yo4b(j) = yo(idremove4(j));
end
for j=1:length(idremove4)
   xo4(idremove4(j),:) = [];
   yo4(idremove4(j)) = [];
end

R4 = 100*length(idr4) / N; 
toc
%*****************************************************************************************************

%Detect outlier spectra via xtrans and remove  (5) ********************
tic
k5a = 0.025;
k5b = 20;
xo5=xo;
yo5=yo;
Xnew=xo;

    %0. Get cal info
    gpr=TMe0.RegressionGP; %existing calibration structure
    % 1. Apply standardization to new data if needed. 
    % function [ypred, ysd] = predictWithEquation(gpr, Xnew, standardize)
    standardize=1;
    muX     = gpr.PredictorLocation;
    sigmaX = gpr.PredictorScale;
    if standardize
        XnewStd = bsxfun(@rdivide,bsxfun(@minus,Xnew,muX'),sigmaX');
    else
        XnewStd = Xnew;
    end

    % 2. Get estimated coefficients.
    alphaHat = gpr.Impl.AlphaHat; %2042x1
    betaHat   = gpr.Impl.BetaHat;   %1x1
    thetaHat  = gpr.Impl.ThetaHat;  %3x1 
    sigmaHat = gpr.Impl.SigmaHat; %1x1

    % 3. Get basis function and kernel function.
    HFcn   = gpr.Impl.HFcn; %structure
    kernel  = gpr.Impl.Kernel;%structure

    % 4. Make function to evaluate kernel products. kfun below can be
    % called like this: KNM = kfun(XN,XM).
    kfun = makeKernelAsFunctionOfXNXM(kernel,thetaHat); %structure

    % 5. Make function to get diagonal of kernel. kfundiag below can be
    % called like this: diagKNN = kfundiag(XN)
    kfundiag = makeDiagKernelAsFunctionOfXN(kernel,thetaHat); %structure

    % 6. Get active set and Cholesky factor for further computation.
    XA = gpr.Impl.ActiveSetX; %2042x58
    L  = gpr.Impl.LFactor; %2042x2042

    % 7. Get predicted values.
    KXnewXA = kfun(XnewStd,XA); %2478x2042

    %x transformed by the model is xs
    xtrans = HFcn(XnewStd)*betaHat / length(alphaHat) + KXnewXA.*alphaHat'; %2478x2042
    ypred  = HFcn(XnewStd)*betaHat + KXnewXA*alphaHat; %2478x1

    % 8. Get standard deviation of new responses.
    gamma = L \ KXnewXA'; %2042x2478
    ysd   = sqrt(max(0,sigmaHat^2 + kfundiag(XnewStd) - sum(gamma.^2,1)')); %2478x1

xb5=xtrans; %xb4=xo4.*b4'; %


% Identify
avgxb5 = mean(xb5);
[N, M] = size(xb5);
for n=1:N
    for m=1:M
        if abs(xb5(n,m) - avgxb5(m)) > abs(k5a*xb5(n,m)); loc5(n,m) = 1; else loc5(n,m) = 0; end
    end
end

loc5=loc5';
S5=sum(loc5);
S5=S5';
S5x=sum(S5);

for n=1:N
      if S5(n) >= k5b oids5(n) = 1; else oids5(n) = 0; end
end

idr5 = find(oids5);

% Remove
xo5=xo;
yo5=yo;
idremove5=flip(union(idr5,idr0)); %remove lines starting from the bottom to avoid messing up

for j=1:length(idremove5)
  xo5b(j,:) = xo(idremove5(j),:);
  yo5b(j) = yo(idremove5(j));
end
for j=1:length(idremove5)
   xo5(idremove5(j),:) = [];
   yo5(idremove5(j)) = [];
end

R5 = 100*length(idr5) / N; 
toc
%*****************************************************************************************************

%Local Outlier Function *******************************************************************************
tic
k6=0.05;

%Train a local outlier factor model by using the lof function. Specify the fraction of outliers (ContaminationFraction) as 0.05, 500 nearest neighbors, and the Mahalanobis distance.
[LOFObj,tf_lof,s_lof] = lof(xo, ContaminationFraction=k6, NumNeighbors=500,Distance="mahalanobis");
%LOFObj is a LocalOutlierFactor object. lof also returns the anomaly indicators (tf_lof) and anomaly scores (s_lof) for the data (feat). lof determines the score threshold value (LOFObj.ScoreThreshold) so that the function detects the specified fraction of observations as outliers.
%Plot a histogram of the score values. Create a vertical line at the score threshold corresponding to the specified fraction.

%figure
%histogram(s_lof,Normalization="probability")
%xline(LOFObj.ScoreThreshold,"k-",  join(["Threshold =" LOFObj.ScoreThreshold]))
%title("Histogram of Anomaly Scores for Local Outlier Factor")

% Identify
idr6=find(tf_lof);
idr6=idr6';

% Remove
xo6=xo;
yo6=yo;
xo6b=[];
yo6b=[];
idremove6=flip(idr6); %remove lines starting from the bottom to avoid messing up

for j=1:length(idremove6)
  xo6(j,:) = xo(idremove6(j),:);
  yo6(j,:) = yo(idremove6(j),:);
end
for j=1:length(idremove6)
   xo6(idremove6(j),:) = [];
   yo6(idremove6(j),:) = [];
end

R6 =100*length(idr6) / Na; 
toc
%*******************************************************************************************

com21 = intersect(idr2,idr1); cp21=100*length(com21)/length(union(idr2,idr1));
com23 = intersect(idr2,idr3); cp23=100*length(com23)/length(union(idr2,idr3));
com24 = intersect(idr2,idr4); cp24=100*length(com24)/length(union(idr2,idr4));
com25 = intersect(idr2,idr5); cp25=100*length(com25)/length(union(idr2,idr5));
com26 = intersect(idr2,idr6); cp26=100*length(com26)/length(union(idr2,idr6));
%*******************************************************************************************

% Remove group*******************************************************************************
xcal=xo;
ycal=yo;
xcalb=xo;
ycalb=yo;
idrt = union(idr1,idr2);
idrt = union(idrt,idr3);
idrt = union(idrt,idr4);
idrt = union(idrt,idr5);
idrt = union(idrt,idr6);
idremovet=flip(idrt); %remove lines starting from the bottom to avoid messing up

for j=1:length(idremovet)
  xcalb(j,:) = xo(idremovet(j),:);
  ycalb(j) = yo(idremovet(j));
end
for j=1:length(idremovet)
   xcal(idremovet(j),:) = [];
   ycal(idremovet(j)) = [];
end

Rt = 100*length(idrt) / N; 
%**********************************************************************************************

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

method1=100*cp21/length(idr1)
method3=100*cp23/length(idr3)
method4=100*cp24/length(idr4)
method5=100*cp25/length(idr5)
method6=100*cp26/length(idr6)

com54 = intersect(idr5,idr6); cp54=100*length(com54)/length(union(idr5,idr4));
com56 = intersect(idr5,idr6); cp56=100*length(com56)/length(union(idr5,idr6));