nikcleju@54: % File: study_analysis_rec_algos nikcleju@54: % Run global experiment to compare algorithms used for analysis-based reconstruction nikcleju@54: % and plot phast transition graphs nikcleju@54: nikcleju@54: clear all nikcleju@54: close all nikcleju@54: nikcleju@54: % ================================= nikcleju@54: % Set up experiment parameters nikcleju@54: %================================== nikcleju@54: % Which form factor, delta and rho we want nikcleju@54: sigmas = 1.2; nikcleju@54: deltas = 0.05:0.05:0.95; nikcleju@54: rhos = 0.05:0.05:0.95; nikcleju@54: % deltas = [0.95]; nikcleju@54: % rhos = [0.1]; nikcleju@54: %deltas = 0.3:0.3:0.9; nikcleju@54: %rhos = 0.3:0.3:0.9; nikcleju@54: nikcleju@54: % Number of vectors to generate each time nikcleju@54: numvects = 100; nikcleju@54: nikcleju@54: % Add noise nikcleju@54: % This is norm(signal)/norm(noise), so power, not energy nikcleju@54: SNRdb = 20; % virtually no noise nikcleju@54: nikcleju@54: % Show progressbar ? (not recommended when running on parallel threads) nikcleju@54: do_progressbar = 0; nikcleju@54: nikcleju@54: % Value of lambda nikcleju@54: lambda = 1e-2; nikcleju@54: nikcleju@54: % What algos to run nikcleju@54: do_abs_ompk = 1; nikcleju@54: do_abs_ompeps = 1; nikcleju@54: do_abs_tst = 1; nikcleju@54: do_abs_bp = 1; nikcleju@54: do_gap = 1; nikcleju@54: do_nesta = 0; nikcleju@54: nikcleju@54: % ================================= nikcleju@54: % Processing the parameters nikcleju@54: %================================== nikcleju@54: % Set up parameter structure nikcleju@54: count = 0; nikcleju@54: for isigma = 1:sigmas nikcleju@54: for idelta = 1:numel(deltas) nikcleju@54: for irho = 1:numel(rhos) nikcleju@54: sigma = sigmas(isigma); nikcleju@54: delta = deltas(idelta); nikcleju@54: rho = rhos(irho); nikcleju@54: nikcleju@54: d = 200; nikcleju@54: p = round(sigma*d); nikcleju@54: m = round(delta*d); nikcleju@54: l = round(d - rho*m); nikcleju@54: nikcleju@54: params(count+1).d = d; nikcleju@54: params(count+1).p = p; nikcleju@54: params(count+1).m = m; nikcleju@54: params(count+1).l = l; nikcleju@54: nikcleju@54: count = count + 1; nikcleju@54: end nikcleju@54: end nikcleju@54: end nikcleju@54: nikcleju@54: % Compute noiselevel from db nikcleju@54: noiselevel = 1 / (10^(SNRdb/10)); nikcleju@54: nikcleju@54: %load study_analysis_init nikcleju@54: nikcleju@54: % Generate an analysis operator Omega nikcleju@54: Omega = Generate_Analysis_Operator(d, p); nikcleju@54: nikcleju@54: % Progressbar nikcleju@54: if do_progressbar nikcleju@54: progressbar('Total', 'Current parameters'); nikcleju@54: end nikcleju@54: nikcleju@54: % Init times nikcleju@54: elapsed_abs_ompk = 0; nikcleju@54: elapsed_abs_ompeps = 0; nikcleju@54: elapsed_abs_tst = 0; nikcleju@54: elapsed_abs_bp = 0; nikcleju@54: elapsed_gap = 0; nikcleju@54: elapsed_nesta = 0; nikcleju@54: nikcleju@54: % Init results structure nikcleju@54: results = []; nikcleju@54: nikcleju@54: % Prepare progressbar reduction variables nikcleju@54: % if do_progressbar nikcleju@54: % incr2 = 1/numvects; nikcleju@54: % incr1 = 1/numvects/count; nikcleju@54: % frac2 = 0; nikcleju@54: % frac1 = 0; nikcleju@54: % end nikcleju@54: nikcleju@54: % ======== nikcleju@54: % Run nikcleju@54: % ======== nikcleju@54: parfor iparam = 1:numel(params) nikcleju@54: nikcleju@54: % Read current parameters nikcleju@54: d = params(iparam).d; nikcleju@54: p = params(iparam).p; nikcleju@54: m = params(iparam).m; nikcleju@54: l = params(iparam).l; nikcleju@54: nikcleju@54: % Init stuff nikcleju@54: xrec_abs_ompk = zeros(d, numvects); nikcleju@54: xrec_abs_ompeps = zeros(d, numvects); nikcleju@54: xrec_abs_tst = zeros(d, numvects); nikcleju@54: xrec_abs_bp = zeros(d, numvects); nikcleju@54: xrec_gap = zeros(d, numvects); nikcleju@54: xrec_nesta = zeros(d, numvects); nikcleju@54: % nikcleju@54: err_abs_ompk = zeros(numvects,1); nikcleju@54: relerr_abs_ompk = zeros(numvects,1); nikcleju@54: err_abs_ompeps = zeros(numvects,1); nikcleju@54: relerr_abs_ompeps = zeros(numvects,1); nikcleju@54: err_abs_tst = zeros(numvects,1); nikcleju@54: relerr_abs_tst = zeros(numvects,1); nikcleju@54: err_abs_bp = zeros(numvects,1); nikcleju@54: relerr_abs_bp = zeros(numvects,1); nikcleju@54: err_gap = zeros(numvects,1); nikcleju@54: relerr_gap = zeros(numvects,1); nikcleju@54: err_nesta = zeros(numvects,1); nikcleju@54: relerr_nesta = zeros(numvects,1); nikcleju@54: nikcleju@54: % Generate data based on parameters nikcleju@54: [x0,y,M,Lambda] = Generate_Data_Known_Omega(Omega, d,p,m,l,noiselevel, numvects,'l0'); nikcleju@54: nikcleju@54: % For every generated signal do nikcleju@54: for iy = 1:size(y,2) nikcleju@54: nikcleju@54: % Compute epsilon nikcleju@54: epsilon = noiselevel * norm(y(:,iy)); nikcleju@54: nikcleju@54: %-------------------------------- nikcleju@54: % Reconstruct (and measure delay) nikcleju@54: % Compute reconstruction error nikcleju@54: %-------------------------------- nikcleju@54: % ABS-OMPk nikcleju@54: if (do_abs_ompk) nikcleju@54: timer_abs_ompk = tic; nikcleju@54: xrec_abs_ompk(:,iy) = ABS_OMPk_approx(y(:,iy), Omega, M, p-l, lambda); nikcleju@54: elapsed_abs_ompk = elapsed_abs_ompk + toc(timer_abs_ompk); nikcleju@54: % nikcleju@54: err_abs_ompk(iy) = norm(x0(:,iy) - xrec_abs_ompk(:,iy)); nikcleju@54: relerr_abs_ompk(iy) = err_abs_ompk(iy) / norm(x0(:,iy)); nikcleju@54: end nikcleju@54: % ABS-OMPeps nikcleju@54: if (do_abs_ompeps) nikcleju@54: timer_abs_ompeps = tic; nikcleju@54: xrec_abs_ompeps(:,iy) = ABS_OMPeps_approx(y(:,iy), Omega, M, epsilon, lambda); nikcleju@54: elapsed_abs_ompeps = elapsed_abs_ompeps + toc(timer_abs_ompeps); nikcleju@54: % nikcleju@54: err_abs_ompeps(iy) = norm(x0(:,iy) - xrec_abs_ompeps(:,iy)); nikcleju@54: relerr_abs_ompeps(iy) = err_abs_ompeps(iy) / norm(x0(:,iy)); nikcleju@54: end nikcleju@54: % ABS-TST nikcleju@54: if (do_abs_tst) nikcleju@54: timer_abs_tst = tic; nikcleju@54: xrec_abs_tst(:,iy) = ABS_TST_approx(y(:,iy), Omega, M, epsilon, lambda); nikcleju@54: elapsed_abs_tst = elapsed_abs_tst + toc(timer_abs_tst); nikcleju@54: % nikcleju@54: err_abs_tst(iy) = norm(x0(:,iy) - xrec_abs_tst(:,iy)); nikcleju@54: relerr_abs_tst(iy) = err_abs_tst(iy) / norm(x0(:,iy)); nikcleju@54: end nikcleju@54: % ABS-BP nikcleju@54: if (do_abs_bp) nikcleju@54: timer_abs_bp = tic; nikcleju@54: xrec_abs_bp(:,iy) = ABS_BP_approx(y(:,iy), Omega, M, epsilon, lambda); nikcleju@54: elapsed_abs_bp = elapsed_abs_bp + toc(timer_abs_bp); nikcleju@54: % nikcleju@54: err_abs_bp(iy) = norm(x0(:,iy) - xrec_abs_bp(:,iy)); nikcleju@54: relerr_abs_bp(iy) = err_abs_bp(iy) / norm(x0(:,iy)); nikcleju@54: end nikcleju@54: % GAP nikcleju@54: if (do_gap) nikcleju@54: gapparams = []; nikcleju@54: gapparams.num_iteration = 40; nikcleju@54: gapparams.greedy_level = 0.9; nikcleju@54: gapparams.stopping_coefficient_size = 1e-4; nikcleju@54: gapparams.l2solver = 'pseudoinverse'; nikcleju@54: gapparams.noise_level = noiselevel; nikcleju@54: timer_gap = tic; nikcleju@54: xrec_gap(:,iy) = GAP(y(:,iy), M, M', Omega, Omega', gapparams, zeros(d,1)); nikcleju@54: elapsed_gap = elapsed_gap + toc(timer_gap); nikcleju@54: % nikcleju@54: err_gap(iy) = norm(x0(:,iy) - xrec_gap(:,iy)); nikcleju@54: relerr_gap(iy) = err_gap(iy) / norm(x0(:,iy)); nikcleju@54: end nikcleju@54: % NESTA nikcleju@54: if (do_nesta) nikcleju@54: try nikcleju@54: timer_nesta = tic; nikcleju@54: xrec_nesta(:,iy) = do_nesta_DemoNonProjector(x0(:,iy), M, Omega', 0); nikcleju@54: elapsed_nesta = elapsed_nesta + toc(timer_nesta); nikcleju@54: catch err nikcleju@54: disp('*****ERROR: NESTA throwed error *****'); nikcleju@54: xrec_nesta(:,iy) = zeros(size(x0(:,iy))); nikcleju@54: end nikcleju@54: % nikcleju@54: err_nesta(iy) = norm(x0(:,iy) - xrec_nesta(:,iy)); nikcleju@54: relerr_nesta(iy) = err_nesta(iy) / norm(x0(:,iy)); nikcleju@54: end nikcleju@54: nikcleju@54: % Update progressbar nikcleju@54: % if do_progressbar nikcleju@54: % %frac2 = iy/numvects; nikcleju@54: % %frac1 = ((iparam-1) + frac2) / count; nikcleju@54: % if norm(frac2 - 1) < 1e-6 nikcleju@54: % frac2 = 0; nikcleju@54: % end nikcleju@54: % frac2 = frac2 + incr2; nikcleju@54: % frac1 = frac1 + incr1; nikcleju@54: % progressbar(frac1, frac2); nikcleju@54: % end nikcleju@54: end nikcleju@54: nikcleju@54: %-------------------------------- nikcleju@54: % Save results in big stucture & display nikcleju@54: %-------------------------------- nikcleju@54: % Save reconstructed signals nikcleju@54: % Save rel & abs errors nikcleju@54: % Display error nikcleju@54: disp(['Simulation no. ' num2str(iparam)]); nikcleju@54: if (do_abs_ompk) nikcleju@54: results(iparam).xrec_abs_ompk = xrec_abs_ompk; nikcleju@54: results(iparam).err_abs_ompk = err_abs_ompk; nikcleju@54: results(iparam).relerr_abs_ompk = relerr_abs_ompk; nikcleju@54: disp([' ABS_OMPk: avg relative error = ' num2str(mean(relerr_abs_ompk))]); nikcleju@54: end nikcleju@54: if (do_abs_ompeps) nikcleju@54: results(iparam).xrec_abs_ompeps = xrec_abs_ompeps; nikcleju@54: results(iparam).err_abs_ompeps = err_abs_ompeps; nikcleju@54: results(iparam).relerr_abs_ompeps = relerr_abs_ompeps; nikcleju@54: disp([' ABS_OMPeps: avg relative error = ' num2str(mean(relerr_abs_ompeps))]); nikcleju@54: end nikcleju@54: if (do_abs_tst) nikcleju@54: results(iparam).xrec_abs_tst = xrec_abs_tst; nikcleju@54: results(iparam).err_abs_tst = err_abs_tst; nikcleju@54: results(iparam).relerr_abs_tst = relerr_abs_tst; nikcleju@54: disp([' ABS_TST: avg relative error = ' num2str(mean(relerr_abs_tst))]); nikcleju@54: end nikcleju@54: if (do_abs_bp) nikcleju@54: results(iparam).xrec_abs_bp = xrec_abs_bp; nikcleju@54: results(iparam).err_abs_bp = err_abs_bp; nikcleju@54: results(iparam).relerr_abs_bp = relerr_abs_bp; nikcleju@54: disp([' ABS_BP: avg relative error = ' num2str(mean(relerr_abs_bp))]); nikcleju@54: end nikcleju@54: if (do_gap) nikcleju@54: results(iparam).xrec_gap = xrec_gap; nikcleju@54: results(iparam).err_gap = err_gap; nikcleju@54: results(iparam).relerr_gap = relerr_gap; nikcleju@54: disp([' GAP: avg relative error = ' num2str(mean(relerr_gap))]); nikcleju@54: end nikcleju@54: if (do_nesta) nikcleju@54: results(iparam).xrec_nesta = xrec_nesta; nikcleju@54: results(iparam).err_nesta = err_nesta; nikcleju@54: results(iparam).relerr_nesta = relerr_nesta; nikcleju@54: disp([' NESTA: avg relative error = ' num2str(mean(relerr_nesta))]); nikcleju@54: end nikcleju@54: end nikcleju@54: nikcleju@54: % ================================= nikcleju@54: % Save nikcleju@54: % ================================= nikcleju@54: save mat/avgerr_SNR20_lbd1e-2 nikcleju@54: nikcleju@54: % ================================= nikcleju@54: % Plot phase transition nikcleju@54: % ================================= nikcleju@54: %-------------------------------- nikcleju@54: % Prepare nikcleju@54: %-------------------------------- nikcleju@54: %d = 200; nikcleju@54: %m = 190; nikcleju@54: %exactthr = d/m * noiselevel; nikcleju@54: %sigma = 1.2; nikcleju@54: iparam = 1; nikcleju@54: for idelta = 1:numel(deltas) nikcleju@54: for irho = 1:numel(rhos) nikcleju@54: % Create exact recovery count matrix nikcleju@54: % nexact_abs_bp (irho, idelta) = sum(results(iparam).relerr_abs_bp < exactthr); nikcleju@54: % nexact_abs_ompk (irho, idelta) = sum(results(iparam).relerr_abs_ompk < exactthr); nikcleju@54: % nexact_abs_ompeps (irho, idelta) = sum(results(iparam).relerr_abs_ompeps < exactthr); nikcleju@54: % nexact_gap (irho, idelta) = sum(results(iparam).relerr_gap < exactthr); nikcleju@54: % nexact_abs_tst (irho, idelta) = sum(results(iparam).relerr_abs_tst < exactthr); nikcleju@54: % % nexact_nesta(irho, idelta) = sum(results(iparam).relerr_nesta < exactthr); nikcleju@54: nikcleju@54: % Get histogram (for a single param set only!) nikcleju@54: % hist_abs_ompk = hist(results(iparam).relerr_abs_ompk, 0.001:0.001:0.1); nikcleju@54: % hist_abs_ompeps = hist(results(iparam).relerr_abs_ompeps, 0.001:0.001:0.1); nikcleju@54: % hist_abs_tst = hist(results(iparam).relerr_abs_tst, 0.001:0.001:0.1); nikcleju@54: % hist_abs_bp = hist(results(iparam).relerr_abs_bp, 0.001:0.001:0.1); nikcleju@54: % hist_gap = hist(results(iparam).relerr_gap, 0.001:0.001:0.1); nikcleju@54: nikcleju@54: % Compute average error value nikcleju@54: if (do_abs_ompk) nikcleju@54: avgerr_abs_ompk(irho, idelta) = 1 - mean(results(iparam).relerr_abs_ompk); nikcleju@54: avgerr_abs_ompk(avgerr_abs_ompk < 0) = 0; nikcleju@54: end nikcleju@54: if (do_abs_ompeps) nikcleju@54: avgerr_abs_ompeps(irho, idelta) = 1 - mean(results(iparam).relerr_abs_ompeps); nikcleju@54: avgerr_abs_ompeps(avgerr_abs_ompeps < 0) = 0; nikcleju@54: end nikcleju@54: if (do_abs_tst) nikcleju@54: avgerr_abs_tst(irho, idelta) = 1 - mean(results(iparam).relerr_abs_tst); nikcleju@54: avgerr_abs_tst(avgerr_abs_tst < 0) = 0; nikcleju@54: end nikcleju@54: if (do_abs_bp) nikcleju@54: avgerr_abs_bp(irho, idelta) = 1 - mean(results(iparam).relerr_abs_bp); nikcleju@54: avgerr_abs_bp(avgerr_abs_bp < 0) = 0; nikcleju@54: end nikcleju@54: if (do_gap) nikcleju@54: avgerr_gap(irho, idelta) = 1 - mean(results(iparam).relerr_gap); nikcleju@54: avgerr_gap(avgerr_gap < 0) = 0; nikcleju@54: end nikcleju@54: if (do_nesta) nikcleju@54: avgerr_nesta(irho, idelta) = 1 - mean(results(iparam).relerr_nesta); nikcleju@54: avgerr_nesta(avgerr_nesta < 0) = 0; nikcleju@54: end nikcleju@54: nikcleju@54: iparam = iparam + 1; nikcleju@54: end nikcleju@54: end nikcleju@54: nikcleju@54: %-------------------------------- nikcleju@54: % Plot nikcleju@54: %-------------------------------- nikcleju@54: show_phasetrans = @show_phasetrans_win; nikcleju@54: iptsetpref('ImshowAxesVisible', 'on'); nikcleju@54: close all nikcleju@54: figbase = 'figs/avgerr_SNR20_lbd1e-2_'; nikcleju@54: do_save = 1; nikcleju@54: % nikcleju@54: if (do_abs_ompk) nikcleju@54: figure; nikcleju@54: %h = show_phasetrans(nexact_abs_ompk, numvects); nikcleju@54: %bar(0.001:0.001:0.1, hist_abs_ompk); nikcleju@54: h = show_phasetrans(avgerr_abs_ompk, 1); nikcleju@54: if do_save nikcleju@54: figname = [figbase 'ABS_OMPk']; nikcleju@54: saveas(h, [figname '.fig']); nikcleju@54: saveas(h, [figname '.png']); nikcleju@54: saveTightFigure(h, [figname '.pdf']); nikcleju@54: end nikcleju@54: end nikcleju@54: % nikcleju@54: if (do_abs_ompeps) nikcleju@54: figure; nikcleju@54: %h = show_phasetrans(nexact_abs_ompeps, numvects); nikcleju@54: %bar(0.001:0.001:0.1, hist_abs_ompeps); nikcleju@54: h = show_phasetrans(avgerr_abs_ompeps, 1); nikcleju@54: if do_save nikcleju@54: figname = [figbase 'ABS_OMPeps']; nikcleju@54: saveas(h, [figname '.fig']); nikcleju@54: saveas(h, [figname '.png']); nikcleju@54: saveTightFigure(h, [figname '.pdf']); nikcleju@54: end nikcleju@54: end nikcleju@54: % nikcleju@54: if (do_abs_tst) nikcleju@54: figure; nikcleju@54: %h = show_phasetrans(nexact_abs_tst, numvects); nikcleju@54: %bar(0.001:0.001:0.1, hist_abs_tst); nikcleju@54: h = show_phasetrans(avgerr_abs_tst, 1); nikcleju@54: if do_save nikcleju@54: figname = [figbase 'ABS_TST']; nikcleju@54: saveas(h, [figname '.fig']); nikcleju@54: saveas(h, [figname '.png']); nikcleju@54: saveTightFigure(h, [figname '.pdf']); nikcleju@54: end nikcleju@54: end nikcleju@54: % nikcleju@54: if (do_abs_bp) nikcleju@54: figure; nikcleju@54: %h = show_phasetrans(nexact_abs_bp, numvects); nikcleju@54: %bar(0.001:0.001:0.1, hist_abs_bp); nikcleju@54: h = show_phasetrans(avgerr_abs_bp, 1); nikcleju@54: if do_save nikcleju@54: figname = [figbase 'ABS_BP']; nikcleju@54: saveas(h, [figname '.fig']); nikcleju@54: saveas(h, [figname '.png']); nikcleju@54: saveTightFigure(h, [figname '.pdf']); nikcleju@54: end nikcleju@54: end nikcleju@54: % nikcleju@54: if (do_gap) nikcleju@54: figure; nikcleju@54: %h = show_phasetrans(nexact_gap, numvects); nikcleju@54: %bar(0.001:0.001:0.1, hist_gap); nikcleju@54: h = show_phasetrans(avgerr_gap, 1); nikcleju@54: if do_save nikcleju@54: figname = [figbase 'GAP']; nikcleju@54: saveas(h, [figname '.fig']); nikcleju@54: saveas(h, [figname '.png']); nikcleju@54: saveTightFigure(h, [figname '.pdf']); nikcleju@54: end nikcleju@54: end nikcleju@54: % nikcleju@54: if (do_nesta) nikcleju@54: figure; nikcleju@54: %h = show_phasetrans(nexact_nesta, numvects); nikcleju@54: %bar(0.001:0.001:0.1, hist_nesta); nikcleju@54: h = show_phasetrans(avgerr_nesta, 1); nikcleju@54: if do_save nikcleju@54: figname = [figbase 'NESTA']; nikcleju@54: saveas(h, [figname '.fig']); nikcleju@54: saveas(h, [figname '.png']); nikcleju@54: saveTightFigure(h, [figname '.pdf']); nikcleju@54: end nikcleju@54: end