script_simulation_adpWBPDN

 partial script for running different solvers for adaptive reweighting

0001 % partial script for running different solvers for adaptive reweighting
0002 
0003 % reweighted setup
0004 adp_wt = 0;  % adaptive selection of the weights?
0005 rwt = 3*(1-adp_wt)+1;        % number of reweighting iterations
0006 
0007 I_ALL = [];
0008 EXP = [0 1; 0 4; 1 1]; 
0009 % 0 1 -- standard BPDN
0010 % 0 x -- iterative reweighting with x reweighting iterations
0011 % 1 1 -- adaptive reweigthing
0012 for exp = 1:3;
0013     % set up random generator
0014     rseed = 2012;
0015     rand('state',rseed);
0016     randn('state',rseed);
0017 
0018     % RandStream.setDefaultStream(RandStream('mt19937ar','seed',rseed));
0019     % seed = rng('shuffle');
0020     % seed = rng(0);
0021 
0022     adp_wt = EXP(exp,1);
0023     rwt = EXP(exp,2);
0024     
0025     fprintf('adp_wt = %d, rwt = %d ',adp_wt,rwt);
0026     
0027     
0028     err_fun = @(z) -20*log10(z);
0029     
0030     %% Main loop
0031     SIM_stack = cell(maxsim,1);
0032     SIM_memory = cell(maxsim,1);
0033     
0034     Ir_stack = [];
0035     for sim = 1:maxsim
0036         %% Generate a random signal
0037         in = []; in.type = sType; in.T = T; in.randgen = 1; in.take_fwt = 1;
0038         [x img wave_struct] = genSignal(N,in);
0039         % gamma_orig = find(abs(x)>0);
0040         [val ind] = sort(abs(x),'descend');
0041         ind_pos = ind(val>0);
0042         gamma_orig = ind_pos(1:min(length(ind_pos),M-1));
0043         
0044         % measurement matrix
0045         in = []; in.type = mType;
0046         AA = genAmat(M,N,in);
0047         if largescale
0048             wType = wave_struct.wType;
0049             FilterBank = wave_struct.FilterBank;
0050             h0 = FilterBank(1,:);
0051             h1 = FilterBank(2,:);
0052             g0 = FilterBank(3,:);
0053             g1 = FilterBank(4,:);
0054             J = wave_struct.J;
0055             sym = wave_struct.sym;
0056             
0057             n = sqrt(N);
0058             rshp1 = @(z) reshape(z,n,n);
0059             rshp2 = @(z) z(:);
0060             
0061             % FT = @(x) rshp2(fwt2(rshp1(x),h0,h1,J,sym));
0062             % iFT = @(x) rshp2(ifwt2(rshp1(x),g0,g1,J,sym));
0063             FT = @(z) z; 
0064             iFT = @(z) z; 
0065             
0066             A = @(x) AA(iFT(x),1);
0067             AT = @(x) FT(AA(x,2));
0068             % measurements
0069             Ax = A(x);
0070             sigma = sqrt(norm(Ax)^2/10^(SNR/10)/M);
0071             %sigma = .05;
0072             e = randn(M,1)*sigma;
0073             
0074             y = Ax+e;
0075             Aty = AT(y);
0076         else
0077             A = AA;
0078             AT = A';
0079             % measurements
0080             Ax = A*x;
0081             sigma = sqrt(norm(Ax)^2/10^(SNR/10)/M);
0082             %sigma = .05;
0083             e = randn(M,1)*sigma;
0084             y = A*x+e;
0085             Aty = AT*y;
0086         end
0087         
0088         maxiter = 2*N;
0089         
0090         %% wSPGL1
0091         % call SPGL1
0092         iter_spgl1 = []; err_spgl1 = []; supp_diff_spgl1 = []; time_spgl1 = [];
0093         % if ~exist('spgSetParms','file'); error('Solver SPGL1 is not found.'); end
0094         
0095         W_new = ones(N,1);
0096         delta = norm(e);
0097         delta = max(sqrt(M)*sigma,1e-1);
0098         x_spgl1 = [];
0099         
0100         for rwt_itr = 1:rwt*0
0101             if adp_wt == 0
0102                 spg_opts = spgSetParms('verbosity',0,'weights',W_new,'iterations',maxiter);
0103                 [x_spgl1,r,g,info] = spgl1(AA,y,0,delta,x_spgl1,spg_opts);
0104                 % [x_spgl1,r,g,info] = spgl1(AA,y,1e5,0,x_spgl1,spg_opts);
0105             else
0106                 spg_opts = wspgSetParms('verbosity',0,'omega',0.3);
0107                 [x_spgl1,r,g,info] = wspgl1(AA,y,0,delta,x_spgl1,spg_opts);
0108             end
0109             
0110             rerr = norm(x_spgl1-x)/norm(x);
0111             iter_spgl1 = [iter_spgl1  (info.nProdA+info.nProdAt)/2];
0112             err_spgl1 = [err_spgl1 rerr];
0113             supp_diff_sgpl1 = [supp_diff_spgl1 length(setxor(gamma_orig,find(x_spgl1)))];
0114             time_spgl1 = [time_spgl1 info.timeTotal];
0115             [alpha beta epsilon] = weight_param(rwt_mode,rwt_itr,x_spgl1,M);
0116             W_new = 1/alpha./(beta*abs(x_spgl1)+epsilon);
0117             % W_new = 1./(abs(x_spgl1)+0.1);
0118             r_spgl1 = A(x_spgl1)-y; Atr_spgl1 = AT(r_spgl1);        
0119         end
0120         
0121         %% SpaRSA
0122         tau = max(1e-3*max(abs(Aty)),sigma*sqrt(log(N)));
0123         % tau = max(abs(Atr_spgl1));
0124         
0125         iter_sparsa = [];
0126         err_sparsa = [];
0127         time_sparsa = [];
0128         x_sparsa = 0;
0129         supp_diff_sparsa = [];
0130         W_new = 1;
0131         
0132         % [alpha beta epsilon] = weight_param(5,rwt_itr,x,M);
0133         % W_new = 1/alpha./(beta*abs(x)+epsilon);
0134         % W_new = abs(Atr_spgl1)/tau;
0135         
0136         for rwt_itr = 1:rwt
0137             psi_function = @(x,tau) soft(x,tau*W_new);
0138             phi_function = @(x) sum(abs(W_new.*x));
0139             tic
0140             [x_SpaRSA,x_debias_SpaRSA_m,obj_SpaRSA_m_cont,...
0141                 times_SpaRSA_m_cont,debias_start_SpaRSA_m,mse_SpaRSA_m,taus_SpaRSA_m, numA, numAt]= ...
0142                 SpaRSA_adpW(y,A,tau,...
0143                 'Monotone',0,...
0144                 'AT',AT,...
0145                 'adp_wt',adp_wt,...
0146                 'Debias',0,...
0147                 'Initialization',x_sparsa,...
0148                 'StopCriterion',2,...
0149                 'ToleranceA',1e-4,...
0150                 'MaxiterA',maxiter,...
0151                 'psi',psi_function,...
0152                 'phi',phi_function,...
0153                 'Verbose',0,...
0154                 'W_new',W_new,...
0155                 'True_x',x,...
0156                 'Continuation',1,...
0157                 'Continuationsteps',-1);
0158             x_sparsa = x_SpaRSA;
0159             time_sparsa = [time_sparsa toc];
0160             
0161             iter_sparsa = [iter_sparsa (numA+numAt)/2];
0162             rerr = norm(x-x_sparsa)/norm(x);
0163             err_sparsa = [err_sparsa rerr];
0164             
0165             supp_diff_sparsa = [supp_diff_sparsa length(setxor(gamma_orig,find(x_sparsa)))];
0166             
0167             [alpha beta epsilon] = weight_param(rwt_mode,rwt_itr,x_sparsa,M);
0168             W_new = 1/alpha./(beta*abs(x_sparsa)+epsilon);
0169             
0170             %         gamma = find(abs(x)>tau*100);
0171             %         W_new = ones(N,1);
0172             %         ewt_b = 2*length(y)*(norm(x,2)/norm(x,1))^2;
0173             %         W_new(gamma) = 1./abs(x(gamma))/ewt_b;
0174         end
0175         % r_sparsa = A(x_sparsa)-y; Atr_sparsa = AT(r_sparsa);
0176         
0177         %% NESTA
0178         %     opts = [];
0179         %     opts.Verbose = 0;
0180         %     opts.tolvar = 1e-5;
0181         %     opts.MaxIntIter = 5;
0182         %     delta = max(sigma*sqrt(M),1e-6);
0183         %     muf = 1e-3;
0184         %
0185         %     iter_nesta = [];
0186         %     err_nesta = [];
0187         %     supp_diff_nesta = [];
0188         %     x_nesta = [];
0189         %     tic;
0190         %     % constrained
0191         %     [x_nesta,niter,resid,outData,optsOut] = NESTA_adpW(A,AT,y,muf,delta,opts);
0192         %     % unconstrained
0193         %     % La = norm( A*A' );
0194         %     % [x_nesta,niter,resid,outData,optsOut] = NESTA_UP(A,[],y,tau,La,muf,opts);
0195         %     time_nesta = toc;
0196         %     iter_nesta = [iter_nesta niter*3];
0197         %     rerr = norm(x-x_nesta)/norm(x);
0198         %     err_nesta = [err_nesta rerr];
0199         %     supp_diff_nesta = [supp_diff_nesta length(setxor(gamma_orig,find(x_nesta)))];
0200         
0201         if largescale == 0
0202             %% Adaptive support & weight selection
0203             in = [];
0204             in.A = A; in.y = y; in.x = x;
0205             in.x_init = zeros(N,1); in.max_rwt = 0;
0206             in.tau = tau;
0207             in.rwt_mode = rwt_mode;
0208             in.delx_mode = delx_mode;
0209             in.debias = 0;
0210             in.verbose = 0;
0211             out = script_rwtBPDN_adaptive(in);
0212             xh_adp = out.x_out;
0213             iter_adp = out.iter;
0214             gamma_adp = out.gamma;
0215             W_new = out.W_new;
0216             xh_adp_init = out.x_init;
0217             err_adp = out.err;
0218             time_adp = out.time;
0219             supp_diff_adp = out.supp_diff;
0220             
0221             %% Oracle rwt BPDN
0222             alpha = 5; beta = 10;
0223             W_new = ones(N,1); W_new(gamma_orig) = 1/alpha./(beta*abs(x(gamma_orig)));
0224             % epsilon = 0.1; W_new = 1/alpha./(beta*(abs(x))+epsilon);
0225             
0226             % To check the accuracy of rwtBPDN_adaptive set
0227             % W_new = out.W_new/tau;
0228             
0229             in = [];
0230             in.tau = tau;
0231             in.maxiter = maxiter;
0232             in.x_orig = x;
0233             in.record = 1;
0234             in.delx_mode = 'mil';
0235             %     in.Te = rwt_itr;
0236             AW = A*diag(1./W_new);
0237             tic;
0238             out_new = BPDN_homotopy_function(AW, y, in); %BPDN
0239             % time_orac = out_new.time;
0240             time_orac = toc;
0241             xh_orac = out_new.x_out.*(1./W_new);
0242             gamma_orac = out_new.gamma;
0243             iter_orac = out_new.iter;
0244             
0245             
0246             %% oracle LS
0247             x_LS = zeros(N,1);
0248             x_LS(gamma_orig) = A(:,gamma_orig)\y;
0249         end
0250                 
0251         %%
0252         if largescale
0253             fprintf('Image-%s-%d, wavelet-%s, adp_wt = %d @ rwt = %d ... ',sType,n,wType,adp_wt, rwt);
0254             
0255             I = ifwt2(reshape(x,n,n),g0,g1,J,sym); Ir_max = max(abs(I(:)));
0256             Ir_sparsa = ifwt2(reshape(x_sparsa,n,n),g0,g1,J,sym);
0257             % Ir_sparsa(Ir_sparsa<0) = 0;
0258             % Ir_sparsa(Ir_sparsa>255) = 255;
0259             
0260             mse_sparsa = norm(Ir_sparsa(:)-I(:))^2/N;
0261             psnr_sparsa = 10*log10(Ir_max*Ir_max/mse_sparsa);
0262             ser_sparsa = norm(I(:))^2/norm(Ir_sparsa(:)-I(:))^2;
0263             fprintf('(PSNR,MSE,SER): SpaRSA-%3.4g,%3.4g,%3.4g  ', psnr_sparsa, mse_sparsa, ser_sparsa);
0264             
0265             % Ir_spgl1 = ifwt2(reshape(x_spgl1,n,n),g0,g1,J,sym);
0266             % mse_spgl1 = norm(Ir_spgl1(:)-I(:))^2/N;
0267             % psnr_spgl1 = 10*log10(Ir_max*Ir_max/mse_spgl1);
0268             % ser_spgl1 = norm(I(:))^2/norm(Ir_spgl1(:)-I(:))^2;
0269             % fprintf('  SPGL1-%3.4g,%3.4g,%3.4g.',psnr_spgl1, mse_spgl1, ser_spgl1)
0270             
0271             fprintf('\n');
0272             
0273             Ir_stack = cat(3,Ir_stack,Ir_sparsa);
0274         end
0275         
0276         %% Save results...
0277         
0278         % fprintf('sim iter. %d: \n',sim);
0279         if largescale            
0280             SIM_stack{sim} = [sim, tau, ...
0281                 psnr_sparsa, mse_sparsa, ser_sparsa, sum(iter_sparsa,2), sum(time_sparsa,2)];
0282             %   psnr_spgl1, mse_spgl1, ser_spgl1, sum(iter_spgl1,2), sum(time_spgl1,2)];
0283             % norm(x-x_nesta)/norm(x), sum(iter_nesta,2), sum(time_nesta,2)];
0284         else
0285             SIM_stack{sim} = [sim, tau, (norm(x-x_LS)/norm(x)), (norm(x-xh_orac)/norm(x)), iter_orac, ...
0286                 (norm(x-xh_adp)/norm(x)), sum(iter_adp,2), sum(time_adp,2), ...
0287                 (norm(x-x_sparsa)/norm(x)), sum(iter_sparsa,2), sum(time_sparsa,2)]; 
0288             %     (norm(x-x_spgl1)/norm(x)), sum(iter_spgl1,2), sum(time_spgl1,2)];
0289             % norm(x-x_nesta)/norm(x), sum(iter_nesta,2), sum(time_nesta,2)];
0290         end
0291         index = 0;
0292         if largescale == 0
0293             % adaptive rwt
0294             index = index+1;
0295             SIM_memory{sim}{index,1} = 'adp';
0296             SIM_memory{sim}{index,2} = iter_adp;
0297             SIM_memory{sim}{index,3} = err_adp;
0298             SIM_memory{sim}{index,4} = supp_diff_adp;
0299             SIM_memory{sim}{index,5} = time_adp;
0300         end
0301         
0302         % sparsa
0303         index = index+1;
0304         SIM_memory{sim}{index,1} = 'sparsa';
0305         SIM_memory{sim}{index,2} = iter_sparsa;
0306         SIM_memory{sim}{index,3} = err_sparsa;
0307         SIM_memory{sim}{index,4} = supp_diff_sparsa;
0308         SIM_memory{sim}{index,5} = time_sparsa;
0309         % spgl1
0310         index = index+1;
0311         SIM_memory{sim}{index,1} = 'spgl1';
0312         SIM_memory{sim}{index,2} = iter_spgl1;
0313         SIM_memory{sim}{index,3} = err_spgl1;
0314         SIM_memory{sim}{index,4} = supp_diff_spgl1;
0315         SIM_memory{sim}{index,5} = time_spgl1;
0316         %     % nesta
0317         %     index = index+1;
0318         %     SIM_memory{sim}{index,1} = 'nesta';
0319         %     SIM_memory{sim}{index,2} = iter_nesta;
0320         %     SIM_memory{sim}{index,3} = err_nesta;
0321         %     SIM_memory{sim}{index,4} = supp_diff_nesta;
0322         %     SIM_memory{sim}{index,5} = time_nesta;
0323         %     % salsa
0324         %     index = index+1;
0325         %     SIM_memory{sim}{6,1} = 'salsa';
0326         %     SIM_memory{sim}{6,2} = iter_salsa;
0327         %     SIM_memory{sim}{6,3} = err_salsa;
0328         %     SIM_memory{sim}{6,4} = supp_diff_salsa;
0329         
0330         % mS = mean(cell2mat(SIM_stack),1);
0331         % fprintf('maxsim %d. oracLS=%3.4g. (err,iter,time): oracWT-%3.4g,%3.4g; adp-%3.4g,%3.4g,%3.4g; rwt-%3.4g,%3.4g,%3.4g; yall1-%3.4g,%3.4g,%3.4g; sparsa-%3.4g,%3.4g,%3.4g; omp-%3.4g,%3.4g; fpc-%3.4g,%3.4g; spgl1-%3.4g,%3.4g. \n', maxsim, mS(2:end));
0332         
0333     end
0334     
0335     
0336     %%
0337     mS = mean(cell2mat(SIM_stack),1);
0338     if largescale == 0
0339         % str1 = 'maxsim %d. tau = %3.4g, oracLS=%3.4g. (err,iter,time): oracWT-%3.4g,%3.4g; adp-%3.4g,%3.4g,%3.4g; sparsa-%3.4g,%3.4g,%3.4g; spgl1-%3.4g,%3.4g,%3.4g; nesta-%3.4g,%3.4g,%3.4g;';
0340         str1 = 'maxsim %d. tau = %3.4g, oracLS=%3.4g. (err,iter,time): oracWT-%3.4g,%3.4g; adp-%3.4g,%3.4g,%3.4g; sparsa-%3.4g,%3.4g,%3.4g; spgl1-%3.4g,%3.4g,%3.4g';
0341     else
0342         % str1 = 'maxsim %d. tau = %3.4g, sparsa-%3.4g,%3.4g,%3.4g; spgl1-%3.4g,%3.4g,%3.4g; nesta-%3.4g,%3.4g,%3.4g;';
0343         str1 = 'maxsim %d. tau = %3.4g,\n~sparsa-%3.4g,%3.4g,%3.4g,%3.4g,%3.4g;'; %\n spgl1-%3.4g,%3.4g,%3.4g';
0344     end
0345     str2 = sprintf([str1,' \n'], maxsim, mS(2:end));
0346     fprintf(str2);       
0347     fprintf('(N,M,SNR)-%d,%d,%d, adp-%d,rwt-%d, PSNR: $%2.4g$dB --- ($%d$) \n',N,M,SNR,adp_wt,rwt,10*log10(Ir_max*Ir_max/mS(4)),round(mS(6)));
0348     
0349     Ir_avg = mean(Ir_sparsa,3);
0350     
0351     figure(100); imagesc(I,[0 255]); colormap gray; axis image off;        
0352     figure(101); imagesc(Ir_avg,[0 255]); colormap gray; axis image off;
0353     figure(102); imagesc(Ir_avg-I,[0 255/10]); colormap gray; axis image off;
0354     % figure(103); imagesc(Ir_spgl1,[0 255]); colormap gray; axis image off;
0355     % figure(104); imagesc(Ir_spgl1-I,[0 255/10]); colormap gray; axis image off;
0356     
0357     if SAVE_RESULTS
0358         figure(100); set(gcf, 'Color', 'w');
0359         eval(sprintf('export_fig %s.pdf',sType));
0360         file_name = sprintf('%s_N%d_M%d_SNR%d_adp%d_rwt%d_%s',sType,N,M,SNR,adp_wt,rwt,wType);
0361         figure(101); set(gcf, 'Color', 'w');
0362         eval(sprintf('export_fig %s.pdf',file_name));
0363         % file_name = sprintf('%s_N%d_M%d_SNR%d_adp%d_rwt%d_PSNR%3.4g_ITER%d_TIME%3.4g_spgl1',sType,N,M,SNR,adp_wt,rwt,psnr_spgl1,round(sum(iter_spgl1,2)),sum(time_spgl1,2));
0364         % figure(103); set(gcf, 'Color', 'w');
0365         % eval(sprintf('export_fig %s.pdf',file_name));
0366     end
0367 
0368     fprintf('\n');    
0369         
0370     I_ALL = [I_ALL Ir_avg];
0371     
0372     figure(200);
0373     % subaxis(2,3,exp,'SpacingVert',0.01,'MT',0.05,'MB',0.05, 'ML',0.01,'MR',0.01);
0374     subplot(2,3,exp)
0375     imagesc(Ir_sparsa, [0 255]); axis off image; colormap gray
0376     % subaxis(2,3,exp+3,'SpacingVert',0.01,'MT',0.05,'MB',0.05, 'ML',0.01,'MR',0.01);
0377     subplot(2,3,exp+3)    
0378     imagesc(abs(I-Ir_sparsa), [0 255/10]); axis off image; colormap gray
0379 end
0380 
0381 if SAVE_RESULTS
0382     figure(200); set(gcf, 'Color', 'w');
0383     file_name = sprintf('%s_N%d_M%d_SNR%d_%s_ALL',sType,N,M,SNR,wType);
0384     eval(sprintf('export_fig %s.pdf',file_name));    
0385 end