Numerical error

Go back to the About page. This link might be useful to keep track of the files created during the preprocessing.

Let us set some global options for all code chunks in this document.

# Set seed for reproducibility
set.seed(1982) 
# Set global options for all code chunks
knitr::opts_chunk$set(
  # Disable messages printed by R code chunks
  message = FALSE,    
  # Disable warnings printed by R code chunks
  warning = FALSE,    
  # Show R code within code chunks in output
  echo = TRUE,        
  # Include both R code and its results in output
  include = TRUE,     
  # Evaluate R code chunks
  eval = TRUE,       
  # Enable caching of R code chunks for faster rendering
  cache = FALSE,      
  # Align figures in the center of the output
  fig.align = "center",
  # Enable retina display for high-resolution figures
  retina = 2,
  # Show errors in the output instead of stopping rendering
  error = TRUE,
  # Do not collapse code and output into a single block
  collapse = FALSE
)
# Start the figure counter
fig_count <- 0
# Define the captioner function
captioner <- function(caption) {
  fig_count <<- fig_count + 1
  paste0("Figure ", fig_count, ": ", caption)
}

Import libraries

library(MetricGraph)
library(Matrix)
library(rSPDE)

library(dplyr)
library(tidyverse)
library(plotly)
library(ggplot2)
library(latex2exp)
library(ggtext)

library(grateful) # Cite all loaded packages
library(here) # here() starts from the home directory

Interval graph

---

Press the Show button below to reveal the code.

---

# library calls
library(MetricGraph)
library(Matrix)
library(rSPDE)

# function 1
gets_graph_interval <- function(n){
  edge <- rbind(c(0,0),c(1,0))
  edges = list(edge)
  graph <- metric_graph$new(edges = edges)
  graph$build_mesh(n = n)
  return(graph)
}

# matern covariance function. Same as in the package
matern.covariance <- function(h, kappa, nu, sigma) {
  if (nu == 1 / 2) {
    C <- sigma^2 * exp(-kappa * abs(h))
  } else {
    C <- (sigma^2 / (2^(nu - 1) * gamma(nu))) *
      ((kappa * abs(h))^nu) * besselK(kappa * abs(h), nu)
  }
  C[h == 0] <- sigma^2
  return(as.matrix(C))
}

# folded.matern.covariance.1d I edited
folded.matern.covariance.1d.local <- function(x, kappa, nu, sigma,
                                              L = 1, N = 10,
                                              boundary = c("neumann",
                                                           "dirichlet", "periodic")) {
  boundary <- tolower(boundary[1])
  if (!(boundary %in% c("neumann", "dirichlet", "periodic"))) {
    stop("The possible boundary conditions are 'neumann',
    'dirichlet' or 'periodic'!")
  }
  addi = t(outer(x, x, "+"))
  diff = t(outer(x, x, "-"))
  s1 <- sapply(-N:N, function(j) { 
    diff + 2 * j * L
  })
  s2 <- sapply(-N:N, function(j) {
    addi + 2 * j * L
  })
  if (boundary == "neumann") {
    C <- rowSums(matern.covariance(h = s1, kappa = kappa,
                                   nu = nu, sigma = sigma) +
                   matern.covariance(h = s2, kappa = kappa,
                                     nu = nu, sigma = sigma))
  } else if (boundary == "dirichlet") {
    C <- rowSums(matern.covariance(h = s1, kappa = kappa,
                                   nu = nu, sigma = sigma) -
                   matern.covariance(h = s2, kappa = kappa,
                                     nu = nu, sigma = sigma))
  } else {
    C <- rowSums(matern.covariance(h = s1,
                                   kappa = kappa, nu = nu, sigma = sigma))
  }
  return(matrix(C, nrow = length(x)))
}

# function 2
gets_true_cov_mat = function(graph, kappa, nu, sigma, N, boundary){
  h = graph$mesh$V[,1] 
  true_cov_mat = folded.matern.covariance.1d.local(x = h, kappa = kappa, nu = nu, sigma = sigma, N = N, boundary = boundary)
  return(true_cov_mat)
}


# parameters
n_vector = c(998, 198, 98, 8)
type_vector = c("covariance", "operator")
type_rational_approximation_vector = c("chebfun", "brasil", "chebfunLB")
rho_vector = c(0.1, 0.5, 1, 2)
m_vector = c(0,1,2,3,4,5)
nu_vector = c(seq(0.1,0.4,by=0.05), seq(0.41,0.59,by=0.01), seq(0.6,1.4,by=0.05), seq(1.41,1.59,by=0.01), seq(1.6,2.4,by=0.05), seq(2.41,2.49,by=0.01))
boundary = "neumann"
sigma = 1
N.folded = 10

Error = list()
for (s in c(1,2,3,4)) { # loop over n_vector
  Error[[as.character(n_vector[s])]] = list()
  n = n_vector[s]
  graph = gets_graph_interval(n = n)
  for (i in c(1)) { # loop over type_vector
    Error[[as.character(n_vector[s])]][[type_vector[i]]] = list()
    type = type_vector[i]
    for (j in c(1,2,3)) { # loop over type_rational_approximation_vector
      Error[[as.character(n_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]] = list()
      type_rational_approximation = type_rational_approximation_vector[j]
      for (k in c(1,2,3,4)) { # loop over rho_vector
        Error[[as.character(n_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]] = list()
        rho = rho_vector[k]
        for (l in 1:length(m_vector)) { # loop over m_vector
          Error[[as.character(n_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]] = matrix(NA, nrow = length(nu_vector), ncol = 2)
          colnames(Error[[as.character(n_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]]) = c("L_inf_error", "L_2_error")
          m = m_vector[l]
          for (r in 1:length(nu_vector)) { # loop over nu_vector
            nu = nu_vector[r]
            
            kappa = sqrt(8*nu)/rho
            tau = sqrt(gamma(nu) / (sigma^2 * kappa^(2*nu) * (4*pi)^(1/2) * gamma(nu + 1/2)))  #sigma = 1, d = 1
            alpha = nu + 1/2
            
            tryCatch({
              # getting true covariance
              true_cov_mat = gets_true_cov_mat(graph = graph,
                                               kappa = kappa,
                                               nu = nu,
                                               sigma = sigma,
                                               N = N.folded,
                                               boundary = boundary)
              
              # getting the approximate covariance matrix
              op = matern.operators(alpha = alpha, 
                                    kappa = kappa, 
                                    tau = tau,
                                    m = m, 
                                    graph = graph,
                                    type = type,
                                    type_rational_approximation = type_rational_approximation)
              
              appr_cov_mat = op$covariance_mesh()
            
            # computing the errors
            L_inf_error = max(abs(true_cov_mat - appr_cov_mat))
            L_2_error = sqrt(as.double(t(graph$mesh$weights)%*%(true_cov_mat - appr_cov_mat)^2%*%graph$mesh$weights))
            
            Error[[as.character(n_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][r,1] = L_inf_error
            Error[[as.character(n_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][r,2] = L_2_error
            }, error = function(err){
              warning(paste("Error occurred at iteration n=", n, ",type=", type, ",type_rational_approximation=", type_rational_approximation, ",rho=", rho, ",m=", m, ",nu=", nu, "Error:", conditionMessage(err)))
              print(paste("Error occurred at iteration n=", n, ",type=", type, ",type_rational_approximation=", type_rational_approximation, ",rho=", rho, ",m=", m, ",nu=", nu, "Error:", conditionMessage(err)))
              Error[[as.character(n_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][r,1] = NA
              Error[[as.character(n_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][r,2] = NA
            })
            print(paste("n=", n, ",type=", type, ",type_rational_approximation=", type_rational_approximation, ",rho=", rho, ",m=", m, ",nu=", nu))
          }
        }
      }
    }
  }
}

Error_interval_DEF = Error
save(Error_interval_DEF, file = here("data_files/Error_interval_DEF.RData"))

Plotting the errors

load(here("data_files/Error_interval_DEF.RData"))

n_vector = c("998", "198", "98", "8") 
type_rational_approximation_vector = c("chebfun", "brasil", "chebfunLB")
nu_max = 53
nu_vector = 0.5 + c(seq(0.1,0.4,by=0.05), seq(0.41,0.59,by=0.01), seq(0.6,1.4,by=0.05), seq(1.41,1.59,by=0.01), seq(1.6,2.4,by=0.05), seq(2.41,2.49,by=0.01))[1:nu_max]

n = "998"
t = "chebfun"


dat = rbind(
data.frame(nu = rep(nu_vector,5), 
           rho = rep(0.1, times = 5*nu_max),
           m = rep(1:5, each = nu_max),
           error = c(Error_interval_DEF[[n]][["covariance"]][[t]][["0.1"]][["1"]][1:nu_max,1],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["0.1"]][["2"]][1:nu_max,1],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["0.1"]][["3"]][1:nu_max,1],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["0.1"]][["4"]][1:nu_max,1],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["0.1"]][["5"]][1:nu_max,1]),
           L = rep("Linf", times = 5*nu_max)
),

data.frame(nu = rep(nu_vector,5), 
           rho = c(rep(0.5, times = 5*nu_max)),
           m = rep(1:5, each = nu_max),
           error = c(Error_interval_DEF[[n]][["covariance"]][[t]][["0.5"]][["1"]][1:nu_max,1],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["0.5"]][["2"]][1:nu_max,1],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["0.5"]][["3"]][1:nu_max,1],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["0.5"]][["4"]][1:nu_max,1],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["0.5"]][["5"]][1:nu_max,1]),
           L = rep("Linf", times = 5*nu_max)
),

data.frame(nu = rep(nu_vector,5), 
           rho = c(rep(1, times = 5*nu_max)),
           m = rep(1:5, each = nu_max),
           error = c(Error_interval_DEF[[n]][["covariance"]][[t]][["1"]][["1"]][1:nu_max,1],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["1"]][["2"]][1:nu_max,1],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["1"]][["3"]][1:nu_max,1],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["1"]][["4"]][1:nu_max,1],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["1"]][["5"]][1:nu_max,1]),
           L = rep("Linf", times = 5*nu_max)
),

data.frame(nu = rep(nu_vector,5), 
           rho = c(rep(2, times = 5*nu_max)),
           m = rep(1:5, each = nu_max),
           error = c(Error_interval_DEF[[n]][["covariance"]][[t]][["2"]][["1"]][1:nu_max,1],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["2"]][["2"]][1:nu_max,1],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["2"]][["3"]][1:nu_max,1],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["2"]][["4"]][1:nu_max,1],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["2"]][["5"]][1:nu_max,1]),
           L = rep("Linf", times = 5*nu_max)
),

## L2

data.frame(nu = rep(nu_vector,5), 
           rho = rep(0.1, times = 5*nu_max),
           m = rep(1:5, each = nu_max),
           error = c(Error_interval_DEF[[n]][["covariance"]][[t]][["0.1"]][["1"]][1:nu_max,2],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["0.1"]][["2"]][1:nu_max,2],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["0.1"]][["3"]][1:nu_max,2],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["0.1"]][["4"]][1:nu_max,2],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["0.1"]][["5"]][1:nu_max,2]),
           L = rep("L2", times = 5*nu_max)
),

data.frame(nu = rep(nu_vector,5), 
           rho = c(rep(0.5, times = 5*nu_max)),
           m = rep(1:5, each = nu_max),
           error = c(Error_interval_DEF[[n]][["covariance"]][[t]][["0.5"]][["1"]][1:nu_max,2],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["0.5"]][["2"]][1:nu_max,2],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["0.5"]][["3"]][1:nu_max,2],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["0.5"]][["4"]][1:nu_max,2],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["0.5"]][["5"]][1:nu_max,2]),
           L = rep("L2", times = 5*nu_max)
),

data.frame(nu = rep(nu_vector,5), 
           rho = c(rep(1, times = 5*nu_max)),
           m = rep(1:5, each = nu_max),
           error = c(Error_interval_DEF[[n]][["covariance"]][[t]][["1"]][["1"]][1:nu_max,2],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["1"]][["2"]][1:nu_max,2],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["1"]][["3"]][1:nu_max,2],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["1"]][["4"]][1:nu_max,2],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["1"]][["5"]][1:nu_max,2]),
           L = rep("L2", times = 5*nu_max)
),

data.frame(nu = rep(nu_vector,5), 
           rho = c(rep(2, times = 5*nu_max)),
           m = rep(1:5, each = nu_max),
           error = c(Error_interval_DEF[[n]][["covariance"]][[t]][["2"]][["1"]][1:nu_max,2],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["2"]][["2"]][1:nu_max,2],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["2"]][["3"]][1:nu_max,2],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["2"]][["4"]][1:nu_max,2],
                     Error_interval_DEF[[n]][["covariance"]][[t]][["2"]][["5"]][1:nu_max,2]),
           L = rep("L2", times = 5*nu_max)
)
)


dat$rho = factor(dat$rho, levels = c("0.1", "0.5", "1", "2"))
dat$L = factor(dat$L, levels = c("L2", "Linf"))
levels(dat$rho) = c("0.1" = latex2exp::TeX("$\\rho = 0.1$"), "0.5" = latex2exp::TeX("$\\rho = 0.5$"), "1" = latex2exp::TeX("$\\rho = 1$"), "2" = latex2exp::TeX("$\\rho = 2$"))
levels(dat$L) = c("L2" = latex2exp::TeX("$L_2(\\Gamma\\times\\Gamma)$ error"), "Linf" = latex2exp::TeX("$L_\\infty(\\Gamma\\times\\Gamma)$ error"))


p <- ggplot(dat, aes(nu, error, colour = as.factor(m))) + 
  facet_grid(L ~ rho, labeller = label_parsed) +
  geom_line() +
  scale_y_log10(n.breaks = 5) +
  scale_x_continuous(n.breaks = 10) +
  theme_bw() +
  theme(panel.spacing = unit(0.3, "cm"), text = element_text(family = "Palatino"),
        plot.title = element_text(size = 16),
        strip.text = element_text(size = 12),        # Panel titles
        axis.title = element_text(size = 14),        # Axis titles
        axis.text = element_text(size = 12),         # Axis text
        legend.title = element_text(size = 14),      # Legend title
        legend.text = element_text(size = 12)) +
  labs(x = bquote(alpha ~ "(smoothness parameter)"), y = "Covariance Error", color = "m") + ggtitle("Interval graph")

p

Figure 1: Covariance error for the interval graph.

ggsave(here("data_files/new_interval_upto1.5.png"), width = 12, height = 6, plot = p, dpi = 300)

Circle graph

---

Press the Show button below to reveal the code.

---

# library calls
library(MetricGraph)
library(Matrix)
library(rSPDE)

# function 1
gets_graph_circle <- function(n){
  r = 1/(pi)
  theta <- seq(from=-pi,to=pi,length.out = 10000)
  edge <- cbind(1+r+r*cos(theta),r*sin(theta))
  edges = list(edge)
  graph <- metric_graph$new(edges = edges)
  graph$set_manual_edge_lengths(edge_lengths = 2)
  graph$build_mesh(n = n)
  return(graph)
}

# matern covariance function. Same as in the package
matern.covariance <- function(h, kappa, nu, sigma) {
  if (nu == 1 / 2) {
    C <- sigma^2 * exp(-kappa * abs(h))
  } else {
    C <- (sigma^2 / (2^(nu - 1) * gamma(nu))) *
      ((kappa * abs(h))^nu) * besselK(kappa * abs(h), nu)
  }
  C[h == 0] <- sigma^2
  return(as.matrix(C))
}

# folded.matern.covariance.1d I edited
folded.matern.covariance.1d.local <- function(x, kappa, nu, sigma, L = 1, N = 10, boundary = c("neumann",
                                                                                               "dirichlet", "periodic")) {
  boundary <- tolower(boundary[1])
  if (!(boundary %in% c("neumann", "dirichlet", "periodic"))) {
    stop("The possible boundary conditions are 'neumann',
    'dirichlet' or 'periodic'!")
  }
  addi = t(outer(x, x, "+"))
  diff = t(outer(x, x, "-"))
  s1 <- sapply(-N:N, function(j) { # s1 is a matrix of size length(h)x(2N+1)
    diff + 2 * j * L
  })
  s2 <- sapply(-N:N, function(j) {
    addi + 2 * j * L
  })
  if (boundary == "neumann") {
    C <- rowSums(matern.covariance(h = s1, kappa = kappa,
                                   nu = nu, sigma = sigma) +
                   matern.covariance(h = s2, kappa = kappa,
                                     nu = nu, sigma = sigma))
  } else if (boundary == "dirichlet") {
    C <- rowSums(matern.covariance(h = s1, kappa = kappa,
                                   nu = nu, sigma = sigma) -
                   matern.covariance(h = s2, kappa = kappa,
                                     nu = nu, sigma = sigma))
  } else {
    C <- rowSums(matern.covariance(h = s1,
                                   kappa = kappa, nu = nu, sigma = sigma))
  }
  return(matrix(C, nrow = length(x)))
}

# function 2
gets_true_cov_mat = function(graph, kappa, nu, sigma, N, boundary){
  h = c(0,graph$get_edge_lengths()[1]*graph$mesh$PtE[,2])
  true_cov_mat = folded.matern.covariance.1d.local(x = h, kappa = kappa, nu = nu, sigma = sigma, N = N, boundary = boundary)
  return(true_cov_mat)
}


# parameters
n_vector = 2*c(998, 198, 98, 8)
type_vector = c("covariance", "operator")
type_rational_approximation_vector = c("chebfun", "brasil", "chebfunLB")
rho_vector = c(0.1, 0.5, 1, 2)
m_vector = c(0,1,2,3,4,5)
nu_vector = c(seq(0.1,0.4,by=0.05), seq(0.41,0.59,by=0.01), seq(0.6,1.4,by=0.05), seq(1.41,1.59,by=0.01), seq(1.6,2.4,by=0.05), seq(2.41,2.49,by=0.01))
boundary = "periodic"
sigma = 1
N.folded = 10

Error = list()
for (s in c(1,2,3,4)) { # loop over n_vector
  Error[[as.character(n_vector[s])]] = list()
  n = n_vector[s]
  graph = gets_graph_circle(n = n)
  for (i in c(1)) { # loop over type_vector
    Error[[as.character(n_vector[s])]][[type_vector[i]]] = list()
    type = type_vector[i]
    for (j in c(1,2,3)) { # loop over type_rational_approximation_vector
      Error[[as.character(n_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]] = list()
      type_rational_approximation = type_rational_approximation_vector[j]
      for (k in c(1,2,3,4)) { # loop over rho_vector
        Error[[as.character(n_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]] = list()
        rho = rho_vector[k]
        for (l in 1:length(m_vector)) { # loop over m_vector
          Error[[as.character(n_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]] = matrix(NA, nrow = length(nu_vector), ncol = 2)
          colnames(Error[[as.character(n_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]]) = c("L_inf_error", "L_2_error")
          m = m_vector[l]
          for (r in 1:length(nu_vector)) { # loop over nu_vector
            nu = nu_vector[r]
            
            kappa = sqrt(8*nu)/rho
            tau = sqrt(gamma(nu) / (sigma^2 * kappa^(2*nu) * (4*pi)^(1/2) * gamma(nu + 1/2)))  #sigma = 1, d = 1
            alpha = nu + 1/2
            
            tryCatch({
              # getting true covariance
              true_cov_mat = gets_true_cov_mat(graph = graph,
                                               kappa = kappa,
                                               nu = nu,
                                               sigma = sigma,
                                               N = N.folded,
                                               boundary = boundary)
              
              # getting the approximate covariance matrix
              op = matern.operators(alpha = alpha, 
                                    kappa = kappa, 
                                    tau = tau,
                                    m = m, 
                                    graph = graph,
                                    type = type,
                                    type_rational_approximation = type_rational_approximation)
              
              appr_cov_mat = op$covariance_mesh()
            
            # computing the errors
            L_inf_error = max(abs(true_cov_mat - appr_cov_mat))
            L_2_error = sqrt(as.double(t(graph$mesh$weights)%*%(true_cov_mat - appr_cov_mat)^2%*%graph$mesh$weights))
            
            Error[[as.character(n_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][r,1] = L_inf_error
            Error[[as.character(n_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][r,2] = L_2_error
            }, error = function(err){
              warning(paste("Error occurred at iteration n=", n, ",type=", type, ",type_rational_approximation=", type_rational_approximation, ",rho=", rho, ",m=", m, ",nu=", nu, "Error:", conditionMessage(err)))
              print(paste("Error occurred at iteration n=", n, ",type=", type, ",type_rational_approximation=", type_rational_approximation, ",rho=", rho, ",m=", m, ",nu=", nu, "Error:", conditionMessage(err)))
              Error[[as.character(n_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][r,1] = NA
              Error[[as.character(n_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][r,2] = NA
            })
            print(paste("n=", n, ",type=", type, ",type_rational_approximation=", type_rational_approximation, ",rho=", rho, ",m=", m, ",nu=", nu))
          }
        }
      }
    }
  }
}
Error_circle_DEF = Error
save(Error_circle_DEF, file = here("data_files/Error_circle_DEF.RData"))

Plotting the errors

load(here("data_files/Error_circle_DEF.RData"))

n_vector = c("1996", "396", "196", "16") 
type_rational_approximation_vector = c("chebfun", "brasil", "chebfunLB")
nu_max = 53
nu_vector = 0.5 + c(seq(0.1,0.4,by=0.05), seq(0.41,0.59,by=0.01), seq(0.6,1.4,by=0.05), seq(1.41,1.59,by=0.01), seq(1.6,2.4,by=0.05), seq(2.41,2.49,by=0.01))[1:nu_max]

n = "1996"
t = "chebfun"
  

dat = rbind(
data.frame(nu = rep(nu_vector,5), 
           rho = rep(0.1, times = 5*nu_max),
           m = rep(1:5, each = nu_max),
           error = c(Error_circle_DEF[[n]][["covariance"]][[t]][["0.1"]][["1"]][1:nu_max,1],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["0.1"]][["2"]][1:nu_max,1],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["0.1"]][["3"]][1:nu_max,1],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["0.1"]][["4"]][1:nu_max,1],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["0.1"]][["5"]][1:nu_max,1]),
           L = rep("Linf", times = 5*nu_max)
),

data.frame(nu = rep(nu_vector,5), 
           rho = c(rep(0.5, times = 5*nu_max)),
           m = rep(1:5, each = nu_max),
           error = c(Error_circle_DEF[[n]][["covariance"]][[t]][["0.5"]][["1"]][1:nu_max,1],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["0.5"]][["2"]][1:nu_max,1],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["0.5"]][["3"]][1:nu_max,1],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["0.5"]][["4"]][1:nu_max,1],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["0.5"]][["5"]][1:nu_max,1]),
           L = rep("Linf", times = 5*nu_max)
),

data.frame(nu = rep(nu_vector,5), 
           rho = c(rep(1, times = 5*nu_max)),
           m = rep(1:5, each = nu_max),
           error = c(Error_circle_DEF[[n]][["covariance"]][[t]][["1"]][["1"]][1:nu_max,1],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["1"]][["2"]][1:nu_max,1],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["1"]][["3"]][1:nu_max,1],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["1"]][["4"]][1:nu_max,1],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["1"]][["5"]][1:nu_max,1]),
           L = rep("Linf", times = 5*nu_max)
),

data.frame(nu = rep(nu_vector,5), 
           rho = c(rep(2, times = 5*nu_max)),
           m = rep(1:5, each = nu_max),
           error = c(Error_circle_DEF[[n]][["covariance"]][[t]][["2"]][["1"]][1:nu_max,1],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["2"]][["2"]][1:nu_max,1],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["2"]][["3"]][1:nu_max,1],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["2"]][["4"]][1:nu_max,1],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["2"]][["5"]][1:nu_max,1]),
           L = rep("Linf", times = 5*nu_max)
),

## L2

data.frame(nu = rep(nu_vector,5), 
           rho = rep(0.1, times = 5*nu_max),
           m = rep(1:5, each = nu_max),
           error = c(Error_circle_DEF[[n]][["covariance"]][[t]][["0.1"]][["1"]][1:nu_max,2],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["0.1"]][["2"]][1:nu_max,2],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["0.1"]][["3"]][1:nu_max,2],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["0.1"]][["4"]][1:nu_max,2],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["0.1"]][["5"]][1:nu_max,2]),
           L = rep("L2", times = 5*nu_max)
),

data.frame(nu = rep(nu_vector,5), 
           rho = c(rep(0.5, times = 5*nu_max)),
           m = rep(1:5, each = nu_max),
           error = c(Error_circle_DEF[[n]][["covariance"]][[t]][["0.5"]][["1"]][1:nu_max,2],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["0.5"]][["2"]][1:nu_max,2],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["0.5"]][["3"]][1:nu_max,2],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["0.5"]][["4"]][1:nu_max,2],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["0.5"]][["5"]][1:nu_max,2]),
           L = rep("L2", times = 5*nu_max)
),

data.frame(nu = rep(nu_vector,5), 
           rho = c(rep(1, times = 5*nu_max)),
           m = rep(1:5, each = nu_max),
           error = c(Error_circle_DEF[[n]][["covariance"]][[t]][["1"]][["1"]][1:nu_max,2],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["1"]][["2"]][1:nu_max,2],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["1"]][["3"]][1:nu_max,2],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["1"]][["4"]][1:nu_max,2],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["1"]][["5"]][1:nu_max,2]),
           L = rep("L2", times = 5*nu_max)
),

data.frame(nu = rep(nu_vector,5), 
           rho = c(rep(2, times = 5*nu_max)),
           m = rep(1:5, each = nu_max),
           error = c(Error_circle_DEF[[n]][["covariance"]][[t]][["2"]][["1"]][1:nu_max,2],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["2"]][["2"]][1:nu_max,2],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["2"]][["3"]][1:nu_max,2],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["2"]][["4"]][1:nu_max,2],
                     Error_circle_DEF[[n]][["covariance"]][[t]][["2"]][["5"]][1:nu_max,2]),
           L = rep("L2", times = 5*nu_max)
)
)


dat$rho = factor(dat$rho, levels = c("0.1", "0.5", "1", "2"))
dat$L = factor(dat$L, levels = c("L2", "Linf"))
levels(dat$rho) = c("0.1" = latex2exp::TeX("$\\rho = 0.1$"), "0.5" = latex2exp::TeX("$\\rho = 0.5$"), "1" = latex2exp::TeX("$\\rho = 1$"), "2" = latex2exp::TeX("$\\rho = 2$"))
levels(dat$L) = c("L2" = latex2exp::TeX("$L_2(\\Gamma\\times\\Gamma)$ error"), "Linf" = latex2exp::TeX("$L_\\infty(\\Gamma\\times\\Gamma)$ error"))

p <- ggplot(dat, aes(nu, error, colour = as.factor(m))) + 
  facet_grid(L ~ rho, labeller = label_parsed) +
  geom_line() +
  scale_y_log10(n.breaks = 5) +
  scale_x_continuous(n.breaks = 10) +
  theme_bw() +
  theme(panel.spacing = unit(0.3, "cm"), text = element_text(family = "Palatino"),
        plot.title = element_text(size = 16),
        strip.text = element_text(size = 12),        # Panel titles
        axis.title = element_text(size = 14),        # Axis titles
        axis.text = element_text(size = 12),         # Axis text
        legend.title = element_text(size = 14),      # Legend title
        legend.text = element_text(size = 12)) +
  labs(x = bquote(alpha ~ "(smoothness parameter)"), y = "Covariance Error", color = "m") + ggtitle("Circle graph")

p

Figure 2: Covariance error for the circle graph.

ggsave(here("data_files/new_circle_upto1.5.png"), width = 12, height = 6, plot = p, dpi = 300)

Tadpole graph

rho = 0.1

---

Press the Show button below to reveal the code.

---

# library calls
library(MetricGraph)
library(Matrix)
library(rSPDE)

# function 1
tadpole.eig <- function(k,graph){
  x1 <- c(0,graph$get_edge_lengths()[1]*graph$mesh$PtE[graph$mesh$PtE[,1]==1,2]) 
  x2 <- c(0,graph$get_edge_lengths()[2]*graph$mesh$PtE[graph$mesh$PtE[,1]==2,2]) 
  
  if(k==0){ 
    f.e1 <- rep(1,length(x1)) 
    f.e2 <- rep(1,length(x2)) 
    f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1]) 
    f = list(phi=f1/sqrt(3)) 
    
  } else {
    f.e1 <- -2*sin(pi*k*1/2)*cos(pi*k*x1/2) 
    f.e2 <- sin(pi*k*x2/2)                  
    
    f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1]) 
    
    if((k %% 2)==1){ 
      f = list(phi=f1/sqrt(3)) 
    } else { 
      f.e1 <- (-1)^{k/2}*cos(pi*k*x1/2)
      f.e2 <- cos(pi*k*x2/2)
      f2 = c(f.e1[1],f.e2[1],f.e1[-1],f.e2[-1]) 
      f <- list(phi=f1,psi=f2/sqrt(3/2))
    }
  }
  
  return(f)
}

# function 2
gets_graph <- function(h){
  edge1 <- rbind(c(0,0),c(1,0))
  theta <- seq(from=-pi,to=pi,length.out = 10000)
  edge2 <- cbind(1+1/pi+cos(theta)/pi,sin(theta)/pi)
  edges = list(edge1, edge2)
  graph <- metric_graph$new(edges = edges)
  graph$set_manual_edge_lengths(edge_lengths = c(1,2))
  graph$build_mesh(h=h)
  return(graph)
}

#function 3
gets_true_cov_mat <- function(graph, kappa, tau, alpha, n.overkill){
  Sigma.kl <- matrix(0,nrow = dim(graph$mesh$V)[1],ncol = dim(graph$mesh$V)[1])
  for(i in 0:n.overkill){
    phi <- tadpole.eig(i,graph)$phi
    Sigma.kl <- Sigma.kl + (1/(kappa^2 + (i*pi/2)^2)^(alpha))*phi%*%t(phi)
    if(i>0 && (i %% 2)==0){ 
      psi <- tadpole.eig(i,graph)$psi
      Sigma.kl <- Sigma.kl + (1/(kappa^2 + (i*pi/2)^2)^(alpha))*psi%*%t(psi)
    }
    
  }
  Sigma.kl <- Sigma.kl/tau^2
  return(Sigma.kl)
}

# parameters
h_vector = c(0.001, 0.005, 0.01, 0.1)
type_vector = c("covariance", "operator")
type_rational_approximation_vector = c("chebfun", "brasil", "chebfunLB")
rho_vector = c(0.1, 0.5, 1, 2)
m_vector = c(0,1,2,3,4,5)
nu_vector = c(seq(0.1,0.4,by=0.05), seq(0.41,0.59,by=0.01), seq(0.6,1.4,by=0.05), seq(1.41,1.59,by=0.01), seq(1.6,2.4,by=0.05), seq(2.41,2.49,by=0.01))
n.overkill = 1000
sigma = 1

Error = list()
for (s in c(1,2,3,4)) { # loop over h_vector
  Error[[as.character(h_vector[s])]] = list()
  h = h_vector[s]
  graph = gets_graph(h = h)
  for (i in c(1)) { # loop over type_vector
    Error[[as.character(h_vector[s])]][[type_vector[i]]] = list()
    type = type_vector[i]
    for (j in c(1,2,3)) { # loop over type_rational_approximation_vector
      Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]] = list()
      type_rational_approximation = type_rational_approximation_vector[j]
      for (k in c(1)) { # loop over rho_vector
        Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]] = list()
        rho = rho_vector[k]
        for (l in 1:length(m_vector)) { # loop over m_vector
          Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]] = matrix(NA, nrow = length(nu_vector), ncol = 2)
          colnames(Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]]) = c("L_inf_error", "L_2_error")
          m = m_vector[l]
          for (n in 1:length(nu_vector)) { # loop over nu_vector
            nu = nu_vector[n]
            
            kappa = sqrt(8*nu)/rho
            tau = sqrt(gamma(nu) / (sigma^2 * kappa^(2*nu) * (4*pi)^(1/2) * gamma(nu + 1/2)))  #sigma = 1, d = 1
            alpha = nu + 1/2
            
            tryCatch({
            # getting true covariance
            true_cov_mat = gets_true_cov_mat(graph = graph,
                                             kappa = kappa,
                                             tau = tau,
                                             alpha = alpha,
                                             n.overkill = n.overkill)
            
            # getting the approximate covariance matrix
            op = matern.operators(alpha = alpha, 
                                  kappa = kappa, 
                                  tau = tau,
                                  m = m, 
                                  graph = graph,
                                  type = type,
                                  type_rational_approximation = type_rational_approximation)
            appr_cov_mat = op$covariance_mesh()
            
            # computing the errors
            L_inf_error = max(abs(true_cov_mat - appr_cov_mat))
            L_2_error = sqrt(as.double(t(graph$mesh$weights)%*%(true_cov_mat - appr_cov_mat)^2%*%graph$mesh$weights))
            
            Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][n,1] = L_inf_error
            Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][n,2] = L_2_error
            }, error = function(err){
              warning(paste("Error occurred at iteration h=", h, ",type=", type, ",type_rational_approximation=", type_rational_approximation, ",rho=", rho, ",m=", m, ",nu=", nu, "Error:", conditionMessage(err)))
              print(paste("Error occurred at iteration h=", h, ",type=", type, ",type_rational_approximation=", type_rational_approximation, ",rho=", rho, ",m=", m, ",nu=", nu, "Error:", conditionMessage(err)))
              Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][n,1] = NA
              Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][n,2] = NA
            })
            print(paste("h=", h, ",type=", type, ",type_rational_approximation=", type_rational_approximation, ",rho=", rho, ",m=", m, ",nu=", nu))
          }
        }
      }
    }
  }
}

Error_tadpole_DEF_rho0.1 = Error
save(Error_tadpole_DEF_rho0.1, file = here("data_files/Error_tadpole_DEF_rho0.1.RData"))

rho = 0.5

---

Press the Show button below to reveal the code.

---

# library calls
library(MetricGraph)
library(Matrix)
library(rSPDE)

# function 1
tadpole.eig <- function(k,graph){
  x1 <- c(0,graph$get_edge_lengths()[1]*graph$mesh$PtE[graph$mesh$PtE[,1]==1,2]) 
  x2 <- c(0,graph$get_edge_lengths()[2]*graph$mesh$PtE[graph$mesh$PtE[,1]==2,2]) 
  
  if(k==0){ 
    f.e1 <- rep(1,length(x1)) 
    f.e2 <- rep(1,length(x2)) 
    f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1]) 
    f = list(phi=f1/sqrt(3)) 
    
  } else {
    f.e1 <- -2*sin(pi*k*1/2)*cos(pi*k*x1/2) 
    f.e2 <- sin(pi*k*x2/2)                  
    
    f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1]) 
    
    if((k %% 2)==1){ 
      f = list(phi=f1/sqrt(3)) 
    } else { 
      f.e1 <- (-1)^{k/2}*cos(pi*k*x1/2)
      f.e2 <- cos(pi*k*x2/2)
      f2 = c(f.e1[1],f.e2[1],f.e1[-1],f.e2[-1]) 
      f <- list(phi=f1,psi=f2/sqrt(3/2))
    }
  }
  
  return(f)
}

# function 2
gets_graph <- function(h){
  edge1 <- rbind(c(0,0),c(1,0))
  theta <- seq(from=-pi,to=pi,length.out = 10000)
  edge2 <- cbind(1+1/pi+cos(theta)/pi,sin(theta)/pi)
  edges = list(edge1, edge2)
  graph <- metric_graph$new(edges = edges)
  graph$set_manual_edge_lengths(edge_lengths = c(1,2))
  graph$build_mesh(h=h)
  return(graph)
}

#function 3
gets_true_cov_mat <- function(graph, kappa, tau, alpha, n.overkill){
  Sigma.kl <- matrix(0,nrow = dim(graph$mesh$V)[1],ncol = dim(graph$mesh$V)[1])
  for(i in 0:n.overkill){
    phi <- tadpole.eig(i,graph)$phi
    Sigma.kl <- Sigma.kl + (1/(kappa^2 + (i*pi/2)^2)^(alpha))*phi%*%t(phi)
    if(i>0 && (i %% 2)==0){ 
      psi <- tadpole.eig(i,graph)$psi
      Sigma.kl <- Sigma.kl + (1/(kappa^2 + (i*pi/2)^2)^(alpha))*psi%*%t(psi)
    }
    
  }
  Sigma.kl <- Sigma.kl/tau^2
  return(Sigma.kl)
}

# parameters
h_vector = c(0.001, 0.005, 0.01, 0.1)
type_vector = c("covariance", "operator")
type_rational_approximation_vector = c("chebfun", "brasil", "chebfunLB")
rho_vector = c(0.1, 0.5, 1, 2)
m_vector = c(0,1,2,3,4,5)
nu_vector = c(seq(0.1,0.4,by=0.05), seq(0.41,0.59,by=0.01), seq(0.6,1.4,by=0.05), seq(1.41,1.59,by=0.01), seq(1.6,2.4,by=0.05), seq(2.41,2.49,by=0.01))
n.overkill = 1000
sigma = 1

Error = list()
for (s in c(1,2,3,4)) { # loop over h_vector
  Error[[as.character(h_vector[s])]] = list()
  h = h_vector[s]
  graph = gets_graph(h = h)
  for (i in c(1)) { # loop over type_vector
    Error[[as.character(h_vector[s])]][[type_vector[i]]] = list()
    type = type_vector[i]
    for (j in c(1,2,3)) { # loop over type_rational_approximation_vector
      Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]] = list()
      type_rational_approximation = type_rational_approximation_vector[j]
      for (k in c(2)) { # loop over rho_vector
        Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]] = list()
        rho = rho_vector[k]
        for (l in 1:length(m_vector)) { # loop over m_vector
          Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]] = matrix(NA, nrow = length(nu_vector), ncol = 2)
          colnames(Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]]) = c("L_inf_error", "L_2_error")
          m = m_vector[l]
          for (n in 1:length(nu_vector)) { # loop over nu_vector
            nu = nu_vector[n]
            
            kappa = sqrt(8*nu)/rho
            tau = sqrt(gamma(nu) / (sigma^2 * kappa^(2*nu) * (4*pi)^(1/2) * gamma(nu + 1/2)))  #sigma = 1, d = 1
            alpha = nu + 1/2
            
            tryCatch({
            # getting true covariance
            true_cov_mat = gets_true_cov_mat(graph = graph,
                                             kappa = kappa,
                                             tau = tau,
                                             alpha = alpha,
                                             n.overkill = n.overkill)
            
            # getting the approximate covariance matrix
            op = matern.operators(alpha = alpha, 
                                  kappa = kappa, 
                                  tau = tau,
                                  m = m, 
                                  graph = graph,
                                  type = type,
                                  type_rational_approximation = type_rational_approximation)
            appr_cov_mat = op$covariance_mesh()
            
            # computing the errors
            L_inf_error = max(abs(true_cov_mat - appr_cov_mat))
            L_2_error = sqrt(as.double(t(graph$mesh$weights)%*%(true_cov_mat - appr_cov_mat)^2%*%graph$mesh$weights))
            
            Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][n,1] = L_inf_error
            Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][n,2] = L_2_error
            }, error = function(err){
              warning(paste("Error occurred at iteration h=", h, ",type=", type, ",type_rational_approximation=", type_rational_approximation, ",rho=", rho, ",m=", m, ",nu=", nu, "Error:", conditionMessage(err)))
              print(paste("Error occurred at iteration h=", h, ",type=", type, ",type_rational_approximation=", type_rational_approximation, ",rho=", rho, ",m=", m, ",nu=", nu, "Error:", conditionMessage(err)))
              Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][n,1] = NA
              Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][n,2] = NA
            })
            print(paste("h=", h, ",type=", type, ",type_rational_approximation=", type_rational_approximation, ",rho=", rho, ",m=", m, ",nu=", nu))
          }
        }
      }
    }
  }
}

Error_tadpole_DEF_rho0.5 = Error
save(Error_tadpole_DEF_rho0.5, file = here("data_files/Error_tadpole_DEF_rho0.5.RData"))

rho = 1

---

Press the Show button below to reveal the code.

---

# library calls
library(MetricGraph)
library(Matrix)
library(rSPDE)

# function 1
tadpole.eig <- function(k,graph){
  x1 <- c(0,graph$get_edge_lengths()[1]*graph$mesh$PtE[graph$mesh$PtE[,1]==1,2]) 
  x2 <- c(0,graph$get_edge_lengths()[2]*graph$mesh$PtE[graph$mesh$PtE[,1]==2,2]) 
  
  if(k==0){ 
    f.e1 <- rep(1,length(x1)) 
    f.e2 <- rep(1,length(x2)) 
    f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1]) 
    f = list(phi=f1/sqrt(3)) 
    
  } else {
    f.e1 <- -2*sin(pi*k*1/2)*cos(pi*k*x1/2) 
    f.e2 <- sin(pi*k*x2/2)                  
    
    f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1]) 
    
    if((k %% 2)==1){ 
      f = list(phi=f1/sqrt(3)) 
    } else { 
      f.e1 <- (-1)^{k/2}*cos(pi*k*x1/2)
      f.e2 <- cos(pi*k*x2/2)
      f2 = c(f.e1[1],f.e2[1],f.e1[-1],f.e2[-1]) 
      f <- list(phi=f1,psi=f2/sqrt(3/2))
    }
  }
  
  return(f)
}

# function 2
gets_graph <- function(h){
  edge1 <- rbind(c(0,0),c(1,0))
  theta <- seq(from=-pi,to=pi,length.out = 10000)
  edge2 <- cbind(1+1/pi+cos(theta)/pi,sin(theta)/pi)
  edges = list(edge1, edge2)
  graph <- metric_graph$new(edges = edges)
  graph$set_manual_edge_lengths(edge_lengths = c(1,2))
  graph$build_mesh(h=h)
  return(graph)
}

#function 3
gets_true_cov_mat <- function(graph, kappa, tau, alpha, n.overkill){
  Sigma.kl <- matrix(0,nrow = dim(graph$mesh$V)[1],ncol = dim(graph$mesh$V)[1])
  for(i in 0:n.overkill){
    phi <- tadpole.eig(i,graph)$phi
    Sigma.kl <- Sigma.kl + (1/(kappa^2 + (i*pi/2)^2)^(alpha))*phi%*%t(phi)
    if(i>0 && (i %% 2)==0){ 
      psi <- tadpole.eig(i,graph)$psi
      Sigma.kl <- Sigma.kl + (1/(kappa^2 + (i*pi/2)^2)^(alpha))*psi%*%t(psi)
    }
    
  }
  Sigma.kl <- Sigma.kl/tau^2
  return(Sigma.kl)
}

# parameters
h_vector = c(0.001, 0.005, 0.01, 0.1)
type_vector = c("covariance", "operator")
type_rational_approximation_vector = c("chebfun", "brasil", "chebfunLB")
rho_vector = c(0.1, 0.5, 1, 2)
m_vector = c(0,1,2,3,4,5)
nu_vector = c(seq(0.1,0.4,by=0.05), seq(0.41,0.59,by=0.01), seq(0.6,1.4,by=0.05), seq(1.41,1.59,by=0.01), seq(1.6,2.4,by=0.05), seq(2.41,2.49,by=0.01)) 
n.overkill = 1000
sigma = 1

Error = list()
for (s in c(1,2,3,4)) { # loop over h_vector
  Error[[as.character(h_vector[s])]] = list()
  h = h_vector[s]
  graph = gets_graph(h = h)
  for (i in c(1)) { # loop over type_vector
    Error[[as.character(h_vector[s])]][[type_vector[i]]] = list()
    type = type_vector[i]
    for (j in c(1,2,3)) { # loop over type_rational_approximation_vector
      Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]] = list()
      type_rational_approximation = type_rational_approximation_vector[j]
      for (k in c(3)) { # loop over rho_vector
        Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]] = list()
        rho = rho_vector[k]
        for (l in 1:length(m_vector)) { # loop over m_vector
          Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]] = matrix(NA, nrow = length(nu_vector), ncol = 2)
          colnames(Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]]) = c("L_inf_error", "L_2_error")
          m = m_vector[l]
          for (n in 1:length(nu_vector)) { # loop over nu_vector
            nu = nu_vector[n]
            
            kappa = sqrt(8*nu)/rho
            tau = sqrt(gamma(nu) / (sigma^2 * kappa^(2*nu) * (4*pi)^(1/2) * gamma(nu + 1/2)))  #sigma = 1, d = 1
            alpha = nu + 1/2
            
            tryCatch({
            # getting true covariance
            true_cov_mat = gets_true_cov_mat(graph = graph,
                                             kappa = kappa,
                                             tau = tau,
                                             alpha = alpha,
                                             n.overkill = n.overkill)
            
            # getting the approximate covariance matrix
            op = matern.operators(alpha = alpha, 
                                  kappa = kappa, 
                                  tau = tau,
                                  m = m, 
                                  graph = graph,
                                  type = type,
                                  type_rational_approximation = type_rational_approximation)
            appr_cov_mat = op$covariance_mesh()
            
            # computing the errors
            L_inf_error = max(abs(true_cov_mat - appr_cov_mat))
            L_2_error = sqrt(as.double(t(graph$mesh$weights)%*%(true_cov_mat - appr_cov_mat)^2%*%graph$mesh$weights))
            
            Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][n,1] = L_inf_error
            Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][n,2] = L_2_error
            }, error = function(err){
              warning(paste("Error occurred at iteration h=", h, ",type=", type, ",type_rational_approximation=", type_rational_approximation, ",rho=", rho, ",m=", m, ",nu=", nu, "Error:", conditionMessage(err)))
              print(paste("Error occurred at iteration h=", h, ",type=", type, ",type_rational_approximation=", type_rational_approximation, ",rho=", rho, ",m=", m, ",nu=", nu, "Error:", conditionMessage(err)))
              Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][n,1] = NA
              Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][n,2] = NA
            })
            print(paste("h=", h, ",type=", type, ",type_rational_approximation=", type_rational_approximation, ",rho=", rho, ",m=", m, ",nu=", nu))
          }
        }
      }
    }
  }
}

Error_tadpole_DEF_rho1 = Error
save(Error_tadpole_DEF_rho1, file = here("data_files/Error_tadpole_DEF_rho1.RData"))

rho = 2

---

Press the Show button below to reveal the code.

---

# library calls
library(MetricGraph)
library(Matrix)
library(rSPDE)

# function 1
tadpole.eig <- function(k,graph){
  x1 <- c(0,graph$get_edge_lengths()[1]*graph$mesh$PtE[graph$mesh$PtE[,1]==1,2]) 
  x2 <- c(0,graph$get_edge_lengths()[2]*graph$mesh$PtE[graph$mesh$PtE[,1]==2,2]) 
  
  if(k==0){ 
    f.e1 <- rep(1,length(x1)) 
    f.e2 <- rep(1,length(x2)) 
    f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1]) 
    f = list(phi=f1/sqrt(3)) 
    
  } else {
    f.e1 <- -2*sin(pi*k*1/2)*cos(pi*k*x1/2) 
    f.e2 <- sin(pi*k*x2/2)                  
    
    f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1]) 
    
    if((k %% 2)==1){ 
      f = list(phi=f1/sqrt(3)) 
    } else { 
      f.e1 <- (-1)^{k/2}*cos(pi*k*x1/2)
      f.e2 <- cos(pi*k*x2/2)
      f2 = c(f.e1[1],f.e2[1],f.e1[-1],f.e2[-1]) 
      f <- list(phi=f1,psi=f2/sqrt(3/2))
    }
  }
  
  return(f)
}

# function 2
gets_graph <- function(h){
  edge1 <- rbind(c(0,0),c(1,0))
  theta <- seq(from=-pi,to=pi,length.out = 10000)
  edge2 <- cbind(1+1/pi+cos(theta)/pi,sin(theta)/pi)
  edges = list(edge1, edge2)
  graph <- metric_graph$new(edges = edges)
  graph$set_manual_edge_lengths(edge_lengths = c(1,2))
  graph$build_mesh(h=h)
  return(graph)
}

#function 3
gets_true_cov_mat <- function(graph, kappa, tau, alpha, n.overkill){
  Sigma.kl <- matrix(0,nrow = dim(graph$mesh$V)[1],ncol = dim(graph$mesh$V)[1])
  for(i in 0:n.overkill){
    phi <- tadpole.eig(i,graph)$phi
    Sigma.kl <- Sigma.kl + (1/(kappa^2 + (i*pi/2)^2)^(alpha))*phi%*%t(phi)
    if(i>0 && (i %% 2)==0){ 
      psi <- tadpole.eig(i,graph)$psi
      Sigma.kl <- Sigma.kl + (1/(kappa^2 + (i*pi/2)^2)^(alpha))*psi%*%t(psi)
    }
    
  }
  Sigma.kl <- Sigma.kl/tau^2
  return(Sigma.kl)
}

# parameters
h_vector = c(0.001, 0.005, 0.01, 0.1)
type_vector = c("covariance", "operator")
type_rational_approximation_vector = c("chebfun", "brasil", "chebfunLB")
rho_vector = c(0.1, 0.5, 1, 2)
m_vector = c(0,1,2,3,4,5)
nu_vector = c(seq(0.1,0.4,by=0.05), seq(0.41,0.59,by=0.01), seq(0.6,1.4,by=0.05), seq(1.41,1.59,by=0.01), seq(1.6,2.4,by=0.05), seq(2.41,2.49,by=0.01)) 
n.overkill = 1000
sigma = 1

Error = list()
for (s in c(1,2,3,4)) { # loop over h_vector
  Error[[as.character(h_vector[s])]] = list()
  h = h_vector[s]
  graph = gets_graph(h = h)
  for (i in c(1)) { # loop over type_vector
    Error[[as.character(h_vector[s])]][[type_vector[i]]] = list()
    type = type_vector[i]
    for (j in c(1,2,3)) { # loop over type_rational_approximation_vector
      Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]] = list()
      type_rational_approximation = type_rational_approximation_vector[j]
      for (k in c(4)) { # loop over rho_vector
        Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]] = list()
        rho = rho_vector[k]
        for (l in 1:length(m_vector)) { # loop over m_vector
          Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]] = matrix(NA, nrow = length(nu_vector), ncol = 2)
          colnames(Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]]) = c("L_inf_error", "L_2_error")
          m = m_vector[l]
          for (n in 1:length(nu_vector)) { # loop over nu_vector
            nu = nu_vector[n]
            
            kappa = sqrt(8*nu)/rho
            tau = sqrt(gamma(nu) / (sigma^2 * kappa^(2*nu) * (4*pi)^(1/2) * gamma(nu + 1/2)))  #sigma = 1, d = 1
            alpha = nu + 1/2
            
            tryCatch({
            # getting true covariance
            true_cov_mat = gets_true_cov_mat(graph = graph,
                                             kappa = kappa,
                                             tau = tau,
                                             alpha = alpha,
                                             n.overkill = n.overkill)
            
            # getting the approximate covariance matrix
            op = matern.operators(alpha = alpha, 
                                  kappa = kappa, 
                                  tau = tau,
                                  m = m, 
                                  graph = graph,
                                  type = type,
                                  type_rational_approximation = type_rational_approximation)
            appr_cov_mat = op$covariance_mesh()
            
            # computing the errors
            L_inf_error = max(abs(true_cov_mat - appr_cov_mat))
            L_2_error = sqrt(as.double(t(graph$mesh$weights)%*%(true_cov_mat - appr_cov_mat)^2%*%graph$mesh$weights))
            
            Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][n,1] = L_inf_error
            Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][n,2] = L_2_error
            }, error = function(err){
              warning(paste("Error occurred at iteration h=", h, ",type=", type, ",type_rational_approximation=", type_rational_approximation, ",rho=", rho, ",m=", m, ",nu=", nu, "Error:", conditionMessage(err)))
              print(paste("Error occurred at iteration h=", h, ",type=", type, ",type_rational_approximation=", type_rational_approximation, ",rho=", rho, ",m=", m, ",nu=", nu, "Error:", conditionMessage(err)))
              Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][n,1] = NA
              Error[[as.character(h_vector[s])]][[type_vector[i]]][[type_rational_approximation_vector[j]]][[as.character(rho_vector[k])]][[as.character(m_vector[l])]][n,2] = NA
            })
            print(paste("h=", h, ",type=", type, ",type_rational_approximation=", type_rational_approximation, ",rho=", rho, ",m=", m, ",nu=", nu))
          }
        }
      }
    }
  }
}

Error_tadpole_DEF_rho2 = Error
save(Error_tadpole_DEF_rho2, file = here("data_files/Error_tadpole_DEF_rho2.RData"))

Plotting the errors

load(here("data_files/Error_tadpole_DEF_rho0.1.RData"))
load(here("data_files/Error_tadpole_DEF_rho0.5.RData"))
load(here("data_files/Error_tadpole_DEF_rho1.RData"))
load(here("data_files/Error_tadpole_DEF_rho2.RData"))

h_vector = c("0.001", "0.005", "0.01", "0.1")
type_rational_approximation_vector = c("chebfun", "brasil", "chebfunLB")
nu_max = 53
nu_vector = 0.5 + c(seq(0.1,0.4,by=0.05), seq(0.41,0.59,by=0.01), seq(0.6,1.4,by=0.05), seq(1.41,1.59,by=0.01), seq(1.6,2.4,by=0.05), seq(2.41,2.49,by=0.01))[1:nu_max]

n = "0.001"
t = "chebfun"


dat = rbind(
data.frame(nu = rep(nu_vector,5), 
           rho = rep(0.1, times = 5*nu_max),
           m = rep(1:5, each = nu_max),
           error = c(Error_tadpole_DEF_rho0.1[[n]][["covariance"]][[t]][["0.1"]][["1"]][1:nu_max,1],
                     Error_tadpole_DEF_rho0.1[[n]][["covariance"]][[t]][["0.1"]][["2"]][1:nu_max,1],
                     Error_tadpole_DEF_rho0.1[[n]][["covariance"]][[t]][["0.1"]][["3"]][1:nu_max,1],
                     Error_tadpole_DEF_rho0.1[[n]][["covariance"]][[t]][["0.1"]][["4"]][1:nu_max,1],
                     Error_tadpole_DEF_rho0.1[[n]][["covariance"]][[t]][["0.1"]][["5"]][1:nu_max,1]),
           L = rep("Linf", times = 5*nu_max)
),

data.frame(nu = rep(nu_vector,5), 
           rho = c(rep(0.5, times = 5*nu_max)),
           m = rep(1:5, each = nu_max),
           error = c(Error_tadpole_DEF_rho0.5[[n]][["covariance"]][[t]][["0.5"]][["1"]][1:nu_max,1],
                     Error_tadpole_DEF_rho0.5[[n]][["covariance"]][[t]][["0.5"]][["2"]][1:nu_max,1],
                     Error_tadpole_DEF_rho0.5[[n]][["covariance"]][[t]][["0.5"]][["3"]][1:nu_max,1],
                     Error_tadpole_DEF_rho0.5[[n]][["covariance"]][[t]][["0.5"]][["4"]][1:nu_max,1],
                     Error_tadpole_DEF_rho0.5[[n]][["covariance"]][[t]][["0.5"]][["5"]][1:nu_max,1]),
           L = rep("Linf", times = 5*nu_max)
),

data.frame(nu = rep(nu_vector,5), 
           rho = c(rep(1, times = 5*nu_max)),
           m = rep(1:5, each = nu_max),
           error = c(Error_tadpole_DEF_rho1[[n]][["covariance"]][[t]][["1"]][["1"]][1:nu_max,1],
                     Error_tadpole_DEF_rho1[[n]][["covariance"]][[t]][["1"]][["2"]][1:nu_max,1],
                     Error_tadpole_DEF_rho1[[n]][["covariance"]][[t]][["1"]][["3"]][1:nu_max,1],
                     Error_tadpole_DEF_rho1[[n]][["covariance"]][[t]][["1"]][["4"]][1:nu_max,1],
                     Error_tadpole_DEF_rho1[[n]][["covariance"]][[t]][["1"]][["5"]][1:nu_max,1]),
           L = rep("Linf", times = 5*nu_max)
),

data.frame(nu = rep(nu_vector,5), 
           rho = c(rep(2, times = 5*nu_max)),
           m = rep(1:5, each = nu_max),
           error = c(Error_tadpole_DEF_rho2[[n]][["covariance"]][[t]][["2"]][["1"]][1:nu_max,1],
                     Error_tadpole_DEF_rho2[[n]][["covariance"]][[t]][["2"]][["2"]][1:nu_max,1],
                     Error_tadpole_DEF_rho2[[n]][["covariance"]][[t]][["2"]][["3"]][1:nu_max,1],
                     Error_tadpole_DEF_rho2[[n]][["covariance"]][[t]][["2"]][["4"]][1:nu_max,1],
                     Error_tadpole_DEF_rho2[[n]][["covariance"]][[t]][["2"]][["5"]][1:nu_max,1]),
           L = rep("Linf", times = 5*nu_max)
),

## L2

data.frame(nu = rep(nu_vector,5), 
           rho = rep(0.1, times = 5*nu_max),
           m = rep(1:5, each = nu_max),
           error = c(Error_tadpole_DEF_rho0.1[[n]][["covariance"]][[t]][["0.1"]][["1"]][1:nu_max,2],
                     Error_tadpole_DEF_rho0.1[[n]][["covariance"]][[t]][["0.1"]][["2"]][1:nu_max,2],
                     Error_tadpole_DEF_rho0.1[[n]][["covariance"]][[t]][["0.1"]][["3"]][1:nu_max,2],
                     Error_tadpole_DEF_rho0.1[[n]][["covariance"]][[t]][["0.1"]][["4"]][1:nu_max,2],
                     Error_tadpole_DEF_rho0.1[[n]][["covariance"]][[t]][["0.1"]][["5"]][1:nu_max,2]),
           L = rep("L2", times = 5*nu_max)
),

data.frame(nu = rep(nu_vector,5), 
           rho = c(rep(0.5, times = 5*nu_max)),
           m = rep(1:5, each = nu_max),
           error = c(Error_tadpole_DEF_rho0.5[[n]][["covariance"]][[t]][["0.5"]][["1"]][1:nu_max,2],
                     Error_tadpole_DEF_rho0.5[[n]][["covariance"]][[t]][["0.5"]][["2"]][1:nu_max,2],
                     Error_tadpole_DEF_rho0.5[[n]][["covariance"]][[t]][["0.5"]][["3"]][1:nu_max,2],
                     Error_tadpole_DEF_rho0.5[[n]][["covariance"]][[t]][["0.5"]][["4"]][1:nu_max,2],
                     Error_tadpole_DEF_rho0.5[[n]][["covariance"]][[t]][["0.5"]][["5"]][1:nu_max,2]),
           L = rep("L2", times = 5*nu_max)
),

data.frame(nu = rep(nu_vector,5), 
           rho = c(rep(1, times = 5*nu_max)),
           m = rep(1:5, each = nu_max),
           error = c(Error_tadpole_DEF_rho1[[n]][["covariance"]][[t]][["1"]][["1"]][1:nu_max,2],
                     Error_tadpole_DEF_rho1[[n]][["covariance"]][[t]][["1"]][["2"]][1:nu_max,2],
                     Error_tadpole_DEF_rho1[[n]][["covariance"]][[t]][["1"]][["3"]][1:nu_max,2],
                     Error_tadpole_DEF_rho1[[n]][["covariance"]][[t]][["1"]][["4"]][1:nu_max,2],
                     Error_tadpole_DEF_rho1[[n]][["covariance"]][[t]][["1"]][["5"]][1:nu_max,2]),
           L = rep("L2", times = 5*nu_max)
),

data.frame(nu = rep(nu_vector,5), 
           rho = c(rep(2, times = 5*nu_max)),
           m = rep(1:5, each = nu_max),
           error = c(Error_tadpole_DEF_rho2[[n]][["covariance"]][[t]][["2"]][["1"]][1:nu_max,2],
                     Error_tadpole_DEF_rho2[[n]][["covariance"]][[t]][["2"]][["2"]][1:nu_max,2],
                     Error_tadpole_DEF_rho2[[n]][["covariance"]][[t]][["2"]][["3"]][1:nu_max,2],
                     Error_tadpole_DEF_rho2[[n]][["covariance"]][[t]][["2"]][["4"]][1:nu_max,2],
                     Error_tadpole_DEF_rho2[[n]][["covariance"]][[t]][["2"]][["5"]][1:nu_max,2]),
           L = rep("L2", times = 5*nu_max)
)
)


dat$rho = factor(dat$rho, levels = c("0.1", "0.5", "1", "2"))
dat$L = factor(dat$L, levels = c("L2", "Linf"))
levels(dat$rho) = c("0.1" = latex2exp::TeX("$\\rho = 0.1$"), "0.5" = latex2exp::TeX("$\\rho = 0.5$"), "1" = latex2exp::TeX("$\\rho = 1$"), "2" = latex2exp::TeX("$\\rho = 2$"))
levels(dat$L) = c("L2" = latex2exp::TeX("$L_2(\\Gamma\\times\\Gamma)$ error"), "Linf" = latex2exp::TeX("$L_\\infty(\\Gamma\\times\\Gamma)$ error"))

p <- ggplot(dat, aes(nu, error, colour = as.factor(m))) + 
  facet_grid(L ~ rho, labeller = label_parsed) +
  geom_line() +
  scale_y_log10(n.breaks = 5) +
  scale_x_continuous(n.breaks = 10) +
  theme_bw() +
  theme(panel.spacing = unit(0.3, "cm"), text = element_text(family = "Palatino"),
        plot.title = element_text(size = 16),
        strip.text = element_text(size = 12),        # Panel titles
        axis.title = element_text(size = 14),        # Axis titles
        axis.text = element_text(size = 12),         # Axis text
        legend.title = element_text(size = 14),      # Legend title
        legend.text = element_text(size = 12)) +
  labs(x = bquote(alpha ~ "(smoothness parameter)"), y = "Covariance Error", color = "m") + ggtitle("Tadpole graph")

p

Figure 3: Covariance error for the tadpole graph.

ggsave(here("data_files/new_tadpole_upto1.5.png"), width = 12, height = 6, plot = p, dpi = 300)

References

cite_packages(output = "paragraph", out.dir = ".")

We used R version 4.4.1 (R Core Team 2024a) and the following R packages: cowplot v. 1.1.3 (Wilke 2024), ggmap v. 4.0.0.900 (Kahle and Wickham 2013), ggpubr v. 0.6.0 (Kassambara 2023), ggtext v. 0.1.2 (Wilke and Wiernik 2022), grid v. 4.4.1 (R Core Team 2024b), here v. 1.0.1 (Müller 2020), htmltools v. 0.5.8.1 (Cheng et al. 2024), INLA v. 24.12.11 (Rue, Martino, and Chopin 2009; Lindgren, Rue, and Lindström 2011; Martins et al. 2013; Lindgren and Rue 2015; De Coninck et al. 2016; Rue et al. 2017; Verbosio et al. 2017; Bakka et al. 2018; Kourounis, Fuchs, and Schenk 2018), inlabru v. 2.12.0.9002 (Yuan et al. 2017; Bachl et al. 2019), knitr v. 1.48 (Xie 2014, 2015, 2024), latex2exp v. 0.9.6 (Meschiari 2022), Matrix v. 1.6.5 (Bates, Maechler, and Jagan 2024), MetricGraph v. 1.4.0.9000 (Bolin, Simas, and Wallin 2023b, 2023a, 2023c, 2024; Bolin et al. 2024), OpenStreetMap v. 0.4.0 (Fellows and JMapViewer library by Jan Peter Stotz 2023), osmdata v. 0.2.5 (Mark Padgham et al. 2017), patchwork v. 1.2.0 (Pedersen 2024), plotly v. 4.10.4 (Sievert 2020), plotrix v. 3.8.4 (J 2006), reshape2 v. 1.4.4 (Wickham 2007), rmarkdown v. 2.28 (Xie, Allaire, and Grolemund 2018; Xie, Dervieux, and Riederer 2020; Allaire et al. 2024), rSPDE v. 2.4.0.9000 (Bolin and Kirchner 2020; Bolin and Simas 2023; Bolin, Simas, and Xiong 2024), scales v. 1.3.0 (Wickham, Pedersen, and Seidel 2023), sf v. 1.0.19 (E. Pebesma 2018; E. Pebesma and Bivand 2023), sp v. 2.1.4 (E. J. Pebesma and Bivand 2005; Bivand, Pebesma, and Gomez-Rubio 2013), tidyverse v. 2.0.0 (Wickham et al. 2019), viridis v. 0.6.4 (Garnier et al. 2023), xaringanExtra v. 0.8.0 (Aden-Buie and Warkentin 2024).

Aden-Buie, Garrick, and Matthew T. Warkentin. 2024. xaringanExtra: Extras and Extensions for “xaringan” Slides. https://CRAN.R-project.org/package=xaringanExtra.

Allaire, JJ, Yihui Xie, Christophe Dervieux, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, et al. 2024. rmarkdown: Dynamic Documents for r. https://github.com/rstudio/rmarkdown.

Bachl, Fabian E., Finn Lindgren, David L. Borchers, and Janine B. Illian. 2019. “inlabru: An R Package for Bayesian Spatial Modelling from Ecological Survey Data.” Methods in Ecology and Evolution 10: 760–66. https://doi.org/10.1111/2041-210X.13168.

Bakka, Haakon, Håvard Rue, Geir-Arne Fuglstad, Andrea I. Riebler, David Bolin, Janine Illian, Elias Krainski, Daniel P. Simpson, and Finn K. Lindgren. 2018. “Spatial Modelling with INLA: A Review.” WIRES (Invited Extended Review) xx (Feb): xx–. http://arxiv.org/abs/1802.06350.

Bates, Douglas, Martin Maechler, and Mikael Jagan. 2024. Matrix: Sparse and Dense Matrix Classes and Methods. https://CRAN.R-project.org/package=Matrix.

Bivand, Roger S., Edzer Pebesma, and Virgilio Gomez-Rubio. 2013. Applied Spatial Data Analysis with R, Second Edition. Springer, NY. https://asdar-book.org/.

Bolin, David, and Kristin Kirchner. 2020. “The Rational SPDE Approach for Gaussian Random Fields with General Smoothness.” Journal of Computational and Graphical Statistics 29 (2): 274–85. https://doi.org/10.1080/10618600.2019.1665537.

Bolin, David, Mihály Kovács, Vivek Kumar, and Alexandre B. Simas. 2024. “Regularity and Numerical Approximation of Fractional Elliptic Differential Equations on Compact Metric Graphs.” Mathematics of Computation 93 (349): 2439–72. https://doi.org/10.1090/mcom/3929.

Bolin, David, and Alexandre B. Simas. 2023. rSPDE: Rational Approximations of Fractional Stochastic Partial Differential Equations. https://CRAN.R-project.org/package=rSPDE.

Bolin, David, Alexandre B. Simas, and Jonas Wallin. 2023a. “Markov Properties of Gaussian Random Fields on Compact Metric Graphs.” arXiv Preprint arXiv:2304.03190. https://doi.org/10.48550/arXiv.2304.03190.

———. 2023b. MetricGraph: Random Fields on Metric Graphs. https://CRAN.R-project.org/package=MetricGraph.

———. 2023c. “Statistical Inference for Gaussian Whittle-Matérn Fields on Metric Graphs.” arXiv Preprint arXiv:2304.10372. https://doi.org/10.48550/arXiv.2304.10372.

———. 2024. “Gaussian Whittle-Matérn Fields on Metric Graphs.” Bernoulli 30 (2): 1611–39. https://doi.org/10.3150/23-BEJ1647.

Bolin, David, Alexandre B. Simas, and Zhen Xiong. 2024. “Covariance-Based Rational Approximations of Fractional SPDEs for Computationally Efficient Bayesian Inference.” Journal of Computational and Graphical Statistics 33 (1): 64–74. https://doi.org/10.1080/10618600.2023.2231051.

Cheng, Joe, Carson Sievert, Barret Schloerke, Winston Chang, Yihui Xie, and Jeff Allen. 2024. htmltools: Tools for HTML. https://CRAN.R-project.org/package=htmltools.

De Coninck, Arne, Bernard De Baets, Drosos Kourounis, Fabio Verbosio, Olaf Schenk, Steven Maenhout, and Jan Fostier. 2016. “Needles: Toward Large-Scale Genomic Prediction with Marker-by-Environment Interaction.” Genetics 203 (1): 543–55. https://doi.org/10.1534/genetics.115.179887.

Fellows, Ian, and using the JMapViewer library by Jan Peter Stotz. 2023. OpenStreetMap: Access to Open Street Map Raster Images. https://CRAN.R-project.org/package=OpenStreetMap.

Garnier, Simon, Ross, Noam, Rudis, Robert, Camargo, et al. 2023. viridis(Lite) - Colorblind-Friendly Color Maps for r. https://doi.org/10.5281/zenodo.4679423.

J, Lemon. 2006. “Plotrix: A Package in the Red Light District of r.” R-News 6 (4): 8–12.

Kahle, David, and Hadley Wickham. 2013. “ggmap: Spatial Visualization with Ggplot2.” The R Journal 5 (1): 144–61. https://journal.r-project.org/archive/2013-1/kahle-wickham.pdf.

Kassambara, Alboukadel. 2023. ggpubr: “ggplot2” Based Publication Ready Plots. https://CRAN.R-project.org/package=ggpubr.

Kourounis, D., A. Fuchs, and O. Schenk. 2018. “Towards the Next Generation of Multiperiod Optimal Power Flow Solvers.” IEEE Transactions on Power Systems PP (99): 1–10. https://doi.org/10.1109/TPWRS.2017.2789187.

Lindgren, Finn, and Håvard Rue. 2015. “Bayesian Spatial Modelling with R-INLA.” Journal of Statistical Software 63 (19): 1–25. http://www.jstatsoft.org/v63/i19/.

Lindgren, Finn, Håvard Rue, and Johan Lindström. 2011. “An Explicit Link Between Gaussian Fields and Gaussian Markov Random Fields: The Stochastic Partial Differential Equation Approach (with Discussion).” Journal of the Royal Statistical Society B 73 (4): 423–98.

Mark Padgham, Bob Rudis, Robin Lovelace, and Maëlle Salmon. 2017. “Osmdata.” Journal of Open Source Software 2 (14): 305. https://doi.org/10.21105/joss.00305.

Martins, Thiago G., Daniel Simpson, Finn Lindgren, and Håvard Rue. 2013. “Bayesian Computing with INLA: New Features.” Computational Statistics and Data Analysis 67: 68–83.

Meschiari, Stefano. 2022. Latex2exp: Use LaTeX Expressions in Plots. https://CRAN.R-project.org/package=latex2exp.

Müller, Kirill. 2020. here: A Simpler Way to Find Your Files. https://CRAN.R-project.org/package=here.

Pebesma, Edzer. 2018. “Simple Features for R: Standardized Support for Spatial Vector Data.” The R Journal 10 (1): 439–46. https://doi.org/10.32614/RJ-2018-009.

Pebesma, Edzer J., and Roger Bivand. 2005. “Classes and Methods for Spatial Data in R.” R News 5 (2): 9–13. https://CRAN.R-project.org/doc/Rnews/.

Pebesma, Edzer, and Roger Bivand. 2023. Spatial Data Science: With applications in R. Chapman and Hall/CRC. https://doi.org/10.1201/9780429459016.

Pedersen, Thomas Lin. 2024. patchwork: The Composer of Plots. https://CRAN.R-project.org/package=patchwork.

R Core Team. 2024a. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

———. 2024b. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Rue, Håvard, Sara Martino, and Nicholas Chopin. 2009. “Approximate Bayesian Inference for Latent Gaussian Models Using Integrated Nested Laplace Approximations (with Discussion).” Journal of the Royal Statistical Society B 71: 319–92.

Rue, Håvard, Andrea I. Riebler, Sigrunn H. Sørbye, Janine B. Illian, Daniel P. Simpson, and Finn K. Lindgren. 2017. “Bayesian Computing with INLA: A Review.” Annual Reviews of Statistics and Its Applications 4 (March): 395–421. http://arxiv.org/abs/1604.00860.

Sievert, Carson. 2020. Interactive Web-Based Data Visualization with r, Plotly, and Shiny. Chapman; Hall/CRC. https://plotly-r.com.

Verbosio, Fabio, Arne De Coninck, Drosos Kourounis, and Olaf Schenk. 2017. “Enhancing the Scalability of Selected Inversion Factorization Algorithms in Genomic Prediction.” Journal of Computational Science 22 (Supplement C): 99–108. https://doi.org/10.1016/j.jocs.2017.08.013.

Wickham, Hadley. 2007. “Reshaping Data with the reshape Package.” Journal of Statistical Software 21 (12): 1–20. http://www.jstatsoft.org/v21/i12/.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.

Wickham, Hadley, Thomas Lin Pedersen, and Dana Seidel. 2023. scales: Scale Functions for Visualization. https://CRAN.R-project.org/package=scales.

Wilke, Claus O. 2024. cowplot: Streamlined Plot Theme and Plot Annotations for “ggplot2”. https://CRAN.R-project.org/package=cowplot.

Wilke, Claus O., and Brenton M. Wiernik. 2022. ggtext: Improved Text Rendering Support for “ggplot2”. https://CRAN.R-project.org/package=ggtext.

Xie, Yihui. 2014. “knitr: A Comprehensive Tool for Reproducible Research in R.” In Implementing Reproducible Computational Research, edited by Victoria Stodden, Friedrich Leisch, and Roger D. Peng. Chapman; Hall/CRC.

———. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.

———. 2024. knitr: A General-Purpose Package for Dynamic Report Generation in r. https://yihui.org/knitr/.

Xie, Yihui, J. J. Allaire, and Garrett Grolemund. 2018. R Markdown: The Definitive Guide. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown.

Xie, Yihui, Christophe Dervieux, and Emily Riederer. 2020. R Markdown Cookbook. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown-cookbook.

Yuan, Yuan, Bachl, Fabian E., Lindgren, Finn, Borchers, et al. 2017. “Point Process Models for Spatio-Temporal Distance Sampling Data from a Large-Scale Survey of Blue Whales.” Ann. Appl. Stat. 11 (4): 2270–97. https://doi.org/10.1214/17-AOAS1078.