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