function [ D, Di, beta, Q ] = OMPSwapping( f, D, Di, beta, Q ) 
% OMPSwapping Swapping based refinement of OMP method 
%
% It interchanges at each step one atom from the atomic decomposition with another atom
% from the dictionary to improve the OMP approximation via adaptive biorthogonalization.
% At each step it modifies the biorthogonal vectors giving rise to the duals of selected
% atoms. The inputs are obtainable from the outputs of  the OMP function
%
% Usage:    [ D, Di, Q, beta ] = OMPSwapping( f, D, Di, Q, beta );
%
% Inputs:
%   D       dictionary, first k functions D(:,1:k) are the selected basis
%   Di      indices of atoms in D with respect to the original dictionary 
%   beta    biorthogonal vectors to D(:,1:k), k=size(beta,2)
%   Q       orthonormal basis spanning the same space as D(:,1:k)
%     
% Outputs:
%   D       updated (re-arranged) dictionary, D(:,1:k) are the selected basis
%   beta    biorthogonal functions to D(:,1:k), k=size(beta,2)  
%   Di      indices of atoms in D with respect to the original dictionary 
%   Q       Q orthonormal basis spanning the same space as D(:,1:k)
% 
% References: 
%   M. Andrle and L. Rebollo-Neira, "A swapping-based refinement of orthogonal matching
%   pursuit strategies", Signal Processing, Vol 86, No 3, pp. 480-495, 2006.
%
% See also OOMPSwapping, OOMPKSwaps, OMPKSwaps, OMPKSwapping

% See   http://www.ncrg.aston.ac.uk/Projects/HNLApprox/
%       http://www.nonlinear-approx.info/
  
tol1=1e-07;%defaul tol1=1e-07;
tic;
%initilization
N=size(D,2);
[L,k]=size(beta);
%delta=1/L;
delta=1;
C=f*beta;
Re=f'-D(:,1:k)*C';
nor1=norm(Re)*sqrt(delta);nor1=nor1*nor1;
fprintf('Initial norm = %g\n\n',nor1);
nor2=0;z=0;

while nor1>nor2
  if nor2~=0 nor1=nor2;else nor2=nor1;end
  dd=zeros(1,k);%clear dd
  for i=1:k  
    dd(i)=beta(:,i)'*beta(:,i);
  end	
  [Cp,p]=min((C.*conj(C))./dd); % p-th atom to be removed
  pp=Di(p);
  [D,Q,beta]=BioBack(D,Q,beta,p);
  dummy=Di(p);Di(p)=[];Di(N)=dummy;
  nor2=nor2+(Cp);
  Re=f-f*Q*Q';
  c=zeros(1,N);cc=zeros(1,N);
  for i=k:N
    cc(i)=abs(Re*D(:,i));
  end
  [max_c,q]=max(cc);qq=Di(q);
  %nor2=nor2-(max_c)^2;
  fprintf('%d  %d~%d\n',pp,q,qq)
    if q~=k
      D(:,[k q])=D(:,[q k]);% swap k-th and q-th columns
      Di([k q])=Di([q k]);
    end	  
  [Q,beta]=BioFor(D,Q,beta,k);
  nor2=nor2- (abs(f*(beta(:,k)/norm(beta(:,k)))))^2;
  if pp==qq fprintf('Stop\n');z=z;break; end
  C=f*beta;
  fprintf('Norm %g\n',sqrt(nor2));
  z=z+1;
  if nor1<=nor2 fprintf('\n nor1<=nor2\n');end
end
fprintf('\n%d refinenments of %d functions lasted %g seconds\n',z,k,toc);
fprintf('New norm  residuum norm (calculated at the end of algorithm) is %g\n',norm(f'-D(:,1:k)*(f*beta)')*sqrt(delta))
fprintf('Residuum norm calculated additively after each iteration: %g\n',sqrt(nor2));
fprintf('\nIf the previous two norms differ much, it means that the algorithm has not been too stable.\n');
ErrorTest(D(:,1:k),beta);
ErrorTest(Q(:,1:k));
%
%Laura Rebollo-Neira (2009)- Based on SWAPPING.m copyright below. 
%
%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.
