diff --git a/+invert/optimize_exp_dist.m b/+optimize/exp_dist.m similarity index 100% rename from +invert/optimize_exp_dist.m rename to +optimize/exp_dist.m diff --git a/+invert/optimize_mart.m b/+optimize/mart.m similarity index 100% rename from +invert/optimize_mart.m rename to +optimize/mart.m diff --git a/+invert/optimize_tikhonov.m b/+optimize/tikhonov.m similarity index 100% rename from +invert/optimize_tikhonov.m rename to +optimize/tikhonov.m diff --git a/+invert/optimize_twomark.m b/+optimize/twomark.m similarity index 100% rename from +invert/optimize_twomark.m rename to +optimize/twomark.m diff --git a/+invert/optimize_twomey.m b/+optimize/twomey.m similarity index 100% rename from +invert/optimize_twomey.m rename to +optimize/twomey.m diff --git a/.gitignore b/.gitignore index 10fac27..6ad8adc 100644 --- a/.gitignore +++ b/.gitignore @@ -9,7 +9,6 @@ main_exper_ua19_salt.m main_jas19_noise.m submodule_update.sh +uncertainty -run_inversions_d.m run_inversions_e.m run_inversions_f.m run_inversions_g.m diff --git a/README.md b/README.md index fd11a8c..854d034 100644 --- a/README.md +++ b/README.md @@ -104,14 +104,20 @@ Contains various functions used to invert the measured data for the desired two-dimensional distribution. This includes implementations of least-squares, Tikhonov regularization, Twomey, Twomey-Markowski (including using the method of [Buckley et al. (2017)][3]), and the multiplicative algebraic reconstruction -technique (MART). Also included are functions that, given the true distribution, -can determine the optimal number of iterations or the optimal regularization -parameter. Development is underyway on the use of an exponential covariance +technique (MART). Development is underyway on the use of an exponential covariance function to correlate pixel values and reduce reconstruction errors. Details on these approaches to inversion are provided in the associated paper, [Sipkens et al. (Submitted)][1]. +#### +optimize + +This package mirrors the content of the +inver package but, +given the true distribution, aims to determine the optimal number of +iterations for the Twomey and MART schemes or the optimal regularization +parameter for the Twomey-Markowski and Tikhonov methods. + + #### +tfer_PMA This is imported from a package distributed with [Sipkens et al. (2019)][2]. diff --git a/run_inversions_a.m b/run_inversions_a.m index 7bc7254..c34e0cd 100644 --- a/run_inversions_a.m +++ b/run_inversions_a.m @@ -29,7 +29,7 @@ %% Tikhonov (0th) implementation disp('Performing Tikhonov (0th) regularization...'); tic; -[x_tk0,lambda_tk0,out_tk0] = invert.optimize_tikhonov(Lb*A,Lb*b,n_x(1),... +[x_tk0,lambda_tk0,out_tk0] = optimize.tikhonov(Lb*A,Lb*b,n_x(1),... [1e-2,1e2],x0,0,[],'interior-point'); t.tk0 = toc; disp('Inversion complete.'); @@ -41,7 +41,7 @@ %% Tikhonov (1st) implementation disp('Performing Tikhonov (1st) regularization...'); tic; -[x_tk1,lambda_tk1,out_tk1] = invert.optimize_tikhonov(Lb*A,Lb*b,n_x(1),... +[x_tk1,lambda_tk1,out_tk1] = optimize.tikhonov(Lb*A,Lb*b,n_x(1),... [1e-6,1e2],x0,1,[],'interior-point'); t.tk1 = toc; disp('Inversion complete.'); @@ -53,7 +53,7 @@ %% Tikhonov (2nd) implementation disp('Performing Tikhonov (2nd) regularization...'); tic; -[x_tk2,lambda_tk2,out_tk2] = invert.optimize_tikhonov(Lb*A,Lb*b,n_x(1),... +[x_tk2,lambda_tk2,out_tk2] = optimize.tikhonov(Lb*A,Lb*b,n_x(1),... [1e-2,1e2],x0,2,[],'interior-point'); t.tk2 = toc; disp('Inversion complete.'); @@ -67,7 +67,7 @@ disp('Performing MART...'); tic; [x_mart,iter_mart,out_mart] = ... - invert.optimize_mart(A,b,x_init,1:300,x0); + optimize.mart(A,b,x_init,1:300,x0); t.MART = toc; disp('Inversion complete.'); disp(' '); @@ -79,7 +79,7 @@ disp('Performing Twomey...'); tic; [x_two,iter_two,out_two] = ... - invert.optimize_twomey(A,b,x_init,1:500,x0); + optimize.twomey(A,b,x_init,1:500,x0); t.two = toc; disp('Completed Twomey.'); @@ -92,7 +92,7 @@ disp('Performing Twomey-Markowski-Buckley...'); tic; [x_two_mh,Sf_two_mh,out_two_mh] = ... - invert.optimize_twomark(A,b,Lb,n_x(1),... + optimize.twomark(A,b,Lb,n_x(1),... x_init,35,[1e2,1e-5],x0,'Buckley'); t.two_mh = toc; diff --git a/run_inversions_d.m b/run_inversions_d.m new file mode 100644 index 0000000..5dc5f43 --- /dev/null +++ b/run_inversions_d.m @@ -0,0 +1,103 @@ + +% RUN_INVERSIONS_D Invert multuple times to determine CPU time. +% Author: Timothy Sipkens, 2019-07-22 +%=========================================================================% + +%% Initial guess for iterative schemes +b_init = b; +b_init(b_init<(1e-5*max(b_init))) = 0; +x_init = interp2(grid_b.edges{2}',grid_b.edges{1}',... + reshape(full(b_init)./(A*ones(size(x0))),grid_b.ne),... + grid_x.elements(:,2),grid_x.elements(:,1)); +x_init(isnan(x_init)) = 0; +x_init(isinf(x_init)) = 0; +x_init = sparse(max(0,x_init)); +chi.init = norm(x0-x_init); +x_init_m = grid_x.marginalize(x_init); + + +for ii=1:20 + +%% Least squares +disp('Performing LS inversion...'); +tic; +x_LSQ = invert.lsq(A,b,'interior-point'); +t.LSQ(ii) = toc; +disp('Inversion complete.'); +disp(' '); + +chi.LSQ = norm(x0-x_LSQ); + + +%% Tikhonov (0th) implementation +disp('Performing Tikhonov (0th) regularization...'); +tic; +x_tk0 = invert.tikhonov(Lb*A,Lb*b,n_x(1),lambda_tk0,0); +t.tk0(ii) = toc; +disp('Inversion complete.'); +disp(' '); + +chi.tk0(ii) = norm(x0-x_tk0); + + +%% Tikhonov (1st) implementation +disp('Performing Tikhonov (1st) regularization...'); +tic; +x_tk1 = invert.tikhonov(Lb*A,Lb*b,n_x(1),lambda_tk1,1); +t.tk1(ii) = toc; +disp('Inversion complete.'); +disp(' '); + +chi.tk1(ii) = norm(x0-x_tk1); + + +%% Tikhonov (2nd) implementation +disp('Performing Tikhonov (2nd) regularization...'); +% lambda_tk2 = 8e1; +tic; +x_tk2 = invert.tikhonov(Lb*A,Lb*b,n_x(1),lambda_tk2,2); +t.tk2(ii) = toc; +disp('Inversion complete.'); +disp(' '); + +chi.tk2(ii) = norm(x0-x_tk2); + + +%% MART, Maximum entropy regularized solution + +disp('Performing MART...'); +tic; +x_mart = invert.mart(A,b,x_init,299); +t.mart(ii) = toc; +disp('Inversion complete.'); +disp(' '); + +chi.mart = norm(x0-x_mart); + + +%% Twomey +disp('Performing Twomey...'); +tic; +x_two = invert.twomey(A,b,x_init,500); +t.two(ii) = toc; +disp('Completed Twomey.'); +disp(' '); + +chi.two = norm(x0-x_two); + + +%% Twomey-Markowski-Buckley +disp('Performing Twomey-Markowski...'); +tic; +x_two_mh = invert.twomark(A,b,Lb,n_x(1),... + x_init,35,'Buckley',1/Sf_two_mh); +t.two_mh(ii) = toc; +disp('Completed Twomey-Markowski.'); + +chi.two_mh = norm(x0-x_two_mh); + + +end + + +