

% Test case 1. Martens data: 
load GlutenStarchNIR
Yall=Matrix';
XX = Yall(13:112,:);
YY = Yall(4,:);

aa=0*a;

na=0;
als=[];
apcr=[];
apls=[];
aloss=[];

% repeated Martens
% 100 times where I take out data each time
for nsis=1:98
X = XX(:,[1:nsis nsis+2:100]); X0 = XX(:,nsis+1);
Y1 = YY(:,[1:nsis nsis+2:100]); Y10 = YY(:,nsis+1);
% Two end cases are handled separately
nsis=0
X = XX(:,[nsis+2:100]); X0 = XX(:,nsis+1);
Y1 = YY(:,[nsis+2:100]); Y10 = YY(:,nsis+1);
nsis=99
X = XX(:,[1:nsis]); X0 = XX(:,nsis+1);
Y1 = YY(:,[1:nsis]); Y10 = YY(:,nsis+1)

% MATLAB code starts here

% Least squares (= PCR with default tol) 
method1='   LS';
Hls = Y1*pinv(X);
res=Hls*X;
res0=Hls*X0;
a11=norm(res-Y1);        
a12=norm(res0-Y10); 
als=[als a12];

% PCR (vary  tol)
% PCR: To find cut-off to find no. of components, semilogy(svd(X))
method6='   PCR';
tol=5.e-4
Hpcr = Y1*pinv(X,tol);
npcr=rank(pinv(X,tol))   % no. of components used in PCR
res=Hpcr*X;
res0=Hpcr*X0;
a61=norm(res-Y1);        
a62=norm(res0-Y10); 
apcr=[apcr a62];

% pls
npls=7
%npls=npcr
[XL,yl,XS,YS,beta,PCTVAR] = plsregress(X',Y1',npls);
% Note that PLS has an additional bias/centering parameter. 
yfit = [ones(size(X',1),1) X']*beta;
yfit0 = [ones(size(X0',1),1) X0']*beta;
a71=norm(yfit-Y1');        
a72=norm(yfit0-Y10'); 
method7='    PLS';
apls=[apls a72];

% Loss method  
%  step 1
Yall = [Y1; X];
[u,s,v]=svd(Y1);
Yall1 = Yall*v;
[r1,c1]=size(Yall); 
[r2,c2]=size(Y1); 
ny=r2;
G1=Yall1(1:ny,1:ny);
Gy=Yall1(ny+1:r1,1:ny);
Xopt = Yall1(ny+1:r1,ny+1:c1); 

% loss method: add artificial noise weights 
noisemag=5.e-4
%noisemag=tol
[rx,cx]=size(X);
noise = noisemag*eye(rx);
Xoptnoise = [Xopt noise];
% Noise Alt2. Add noise proportional to variation in each output
%[rx,cx]=size(Xopt); 
%mag=[];
%for i = 1:rx
%mag=[mag norm(Xopt(i,:))]; % the magnitudes are about 0.8
%end
%noise = noisemag*diag(mag);
%Xoptnoise = [Xopt noise];

% step 2- with artificial noise
% semilogy(svd(X)) 
% semilogy(svd([X X0]) 
[Hopt,loss]=soc_avg(Gy,Xoptnoise);
D=Hopt*Gy*inv(G1); % 
Hopt9=inv(D)*Hopt;
resb=Hopt9*X;
resb0=Hopt9*X0 ;
a91=norm(resb-Y1) ;    
a92=norm(resb0-Y10) ;
method9='    loss w/noise';   
aloss=[aloss a92];

% Summary of results
methods=[method1 method6 method7 method9]
a =[a11 a12; a61 a62; a71 a72; ; a91 a92]'

% This is for repeated runs:
na=na+1
aa=aa+a;
aaa=aa/na

end


