MCMC-Laplace.h

Go to the documentation of this file.

 5,
    /* fe_degree = */ 1,
    "");
  LogLikelihood::Gaussian        log_likelihood(exact_solution, 0.05);
  LogPrior::LogGaussian          log_prior(0, 2);
 
  const unsigned int n_theta = 64;
  for (unsigned int test=0; test<10; ++test)
    {
      std::cout << "Generating output for test " << test << std::endl;
 
      // For each test, read the input file...
      std::ifstream test_input ("../testing/input." + std::to_string(test) + ".txt");
      Assert (test_input, ExcIO());
 
      Vector<double> theta(n_theta);
      for (unsigned int i=0; i<n_theta; ++i)
        test_input >> theta[i];
 
      // ...run it through the forward simulator to get the
      // z vector (which we then output to a file)...
      const Vector<double> z = laplace_problem.evaluate(theta);
      std::ofstream test_output_z ("output." + std::to_string(test) + ".z.txt");
      z.print(test_output_z, 16);
 
      // ...and then also evaluate prior and likelihood for these
      // theta vectors:
      std::ofstream test_output_likelihood ("output." + std::to_string(test) + ".loglikelihood.txt");
      test_output_likelihood.precision(12);
      test_output_likelihood << log_likelihood.log_likelihood(z) << std::endl;
 
      std::ofstream test_output_prior ("output." + std::to_string(test) + ".logprior.txt");
      test_output_prior.precision(12);
      test_output_prior << log_prior.log_prior(theta) << std::endl;
    }
}
@endcode
This code reads in each of the input files (assuming that the executable is located in a
build directory parallel to the `testing/` directory) and outputs its results into the
current directory. The inputs you get from a modified program should then be compared
against the ones stored in the `testing/` directory. They should match to several digits.
 
 
 
An alternative implementation in Matlab
---------------------------------------
 
To facilitate experiments, this directory also contains alternative
implementations of the benchmark. The first one was written by David
Aristoff in Matlab and can be found in the `Matlab/` directory. As is
common in Matlab programs, each function is provided in its own
file. We have verified that the program generates the same results (to
12 or more digits) as the C++ program, using the tests mentioned in
the previous section.
 
 
 
An alternative implementation in Python
---------------------------------------
 
Another alternative, written by Wolfgang Bangerth, can be found in the
`Python/` directory and uses Python3. As is common for Python
programs, the whole program is provided in one file. As for the Matlab
version, we have verified that the program generates the same results
(to 12 or more digits) as the C++ program, using the tests mentioned
in the previous section. In fact, if you just execute the program
as-is, it runs diagnostics that output the errors.
 
This Python version is essentially a literal translation of the Matlab
code. It is not nearly as efficient (taking around 5 times as much
time for each function evaluation) and could probably be optimized
substantially if desired. A good starting point is the insertion of
the local elements into the global matrix in the line
@code{.sh}
            A[np.ix_(dof,dof)] += theta_loc * A_loc
@endcode
that takes up a substantial fraction of the overall run time.
 
 
<a name="ann-Matlab/exact_values.m"></a>
<h1>Annotated version of Matlab/exact_values.m</h1>
@code{.m}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% 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];
@endcode
 
 
<a name="ann-Matlab/forward_solver.m"></a>
<h1>Annotated version of Matlab/forward_solver.m</h1>
@code{.m}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%% 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
        theta_loc = 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) + theta_loc*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 = M*U;
@endcode
 
 
<a name="ann-Matlab/get_statistics.m"></a>
<h1>Annotated version of Matlab/get_statistics.m</h1>
@code{.m}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% 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
%autocovar = 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;
@endcode
 
 
<a name="ann-Matlab/log_probability.m"></a>
<h1>Annotated version of Matlab/log_probability.m</h1>
@code{.m}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% 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;
@endcode
 
 
<a name="ann-Matlab/main.m"></a>
<h1>Annotated version of Matlab/main.m</h1>
@code{.m}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% 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,autocovar] = 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'])
@endcode
 
 
<a name="ann-Matlab/plot_solution.m"></a>
<h1>Annotated version of Matlab/plot_solution.m</h1>
@code{.m}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% 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
        theta_loc = 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) + theta_loc*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)
@endcode
 
 
<a name="ann-Matlab/precomputations.m"></a>
<h1>Annotated version of Matlab/precomputations.m</h1>
@code{.m}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% 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, theta_hat
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
@endcode
 
 
<a name="ann-Python/forward_solver.py"></a>
<h1>Annotated version of Python/forward_solver.py</h1>
@code{.py}
import numpy as np
import scipy.sparse
from scipy.sparse.linalg import spsolve
import time
 
###########################################################################
############ list of "exact" measurement values, z_hat ####################
###########################################################################
 
z_hat = np.array(
     [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])
 
 
###########################################################################
####### do all precomputations necessary for MCMC simulations #############
###########################################################################
 
# Define the mesh width
h = 1/32
 
# Define characteristic function of unit square
def heaviside(x) :
    if x<0 :
        return 0
    else :
        return 1
    
def S(x,y) :
    return heaviside(x)*heaviside(y) * (1-heaviside(x-h))*(1-heaviside(y-h));
 
# Define tent function on the domain [0,2h]x[0,2h]
def phi(x,y) :
    return ((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 conversion function for dof's from 2D to scalar label, and
# its inverse
def ij_to_dof_index(i,j) :
    return 33*j+i
 
def inv_ij_to_dof_index(k) :
    return [k-33*int(k/33),int(k/33)]
 
 
# Construct measurement matrix, M, for measurements
xs = np.arange(1./14,13./14,1./14);    #measurement points
 
M = np.zeros((13,13,33**2));
for k in range(33**2) :
    c = inv_ij_to_dof_index(k)
    for i in range(13) :
        for j in range(13) :
            M[i,j,k] = phi(xs[i]-h*c[0], xs[j]-h*c[1])
M = M.reshape((13**2, 33**2))
M = scipy.sparse.csr_matrix(M);
 
# Construct local overlap matrix, A_loc, and identity matrix Id
A_loc = np.array([[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 = np.eye(33**2,33**2)
 
# Locate boundary labels
boundaries = ([ij_to_dof_index(i,0) for i in range(33)] +
              [ij_to_dof_index(i,32) for i in range(33)] +
              [ij_to_dof_index(0,j+1) for j in range(31)] +
              [ij_to_dof_index(32,j+1) for j in range(31)])
 
# Define RHS of FEM linear system, AU = b
b = np.ones(33**2)*10*h**2
b[boundaries] = 0    #enforce boundary conditions on b
 
 
 
 
 
###########################################################################
###################### forward solver function ############################
###########################################################################
 
def forward_solver(theta) :
    # Initialize matrix A for FEM linear solve, AU = b
    A = np.zeros((33**2,33**2))
 
    # Build A by summing over contribution from each cell
    for i in range(32) :
        for j in range (32) :
            # Find local coefficient in 8x8 grid
            theta_loc = theta[int(i/4)+int(j/4)*8]
 
            # Update A by including contribution from cell (i,j)
            dof = [ij_to_dof_index(i,j),
                   ij_to_dof_index(i,j+1),
                   ij_to_dof_index(i+1,j+1),
                   ij_to_dof_index(i+1,j)]
            A[np.ix_(dof,dof)] += theta_loc * A_loc
 
    # Enforce boundary condition: Zero out rows and columns, then
    # put a one back into the diagonal entries.
    A[boundaries,:] = 0
    A[:,boundaries] = 0
    A[boundaries,boundaries] = 1
 
    # Solve linear equation for coefficients, U, and then
    # get the Z vector by multiplying by the measurement matrix
    u = spsolve(scipy.sparse.csr_matrix(A), b)
    
    z = M * u
    
    return z
 
 
 
 
###########################################################################
################# compute log probability, log pi #########################
###########################################################################
 
def log_likelihood(theta) :
    z = forward_solver(theta)
    misfit = z - z_hat
    sig = 0.05             #likelihood standard deviation
    return -np.dot(misfit,misfit)/(2*sig**2)
 
def log_prior(theta) :
    sig_pr = 2             #prior (log) standard deviation
    return -np.linalg.norm(np.log(theta))**2/(2*sig_pr**2)
 
def log_posterior(theta) :
    return log_likelihood(theta) + log_prior(theta)
 
 
 
###########################################################################
############# A function to test against known output #####################
###########################################################################
 
 
def verify_against_stored_tests() :
    for i in range(10) :
        print ("Verifying against data set", i)
 
        # Read the input vector
        f_input = open ("../testing/input.{}.txt".format(i), 'r')
        theta = np.fromfile(f_input, count=64, sep=" ")
 
        # Then compute both the forward solution and its statistics.
        # This is not efficiently written here (it calls the forward
        # solver twice), but we don't care about efficiency here since
        # we are only computing with ten samples
        this_z              = forward_solver(theta)
        this_log_likelihood = log_likelihood(theta)
        this_log_prior      = log_prior(theta)
 
        # Then also read the reference output generated by the C++ program:
        f_output_z = open ("../testing/output.{}.z.txt".format(i), 'r')
        f_output_likelihood = open ("../testing/output.{}.loglikelihood.txt".format(i), 'r')
        f_output_prior = open ("../testing/output.{}.logprior.txt".format(i), 'r')
 
        reference_z              = np.fromfile(f_output_z, count=13**2, sep=" ")
        reference_log_likelihood = float(f_output_likelihood.read())
        reference_log_prior = float(f_output_prior.read())
 
        print ("  || z-z_ref ||  : ",
               np.linalg.norm(this_z - reference_z))
        print ("  log likelihood : ",
               "Python value=", this_log_likelihood,
               "(C++ reference value=", reference_log_likelihood,
               ", error=", abs(this_log_likelihood - reference_log_likelihood),
               ")")
        print ("  log prior      : ",
               "Python value=", this_log_prior,
               "(C++ reference value=", reference_log_prior,
               ", error=", abs(this_log_prior - reference_log_prior),
               ")")
 
 
def time_forward_solver() :
    begin = time.time()
 
    n_runs = 100
    for i in range(n_runs) :
        # Create a random vector (with entries between 0 and 1), scale
        # it by a factor of 4, subtract 2, then take the exponential
        # of each entry to get random entries between e^{-2} and
        # e^{+2}
        theta = np.exp(np.random.rand(64) * 4 - 2)
        z = forward_solver(theta)
    end = time.time()
    print ("Time per forward evaluation:", (end-begin)/n_runs)
 
verify_against_stored_tests()
time_forward_solver()
@endcode
 
 
<a name="ann-mcmc-laplace.cc"></a>
<h1>Annotated version of mcmc-laplace.cc</h1>
 * 
 * 
 * 
 * 
 * @code
 *   /* ---------------------------------------------------------------------
 *    *
 *    * Copyright (C) 2019 by the deal.II authors and Wolfgang Bangerth.
 *    *
 *    * This file is part of the deal.II library.
 *    *
 *    * The deal.II library is free software; you can use it, redistribute
 *    * it, and/or modify it under the terms of the GNU Lesser General
 *    * Public License as published by the Free Software Foundation; either
 *    * version 2.1 of the License, or (at your option) any later version.
 *    * The full text of the license can be found in the file LICENSE.md at
 *    * the top level directory of deal.II.
 *    *
 *    * ---------------------------------------------------------------------
 *   
 *    *
 *    * Author: Wolfgang Bangerth, Colorado State University, 2019.
 *    */
 *   
 *   
 *   
 *   #include <deal.II/base/timer.h>
 *   #include <deal.II/base/multithread_info.h>
 *   #include <deal.II/grid/tria.h>
 *   #include <deal.II/dofs/dof_handler.h>
 *   #include <deal.II/grid/grid_generator.h>
 *   #include <deal.II/grid/tria_accessor.h>
 *   #include <deal.II/grid/tria_iterator.h>
 *   #include <deal.II/dofs/dof_accessor.h>
 *   #include <deal.II/fe/fe_q.h>
 *   #include <deal.II/dofs/dof_tools.h>
 *   #include <deal.II/fe/fe_values.h>
 *   #include <deal.II/base/quadrature_lib.h>
 *   #include <deal.II/base/function.h>
 *   #include <deal.II/numerics/vector_tools.h>
 *   #include <deal.II/numerics/matrix_tools.h>
 *   #include <deal.II/lac/vector.h>
 *   #include <deal.II/lac/full_matrix.h>
 *   #include <deal.II/lac/sparse_matrix.h>
 *   #include <deal.II/lac/dynamic_sparsity_pattern.h>
 *   #include <deal.II/lac/solver_cg.h>
 *   #include <deal.II/lac/precondition.h>
 *   #include <deal.II/lac/sparse_ilu.h>
 *   
 *   #include <deal.II/numerics/data_out.h>
 *   
 *   #include <fstream>
 *   #include <iostream>
 *   #include <random>
 *   
 *   #include <deal.II/base/logstream.h>
 *   
 *   using namespace dealii;
 *   
 *   
 * @endcode
 * 
 * The following is a namespace in which we define the solver of the PDE.
 * The main class implements an abstract `Interface` class declared at
 * the top, which provides for an `evaluate()` function that, given
 * a coefficient vector, solves the PDE discussed in the Readme file
 * and then evaluates the solution at the 169 mentioned points.
 * 
 
 * 
 * The solver follows the basic layout of @ref step_4 "step-4", though it precomputes
 * a number of things in the `setup_system()` function, such as the
 * evaluation of the matrix that corresponds to the point evaluations,
 * as well as the local contributions to matrix and right hand side.
 * 
 
 * 
 * Rather than commenting on everything in detail, in the following
 * we will only document those things that are not already clear from
 * @ref step_4 "step-4" and a small number of other tutorial programs.
 * 
 * @code
 *   namespace ForwardSimulator
 *   {
 *     class Interface
 *     {
 *     public:
 *       virtual Vector<double> evaluate(const Vector<double> &coefficients) = 0;
 *   
 *       virtual ~Interface() = default;
 *     };
 *   
 *   
 *   
 *     template <int dim>
 *     class PoissonSolver : public Interface
 *     {
 *     public:
 *       PoissonSolver(const unsigned int global_refinements,
 *                     const unsigned int fe_degree,
 *                     const std::string &dataset_name);
 *       virtual Vector<double>
 *       evaluate(const Vector<double> &coefficients) override;
 *   
 *     private:
 *       void make_grid(const unsigned int global_refinements);
 *       void setup_system();
 *       void assemble_system(const Vector<double> &coefficients);
 *       void solve();
 *       void output_results(const Vector<double> &coefficients) const;
 *   
 *       Triangulation<dim>        triangulation;
 *       FE_Q<dim>                 fe;
 *       DoFHandler<dim>           dof_handler;
 *   
 *       FullMatrix<double>        cell_matrix;
 *       Vector<double>            cell_rhs;
 *       std::map<types::global_dof_index,double> boundary_values;
 *   
 *       SparsityPattern           sparsity_pattern;
 *       SparseMatrix<double>      system_matrix;
 *   
 *       Vector<double>            solution;
 *       Vector<double>            system_rhs;
 *   
 *       std::vector<Point<dim>>   measurement_points;
 *   
 *       SparsityPattern           measurement_sparsity;
 *       SparseMatrix<double>      measurement_matrix;
 *   
 *       TimerOutput  timer;
 *       unsigned int nth_evaluation;
 *   
 *       const std::string &dataset_name;
 *     };
 *   
 *   
 *   
 *     template <int dim>
 *     PoissonSolver<dim>::PoissonSolver(const unsigned int global_refinements,
 *                                       const unsigned int fe_degree,
 *                                       const std::string &dataset_name)
 *       : fe(fe_degree)
 *       , dof_handler(triangulation)
 *       , timer(std::cout, TimerOutput::summary, TimerOutput::cpu_times)
 *       , nth_evaluation(0)
 *       , dataset_name(dataset_name)
 *     {
 *       make_grid(global_refinements);
 *       setup_system();
 *     }
 *   
 *   
 *   
 *     template <int dim>
 *     void PoissonSolver<dim>::make_grid(const unsigned int global_refinements)
 *     {
 *       Assert(global_refinements >= 3,
 *              ExcMessage("This program makes the assumption that the mesh for the "
 *                         "solution of the PDE is at least as fine as the one used "
 *                         "in the definition of the coefficient."));
 *       GridGenerator::hyper_cube(triangulation, 0, 1);
 *       triangulation.refine_global(global_refinements);
 *   
 *       std::cout << "   Number of active cells: " << triangulation.n_active_cells()
 *                 << std::endl;
 *     }
 *   
 *   
 *   
 *     template <int dim>
 *     void PoissonSolver<dim>::setup_system()
 *     {
 * @endcode
 * 
 * First define the finite element space:
 * 
 * @code
 *       dof_handler.distribute_dofs(fe);
 *   
 *       std::cout << "   Number of degrees of freedom: " << dof_handler.n_dofs()
 *                 << std::endl;
 *   
 * @endcode
 * 
 * Then set up the main data structures that will hold the discrete problem:
 * 
 * @code
 *       {
 *         DynamicSparsityPattern dsp(dof_handler.n_dofs());
 *         DoFTools::make_sparsity_pattern(dof_handler, dsp);
 *         sparsity_pattern.copy_from(dsp);
 *   
 *         system_matrix.reinit(sparsity_pattern);
 *   
 *         solution.reinit(dof_handler.n_dofs());
 *         system_rhs.reinit(dof_handler.n_dofs());
 *       }
 *   
 * @endcode
 * 
 * And then define the tools to do point evaluation. We choose
 * a set of 13x13 points evenly distributed across the domain:
 * 
 * @code
 *       {
 *         const unsigned int n_points_per_direction = 13;
 *         const double       dx = 1. / (n_points_per_direction + 1);
 *   
 *         for (unsigned int x = 1; x <= n_points_per_direction; ++x)
 *           for (unsigned int y = 1; y <= n_points_per_direction; ++y)
 *             measurement_points.emplace_back(x * dx, y * dx);
 *   
 * @endcode
 * 
 * First build a full matrix of the evaluation process. We do this
 * even though the matrix is really sparse -- but we don't know
 * which entries are nonzero. Later, the `copy_from()` function
 * calls build a sparsity pattern and a sparse matrix from
 * the dense matrix.
 * 
 * @code
 *         Vector<double>     weights(dof_handler.n_dofs());
 *         FullMatrix<double> full_measurement_matrix(n_points_per_direction *
 *                                                      n_points_per_direction,
 *                                                    dof_handler.n_dofs());
 *   
 *         for (unsigned int index = 0; index < measurement_points.size(); ++index)
 *           {
 *             VectorTools::create_point_source_vector(dof_handler,
 *                                                     measurement_points[index],
 *                                                     weights);
 *             for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
 *               full_measurement_matrix(index, i) = weights(i);
 *           }
 *   
 *         measurement_sparsity.copy_from(full_measurement_matrix);
 *         measurement_matrix.reinit(measurement_sparsity);
 *         measurement_matrix.copy_from(full_measurement_matrix);
 *       }
 *   
 * @endcode
 * 
 * Next build the mapping from cell to the index in the 64-element
 * coefficient vector:
 * 
 * @code
 *       for (const auto &cell : triangulation.active_cell_iterators())
 *         {
 *           const unsigned int i = std::floor(cell->center()[0] * 8);
 *           const unsigned int j = std::floor(cell->center()[1] * 8);
 *   
 *           const unsigned int index = i + 8 * j;
 *   
 *           cell->set_user_index(index);
 *         }
 *   
 * @endcode
 * 
 * Finally prebuild the building blocks of the linear system as
 * discussed in the Readme file :
 * 
 * @code
 *       {
 *         const unsigned int dofs_per_cell = fe.dofs_per_cell;
 *   
 *         cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
 *         cell_rhs.reinit(dofs_per_cell);
 *   
 *         const QGauss<dim>  quadrature_formula(fe.degree+1);
 *         const unsigned int n_q_points = quadrature_formula.size();
 *   
 *         FEValues<dim> fe_values(fe,
 *                                 quadrature_formula,
 *                                 update_values | update_gradients |
 *                                   update_JxW_values);
 *   
 *         fe_values.reinit(dof_handler.begin_active());
 *   
 *         for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             {
 *               for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                 cell_matrix(i, j) +=
 *                   (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
 *                    fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
 *                    fe_values.JxW(q_index));           // dx
 *   
 *               cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
 *                               10.0 *                              // f(x_q)
 *                               fe_values.JxW(q_index));            // dx
 *             }
 *   
 *         VectorTools::interpolate_boundary_values(dof_handler,
 *                                                  0,
 *                                                  Functions::ZeroFunction<dim>(),
 *                                                  boundary_values);
 *       }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * Given that we have pre-built the matrix and right hand side contributions
 * for a (representative) cell, the function that assembles the matrix is
 * pretty short and straightforward:
 * 
 * @code
 *     template <int dim>
 *     void PoissonSolver<dim>::assemble_system(const Vector<double> &coefficients)
 *     {
 *       Assert(coefficients.size() == 64, ExcInternalError());
 *   
 *       system_matrix = 0;
 *       system_rhs    = 0;
 *   
 *       const unsigned int dofs_per_cell = fe.dofs_per_cell;
 *   
 *       std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *   
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           const double coefficient = coefficients(cell->user_index());
 *   
 *           cell->get_dof_indices(local_dof_indices);
 *           for (unsigned int i = 0; i < dofs_per_cell; ++i)
 *             {
 *               for (unsigned int j = 0; j < dofs_per_cell; ++j)
 *                 system_matrix.add(local_dof_indices[i],
 *                                   local_dof_indices[j],
 *                                   coefficient * cell_matrix(i, j));
 *   
 *               system_rhs(local_dof_indices[i]) += cell_rhs(i);
 *             }
 *         }
 *   
 *       MatrixTools::apply_boundary_values(boundary_values,
 *                                          system_matrix,
 *                                          solution,
 *                                          system_rhs);
 *     }
 *   
 *   
 * @endcode
 * 
 * The same is true for the function that solves the linear system:
 * 
 * @code
 *     template <int dim>
 *     void PoissonSolver<dim>::solve()
 *     {
 *       SparseILU<double> ilu;
 *       ilu.initialize(system_matrix);
 *       SolverControl control(100, 1e-10*system_rhs.l2_norm());
 *       SolverCG<> solver(control);
 *       solver.solve(system_matrix, solution, system_rhs, ilu);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * The following function outputs graphical data for the most recently
 * used coefficient and corresponding solution of the PDE. Collecting
 * the coefficient values requires translating from the 64-element
 * coefficient vector and the cells that correspond to each of these
 * elements. The rest remains pretty obvious, with the exception
 * of including the number of the current sample into the file name.
 * 
 * @code
 *     template <int dim>
 *     void
 *     PoissonSolver<dim>::output_results(const Vector<double> &coefficients) const
 *     {
 *       Vector<float> coefficient_values(triangulation.n_active_cells());
 *       for (const auto &cell : triangulation.active_cell_iterators())
 *         coefficient_values[cell->active_cell_index()] =
 *           coefficients(cell->user_index());
 *   
 *       DataOut<dim> data_out;
 *   
 *       data_out.attach_dof_handler(dof_handler);
 *       data_out.add_data_vector(solution, "solution");
 *       data_out.add_data_vector(coefficient_values, "coefficient");
 *   
 *       data_out.build_patches();
 *   
 *       std::ofstream output("solution-" +
 *                            Utilities::int_to_string(nth_evaluation, 10) + ".vtu");
 *       data_out.write_vtu(output);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * The following is the main function of this class: Given a coefficient
 * vector, it assembles the linear system, solves it, and then evaluates
 * the solution at the measurement points by applying the measurement
 * matrix to the solution vector. That vector of "measured" values
 * is then returned.
 *   
 
 * 
 * The function will also output the solution in a graphical format
 * if you un-comment the corresponding statement in the third
 * code block. However, you may end up with a very large amount
 * of data: This code is producing, at the minimum, 10,000 samples
 * and creating output for each one of them is surely more data
 * than you ever want to see!
 *   
 
 * 
 * At the end of the function, we output some timing information
 * every 10,000 samples.
 * 
 * @code
 *     template <int dim>
 *     Vector<double>
 *     PoissonSolver<dim>::evaluate(const Vector<double> &coefficients)
 *     {
 *       {
 *         TimerOutput::Scope section(timer, "Building linear systems");
 *         assemble_system(coefficients);
 *       }
 *   
 *       {
 *         TimerOutput::Scope section(timer, "Solving linear systems");
 *         solve();
 *       }
 *   
 *       Vector<double> measurements(measurement_matrix.m());
 *       {
 *         TimerOutput::Scope section(timer, "Postprocessing");
 *   
 *         measurement_matrix.vmult(measurements, solution);
 *         Assert(measurements.size() == measurement_points.size(),
 *                ExcInternalError());
 *   
 *         /*  output_results(coefficients);  */
 *       }
 *   
 *       ++nth_evaluation;
 *       if (nth_evaluation % 10000 == 0)
 *         timer.print_summary();
 *   
 *       return measurements;
 *     }
 *   } // namespace ForwardSimulator
 *   
 *   
 * @endcode
 * 
 * The following namespaces define the statistical properties of the Bayesian
 * inverse problem. The first is about the definition of the measurement
 * statistics (the "likelihood"), which we here assume to be a normal
 * distribution @f$N(\mu,\sigma I)@f$ with mean value @f$\mu@f$ given by the
 * actual measurement vector (passed as an argument to the constructor
 * of the `Gaussian` class and standard deviation @f$\sigma@f$.
 * 
 
 * 
 * For reasons of numerical accuracy, it is useful to not return the
 * actual likelihood, but its logarithm. This is because these
 * values can be very small, occasionally on the order of @f$e^{-100}@f$,
 * for which it becomes very difficult to compute accurate
 * values.
 * 
 * @code
 *   namespace LogLikelihood
 *   {
 *     class Interface
 *     {
 *     public:
 *       virtual double log_likelihood(const Vector<double> &x) const = 0;
 *   
 *       virtual ~Interface() = default;
 *     };
 *   
 *   
 *     class Gaussian : public Interface
 *     {
 *     public:
 *       Gaussian(const Vector<double> &mu, const double sigma);
 *   
 *       virtual double log_likelihood(const Vector<double> &x) const override;
 *   
 *     private:
 *       const Vector<double> mu;
 *       const double         sigma;
 *     };
 *   
 *     Gaussian::Gaussian(const Vector<double> &mu, const double sigma)
 *       : mu(mu)
 *       , sigma(sigma)
 *     {}
 *   
 *   
 *     double Gaussian::log_likelihood(const Vector<double> &x) const
 *     {
 *       Vector<double> x_minus_mu = x;
 *       x_minus_mu -= mu;
 *   
 *       return -x_minus_mu.norm_sqr() / (2 * sigma * sigma);
 *     }
 *   } // namespace LogLikelihood
 *   
 *   
 * @endcode
 * 
 * Next up is the "prior" imposed on the coefficients. We assume
 * that the logarithms of the entries of the coefficient vector
 * are all distributed as a Gaussian with given mean and standard
 * deviation. If the logarithms of the coefficients are normally
 * distributed, then this implies in particular that the coefficients
 * can only be positive, which is a useful property to ensure the
 * well-posedness of the forward problem.
 * 
 
 * 
 * For the same reasons as for the likelihood above, the interface
 * for the prior asks for returning the *logarithm* of the prior,
 * instead of the prior probability itself.
 * 
 * @code
 *   namespace LogPrior
 *   {
 *     class Interface
 *     {
 *     public:
 *       virtual double log_prior(const Vector<double> &x) const = 0;
 *   
 *       virtual ~Interface() = default;
 *     };
 *   
 *   
 *     class LogGaussian : public Interface
 *     {
 *     public:
 *       LogGaussian(const double mu, const double sigma);
 *   
 *       virtual double log_prior(const Vector<double> &x) const override;
 *   
 *     private:
 *       const double mu;
 *       const double sigma;
 *     };
 *   
 *     LogGaussian::LogGaussian(const double mu, const double sigma)
 *       : mu(mu)
 *       , sigma(sigma)
 *     {}
 *   
 *   
 *     double LogGaussian::log_prior(const Vector<double> &x) const
 *     {
 *       double log_of_product = 0;
 *   
 *       for (const auto &el : x)
 *         log_of_product +=
 *           -(std::log(el) - mu) * (std::log(el) - mu) / (2 * sigma * sigma);
 *   
 *       return log_of_product;
 *     }
 *   } // namespace LogPrior
 *   
 *   
 *   
 * @endcode
 * 
 * The Metropolis-Hastings algorithm requires a method to create a new sample
 * given a previous sample. We do this by perturbing the current (coefficient)
 * sample randomly using a Gaussian distribution centered at the current
 * sample. To ensure that the samples' individual entries all remain
 * positive, we use a Gaussian distribution in logarithm space -- in other
 * words, instead of *adding* a small perturbation with mean value zero,
 * we *multiply* the entries of the current sample by a factor that
 * is the exponential of a random number with mean zero. (Because the
 * exponential of zero is one, this means that the most likely factors
 * to multiply the existing sample entries by are close to one. And
 * because the exponential of a number is always positive, we never
 * get negative samples this way.)
 * 
 
 * 
 * But the Metropolis-Hastings sampler doesn't just need a perturbed
 * sample @f$y@f$ location given the current sample location @f$x@f$. It also
 * needs to know the ratio of the probability of reaching @f$y@f$ from
 * @f$x@f$, divided by the probability of reaching @f$x@f$ from @f$y@f$. If we
 * were to use a symmetric proposal distribution (e.g., a Gaussian
 * distribution centered at @f$x@f$ with a width independent of @f$x@f$), then
 * these two probabilities would be the same, and the ratio one. But
 * that's not the case for the Gaussian in log space. It's not
 * terribly difficult to verify that in that case, for a single
 * component the ratio of these probabilities is @f$y_i/x_i@f$, and
 * consequently for all components of the vector together, the
 * probability is the product of these ratios.
 * 
 * @code
 *   namespace ProposalGenerator
 *   {
 *     class Interface
 *     {
 *     public:
 *       virtual
 *       std::pair<Vector<double>,double>
 *       perturb(const Vector<double> &current_sample) const = 0;
 *   
 *       virtual ~Interface() = default;
 *     };
 *   
 *   
 *     class LogGaussian : public Interface
 *     {
 *     public:
 *       LogGaussian(const unsigned int random_seed, const double log_sigma);
 *   
 *       virtual
 *       std::pair<Vector<double>,double>
 *       perturb(const Vector<double> &current_sample) const override;
 *   
 *     private:
 *       const double         log_sigma;
 *       mutable std::mt19937 random_number_generator;
 *     };
 *   
 *   
 *   
 *     LogGaussian::LogGaussian(const unsigned int random_seed,
 *                              const double       log_sigma)
 *       : log_sigma(log_sigma)
 *     {
 *       random_number_generator.seed(random_seed);
 *     }
 *   
 *   
 *     std::pair<Vector<double>,double>
 *     LogGaussian::perturb(const Vector<double> &current_sample) const
 *     {
 *       Vector<double> new_sample = current_sample;
 *       double         product_of_ratios = 1;
 *       for (auto &x : new_sample)
 *         {
 *           const double rnd = std::normal_distribution<>(0, log_sigma)(random_number_generator);
 *           const double exp_rnd = std::exp(rnd);
 *           x *= exp_rnd;
 *           product_of_ratios *= exp_rnd;
 *         }
 *   
 *       return {new_sample, product_of_ratios};
 *     }
 *   
 *   } // namespace ProposalGenerator
 *   
 *   
 * @endcode
 * 
 * The last main class is the Metropolis-Hastings sampler itself.
 * If you understand the algorithm behind this method, then
 * the following implementation should not be too difficult
 * to read. The only thing of relevance is that descriptions
 * of the algorithm typically ask whether the *ratio* of two
 * probabilities (the "posterior" probabilities of the current
 * and the previous samples, where the "posterior" is the product of the
 * likelihood and the prior probability) is larger or smaller than a
 * randomly drawn number. But because our interfaces return the
 * *logarithms* of these probabilities, we now need to take
 * the ratio of appropriate exponentials -- which is made numerically
 * more stable by considering the exponential of the difference of
 * the log probabilities. The only other slight complication is that
 * we need to multiply this ratio by the ratio of proposal probabilities
 * since we use a non-symmetric proposal distribution.
 * 
 
 * 
 * Finally, we note that the output is generated with 7 digits of
 * accuracy. (The C++ default is 6 digits.) We do this because,
 * as shown in the paper, we can determine the mean value of the
 * probability distribution we are sampling here to at least six
 * digits of accuracy, and do not want to be limited by the precision
 * of the output.
 * 
 * @code
 *   namespace Sampler
 *   {
 *     class MetropolisHastings
 *     {
 *     public:
 *       MetropolisHastings(ForwardSimulator::Interface &       simulator,
 *                          const LogLikelihood::Interface &    likelihood,
 *                          const LogPrior::Interface &         prior,
 *                          const ProposalGenerator::Interface &proposal_generator,
 *                          const unsigned int                  random_seed,
 *                          const std::string &                 dataset_name);
 *   
 *       void sample(const Vector<double> &starting_guess,
 *                   const unsigned int    n_samples);
 *   
 *     private:
 *       ForwardSimulator::Interface &       simulator;
 *       const LogLikelihood::Interface &    likelihood;
 *       const LogPrior::Interface &         prior;
 *       const ProposalGenerator::Interface &proposal_generator;
 *   
 *       std::mt19937 random_number_generator;
 *   
 *       unsigned int sample_number;
 *       unsigned int accepted_sample_number;
 *   
 *       std::ofstream output_file;
 *   
 *       void write_sample(const Vector<double> &current_sample,
 *                         const double          current_log_likelihood);
 *     };
 *   
 *   
 *     MetropolisHastings::MetropolisHastings(
 *       ForwardSimulator::Interface &       simulator,
 *       const LogLikelihood::Interface &    likelihood,
 *       const LogPrior::Interface &         prior,
 *       const ProposalGenerator::Interface &proposal_generator,
 *       const unsigned int                  random_seed,
 *       const std::string &                 dataset_name)
 *       : simulator(simulator)
 *       , likelihood(likelihood)
 *       , prior(prior)
 *       , proposal_generator(proposal_generator)
 *       , sample_number(0)
 *       , accepted_sample_number(0)
 *     {
 *       output_file.open("samples-" + dataset_name + ".txt");
 *       output_file.precision(7);
 *   
 *       random_number_generator.seed(random_seed);
 *     }
 *   
 *   
 *     void MetropolisHastings::sample(const Vector<double> &starting_guess,
 *                                     const unsigned int    n_samples)
 *     {
 *       std::uniform_real_distribution<> uniform_distribution(0, 1);
 *   
 *       Vector<double> current_sample = starting_guess;
 *       double         current_log_posterior =
 *         (likelihood.log_likelihood(simulator.evaluate(current_sample)) +
 *          prior.log_prior(current_sample));
 *   
 *       ++sample_number;
 *       ++accepted_sample_number;
 *       write_sample(current_sample, current_log_posterior);
 *   
 *       for (unsigned int k = 1; k < n_samples; ++k, ++sample_number)
 *         {
 *           std::pair<Vector<double>,double>
 *             perturbation = proposal_generator.perturb(current_sample);
 *           const Vector<double> trial_sample                   = std::move (perturbation.first);
 *           const double         perturbation_probability_ratio = perturbation.second;
 *   
 *           const double trial_log_posterior =
 *             (likelihood.log_likelihood(simulator.evaluate(trial_sample)) +
 *              prior.log_prior(trial_sample));
 *   
 *           if (std::exp(trial_log_posterior - current_log_posterior) * perturbation_probability_ratio
 *               >=
 *               uniform_distribution(random_number_generator))
 *             {
 *               current_sample        = trial_sample;
 *               current_log_posterior = trial_log_posterior;
 *   
 *               ++accepted_sample_number;
 *             }
 *   
 *           write_sample(current_sample, current_log_posterior);
 *         }
 *     }
 *   
 *   
 *   
 *     void MetropolisHastings::write_sample(const Vector<double> &current_sample,
 *                                           const double current_log_posterior)
 *     {
 *       output_file << current_log_posterior << '\t';
 *       output_file << accepted_sample_number << '\t';
 *       for (const auto &x : current_sample)
 *         output_file << x << ' ';
 *       output_file << '\n';
 *   
 *       output_file.flush();
 *     }
 *   } // namespace Sampler
 *   
 *   
 * @endcode
 * 
 * The final function is `main()`, which simply puts all of these pieces
 * together into one. The "exact solution", i.e., the "measurement values"
 * we use for this program are tabulated to make it easier for other
 * people to use in their own implementations of this benchmark. These
 * values created using the same main class above, but using 8 mesh
 * refinements and using a Q3 element -- i.e., using a much more accurate
 * method than the one we use in the forward simulator for generating
 * samples below (which uses 5 global mesh refinement steps and a Q1
 * element). If you wanted to regenerate this set of numbers, then
 * the following code snippet would do that:
 * <div class=CodeFragmentInTutorialComment>
 * @code
 *  /* Set the exact coefficient: */
 *  Vector<double> exact_coefficients(64);
 *  for (auto &el : exact_coefficients)
 *    el = 1.;
 *  exact_coefficients(9) = exact_coefficients(10) = exact_coefficients(17) =
 *    exact_coefficients(18)                       = 0.1;
 *  exact_coefficients(45) = exact_coefficients(46) = exact_coefficients(53) =
 *    exact_coefficients(54)                        = 10.;
 *
 
 *  /* Compute the "correct" solution vector: */
 *  const Vector<double> exact_solution =
 *    ForwardSimulator::PoissonSolver<2>(/* global_refinements = */ 8,
 *                                       /* fe_degree = */ 3,
 *                                       /* prefix = */ "exact")
 *      .evaluate(exact_coefficients);
 * @endcode
 * </div>
 * 
 * @code
 *   int main()
 *   {
 *     const bool testing = true;
 *   
 * @endcode
 * 
 * Run with one thread, so as to not step on other processes
 * doing the same at the same time. It turns out that the problem
 * is also so small that running with more than one thread
 * *increases* the runtime.
 * 
 * @code
 *     MultithreadInfo::set_thread_limit(1);
 *   
 *     const unsigned int random_seed  = (testing ? 1U : std::random_device()());
 *     const std::string  dataset_name = Utilities::to_string(random_seed, 10);
 *   
 *     const Vector<double> exact_solution(
 *       {   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                         });
 *   
 * @endcode
 * 
 * Now run the forward simulator for samples:
 * 
 * @code
 *     ForwardSimulator::PoissonSolver<2> laplace_problem(
 *       /* global_refinements = */ 5,
 *       /* fe_degree = */ 1,
 *       dataset_name);
 *     LogLikelihood::Gaussian        log_likelihood(exact_solution, 0.05);
 *     LogPrior::LogGaussian          log_prior(0, 2);
 *     ProposalGenerator::LogGaussian proposal_generator(
 *       random_seed, 0.09); /* so that the acceptance ratio is ~0.24 */
 *     Sampler::MetropolisHastings sampler(laplace_problem,
 *                                         log_likelihood,
 *                                         log_prior,
 *                                         proposal_generator,
 *                                         random_seed,
 *                                         dataset_name);
 *   
 *     Vector<double> starting_coefficients(64);
 *     for (auto &el : starting_coefficients)
 *       el = 1.;
 *     sampler.sample(starting_coefficients,
 *                    (testing ? 250 * 40 /* takes 10 seconds */
 *                               :
 *                               100000000 /* takes 1.5 days */
 *                     ));
 *   }
 * @endcode
 
 
*/