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