% TestVS  test of Vector Selection algorithms

%----------------------------------------------------------------------
% Copyright (c) 2002-2003.  Karl Skretting.  All rights reserved.
% Hogskolen in Stavanger (Stavanger University), Signal Processing Group
% Mail:  karl.skretting@tn.his.no   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:
% Ver. 1.0  19.09.2003  KS: m-file made 
%----------------------------------------------------------------------

clear all;
global F FF

v=version;
if (v(1)=='6'); NoFlops=1; else NoFlops=0; end;

N=20;
K=50;
S=8;
% K=34;S=4;   % makes full search possible

% randn('state',208);
randn('state',sum(100*clock));
F=randn(N,K);
F=NormalizeF(F);
FF=F'*F;
x=randn(N,1);
disp(['TestVS using size(F)=',int2str(N),'x',int2str(K),...
      '  S=',int2str(S),', and norm(x)=',num2str(norm(x)),'.']);

if (exist('VSmp2')==2)  
    if (NoFlops==0); flops(0); end;
    w=VSmp2(x,S,S);    % as few iterations as possible
    if (NoFlops==0); FC=flops; else; FC=0; end;  
    xr=F*w;
    xe=x-xr;
    disp(['VSmp2(,S)   gives norm of error =',num2str(norm(xe)),...
          ', norm of weighs=',num2str(norm(w)),'  (flops=',int2str(FC),')']);
    if (NoFlops==0); flops(0); end;
    w=VSmp2(x,S);
    if (NoFlops==0); FC=flops; else; FC=0; end;  
    xr=F*w;
    xe=x-xr;
    disp(['VSmp2       gives norm of error =',num2str(norm(xe)),...
          ', norm of weighs=',num2str(norm(w)),'  (flops=',int2str(FC),')']);
    if (NoFlops==0); flops(0); end;
    w=VSmp2(x,S,S*6);
    if (NoFlops==0); FC=flops; else; FC=0; end;  
    xr=F*w;
    xe=x-xr;
    disp(['VSmp2(,S*6) gives norm of error =',num2str(norm(xe)),...
          ', norm of weighs=',num2str(norm(w)),'  (flops=',int2str(FC),')']);
end

if (exist('VSmp2')==2) & (exist('VSab2')==2) 
    if (NoFlops==0); flops(0); end;
    for i=1:10 
        w=VSab2(x,w); 
    end
    if (NoFlops==0); FC=flops; else; FC=0; end;  
    xr=F*w;
    xe=x-xr;
    disp(['VSab2       gives norm of error =',num2str(norm(xe)),...
          ', norm of weighs=',num2str(norm(w)),'  (flops=',int2str(FC),')']);
end

if (exist('VSomp')==2)  
    if (NoFlops==0); flops(0); end;
    w=VSomp(x,S);
    if (NoFlops==0); FC=flops; else; FC=0; end;  
    xr=F*w;
    xe=x-xr;
    disp(['VSomp       gives norm of error =',num2str(norm(xe)),...
          ', norm of weighs=',num2str(norm(w)),'  (flops=',int2str(FC),')']);
end

if (exist('VSompJHH')==2)  % another implementation of OMP
    if (NoFlops==0); flops(0); end;
    w=VSompJHH(x,S);
    if (NoFlops==0); FC=flops; else; FC=0; end;  
    xr=F*w;
    xe=x-xr;
    disp(['VSompJHH    gives norm of error =',num2str(norm(xe)),...
          ', norm of weighs=',num2str(norm(w)),'  (flops=',int2str(FC),')']);
end

if (exist('VSormp')==2)  
    if (NoFlops==0); flops(0); end;
    w=VSormp(x,S);
    if (NoFlops==0); FC=flops; else; FC=0; end;  
    xr=F*w;
    xe=x-xr;
    disp(['VSormp      gives norm of error =',num2str(norm(xe)),...
          ', norm of weighs=',num2str(norm(w)),'  (flops=',int2str(FC),')']);
end

if (exist('VSormp2')==2)   % another implementation of ORMP
    if (NoFlops==0); flops(0); end;
    w=VSormp2(x,S);
    if (NoFlops==0); FC=flops; else; FC=0; end;  
    xr=F*w;
    xe=x-xr;
    disp(['VSormp2     gives norm of error =',num2str(norm(xe)),...
          ', norm of weighs=',num2str(norm(w)),'  (flops=',int2str(FC),')']);
end

if (exist('VSfomp2')==2)   % another implementation of ORMP
    if (NoFlops==0); flops(0); end;
    w=VSfomp2(x,S);
    if (NoFlops==0); FC=flops; else; FC=0; end;  
    xr=F*w;
    xe=x-xr;
    disp(['VSfomp2     gives norm of error =',num2str(norm(xe)),...
          ', norm of weighs=',num2str(norm(w)),'  (flops=',int2str(FC),')']);
end

if (exist('VSps')==2)   % partial search, recursive (and slow) algorithm
    if (NoFlops==0); flops(0); end;
    P=ones(1,S);P(1:3)=[5,4,2];
    w=VSps(F,x,S,P);
    if (NoFlops==0); FC=flops; else; FC=0; end;  
    xr=F*w;
    xe=x-xr;
    disp(['VSps 5-4-2  gives norm of error =',num2str(norm(xe)),...
          ', norm of weighs=',num2str(norm(w)),'  (flops=',int2str(FC),')']);
end

if (exist('VSps2')==2)   % partial search, (effective algorithm)
    M=4; 
    if (NoFlops==0); flops(0); end;
    w=VSps2(x,S,M);
    if (NoFlops==0); FC=flops; else; FC=0; end;  
    xr=F*w;
    xe=x-xr;
    disp(['VSps2 M=4   gives norm of error =',num2str(norm(xe)),...
          ', norm of weighs=',num2str(norm(w)),'  (flops=',int2str(FC),')']);
    M=40; 
    if (NoFlops==0); flops(0); end;
    w=VSps2(x,S,M);
    if (NoFlops==0); FC=flops; else; FC=0; end;  
    xr=F*w;
    xe=x-xr;
    disp(['VSps2 M=40  gives norm of error =',num2str(norm(xe)),...
          ', norm of weighs=',num2str(norm(w)),'  (flops=',int2str(FC),')']);
    M=400; 
    if (NoFlops==0); flops(0); end;
    w=VSps2(x,S,M);
    if (NoFlops==0); FC=flops; else; FC=0; end;  
    xr=F*w;
    xe=x-xr;
    disp(['VSps2 M=400 gives norm of error =',num2str(norm(xe)),...
          ', norm of weighs=',num2str(norm(w)),'  (flops=',int2str(FC),')']);
end;

if (exist('VSfs')==2) & (nchoosek(K,S)<50000)  % full search
    if (NoFlops==0); flops(0); end;
    w=VSfs(F,x,S);
    if (NoFlops==0); FC=flops; else; FC=0; end;  
    xr=F*w;
    xe=x-xr;
    disp(['VSfs        gives norm of error =',num2str(norm(xe)),...
          ', norm of weighs=',num2str(norm(w)),'  (flops=',int2str(FC),')']);
end

disp(' ');
return