function [ re, resid, D, Di, beta, C, Q ] = KSwapping( f, D, Di, beta, C, Q, swi, sws, tol )
% KSwapping extends Swapping to considering K swaps
% 
% Given an initial approximation of f, it improves upon the approximation by interchanging
% swi-pairs of atoms (swi from the approximation and  swi from the dictionary) then
% (swi+1)-atoms and (swi+1)-atoms, (swi+2)-atoms and (swi+2)-atoms and so on up to
% sws-atoms, unless the desired precision tol has been reached. (See Ref[1]). If the
% number of atoms involved in the swapping process is equal to sws and the stopping
% criterion based on the precision tol has not been reached the function returns the value
% re=0. Note: The inputs are obtainable by running OOMP first (see the example)
%
% Usage:    [ re, resid, D, Di, beta, C, Q ] = KSwapping( f, D, Di, Q, beta, C, swi,...
%               sws, tol );
%           [ re, resid, D, Di, beta, C, Q ] = KSwapping( f, D, Di, Q, beta, C );
%
% Inputs:
%   f       signal to be decomposed   
%   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 functions to D(:,1:k), k=size(beta,2)
%   C       coefficients in the expansion  
%   Q       Q(:,1:k) orthonormal basis spanning the same space as D(:,1:k), Q(:,k+1:end)
%           unselected atoms subtracted by their component in D(:,1:k) 
%   swi     minimum number of atoms to be swapped, default swi=1 
%   sws     maximum number of atoms to be swapped
%   tol     tolerance for the approximation, default tol= 1e-8
%
% Outputs:
%   re      convergence indicator: re=1 if  the method converges within the given tol and
%           re=0 otherwise 
%   resid   vector of length sws to store the residuals at each swapping. The first
%           component of resid is the residual when swi atoms are swapped, the second
%           component is the residual when (swi+1) atoms are swapped and so on. If the
%           swapping is started from swi atoms, resid is of length sws-swi+1
%   D       updated (re-arranged) dictionary, D(:,1:k) is the selected basis
%   Di      indices of atoms in D with respect to the original dictionary  
%   beta    biorthogonal vectors to D(:,1:k) 
%   C       coefficients in the expansion
%   Q       Q(:,1:k) orthonormal basis spanning the same space as D(:,1:k), Q(:,k+1:end)
%           unselected atoms subtracted by their component in D(:,1:k)  
%   
% References:
%  [1] M. Andrle and L. Rebollo-Neira, "Improvement of Orthogonal Matching Pursuit
%  strategies by Backward and Forward movements," in Proc. of the 31st International
%  Conference on Acoustics, Speech, and Signal Processing (ICASSP'06) 
%  [2] M. Andrle and L. Rebollo-Neira, "A swapping-based refinement of orthogonal matching
%  pursuit strategies", Signal Processing, Vol (86,3), 480-495 (2006) 
%  [3] L. Rebollo-Neira, "Measurements design and phenomena discrimination", J. Phys. A:
%  Math. Theor. 42 (2009)
%
% See also OOMPKSwaps OMPKSwaps OBOMPKSwaps Swapping OOMP

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

%=====================================================================
% For an example on how to use the routine see  exa_OOMP_KSwaps
%=====================================================================

  name='KSWAPPING';
  %atoms having smaller norm than tol1 are supposed to be zero ones
  tol1=1e-7;
  %initilization
  if  nargin<9;  tol=1e-8 ;end
  if ((nargin<8) | (isempty(sws)==1)) ;sws=[];end
  if ((nargin<7) | (isempty(swi)==1)) ;swi=1;end
  if (sws <swi) error('swi should be < or = sws');end
  N=size(D,2); re=0;
  [L,k]=size(beta); delta=1;
  D0=D;Di0=Di;Q0=Q;beta0=beta;C0=C;k0=k;
  %start the swapping loop
  for s=swi:sws
    tic;
      if s>1; D=D0;Di=Di0;Q=Q0;beta=beta0;C=C0;k=k0;end
      nor1=norm(f'-D(:,1:k)*C')*sqrt(delta);  
      fprintf('Initial norm = %g\n',nor1);
      nor_old=nor1;
      nor1=nor1*nor1;
      nor2=0;z=0;
      %swap until improvement is possible
      while nor1>nor2
        if nor2~=0 nor1=nor2;else nor2=nor1;end
        backind=zeros(1,s);
        forwind=zeros(1,s);
        %============================
        %start the backward loop
        %============================
        for i=1:s
          dd=zeros(1,k);%clear dd
          for j=1:k  
            dd(j)=beta(:,j)'*beta(:,j);
          end	
          [Cp,p]=min((C.*conj(C))./dd);pp=Di(p);% p-th atom to be removed
          [D,Q,beta]=BioDictDel(D,Q,beta,p);
          dummy=Di(p);Di(p)=[];Di(N)=dummy;
          nor2=nor2+(Cp);
          backind(i)=pp;
          C=f*beta;
          k=k-1;
        end
        %============================================
        % start the forward loop
        %============================================
          for i=1:s
            c=zeros(1,N);cc=zeros(1,N);
            l=k+1;
            for j=l:N
              c(j)=norm(Q(:,j));
              if c(j)>tol1 cc(j)=abs(f*Q(:,j))/c(j);end;
            end
            [max_c,q]=max(cc);qq=Di(q);%q-th atom to be included
            %=======================================
            %if the contribution of the new atom is very small break it
            %=======================================
            if q~=l
              Q(:,[l q])=Q(:,[q l]); % swap l-th and q-th columns
              D(:,[l q])=D(:,[q l]);
              Di([l q])=Di([q l]);
            end
            [Q,beta]=BioDictIns(D,Q,beta,l);
            nor2=nor2-(abs(f*(beta(:,l)/norm(beta(:,l)))))^2;
            forwind(i)=qq;k=k+1;
            C=f*beta;res=norm(f'-D(:,1:k)*C');
            %=========================================================
            %if the required tolerance reached inside the forward loop,
            %break it
            %=========================================================
            if(res<tol) fprintf('\n Req. tol reached at forward step %d,error=%g\n',i,res); 
            re=1;
            break;
          end;
          end 
          if (re==1) break;end;
          if sort(backind)==sort(forwind) fprintf('Stop\n');break; end
          z=z+1;
          if nor1<=nor2 fprintf('\n nor1<=nor2\n');end  
          %fprintf('%d  %d~%d\n',pp,q,qq); 
          fprintf('z=%d, Norm %g\n',z, sqrt(nor2));
      end
      fprintf('\n%d Refinenments of %d functions lasted %g seconds (s=%d)\n',z,k,toc,s);
      fprintf('The original residuum norm was %g\n',nor_old); 
      fprintf('Residuum norm calculated additively after each iteration: %g\n',sqrt(nor2)*sqrt(delta));
      C=f*beta;res=norm(f'-D(:,1:k)*C');
      resid(s-swi+1)=res;
      fprintf('The final residuum norm was %g\n',res);  
      fprintf(['\nIf the previous two norms differ much, the ' ...
               'algorithm has not been too stable.\n']); 
      % biorthogonality
      ErrorTest(D(:,1:k),beta);  
      %orthogonality
      ErrorTest(Q(:,1:k));
      if (re==1) break;end;
      %==========================================================
      % check whether tolerance is reached or not
      %==========================================================
      if (res<tol) fprintf('Req. tol reached, error=%g\n',res); 
        re=1;
        break;
      end;
    end

%This function generalizes SWAPPING.m  copyright below. 
%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.
