
%2) Removes spectra based on residuals determined by Exp-GPR regression
%*****************************************************************************************************
%Mike Sessions: 1, 2, 7, 9, 12. Prediction: 7. Calibration: 1,2,9,12 
%Pam Sessions: (3), 6, 10, 14, 17, 20. Prediciton: 14. Calibration: 6, 10, 17, 20
%Jake Sessions: 4,13,15,18, 21. Prediction: 15. Calibration: 4,13,18,21
%Univ Sessions: all above, Prediciton 7, 14, 15. Calibration: 1,2,4,6,9,10,12,13,14,17,18,21
%using the 25hz 64pixel normalized PDS, i suggest the following for selecting outliers for removal from presentation to show expected performance with higher SNR.
%spe_d1r0c0s0122bPDn.mat (x1pdnn) and sgl_y2b.mat (onetouch/verio)
%x0119(:,[1:28 35:64]);
close all
clear;

% Definitions spectra and glucose reading arrays to operate on ********************
load('spe_d1r0c0s0122bPDn', 'x1pdPam','x1pd14');
load('sgl_y2b');
xo=x1pdPam; 
yo=y2Pam;
idtest =[];
%***********************************************************************************************

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

xo2=xo;
yo2=yo;

%*********************************************************************************************
tic
%perform calibration to determine residuals
k2 = 0.15; %k2 = 25; %select mg/dL threshold value for residuals
[trainedModel, validationRMSE] = trainRegressionModelEGPR64(xo2, yo2);
partitionedModel = crossval(trainedModel.RegressionGP, 'KFold', 5);
zoe0 = kfoldPredict(partitionedModel);
residualoe0= yo2-zoe0;
TMe640=trainedModel; %save model
resper = 100*k2 / max(residualoe0);

% Identify
[I J] = size(residualoe0);
for i=1:I
  if  zoe0(i) - yo2(i)*(1+k2)  >= 0  oids2(i) = 1; else oids2(i) = 0; end
 end
idr2=find(oids2);

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

%*********************************************************************************************
tic
xo3=xo2;
yo3=yo2;

%perform calibration to determine residuals
k3 = 0.1; %k2 = 25; %select mg/dL threshold value for residuals
[trainedModel, validationRMSE] = trainRegressionModelEGPR64(xo3, yo3);
partitionedModel = crossval(trainedModel.RegressionGP, 'KFold', 5);
zoe0 = kfoldPredict(partitionedModel);
residualoe0= yo3-zoe0;
TMe640=trainedModel; %save model
resper = 100*k3 / max(residualoe0);

% Identify
[I J] = size(residualoe0);
for i=1:I
  if  zoe0(i) - yo3(i)*(1+k3)  >= 0  oids3(i) = 1; else oids3(i) = 0; end
 end
idr3=find(oids3);

% Remove
idremove3=flip(idr3);
idremove3=idremove3';
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
yo3b=yo3b';
toc
%***************************************************************************************

xcal1=xo; 
ycal1=yo; 
xcal2=xo2;
ycal2=yo2;
xcal3=xo3;
ycal3=yo3;
%*******************************************************************************************

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

