Description of Simulation_DynamicX_BPDN

0001 % BPDN update with dynamic change in x
0002 % Author: Salman Asif
0003 % Created: November 2008
0004 % Modified: February 2009
0005 
0006 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0007 % In the submitted paper
0008 % 03/02/09
0009 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0010 % DynamicX_..._Blocks_final_new
0011 % DynamicX_..._PcwPoly_final_new
0012 % DynamicX_..._Houste_final_new_Haar
0013 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0014 
0015 % Need WaveLab environment in your MATLAB path to run this simulation.
0016 
0017 clear
0018 close all
0019 % Model of the signal
0020 model_sig = 1; % Test only 1 and 2
0021 
0022 sim_runs = 10; % Use large number of simulations for evaluations
0023 if model_sig == 2
0024     sim_runs = 255; % sim_runs should be the size of image for slices of image
0025 end
0026 
0027 stack_time=zeros(sim_runs,6);
0028 stack_iter=zeros(sim_runs,9);
0029 
0030 for sim_iter = 1:sim_runs
0031 
0032     J = 1;
0033     if model_sig == 1
0034         % Use wavelets on randomly perturbed standard signals from WaveLab
0035         %  Name string: 'Blocks', 'Piece-Polynomial' (Piece-Wise 3rd degree polynomial)
0036 
0037         % signal length
0038         N = 2048;
0039         % number of observations to make
0040         m = 1024;
0041 
0042         Name = 'Blocks';
0043         sig = MakeSignal(Name,N)';
0044         qmf = MakeONFilter('Haar');
0045         
0046         Name = 'Piece-Polynomial';
0047         sig = MakeSignal(Name,N)';
0048         qmf = MakeONFilter('Daubechies',8);
0049 
0050         wc_sig  = FWT_PO(sig,J,qmf);
0051         x = wc_sig;
0052     elseif model_sig == 2
0053         % Use strips from some image
0054         
0055         % I = imread('cameraman.tif');
0056         % I = imread('pirate.tif')';
0057         % I = imread('barb.gif');
0058         % I = (imread('Lena.jpg'));
0059         % load boats;
0060         % I = boats;
0061         % load blocks
0062         % I = blocks;
0063         I = rgb2gray(imread('house.tiff'));
0064         
0065         N = size(I,1);
0066         m = N/2;
0067         sig = double(I(:,sim_iter))/256;
0068         qmf = MakeONFilter('Haar');
0069         % qmf = MakeONFilter('Daubechies',4);
0070         wc_sig  = FWT_PO(sig,J,qmf);
0071         x = wc_sig;
0072     else
0073         disp('Do it yourself');
0074     end
0075 
0076 
0077     % measurement matrix
0078     G = randn(m,N)/sqrt(m);
0079     % A = (orth(G'))';
0080     A = G;
0081 
0082     % observations
0083     sigma = .00;
0084     e = randn(m,1)*sigma;
0085     y = A*x + e;
0086     %figure(101); clf; plot(y); hold on; plot(e,'r');
0087 
0088     % regularization parameter
0089     tau = 0.005*max(abs(A'*y));    
0090     % tau = sigma * sqrt(log(N)*2); % BPDN
0091 
0092     % Solve for measurements y
0093     tic
0094     [xp, gamma_x, iter_xp, t_homotopy1] = BPDN_homotopy_function(A, y, tau,4*N); %BPDN
0095     t_homotopy1 = toc;
0096     tolA_h = tau*sum(abs(xp))+1/2*(norm(A*xp-y))^2;
0097 
0098     % GPSR
0099     pdg_scale = 1e-8;
0100     t_BB_mono = 0;
0101     tolA = tolA_h+pdg_scale;
0102     stopCri=4;
0103     debias = 0;
0104     tic;
0105     [x_BB_mono,x_debias_BB_mono,obj_BB_mono,...
0106         times_BB_mono,debias_start_BB_mono,mse]= ...
0107         GPSR_BB(y,A,tau,...
0108         'Debias',debias,...
0109         'Monotone',1,...
0110         'Initialization',0,...
0111         'Continuation',0,...
0112         'StopCriterion',stopCri,...
0113         'Verbose', 0,...
0114         'MaxiterA', 1000,...
0115         'ToleranceA',tolA,...
0116         'ToleranceD',0.001);
0117     t_BB_mono = toc;
0118 
0119     opts.f_value_tol = tolA;
0120     M = [];
0121     tic
0122     n = N;
0123     tic;
0124     [x_FPC, Out_FPC] = FPC_AS(n,A,y,tau,M,opts); % call FPC_AS
0125     t_FPC = toc;
0126     
0127     % Model the change in signal
0128     if model_sig == 1
0129         % Perturb the signal slightly at random
0130         sig_mod = DynamicX_varying_model(Name,N)';
0131         wc_mod  = FWT_PO(sig_mod,J,qmf);
0132         xt = wc_mod;
0133         
0134     elseif model_sig == 2
0135         % Use strips from some image
0136         pm1 = randsrc;
0137         sig_mod = double(I(:, sim_iter+1))/256;
0138         wc_mod  = FWT_PO(sig_mod,J,qmf);
0139         xt = wc_mod;
0140     else
0141         disp('hehehe...');
0142     end
0143     dx = xt-x;
0144     yt = A*xt;
0145 
0146     pk_old = A'*(A*xp-y);
0147     AtAgx = A(:,gamma_x)'*A(:,gamma_x);
0148     iAtAgx = inv(AtAgx);
0149 
0150     if norm(y-yt)>2*eps
0151         tic;
0152         [xp_h, gamma_xh, iter_xp_update, t_homotopy_update] = DynamicX_BPDN_function(A, AtAgx, iAtAgx, y, yt, xp, gamma_x, pk_old, tau, 4*N);
0153         t_homotopy_update = toc;
0154     else
0155         xp_h = xp;
0156         gamma_xh = gamma_x;
0157         iter_xp_update = 0;
0158         t_homotopy_update = 0;
0159     end
0160     tolA_h = tau*sum(abs(xp_h))+1/2*(norm(A*xp_h-yt))^2;
0161 
0162     % GPSR
0163     tolA = tolA_h+pdg_scale;
0164     stopCri=4;
0165     debias = 0;
0166     tic;
0167     % PUTTING CONTINUATION IN HERE MAKES THINGS WORSE
0168     [x_BB_mono_update,x_debias_BB_mono,obj_BB_mono,...
0169         times_BB_mono,debias_start_BB_mono,mse]= ...
0170         GPSR_BB(yt,A,tau,...
0171         'Debias',debias,...
0172         'Monotone',1,...
0173         'Initialization',x_BB_mono,...
0174         'StopCriterion',stopCri,...
0175         'ToleranceA',tolA,...
0176         'MaxiterA', 1000,...
0177         'Verbose', 0,...
0178         'ToleranceD',0.001);
0179     t_BB_mono_update = toc;
0180 
0181     opts.x0 = x_FPC;
0182     opts.f_value_tol = tolA;
0183     M = [];
0184     tic
0185     [x_FPC_update, Out_FPC_update] = FPC_AS(n,A,yt,tau,M,opts); % call FPC_AS
0186     t_FPC_update = toc;
0187     
0188     stack_time(sim_iter,:) = [t_homotopy1 t_BB_mono t_FPC t_homotopy_update t_BB_mono_update t_FPC_update];
0189     % For old DynamicX_GPSR_Blocks_final or similar
0190     %stack_iter(sim_iter,:) = [iter_xp iter_xp_update norm(dx) norm(xt-xp_h)/norm(xt) norm(xt-x_BB_mono_update)/norm(xt) tau length(gamma_x)];
0191 
0192     % For new mat files
0193     stack_iter(sim_iter,:) = [iter_xp iter_xp_update length(times_BB_mono) (Out_FPC_update.nProdA+Out_FPC_update.nProdAt)/2 norm(xt-xp_h)/norm(xt) norm(xp_h-x_BB_mono_update)/norm(xt) norm(xt-x_FPC_update)/norm(xt) tau length(gamma_x)];
0194     % I_rec(:,sim_iter) = IWT_PO(xp,J,qmf);
0195 end
0196 
0197 % Generate table for performance
0198 disp(['       DynamicX          BPDN homotopy           GPSR            FPC_AS']);
0199 disp(['   nProdAtA    CPU     nProdAtA    CPU    nProdAtA   CPU    nProdAtA   CPU']);
0200 outer_iterations = 1;
0201 table_stats = zeros(1,8);
0202 table_stats(:,3) = mean(reshape(stack_iter(:,1),sim_runs,outer_iterations));
0203 table_stats(:,1) = mean(reshape(stack_iter(:,2),sim_runs,outer_iterations));
0204 table_stats(:,5) = mean(reshape(stack_iter(:,3),sim_runs,outer_iterations));
0205 table_stats(:,7) = mean(reshape(stack_iter(:,4),sim_runs,outer_iterations));
0206 table_stats(:,4) = mean(reshape(stack_time(:,1),sim_runs,outer_iterations));
0207 table_stats(:,2) = mean(reshape(stack_time(:,4),sim_runs,outer_iterations));
0208 table_stats(:,6) = mean(reshape(stack_time(:,5),sim_runs,outer_iterations));
0209 table_stats(:,8) = mean(reshape(stack_time(:,6),sim_runs,outer_iterations));
0210 table_stats
Simulation_DynamicX_BPDN_Wavelet

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
