%
%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('spe_d1r0c0s0122ban');
%load('sgl_d1r0c0s0122b');
load('dat_d1r0c0s0119cdc');
xo=x1dc2022n; %x0119(:,[1:28 35:64]);
%yo=yo2b;
yo=y2022;
idtest=[];%idtest=45*[1:50]; %these test spectra to remove before bad spectra removal
%***********************************************************************************************

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

%Detect outlier spectra via xtrans and remove  (5) ********************
tic
k5a = 4;
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.01;

%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) / N; 
toc
%*******************************************************************************************

% Remove group*******************************************************************************
xprd=xo;
yprd=yo;
xprdb=xo;
yprdb=yo;

idrt = union(idr5,idr6);
idremovet=flip(idrt); %remove lines starting from the bottom to avoid messing up

for j=1:length(idremovet)
  xprdb(j,:) = xo(idremovet(j),:);
  yprdb(j) = yo(idremovet(j));
end
for j=1:length(idremovet)
   xprd(idremovet(j),:) = [];
   yprd(idremovet(j)) = [];
end

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

%Subtracts min value and divides by value of first pixel
so = xprd;
figure;plot(so')

minX1 = min(so);
minX2 = min(minX1);
xoa = so - minX2;
figure;plot(xoa')

[N M] = size(so);
xtemp = xoa(:,1).*ones(N,M);
xob= xoa./xtemp;
figure;plot(xob');

xprdn = xob;
%*******************************************************************************************

% Save data and parameters
save('dat_d1r0c0s2022cdc');
%*******************************************************************************************

