idamnjanovic@3: function [A, resF]=SMALL_cgp(Dict,X, m,  maxNumCoef, errorGoal, varargin) 
ivan@128: %%  Partial Conjugate Gradient Pursuit Multiple vectors sparse representation
ivan@128: %   
ivan@128: %   Sparse coding of a group of signals based on a given 
ivan@128: %   dictionary and specified number of atoms to use. 
ivan@128: %   input arguments: Dict - the dictionary
ivan@128: %                  X - the signals to represent
ivan@128: %                  m - number of atoms in Dictionary
ivan@128: %                  errorGoal - the maximal allowed representation error for
ivan@128: %                  each signal.
ivan@128: %
ivan@128: %   optional:        if Dict is function handle then Transpose Dictionary
ivan@128: %                  handle needs to be specified.
ivan@128: %
ivan@128: %   output arguments: A - sparse coefficient matrix.
ivan@128: 
idamnjanovic@22: %
idamnjanovic@22: %   Centre for Digital Music, Queen Mary, University of London.
idamnjanovic@22: %   This file copyright 2009 Ivan Damnjanovic.
idamnjanovic@22: %
idamnjanovic@22: %   This program is free software; you can redistribute it and/or
idamnjanovic@22: %   modify it under the terms of the GNU General Public License as
idamnjanovic@22: %   published by the Free Software Foundation; either version 2 of the
idamnjanovic@22: %   License, or (at your option) any later version.  See the file
idamnjanovic@22: %   COPYING included with this distribution for more information.
idamnjanovic@1: %
idamnjanovic@1: %%
idamnjanovic@1: 
idamnjanovic@1: % This Dictionary check is based on Thomas Blumensath work in sparsify 0_4 greedy solvers 
idamnjanovic@1: if     isa(Dict,'float')      
idamnjanovic@1:             D =@(z) Dict*z;  
idamnjanovic@1:             Dt =@(z) Dict'*z;
idamnjanovic@1: elseif isobject(Dict) 
idamnjanovic@1:             D =@(z) Dict*z;  
idamnjanovic@1:             Dt =@(z) Dict'*z;
idamnjanovic@1: elseif isa(Dict,'function_handle') 
idamnjanovic@1:     try
idamnjanovic@1:         DictT=varargin{1};
idamnjanovic@1:         if isa(DictT,'function_handle'); 
idamnjanovic@73:                 D=Dict; 
idamnjanovic@73:                 Dt=DictT;
idamnjanovic@1:         else
idamnjanovic@1:                 error('If Dictionary is a function handle,Transpose Dictionary also needs to be a function handle. ');
idamnjanovic@1:         end
idamnjanovic@1:     catch
idamnjanovic@1:         error('If Dictionary is a function handle, Transpose Dictionary  needs to be specified. Exiting.');
idamnjanovic@1:     end
idamnjanovic@1: else
idamnjanovic@1:     error('Dictionary is of unsupported type. Use explicit matrix, function_handle or object. Exiting.');
idamnjanovic@1: end
idamnjanovic@1: %%
idamnjanovic@62: %positivity=1;
idamnjanovic@1: [n,P]=size(X);
idamnjanovic@1: 
idamnjanovic@1: 
idamnjanovic@1: 
idamnjanovic@1: A = sparse(m,size(X,2));
idamnjanovic@1: 
idamnjanovic@1: for k=1:1:P,
idamnjanovic@1:     
idamnjanovic@1:     x = X(:,k);
idamnjanovic@1:     residual = x;
idamnjanovic@1: 	indx =[];
idamnjanovic@1: 	j=0;
idamnjanovic@1:    
idamnjanovic@1: 	
idamnjanovic@1:     currResNorm = residual'*residual/n;
idamnjanovic@3:     errGoal=errorGoal*currResNorm;
idamnjanovic@1:     a = zeros(m,1);
idamnjanovic@1:     p = zeros(m,1);
idamnjanovic@1:     
idamnjanovic@1:     while  j < maxNumCoef,
idamnjanovic@1:         
idamnjanovic@1:         j = j+1;
idamnjanovic@1:         dir=Dt(residual);
idamnjanovic@22:         if exist('positivity','var')&&(positivity==1)
idamnjanovic@22:             [tmp__, pos]=max(dir);
idamnjanovic@22:         else
idamnjanovic@22:             [tmp__, pos]=max(abs(dir));
idamnjanovic@22:         end
idamnjanovic@1:         indx(j)=pos;
idamnjanovic@1:         
idamnjanovic@1:         p(indx)=dir(indx);
idamnjanovic@1:         Dp=D(p);
idamnjanovic@1:         pDDp=Dp'*Dp;
idamnjanovic@1:         
idamnjanovic@1:         if (j==1) 
idamnjanovic@1:             beta=0;
idamnjanovic@1:         else
idamnjanovic@1:             beta=-Dp'*residual/pDDp;
idamnjanovic@1:         end
idamnjanovic@1:         
idamnjanovic@1:         alfa=residual'*Dp/pDDp;
idamnjanovic@1:         a=a+alfa*p;
idamnjanovic@1:         p(indx)=dir(indx)+beta*p(indx);
idamnjanovic@1:         
idamnjanovic@1:         residual=residual-alfa*Dp;
idamnjanovic@1:     
idamnjanovic@1: 		currResNorm = residual'*residual/n;
idamnjanovic@3:         if currResNorm<errGoal
idamnjanovic@1:             fprintf('\nFound exact representation! \n');
idamnjanovic@1:             break;
idamnjanovic@1:         end
idamnjanovic@1:    end;
idamnjanovic@1:    if (~isempty(indx))
idamnjanovic@3:        resF(k)=currResNorm;
idamnjanovic@1:        A(indx,k)=a(indx);
idamnjanovic@1:    end
idamnjanovic@1: end;
idamnjanovic@1: return;
idamnjanovic@1: 
idamnjanovic@1: 
idamnjanovic@1: