
% This is file matlab-test-cases.m
% Load data from C:\Documents and Settings\skoge\My Documents\MATLAB 

% -------------------------------------------------
% Test case 1. Martens data: 
load GlutenStarchNIR
%    X: Variablene  13-112 er NIR absorbanser (log(1/T) fra 850-1048 nm). 
%    y1: Variabel 4 er kvantitiv analytt-konsentrasjon (gluten),
% Matrix: 100 rows  112 cols: regular MATLAB matrix
% VarLabels: 112 rows  17 cols: MATLAB character string
% ObjLabels: 100 rows  7 cols: MATLAB character string
Yall=Matrix';
XX = Yall(13:112,:);
YY = Yall(4,:);
X = XX(:,1:50);
X0= XX(:,51:100);
Y1= YY(:,1:50);
Y10=YY(:,51:100);

% 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)

% Jump to.. MATLAB code starts here

% -------------------------------------------------
% Test case 2.
% Bjørn Alsberg data (epost 15.11 2010)
load wheat_spectra
%Your variables are:
%X        Y        idxcal   idxval   moist    protein  
Xwheat=X'; Ywheat=Y';  % 2 y's, 701 x's, 100 data set
X =Xwheat(:,idxcal); Y1 =Ywheat(:,idxcal); % 50 calibration sets
X0=Xwheat(:,idxval); Y10=Ywheat(:,idxval); % 50 validation sets 
%X =Xwheat(:,idxval); Y1 =Ywheat(:,idxval); % 50 calibration sets - switched
%X0=Xwheat(:,idxcal); Y10=Ywheat(:,idxcal); % 50 validation sets - switched


% Jump to.. MATLAB code starts here

% -------------------------------------------------
% Test case 3. Own example

G1 = [1 0; 0 1] 
Gd1= [0 0; 0 0]
gy1 =[0.2  0  3  4  5   6 3] 
gy2 =[0  0.2  1  3  4   8 9]
Gy = [gy1',gy2']
gd1 =[0  0  4  5   6  8  -9]
gd2 =[0  0  -3  -5   5  9 18]
Gyd = [gd1',gd2']

Y10 = [G1  Gd1]
X0 = [Gy Gyd]

Y100=Y10;
na=0; aa=a*0;

% Run repeatedly from here for several cases

Y10=Y100;

% 8 data sets
Noise = 0.5*randn(7,8)
X =  [X0 -X0] + Noise
Y1 = [Y10 -Y10]  % use Y10 and -Y10 to get centered data

%%32 data sets
%X =  [X0 -X0] 
%Noise = 0.5*randn(7,32)
%X = [X X X X] + Noise
%Y1 = [Y10 -Y10]  % use Y10 and -Y10 to get centered data
%Y1 = [Y1 Y1 Y1 Y1] 
%%NoiseY= 0.0*randn(2,32)  % with 2 y1's
%%Y1 = [Y1 Y1 Y1 Y1] + NoiseY
% 100 times where I take out data each time

%--------------------------------------------------
% 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); 

% PCR (vary  tol)
% PCR: To find cut-off to find no. of components, semilogy(svd(X))
method6='   PCR';
tol=1.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); 

% Weighted least squares (to get relative noise assumption rather than additive noise)
method5='   weightedLS';
[rx,cx]=size(X); 
mag=[];
for i = 1:rx
mag=[mag norm(X(i,:))]; % the magnitudes are about 0.8
end
Xs = inv(diag(mag))*X;
Hlss1 = Y1*pinv(Xs);
Hlss = Hlss1*inv(diag(mag)); 
ress=Hlss*X;
ress0=Hlss*X0;
a51=norm(ress-Y1) ;       
a52=norm(ress0-Y10) ;
 
% pls
npls=35
%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';

% 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); 

% step 2, 
%Hopt = (pinv(Xopt*Xopt')*Gy)'; %analytical expression 
[Hopt,loss]=soc_avg(Gy,Xopt);
D=Hopt*Gy*inv(G1); % 
Hopt4=inv(D)*Hopt;
resb=Hopt4*X;
resb0=Hopt4*X0 ;
a41=norm(resb-Y1) ;    
a42=norm(resb0-Y10) ;
method4='    loss '; 

% NEW modified loss method for case with noise on Y: redo step 1.
% New estimate of Y1.
Y1smooth = Hopt4*X; 
% then redo step 1 and 2
Yallsmooth = [Y1smooth; X];
[u,s,v]=svd(Y1smooth);
Yallsmooth = Yall*v;
G1smooth=Yallsmooth(1:ny,1:ny);
Gysmooth=Yallsmooth(ny+1:r1,1:ny);
Xsmooth = Yallsmooth(ny+1:r1,ny+1:c1); 
% step 2
[Hopt,loss]=soc_avg(Gysmooth,Xsmooth);
D=Hopt*Gysmooth*inv(G1smooth); % 
Hopt10=inv(D)*Hopt;
resa=Hopt10*X;
resa0=Hopt10*X0 ;
a101=norm(resa-Y1) ;    
a102=norm(resa0-Y10) ;   
method10='   loss w/smooth';

% loss method: add artificial noise weights 
noisemag=2.e-4
%noisemag=tol
% Noise Alt.1 just additive noise:
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';   

% Summary of results
methods=[method1 method5 method6 method7 method4 method10 method9]
a =[a11 a12; a51 a52; a61 a62; a71 a72; ; a41 a42; a101 a102; a91 a92]'

% For repeated case 1 and case 3:
na=na+1
aa=aa+a;
aaa=aa/na
% For repeated case 1
end


