From 38bfe22200c09ba295b3087e8f0faaad45d568ce Mon Sep 17 00:00:00 2001
From: "Timothy A. Sipkens" <tsipkens@gmail.com>
Date: Fri, 28 Feb 2020 16:34:23 -0800
Subject: [PATCH 1/2] Update .gitignore

---
 .gitignore | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/.gitignore b/.gitignore
index 8f9f3ae..c9b725c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,7 +7,8 @@ results/
 +uncertainty
 *.asv
 
-main_exper_fn18.m
+main_exper_fn18_invert.m
+main_exper_fn18_post.m
 main_exper_ua19_salt.m
 main_exper_gra_test.m
 main_exper_gra.m

From 7f2610b849b92a76073e64a907a4eba886d9791b Mon Sep 17 00:00:00 2001
From: "Timothy A. Sipkens" <tsipkens@gmail.com>
Date: Fri, 28 Feb 2020 17:16:56 -0800
Subject: [PATCH 2/2] Added tikhonov_op2d_bf.m

This is a contribution from @ArashNaseri, with minor nomenclaute updates.
+ Bug fixes in a couple other optimize functions.
---
 +optimize/exp_dist_op.m      |  2 +-
 +optimize/tikhonov_op.m      |  4 +-
 +optimize/tikhonov_op2d_bf.m | 75 ++++++++++++++++++++++++++++++++++++
 3 files changed, 78 insertions(+), 3 deletions(-)
 create mode 100644 +optimize/tikhonov_op2d_bf.m

diff --git a/+optimize/exp_dist_op.m b/+optimize/exp_dist_op.m
index 4c85871..c8ff6e4 100644
--- a/+optimize/exp_dist_op.m
+++ b/+optimize/exp_dist_op.m
@@ -46,7 +46,7 @@
 
 disp('Optimizing exponential distance regularization:');
 tools.textbar(0);
-for ii=length(lambda):-1:1
+for ii=length(lambda):-1:1 % reverse loop to pre-allocate
     %-- Store case parameters ----------------------%
     output(ii).lambda = lambda(ii);
     output(ii).lm = sqrt(Gd(1,1));
diff --git a/+optimize/tikhonov_op.m b/+optimize/tikhonov_op.m
index 412a26c..2f41317 100644
--- a/+optimize/tikhonov_op.m
+++ b/+optimize/tikhonov_op.m
@@ -46,7 +46,7 @@
 
 disp('Optimizing Tikhonov regularization:');
 tools.textbar(0);
-for ii=length(lambda):-1:1
+for ii=length(lambda):-1:1 % reverse loop to pre-allocate
     output(ii).lambda = lambda(ii); % store regularization parameter
     
     %-- Perform inversion --%
@@ -68,7 +68,7 @@
 if ~isempty(x_ex) % if exact solution is supplied
     [~,ind_min] = min([output.chi]);
 else
-    ind_min = [];
+    [~,ind_min] = max([output.B]);
 end
 lambda = output(ind_min).lambda;
 x = output(ind_min).x;
diff --git a/+optimize/tikhonov_op2d_bf.m b/+optimize/tikhonov_op2d_bf.m
new file mode 100644
index 0000000..229610a
--- /dev/null
+++ b/+optimize/tikhonov_op2d_bf.m
@@ -0,0 +1,75 @@
+
+% TIKHONOV_OP2D_BF  Finds optimal lambda and alhpa for Tikhonov solver.
+% Author: Arash Naseri, Timothy Sipkens, 2020-02-28
+%=========================================================================%
+
+function [x,lambda,alpha,out,chi] = tikhonov_op2d_bf(A,b,C,d,span1,span2,order,n,x_ex,xi,solver)
+
+
+%-- Parse inputs ---------------------------------------------------------%
+if ~exist('order','var'); order = []; end
+if ~exist('xi','var'); xi = []; end
+if ~exist('x_ex','var'); x_ex = []; end
+if ~exist('solver','var'); solver = []; end
+%-------------------------------------------------------------------------%
+
+
+%-- Compute credence, fit, and Bayes factor ------------------------------%
+% Initially meshing the domain of (lambda, alpha ) to roughly find the 
+% location of global extremum of B
+lambda = logspace(log10(span1(1)),log10(span1(1)),3);
+alpha =  logspace(log10(span2(1)),log10(span2(1)),3);
+[lambda_mat,alpha_mat] = meshgrid(lambda,alpha);
+param = [lambda_mat(:),alpha_mat(:)]; % set of lambda and alpha to consider
+x_length = size(A,2);
+
+Lpr0 = invert.tikhonov_lpr(order,n,x_length); % get Tikhonov matrix
+
+tools.textbar(0);
+for ii=length(param):-1:1 % reverse loop to pre-allocate
+    out(ii).lambda = param(ii,1); % store regularization parameter
+    out(ii).alpha = param(ii,2); % store regularization parameter
+    
+    %-- Perform inversion --%
+    [out(ii).x,~,Lpr0] = invert.tikhonov(...
+        [param(ii,2).*A;C],[param(ii,2).*b;d],param(ii,1),Lpr0,[],xi,solver);
+    %-- Store ||Ax-b|| and Euclidean error --%
+    if ~isempty(x_ex); out(ii).chi = norm(out(ii).x-x_ex); end
+     out(ii).Axb = norm(A*out(ii).x-b);
+    
+    %-- Compute credence, fit, and Bayes factor --%
+    out(ii).x = invert.tikhonov([param(ii,2)*A;C],[param(ii,2)*b;d],...
+        param(ii,1),Lpr0,[],xi,solver);
+    [out(ii).B,out(ii).F,out(ii).C] = ...
+        optimize.bayesf_b([param(ii,2)*A;C],[param(ii,2)*b;d],...
+        out(ii).x,Lpr0,param(ii,1));
+    tools.textbar((length(param)-ii+1)/length(param));
+end
+
+%-- Record a rough estimate of the solution --%
+[~,ind_min] = max([output.B]); % get optimal w.r.t. Bayes factor
+out(1).ind_min = ind_min;
+%-------------------------------------------------------------------------%
+
+
+%-- Add fminsearch step to optimize parameter set --------------------%
+disp('Optimizing Tikhonov regularization:');
+fun = @(lambda) log(-1*optimize.bayesf_b([lambda(2)*A;C],[lambda(2)*b;d],invert.tikhonov...
+    ([lambda(2)*A;C],[lambda(2)*b;d],lambda(1),Lpr0,[],xi,solver),Lpr0,lambda(1)));
+
+y0 = [out(ind_min).lambda out(ind_min).alpha]; % initial guess for fminsearch
+options = optimset('TolFun',10^-8,'TolX',10^-8,'Display','iter');
+y1 = fminsearch(fun,y0,options); % get optimal lambda and alpha
+
+lambda = y1(1); % assign output variables
+alpha = y1(2);
+disp('Complete.');
+%---------------------------------------------------------------------%
+
+
+x = invert.tikhonov(...
+    [alpha*A;C],[alpha*b;d],lambda,Lpr0,[],xi,solver);
+chi = norm(x-x_ex);
+
+
+end