```
# default_exp twoDim
```


```
#hide
import matplotlib.pyplot as plt
import seaborn as sns
import matplotlib.cm as cm
import matplotlib.colors as colors
import matplotlib.colorbar as colorbar
import matplotlib.tri as tri
plt.rcParams['figure.figsize'] = (10,6)
import sympy; sympy.init_printing()
# code for displaying matrices nicely
def display_matrix(m):
    display(sympy.Matrix(m))
```

# twoDim
> Code for a 2-D problem.


```
#hide
from nbdev.showdoc import *
```

# 2 dimensional case (PDE)

We consider the following 2-D problem:

$$\nabla\cdot\left(\kappa(x)\nabla u(x)\right)=f(x) \quad\forall x\in D=[0,1]^{2}$$
$$u(x)=0\quad\forall x\in\partial D$$
where here $f$ is again a random forcing term, assumed to be a GP in this work.

## Variational formulation

The variational formulation is given by:

$$a(u,v)=L(v)$$

where:

$$a(u,v)=\int_{D}\nabla u\cdot\left(\kappa\nabla u\right)dx$$

and

$$L(v)=\int_{D}fvdx$$

We will make the following choices for $\kappa,f$:

$$\kappa(x)=1$$

$$f\sim\mathcal{G}\mathcal{P}(\bar{f},k_{f})$$

$$\bar{f}(x)=1$$

$$ k_{f}(x,y) = \sigma_f^{2}\exp\left(-\frac{\|x-y\|^2}{2l_f^2}\right)$$

$$ \sigma_{f} = 0.1$$

$$ l_f = 0.4 $$

where $\|\cdot\|$ is the usual Euclidean norm.

Since we do not have access to a suitable Green's function for this problem, we will have to estimate the rate of convergence of the statFEM prior and posterior by comparing them on a sequence of refined meshes. More details on this will follow later. Thus, we need similar code as for the 1-D problem.

## statFEM prior mean

We will again utilise FEniCS to obtain the statFEM prior mean. For this purpose, we create a function `mean_assembler` which will assemble the mean for the statFEM prior.


```
#export
from dolfin import *
import numpy as np
from scipy import integrate
from scipy.spatial.distance import cdist
from scipy.linalg import sqrtm
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve
from scipy.interpolate import interp1d
from joblib import Parallel, delayed
import multiprocessing

# code to assemble the mean for a given mesh size
def mean_assembler(h,f_bar):
    "This function assembles the mean for the statFEM prior for our 2-D problem"
    # get size of the grid
    J = int(np.round(1/h))

    # set up the mesh and function space for FEM
    mesh = UnitSquareMesh(J,J)
    V = FunctionSpace(mesh,'Lagrange',1)

    # set up boundary condition
    def boundary(x, on_boundary):
        return on_boundary

    bc = DirichletBC(V, 0.0, boundary)
    # set up the functions κ and f
    κ = Constant(1.0)
    f = f_bar

    # set up the bilinear form for the variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    a = inner(κ*grad(u),grad(v))*dx

    # set up the linear form
    L = f*v*dx

    # solve the variational problem
    μ = Function(V)
    solve(a == L, μ, bc)

    return μ
```

`mean_assembler` takes in the mesh size `h` and the mean function `f_bar` for the forcing and computes the mean of the approximate statFEM prior, returning this as a FEniCS function.

> Important: `mean_assembler` requires `f_bar` to be represented as a FEniCS function/expression/constant.

Let's check that this is working:


```
h = 0.1
f_bar = Constant(1.0)
μ = mean_assembler(h,f_bar)
μ
```


```
# check the type of μ
assert type(μ) == function.function.Function
```

Let's plot $\mu$:


```
#hide_input
# use FEniCS to plot μ
plot(μ)
plt.xlabel(r'$x$')
plt.ylabel(r'$y$')
plt.title(r'Plot of statFEM mean for $h=%.2f$'%h)
plt.show()
```

## statFEM prior covariance

We will also utilise FEniCS again to obtain an approximation of our statFEM covariance function.

The statFEM covariance can be approximated as follows:

$$c_u^{\text{FEM}}(x,y)\approx\sum_{i,j=1}^{J}\varphi_{i}(x)Q_{ij}\varphi_{j}(y)$$

where $Q=A^{-1}MC_{f}M^{T}A^{-T}$ and where the $\{\varphi_{i}\}_{i=1}^{J}$ are the FE basis functions corresponding to the interior nodes of our domain.

with $C_f$ being the kernel matrix of $f$ (evaluated on the FEM grid).

As we will be comparing the statFEM covariance functions for finer and finer FE mesh sizes we will need to be able to assemble the statFEM covariance function on a grid. As discussed in <a href="/statFEM/oneDim.html#Difference-between-true-prior-covariance-and-statFEM-prior-covariance"><code>oneDim</code></a>, we can assemble such covariance matrices in a very efficient manner. The code remains largely the same as in the 1-D case and so we do not go into as much detail here.

We start by creating a function `kernMat` which assembles the covariance matrix corresponding to a covariance function `k` on a grid `grid`.


```
#export
def kernMat(k,grid,parallel=True,translation_inv=False):
    "Function to compute the covariance matrix $K$ corresponding to the covariance kernel $k$ on a grid. This matrix has $ij$-th entry $K_{ij}=k(x_i,x_j)$ where $x_i$ is the $i$-th point of the grid."
    # get the length of the grid
    n = grid.shape[0]
    # preallocate an n x n array of zeros to hold the cov matrix
    K = np.zeros((n,n))

    # check if the cov matrix should be computed in parallel
    if parallel:
        # compute the cov matrix in parallel by computing the upper triangular part column by column
        # set up function to compute the ith column of the upper triangular part:
        def processInput(i):
            return np.array([k(grid[i,:],grid[j,:]) for j in range(i,n)])
        
        # get the number of cpu cores present and compute the upper triangular columns in parallel
        num_cores = multiprocessing.cpu_count()
        results = Parallel(n_jobs=num_cores)(delayed(processInput)(i) for i in range(n))

        # store the results in the appropriate positions in K
        #for (i,v) in enumerate(results[0:n-1]):
        for (i,v) in enumerate(results): # is this correct???
            K[i,i:] = v
        
        # only the upper triangular part has been formed, so use the symmetry of the cov mat to get full K:
        K = K + K.T - np.diag(K.diagonal())
        return K
    elif translation_inv:
        # reshape grid so that it has correct dimensions
        grid = grid.reshape(n,-1)
        
        # compute the distance matrix D
        D = cdist(grid,grid)

        # evaluate the kernel function using D
        K = k(D)
        return K
    else:
        # compute the cov mat using a nested for loop
        for i in range(n):
            for j in range(i,n):
                K[i,j] = k(grid[i,:],grid[j,:])
        K = K + K.T - np.diag(K.diagonal())
        return K    
```

> Note: This function takes in two optional boolean arguments `parallel` and `translation_inv`. The first of these specifies whether or not the cov matrix should be computed in parallel and the second specifies whether or not the cov kernel is translation invariant. If it is, the covariance matrix is computed more efficiently using the `cdist` function from scipy.

Let's quickly test if this function is working, by computing the cov matrix for white noise, which has kernel function $k(x,y)=\delta(x-y)$. For a grid of length $N$ this should be the $N\times N$ identity matrix.


```
# set up the kernel function
# set up tolerance for comparison
tol = 1e-16
def k(x,y):
    if (np.abs(x-y) < tol).all():
        # x == y within the tolerance
        return 1.0
    else:
        # x != y within the tolerance
        return 0.0

# set up grid
n = 21
x_range = np.linspace(0,1,n)
grid = np.array([[x,y] for x in x_range for y in x_range])
N = len(grid) # get length of grid (N=n^2)
K = kernMat(k,grid,True,False) # parallel mode

# check that this is the N x N identity matrix
assert (K == np.eye(N)).all()
```

We now create a function `BigPhiMat` to utilise FEniCS to efficiently compute the matrix $\boldsymbol{\Phi}$ defined above.


```
#export
def BigPhiMat(J,grid):
    "Function to compute the $\Phi$ matrix."
    # create the FE mesh and function space
    mesh = UnitSquareMesh(J,J)
    V = FunctionSpace(mesh,'Lagrange',1)
    # get the tree for the mesh
    tree = mesh.bounding_box_tree()
    # set up a function to compute the ith column of Phi corresponding to the ith grid point
    def Φ(i):
        x = grid[i]
        cell_index = tree.compute_first_entity_collision(Point(*x))
        cell = Cell(mesh,cell_index)
        cell_global_dofs = V.dofmap().cell_dofs(cell_index)
        vertex_coordinates = cell.get_vertex_coordinates()
        cell_orientation = cell.orientation()
        data = V.element().evaluate_basis_all(x,vertex_coordinates,cell_orientation)
        return (data,cell_global_dofs,i*np.ones_like(cell_global_dofs))
    # compute all the columns of Phi using the function above
    res = [Φ(i) for i in range(len(grid))]
    # assemble the sparse matrix Phi using the results
    data = np.hstack([res[i][0] for i in range(len(grid))])
    row = np.hstack([res[i][1] for i in range(len(grid))])
    col = np.hstack([res[i][2] for i in range(len(grid))])
    return csr_matrix((data,(row,col)),shape=(V.dim(),len(grid)))
```

`BigPhiMat` takes in two arguments: `J`, which controls the FE mesh size ($h=1/J^{2}$), and `grid` which is the grid in the definition of $\boldsymbol{\Phi}$. `BigPhiMat` returns $\boldsymbol{\Phi}$ as a sparse `csr_matrix` for memory efficiency.

> Note: Note that since FEniCS works with the FE functions corresponding to all the FE dofs and our statFEM cov matix only uses the FE functions corresponding to non-boundary dofs we need to account for this in the code. See the source code for `BigPhiMat` to see how this is done.

We now create a function `cov_asssembler` which assembles the approximate FEM covariance matrix on the grid.


```
#export
# function to assemble the fem covariance
def cov_assembler(J,k_f,grid,parallel,translation_inv):
    "Function to assemble the approximate FEM covariance matrix on the reference grid."
    
    # set up mesh and function space
    mesh = UnitSquareMesh(J,J)
    V = FunctionSpace(mesh,'Lagrange',1)
    
    # set up FE grid
    x_grid = V.tabulate_dof_coordinates()
    
    # set up boundary condition
    def boundary(x, on_boundary):
        return on_boundary
    
    bc = DirichletBC(V, 0.0, boundary)
    
    # get the boundary and interior dofs
    bc_dofs = bc.get_boundary_values().keys()
    first, last = V.dofmap().ownership_range()
    all_dofs = range(last - first)
    interior_dofs = list(set(all_dofs) - set(bc_dofs))
    bc_dofs = list(set(bc_dofs))
    
    # set up the function κ
    κ = Constant(1.0)

    # get the mass and stiffness matrices as sparse csr_matrices
    u = TrialFunction(V)
    v = TestFunction(V)
    
    mass_form = u*v*dx
    a = inner(κ*grad(u),grad(v))*dx
    
    M = assemble(mass_form)
    A = assemble(a)
    M = as_backend_type(M).mat()
    A = as_backend_type(A).mat()
    M = csr_matrix(M.getValuesCSR()[::-1],shape=M.size)
    A = csr_matrix(A.getValuesCSR()[::-1],shape=A.size)
    
    # extract the submatrices corresponding to the interior dofs
    M = M[interior_dofs,:][:,interior_dofs]
    A = A[interior_dofs,:][:,interior_dofs]
    
    # get the forcing cov matrix on the interior nodes of the grid
    Σ_int = kernMat(k_f,x_grid[interior_dofs],parallel,translation_inv)
    
    # form the matrix Q in the defintion of the approximate FEM cov mat
    # Note: overwrite Σ_int for memory efficiency.
#     Σ_int = M @ Σ_int @ M.T
    Σ_int = Σ_int @ M.T
    Σ_int = M @ Σ_int

    Σ_int = spsolve(A,Σ_int)
    Σ_int = spsolve(A,Σ_int.T).T
    
    # ensure Σ_int is symmetric
    Σ_int = 0.5*(Σ_int + Σ_int.T)
    
    # get big phi matrix on the grid (extracting only the rows corresponding to the
    # interior dofs)
    Phi = BigPhiMat(J,grid)[interior_dofs,:]
    #print("Computed Phi")
    
    # assemble cov mat on grid using Phi and Σ_int
    Σ = Phi.T @ Σ_int @ Phi
    
    # ensure Σ is symmetric and return
    Σ = 0.5*(Σ + Σ.T)
    return Σ
```

`cov_assembler` takes in several arguments which are explained below:

- `J`: controls the FE mesh size ($h=1/J^{2})$
- `k_f`: the covariance function for the forcing $f$
- `grid`: the reference grid where the FEM cov matrix should be computed on
- `parallel`: boolean argument indicating whether the intermediate computation of $C_f$ should be done in parallel 
- `translation_inv`: boolean argument indicating whether the intermediate computation of $C_f$ should be computed assuming `k_f` is translation invariant or not

As a quick demonstration that the code is working, we will the statFEM cov matrix for a relatively coarse grid.


```
# set up kernel function for forcing
f_bar = Constant(1.0)

l_f = 0.4
σ_f = 0.1

def k_f(x):
    return (σ_f**2)*np.exp(-(x**2)/(2*(l_f**2)))

# set up grid
n = 21
x_range = np.linspace(0,1,n)
grid = np.array([[x,y] for x in x_range for y in x_range])

# get the statFEM grid for a particular choice of J
J = 10
Σ = cov_assembler(J,k_f,grid,False,True)
```

Let's plot a heatmap of the statFEM cov matrix:


```
#hide_input
sns.heatmap(Σ,cbar=True,
              annot=False,
              xticklabels=False,
              yticklabels=False,
              cmap=cm.viridis)
plt.title('Heat map of statFEM covariance matrix')
plt.show()
```

> Note: that the banded structure in the above statFEM covariance matrix is due to the internal ordering of the FE grid in FEniCS.

## statFEM posterior mean

The statFEM posterior from incorporating sensor readings has the same form as given in <a href="/statFEM/oneDim.html#Posterior-from-incorporating-sensor-readings"><code>oneDim</code></a>. We will thus require very similar code as to the 1-D case. We start by creating a function `m_post` which evaluates the posterior mean at a given point.


```
#export
def m_post(x,m,c,v,Y,B):
    "This function evaluates the posterior mean at the point $x$."
    m_vect = np.array([m(y_i) for y_i in Y]).flatten()
    c_vect = c(x).flatten()
    
    # compute the update term
    update = c_vect @ np.linalg.solve(B,m_vect-v)
    
    # return m_post
    return (m(x) - update)
```

`m_post` takes in several arguments which are explained below:

- `x`: point where the posterior mean will be evaluated
- `m`: function which computes the prior mean at a given point y
- `c`: function which returns the vector (c(x,y)) for y in Y (note: c is the prior covariance function)
- `v`: vector of noisy sensor readings
- `Y`: vector of sensor locations
- `B`: the matrix $\epsilon^{2}I+C_Y$ to be inverted in order to obtain the posterior

We now require code to generate samples from a GP with mean $m$ and cov function $k$ on a grid. We write the function `sample_gp` for this purpose.


```
#export
def sample_gp(n_sim,m,k,grid,par,trans,tol=1e-9):
    "Function to sample a GP with mean $m$ and cov $k$ on a grid."
    
    # get length of grid
    d = len(grid)
    
    # construct mean vector
    μ = np.array([m(x) for x in grid]).reshape(d,1)
    
    # construct covariance matrix
    Σ = kernMat(k,grid,parallel=par,translation_inv=trans)
    
    # construct the cholesky decomposition Σ = GG^T
    # we add a small diagonal perturbation to Σ to ensure it 
    # strictly positive definite
    G = np.linalg.cholesky(Σ + tol * np.eye(d))
    
    # draw iid standard normal random vectors
    Z = np.random.normal(size=(d,n_sim))
    
    # construct samples from GP(m,k)
    Y = G@Z + np.tile(μ,n_sim)
    
    # return the sampled fields
    return Y
```

`sample_gp` takes in several arguments which are explained below:

- `n_sim`: number of trajectories to be sampled
- `m`: mean function for the GP
- `k`: cov function for the GP
- `grid`: grid of points on which to sample the GP
- `par`: boolean argument indicating whether the computation of the cov matrix should be done in parallel
- `trans`: boolean argument indicating whether the computation of the cov matrix should be computed assuming `k` is translation invariant or not
- `tol`: controls the size of the tiny diagonal perturbation added to cov matrix to ensure it is strictly positive definite (defaults to `1e-9`)

As a quick demonstration that the code is working lets generate 2 realisations of white noise, using the kernel `k` from one of the previous tests and plot a heatmap of these random fields side-by-side.


```
#hide_input
n = 41
x_range = np.linspace(0,1,n)
grid = np.array([[x,y] for x in x_range for y in x_range])

# set up mean
def m(x):
    return 0.0

np.random.seed(23534)
samples = sample_gp(2,m,k,grid,True,False)

sample_1 = samples[:,0].flatten()
sample_2 = samples[:,1].flatten()

vmin = min(sample_1.min(),sample_2.min())
vmax = max(sample_1.max(),sample_2.max())
cmap = cm.jet
norm = colors.Normalize(vmin=vmin,vmax=vmax)

x = grid[:,0].flatten()
y = grid[:,1].flatten()
triang = tri.Triangulation(x,y)

plt.rcParams['figure.figsize'] = (12,6)
fig, axs = plt.subplots(ncols=3, gridspec_kw=dict(width_ratios=[4,4,0.2]))
axs[0].tricontourf(triang,sample_1.flatten(),cmap=cmap)
axs[1].tricontourf(triang,sample_2.flatten(),cmap=cmap)
cb = colorbar.ColorbarBase(axs[2],cmap=cmap,norm=norm)
fig.suptitle('Realisations of white-noise fields')
plt.show()
```

Let's also quickly generate 2 realisations for the kernel `k_f` above:


```
#hide_input
np.random.seed(534)
samples = sample_gp(2,m,k_f,grid,False,True)

sample_1 = samples[:,0].flatten()
sample_2 = samples[:,1].flatten()

vmin = min(sample_1.min(),sample_2.min())
vmax = max(sample_1.max(),sample_2.max())
cmap = cm.jet
norm = colors.Normalize(vmin=vmin,vmax=vmax)

x = grid[:,0].flatten()
y = grid[:,1].flatten()
triang = tri.Triangulation(x,y)

plt.rcParams['figure.figsize'] = (12,6)
fig, axs = plt.subplots(ncols=3, gridspec_kw=dict(width_ratios=[4,4,0.2]))
axs[0].tricontourf(triang,sample_1.flatten(),cmap=cmap)
axs[1].tricontourf(triang,sample_2.flatten(),cmap=cmap)
cb = colorbar.ColorbarBase(axs[2],cmap=cmap,norm=norm)
fig.suptitle(r'Realisations of random fields with covariance $k_f$')
plt.show()
```

The next bit of code we require is code to generate noisy sensor readings from our system. We write the function `gen_sensor` for this purpose.


```
#export
def gen_sensor(ϵ,m,k,Y,J,par,trans,tol=1e-9,require=False):
    "Function to generate noisy sensor readings of the solution u on a sensor grid Y."
    
    # get number of sensors from the sensor grid Y
    s = len(Y)
    
    # create FEM space and grid
    mesh = UnitSquareMesh(J,J)
    V = FunctionSpace(mesh,'Lagrange',1)
    grid = V.tabulate_dof_coordinates()
    
    # sample a single f on the grid
    f_sim = sample_gp(1,m,k,grid,par=par,trans=trans,tol=tol)
    
    # set up a FEM function for this realisation
    f = Function(V)
    f.vector().set_local(f_sim.flatten())
    
    # use FENICS to find the corresponding solution u
    # set up boundary condition
    def boundary(x, on_boundary):
        return on_boundary

    bc = DirichletBC(V, 0.0, boundary)
    
    # set up the function κ
    κ = Constant(1.0)
    
    # set up the bilinear form for the variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    a = inner(κ*grad(u),grad(v))*dx
    
    # set up the linear form
    L = f*v*dx
    
    # solve the variational problem
    u_sol = Function(V)
    solve(a == L, u_sol, bc)
    
    # get solution on grid Y:
    u_Y = np.array([u_sol(y_i) for y_i in Y])
    
    # add N(0,ϵ^2) to each evaluation point
    u_S = u_Y + ϵ*np.random.normal(size=s)
    
    if require:
        return u_S, f_sim, u_sol
    else:
        return u_S
```

`gen_sensor` takes in several arguments which are explained below:

- `ϵ`: controls the amount of sensor noise
- `m`: mean function for the forcing f
- `k`: cov function for the forcing f
- `Y`: vector of sensor locations
- `J`: controls the FE mesh size ($h=1/J^{2}$)
- `par`: boolean argument indicating whether the computation of the forcing cov matrix should be done in parallel
- `trans`: boolean argument indicating whether the computation of the forcing cov matrix should be computed assuming `k` is translation invariant or not
- `tol`: controls the size of the tiny diagonal perturbation added to forcing cov matrix to ensure it is strictly positive definite (defaults to `1e-9`)
- `require` : boolean argument indicating whether or not to also return the realisation of the forcing `f_sim` and the FEniCS solution `u_sol` (defaults to `False`)

> Warning: Since we do not have access to the true solution we must use FEniCS to get the solution for our system. Thus, one must choose a small enough `J` in `gen_sensor` above to ensure we get realistic noisy sensor readings.

Let's demonstrate that this code is working, by generating $s=25$ sensor observations with the sensors equally space in the domain $D$.


```
# set up mean function for forcing
def m_f(x):
    return 1.0

# set up sensor grid and sensor noise level
ϵ = 0.2
s = 25
s_sqrt = int(np.round(np.sqrt(s)))
Y_range = np.linspace(0.01,0.99,s_sqrt)
Y = np.array([[x,y] for x in Y_range for y in Y_range])
J_fine = 100 # FE mesh size to compute solution on

# generate the sensor observations
np.random.seed(235)
v_dat = gen_sensor(ϵ,m_f,k_f,Y,J_fine,False,True)
```


```
#export
class MyExpression(UserExpression):
    "Class to allow users to user their own functions to create a FEniCS UserExpression."
    def eval(self, value, x):
        value[0] = self.f(x)
    def value_shape(self):
        return ()
```


```
show_doc(MyExpression,title_level=4)
```


<h4 id="MyExpression" class="doc_header"><code>class</code> <code>MyExpression</code><a href="" class="source_link" style="float:right">[source]</a></h4>

> <code>MyExpression</code>(**\*`args`**, **\*\*`kwargs`**) :: `UserExpression`

Class to allow users to user their own functions to create a FEniCS UserExpression.


We now require code which will create the matrix $C_Y,h$ and the function $\mathbf{c}^{(h)}$ required for the statFEM posterior mean. We will create the function `fem_cov_assembler_post` for this purpose. 


```
#export
def fem_cov_assembler_post(J,k_f,Y,parallel,translation_inv):
    "Function to create the matrix $C_{Y,h}$ and the vector function $c^{(h)}$ required for the statFEM posterior mean."
    
    # set up mesh and function space
    mesh = UnitSquareMesh(J,J)
    V = FunctionSpace(mesh,'Lagrange',1)
    tree = mesh.bounding_box_tree()
    
    # set up grid
    x_grid = V.tabulate_dof_coordinates()
    
    # set up boundary condition
    def boundary(x, on_boundary):
        return on_boundary
    
    bc = DirichletBC(V, 0.0, boundary)
    
    # get the boundary and interior dofs
    bc_dofs = bc.get_boundary_values().keys()
    first, last = V.dofmap().ownership_range()
    all_dofs = range(last - first)
    interior_dofs = list(set(all_dofs) - set(bc_dofs))
    bc_dofs = list(set(bc_dofs))
    
    # set up the function κ
    κ = Constant(1.0)

    # get the mass and stiffness matrices
    u = TrialFunction(V)
    v = TestFunction(V)
    
    mass_form = u*v*dx
    a = inner(κ*grad(u),grad(v))*dx
    
    M = assemble(mass_form)
    A = assemble(a)
    M = as_backend_type(M).mat()
    A = as_backend_type(A).mat()
    M = csr_matrix(M.getValuesCSR()[::-1],shape=M.size)
    A = csr_matrix(A.getValuesCSR()[::-1],shape=A.size)
    
    # extract the submatrices corresponding to the interior dofs
    M = M[interior_dofs,:][:,interior_dofs]
    A = A[interior_dofs,:][:,interior_dofs]
    
    # get the forcing cov matrix on the interior nodes of the grid
    Σ_int = kernMat(k_f,x_grid[interior_dofs],parallel,translation_inv)

    # form the matrix Q in the defintion of the approximate FEM cov mat
    # Note: overwrite Σ_int for memory efficiency
    Σ_int = M @ Σ_int @ M.T

    Σ_int = spsolve(A,Σ_int)
    Σ_int = spsolve(A,Σ_int.T).T
    
    # ensure Σ_int is symmetric
    Σ_int = 0.5*(Σ_int + Σ_int.T)
    
    # get big phi matrix on the sensor grid (only need the interior dofs)
    Phi = BigPhiMat(J,Y)[interior_dofs,:]

    # assemble the FEM cov mat on the sensor grid and ensure it is symmetric
    Σ_s = Phi.T @ Σ_int @ Phi
    Σ_s = 0.5*(Σ_s + Σ_s.T)
    
    # set up function to yield the vector (c(x,y)) for y in Y
    def Φ(x):
        cell_index = tree.compute_first_entity_collision(Point(*x))
        cell_global_dofs = V.dofmap().cell_dofs(cell_index)
        cell = Cell(mesh, cell_index)
        vertex_coordinates = cell.get_vertex_coordinates()
        cell_orientation = cell.orientation()
        data = V.element().evaluate_basis_all(x,vertex_coordinates,cell_orientation)
        col = np.zeros_like(cell_global_dofs)
        res = csr_matrix((data,(cell_global_dofs,col)),shape=(V.dim(),1))[interior_dofs,:]
        return res
    
    def c_fem(x):
        return Φ(x).T @ Σ_int @ Phi
    
    #return Σ_s and c_fem
    return Σ_s, c_fem
```

`fem_cov_assembler_post` takes in several arguments which are explained below:

- `J`: controls the FE mesh size ($h=1/J^2$)
- `k_f`: the covariance function for the forcing $f$
- `Y`: vector of sensor locations
- `parallel`: boolean argument indicating whether the computation of the forcing cov mat should be done in parallel
- `translation_inv`: boolean argument indicating whether the computation of the forcing cov mat should be computed assuming `k_f` is translation invariant or not

With all of this code in place we can now finally write the function `m_post_fem_assmebler` which will assemble the statFEM posterior mean function.


```
#export
def m_post_fem_assembler(J,f_bar,k_f,ϵ,Y,v_dat,par=False,trans=True):
    "Function to assemble the statFEM posterior mean function."
    
    # get number of sensors
    s = len(Y)
    
    # set up mesh and function space
    mesh = UnitSquareMesh(J,J)
    V = FunctionSpace(mesh,'Lagrange',1)
    
    # set up boundary condition
    def boundary(x, on_boundary):
        return on_boundary
    
    bc = DirichletBC(V, 0.0, boundary)
    
    # set up the functions κ and f
    κ = Constant(1.0)
    f = f_bar
    
    # set up the bilinear form for the variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    a = inner(κ*grad(u),grad(v))*dx
    
    # set up linear form
    L = f*v*dx
    
    # solve the variational problem
    μ_fem = Function(V)
    solve(a == L, μ_fem, bc)
    
    # use fem_cov_assembler_post to obtain cov mat on sensor grid and function to compute vector
    # (c(x,y)) for y in Y
    C_fem_s, c_fem = fem_cov_assembler_post(J,k_f,Y,parallel=par,translation_inv=trans)
    
    # form B_fem_s by adding noise contribution
    C_fem_s += (ϵ**2)*np.eye(s)
    
    # assemble function to compute posterior mean and return
    def m_post_fem(x):
        return m_post(x,μ_fem,c_fem,v_dat,Y,C_fem_s)
    
    return m_post_fem
```

`m_post_fem_assembler` takes in several arguments which are explained below:

- `J`: controls the FE mesh size ($h=1/J^{2}$)
- `f_bar`: the mean function for the forcing $f$
- `k_f`: the covariance function for the forcing $f$
- `ϵ`: controls the amount of sensor noise
- `Y`: vector of sensor locations
- `v_dat`: vector of noisy sensor observations
- `par`: boolean argument passed to `fem_cov_assembler_post`'s argument `parallel` (defaults to `False`)
- `trans`: boolean argument passed to `fem_cov_assembler_post`'s argument `translation_inv` (defaults to `True`)

Let's quickly check that this function is working.


```
J = 20
f_bar = Constant(1.0)
m_post_fem = m_post_fem_assembler(J,f_bar,k_f,ϵ,Y,v_dat)
# compute posterior mean at a location x in D
x = np.array([0.3,0.1])
m_post_fem(x)
```

## statFEM posterior covariance

The form of the statFEM posterior covariance remains the same as given in <a href="/statFEM/oneDim.html#Posterior-covariance"><code>oneDim</code></a>. Thus, we require very similar code as to the 1-D case. We start by creating a function `c_post` which evaluates the posterior covariance at a given point.


```
#export
def c_post(x,y,c,Y,B):
    "This function evaluates the posterior covariance at $(x,y)$"
    
    # compute vectors c_x and c_y:
    c_x = np.array([c(x,y_i) for y_i in Y])
    c_y = np.array([c(y_i,y) for y_i in Y])
    
    # compute update term
    update = c_x @ np.linalg.solve(B,c_y)
    
    # return c_post
    return (c(x,y) - update)
```

`c_post` takes in several arguments which are explained below:

- `x`,`y`: points to evaluate the covariance at
- `c`: function which returns the prior covariance at any given pair $(x,y)$
- `Y`: vector of sensor locations
- `B`: the matrix $\epsilon^{2}I+C_{Y}$ to be inverted in order to obtain the posterior

To compare the statFEM covariance matrices for finer and finer FE mesh sizes we will require some more code. First we create a function `post_fem_cov_assembler` which helps us to quickly assemble the statFEM posterior covariance matrix as explained in <a href="/statFEM/oneDim.html#Difference-between-posterior-covariances"><code>oneDim</code></a>.


```
#export
def post_fem_cov_assembler(J,k_f,grid,Y,parallel,translation_inv):
    "Function which assembles the matrices $Σ_X$ ,$Σ_{XY}$, and $Σ_Y$ required for the statFEM posterior covariance."

    # set up mesh and function space
    mesh = UnitSquareMesh(J,J)
    V = FunctionSpace(mesh,'Lagrange',1)
    
    # set up grid
    x_grid = V.tabulate_dof_coordinates()
    
    # set up boundary condition
    def boundary(x, on_boundary):
        return on_boundary
    
    bc = DirichletBC(V, 0.0, boundary)
    
    # get the boundary and interior dofs
    bc_dofs = bc.get_boundary_values().keys()
    first, last = V.dofmap().ownership_range()
    all_dofs = range(last - first)
    interior_dofs = list(set(all_dofs) - set(bc_dofs))
    bc_dofs = list(set(bc_dofs))
    
    # set up the function κ
    κ = Constant(1.0)

    # get the mass and stiffness matrices
    u = TrialFunction(V)
    v = TestFunction(V)
    
    mass_form = u*v*dx
    a = inner(κ*grad(u),grad(v))*dx
    
    M = assemble(mass_form)
    A = assemble(a)
    M = as_backend_type(M).mat()
    A = as_backend_type(A).mat()
    M = csr_matrix(M.getValuesCSR()[::-1],shape=M.size)
    A = csr_matrix(A.getValuesCSR()[::-1],shape=A.size)
    
    # extract the submatrices corresponding to the interior dofs
    M = M[interior_dofs,:][:,interior_dofs]
    A = A[interior_dofs,:][:,interior_dofs]
    
    # get the forcing cov matrix on the interior nodes of the grid
    Σ_int = kernMat(k_f,x_grid[interior_dofs],parallel,translation_inv)

    # form the matrix Q in the defintion of the approximate FEM cov mat
    # Note: overwrite Σ_int for memory efficiency
    Σ_int = M @ Σ_int @ M.T

    Σ_int = spsolve(A,Σ_int)
    Σ_int = spsolve(A,Σ_int.T).T
    
    # ensure Σ_int is symmetric
    Σ_int = 0.5*(Σ_int + Σ_int.T)
    
    # get big phi matrix on the grid (only need the interior nodes)
    Phi_grid = BigPhiMat(J,grid)[interior_dofs,:]
    
    # get big phi matrix on the sensor grid (only need the interior nodes)
    Phi_Y = BigPhiMat(J,Y)[interior_dofs,:]
    
    # assemble the FEM cov mat on the sensor grid using Σ_int and Phi_Y
    Σ_Y = Phi_Y.T @ Σ_int @ Phi_Y
    
    # assemble the FEM cov mat on the grid using Σ_int and Phi_grid
    Σ_X = Phi_grid.T @ Σ_int @ Phi_grid
    
    # assemble cross term matrix (with ijth entry c(x_i,y_j))
    Σ_XY = Phi_grid.T @ Σ_int @ Phi_Y
    
    # return these sigma matrices
    return Σ_Y, Σ_X, Σ_XY
```

`post_fem_cov_assembler` takes in several arguments which are explained below:

- `J`: controls the FE mesh size ($h=1/J^2$)
- `k_f`: the covariance function for the forcing $f$
- `grid`: the fixed reference grid $\{x_{i}\}_{i=1}^{N}$ on which to assemble the posterior cov mat
- `Y`: vector of sensor locations.
- `parallel`: boolean argument indicating whether the computation of the forcing cov mat should be done in parallel
- `translation_inv`: boolean argument indicating whether the computation of the forcing cov mat should be computed assuming `k_f` is translation invariant or not

Finally, we create the function `c_post_fem_assembler` which assembles the statFEM posterior cov mat on the reference grid using the matrices `post_fem_cov_assembler` returns.


```
#export
def c_post_fem_assembler(J,k_f,grid,Y,ϵ,par,trans):
    "Function to assemble the statFEM posterior cov mat on a reference grid specified by grid."
    # use post_fem_cov_assembler to get the sigma matrices needed for posterior cov mat
    Σ_Y, Σ_X, Σ_XY = post_fem_cov_assembler(J,k_f,grid,Y,parallel=par,translation_inv=trans)
    
    # create the matrix B (store in Σ_Y for memory efficiency)
    s = len(Y) # number of sensor points
    Σ_Y += (ϵ**2)*np.eye(s)
    
    #form the posterior cov matrix (store in Σ_X for memory efficiency)
    Σ_X -= Σ_XY @ np.linalg.solve(Σ_Y,Σ_XY.T)
    
    return Σ_X
```


```
#hide
from nbdev.export import notebook2script; notebook2script()
```

    Converted 00_oneDim.ipynb.
    Converted 01_twoDim.ipynb.
    Converted index.ipynb.
    Converted oneDim_prior_results.ipynb.


e7a53fa05d7b47d06f89fd73dc64a851b09c016c
ipynb
Jupyter Notebook
01_twoDim.ipynb
YanniPapandreou/statFEM
189ddbb9c2f5a363d6e7e2f62a893cb3706e45bb
Apache-2.0
2022-02-04T09:26:33.000Z
2022-02-04T09:26:33.000Z
01_twoDim.ipynb
YanniPapandreou/statFEM
189ddbb9c2f5a363d6e7e2f62a893cb3706e45bb
Apache-2.0
01_twoDim.ipynb
YanniPapandreou/statFEM
189ddbb9c2f5a363d6e7e2f62a893cb3706e45bb
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES