---
title: "The SPDE model with transparent barriers"
author: "Elias T Krainski"
date: "October-2024"
output: 
  - rmarkdown::html_vignette
  - rmarkdown::pdf_document
bibliography: web/INLAspacetime.bib
vignette: >
  %\VignetteIndexEntry{The SPDE model with transparent barriers}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# The transparent barrier model

The transparent barrier model consider a domain
$\Omega$ which is partitioned into $k$
sub-domains, $\Omega_d$ for $d\in\{1,\ldots,k\}$, 
where $\cup_{d=1}^k\Omega_d=\Omega$.
The SPDE model assumes a common marginal variance
and a particular to each $\Omega_d$, $r_d$.

From @bakka2019barrier, the precision matrix is 
\[
\mathbf{Q} = \frac{1}{\sigma^2}\mathbf{R}_r\mathbf{\tilde{C}}^{-1}\mathbf{R}_r
\textrm{ for } 
 \mathbf{R}_r = \mathbf{C} +
 \frac{1}{8}\sum_{d=1}^kr_d^2\mathbf{G}_d ,
\;\;\;
\mathbf{\tilde{C}}_r = 
\frac{\pi}{2}\sum_{d=1}^kr_d^2\mathbf{\tilde{C}}_d
\]
where $\sigma^2$ is the marginal variance. 
The Finite Element Method - FEM matrices:
$\mathbf{C}$, defined as
\[ \mathbf{C}_{i,j} = \langle \psi_i, \psi_j \rangle = 
  \int_\Omega \psi_i(\mathbf{s}) \psi_j(\mathbf{s}) \partial \mathbf{s},\]
computed over the whole domain,
while $\mathbf{G}_d$ and $\mathbf{\tilde{C}}_d$ 
are defined as a pair of matrices for each subdomain
\[ (\mathbf{G}_d)_{i,j} = \langle 1_{\Omega_d} \nabla \psi_i, \nabla \psi_j \rangle = 
  \int_{\Omega_d} \nabla \psi_i(\mathbf{s}) \nabla \psi_j(\mathbf{s}) \partial \mathbf{s}\; \textrm{ and }\;
 (\mathbf{\tilde{C}}_d)_{i,i} = \langle 1_{\Omega_d} \psi_i, 1 \rangle = 
  \int_{\Omega_d} \psi_i(\mathbf{s}) \partial \mathbf{s} . \]

In the case when $r = r_1 = r_2 = \ldots = r_k$ we have
$\mathbf{R}_r = \mathbf{C}+\frac{r^2}{8}\mathbf{G}$ 
and $\mathbf{\tilde{C}}_r = \frac{\pi r^2}{2}\mathbf{\tilde{C}}$
giving 
\[ \mathbf{Q} = \frac{2}{\pi\sigma^2}(
\frac{1}{r^2}\mathbf{C}\mathbf{\tilde{C}}^{-1}\mathbf{C} + 
\frac{1}{8}\mathbf{C}\mathbf{\tilde{C}}^{-1}\mathbf{G} +
\frac{1}{8}\mathbf{G}\mathbf{\tilde{C}}^{-1}\mathbf{C} +
\frac{r^2}{64}\mathbf{G}\mathbf{\tilde{C}}^{-1}\mathbf{G} 
) \]
which coincides with the stationary case in @lindgren2015bayesian,
when using $\tilde{\mathbf{C}}$ in place of $\mathbf{C}$.

# Implementation

We assume $r_d = p_d r$ for known $p_1,\ldots,p_k$ constants. 
This gives
\[ \mathbf{\tilde{C}}_r = 
\frac{\pi r^2}{2}\sum_{d=1}^kp_d^2\mathbf{\tilde{C}}_d =
\frac{\pi r^2}{2} \mathbf{\tilde{C}}_{p_1,\ldots,p_k}
\textrm{ and }
\frac{1}{8}\sum_{d=1}^kr_d^2\mathbf{G}_d =
\frac{r^2}{8}\sum_{d=1}^kp_d^2\mathbf{\tilde{G}}_d = 
\frac{r^2}{8}\mathbf{\tilde{G}}_{p_1,\ldots,p_k}
\]
where $\mathbf{\tilde{C}}_{p_1,\ldots,p_k}$
and $\mathbf{\tilde{G}}_{p_1,\ldots,p_k}$ 
are pre-computed matrices.


# References