function r=MCsimF(X,Fin,S,VSalg);
% MCsimF    Monte Carlo simulation for the frame F, selecting S vectors
%           returning the norm of the error for each 
% training vector (or each simulation). The data to test is given in the
% X matrix, or X may indicate which kind of random data to use.
% This function is quite similar to VSblock, but does not return the weights.
%
% r=MCsimF(X,F,S);
% r=MCsimF(X,F,S,VSalg);
%-----------------------------------------------------------------------------------
% arguments:
%   X     - the data, a matrix of siz NxL
%           or if X is a 2x1 (or 1x2) vector then X(2) is L and
%           X(1)==1  for random distribution within unit ball (radius==1) 
%           X(1)==2  for random distribution within cube, uniform in range [-1,1] 
%           X(1)==3  random distribution, each variable is Gaussian. 
%                    Note, this is uniform on the unit ball!
%   F     - the normalized dictionary (or F matrix), size NxK 
%   S     - number of vectors to select or non-zero weights, 1<=S<=N 
%   VSalg - Which vector selection algorithm to use, alternatives are:
%           VSfs (default), VSfomp2
%   r     - norm of error for each try, size 1xL
%-----------------------------------------------------------------------------------

%----------------------------------------------------------------------
% Copyright (c) 2001.  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  06.12.2001  KS: function made 
% Ver. 1.1  19.12.2001  KS: changed arguments 
% Ver. 1.2  04.12.2002  KS: moved from ..\Frames\ to ..\FrameTools
%----------------------------------------------------------------------

Mfile='MCsimF';

if nargin<4
   VSalg='VSfs';      % full search
end
if nargin<3
   error([Mfile,': wrong number of arguments, see help.']);
end
if strcmp(VSalg,'VSfomp2')
   global F FF
   FF=F'*F;
end
F=Fin;
clear Fin;

[N,K]=size(F);
if prod(size(X))==2
   Xtype=X(1);
   L=X(2);
else
   [n,L]=size(X);
   if n~=N
      error([Mfile,': wrong size of argument 1 (X), see help.']);
   end
   Xtype=0;       % as argument X
end   

if Xtype==1        % within unit ball
   X=zeros(N,L);   % none made so far
   l=0;
   while l<L
      x=2*rand(N,1)-1;    % uniform within unit cube (-1, 1)
      xx=x'*x;
      if (xx<1) % within unit ball 
         l=l+1;
         X(:,l)=x/sqrt(xx);     % on unit ball
      end  
   end
elseif Xtype==2    % Uniform
   X=2*rand(N,L)-1;
   X=X./(ones(N,1)*sqrt(sum(X.*X)));     % normalize each column vector
elseif Xtype==3    % Gaussian
   X=randn(N,L);
   X=X./(ones(N,1)*sqrt(sum(X.*X)));     % normalize each column vector
else
   X=X./(ones(N,1)*sqrt(sum(X.*X)));     % normalize each column vector
end

r=zeros(1,L);
normw=zeros(1,L);

%
if strcmp(VSalg,'VSfs');
   if S==1
      c=F'*X;       % cosine of angles
      normw=max(abs(c));
      r=sqrt(1-normw.*normw);
      [temp,i]=max(r);
      %disp(X(:,i));
   elseif S==2
      for i1=1:(K-1)
         ym=zeros(K-i1,L);
         for i2=(i1+1):K
            [Qm,Rm]=qr(F(:,[i1,i2]),0);  
            Qm=Qm(:,find(diag(Rm)));
            Ym=Qm'*X;  
            ym(i2-i1,:)=sum(Ym.*Ym); 
         end
         r=max([r;ym]);    % actually best y'*y = 1-r'*r
      end
      r=sqrt(1-r);
   elseif S==3
      for i1=1:(K-2)
         for i2=(i1+1):(K-1)
            ym=zeros(K-i2,L);
            for i3=(i2+1):K
               [Qm,Rm]=qr(F(:,[i1,i2,i3]),0);  
               Qm=Qm(:,find(diag(Rm)));
               Ym=Qm'*X;
               ym(i3-i2,:)=sum(Ym.*Ym); 
            end
            r=max([r;ym]);    % actually best y'*y = 1-r'*r
         end
      end
      r=sqrt(1-r);
   elseif S==4
      for i1=1:(K-3)
         for i2=(i1+1):(K-2)
            for i3=(i2+1):(K-1)
               ym=zeros(K-i3,L);
               for i4=(i3+1):K
                  [Qm,Rm]=qr(F(:,[i1,i2,i3,i4]),0);  
                  Qm=Qm(:,find(diag(Rm)));
                  Ym=Qm'*X;
                  ym(i4-i3,:)=sum(Ym.*Ym); 
               end
               r=max([r;ym]);    % actually best y'*y = 1-r'*r
            end
         end
      end
      r=sqrt(1-r);
   elseif S==5
      for i1=1:(K-4)
         for i2=(i1+1):(K-3)
            for i3=(i2+1):(K-2)
               for i4=(i3+1):(K-1)
                  ym=zeros(K-i4,L);
                  for i5=(i4+1):K
                     [Qm,Rm]=qr(F(:,[i1,i2,i3,i4,i5]),0);  
                     Qm=Qm(:,find(diag(Rm)));
                     Ym=Qm'*X;
                     ym(i5-i4,:)=sum(Ym.*Ym); 
                  end
                  r=max([r;ym]);    % actually best y'*y = 1-r'*r
               end
            end
         end
      end
      r=sqrt(1-r);
   else   % S>=6
      for l=1:L
         x=X(:,l);
         w=VSfs(F,x,S);     % do full search
         rm=x-F*w;
         r(l)=sqrt(rm'*rm);
      end
   end
elseif strcmp(VSalg,'VSfomp2');
   for l=1:L
      x=X(:,l);
      w=VSfomp2(x,S);     
      rm=x-F*w;
      r(l)=sqrt(rm'*rm);
   end
end

return


% test of Gaussian distribution, this shows that it is uniform on the unit ball
j=sqrt(-1);
N=2;L=20000;
X=randn(N,L);
X=X./(ones(N,1)*sqrt(sum(X.*X)));     % normalize each column vector
theta=zeros(1,L);
theta=angle(X(1,:)+j*X(2,:));
theta=sort(theta);
% pdffun=PDFpoly(theta,100,[-pi,pi],0); 
pdffun=PDFpoly(theta,200,[-pi,-pi/2,0,pi/2,pi],[3,3,3,3],[2,2,2]);