function [ fe, c, lm ] = RegObliProj( V, WC, f, err, opt, No )
% RegObliProj regularizes the oblique projection of a given noisy signal by  truncation
% of singular valueas  
%
% The regularization tries to fulfill:
%
% ||P_W f' -  P_W E_{VWC}f'||\le ||  P_W err'||, 
%
% 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 ] = RegObliProj( D, WC, f, err, opt, No );
%           [ fe, c, lm ] = RegObliProj( 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 (for details see there) default opt=[1,2] (with
%           tol_orth= 1e-7 for DREOr) 
%   err     (optional) error of each point of f, or tolerance for the error's norm default
%           err=ones(size(f))*5e-7 
%   No      maximum number of eigenvalues to consider default all
%   []      can be used for err, opt, or No
%
% 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   http://www.ncrg.aston.ac.uk/Projects/HNLApprox/
%       http://www.nonlinear-approx.info/
   ND=size(V,2);
   delta=1;   
%=====================================================================
% setting of default tolerance parameter
%=====================================================================
   if ((nargin<6) | (isempty(No)==1)); No=ND;    end
   if ((nargin<5) | (isempty(opt)==1)); opt=[1,2];  end
   if ((nargin<4) | (isempty(err)==1)); err=ones(1,length(f))*5e-7; end
%
   if max(size(err))==1 tol=err;
   else tol=norm(err);
   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
%==============================================================
Wperp=NormDict(WC);
%
[Q,pp]=OrthProj(Wperp,opt_ort,tol_ort);
%
U=V-pp*V;
%
% project the signal and error onto the space spanned by the columns of U
%
fw=f'-pp*f';
fw2=f';
%
if max(size(err))==1 tol=err;
    else tol=norm(err' -pp*err');
end
%
opt1=opt(1);
opt2=opt(2);
%
if opt2==1
     switch opt1;
     case 1
      G=U'*V;
      ctil=U'*fw2;
      [U1,S1,V1]=svd(G,0);s=diag(S1);
      sd=length(s);
      c=V1(:,1)*U1(:,1)'/s(1)*ctil;
      sm=min(sd,No);
      for i=2:sm;
       gi=V1(:,i)*U1(:,i)'/s(i);
       c=c+gi*ctil;
       fwp=U*c;
       er=norm(fwp-fw);
        if (tol~= 0) & (er*sqrt(delta) < tol) break;end;
      end
      lm=i;
      case 2
       G=U'*U;
       ctil=U'*fw2;
       [U1,S1,V1]=svd(G,0);s=diag(S1);
       sd=length(s);
       c=V1(:,1)*U1(:,1)'/s(1)*ctil;
       sm=min(sd,No);
       for i=2:sm;
        gi=V1(:,i)*U1(:,i)'/s(i);
        c=c+gi*ctil;
        fwp=U*c;
        er=norm(fwp-fw);
         if (tol~= 0) & (er*sqrt(delta) < tol) break;end;
       end
       lm=i;
     case 3
      [U1,S1,V1]=svd(U,0);s=diag(S1);
      sd=length(s);
      c=V1(:,1)*U1(:,1)'/s(1)*fw2;
      sm=min(sd,No);
      for i=2:sm;
       gi=V1(:,i)*U1(:,i)'/s(i);
       c=c+gi*fw2;
       fwp=U*c;
       er=norm(fwp-fw);
        if (tol~= 0) & (er*sqrt(delta) < tol) break;end;
      end
      lm=i;
      otherwise
      disp('wrong option 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;
   ctil=or'*fw2;
   [U1,S1,V1]=svd(G,0);s=diag(S1);
   sd=length(s);
   c=V1(:,1)*U1(:,1)'/s(1)*ctil;
   sm=min(sd,No);
   for i=2:sm;
    gi=V1(:,i)*U1(:,i)'/s(i);
    c=c+gi*ctil;
    fwp=U*c;
    er=norm(fwp-fw);
    if (tol~= 0) & (er*sqrt(delta) < tol) break;end;
   end
   lm=i;
 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;
      ctil=or'*fw2;
      [U1,S1,V1]=svd(G,0);s=diag(S1);
      sd=length(s);
      c=V1(:,1)*U1(:,1)'/s(1)*ctil;
      sm=min(sd,No);
      for i=2:sm;
       gi=V1(:,i)*U1(:,i)'/s(i);
       c=c+gi*ctil;
       fwp=U*c;
       er=norm(fwp-fw);
        if (tol~= 0) & (er*sqrt(delta) < tol) break;end;
      end
      lm=i;
    end
  end
end
fe=V*c;
%
%Laura REBOLLO-NEIRA 2008
%
