view DL/RLS-DLA/SMALL_rlsdla 05032011.m @ 77:62f20b91d870

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

function Dictionary = SMALL_rlsdla(X, params)





global CODE_SPARSITY CODE_ERROR codemode
global MEM_LOW MEM_NORMAL MEM_HIGH memusage
global ompfunc ompparams exactsvd

CODE_SPARSITY = 1;
CODE_ERROR = 2;

MEM_LOW = 1;
MEM_NORMAL = 2;
MEM_HIGH = 3;


% p = randperm(size(X,2));

        % coding mode %
%X_norm=sqrt(sum(X.^2, 1));
% X_norm_1=sum(abs(X));
% X_norm_inf=max(abs(X));
%[X_norm_sort, p]=sort(X_norm);%, 'descend');
% [X_norm_sort1, p5]=sort(X_norm_1);%, 'descend');

%         if (isfield(params,'codemode'))
%           switch lower(params.codemode)
%             case 'sparsity'
%               codemode = CODE_SPARSITY;
%               thresh = params.Tdata;
%             case 'error'
%               codemode = CODE_ERROR;
%               thresh = params.Edata;
%             otherwise
%               error('Invalid coding mode specified');
%           end
%         elseif (isfield(params,'Tdata'))
%           codemode = CODE_SPARSITY;
%           thresh = params.Tdata;
%         elseif (isfield(params,'Edata'))
%           codemode = CODE_ERROR;
%           thresh = params.Edata;
% 
%         else
%           error('Data sparse-coding target not specified');
%         end

thresh = params.Edata;
% max number of atoms %

% if (codemode==CODE_ERROR && isfield(params,'maxatoms'))
%   ompparams{end+1} = 'maxatoms';
%   ompparams{end+1} = params.maxatoms;
% end


% memory usage %

if (isfield(params,'memusage'))
  switch lower(params.memusage)
    case 'low'
      memusage = MEM_LOW;
    case 'normal'
      memusage = MEM_NORMAL;
    case 'high'
      memusage = MEM_HIGH;
    otherwise
      error('Invalid memory usage mode');
  end
else
  memusage = MEM_NORMAL;
end


% iteration count %

if (isfield(params,'iternum'))
  iternum = params.iternum;
else
  iternum = 10;
end


% omp function %

if (codemode == CODE_SPARSITY)
  ompfunc = @omp;
else
  ompfunc = @omp2;
end


% % status messages %
% 
% printiter = 0;
% printreplaced = 0;
% printerr = 0;
% printgerr = 0;
% 
% verbose = 't';
% msgdelta = -1;
% 

% 
% for i = 1:length(verbose)
%   switch lower(verbose(i))
%     case 'i'
%       printiter = 1;
%     case 'r'
%       printiter = 1;
%       printreplaced = 1;
%     case 't'
%       printiter = 1;
%       printerr = 1;
%       if (isfield(params,'testdata'))
%         printgerr = 1;
%       end
%   end
% end
% 
% if (msgdelta<=0 || isempty(verbose))
%   msgdelta = -1; 
% end
% 
% ompparams{end+1} = 'messages';
% ompparams{end+1} = msgdelta;
% 
% 
% 
% % compute error flag %
% 
% comperr = (nargout>=3 || printerr);
% 
% 
% % validation flag %
% 
% testgen = 0;
% if (isfield(params,'testdata'))
%   testdata = params.testdata;
%   if (nargout>=4 || printgerr)
%     testgen = 1;
%   end
% end

% 
% % data norms %
% 
% XtX = []; XtXg = [];
% if (codemode==CODE_ERROR && memusage==MEM_HIGH)
%   XtX = colnorms_squared(data);
%   if (testgen)
%     XtXg = colnorms_squared(testdata);
%   end
% end


% mutual incoherence limit %

if (isfield(params,'muthresh'))
  muthresh = params.muthresh;
else
  muthresh = 0.99;
end
if (muthresh < 0)
  error('invalid muthresh value, must be non-negative');
end





% determine dictionary size %

if (isfield(params,'initdict'))
  if (any(size(params.initdict)==1) && all(iswhole(params.initdict(:))))
    dictsize = length(params.initdict);
  else
    dictsize = size(params.initdict,2);
  end
end
if (isfield(params,'dictsize'))    % this superceedes the size determined by initdict
  dictsize = params.dictsize;
end

if (size(X,2) < dictsize)
  error('Number of training signals is smaller than number of atoms to train');
end


% initialize the dictionary %

if (isfield(params,'initdict'))
  if (any(size(params.initdict)==1) && all(iswhole(params.initdict(:))))
    D = X(:,params.initdict(1:dictsize));
  else
    if (size(params.initdict,1)~=size(X,1) || size(params.initdict,2)<dictsize)
      error('Invalid initial dictionary');
    end
    D = params.initdict(:,1:dictsize);
  end
else
  data_ids = find(colnorms_squared(X) > 1e-6);   % ensure no zero data elements are chosen
  perm = randperm(length(data_ids));
  D = X(:,data_ids(perm(1:dictsize)));
end


% normalize the dictionary %

% D = normcols(D);
% DtD=D'*D;

err = zeros(1,iternum);
gerr = zeros(1,iternum);

if (codemode == CODE_SPARSITY)
  errstr = 'RMSE';
else
  errstr = 'mean atomnum';
end
%X(:,p(X_norm_sort<thresh))=0;
% if (iternum==4)
%     X_im=col2imstep(X, [256 256], [8 8]);
% else
%     X_im=col2imstep(X, [512 512], [8 8]);
% end
% figure(10); imshow(X_im);

%p1=p(cumsum(X_norm_sort)./[1:size(X_norm_sort,2)]>thresh);
%p1=p(X_norm_sort>thresh);
% X(:,setxor(p1,1:end))=0;
% X_im=col2imstep(X, [256 256], [8 8]);
% figure(10); imshow(X_im);
% if iternum==2
% D(:,1)=D(:,2);
% end
%p1=p1(p2(1:40000));
%end-min(40000, end)+1:end));%1:min(40000, end)));
%p1 = randperm(size(data,2));%size(data,2)
%data=data(:,p1);

C=(100000*thresh)*eye(dictsize);
% figure(11);
w=zeros(dictsize,1);
replaced=zeros(dictsize,1);
u=zeros(dictsize,1);
%   dictimg = showdict(D,[8 8],round(sqrt(size(D,2))),round(sqrt(size(D,2))),'lines','highcontrast');
%   figure(11);imshow(imresize(dictimg,2,'nearest'));
%  pause(1);
lambda=0.9997%0.99986;%3+0.0001*params.linc;
for j=1:1
    %data=X;
if size(X,2)>40000
    p2 = randperm(size(X,2));
    
    p2=sort(p2(1:40000));%min(floor(size(p1,2)/2),40000)));
    size(p2,2)
    data=X(:,p2);
elseif size(X,2)>0
    %p2 = randperm(size(p1,2));
    size(X,2)
    data=X;
else
    break;
end
% figure(1);
% plot(sqrt(sum(data.^2, 1)));
% a=size(data,2)/4;
% lambda0=0.99;%1-16/numS+iternum*0.0001-0.0002
%C(1,1)=0;
modi=1000;
for i = 1:size(data,2)
%     if norm(data(:,i))>thresh
    %      par.multA= @(x,par)  multMatr(D,x);     % user function   y=Ax
    %      par.multAt=@(x,par)  multMatrAdj(D,x);  % user function  y=A'*x
    %      par.y=data(:,i);
     %  w=SolveFISTA(D,data(:,i),'lambda',0.5*thresh);
     % w=sesoptn(zeros(dictsize,1),par.func_u, par.func_x, par.multA, par.multAt,options,par);
   %w = SMALL_chol(D,data(:,i), 256,32, thresh);%
   %w = sparsecode(data(:,i), D, [], [], thresh);
   w = omp2mex(D,data(:,i),[],[],[],thresh,0,-1,-1,0);

   %w(find(w<1))=0; 
   %^2;
%   lambda(i)=1-0.001/(1+i/a);
% if i<a
%     lambda(i)=1-0.001*(1-(i/a));
% else
%     lambda(i)=1;
% end
%    param.lambda=thresh;
%    param.mode=2;
%    param.L=32;
%    w=mexLasso(data(:,i), D, param);
    spind=find(w);
    %replaced(spind)=replaced(spind)+1;
    %-0.001*(1/2)^(i/a);
%   w_sp(i)=nnz(w);
    residual = data(:,i) - D * w;
   %if ~isempty(spind) 
    %i
   if (j==1) 
    C = C *(1/ lambda);
   end
    u = C(:,spind) * w(spind);
  
   %spindu=find(u);
   % v = D' * residual;
    
    alfa = 1/(1 + w' * u);
    
    D = D + (alfa * residual) * u';
    
    %uut=;
    C = C - (alfa * u)* u';
   % lambda=(19*lambda+1)/20;
   % DtD = DtD + alfa * ( v*u' + u*v') + alfa^2 * (residual'*residual) * uut;
   
%    if (mod(i,modi)==0)
%        Ximd=zeros(size(X));
%        Ximd(:,p1((i-modi+1:i)))=data(:,i-modi+1:i);
%        
%        if (iternum==4)
%             X_ima(:,:,1)=col2imstep(Ximd, [256 256], [8 8]);
%             X_ima(:,:,2)=col2imstep(X, [256 256], [8 8]);
%             X_ima(:,:,3)=zeros(256,256);
%        else
%             X_ima(:,:,1)=col2imstep(Ximd, [512 512], [8 8]);
%             X_ima(:,:,2)=col2imstep(X, [512 512], [8 8]);
%             X_ima(:,:,3)=zeros(512,512);
%        end
%        
%        dictimg1=dictimg;
%        dictimg = showdict(D,[8 8],...
%         round(sqrt(size(D,2))),round(sqrt(size(D,2))),'lines','highcontrast');
%         dictimg1=(dictimg-dictimg1);
%         
%         figure(2);
%         subplot(2,2,1); imshow(X_ima); title(sprintf('%d',i));
%         subplot(2,2,3); imshow(imresize(dictimg,2,'nearest'));
%         subplot(2,2,4);    imshow(imresize(dictimg1,2,'nearest'));
%         subplot(2,2,2);imshow(C*(255/max(max(C))));
%         pause(0.02);
%         if (i>=35000) 
%             modi=100;
%             pause 
%         end;
%     end
%    end
end
%p1=p1(setxor(p2,1:end));
%[D,cleared_atoms] = cleardict(D,X,muthresh,p1,replaced);
%replaced=zeros(dictsize,1);
% W=sparsecode(data, D, [], [], thresh);
% data=D*W;
%lambda=lambda+0.0002
end
%Gamma=mexLasso(data, D, param);
%err=compute_err(D,Gamma, data);
%[y,i]=max(err);
%D(:,1)=data(:,i)/norm(data(:,i));

Dictionary = D;%D(:,p);
% figure(3);
% plot(lambda);
% mean(lambda);
% figure(4+j);plot(w_sp);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             sparsecode               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Gamma = sparsecode(data,D,XtX,G,thresh)

global CODE_SPARSITY codemode
global MEM_HIGH memusage
global ompfunc ompparams

if (memusage < MEM_HIGH)
  Gamma = ompfunc(D,data,G,thresh,ompparams{:});
  
else  % memusage is high
  
  if (codemode == CODE_SPARSITY)
    Gamma = ompfunc(D'*data,G,thresh,ompparams{:});
    
  else
    Gamma = ompfunc(D, data, G, thresh,ompparams{:});
  end
  
end

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             compute_err              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function err = compute_err(D,Gamma,data)
  
global CODE_SPARSITY codemode

if (codemode == CODE_SPARSITY)
  err = sqrt(sum(reperror2(data,D,Gamma))/numel(data));
else
  err = nnz(Gamma)/size(data,2);
end

end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           cleardict                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [D,cleared_atoms] = cleardict(D,X,muthresh,unused_sigs,replaced_atoms)

use_thresh = 4;  % at least this number of samples must use the atom to be kept

dictsize = size(D,2);

% compute error in blocks to conserve memory
% err = zeros(1,size(X,2));
% blocks = [1:3000:size(X,2) size(X,2)+1];
% for i = 1:length(blocks)-1
%   err(blocks(i):blocks(i+1)-1) = sum((X(:,blocks(i):blocks(i+1)-1)-D*Gamma(:,blocks(i):blocks(i+1)-1)).^2);
% end

cleared_atoms = 0;
usecount = replaced_atoms;%sum(abs(Gamma)>1e-7, 2);

for j = 1:dictsize
  
  % compute G(:,j)
  Gj = D'*D(:,j);
  Gj(j) = 0;
  
  % replace atom
  if ( (max(Gj.^2)>muthresh^2 || usecount(j)<use_thresh) && ~replaced_atoms(j) )
%     [y,i] = max(err(unused_sigs));
    D(:,j) = X(:,unused_sigs(end)) / norm(X(:,unused_sigs(end)));
    unused_sigs = unused_sigs([1:end-1]);
    cleared_atoms = cleared_atoms+1;
  end
end

end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%            misc functions            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function err2 = reperror2(X,D,Gamma)

% compute in blocks to conserve memory
err2 = zeros(1,size(X,2));
blocksize = 2000;
for i = 1:blocksize:size(X,2)
  blockids = i : min(i+blocksize-1,size(X,2));
  err2(blockids) = sum((X(:,blockids) - D*Gamma(:,blockids)).^2);
end

end


function Y = colnorms_squared(X)

% compute in blocks to conserve memory
Y = zeros(1,size(X,2));
blocksize = 2000;
for i = 1:blocksize:size(X,2)
  blockids = i : min(i+blocksize-1,size(X,2));
  Y(blockids) = sum(X(:,blockids).^2);
end

end