%
%2) Removes spectra based on residuals determined by Exp-GPR regression
%*****************************************************************************************************

close all
clear;

% Definitions spectra and glucose reading arrays to operate on ********************
load('spe_d1r0c0s030610141720caln');
load('sgl_d1r0c0s030610141720b');
%load('dat_d1r0c0s0119cad','TMe0');
%load('spe_d1r0c0s0122baln','x1ad0119n');
%load('sgl_d1r0c0s0122b','y0119');
xo=x1dc030610141720n; %x0119(:,[1:28 35:64]);
yo=y1030610141720;
%xo=x1ad0119n; %x0119(:,[1:28 35:64]);
%yo=y0119;
idtest=28*[1:30]; %these test spectra to remove before bad spectra removal
%***********************************************************************************************

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

xo2=xo;
yo2=yo;

idremovet=idtest;
for j=1:length(idremovet)
  xot(j,:) = xo(idremovet(j),:);
  yot(j) = yo(idremovet(j));
end
for j=1:length(idremovet)
   xo2(idremovet(j),:) = [];
   yo2(idremovet(j)) = [];
end
yot=yot';
%*********************************************************************************

%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] = trainRegressionModelSEGPR(xo2, yo2);
partitionedModel = crossval(trainedModel.RegressionGP, 'KFold', 5);
zoe0 = kfoldPredict(partitionedModel);
residualoe0= yo2-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=flip(idr2);
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
yo2b=yo2b';
xcal=xo2;
ycal=yo2;
R2 = 100*length(idr2) / N; 
toc
%*******************************************************************************************

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

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

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

xcaln = xob;
%*******************************************************************************************

% Save data and parameters
%save('cal_model_res','TMe0');
save('dat_d1r0c0s0119cdc');
%*******************************************************************************************

