annotate DL/RLS-DLA/SMALL_rlsdla1.m @ 45:b9465d2bb3b0

(none)
author idamnjanovic
date Mon, 14 Mar 2011 15:42:52 +0000
parents 6416fc12f2b8
children
rev   line source
idamnjanovic@40 1 function Dictionary = SMALL_rlsdla1(X, params)
idamnjanovic@40 2
idamnjanovic@40 3
idamnjanovic@40 4
idamnjanovic@40 5
idamnjanovic@40 6
idamnjanovic@40 7 global CODE_SPARSITY CODE_ERROR codemode
idamnjanovic@40 8 global MEM_LOW MEM_NORMAL MEM_HIGH memusage
idamnjanovic@40 9 global ompfunc ompparams exactsvd
idamnjanovic@40 10
idamnjanovic@40 11 CODE_SPARSITY = 1;
idamnjanovic@40 12 CODE_ERROR = 2;
idamnjanovic@40 13
idamnjanovic@40 14 MEM_LOW = 1;
idamnjanovic@40 15 MEM_NORMAL = 2;
idamnjanovic@40 16 MEM_HIGH = 3;
idamnjanovic@40 17
idamnjanovic@40 18
idamnjanovic@40 19 % p = randperm(size(X,2));
idamnjanovic@40 20
idamnjanovic@40 21 % coding mode %
idamnjanovic@40 22 X_norm=sqrt(sum(X.^2, 1));
idamnjanovic@40 23 % X_norm_1=sum(abs(X));
idamnjanovic@40 24 % X_norm_inf=max(abs(X));
idamnjanovic@40 25 [X_norm_sort, p]=sort(X_norm);%, 'descend');
idamnjanovic@40 26 % [X_norm_sort1, p5]=sort(X_norm_1);%, 'descend');
idamnjanovic@40 27
idamnjanovic@40 28 % if (isfield(params,'codemode'))
idamnjanovic@40 29 % switch lower(params.codemode)
idamnjanovic@40 30 % case 'sparsity'
idamnjanovic@40 31 % codemode = CODE_SPARSITY;
idamnjanovic@40 32 % thresh = params.Tdata;
idamnjanovic@40 33 % case 'error'
idamnjanovic@40 34 % codemode = CODE_ERROR;
idamnjanovic@40 35 % thresh = params.Edata;
idamnjanovic@40 36 % otherwise
idamnjanovic@40 37 % error('Invalid coding mode specified');
idamnjanovic@40 38 % end
idamnjanovic@40 39 % elseif (isfield(params,'Tdata'))
idamnjanovic@40 40 % codemode = CODE_SPARSITY;
idamnjanovic@40 41 % thresh = params.Tdata;
idamnjanovic@40 42 % elseif (isfield(params,'Edata'))
idamnjanovic@40 43 % codemode = CODE_ERROR;
idamnjanovic@40 44 % thresh = params.Edata;
idamnjanovic@40 45 %
idamnjanovic@40 46 % else
idamnjanovic@40 47 % error('Data sparse-coding target not specified');
idamnjanovic@40 48 % end
idamnjanovic@40 49
idamnjanovic@40 50 thresh = params.Edata;
idamnjanovic@40 51 % max number of atoms %
idamnjanovic@40 52
idamnjanovic@40 53 % if (codemode==CODE_ERROR && isfield(params,'maxatoms'))
idamnjanovic@40 54 % ompparams{end+1} = 'maxatoms';
idamnjanovic@40 55 % ompparams{end+1} = params.maxatoms;
idamnjanovic@40 56 % end
idamnjanovic@40 57
idamnjanovic@40 58
idamnjanovic@40 59 % memory usage %
idamnjanovic@40 60
idamnjanovic@40 61 if (isfield(params,'memusage'))
idamnjanovic@40 62 switch lower(params.memusage)
idamnjanovic@40 63 case 'low'
idamnjanovic@40 64 memusage = MEM_LOW;
idamnjanovic@40 65 case 'normal'
idamnjanovic@40 66 memusage = MEM_NORMAL;
idamnjanovic@40 67 case 'high'
idamnjanovic@40 68 memusage = MEM_HIGH;
idamnjanovic@40 69 otherwise
idamnjanovic@40 70 error('Invalid memory usage mode');
idamnjanovic@40 71 end
idamnjanovic@40 72 else
idamnjanovic@40 73 memusage = MEM_NORMAL;
idamnjanovic@40 74 end
idamnjanovic@40 75
idamnjanovic@40 76
idamnjanovic@40 77 % iteration count %
idamnjanovic@40 78
idamnjanovic@40 79 if (isfield(params,'iternum'))
idamnjanovic@40 80 iternum = params.iternum;
idamnjanovic@40 81 else
idamnjanovic@40 82 iternum = 10;
idamnjanovic@40 83 end
idamnjanovic@40 84
idamnjanovic@40 85
idamnjanovic@40 86 % omp function %
idamnjanovic@40 87
idamnjanovic@40 88 if (codemode == CODE_SPARSITY)
idamnjanovic@40 89 ompfunc = @omp;
idamnjanovic@40 90 else
idamnjanovic@40 91 ompfunc = @omp2;
idamnjanovic@40 92 end
idamnjanovic@40 93
idamnjanovic@40 94
idamnjanovic@40 95 % % status messages %
idamnjanovic@40 96 %
idamnjanovic@40 97 % printiter = 0;
idamnjanovic@40 98 % printreplaced = 0;
idamnjanovic@40 99 % printerr = 0;
idamnjanovic@40 100 % printgerr = 0;
idamnjanovic@40 101 %
idamnjanovic@40 102 % verbose = 't';
idamnjanovic@40 103 % msgdelta = -1;
idamnjanovic@40 104 %
idamnjanovic@40 105
idamnjanovic@40 106 %
idamnjanovic@40 107 % for i = 1:length(verbose)
idamnjanovic@40 108 % switch lower(verbose(i))
idamnjanovic@40 109 % case 'i'
idamnjanovic@40 110 % printiter = 1;
idamnjanovic@40 111 % case 'r'
idamnjanovic@40 112 % printiter = 1;
idamnjanovic@40 113 % printreplaced = 1;
idamnjanovic@40 114 % case 't'
idamnjanovic@40 115 % printiter = 1;
idamnjanovic@40 116 % printerr = 1;
idamnjanovic@40 117 % if (isfield(params,'testdata'))
idamnjanovic@40 118 % printgerr = 1;
idamnjanovic@40 119 % end
idamnjanovic@40 120 % end
idamnjanovic@40 121 % end
idamnjanovic@40 122 %
idamnjanovic@40 123 % if (msgdelta<=0 || isempty(verbose))
idamnjanovic@40 124 % msgdelta = -1;
idamnjanovic@40 125 % end
idamnjanovic@40 126 %
idamnjanovic@40 127 % ompparams{end+1} = 'messages';
idamnjanovic@40 128 % ompparams{end+1} = msgdelta;
idamnjanovic@40 129 %
idamnjanovic@40 130 %
idamnjanovic@40 131 %
idamnjanovic@40 132 % % compute error flag %
idamnjanovic@40 133 %
idamnjanovic@40 134 % comperr = (nargout>=3 || printerr);
idamnjanovic@40 135 %
idamnjanovic@40 136 %
idamnjanovic@40 137 % % validation flag %
idamnjanovic@40 138 %
idamnjanovic@40 139 % testgen = 0;
idamnjanovic@40 140 % if (isfield(params,'testdata'))
idamnjanovic@40 141 % testdata = params.testdata;
idamnjanovic@40 142 % if (nargout>=4 || printgerr)
idamnjanovic@40 143 % testgen = 1;
idamnjanovic@40 144 % end
idamnjanovic@40 145 % end
idamnjanovic@40 146
idamnjanovic@40 147 %
idamnjanovic@40 148 % % data norms %
idamnjanovic@40 149 %
idamnjanovic@40 150 % XtX = []; XtXg = [];
idamnjanovic@40 151 % if (codemode==CODE_ERROR && memusage==MEM_HIGH)
idamnjanovic@40 152 % XtX = colnorms_squared(data);
idamnjanovic@40 153 % if (testgen)
idamnjanovic@40 154 % XtXg = colnorms_squared(testdata);
idamnjanovic@40 155 % end
idamnjanovic@40 156 % end
idamnjanovic@40 157
idamnjanovic@40 158
idamnjanovic@40 159 % mutual incoherence limit %
idamnjanovic@40 160
idamnjanovic@40 161 if (isfield(params,'muthresh'))
idamnjanovic@40 162 muthresh = params.muthresh;
idamnjanovic@40 163 else
idamnjanovic@40 164 muthresh = 0.99;
idamnjanovic@40 165 end
idamnjanovic@40 166 if (muthresh < 0)
idamnjanovic@40 167 error('invalid muthresh value, must be non-negative');
idamnjanovic@40 168 end
idamnjanovic@40 169
idamnjanovic@40 170
idamnjanovic@40 171
idamnjanovic@40 172
idamnjanovic@40 173
idamnjanovic@40 174 % determine dictionary size %
idamnjanovic@40 175
idamnjanovic@40 176 if (isfield(params,'initdict'))
idamnjanovic@40 177 if (any(size(params.initdict)==1) && all(iswhole(params.initdict(:))))
idamnjanovic@40 178 dictsize = length(params.initdict);
idamnjanovic@40 179 else
idamnjanovic@40 180 dictsize = size(params.initdict,2);
idamnjanovic@40 181 end
idamnjanovic@40 182 end
idamnjanovic@40 183 if (isfield(params,'dictsize')) % this superceedes the size determined by initdict
idamnjanovic@40 184 dictsize = params.dictsize;
idamnjanovic@40 185 end
idamnjanovic@40 186
idamnjanovic@40 187 if (size(X,2) < dictsize)
idamnjanovic@40 188 error('Number of training signals is smaller than number of atoms to train');
idamnjanovic@40 189 end
idamnjanovic@40 190
idamnjanovic@40 191
idamnjanovic@40 192 % initialize the dictionary %
idamnjanovic@40 193
idamnjanovic@40 194 if (isfield(params,'initdict'))
idamnjanovic@40 195 if (any(size(params.initdict)==1) && all(iswhole(params.initdict(:))))
idamnjanovic@40 196 D1 = X(:,params.initdict(1:dictsize));
idamnjanovic@40 197 D2 = X(:,params.initdict(1:dictsize));
idamnjanovic@40 198 else
idamnjanovic@40 199 if (size(params.initdict,1)~=size(X,1) || size(params.initdict,2)<dictsize)
idamnjanovic@40 200 error('Invalid initial dictionary');
idamnjanovic@40 201 end
idamnjanovic@40 202 D1 = params.initdict(:,1:dictsize);
idamnjanovic@40 203 D2 = params.initdict(:,1:dictsize);
idamnjanovic@40 204 end
idamnjanovic@40 205 else
idamnjanovic@40 206 data_ids = find(colnorms_squared(X) > 1e-6); % ensure no zero data elements are chosen
idamnjanovic@40 207 perm = randperm(length(data_ids));
idamnjanovic@40 208 D = X(:,data_ids(perm(1:dictsize)));
idamnjanovic@40 209 end
idamnjanovic@40 210
idamnjanovic@40 211
idamnjanovic@40 212 % normalize the dictionary %
idamnjanovic@40 213
idamnjanovic@40 214 % D = normcols(D);
idamnjanovic@40 215 % DtD=D'*D;
idamnjanovic@40 216
idamnjanovic@40 217 err = zeros(1,iternum);
idamnjanovic@40 218 gerr = zeros(1,iternum);
idamnjanovic@40 219
idamnjanovic@40 220 if (codemode == CODE_SPARSITY)
idamnjanovic@40 221 errstr = 'RMSE';
idamnjanovic@40 222 else
idamnjanovic@40 223 errstr = 'mean atomnum';
idamnjanovic@40 224 end
idamnjanovic@40 225 X(:,p(X_norm_sort<thresh))=0;
idamnjanovic@40 226 % if (iternum==4)
idamnjanovic@40 227 % X_im=col2imstep(X, [256 256], [8 8]);
idamnjanovic@40 228 % else
idamnjanovic@40 229 % X_im=col2imstep(X, [512 512], [8 8]);
idamnjanovic@40 230 % end
idamnjanovic@40 231 % figure(10); imshow(X_im);
idamnjanovic@40 232 p1=p(X_norm_sort>thresh);
idamnjanovic@40 233
idamnjanovic@40 234 %p1=p1(p2(1:40000));
idamnjanovic@40 235 %end-min(40000, end)+1:end));%1:min(40000, end)));
idamnjanovic@40 236 %p1 = randperm(size(data,2));%size(data,2)
idamnjanovic@40 237 %data=data(:,p1);
idamnjanovic@40 238
idamnjanovic@40 239 C1=(100000*thresh)*eye(dictsize);
idamnjanovic@40 240 C2=(100000*thresh)*eye(dictsize);
idamnjanovic@40 241 % figure(11);
idamnjanovic@40 242 w=zeros(dictsize,1);
idamnjanovic@40 243 replaced=zeros(dictsize,1);
idamnjanovic@40 244 u=zeros(dictsize,1);
idamnjanovic@40 245 dictimg = showdict(D1,[8 8],round(sqrt(size(D1,2))),round(sqrt(size(D1,2))),'lines','highcontrast');
idamnjanovic@40 246 % h=imshow(imresize(dictimg,2,'nearest'));
idamnjanovic@40 247 lambda=0.9998
idamnjanovic@40 248 for j=1:3
idamnjanovic@40 249 if size(p1,2)>20000
idamnjanovic@40 250 p2 = randperm(floor(size(p1,2)/2));
idamnjanovic@40 251 p2=sort(p2(1:20000));
idamnjanovic@40 252 data1=X(:,p1(p2));
idamnjanovic@40 253 data2=X(:,p1(floor(size(p1,2)/2)+p2));
idamnjanovic@40 254 elseif size(p1,2)>0
idamnjanovic@40 255 data=X(:,p1);
idamnjanovic@40 256 else
idamnjanovic@40 257 break;
idamnjanovic@40 258 end
idamnjanovic@40 259 % figure(1);
idamnjanovic@40 260 % plot(sqrt(sum(data.^2, 1)));
idamnjanovic@40 261 % a=size(data,2)/4;
idamnjanovic@40 262 % lambda0=0.99;%1-16/numS+iternum*0.0001-0.0002
idamnjanovic@40 263 C1(1,1)=0;
idamnjanovic@40 264 C2(1,1)=0;
idamnjanovic@40 265 for i = 1:size(data1,2)
idamnjanovic@40 266 % if norm(data(:,i))>thresh
idamnjanovic@40 267 % par.multA= @(x,par) multMatr(D,x); % user function y=Ax
idamnjanovic@40 268 % par.multAt=@(x,par) multMatrAdj(D,x); % user function y=A'*x
idamnjanovic@40 269 % par.y=data(:,i);
idamnjanovic@40 270 % w=SolveFISTA(D,data(:,i),'lambda',0.5*thresh);
idamnjanovic@40 271 % w=sesoptn(zeros(dictsize,1),par.func_u, par.func_x, par.multA, par.multAt,options,par);
idamnjanovic@40 272 %w = SMALL_chol(D,data(:,i), 256,32, thresh);%
idamnjanovic@40 273 %w = sparsecode(data(:,i), D, [], [], thresh);
idamnjanovic@40 274 w1 = omp2mex(D1,data1(:,i),[],[],[],thresh,0,-1,-1,0);
idamnjanovic@40 275 w2 = omp2mex(D2,data2(:,i),[],[],[],thresh,0,-1,-1,0);
idamnjanovic@40 276 %w(find(w<1))=0;
idamnjanovic@40 277 %^2;
idamnjanovic@40 278 % lambda(i)=1-0.001/(1+i/a);
idamnjanovic@40 279 % if i<a
idamnjanovic@40 280 % lambda(i)=1-0.001*(1-(i/a));
idamnjanovic@40 281 % else
idamnjanovic@40 282 % lambda(i)=1;
idamnjanovic@40 283 % end
idamnjanovic@40 284 % param.lambda=thresh;
idamnjanovic@40 285 % param.mode=2;
idamnjanovic@40 286 % param.L=32;
idamnjanovic@40 287 % w=mexLasso(data(:,i), D, param);
idamnjanovic@40 288 spind1=find(w1);
idamnjanovic@40 289 spind2=find(w2);
idamnjanovic@40 290
idamnjanovic@40 291 %replaced(spind)=replaced(spind)+1;
idamnjanovic@40 292 %-0.001*(1/2)^(i/a);
idamnjanovic@40 293 % w_sp(i)=nnz(w);
idamnjanovic@40 294 residual1 = data1(:,i) - D1 * w1;
idamnjanovic@40 295 residual2 = data2(:,i) - D2 * w2;
idamnjanovic@40 296 %if ~isempty(spind)
idamnjanovic@40 297 %i
idamnjanovic@40 298
idamnjanovic@40 299 C1 = C1 *(1/ lambda);
idamnjanovic@40 300 C2 = C2 *(1/ lambda);
idamnjanovic@40 301 u1 = C1(:,spind1) * w1(spind1);
idamnjanovic@40 302 u2 = C2(:,spind2) * w2(spind2);
idamnjanovic@40 303 %spindu=find(u);
idamnjanovic@40 304 % v = D' * residual;
idamnjanovic@40 305
idamnjanovic@40 306 alfa1 = 1/(1 + w1' * u1);
idamnjanovic@40 307 alfa2 = 1/(1 + w2' * u2);
idamnjanovic@40 308 D1 = D1 + (alfa1 * residual1) * u1';
idamnjanovic@40 309 D2 = D2 + (alfa2 * residual2) * u2';
idamnjanovic@40 310 %uut=;
idamnjanovic@40 311 C1 = C1 - (alfa1 * u1)* u1';
idamnjanovic@40 312 C2 = C2 - (alfa2 * u2)* u2';
idamnjanovic@40 313 % lambda=(19*lambda+1)/20;
idamnjanovic@40 314 % DtD = DtD + alfa * ( v*u' + u*v') + alfa^2 * (residual'*residual) * uut;
idamnjanovic@40 315 % modi=5000;
idamnjanovic@40 316 % if (mod(i,modi)==0)
idamnjanovic@40 317 % Ximd=zeros(size(X));
idamnjanovic@40 318 % Ximd(:,p((i-modi+1:i)))=data(:,i-modi+1:i);
idamnjanovic@40 319 %
idamnjanovic@40 320 % if (iternum==4)
idamnjanovic@40 321 % X_ima=col2imstep(Ximd, [256 256], [8 8]);
idamnjanovic@40 322 % else
idamnjanovic@40 323 % X_ima=col2imstep(Ximd, [512 512], [8 8]);
idamnjanovic@40 324 % end
idamnjanovic@40 325 % dictimg1=dictimg;
idamnjanovic@40 326 % dictimg = showdict(D,[8 8],...
idamnjanovic@40 327 % round(sqrt(size(D,2))),round(sqrt(size(D,2))),'lines','highcontrast');
idamnjanovic@40 328 % dictimg1=(dictimg-dictimg1)*255;
idamnjanovic@40 329 %
idamnjanovic@40 330 % figure(2);
idamnjanovic@40 331 % subplot(2,2,1); imshow(X_ima);
idamnjanovic@40 332 % subplot(2,2,3); imshow(imresize(dictimg,2,'nearest'));
idamnjanovic@40 333 % subplot(2,2,4); imshow(imresize(dictimg1,2,'nearest'));
idamnjanovic@40 334 % subplot(2,2,2);imshow(C*(255/max(max(C))));
idamnjanovic@40 335 % pause(0.02);
idamnjanovic@40 336 % end
idamnjanovic@40 337 % end
idamnjanovic@40 338 end
idamnjanovic@40 339 %p1=p1(setxor(p2,1:end));
idamnjanovic@40 340 %[D,cleared_atoms] = cleardict(D,X,muthresh,p1,replaced);
idamnjanovic@40 341 %replaced=zeros(dictsize,1);
idamnjanovic@40 342 % W=sparsecode(data, D, [], [], thresh);
idamnjanovic@40 343 % data=D*W;
idamnjanovic@40 344 lambda=lambda+0.0001
idamnjanovic@40 345 end
idamnjanovic@40 346 %Gamma=mexLasso(data, D, param);
idamnjanovic@40 347 %err=compute_err(D,Gamma, data);
idamnjanovic@40 348 %[y,i]=max(err);
idamnjanovic@40 349 %D(:,1)=data(:,i)/norm(data(:,i));
idamnjanovic@40 350 % D=normcols(D);
idamnjanovic@40 351 % D_norm=sqrt(sum(D.^2, 1));
idamnjanovic@40 352 % D_norm_1=sum(abs(D));
idamnjanovic@40 353 % X_norm_1=sum(abs(X));
idamnjanovic@40 354 % X_norm_inf=max(abs(X));
idamnjanovic@40 355 % [D_norm_sort, p]=sort(D_norm_1, 'descend');
idamnjanovic@40 356 Dictionary =[D1 D2];
idamnjanovic@40 357 % figure(3);
idamnjanovic@40 358 % plot(lambda);
idamnjanovic@40 359 % mean(lambda);
idamnjanovic@40 360 % figure(4+j);plot(w_sp);
idamnjanovic@40 361 end
idamnjanovic@40 362
idamnjanovic@40 363 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
idamnjanovic@40 364 % sparsecode %
idamnjanovic@40 365 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
idamnjanovic@40 366
idamnjanovic@40 367 function Gamma = sparsecode(data,D,XtX,G,thresh)
idamnjanovic@40 368
idamnjanovic@40 369 global CODE_SPARSITY codemode
idamnjanovic@40 370 global MEM_HIGH memusage
idamnjanovic@40 371 global ompfunc ompparams
idamnjanovic@40 372
idamnjanovic@40 373 if (memusage < MEM_HIGH)
idamnjanovic@40 374 Gamma = ompfunc(D,data,G,thresh,ompparams{:});
idamnjanovic@40 375
idamnjanovic@40 376 else % memusage is high
idamnjanovic@40 377
idamnjanovic@40 378 if (codemode == CODE_SPARSITY)
idamnjanovic@40 379 Gamma = ompfunc(D'*data,G,thresh,ompparams{:});
idamnjanovic@40 380
idamnjanovic@40 381 else
idamnjanovic@40 382 Gamma = ompfunc(D, data, G, thresh,ompparams{:});
idamnjanovic@40 383 end
idamnjanovic@40 384
idamnjanovic@40 385 end
idamnjanovic@40 386
idamnjanovic@40 387 end
idamnjanovic@40 388
idamnjanovic@40 389
idamnjanovic@40 390 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
idamnjanovic@40 391 % compute_err %
idamnjanovic@40 392 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
idamnjanovic@40 393
idamnjanovic@40 394
idamnjanovic@40 395 function err = compute_err(D,Gamma,data)
idamnjanovic@40 396
idamnjanovic@40 397 global CODE_SPARSITY codemode
idamnjanovic@40 398
idamnjanovic@40 399 if (codemode == CODE_SPARSITY)
idamnjanovic@40 400 err = sqrt(sum(reperror2(data,D,Gamma))/numel(data));
idamnjanovic@40 401 else
idamnjanovic@40 402 err = nnz(Gamma)/size(data,2);
idamnjanovic@40 403 end
idamnjanovic@40 404
idamnjanovic@40 405 end
idamnjanovic@40 406
idamnjanovic@40 407
idamnjanovic@40 408
idamnjanovic@40 409 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
idamnjanovic@40 410 % cleardict %
idamnjanovic@40 411 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
idamnjanovic@40 412
idamnjanovic@40 413
idamnjanovic@40 414 function [D,cleared_atoms] = cleardict(D,X,muthresh,unused_sigs,replaced_atoms)
idamnjanovic@40 415
idamnjanovic@40 416 use_thresh = 4; % at least this number of samples must use the atom to be kept
idamnjanovic@40 417
idamnjanovic@40 418 dictsize = size(D,2);
idamnjanovic@40 419
idamnjanovic@40 420 % compute error in blocks to conserve memory
idamnjanovic@40 421 % err = zeros(1,size(X,2));
idamnjanovic@40 422 % blocks = [1:3000:size(X,2) size(X,2)+1];
idamnjanovic@40 423 % for i = 1:length(blocks)-1
idamnjanovic@40 424 % err(blocks(i):blocks(i+1)-1) = sum((X(:,blocks(i):blocks(i+1)-1)-D*Gamma(:,blocks(i):blocks(i+1)-1)).^2);
idamnjanovic@40 425 % end
idamnjanovic@40 426
idamnjanovic@40 427 cleared_atoms = 0;
idamnjanovic@40 428 usecount = replaced_atoms;%sum(abs(Gamma)>1e-7, 2);
idamnjanovic@40 429
idamnjanovic@40 430 for j = 1:dictsize
idamnjanovic@40 431
idamnjanovic@40 432 % compute G(:,j)
idamnjanovic@40 433 Gj = D'*D(:,j);
idamnjanovic@40 434 Gj(j) = 0;
idamnjanovic@40 435
idamnjanovic@40 436 % replace atom
idamnjanovic@40 437 if ( (max(Gj.^2)>muthresh^2 || usecount(j)<use_thresh) && ~replaced_atoms(j) )
idamnjanovic@40 438 % [y,i] = max(err(unused_sigs));
idamnjanovic@40 439 D(:,j) = X(:,unused_sigs(end)) / norm(X(:,unused_sigs(end)));
idamnjanovic@40 440 unused_sigs = unused_sigs([1:end-1]);
idamnjanovic@40 441 cleared_atoms = cleared_atoms+1;
idamnjanovic@40 442 end
idamnjanovic@40 443 end
idamnjanovic@40 444
idamnjanovic@40 445 end
idamnjanovic@40 446
idamnjanovic@40 447
idamnjanovic@40 448
idamnjanovic@40 449 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
idamnjanovic@40 450 % misc functions %
idamnjanovic@40 451 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
idamnjanovic@40 452
idamnjanovic@40 453
idamnjanovic@40 454 function err2 = reperror2(X,D,Gamma)
idamnjanovic@40 455
idamnjanovic@40 456 % compute in blocks to conserve memory
idamnjanovic@40 457 err2 = zeros(1,size(X,2));
idamnjanovic@40 458 blocksize = 2000;
idamnjanovic@40 459 for i = 1:blocksize:size(X,2)
idamnjanovic@40 460 blockids = i : min(i+blocksize-1,size(X,2));
idamnjanovic@40 461 err2(blockids) = sum((X(:,blockids) - D*Gamma(:,blockids)).^2);
idamnjanovic@40 462 end
idamnjanovic@40 463
idamnjanovic@40 464 end
idamnjanovic@40 465
idamnjanovic@40 466
idamnjanovic@40 467 function Y = colnorms_squared(X)
idamnjanovic@40 468
idamnjanovic@40 469 % compute in blocks to conserve memory
idamnjanovic@40 470 Y = zeros(1,size(X,2));
idamnjanovic@40 471 blocksize = 2000;
idamnjanovic@40 472 for i = 1:blocksize:size(X,2)
idamnjanovic@40 473 blockids = i : min(i+blocksize-1,size(X,2));
idamnjanovic@40 474 Y(blockids) = sum(X(:,blockids).^2);
idamnjanovic@40 475 end
idamnjanovic@40 476
idamnjanovic@40 477 end
idamnjanovic@40 478
idamnjanovic@40 479