
load protosessions session

outlier=struct([]);
tic
for sesid=1:length(session)
% sesid= 10; % Session ID******************************
% k = 1; % Ref Meter ID
kfactor1  = 1.5; % for Ralf Cook distance on spectra
kfactor2  = 1.5; % for Ralf Leverage distance on spectra
kfactor3  = 3.0; % for Quartile distance on spectra
kfactor4  = 1.5; % for Cook distance on pairs
kfactor5  = 1.5; % for Mahal distance on pairs
kfactor6  = 1.5; % for AbsDif distance on pair
kfactor7  = 1.5; % for Ralf res distance on spectra
kfactor8  = 1.5; % for Ralf lc_res distance on spectra
kfactor9  = 1.5; % for Ralf stud_res distance on spectra
kfactor10 = 1.5; % for Ralf Lof distance on spectraz
kfactor11 = 1.5; % for Ralf Cook distance on loop spectra
kfactor12 = 1.5; % for Ralf Leverage distance on loop spectra
kfactor13 = 1.5; % for Ralf res distance on loop spectra
kfactor14 = 1.5; % for Ralf lc_res distance on loop spectra
kfactor15 = 1.5; % for Ralf stud_res distance on loop spectra
kfactor16 = 1.5; % for Ralf Lof distance on loop spectra

x=session(sesid).sitabs(:,1:64);%this gets spectra

dnuid = [ ]; % no siiting exclusions   

% %dnuid = [6 7 16 18] % exclude for sesid=1 k=2 Mike Cook improves no xnew
% %dnuid = [6 7 9 12 16 23]  % exclude for sesid=1 k=2 Mike AbsDif improves morre no xnew
% 
% %dnuid = [1 5 8 11] % exclude for sesid=2 k=2 Mike Cook improves no xnew
% %dnuid = [1 5 11] xclude for sesid=2 k=2 Mike AbsDif improves moreno xnew
% 
% %dnuid = [1 2 3 12 26 28] % exclude for sesid=3 k=1 Pam Cook improves less no xnew
% %dnuid = [1 3 12 26 27 28 34] % exclude for sesid=3 k=1 Pam AbsDif improves more no xnew
% %dnuid = [1 3 12 26 31 36] % exclude for sesid=3 k=2 Pam Cook improves no xnew
% %dnuid = [1 3 12 21 26 27 31 36] % exclude for sesid=3 k=2 Pam AbsDif improves a bit more no xnew
% 
% %dnuid = [8 13 19] % exclude for sesid=4 k=2 new Jake Cook improves no xnew
% %dnuid = [7 13 17 19] % exclude for sesid=4 k=2 Jake AbsDif improves a bit more no xnew
% %dnuid = [1 2 3 19] % exclude for sesid=4 k=3 new Jake Cook improves no xnew
% %dnuid = [1 7 19] % exclude for sesid=4 k=3 Jake AbsDif a bit less no xnew
% 
% %dnuid = [3 4 16] % exclude for sesid=5 k=2 Anna Cook improves no xnew
% %dnuid = [2 4 12 16] % exclude for sesid=5 k=2 Anna AbsDIf makes worse no xnew
% %dnuid = [3 4 16] % exclude for sesid=5 k=2 Anna AbsDif improves a bit with xnew
% %dnuid = [3 10] % exclude for sesid=5 k=3 Anna Cook improves no xnew
% %dnuid = [2 4 8 10 16] % exclude for sesid=5 k=3 Anna AbsDIf makes worse no xnew********
% 
% %dnuid = [18 20 21 27 30 32] % exclude for sesid=6 k=1 Pam Cook improves a bit but high MARD no xnew
% %dnuid = [5 6 19 20 21 30 32] % exclude for sesid=6 k=1 Pam AbsDIf improves a lot but MARD > 15 no xnew
% %dnuid = [5 6 19 20 21 30 32] % exclude for sesid=6 k=1 Pam AbsDIf improves more but still high MARD with xnew
% %dnuid = [14 21 30] % exclude for sesid=6 k=2 Pam Cook improves but high MARD no xnew
% %dnuid = [5 19 20 21 27 30 32] % exclude for sesid=6 k=2 Pam AbsDif improves more but still high MARD no xnew and bit more with xnew
% %dnuid = [20 21 30 32] % exclude for sesid=6 k=3 Pam Cook improves some but high MARD no xnew
% %dnuid = [4 19 20 21 30 32] % exclude for sesid=6 k=3 Pam AbsDIf improves morebutstill high MARD  no xnew more with xnew but still high MARD
% %dnuid = [15 16 21 26 32] % exclude for sesid=6 k=3 Pam Cook improves some with xnew first 
% %dnuid = [3 18 20 21 26 30 31 32] % exclude for sesid=6 k=3 Pam AbsDiff improves more with xnew first
% %dnuid = [2 14 19 21 22 26 27] % exclude for sesid=6 k=3 Pam Mahal makes worse with xnew f
% 
% %dnuid = [3 8 14 16] % exclude for sesid=7 k=4 Mike Cook improves no xnew
% %dnuid = [3 8 9 16] % exclude for sesid=7 k=4 Mike AbsDif improves a lot no xnew
% %dnuid = [3 8 14 16] % exclude for sesid=7 k=4 Mike Cook improves with xnew first
% %dnuid = [3 8 9 11 16] % exclude for sesid=7 k=4 Mike AbsDif improves with xnew first
% %dnuid = [7 8 13] % exclude for sesid=7 k=5 Mike Cook improves just a bit no xnew
% %dnuid = [1 3 4 7 8] % exclude for sesid=7 k=5 Mike AbsDif improves a lot more no xnew
% 
% %dnuid = [8 10 11] % exclude for sesid=8 k=4 Deb Cook improves no xnew
% %dnuid = [1 11] % exclude for sesid=8 k=4 Deb AbsDif makes worse no xnew****
% %dnuid = [4 8 11] % exclude for sesid=8 k=5 Deb Cook makes worse no xnew 
% %dnuid = [1 8 11] % exclude for sesid=8 k=5 Deb AbsDif improves no xnew
% 
% %dnuid = [1 2] % exclude for ssesid=9 k=3 Mike Cook improves no xnew
% %dnuid = [2 4 7] % exclude for ssesid=9 k=3 Mike AbsDif slight improvement no xnew
% %dnuid = [1 5 8 11] % exclude for ssesid=9 k=4 Mike Cook improves no xnew
% %dnuid = [5 9 11] % exclude for ssesid=9 k=4 Mike AbsDif slight improvement no xnew
% %dnuid = [1 2 5] % exclude for ssesid=9 k=5 Mike Cook makes worse no xnew
% %dnuid = [1 4] % exclude for ssesid=9 k=5 Mike AbsDif improves a bit no xnew
% 
% %dnuid = [1 9 15 23 31 32 33 34] % exclude for sesid=10 k=1 Pam Cook improves no xnew
% %dnuid = [1 7 9 15 19 21 23 29 31 33 34] % exclude for sesid=10 k=1 Pam AbsDif improves no xnew
% %dnuid = [9 23 29 33 34] % exclude for sesid=10 k=4 Pam Cook improves no xnew
% %dnuid = [9 13 15 19 22 23 29 34] % exclude for sesid=10 k=4 Pam AbsDif improves less no xnew****
% %dnuid = [9 15 19 22 23 29 31 32 33 34] % exclude for sesid=10 k=5 Pam Cook makes worse no xnew
% %dnuid = [9 15 19 22 23 29 31 32 33 34] % exclude for sesid=10 k=5 Pam AbsDif improves no xnew
% 
% %dnuid = [5 6 9 14] % exclude for sesid=11 k=4 no xnew Deb Cook improves less no xnew
% %dnuid = [6 9 13 14] % exclude for sesid=11 k=4 no xnew Deb AbsDif improves more no xnew
% %dnuid = [5 6 14] % exclude for sesid=11 k=5 Deb Cook makes worse noxnew******
% %dnuid = [6 9 14] % exclude for sesid=11 k=5 Deb AbsDif makes worse no xnew *******
% 
% %dnuid = [3 13 18] % exclude for sesid=12 k=4 MikeCook improves a bit no xnew
% %dnuid = [2 3 10 13 18] % exclude for sesid=12 k=4 Mike AbsDif improves no xnew
% %dnuid = [3 6 18] % exclude for sesid=12 k=5 Mike Cook improves a bit no xnew
% %dnuid = [3 6 10] % exclude for sesid=12 k=5 Mike AbsDif improves more but MARD>15 no xnew
% 
% %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
g=session(sesid).glu;%this has the refmeter glucose values
session(sesid).refmeteruse; %this shows which ref meters are used [cgm, ultra, lite, oxi, verio]
st=session(sesid).spectimes;%these are the time staps for each spectrum
t=avp(st,6);%this averages the time stamps in chunks of 6

refmetercount=length(session(sesid).refmeteruse);
refmeterids=1:refmetercount;

for k=refmeterids(session(sesid).refmeteruse) %cycles through ref meters

sg=g(:,k);%start with reference values from the CGM
ax=x;
ax(:,69)=sg;%make the last col of the spectral array the glucose values
% clear loop Rlev Rcook
% loop=[];
% Rlev=[];
% RCook=[];

[loop,Rres_loop,Rlev_loop,Rlc_res_loop,Rstud_res_loop,Rcook_loop,RLof_loop]=...
    leaveoneoutpreds(ax(:,[1:28,35:64]),ax(:,69),1,ax(:,[1:28,35:64]),ax(:,69),1,1);
Rstud_res_loop_real= abs(Rstud_res_loop);
Rstud_res_loop_imag=imag(Rstud_res_loop);
% for i=1:size(ax,1)/6
%     % cid=logical(ones(size(ax,1),1));
%     % cid(i)=0;
%     cid=logical(ones(6,1)*ones(1,size(ax,1)/6));%make a set of logicals for all entries in the array
%     cid(:,i)=0;%set the ones corresponding to sitting i to zero so they're left out of the calibration
%     ccid=cid;
%     %***********************************
%     ccid(:,dnuid)=0; % do not use these sittings
%     %************************************
%     cid=cid(:);
%     ccid=ccid(:);
% 
%     %the code below can be used for the rectangular distribution of the
%     %glucose information
%     % [hr,rcal]=rba(ax(cid,1:64),ax(cid,end),500,50);
%     % [hr,rcal]=rba(ax(ccid,1:64),ax(ccid,end),1000,500);
%     % [hr,rcal]=rba(ax(ccid,1:64),ax(ccid,end),1000,50);
% 
%     %the code below performs a partial least squares regression and slope
%     %intercept corrects the bvector using the orthogonal regression slope
%     %of the self predicted results
%     % [mbr2bbcr,ysipreds,bvects]=jpredosic(ax(ccid,[1:28,35:64]),ax(ccid,69),1,ax(:,[1:28,35:64]),ax(:,69),1);
%     [~,ysipreds,~,lev,~,~,~,Cook,~]=jpredosicwplsout(ax(ccid,[1:28,35:64]),ax(ccid,69),1,ax(:,[1:28,35:64]),ax(:,69),1);
%     % [mbr2bbcr,ysipreds,bvects]=jpredosic(ax(ccid,[1:64]),ax(ccid,66),1,ax(:,[1:64]),ax(:,66),1);
%     % [mbr2bbcr,ysipreds,bvects]=jpredosic(rcal(:,[1:28,35:64]),hr,1,ax(:,[1:28,35:64]),ax(:,66),1);
%     % loop(i,:)=ysipreds(i,:);
%     loop=[loop;ysipreds(~cid,:)];%all of the spectra hag glucose values predicted, but save only the ones left out
%     Rlev=[Rlev;lev(~cid,:)];
%     RCook=[RCook;Cook(~cid,:)];   
% end %i
% 
% %figure
% %plot((t-t(1))*24*60,  (avp([loop(:,end),  loop(:,end-2)],  6)),'o:');grid
% %xlabel('Time')
% %ylabel('Blood Glucose Concentration (mg/dL)')
% %title(['Proto Measurements: ',session(sesid).user,', ',char(session(sesid).date),', ',char(session(sesid).refmeternames{k}),' as Reference Values'])
% %legend('Ref Meter Spline','ProtoData')
% %ylim([50 300]);
% 
% %end %k
% 
% w=avp([loop(:,end),  loop(:,end-2)],  6);
% bid = logical(ones(length(w),1)); bid(dnuid) = 0;  % exclude removed sittings
% figure; plot (w,'o:'); title(['Proto Measurements: ',session(sesid).user,', ',char(session(sesid).date),', ',char(session(sesid).refmeternames{k}),' as Reference Values']);grid
% 
% % Spectral data outlier detection using Ralf's functions*************************
% Ry= session(sesid).glu(:,k);
% % Rx=x;%Rx=x(:, 1:64);
% Rx=x(:,[1:28,35:64]);%Rx=x(:, 1:64);
% RuseR=57;
% [Rstat,Ru,Rb,Rv,Rxav,Ryav] = pls2(Rx,Ry,RuseR);
% RouteR = 20 %[1:28 35:64]; %ranks to be used
% [Rres,Rlev,lc_Rres,Rstud_Rres,RCook,Rlof,Rcrit] = plsout(Ru,Rb,Rv,Rxav,Ry,RouteR);
% Outlier_RCook = find(RCook>kfactor1*mean(RCook));
% %figure; plot(RCook); title('Ralf Cook distance on spectra')
% Outlier_RLev = find(Rlev>kfactor2*mean(Rlev));
% %figure; plot(Rlev); title('Ralf Leverage distance on spectra')
% %u = w(ecid,1); %ref
% 
% %*******************************************************************
% 
% % Jerry's ***********************************************************
% % Use this command line for loop to look at the predictions for each rank.
% % This helps to determine if all of the mathematically calculable ranks are needed (ie the min of m x n
% % of the cal matrix is the max rank that can be calculated).  Also, since the pls routine mean centers
% % he data, the last rank may be numerically deficient (ie, all of the relevant info is squeezed out of the data to below the bidiagonalization threshold).  If you see the warning message "bidiagonalization stop: Minute lu or lw at R = XX", rank XX and above may be useless.
% %  for j=1:57;plot((st-st(1))*24*60,(avp([loop(:,end),loop(:,j)],1)),'o:');grid;title(int2str(j));shg;pause;end % for individual sitting spectra
%  % for j=1:57;plot((t-t(1))*24*60,(avp([loop(:,end),loop(:,j)],6)),'o:');grid;title(int2str(j));shg;pause;end % for averaged sitting spectra
[Rres,Rlev,Rlc_res,Rstud_res,Rcook,RLof]=plsoutallranks(ax(:,[1:28,35:64]),ax(:,69));
% %******************************************************************
% 
% % Spectral data outlier detection using Quartiles with replacements in xnew
fillmethod = "spline";%"pchip"
findmethod ="quartiles";% "gesd" %"grubbs" %"median"%"mean"%
xnew = filloutliers(x,fillmethod,findmethod);
xdelta = abs(xnew - x);
%figure; plot(xdelta'); title('Quartiles distance on spectra')
Outlier_Quart = find((xdelta)>kfactor3*mean(xdelta));
xnew(:,65:69)=ax(:,65:69);

% 
% 
% % Paired data outlier detection
% u = w(:,1); %ref
% v = w(:,2); %proto
% dif = u - v;
% AD = abs(dif);
% 
% % Cook distance on pairs
% mdl = fitlm(v,u); %  LM = fitlm(X,Y) fits a linear regression model using Y a response variable
% % and X as predictor variable, and returns the linear model LM.
% cookd = mdl.Diagnostics.CooksDistance;
% Outlier_Cook = find((cookd)>kfactor4*mean(cookd));
% %figure; plotDiagnostics(mdl,'cookd');  title('Cook distance on paired')
% %figure; plot(cookd);  title('Cook distance on paired')
% 
% % Mahal distance on pairs
% mahad =  mahal(v,u); %
% Outlier_Maha =  find((mahad)>kfactor5*mean(mahad));
% %figure; plotDiagnostics(mdl,'cookd');  title('Mahal distance on paired')
% %figure; plot(mahad),  title('Mahal distance on paired')
% 
% % Absolute difference distance on pairs
% Outlier_AbsDif = find(AD>kfactor6*mean(AD));
% %figure; plot(AD); title('AbsDif distance on paired')
% 
% %Statistics on paired data
% ecid = logical(ones(length(w),1));  ecid(dnuid) = 0;  % exclude removed sittings
% u = w(ecid,1); %ref
% v = w(ecid,2); %proto
% dif = u - v;
% Stat_SEP1  = std(dif);
% Stat_MARD1 = 100*mean(abs(dif)./v);
% Tname = (['Proto Measurements: ',session(sesid).user,', ',char(session(sesid).date),', ',char(session(sesid).refmeternames{k}),' as Reference Values']);
% [STAT_Total, STAT_Clark_Percent] = clarkeEGC(u,v, Tname); %total points in each section %percentage of points in eah section

[cookd,mahad,AD,SMTP]=getPPstats(avp(loop,6));
% Outlier_AD             =  find((AD(:,end))           >kfactor6*mean(AD(:,end)));
% Outlier_Cook           =  find((cookd(:,end))        >kfactor4*mean(cookd(:,end)));
% Outlier_Maha           =  find((mahad(:,end))        >kfactor5*mean(mahad(:,end)));
% Outlier_Rres           =  find(Rres(:,end)           >kfactor7*mean(Rres(:,end)));
% Outlier_Rlev           =  find(Rlev(:,end)           >kfactor2*mean(Rlev(:,end)));
% Outlier_Rlc_res        =  find(Rlc_res(:,end)        >kfactor8*mean(Rlc_res(:,end)));
% Outlier_Rstud_res      =  find(Rstud_res(:,end)      >kfactor9*mean(Rstud_res(:,end)));
% Outlier_Rcook          =  find(Rcook(:,end)          >kfactor1*mean(Rcook(:,end)));
% Outlier_RLof           =  find(RLof(:,end)           >kfactor10*mean(RLof(:,end)));
% Outlier_Rres_loop      =  find(Rres_loop(:,end)      >kfactor13*mean(Rres_loop(:,end)));
% Outlier_Rlev_loop      =  find(Rlev_loop(:,end)      >kfactor12*mean(Rlev_loop(:,end)));
% Outlier_Rlc_res_loop   =  find(Rlc_res_loop(:,end)   >kfactor14*mean(Rlc_res_loop(:,end)));
% Outlier_Rstud_res_loop =  find(Rstud_res_loop_real(:,end) >kfactor15*mean(Rstud_res_loop_real(:,end)));
% Outlier_Rcook_loop     =  find(Rcook_loop(:,end)     >kfactor11*mean(Rcook_loop(:,end)));
% Outlier_RLof_loop      =  find(RLof_loop(:,end)      >kfactor16*mean(RLof_loop(:,end)));

% OutlierRank=57;
% OutlierRank=20;
OutlierRank=8;

Outlier_AD             =  find((AD(:,OutlierRank))                     >kfactor6*mean(AD(:,OutlierRank)));
Outlier_Cook           =  find((cookd(:,OutlierRank))                  >kfactor4*mean(cookd(:,OutlierRank)));
Outlier_Maha           =  find((mahad(:,OutlierRank))                  >kfactor5*mean(mahad(:,OutlierRank)));
Outlier_Rres           =  find(abs(Rres(:,OutlierRank))                >kfactor7*mean(abs(Rres(:,OutlierRank))));
Outlier_Rlev           =  find(Rlev(:,OutlierRank)                     >kfactor2*mean(Rlev(:,OutlierRank)));
Outlier_Rlc_res        =  find(abs(Rlc_res(:,OutlierRank))             >kfactor8*mean(abs(Rlc_res(:,OutlierRank))));
Outlier_Rstud_res      =  find(abs(Rstud_res(:,OutlierRank))           >kfactor9*mean(abs(Rstud_res(:,OutlierRank))));
Outlier_Rcook          =  find(Rcook(:,OutlierRank)                    >kfactor1*mean(Rcook(:,OutlierRank)));
Outlier_RLof           =  find(RLof(:,OutlierRank)                     >kfactor10*mean(RLof(:,OutlierRank)));
Outlier_Rres_loop      =  find(abs(Rres_loop(:,OutlierRank))           >kfactor13*mean(abs(Rres_loop(:,OutlierRank))));
Outlier_Rlev_loop      =  find(Rlev_loop(:,OutlierRank)                >kfactor12*mean(Rlev_loop(:,OutlierRank)));
Outlier_Rlc_res_loop   =  find(abs(Rlc_res_loop(:,OutlierRank))        >kfactor14*mean(abs(Rlc_res_loop(:,OutlierRank))));
Outlier_Rstud_res_loop =  find(abs(Rstud_res_loop_real(:,OutlierRank)) >kfactor15*mean(abs(Rstud_res_loop_real(:,OutlierRank))));
Outlier_Rcook_loop     =  find(Rcook_loop(:,OutlierRank)               >kfactor11*mean(Rcook_loop(:,OutlierRank)));
Outlier_RLof_loop      =  find(RLof_loop(:,OutlierRank)                >kfactor16*mean(RLof_loop(:,OutlierRank)));
% did1  = true(size(x,1),1); did1(Outlier_Quart) = 0; xnew1=ax; xnew1(~did1,:)=nan; %sg=g(did1,k);%start with reference values from the CGM  % Use to remove outliers by Quartiles
did2  = true(size(x,1),1); did2(Outlier_Rcook)           = 0; xnew2 =ax; xnew2(~did2,:)  =nan; %sg=g(did2,k);%start with reference values from the CGM % Use to remove outliers by Cook
did3  = true(size(x,1),1); did3(Outlier_Rlev)            = 0; xnew3 =ax; xnew3(~did3,:)  =nan; %sg=g(did3,k);%start with reference values from the CGM % Use to remove outliers by Leverag
did4  = true(size(x,1),1); did4(Outlier_AD)              = 0; xnew4 =ax; xnew4(~did4,:)  =nan; %sg=g(did3,k);%start with reference values from the CGM % Use to remove outliers by Leverag
did5  = true(size(x,1),1); did5(Outlier_Cook)            = 0; xnew5 =ax; xnew5(~did5,:)  =nan; %sg=g(did3,k);%start with reference values from the CGM % Use to remove outliers by Leverag
did6  = true(size(x,1),1); did6(Outlier_Maha)            = 0; xnew6 =ax; xnew6(~did6,:)  =nan; %sg=g(did3,k);%start with reference values from the CGM % Use to remove outliers by Leverag
did7  = true(size(x,1),1); did7(Outlier_Rres)            = 0; xnew7 =ax; xnew7(~did7,:)  =nan; %sg=g(did3,k);%start with reference values from the CGM % Use to remove outliers by Leverag
did8  = true(size(x,1),1); did8(Outlier_Rlc_res)         = 0; xnew8 =ax; xnew8(~did8,:)  =nan; %sg=g(did3,k);%start with reference values from the CGM % Use to remove outliers by Leverag
did9  = true(size(x,1),1); did9(Outlier_Rstud_res)       = 0; xnew9 =ax; xnew9(~did9,:)  =nan; %sg=g(did3,k);%start with reference values from the CGM % Use to remove outliers by Leverag
did10 = true(size(x,1),1); did10(Outlier_RLof)           = 0; xnew10=ax; xnew10(~did10,:)=nan; %sg=g(did3,k);%start with reference values from the CGM % Use to remove outliers by Leverag
did11 = true(size(x,1),1); did11(Outlier_Rcook_loop)     = 0; xnew11=ax; xnew11(~did11,:)=nan; %sg=g(did3,k);%start with reference values from the CGM % Use to remove outliers by Leverag
did12 = true(size(x,1),1); did12(Outlier_Rlc_res_loop)   = 0; xnew12=ax; xnew12(~did12,:)=nan; %sg=g(did3,k);%start with reference values from the CGM % Use to remove outliers by Leverag
did13 = true(size(x,1),1); did13(Outlier_Rstud_res_loop) = 0; xnew13=ax; xnew13(~did13,:)=nan; %sg=g(did3,k);%start with reference values from the CGM % Use to remove outliers by Leverag
did14 = true(size(x,1),1); did14(Outlier_Rlev_loop)      = 0; xnew14=ax; xnew14(~did14,:)=nan; %sg=g(did3,k);%start with reference values from the CGM % Use to remove outliers by Leverag
did15 = true(size(x,1),1); did15(Outlier_RLof_loop)      = 0; xnew15=ax; xnew15(~did15,:)=nan; %sg=g(did3,k);%start with reference values from the CGM % Use to remove outliers by Leverag
did16 = true(size(x,1),1); did16(Outlier_Rres_loop)      = 0; xnew16=ax; xnew16(~did16,:)=nan; %sg=g(did3,k);%start with reference values from the CGM % Use to remove outliers by Leverag

    [Qres       , T2]        = getQT(ax   (:    ,[1:28,35:64]) ,ax   (:    ,69));
    [Qres_xnew  , T2_xnew]   = getQT(xnew (:    ,[1:28,35:64]) ,xnew (:    ,69));
    [Qres_xnew2 , T2_xnew2]  = getQT(ax   (did2 ,[1:28,35:64]) ,ax   (did2 ,69));
    [Qres_xnew3 , T2_xnew3]  = getQT(ax   (did3 ,[1:28,35:64]) ,ax   (did3 ,69));
    [Qres_xnew4 , T2_xnew4]  = getQT(ax   (did4 ,[1:28,35:64]) ,ax   (did4 ,69));
    [Qres_xnew5 , T2_xnew5]  = getQT(ax   (did5 ,[1:28,35:64]) ,ax   (did5 ,69));
    [Qres_xnew6 , T2_xnew6]  = getQT(ax   (did6 ,[1:28,35:64]) ,ax   (did6 ,69));
    [Qres_xnew7 , T2_xnew7]  = getQT(ax   (did7 ,[1:28,35:64]) ,ax   (did7 ,69));
    [Qres_xnew8 , T2_xnew8]  = getQT(ax   (did8 ,[1:28,35:64]) ,ax   (did8 ,69));
    [Qres_xnew9 , T2_xnew9]  = getQT(ax   (did9 ,[1:28,35:64]) ,ax   (did9 ,69));
    [Qres_xnew10, T2_xnew10] = getQT(ax   (did10,[1:28,35:64]) ,ax   (did10,69));
    [Qres_xnew11, T2_xnew11] = getQT(ax   (did11,[1:28,35:64]) ,ax   (did11,69));
    [Qres_xnew12, T2_xnew12] = getQT(ax   (did12,[1:28,35:64]) ,ax   (did12,69));
    [Qres_xnew13, T2_xnew13] = getQT(ax   (did13,[1:28,35:64]) ,ax   (did13,69));
    [Qres_xnew14, T2_xnew14] = getQT(ax   (did14,[1:28,35:64]) ,ax   (did14,69));
    [Qres_xnew15, T2_xnew15] = getQT(ax   (did15,[1:28,35:64]) ,ax   (did15,69));
    [Qres_xnew16, T2_xnew16] = getQT(ax   (did16,[1:28,35:64]) ,ax   (did16,69));


[loop_xnew ] =leaveoneoutpreds(xnew (: ,[1:28,35:64]),xnew (: ,69),1,ax(:,[1:28,35:64]),ax(:,69),1,0);
% [loop_xnew1] =leaveoneoutpreds(xnew1(:,[1:28,35:64]),xnew1(:,69),1,ax(:,[1:28,35:64]),ax(:,69),1,0);
[loop_xnew2] =leaveoneoutpreds(xnew2(: ,[1:28,35:64]),xnew2(: ,69),1,ax(:,[1:28,35:64]),ax(:,69),1,0);
[loop_xnew3] =leaveoneoutpreds(xnew3(: ,[1:28,35:64]),xnew3(: ,69),1,ax(:,[1:28,35:64]),ax(:,69),1,0);
[loop_xnew4] =leaveoneoutpreds(xnew4(: ,[1:28,35:64]),xnew4(: ,69),1,ax(:,[1:28,35:64]),ax(:,69),1,0);
[loop_xnew5] =leaveoneoutpreds(xnew5(: ,[1:28,35:64]),xnew5(: ,69),1,ax(:,[1:28,35:64]),ax(:,69),1,0);
[loop_xnew6] =leaveoneoutpreds(xnew6(: ,[1:28,35:64]),xnew6(: ,69),1,ax(:,[1:28,35:64]),ax(:,69),1,0);
[loop_xnew7] =leaveoneoutpreds(xnew7(: ,[1:28,35:64]),xnew7(: ,69),1,ax(:,[1:28,35:64]),ax(:,69),1,0);
[loop_xnew8] =leaveoneoutpreds(xnew8(: ,[1:28,35:64]),xnew8(: ,69),1,ax(:,[1:28,35:64]),ax(:,69),1,0);
[loop_xnew9] =leaveoneoutpreds(xnew9(: ,[1:28,35:64]),xnew9(: ,69),1,ax(:,[1:28,35:64]),ax(:,69),1,0);
[loop_xnew10]=leaveoneoutpreds(xnew10(:,[1:28,35:64]),xnew10(:,69),1,ax(:,[1:28,35:64]),ax(:,69),1,0);
[loop_xnew11]=leaveoneoutpreds(xnew11(:,[1:28,35:64]),xnew11(:,69),1,ax(:,[1:28,35:64]),ax(:,69),1,0);
[loop_xnew12]=leaveoneoutpreds(xnew12(:,[1:28,35:64]),xnew12(:,69),1,ax(:,[1:28,35:64]),ax(:,69),1,0);
[loop_xnew13]=leaveoneoutpreds(xnew13(:,[1:28,35:64]),xnew13(:,69),1,ax(:,[1:28,35:64]),ax(:,69),1,0);
[loop_xnew14]=leaveoneoutpreds(xnew14(:,[1:28,35:64]),xnew14(:,69),1,ax(:,[1:28,35:64]),ax(:,69),1,0);
[loop_xnew15]=leaveoneoutpreds(xnew15(:,[1:28,35:64]),xnew15(:,69),1,ax(:,[1:28,35:64]),ax(:,69),1,0);
[loop_xnew16]=leaveoneoutpreds(xnew16(:,[1:28,35:64]),xnew16(:,69),1,ax(:,[1:28,35:64]),ax(:,69),1,0);

[cookd_xnew  ,mahad_xnew   ,AD_xnew ,SMTP_xnew]  =getPPstats(avp(loop_xnew  ,6));
[cookd_xnew2 ,mahad_xnew2  ,AD_xnew2,SMTP_xnew2] =getPPstats(avp(loop_xnew2 ,6));
[cookd_xnew3 ,mahad_xnew3  ,AD_xnew3,SMTP_xnew3] =getPPstats(avp(loop_xnew3 ,6));
[cookd_xnew4 ,mahad_xnew4  ,AD_xnew4,SMTP_xnew4] =getPPstats(avp(loop_xnew4 ,6));
[cookd_xnew5 ,mahad_xnew5  ,AD_xnew5,SMTP_xnew5] =getPPstats(avp(loop_xnew5 ,6));
[cookd_xnew6 ,mahad_xnew6  ,AD_xnew6,SMTP_xnew6] =getPPstats(avp(loop_xnew6 ,6));
[cookd_xnew7 ,mahad_xnew7  ,AD_xnew7,SMTP_xnew7] =getPPstats(avp(loop_xnew7 ,6));
[cookd_xnew8 ,mahad_xnew8  ,AD_xnew8,SMTP_xnew8] =getPPstats(avp(loop_xnew8 ,6));
[cookd_xnew9 ,mahad_xnew9  ,AD_xnew9,SMTP_xnew9] =getPPstats(avp(loop_xnew9 ,6));
[cookd_xnew10,mahad_xnew10,AD_xnew10,SMTP_xnew10]=getPPstats(avp(loop_xnew10,6));
[cookd_xnew11,mahad_xnew11,AD_xnew11,SMTP_xnew11]=getPPstats(avp(loop_xnew11,6));
[cookd_xnew12,mahad_xnew12,AD_xnew12,SMTP_xnew12]=getPPstats(avp(loop_xnew12,6));
[cookd_xnew13,mahad_xnew13,AD_xnew13,SMTP_xnew13]=getPPstats(avp(loop_xnew13,6));
[cookd_xnew14,mahad_xnew14,AD_xnew14,SMTP_xnew14]=getPPstats(avp(loop_xnew14,6));
[cookd_xnew15,mahad_xnew15,AD_xnew15,SMTP_xnew15]=getPPstats(avp(loop_xnew15,6));
[cookd_xnew16,mahad_xnew16,AD_xnew16,SMTP_xnew16]=getPPstats(avp(loop_xnew16,6));


outlier(sesid,k).loop        =loop;
outlier(sesid,k).loop_xnew   =loop_xnew;
outlier(sesid,k).loop_xnew2  =loop_xnew2;
outlier(sesid,k).loop_xnew3  =loop_xnew3;
outlier(sesid,k).loop_xnew4  =loop_xnew4;
outlier(sesid,k).loop_xnew5  =loop_xnew5;
outlier(sesid,k).loop_xnew6  =loop_xnew6;
outlier(sesid,k).loop_xnew7  =loop_xnew7;
outlier(sesid,k).loop_xnew8  =loop_xnew8;
outlier(sesid,k).loop_xnew9  =loop_xnew9;
outlier(sesid,k).loop_xnew10 =loop_xnew10;
outlier(sesid,k).loop_xnew11 =loop_xnew11;
outlier(sesid,k).loop_xnew12 =loop_xnew12;
outlier(sesid,k).loop_xnew13 =loop_xnew13;
outlier(sesid,k).loop_xnew14 =loop_xnew14;
outlier(sesid,k).loop_xnew15 =loop_xnew15;
outlier(sesid,k).loop_xnew16 =loop_xnew16;
outlier(sesid,k).did2  =did2;
outlier(sesid,k).did3  =did3;
outlier(sesid,k).did4  =did4;
outlier(sesid,k).did5  =did5;
outlier(sesid,k).did6  =did6;
outlier(sesid,k).did7  =did7;
outlier(sesid,k).did8  =did8;
outlier(sesid,k).did9  =did9;
outlier(sesid,k).did10 =did10;
outlier(sesid,k).did11 =did11;
outlier(sesid,k).did12 =did12;
outlier(sesid,k).did13 =did13;
outlier(sesid,k).did14 =did14;
outlier(sesid,k).did15 =did15;
outlier(sesid,k).did16 =did16;
outlier(sesid,k).cookd        =cookd;
outlier(sesid,k).cookd_xnew   =cookd_xnew;
outlier(sesid,k).cookd_xnew2  =cookd_xnew2;
outlier(sesid,k).cookd_xnew3  =cookd_xnew3;
outlier(sesid,k).cookd_xnew4  =cookd_xnew4;
outlier(sesid,k).cookd_xnew5  =cookd_xnew5;
outlier(sesid,k).cookd_xnew6  =cookd_xnew6;
outlier(sesid,k).cookd_xnew7  =cookd_xnew7;
outlier(sesid,k).cookd_xnew8  =cookd_xnew8;
outlier(sesid,k).cookd_xnew9  =cookd_xnew9;
outlier(sesid,k).cookd_xnew10 =cookd_xnew10;
outlier(sesid,k).cookd_xnew11 =cookd_xnew11;
outlier(sesid,k).cookd_xnew12 =cookd_xnew12;
outlier(sesid,k).cookd_xnew13 =cookd_xnew13;
outlier(sesid,k).cookd_xnew14 =cookd_xnew14;
outlier(sesid,k).cookd_xnew15 =cookd_xnew15;
outlier(sesid,k).cookd_xnew16 =cookd_xnew16;
outlier(sesid,k).mahad        =mahad;
outlier(sesid,k).mahad_xnew   =mahad_xnew;
outlier(sesid,k).mahad_xnew2  =mahad_xnew2;
outlier(sesid,k).mahad_xnew3  =mahad_xnew3;
outlier(sesid,k).mahad_xnew4  =mahad_xnew4;
outlier(sesid,k).mahad_xnew5  =mahad_xnew5;
outlier(sesid,k).mahad_xnew6  =mahad_xnew6;
outlier(sesid,k).mahad_xnew7  =mahad_xnew7;
outlier(sesid,k).mahad_xnew8  =mahad_xnew8;
outlier(sesid,k).mahad_xnew9  =mahad_xnew9;
outlier(sesid,k).mahad_xnew10 =mahad_xnew10;
outlier(sesid,k).mahad_xnew11 =mahad_xnew11;
outlier(sesid,k).mahad_xnew12 =mahad_xnew12;
outlier(sesid,k).mahad_xnew13 =mahad_xnew13;
outlier(sesid,k).mahad_xnew14 =mahad_xnew14;
outlier(sesid,k).mahad_xnew15 =mahad_xnew15;
outlier(sesid,k).mahad_xnew16 =mahad_xnew16;
outlier(sesid,k).AD        =AD;
outlier(sesid,k).AD_xnew   =AD_xnew;
outlier(sesid,k).AD_xnew2  =AD_xnew2;
outlier(sesid,k).AD_xnew3  =AD_xnew3;
outlier(sesid,k).AD_xnew4  =AD_xnew4;
outlier(sesid,k).AD_xnew5  =AD_xnew5;
outlier(sesid,k).AD_xnew6  =AD_xnew6;
outlier(sesid,k).AD_xnew7  =AD_xnew7;
outlier(sesid,k).AD_xnew8  =AD_xnew8;
outlier(sesid,k).AD_xnew9  =AD_xnew9;
outlier(sesid,k).AD_xnew10 =AD_xnew10;
outlier(sesid,k).AD_xnew11 =AD_xnew11;
outlier(sesid,k).AD_xnew12 =AD_xnew12;
outlier(sesid,k).AD_xnew13 =AD_xnew13;
outlier(sesid,k).AD_xnew14 =AD_xnew14;
outlier(sesid,k).AD_xnew15 =AD_xnew15;
outlier(sesid,k).AD_xnew16 =AD_xnew16;
outlier(sesid,k).SMTP        =SMTP;
outlier(sesid,k).SMTP_xnew   =SMTP_xnew;
outlier(sesid,k).SMTP_xnew2  =SMTP_xnew2;
outlier(sesid,k).SMTP_xnew3  =SMTP_xnew3;
outlier(sesid,k).SMTP_xnew4  =SMTP_xnew4;
outlier(sesid,k).SMTP_xnew5  =SMTP_xnew5;
outlier(sesid,k).SMTP_xnew6  =SMTP_xnew6;
outlier(sesid,k).SMTP_xnew7  =SMTP_xnew7;
outlier(sesid,k).SMTP_xnew8  =SMTP_xnew8;
outlier(sesid,k).SMTP_xnew9  =SMTP_xnew9;
outlier(sesid,k).SMTP_xnew10 =SMTP_xnew10;
outlier(sesid,k).SMTP_xnew11 =SMTP_xnew11;
outlier(sesid,k).SMTP_xnew12 =SMTP_xnew12;
outlier(sesid,k).SMTP_xnew13 =SMTP_xnew13;
outlier(sesid,k).SMTP_xnew14 =SMTP_xnew14;
outlier(sesid,k).SMTP_xnew15 =SMTP_xnew15;
outlier(sesid,k).SMTP_xnew16 =SMTP_xnew16;
outlier(sesid,k).Qres        =Qres;
outlier(sesid,k).Qres_xnew   =Qres_xnew;
outlier(sesid,k).Qres_xnew2  =Qres_xnew2;
outlier(sesid,k).Qres_xnew3  =Qres_xnew3;
outlier(sesid,k).Qres_xnew4  =Qres_xnew4;
outlier(sesid,k).Qres_xnew5  =Qres_xnew5;
outlier(sesid,k).Qres_xnew6  =Qres_xnew6;
outlier(sesid,k).Qres_xnew7  =Qres_xnew7;
outlier(sesid,k).Qres_xnew8  =Qres_xnew8;
outlier(sesid,k).Qres_xnew9  =Qres_xnew9;
outlier(sesid,k).Qres_xnew10 =Qres_xnew10;
outlier(sesid,k).Qres_xnew11 =Qres_xnew11;
outlier(sesid,k).Qres_xnew12 =Qres_xnew12;
outlier(sesid,k).Qres_xnew13 =Qres_xnew13;
outlier(sesid,k).Qres_xnew14 =Qres_xnew14;
outlier(sesid,k).Qres_xnew15 =Qres_xnew15;
outlier(sesid,k).Qres_xnew16 =Qres_xnew16;
outlier(sesid,k).T2        =T2;
outlier(sesid,k).T2_xnew   =T2_xnew;
outlier(sesid,k).T2_xnew2  =T2_xnew2;
outlier(sesid,k).T2_xnew3  =T2_xnew3;
outlier(sesid,k).T2_xnew4  =T2_xnew4;
outlier(sesid,k).T2_xnew5  =T2_xnew5;
outlier(sesid,k).T2_xnew6  =T2_xnew6;
outlier(sesid,k).T2_xnew7  =T2_xnew7;
outlier(sesid,k).T2_xnew8  =T2_xnew8;
outlier(sesid,k).T2_xnew9  =T2_xnew9;
outlier(sesid,k).T2_xnew10 =T2_xnew10;
outlier(sesid,k).T2_xnew11 =T2_xnew11;
outlier(sesid,k).T2_xnew12 =T2_xnew12;
outlier(sesid,k).T2_xnew13 =T2_xnew13;
outlier(sesid,k).T2_xnew14 =T2_xnew14;
outlier(sesid,k).T2_xnew15 =T2_xnew15;
outlier(sesid,k).T2_xnew16 =T2_xnew16;
outlier(sesid,k).Outlier_AD=Outlier_AD;
outlier(sesid,k).Outlier_Cook=Outlier_Cook;
outlier(sesid,k).Outlier_Maha=Outlier_Maha;
outlier(sesid,k).Outlier_Rres=Outlier_Rres;
outlier(sesid,k).Outlier_Rlev=Outlier_Rlev;
outlier(sesid,k).Outlier_Rlc_res=Outlier_Rlc_res;
outlier(sesid,k).Outlier_Rstud_res=Outlier_Rstud_res;
outlier(sesid,k).Outlier_Rcook=Outlier_Rcook;
outlier(sesid,k).Outlier_RLof=Outlier_RLof;
outlier(sesid,k).Outlier_Rres_loop=Outlier_Rres_loop;
outlier(sesid,k).Outlier_Rlev_loop=Outlier_Rlev_loop;
outlier(sesid,k).Outlier_Rlc_res_loop=Outlier_Rlc_res_loop;
outlier(sesid,k).Outlier_Rstud_res_loop=Outlier_Rstud_res_loop;
outlier(sesid,k).Outlier_Rcook_loop=Outlier_Rcook_loop;
outlier(sesid,k).Outlier_RLof_loop=Outlier_RLof_loop;
outlier(sesid,k).Outlier_Quart=Outlier_Quart;
outlier(sesid,k).Rres=Rres;
outlier(sesid,k).Rlev=Rlev;
outlier(sesid,k).Rlc_res=Rlc_res;
outlier(sesid,k).Rstud_res=Rstud_res;
outlier(sesid,k).Rcook=Rcook;
outlier(sesid,k).RLof=RLof;
outlier(sesid,k).Rres_loop=Rres_loop;
outlier(sesid,k).Rlev_loop=Rlev_loop;
outlier(sesid,k).Rlc_res_loop=Rlc_res_loop;
outlier(sesid,k).Rstud_res_loop=Rstud_res_loop;
outlier(sesid,k).Rcook_loop=Rcook_loop;
outlier(sesid,k).RLof_loop=RLof_loop;
outlier(sesid,k).kfactor1  =kfactor1;
outlier(sesid,k).kfactor2  =kfactor2;
outlier(sesid,k).kfactor3  =kfactor3;
outlier(sesid,k).kfactor4  =kfactor4;
outlier(sesid,k).kfactor5  =kfactor5;
outlier(sesid,k).kfactor6  =kfactor6;
outlier(sesid,k).kfactor7  =kfactor7;
outlier(sesid,k).kfactor8  =kfactor8;
outlier(sesid,k).kfactor9  =kfactor9;
outlier(sesid,k).kfactor10 =kfactor10;
outlier(sesid,k).kfactor11 =kfactor11;
outlier(sesid,k).kfactor12 =kfactor12;
outlier(sesid,k).kfactor13 =kfactor13;
outlier(sesid,k).kfactor14 =kfactor14;
outlier(sesid,k).kfactor15 =kfactor15;
outlier(sesid,k).kfactor16 =kfactor16;

disp([sesid,k,toc])

end
end
% save outlieronallsessions outlier
% save outlieronallsessionsOutlierRank20 outlier
save outlieronallsessionsOutlierRank8 outlier

% Enable a line to check effect of removing outliers (cmpare MARD2 to MARD1)
%x=xnew; % xnew for replaced specta by outliers in pairs
%x=xnew1; % xnew for removed specta by Quartiles RESULTS TO SAME MARD
%x=xnew2; % xnew for removed specta by Cook GIVES ERROR
%x=xnew3; % xnew for removed specta by Leverage GIVES ERROR

%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
% g=session(sesid).glu;%this has the refmeter glucose values
% session(sesid).refmeteruse; %this shows which ref meters are used [cgm, ultra, lite, oxi, verio]
% st=session(sesid).spectimes;%these are the time staps for each spectrum
% t=avp(st,6);%this averages the time stamps in chunks of 6
% 
% refmetercount=length(session(sesid).refmeteruse);
% refmeterids=1:refmetercount;
% 
% %for k=refmeterids(session(sesid).refmeteruse) %cycles through ref meters
% 
% sg=g(:,k);%start with reference values from the CGM
% ax=x;
% ax(:,69)=sg;%make the last col of the spectral array the glucose values
% clear loop
% loop=[];
% 
% for i=1:size(ax,1)/6
%     % cid=logical(ones(size(ax,1),1));
%     % cid(i)=0;
%     cid=logical(ones(6,1)*ones(1,size(ax,1)/6));%make a set of logicals for all entries in the array
%     cid(:,i)=0;%set the ones corresponding to sitting i to zero so they're left out of the calibration
%     ccid=cid;
%     %***********************************
%     ccid(:,dnuid)=0; % do not use these sittings
%     %************************************
%     cid=cid(:);
%     ccid=ccid(:);
% 
%     %the code below can be used for the rectangular distribution of the
%     %glucose information
%     % [hr,rcal]=rba(ax(cid,1:64),ax(cid,end),500,50);
%     % [hr,rcal]=rba(ax(ccid,1:64),ax(ccid,end),1000,500);
%     % [hr,rcal]=rba(ax(ccid,1:64),ax(ccid,end),1000,50);
% 
%     %the code below performs a partial least squares regression and slope
%     %intercept corrects the bvector using the orthogonal regression slope
%     %of the self predicted results
%     [mbr2bbcr,ysipreds,bvects]=jpredosic(ax(ccid,[1:28,35:64]),ax(ccid,69),1,ax(:,[1:28,35:64]),ax(:,69),1);
%     % [mbr2bbcr,ysipreds,bvects]=jpredosic(ax(ccid,[1:64]),ax(ccid,66),1,ax(:,[1:64]),ax(:,66),1);
%     % [mbr2bbcr,ysipreds,bvects]=jpredosic(rcal(:,[1:28,35:64]),hr,1,ax(:,[1:28,35:64]),ax(:,66),1);
%     % loop(i,:)=ysipreds(i,:);
%     loop=[loop;ysipreds(~cid,:)];%all of the spectra hag glucose values predicted, but save only the ones left out
% end %i
% 
% %figure
% %plot((t-t(1))*24*60,  (avp([loop(:,end),  loop(:,end-2)],  6)),'o:');grid
% %xlabel('Time')
% %ylabel('Blood Glucose Concentration (mg/dL)')
% %title(['Proto Measurements: ',session(sesid).user,', ',char(session(sesid).date),', ',char(session(sesid).refmeternames{k}),' as Reference Values'])
% %legend('Ref Meter Spline','ProtoData')
% %ylim([50 300]);
% 
% %end %k
% 
% w=avp([loop(:,end),  loop(:,end-2)],  6);
% bid = logical(ones(length(w),1)); bid(dnuid) = 0;  % exclude removed sittings
% figure; plot (w,'o:'); title(['Proto Measurements: ',session(sesid).user,', ',char(session(sesid).date),', ',char(session(sesid).refmeternames{k}),' as Reference Values']);grid
% 
% % Spectral data outlier detection using Ralf's functions*************************
% Ry= session(sesid).glu(:,k);
% Rx=x;%Rx=x(:, 1:64);
% RuseR=58;
% [Rstat,Ru,Rb,Rv,Rxav,Ryav] = pls2(Rx,Ry,RuseR);
% RouteR = 20 %[1:28 35:64]; %ranks to be used
% [Rres,Rlev,lc_Rres,Rstud_Rres,RCook,Rlof,Rcrit] = plsout(Ru,Rb,Rv,Rxav,Ry,RouteR);
% Outlier_RCook = find(RCook>kfactor1*mean(RCook));
% %figure; plot(RCook); title('Ralf Cook distance on spectra')
% Outlier_RLev = find(Rlev>kfactor2*mean(Rlev));
% %figure; plot(Rlev); title('Ralf Leverage distance on spectra')
% %u = w(ecid,1); %ref
% 
% %*******************************************************************
% 
% % Jerry's ***********************************************************
% % Use this command line for loop to look at the predictions for each rank.
% % This helps to determine if all of the mathematically calculable ranks are needed (ie the min of m x n
% % of the cal matrix is the max rank that can be calculated).  Also, since the pls routine mean centers
% % he data, the last rank may be numerically deficient (ie, all of the relevant info is squeezed out of the data to below the bidiagonalization threshold).  If you see the warning message "bidiagonalization stop: Minute lu or lw at R = XX", rank XX and above may be useless.
% %  for j=1:57;plot((st-st(1))*24*60,(avp([loop(:,end),loop(:,j)],1)),'o:');grid;title(int2str(j));shg;pause;end % for individual sitting spectra
%  % for j=1:57;plot((t-t(1))*24*60,(avp([loop(:,end),loop(:,j)],6)),'o:');grid;title(int2str(j));shg;pause;end % for averaged sitting spectra
% %******************************************************************
% 
% % Spectral data outlier detection using Quartiles with replacements in xnew
% fillmethod = "spline";%"pchip"
% findmethod ="quartiles";% "gesd" %"grubbs" %"median"%"mean"%
% xnew = filloutliers(x,fillmethod,findmethod);
% xdelta = abs(xnew - x);
% %figure; plot(xdelta'); title('Quartiles distance on spectra')
% Outlier_Quart = find((xdelta)>kfactor3*mean(xdelta));
% 
% 
% % Paired data outlier detection
% u = w(:,1); %ref
% v = w(:,2); %proto
% dif = u - v;
% AD = abs(dif);
% 
% % Cook distance on pairs
% mdl = fitlm(v,u); %  LM = fitlm(X,Y) fits a linear regression model using Y a response variable
% % and X as predictor variable, and returns the linear model LM.
% cookd = mdl.Diagnostics.CooksDistance;
% Outlier_Cook = find((cookd)>kfactor4*mean(cookd));
% %figure; plotDiagnostics(mdl,'cookd');  title('Cook distance on paired')
% %figure; plot(cookd);  title('Cook distance on paired')
% 
% % Mahal distance on pairs
% mahad =  mahal(v,u); %
% Outlier_Maha =  find((mahad)>kfactor5*mean(mahad));
% %figure; plotDiagnostics(mdl,'cookd');  title('Mahal distance on paired')
% %figure; plot(mahad),  title('Mahal distance on paired')
% 
% % Absolute difference distance on pairs
% Outlier_AbsDif = find(AD>kfactor6*mean(AD));
% %figure; plot(AD); title('AbsDif distance on paired')
% 
% %Statistics on paired data
% ecid = logical(ones(length(w),1));  ecid(dnuid) = 0;  % exclude removed sittings
% u = w(ecid,1); %ref
% v = w(ecid,2); %proto
% dif = u - v;
% Stat_SEP2  = std(dif);
% Stat_MARD2 = 100*mean(abs(dif)./v);
% Tname = (['Proto Measurements: ',session(sesid).user,', ',char(session(sesid).date),', ',char(session(sesid).refmeternames{k}),' as Reference Values']);
% [STAT_Total, STAT_Clark_Percent] = clarkeEGC(u,v, Tname); %total points in each section %percentage of points in eah section
