function  VFSwapping( f, D )
% VFSwapping Swapping based refinement of OMP methods (inner product implementation)
%
% It interchange at each step one atom from the atomic decomposition for another atom from
% the dictionary to improve the signal approximation.  The process is carried out until
% the approximation error will not increase. The routine is not modifying biorthogonal
% functions or unselected atoms subtracted by the component in the selected space but only
% several related inner product matrices. See Alg. 4b in the reference below. 
%  
% Usage:    VFSwapping( f, D );
%
% Inputs:
%  f        signal  
%  D        dictionary  
%  
% These inputs and all outputs are realized through global variables AL LL LF J IJ I
%   I       set of all indices
%   J       set of indices (subset of I) which are not in the basis
%   IJ      set of indices (subset of I) which are in the basis
%   AL      inner products between dictionary D and lambdas
%   LL      inner products between lambda(IJ), between biorthogonal functions to D(:,IJ)
%   LF      inner products between lambda and signal f 
%
% Note:     lambda(IJ) means biorthogonal functions to D(IJ)
%           lambda(J) means D(J) subtracted by their component in selected basis D(IJ) 
%
% References: 
%   M. Andrle and L. Rebollo-Neira, "A swapping-based refinement of orthogonal matching 
%   pursuit strategies", Signal Processing, Vol (86,3), 480-495 (2006).
%
% See also Swapping, KSwapping.  

% See   http://www.ncrg.aston.ac.uk/Projects/HNLApprox/
%       http://www.nonlinear-approx.info/
  
global  AL LL  LF J IJ I;
tol1=1e-7;
tol2=1e-10;

tic;
%initilization
N=length(I); %number of all atoms
k=length(IJ); %number of selected atoms
delta=1;

nor1=norm(f'-D(:,IJ)*LF(IJ)')*sqrt(delta);
nor_old=nor1;
nor1=nor1*nor1;
fprintf('Initial norm = %g\n\n',nor_old);
nor2=0;z=0;

while nor1>nor2
  if nor2~=0 nor1=nor2;else nor2=nor1;end
  %clear q;
  [min_q,j]=min((LF(IJ).*conj(LF(IJ)))./diag(LL).'); % lj-th atom to be removed
  lj=IJ(j);
  vfastbackward(lj);
  nor2=nor2+(min_q);
  fprintf('Norm %g\n',sqrt(nor2))
  e=zeros(1,N);  
  for i=J
    c=sqrt(AL(i,i));
    if c>tol1 e(i)=abs(LF(i))/c;end
  end
  [max_e,lk]=max(e);
  %nor2=nor2-(max_e)^2;
  fprintf('%d  %d\n',lj,lk)
  %re-orthogonalization
  for i=IJ
    AL(:,lk)=AL(:,lk)-AL(i,lk)*AL(:,i);
    LF(lk)=LF(:,lk)-conj(AL(i,lk))*LF(i);
  end
  %end of re-orthogonalization
  nor2=nor2-(LF(lk)/sqrt(AL(lk,lk)))^2;
  vfastforward(lk);
  fprintf('%g\n',norm(f'-D(:,IJ)*LF(IJ)'));
  fprintf('Norm %g\n',sqrt(nor2));
  if lj==lk fprintf('Stop\n');break; end
  z=z+1;
  if nor1<=nor2 fprintf('\n nor1<=nor2\n');end
end
k=size(LL,1);
IJ=setdiff(I,J);
fprintf('\n%d refinenments of %d functions lasted %g seconds\n',z,length(IJ),toc);
fprintf('The original residuum norm was %g\n',nor_old);
fprintf('New residuum norm (calculated at the end of algorithm) is %1.10f\n',norm(f'-D(:,IJ)*LF(IJ)')*sqrt(delta));
fprintf('Residuum norm calculated additively after each iteration: %1.10f\n',sqrt(nor2)*sqrt(delta));
fprintf('\nIf the previous two norms differ much, it means that the algorithm has not been too stable.\n');
fprintf('Biorthogonality: %g\n',norm(eye(k)-AL([IJ],[IJ])));


%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.
