view solvers/SMALL_cgp.m @ 77:62f20b91d870

add routines from sparco problems privite folder to {root}\util
some changes to ksvd vs rlsdla image denoising example

function [A, resF]=SMALL_cgp(Dict,X, m,  maxNumCoef, errorGoal, varargin) 
%%% Conjugate Gradient Pursuit Multiple vectors sparse representation
%
%   Centre for Digital Music, Queen Mary, University of London.
%   This file copyright 2009 Ivan Damnjanovic.
%
%   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.  See the file
%   COPYING included with this distribution for more information.
%   
% Sparse coding of a group of signals based on a given 
% dictionary and specified number of atoms to use. 
% input arguments: Dict - the dictionary
%                  X - the signals to represent
%                  m - number of atoms in Dictionary
%                  errorGoal - the maximal allowed representation error for
%                  each signal.
%
% optional:        if Dict is function handle then Transpose Dictionary
%                  handle needs to be specified.
%
% output arguments: A - sparse coefficient matrix.
%
%%

% This Dictionary check is based on Thomas Blumensath work in sparsify 0_4 greedy solvers 
if     isa(Dict,'float')      
            D =@(z) Dict*z;  
            Dt =@(z) Dict'*z;
elseif isobject(Dict) 
            D =@(z) Dict*z;  
            Dt =@(z) Dict'*z;
elseif isa(Dict,'function_handle') 
    try
        DictT=varargin{1};
        if isa(DictT,'function_handle'); 
                D=Dict; 
                Dt=DictT;
        else
                error('If Dictionary is a function handle,Transpose Dictionary also needs to be a function handle. ');
        end
    catch
        error('If Dictionary is a function handle, Transpose Dictionary  needs to be specified. Exiting.');
    end
else
    error('Dictionary is of unsupported type. Use explicit matrix, function_handle or object. Exiting.');
end
%%
%positivity=1;
[n,P]=size(X);



A = sparse(m,size(X,2));

for k=1:1:P,
    
    x = X(:,k);
    residual = x;
	indx =[];
	j=0;
   
	
    currResNorm = residual'*residual/n;
    errGoal=errorGoal*currResNorm;
    a = zeros(m,1);
    p = zeros(m,1);
    
    while  j < maxNumCoef,
        
        j = j+1;
        dir=Dt(residual);
        if exist('positivity','var')&&(positivity==1)
            [tmp__, pos]=max(dir);
        else
            [tmp__, pos]=max(abs(dir));
        end
        indx(j)=pos;
        
        p(indx)=dir(indx);
        Dp=D(p);
        pDDp=Dp'*Dp;
        
        if (j==1) 
            beta=0;
        else
            beta=-Dp'*residual/pDDp;
        end
        
        alfa=residual'*Dp/pDDp;
        a=a+alfa*p;
        p(indx)=dir(indx)+beta*p(indx);
        
        residual=residual-alfa*Dp;
    
		currResNorm = residual'*residual/n;
        if currResNorm<errGoal
            fprintf('\nFound exact representation! \n');
            break;
        end
   end;
   if (~isempty(indx))
       resF(k)=currResNorm;
       A(indx,k)=a(indx);
   end
end;
return;