function [ D, Di, beta, c, Q ] = BOOMPQ( f, D, Di, beta, Q, tol, No )
% BOOMPQ Backward-Optimized Orthogonal Matching and gives the orthonormal 
% basis
%
% Using the Least Square criterion at each step it eliminates one function from a given
% basis to have best possible representation in the reduced space. It also modifies the
% corresponding biorthogonal functions and orthonormal basis (obtained by modified
% Gram-Schmidt).
%  
% Usage:    [ D, Di, Q, beta, c ] = BOOMPQ( f, D, Di, Q, beta, toli, No);
%           [ D, Di, Q, beta, c ] = BOOMPQ( f, D, Di, Q, beta, tol );
%
% Inputs:  
%   f       analyzing signal D dictionary D(:,1:N) is the selected basis, N=size(beta,2)
%   Di      indices of atoms Q Q(:,1:N) orthonormal set spanning D(:,1:N), Q(:,N+1:end) 
%           the rest of the dictionary orthogonalized w.r.t Q(:,1:N) 
%   beta    set of biorthogonal functions to D(:,1:N) tol tolerance, difference between
%           signal and approx, (optional, default tol=1.0e-2) 
%   No      (optional) desired number of atoms in the decomposition if you want really No
%           atoms set tol=0, it speeds the process 
%  
% Outputs: 
%   D       rearranged dictionary
%           D(:,1:s) reduced basis, s=size(beta,2)  
%   Di      indices of atoms in new D w.r.t. to original D Q Q(:,1:s) orthonormal set 
%           spanning D(:,1:N), Q(:,s+1:end) the rest of the dictionary orthogonalized w.r.t
%           Q(:,1:s)  
%   beta    biorthogonal functions to new D(:,1:s) c coefficients of the atomic
%           decomposition  
%  
% References:
%   M. Andrle, L. Rebollo-Neira, and E. Sagianos, "Backward-Optimized Orthogonal Matching
%   Pursuit Approach", IEEE Signal Processing Letters, Vol (11,9), 705-708 (2004). 

% See   http://www.ncrg.aston.ac.uk/Projects/HNLApprox/
%       http://www.nonlinear-approx.info/

name='BOOMPQ';  
tic;
[L,N]=size(beta); 
Ni=size(Di,2);
if nargin==5; tol=1.0e-2;end; 
if nargin==6; No=0;end;
if No>N; error('%s: You want to reduce %d atoms into %d atoms!',name,N,No);end
%initilization
delta=1;% 1/L;
C=f*beta;
c=C;   
nor1=norm(f'-D(:,1:N)*C')*sqrt(delta);
nor2=nor1;

if (nor1>tol) & (No==0)
  fprintf('\n%s: Approximation is too coarse, I cannot satisfy the tolerance condition\n',name);
else
  for n=1:N-No
    s=size(beta,2);
    dd=zeros(1,s);%clear dd
    for i=1:s  
      dd(i)=beta(:,i)'*beta(:,i);
    end	
    [Cp,p]=min((C.*conj(C))./dd); % p-th atom to be removed
    temp=(beta'*beta(:,p))*C(p)/dd(p);
    C=C-temp';
    %stopping criterion
    if (tol~=0) & (norm(f'-D(:,1:s)*C')*sqrt(delta)>tol); break; end
    [D,Q,beta]=BioDictDel(D,Q,beta,p);
    dummy=Di(p);Di(p)=[];Di(Ni)=dummy;
    C(p)=[];c=C;
  end
  s=size(beta,2);
  nor2=norm(f'-D(:,1:s)*c')*sqrt(delta);
end  

fprintf('\n%s reduces the space of %d functions into a space of %d functions\n',name,N,size(beta,2));
fprintf('Previous norm ||f-f_approx|| was %g, new norm is %g\n',nor1,nor2);
if No~=0 
  fprintf('You have required %d atoms in decomposition\n',No);
else  
  fprintf('Required precision was tol=%g\n',tol);
end  
fprintf('%s ran for %f seconds\n',name,toc);

%Last modification Laura REBOLLO-NEIRA (2009)
%
%Copyright (C) 2006 Miroslav ANDRLE and Laura REBOLLO-NEIRA
%
%This program is free software; you can redistribute it and/or modify it under the terms 
%of the GNU General Public License as published by the Free Software Foundation; either 
%version 2 of the License, or (at your option) any later version.
%
%This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
%without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
%See the GNU General Public License for more details.
%
%You should have received a copy of the GNU General Public License along with this program;
%if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
%Boston, MA  02110-1301, USA.
