Created opimtize package

This involved transferinf and renaming functions from the +invert package. The README was updated accordingly.
`run_inversion_d.m` was also added back to the program.

.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

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].

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;
 

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
+
+
+