function w=VSps(F,x,S,P)% VSps      Vector Selection algorithm doing Partial Search
%
% w=VSps(F,x,S,P);
%-----------------------------------------------------------------------------------
% arguments:
%   F    - the normalized dictionary (or F matrix), size NxK 
%   x    - the vector to approximate, size Nx1
%   S    - number of vectors to select or non-zero weights, 1<=S<=N 
%   P    - how many branches to consider at each level, size 1xS
%          total number of combinations examined is prod(P)
%          If P=K:(-1):(K+1-S); full search is done, VSfs will usually work faster
%          since combinations that only differ in order is only done once.
%          If P=ones(1,S); this function should be equal to ORMP (VSfomp2)
%          When last vector is found all is investigated, so P(S)=K-S+1;
%   w    - the weights where at most S of the elements are non-zero, size Kx1
%-----------------------------------------------------------------------------------
% Note that F is not global in this function as in many other VSxxx functions


%----------------------------------------------------------------------
% Copyright (c) 2001-2002.  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  20.12.2001  KS: function made 
% Ver. 1.1  25.11.2002  KS: moved from ..\Frames to ..\FrameTools%----------------------------------------------------------------------

Mfile='VSps';
global gPDisplay=0;

if nargin<4
   error([Mfile,': wrong number of arguments, see help.']);
end

[N,K]=size(F);
x=x(:);
if length(x)~=N
   error([Mfile,': size of x and F do not match, see help.']);
end

S=floor(S(1));
if S<1; S=1; end;
if S>N; S=N; end;

P(S)=1;                       % even if all are tested (like ORMP)
P=min([P;K:(-1):(K+1-S)]);
P=max([P;ones(1,S)]);
gP=fliplr(P);if Display
    M=prod(P);      % the number of combinations to do
    disp([Mfile,': N=',int2str(N),', K=',int2str(K),', S=',int2str(S),', M=',int2str(M)]);end
if S>16
   disp([Mfile,': number of vectors to select is S=',int2str(S),...
         ', and that is too large.']);
   w=F(:,1:S)\x;    
   return
end

% start the search
[rb,Ib,Qb,pb]=rfun(F,x,S);
% and find the solution
wi=(Qb'*F(:,Ib))\(Qb'*x);
w=zeros(K,1);
w(Ib)=wi;if Display
    t1=' ';
    for s=1:S; t1=[t1,int2str(pb(s)),':',int2str(P(s)),' ']; end;
    disp(['best comb. was ',t1]);end

return;

% the recursive function
% F is size Nx(K-S+s) and normalized (and orthogonal to previously selected)
%   where S is total number of vectors to select 
%   and s is remaining number of vectors to select. 
% Qb is the Q-part for the QR decomposition of the selected frame vectors
% Ib is the indexes for the selected frame vectors
% rb is the norm squared for the residual, ||r||^2
function [rb,Ib,Qb,pb]=rfun(F,r,s);
global gP  % the number of combinations to try at each level is global
c=r'*F;    % the inner products, 1x(K-S+s)
if s==1    % only select one vector
   [t,Ib] = max(abs(c));
   Qb=F(:,Ib);
   rb=r-c(Ib)*Qb;
   rb=(rb'*rb);     % ||r||^2
   pb=1;
else
   [t,I]=sort(-abs(c));   % c(I) is sorted, c(I(1)) is best
   rb=0;
   for p=1:gP(s)
      J=setdiff(1:size(F,2),I(p));
      f=F(:,I(p));   % the selected frame vector
      F1=F(:,J);     % F1 is the new frame (without the selected frame vector)
      F1=F1-f*(f'*F1);                      % orthogonalize F1 to f
      t=sum(F1.*F1);t(find(t==0))=1;
      F1=F1./(ones(size(F1,1),1)*sqrt(t));  % normalize F1 
      [rc,Ic,Qc,pc]=rfun(F1,r-c(I(p))*f,s-1);  % the recursive call
      if (p==1) | (rc<rb)
         rb=rc;Ib=[I(p),J(Ic)];Qb=[f,Qc];pb=[p,pc];
      end
   end
end
return            

% test of function
%clear all
%N=8;K=16;S=5;
% want a uniform distribution within the unit ball 
%randn('state',200);   % ex 100, 200, 400
%F=randn(N,K);
%F=NormalizeF(F);
%x=randn(N,1);
%flops(0)
%wfs=VSfs(F,x,S);
%disp(['Flops for full search    : ',int2str(flops)]);
% here partial search
%P=[4,3,1,1,K-S+1];       % best for state 400
%P=[5,9,5,2,K-S+1];       % best for state 200
%P=[1,2,1,1,K-S+1];       % best for state 100
%P=[4,3,2,2,K-S+1];
%P=[4,4,2,2,1];
%flops(0)
%wps=VSps(F,x,S,P);
%disp(['Flops for partial search : ',int2str(flops)]);
% 
%flops(0);
%wfomp=VSblock(x,F,S/N,'VSfomp2',1);
%disp(['Flops for ORMP (VSfomp2) : ',int2str(flops)]);
%wfomp=full(wfomp);
%disp('The weights are (Full Search, Partial Search, VSfomp2):');
%disp([wfs,wps,wfomp]);
%disp('The errors are (Full Search, Partial Search, VSfomp2):');
%disp([norm(x-F*wfs),norm(x-F*wps),norm(x-F*wfomp)]);