function [ fe, c, lm ] = OptObliProj( V, WC, f, opt, eps )
% OptObliProj regularizes the oblique projector by truncation of singular values 
%
% The projection of f onto span of V  along span of WC is regularized by truncation of
% singular values- the number of singular values is decided by minimizing:
%
%                  ||P_W f' -  P_W E_{VWC}f'|| 
%
% where P_W is the orthogonal projector onto W (the orthogonal complement of span WC) and
% E_{V,WC}) is the oblique projector onto span V along span WC.
% 
% Usage:    [ fe, c, lm ] = OptObliProj( D, WC, f, opt, eps );
%           [ fe, c, lm ] = OptObliProj( D, WC, f );
% Inputs:
%  	V       matrix, the columns of which span the space to project onto
%  	WC      matrix, the columns of which span the space to project along 
%  	f       vector, with signal to be projected
%   opt     equivalent role as in ObliProj and (for details see there) but some of them
%           may not always be good here default opt=[3,1]. 
%   eps     minimum eigenvalue to be considered nonzero 
%   []      can be used for eps or/and opt
%
% Outputs:
%   fe      regularized oblique projection of f
%   c       coefficients in the decomposition fe=V*c 
%   lm      resulting number of nonzero eigenvalues to calculate fe
%
% See  also ObliProj RegObliProj OrthProj

% See   http://www.ncrg.aston.ac.uk/Projects/HNLApprox/
%       http://www.nonlinear-approx.info/
   ND=size(V,2);
%=====================================================================
%  setting of defaults
%=====================================================================
   if (nargin<5) |(isempty(eps==1));eps=1e-10;  end
   if (nargin<4) |(isempty(opt==1));opt=[3,1];  end
%==============================================================
%  parameters to call ORTH_PROJ see that functions to change them
%============================================================== 
   tol_ort=1e-7;
   opt_ort=2; %calculate orth_proj for option 2 in DREORP
%==============================================================
delta=1; 
Wperp=NormDict(WC);
%
[pp,Q]=OrthProj(Wperp,opt_ort,tol_ort);
%
U=V-pp*V;
%
% project the signal onto the space spanned by the columns of U
%
fw=f'-pp*f';
fw2=f';
%fw2=fw;
%
opt1=opt(1);
opt2=opt(2);
if opt2==1
     switch opt1;
     case 1
      G=U'*V;
      [U1,S1,V1]=svd(G,0);s=diag(S1);
       %decide the truncation of singular values
        lr=sum(s>eps);
	ctil=U'*fw2;
	ctil2=U'*fw2;
        cc(:,1)=V1(:,1)*U1(:,1)'/s(1)*ctil2;
        cw=V1(:,1)*U1(:,1)'/s(1)*ctil;
	d(1)=norm(fw-U*cw);
	  for i=2:lr;
	  cc(:,i)= cc(:,i-1)+V1(:,i)*U1(:,i)'/s(i)*ctil2;
	  cw=cw+V1(:,i)*U1(:,i)'/s(i)*ctil;
	  d(i)=norm(fw-U*cw);
	  end
	[md,lm]=min(d);
	c=cc(:,lm);
      case 2
       G=U'*U;
       [U1,S1,V1]=svd(G,0);s=diag(S1);
       lr=sum(s>eps);
       %decide the truncation of singular values
       ctil=U'*fw2;
       ctil2=U'*fw2;
       cc(:,1)=V1(:,1)*U1(:,1)'/s(1)*ctil2;
       cw=V1(:,1)*U1(:,1)'/s(1)*ctil;
       d(1)=norm(fw-U*cw);
	 for i=2:lr;
	   cc(:,i)= cc(:,i-1)+V1(:,i)*U1(:,i)'/s(i)*ctil2;
	   cw=cw+V1(:,i)*U1(:,i)'/s(i)*ctil;
	   d(i)=norm(fw-U*cw);
	 end
	  [md,lm]=min(d);
	   c=cc(:,lm);
       case 3
       [U1,S1,V1]=svd(U,0);s=diag(S1);
       lr=sum(s>eps);
       %compute the distance  between orthogonal projections
       for k=1:lr
          d(k)=norm(fw-U1(:,1:k)*(U1(:,1:k)'*fw2));
       end
       %decide the truncation of singular values
       [md,lm]=min(d);
       W=V1(:,1)*U1(:,1)'/s(1);
       for i=2:lm;
       W=W+V1(:,i)*U1(:,i)'/s(i);
       end
       c=W*fw2;
      otherwise
      disp('wrong option to compute the pseudoinverse')
      end
else
  if opt1==1 & opt2~=1
   switch opt2
    case 2
      or=DREOr(U,tol_ort);
   case 3
     [or,rr]=qr(U,0);
   case 4
   or=orth(U);
   otherwise
   disp('wrong option to orthogonalize vectors U')
  end
   G=or'*V;
   [U1,S1,V1]=svd(G,0);s=diag(S1);
       lr=sum(s>eps);
         %decide the truncation of singular values
   ctil=or'*fw2;
   cc(:,1)=V1(:,1)*U1(:,1)'/s(1)*ctil;
   d(1)=norm(fw-U*cc(:,1));
   for i=2:lr;
    cc(:,i)= cc(:,i-1)+V1(:,i)*U1(:,i)'/s(i)*ctil;
    d(i)=norm(fw-U*cc(:,i));
   end 
    [md,lm]=min(d);
    c=cc(:,lm);
 else
    if opt1==2 & opt2~=1
     switch opt2
      case 2
      or=DREOr(U,tol_ort);
      case 3
      [or,rr]=qr(U,0);
      case 4
      or=orth(U);
      otherwise
      disp('wrong option to orthogonalize vectors U')
     end
      G=or'*U;
     [U1,S1,V1]=svd(G,0);s=diag(S1);
       lr=sum(s>eps);
     %decide the truncation of singular values
      ctil=or'*fw2;
      cc(:,1)=V1(:,1)*U1(:,1)'/s(1)*ctil;
      d(1)=norm(fw-U*cc(:,1));
      for i=2:lr;
        cc(:,i)= cc(:,i-1)+V1(:,i)*U1(:,i)'/s(i)*ctil;
        d(i)=norm(fw-U*cc(:,i));
      end
      [md,lm]=min(d); 
      c=cc(:,lm);
   end
end
end
fe=V*c;
%
%Laura REBOLLO-NEIRA 2008
%
