From cac34eb90c7966a30975453ed2f00b18fe37d54f Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Wed, 16 Feb 2022 12:05:22 -0700
Subject: [PATCH] Also add the Matlab version.

---
 MCMC-Laplace/Matlab/README.md         |   9 ++
 MCMC-Laplace/Matlab/exact_values.m    | 173 ++++++++++++++++++++++++++
 MCMC-Laplace/Matlab/forward_solver.m  |  48 +++++++
 MCMC-Laplace/Matlab/get_statistics.m  |  46 +++++++
 MCMC-Laplace/Matlab/log_probability.m |  24 ++++
 MCMC-Laplace/Matlab/main.m            |  74 +++++++++++
 MCMC-Laplace/Matlab/plot_solution.m   |  52 ++++++++
 MCMC-Laplace/Matlab/precomputations.m |  67 ++++++++++
 8 files changed, 493 insertions(+)
 create mode 100644 MCMC-Laplace/Matlab/README.md
 create mode 100644 MCMC-Laplace/Matlab/exact_values.m
 create mode 100644 MCMC-Laplace/Matlab/forward_solver.m
 create mode 100644 MCMC-Laplace/Matlab/get_statistics.m
 create mode 100644 MCMC-Laplace/Matlab/log_probability.m
 create mode 100644 MCMC-Laplace/Matlab/main.m
 create mode 100644 MCMC-Laplace/Matlab/plot_solution.m
 create mode 100644 MCMC-Laplace/Matlab/precomputations.m

diff --git a/MCMC-Laplace/Matlab/README.md b/MCMC-Laplace/Matlab/README.md
new file mode 100644
index 0000000..ae3dc37
--- /dev/null
+++ b/MCMC-Laplace/Matlab/README.md
@@ -0,0 +1,9 @@
+# Bayesian inversion benchmark: Matlab version
+ 
+This code implements a Matlab version of the benchmark in D. Aristoff and W. Bangerth, "A benchmark for the Bayesian inversion of coefficients in partial differential equations," arXiv:2102.07263. The article is available on the web at https://arxiv.org/abs/2102.07263.
+
+The benchmark is intended for comparing advanced methods for Bayesian inverse problems. The benchmark is designed to be more complex than simple testcases like Gaussian mixtures, but not so complex as to be intractable or unrepeatable. The benchmark concerns the determination of a spatially variable coefficient, discretized by 64 values, in a Poisson equation, based on point measurements of the solution. 
+
+This code runs a basic Metropolis-Hastings sampler for the posterior distribution. To run it, put all the .m and the .mat files from this repository in your current Matlab directory, and type "main" into the Matlab command window. This starts the main program, main.m, which lets you can choose a number of independent Markov chains, the length of each chain, the lag time between recording samples, and the number of parallel workers to use. Other parameters for the simulations can be found within precomputations.m.
+
+The posterior distribution has argument theta, which is a vector of length 64 representing the spatially variable coefficient in the Poisson equation. This vector represents the values of the coefficient in a 8x8 grid. In this code, theta is represented by an 8x8 matrix. To convert a vector, theta, of length 64 -- with entries associated to coefficients as described in the article above -- into an equivalent 8x8 matrix suitable for this code, use the Matlab command "theta = reshape(theta,[8 8])';". The matrix of exact coefficients can be obtained using the command "theta = ones(8,8); theta(2:3,2:3) = 0.1; theta(6:7,6:7) = 10;".
diff --git a/MCMC-Laplace/Matlab/exact_values.m b/MCMC-Laplace/Matlab/exact_values.m
new file mode 100644
index 0000000..d66e47f
--- /dev/null
+++ b/MCMC-Laplace/Matlab/exact_values.m
@@ -0,0 +1,173 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%% list of "exact" measurement values, z_hat %%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+z_hat = [0.06076511762259369;
+         0.09601910120848481;
+         0.1238852517838584; 
+         0.1495184117375201;
+         0.1841596127549784;  
+         0.2174525028261122;
+         0.2250996160898698;  
+         0.2197954769002993;
+         0.2074695698370926;  
+         0.1889996477663016;
+         0.1632722532153726;  
+         0.1276782480038186;
+         0.07711845915789312; 
+         0.09601910120848552;
+         0.2000589533367983;  
+         0.3385592591951766;
+         0.3934300024647806;  
+         0.4040223892461541;
+         0.4122329537843092;  
+         0.4100480091545554;
+         0.3949151637189968;  
+         0.3697873264791232;
+         0.33401826235924;    
+         0.2850397806663382;
+         0.2184260032478671;  
+         0.1271121156350957;
+         0.1238852517838611;  
+         0.3385592591951819;
+         0.7119285162766475;  
+         0.8175712861756428;
+         0.6836254116578105;  
+         0.5779452419831157;
+         0.5555615956136897;  
+         0.5285181561736719;
+         0.491439702849224;   
+         0.4409367494853282;
+         0.3730060082060772;  
+         0.2821694983395214;
+         0.1610176733857739;  
+         0.1495184117375257;
+         0.3934300024647929;  
+         0.8175712861756562;
+         0.9439154625527653;  
+         0.8015904115095128;
+         0.6859683749254024;  
+         0.6561235366960599;
+         0.6213197201867315;   
+         0.5753611315000049;
+         0.5140091754526823;  
+         0.4325325506354165;
+         0.3248315148915482;  
+         0.1834600412730086;
+         0.1841596127549917;  
+         0.4040223892461832;
+         0.6836254116578439;  
+         0.8015904115095396;
+         0.7870119561144977;  
+         0.7373108331395808;
+         0.7116558878070463;  
+         0.6745179049094283;
+         0.6235300574156917;  
+         0.5559332704045935;
+         0.4670304994474178;  
+         0.3499809143811;
+         0.19688263746294;    
+         0.2174525028261253;
+         0.4122329537843404;  
+         0.5779452419831566;
+         0.6859683749254372;  
+         0.7373108331396063;
+         0.7458811983178246;  
+         0.7278968022406559;
+         0.6904793535357751;  
+         0.6369176452710288;
+         0.5677443693743215;  
+         0.4784738764865867;
+         0.3602190632823262;  
+         0.2031792054737325;
+         0.2250996160898818;  
+         0.4100480091545787;
+         0.5555615956137137;  
+         0.6561235366960938;
+         0.7116558878070715;  
+         0.727896802240657;
+         0.7121928678670187;  
+         0.6712187391428729;
+         0.6139157775591492;  
+         0.5478251665295381;
+         0.4677122687599031;  
+         0.3587654911000848;
+         0.2050734291675918;  
+         0.2197954769003094;
+         0.3949151637190157;  
+         0.5285181561736911;
+         0.6213197201867471;  
+         0.6745179049094407;
+         0.690479353535786;   
+         0.6712187391428787;
+         0.6178408289359514;  
+         0.5453605027237883;
+         0.489575966490909;   
+         0.4341716881061278;
+         0.3534389974779456;  
+         0.2083227496961347;
+         0.207469569837099;   
+         0.3697873264791366;
+         0.4914397028492412;  
+         0.5753611315000203;
+         0.6235300574157017;  
+         0.6369176452710497;
+         0.6139157775591579;  
+         0.5453605027237935;
+         0.4336604929612851;  
+         0.4109641743019312;
+         0.3881864790111245;  
+         0.3642640090182592;
+         0.2179599909280145;  
+         0.1889996477663011;
+         0.3340182623592461;  
+         0.4409367494853381;
+         0.5140091754526943;  
+         0.5559332704045969;
+         0.5677443693743304;  
+         0.5478251665295453;
+         0.4895759664908982;  
+         0.4109641743019171;
+         0.395727260284338;   
+         0.3778949322004734;
+         0.3596268271857124;  
+         0.2191250268948948;
+         0.1632722532153683;  
+         0.2850397806663325;
+         0.373006008206081;   
+         0.4325325506354207;
+         0.4670304994474315;  
+         0.4784738764866023;
+         0.4677122687599041;  
+         0.4341716881061055;
+         0.388186479011099;   
+         0.3778949322004602;
+         0.3633362567187364;  
+         0.3464457261905399;
+         0.2096362321365655;  
+         0.1276782480038148;
+         0.2184260032478634;  
+         0.2821694983395252;
+         0.3248315148915535;  
+         0.3499809143811097;
+         0.3602190632823333;  
+         0.3587654911000799;
+         0.3534389974779268;  
+         0.3642640090182283;
+         0.35962682718569;    
+         0.3464457261905295;
+         0.3260728953424643;  
+         0.180670595355394;
+         0.07711845915789244; 
+         0.1271121156350963;
+         0.1610176733857757;  
+         0.1834600412730144;
+         0.1968826374629443;  
+         0.2031792054737354;
+         0.2050734291675885;  
+         0.2083227496961245;
+         0.2179599909279998;  
+         0.2191250268948822;
+         0.2096362321365551;  
+         0.1806705953553887;
+         0.1067965550010013];
diff --git a/MCMC-Laplace/Matlab/forward_solver.m b/MCMC-Laplace/Matlab/forward_solver.m
new file mode 100644
index 0000000..833532b
--- /dev/null
+++ b/MCMC-Laplace/Matlab/forward_solver.m
@@ -0,0 +1,48 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%% forward solver function %%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%INPUTS: 
+%theta = current 8x8 parameter matrix
+%lbl = cell labeling function
+%A_loc = matrix of local contributions to A
+%Id = Identity matrix of size 128x128
+%boundaries = labels of boundary cells
+%b = right hand side of linear system (AU = b)
+%M = measurement matrix
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%OUTPUTS:
+%z = vector of measurements
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function z = forward_solver(theta,lbl,A_loc,Id,boundaries,b,M)
+
+%initialize matrix A for FEM linear solve, AU = b
+A = zeros(33^2,33^2);
+
+%loop over cells to build A
+for i=0:31
+    for j=0:31    %build A by summing over contribution from each cell
+
+        %find local coefficient in 8x8 grid
+        thet = theta(floor(i/4)+1,floor(j/4)+1);
+
+        %update A by including contribution from cell (i,j)
+        dof = [lbl(i,j),lbl(i,j+1),lbl(i+1,j+1),lbl(i+1,j)];
+        A(dof,dof) = A(dof,dof) + thet*A_loc;
+    end
+end
+
+%enforce boundary condition
+A(boundaries,:) = Id(boundaries,:);
+A(:,boundaries) = Id(:,boundaries);
+
+%sparsify A
+A = sparse(A);
+
+%solve linear equation for coefficients, U
+U = A\b;
+
+%get new z values
+z = reshape((M*U)',[13^2 1]);
+
+        
diff --git a/MCMC-Laplace/Matlab/get_statistics.m b/MCMC-Laplace/Matlab/get_statistics.m
new file mode 100644
index 0000000..ef3c48a
--- /dev/null
+++ b/MCMC-Laplace/Matlab/get_statistics.m
@@ -0,0 +1,46 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%% compute statistics on data set %%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%INPUTS: 
+%data = tensor of theta samples from each lag time and chain
+%theta_means = means of theta over each independent chain
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%OUTPUTS:
+%theta_mean = overall mean of chains
+%covars = covariance matrices of each independent chain
+%autocov = mean of autocovariance matrix over all the chains
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [theta_mean,covars,autocovar] = get_statistics(data,theta_means);
+
+%compute overall mean of data, and get size of data matrix
+theta_mean = mean(theta_means,3);
+[~,~,L,N] = size(data);
+
+%initialize covariance matrices and mean autocovariance matrix
+covars = zeros(64,64,N);
+autocovar = zeros(64,64,2*L-1);
+
+%compute covariance matrices and mean autocovariance matrix
+for n=1:N                            %loop over independent Markov chains
+    
+    %get data from chain n
+    data_ = reshape(permute(data(:,:,:,n),[3 2 1]),[L 64]);
+
+    %compute autocovariance matrix of chain n
+    mat = xcov(data_,'unbiased');
+
+    %store covariance matrix of chain n
+    covars(:,:,n) = reshape(mat(L,:),64,64);
+
+    %update mean autocovariance matrix
+    autocovar = autocovar + reshape(mat',[64 64 2*L-1]);
+end
+
+%compute mean of autocovariance matrix
+autocovar = autocovar(1:64,1:64,L:2*L-1)/N;
+
+
+
+
+
diff --git a/MCMC-Laplace/Matlab/log_probability.m b/MCMC-Laplace/Matlab/log_probability.m
new file mode 100644
index 0000000..62ed698
--- /dev/null
+++ b/MCMC-Laplace/Matlab/log_probability.m
@@ -0,0 +1,24 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%% compute log probability, log pi %%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%INPUTS: 
+%theta = current 8x8 parameter matrix
+%z = current vector of measurements
+%z_hat = vector of "exact" measurements
+%sig = standard deviation parameter in likelihood
+%sig_pr = standard deviation parameter in prior
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%OUTPUTS:
+%log_pi = logarithm of posterior probability
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function log_pi = log_probability(theta,z,z_hat,sig,sig_pr)
+
+%compute log likelihood
+log_L = -sum((z-z_hat).^2)/(2*sig^2);
+
+%compute log prior
+log_pi_pr = -sum(log(theta).^2,'all')/(2*sig_pr^2);
+
+%compute log posterior
+log_pi = log_L+log_pi_pr;
diff --git a/MCMC-Laplace/Matlab/main.m b/MCMC-Laplace/Matlab/main.m
new file mode 100644
index 0000000..f4c4786
--- /dev/null
+++ b/MCMC-Laplace/Matlab/main.m
@@ -0,0 +1,74 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%% Run MCMC sampler to estimate posterior distribution %%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%define number of chains, chain length, and lag time
+N = input('number of independent Markov chains: ');
+N_L = input('length of each Markov chain: ');
+lag = input('lag time for measurements: ');
+workers = input('number of parallel workers: ');
+L = N_L/lag;
+
+%open Matlab parallel pool
+parpool(workers)
+
+%load precomputations
+load precomputations.mat
+
+%define lag time and data matrix
+data = zeros(8,8,L,N);   %data matrix of samples at lag times
+theta_means = zeros(8,8,N);      %overall mean of theta
+
+tic
+
+parfor n=1:N
+
+    %set initial theta, theta mean, and z values of chain
+    theta = theta0;
+    theta_mean = 0;
+    z = z0;
+
+    for m=1:L
+
+        for l=1:lag
+
+            %define proposal, theta_tilde
+            xi = normrnd(0,sig_prop,[8 8]);
+            theta_tilde = theta.*exp(xi);
+        
+            %compute new z values
+            z_tilde = forward_solver_(theta_tilde);
+        
+            %compute posterior log probability of theta_tilde
+            log_pi_tilde = log_probability_(theta_tilde,z_tilde);
+            log_pi = log_probability_(theta,z);
+        
+            %compute acceptance probability; accept proposal appropriately
+            accept = exp(log_pi_tilde-log_pi) ... 
+                    *prod(theta_tilde./theta,'all');
+            if rand<accept
+                theta = theta_tilde;   %accept new theta values
+                z = z_tilde;           %record associated measurements
+            end
+            
+            %update mean of theta
+            theta_mean = theta_mean + theta;
+        
+        end
+    
+        %update data matrix
+        data(:,:,m,n) = theta;
+
+    end
+
+    %update theta means
+    theta_means(:,:,n) = theta_mean/N_L;
+end
+
+toc
+
+%compute statistics on data set
+[theta_mean,covars,autocovars] = get_statistics(data,theta_means);
+
+%save data to Matlab workspace, labeled by N and N_L
+save (['data_N_' num2str(N) '_N_L_ ' num2str(N_L) '.mat'])
diff --git a/MCMC-Laplace/Matlab/plot_solution.m b/MCMC-Laplace/Matlab/plot_solution.m
new file mode 100644
index 0000000..fb353f4
--- /dev/null
+++ b/MCMC-Laplace/Matlab/plot_solution.m
@@ -0,0 +1,52 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%% plots solution, u, to Poisson equation %%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%% associated to theta and a 32x32 mesh %%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%input matrix theta for plotting (e.g., theta = theta_hat)
+theta = input('choose theta for plotting u: theta = ');
+
+%construct mass matrix, M_plot, for plotting
+xsp = 0:0.02:1;
+n = length(xsp);
+Mp = zeros(n,n,33^2);
+for k=1:33^2
+    c = inv_lbl(k);
+    for i=1:n
+        for j=1:n
+            Mp(i,j,k) = phi(xsp(i)-h*c(1),xsp(j)-h*c(2));
+        end
+    end
+end
+Mp = reshape(Mp,[n^2 33^2]);
+
+%run forward solver on mean of theta
+A = zeros(33^2,33^2);
+for i=0:31
+    for j=0:31    %build A by summing over contribution from each cell
+
+        %find local coefficient in 8x8 grid
+        thet = theta(floor(i/4)+1,floor(j/4)+1);
+
+        %update A by including contribution from cell (i,j)
+        dof = [lbl(i,j),lbl(i,j+1),lbl(i+1,j+1),lbl(i+1,j)];
+        A(dof,dof) = A(dof,dof) + thet*A_loc;
+    end
+end
+
+%enforce boundary condition
+A(boundaries,:) = Id(boundaries,:);
+
+%sparsify A
+A = sparse(A);
+
+%solve linear equation for coefficients, U
+U = A\b;
+
+%close all current plots
+close all
+
+%plot solution
+figure
+zs = reshape(Mp*U,[n n]);
+surf(xsp,xsp,zs)
\ No newline at end of file
diff --git a/MCMC-Laplace/Matlab/precomputations.m b/MCMC-Laplace/Matlab/precomputations.m
new file mode 100644
index 0000000..76cec1b
--- /dev/null
+++ b/MCMC-Laplace/Matlab/precomputations.m
@@ -0,0 +1,67 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%% do all precomputations necessary for MCMC simulations %%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%define mesh width
+h = 1/32;   
+
+%define characteristic function of unit square
+S = @(x,y) heaviside(x).*heaviside(y) ...
+           .*(1-heaviside(x-h)).*(1-heaviside(y-h));
+
+%define tent function on the domain [-h,h]x[-h,h]
+phi = @(x,y) ((x+h).*(y+h).*S(x+h,y+h) + (h-x).*(h-y).*S(x,y) ... 
+          + (x+h).*(h-y).*S(x+h,y) + (h-x).*(y+h).*S(x,y+h))/h^2;
+
+%define function that converts from (i,j) to dof, and its inverse
+lbl = @(i,j) 33*j+i+1;
+inv_lbl = @(k) [k-1-33*floor((k-1)/33),floor((k-1)/33)];
+
+%construct measurement matrix, M
+xs = 1/14:1/14:13/14;    %measurement points
+M = zeros(13,13,33^2);
+for k=1:33^2
+    c = inv_lbl(k);
+    for i=1:13
+        for j=1:13
+            M(i,j,k) = phi(xs(i)-h*c(1),xs(j)-h*c(2));
+        end
+    end
+end
+M = reshape(M,[13^2 33^2]);
+
+%construct exact coefficient matrix, thetas0
+theta_hat = ones(8,8);
+theta_hat(2:3,2:3) = 0.1;
+theta_hat(6:7,6:7) = 10;
+
+%construct local overlap matrix, A_loc, and identity matrix Id
+A_loc = [2/3  -1/6  -1/3  -1/6;
+          -1/6  2/3  -1/6  -1/3;
+          -1/3 -1/6   2/3  -1/6;
+          -1/6 -1/3  -1/6   2/3];
+Id = eye(33^2);
+
+%locate boundary labels
+boundaries = [lbl(0:1:32,0),lbl(0:1:32,32),lbl(0,1:1:31),lbl(32,1:1:31)];
+
+%define RHS of FEM linear system, AU = b
+b = ones(33^2,1)*10*h^2;
+b(boundaries) = zeros(128,1);    %enforce boundary conditions on b
+
+%load exact z_hat values
+exact_values
+
+%set global parameters and functions for simulation
+sig = 0.05;             %likelihood standard deviation
+sig_pr = 2;             %prior (log) standard deviation
+sig_prop = 0.0725;      %proposal (log) standard deviation
+theta0 = ones(8,8);     %initial theta values
+forward_solver_ = @(theta) ... 
+                  forward_solver(theta,lbl,A_loc,Id,boundaries,b,M);
+log_probability_ = @(theta,z) log_probability(theta,z,z_hat,sig,sig_pr);
+
+%find initial z values
+z0 = forward_solver_(theta0);
+
+save precomputations.mat
-- 
2.39.5