comparison solvers/SMALL_pcgp.m @ 163:855025f4c779 ivand_dev

renaiming small_cgp to small_pcgp

+function [A, resF]=SMALL_pcgp(Dict,X, m,  maxNumCoef, errorGoal, varargin) 
+%%  Partial Conjugate Gradient Pursuit Multiple vectors sparse representation
+%   
+%   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.
+
+%
+%   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.
+%
+%%
+
+% 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;
+
+
+