R/fractional.operators.R

In rSPDE: Rational Approximations of Fractional Stochastic Partial Differential Equations

## fractional.operators.R
##
##   Copyright (C) 2018, 2019, David Bolin
##
##   This program is free software: you can redistribute it and/or modify
##   it under the terms of the GNU General Public License as published by
##   the Free Software Foundation, either version 3 of the License, or
##   (at your option) any later version.
##
##   This program is distributed in the hope that it will be useful,
##   but WITHOUT ANY WARRANTY; without even the implied warranty of
##   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##   GNU General Public License for more details.
##
##   You should have received a copy of the GNU General Public License
##   along with this program.  If not, see <http://www.gnu.org/licenses/>.

#' Rational approximations of fractional operators
#'
#' `fractional.operators` is used for computing an approximation,
#' which can be used for inference and simulation, of the fractional SPDE
#' \deqn{L^\beta (\tau u(s)) = W.}
#' Here \eqn{L} is a differential operator, \eqn{\beta>0} is
#' the fractional power, \eqn{\tau} is a positive scalar or vector that
#' scales the variance of the solution \eqn{u}, and \eqn{W} is white noise.
#'
#' @param L A finite element discretization of the operator \eqn{L}{L}.
#' @param beta The positive fractional power.
#' @param C The mass matrix of the finite element discretization.
#' @param scale.factor A constant \eqn{c} is a lower bound for the the smallest
#' eigenvalue of the non-discretized operator \eqn{L}{L}.
#' @param m The order of the rational approximation, which needs to be a
#' positive integer. The default value is 1.
#' Higer values gives a more accurate approximation, which are more
#' computationally expensive to use for inference. Currently, the largest value
#' of m that is implemented is 4.
#' @param tau The constant or vector that scales the variance of the solution.
#' The default value is 1.
#'
#' @return `fractional.operators` returns an object of class "rSPDEobj".
#' This object contains the following quantities:
#' \item{Pl}{The operator \eqn{P_l}.}
#' \item{Pr}{The operator \eqn{P_r}.}
#' \item{C}{The mass lumped mass matrix.}
#' \item{Ci}{The inverse of `C`.}
#' \item{m}{The order of the rational approximation.}
#' \item{beta}{The fractional power.}
#' \item{type}{String indicating the type of approximation.}
#' \item{Q}{The matrix `t(Pl) %*% solve(C,Pl)`.}
#' \item{type}{String indicating the type of approximation.}
#' \item{Pl.factors}{List with elements that can be used to assemble \eqn{P_l}.}
#' \item{Pr.factors}{List with elements that can be used to assemble \eqn{P_r}.}
#' @export
#' @seealso [matern.operators()], [spde.matern.operators()],
#' [matern.operators()]
#' @details The approximation is based on a rational approximation of
#' the fractional operator, resulting in an
#' approximate model on the form \deqn{P_l u(s) = P_r W,}
#' where \eqn{P_j = p_j(L)} are non-fractional operators defined in terms of
#' polynomials \eqn{p_j} for \eqn{j=l,r}. The order of \eqn{p_r} is given by
#' `m` and the order of \eqn{p_l} is \eqn{m + m_\beta}
#' where \eqn{m_\beta} is the integer part of \eqn{\beta} if \eqn{\beta>1} and
#' \eqn{m_\beta = 1} otherwise.
#'
#' The discrete approximation can be written as \eqn{u = P_r x} where
#' \eqn{x \sim N(0,Q^{-1})}{x ~ N(0,Q^{-1})} and \eqn{Q = P_l^T C^{-1} P_l}.
#' Note that the matrices \eqn{P_r} and \eqn{Q} may be be ill-conditioned
#' for \eqn{m>1}. In this case, the methods in [operator.operations()]
#' should be used for operations involving the matrices, since these methods
#' are more numerically stable.
#'
#' @examples
#' # Compute rational approximation of a Gaussian process with a
#' # Matern covariance function on R
#' kappa <- 10
#' sigma <- 1
#' nu <- 0.8
#'
#' # create mass and stiffness matrices for a FEM discretization
#' x <- seq(from = 0, to = 1, length.out = 101)
#' fem <- rSPDE.fem1d(x)
#'
#' # compute rational approximation of covariance function at 0.5
#' tau <- sqrt(gamma(nu) / (sigma^2 * kappa^(2 * nu) *
#'   (4 * pi)^(1 / 2) * gamma(nu + 1 / 2)))
#' op <- fractional.operators(
#'   L = fem$G + kappa^2 * fem$C, beta = (nu + 1 / 2) / 2,
#'   C = fem$C, scale.factor = kappa^2, tau = tau
#' )
#'
#' v <- t(rSPDE.A1d(x, 0.5))
#' c.approx <- Sigma.mult(op, v)
#'
#' # plot the result and compare with the true Matern covariance
#' plot(x, matern.covariance(abs(x - 0.5), kappa, nu, sigma),
#'   type = "l", ylab = "C(h)",
#'   xlab = "h", main = "Matern covariance and rational approximations"
#' )
#' lines(x, c.approx, col = 2)
#'
fractional.operators <- function(L,
                                 beta,
                                 C,
                                 scale.factor,
                                 m = 1,
                                 tau = 1) {
  if (min(tau) < 0) {
    stop("tau should be positive")
  }
  if ((m %% 1) != 0 || m < 0) {
    stop("m must be a positive integer")
  }
  if (scale.factor <= 0) {
    stop("the scaling factor must be positive")
  }

  C <- Matrix::Diagonal(dim(C)[1], rowSums(C))
  Ci <- Matrix::Diagonal(dim(C)[1], 1 / rowSums(C))
  I <- Matrix::Diagonal(dim(C)[1])
  L <- L / scale.factor
  CiL <- Ci %*% L
  if (beta %% 1 == 0) { # not fractional
    Pr <- I
    Pl <- L
    if (beta > 1) {
      for (i in 1:(beta - 1)) {
        Pl <- Pl %*% CiL
      }
    }
    Pl.roots <- Pl.factors <- NULL
    Pl.k <- beta
    Pl.scaling <- scale.factor^beta
    Pr.roots <- Pr.factors <- NULL
  } else {
    roots <- get.roots(m, beta)
    Pl.roots <- roots$rb
    Pr.roots <- roots$rc
    m_beta <- max(1, floor(beta))
    Pr.factors <- Pl.factors <- list()
    # construct Pl
    Pl <- I - CiL * roots$rb[1]
    Pl.factors[[1]] <- I - CiL * roots$rb[1]
    if (length(roots$rb) > 1) {
      for (i in 2:length(roots$rb)) {
        Pl <- Pl %*% (I - CiL * roots$rb[i])
        Pl.factors[[i]] <- I - CiL * roots$rb[i]
      }
    }
    Lp <- C
    if (m_beta > 1) {
      for (i in 1:(m_beta - 1)) {
        Lp <- Lp %*% CiL
      }
    }
    Pl.k <- m_beta - 1
    Pl.scaling <- scale.factor^beta / roots$factor
    Pl <- Lp %*% Pl

    # construct Pr
    Pr <- I - CiL * roots$rc[1]
    Pr.factors[[1]] <- I - CiL * roots$rc[1]
    if (length(roots$rc) > 1) {
      for (i in 2:length(roots$rc)) {
        Pr <- Pr %*% (I - CiL * roots$rc[i])
        Pr.factors[[i]] <- I - CiL * roots$rc[i]
      }
    }
  }
  Pl <- Pl * Pl.scaling
  Phi <- Matrix::Diagonal(dim(C)[1], 1 / tau)
  output <- list(
    Q = t(Pl) %*% Ci %*% Pl,
    Pl = Pl,
    Pr = Phi %*% Pr,
    Ci = Ci,
    C = C,
    CiL = CiL,
    L = L,
    Pl.factors = list(
      scaling = Pl.scaling,
      roots = Pl.roots,
      factor = Pl.factors,
      k = Pl.k
    ),
    Pr.factors = list(
      roots = Pr.roots,
      factor = Pr.factors,
      Phi = Phi
    ),
    m = m,
    beta = beta,
    type = "fractional approximation"
  )
  class(output) <- "rSPDEobj"
  return(output)
}

#' Rational approximations of stationary Gaussian Matern random fields
#'
#' `matern.operators` is used for computing a rational SPDE approximation
#' of a stationary Gaussian random fields on \eqn{R^d} with a Matern covariance
#' function
#' \deqn{C(h) = \frac{\sigma^2}{2^{\nu-1}\Gamma(\nu)}
#' (\kappa h)^\nu K_\nu(\kappa h)}{C(h) = (\sigma^2/(2^{\nu-1}\Gamma(\nu))
#' (\kappa h)^\nu K_\nu(\kappa h).}
#'
#' @param kappa Parameter kappa of the SPDE representation. If `NULL`, the range parameter will be used. If the range is also `NULL`, a starting value based on the mesh will be supplied.
#' @param tau Parameter tau of the SPDE representation. If both sigma and tau are `NULL`, a starting value based on the mesh will be supplied.
#' @param alpha Parameter alpha of the SPDE representation. If `alpha` is `NULL`, a starting value will be supplied.
#' @param range Range parameter of the covariance function. Used if `parameterization` is `matern`. If range is `NULL`, a starting value based on the mesh will be supplied.
#' @param sigma Standard deviation of the covariance function. Used if `parameterization` is `matern`. If `NULL`, tau will be used. If tau is also `NULL`, a starting value based on the mesh will be supplied.
#' @param nu Shape parameter of the covariance function. Used if `parameterization` is `matern`. If `NULL`, a starting value will be supplied.
#' @param G The stiffness matrix of a finite element discretization of the
#' domain of interest. Does not need to be given if either `mesh` or `graph` is supplied.
#' @param C The mass matrix of a finite element discretization of the domain
#' of interest. Does not need to be given if either `mesh` or `graph` is supplied.
#' @param mesh An optional fmesher mesh. Replaces `d`, `C` and `G`.
#' @param graph An optional `metric_graph` object. Replaces `d`, `C` and `G`.
#' @param range_mesh The range of the mesh. Will be used to provide starting values for the parameters. Will be used if `mesh` and `graph` are `NULL`, and if one of the parameters (kappa or tau for spde parameterization, or sigma or range for matern parameterization) are not provided.
#' @param loc_mesh The mesh locations used to construct the matrices C and G. This option should be provided if one wants to use the `rspde_lme()` function and will not provide neither graph nor mesh. Only works for 1d data. Does not work for metric graphs. For metric graphs you should supply the graph using the `graph` argument.
#' @param d The dimension of the domain. Does not need to be given if either
#' `mesh` or `graph` is provided.
#' @param m The order of the rational approximation, which needs to be a
#' positive integer. The default value is 1.
#' @param type The type of the rational approximation. The options are
#' "covariance" and "operator". The default is "covariance".
#' @param parameterization Which parameterization to use? `matern` uses range, std. deviation and nu (smoothness). `spde` uses kappa, tau and alpha. The default is `spde`.
#' @param compute_higher_order Logical. Should the higher order finite
#' element matrices be computed?
#' @param return_block_list Logical. For `type = "covariance"`,
#' should the block parts of the precision matrix be returned
#' separately as a list?
#' @param type_rational_approximation Which type of rational
#' approximation should be used? The current types are
#' "brasil", "chebfun" or "chebfunLB".
#' @param compute_logdet Should log determinants be computed while building the model? (For covariance-based models)
#' @return If `type` is "covariance", then `matern.operators`
#' returns an object of class "CBrSPDEobj".
#' This object is a list containing the
#' following quantities:
#' \item{C}{The mass lumped mass matrix.}
#' \item{Ci}{The inverse of `C`.}
#' \item{GCi}{The stiffness matrix G times `Ci`}
#' \item{Gk}{The stiffness matrix G along with the higher-order
#' FEM-related matrices G2, G3, etc.}
#' \item{fem_mesh_matrices}{A list containing the mass
#' lumped mass matrix, the stiffness matrix and
#' the higher-order FEM-related matrices.}
#' \item{m}{The order of the rational approximation.}
#' \item{alpha}{The fractional power of the precision operator.}
#' \item{type}{String indicating the type of approximation.}
#' \item{d}{The dimension of the domain.}
#' \item{nu}{Shape parameter of the covariance function.}
#' \item{kappa}{Range parameter of the covariance function}
#' \item{tau}{Scale parameter of the covariance function.}
#' \item{sigma}{Standard deviation of the covariance function.}
#' \item{type}{String indicating the type of approximation.}
#'
#' If `type` is "operator", then `matern.operators`
#'  returns an object of class "rSPDEobj". This object contains the
#' quantities listed in the output of [fractional.operators()],
#' the `G` matrix, the dimension of the domain, as well as the
#' parameters of the covariance function.
#'
#' @details If `type` is "covariance", we use the
#' covariance-based rational approximation of the fractional operator.
#' In the SPDE approach, we model \eqn{u} as the solution of the following SPDE:
#' \deqn{L^{\alpha/2}(\tau u) = \mathcal{W},}
#' where
#' \eqn{L  = -\Delta +\kappa^2 I} and \eqn{\mathcal{W}} is the standard
#' Gaussian white noise. The covariance operator of \eqn{u} is given
#' by \eqn{L^{-\alpha}}. Now, let \eqn{L_h} be a finite-element
#' approximation of \eqn{L}. We can use
#' a rational approximation of order \eqn{m} on \eqn{L_h^{-\alpha}} to
#' obtain the following approximation:
#' \deqn{L_{h,m}^{-\alpha} = L_h^{-m_\alpha} p(L_h^{-1})q(L_h^{-1})^{-1},}
#' where \eqn{m_\alpha = \lfloor \alpha\rfloor}, \eqn{p} and \eqn{q} are
#' polynomials arising from such rational approximation.
#' From this approximation we construct an approximate precision
#' matrix for \eqn{u}.
#'
#' If `type` is "operator", the approximation is based on a
#' rational approximation of the fractional operator
#' \eqn{(\kappa^2 -\Delta)^\beta}, where \eqn{\beta = (\nu + d/2)/2}.
#' This results in an approximate model of the form \deqn{P_l u(s) = P_r W,}
#' where \eqn{P_j = p_j(L)} are non-fractional operators defined in terms
#' of polynomials \eqn{p_j} for \eqn{j=l,r}. The order of \eqn{p_r} is given
#' by `m` and the order of \eqn{p_l} is \eqn{m + m_\beta}
#' where \eqn{m_\beta} is the integer part of \eqn{\beta} if \eqn{\beta>1} and
#' \eqn{m_\beta = 1} otherwise.
#'
#' The discrete approximation can be written as \eqn{u = P_r x} where
#' \eqn{x \sim N(0,Q^{-1})}{x ~ N(0,Q^{-1})} and \eqn{Q = P_l^T C^{-1} P_l}.
#' Note that the matrices \eqn{P_r} and \eqn{Q} may be be
#' ill-conditioned for \eqn{m>1}. In this case, the methods in
#' [operator.operations()] should be used for operations involving
#' the matrices, since these methods are more numerically stable.
#' @export
#' @seealso [fractional.operators()],
#' [spde.matern.operators()],
#' [matern.operators()]
#'
#' @examples
#' # Compute the covariance-based rational approximation of a
#' # Gaussian process with a Matern covariance function on R
#' kappa <- 10
#' sigma <- 1
#' nu <- 0.8
#' range <- sqrt(8 * nu) / kappa
#'
#' # create mass and stiffness matrices for a FEM discretization
#' nobs <- 101
#' x <- seq(from = 0, to = 1, length.out = 101)
#' fem <- rSPDE.fem1d(x)
#'
#' # compute rational approximation of covariance function at 0.5
#' op_cov <- matern.operators(
#'   loc_mesh = x, nu = nu,
#'   range = range, sigma = sigma, d = 1, m = 2,
#'   parameterization = "matern"
#' )
#'
#' v <- t(rSPDE.A1d(x, 0.5))
#' # Compute the precision matrix
#' Q <- op_cov$Q
#' # A matrix here is the identity matrix
#' A <- Diagonal(nobs)
#' # We need to concatenate 3 A's since we are doing a covariance-based rational
#' # approximation of order 2
#' Abar <- cbind(A, A, A)
#' w <- rbind(v, v, v)
#' # The approximate covariance function:
#' c_cov.approx <- (Abar) %*% solve(Q, w)
#' c.true <- folded.matern.covariance.1d(
#'   rep(0.5, length(x)),
#'   abs(x), kappa, nu, sigma
#' )
#'
#' # plot the result and compare with the true Matern covariance
#' plot(x, c.true,
#'   type = "l", ylab = "C(h)",
#'   xlab = "h", main = "Matern covariance and rational approximations"
#' )
#' lines(x, c_cov.approx, col = 2)
#'
#'
#' # Compute the operator-based rational approximation of a Gaussian
#' # process with a Matern covariance function on R
#' kappa <- 10
#' sigma <- 1
#' nu <- 0.8
#' range <- sqrt(8 * nu) / kappa
#'
#' # create mass and stiffness matrices for a FEM discretization
#' x <- seq(from = 0, to = 1, length.out = 101)
#' fem <- rSPDE.fem1d(x)
#'
#' # compute rational approximation of covariance function at 0.5
#' op <- matern.operators(
#'   range = range, sigma = sigma, nu = nu,
#'   loc_mesh = x, d = 1,
#'   type = "operator",
#'   parameterization = "matern"
#' )
#'
#' v <- t(rSPDE.A1d(x, 0.5))
#' c.approx <- Sigma.mult(op, v)
#' c.true <- folded.matern.covariance.1d(
#'   rep(0.5, length(x)),
#'   abs(x), kappa, nu, sigma
#' )
#'
#' # plot the result and compare with the true Matern covariance
#' plot(x, c.true,
#'   type = "l", ylab = "C(h)",
#'   xlab = "h", main = "Matern covariance and rational approximation"
#' )
#' lines(x, c.approx, col = 2)
matern.operators <- function(kappa = NULL,
                             tau = NULL,
                             alpha = NULL,
                             sigma = NULL,
                             range = NULL,
                             nu = NULL,
                             G = NULL,
                             C = NULL,
                             d = NULL,
                             mesh = NULL,
                             graph = NULL,
                             range_mesh = NULL,
                             loc_mesh = NULL,
                             m = 1,
                             type = c("covariance", "operator"),
                             parameterization = c("spde", "matern"),
                             compute_higher_order = FALSE,
                             return_block_list = FALSE,
                             type_rational_approximation = c(
                               "brasil",
                               "chebfun",
                               "chebfunLB"
                             ),
                             compute_logdet = FALSE) {
  type <- type[[1]]

  if (!type %in% c("covariance", "operator")) {
    stop("The type should be 'covariance' or 'operator'!")
  }
  if (is.null(d) && is.null(mesh) && is.null(graph)) {
    stop("You should give either the dimension d, the mesh or graph!")
  }

  mesh_1d <- NULL

  if ((is.null(C) || is.null(G)) && is.null(mesh) && is.null(graph) && (is.null(loc_mesh) || d != 1)) {
    stop("You should either provide mesh, graph, or provide both C *and* G!")
  }

  mesh_1d <- NULL

  if (!is.null(loc_mesh) && d != 1) {
    stop("loc_mesh only works with dimension 1.")
  } else if (!is.null(loc_mesh)) {
    # mesh_1d <- fmesher::fm_mesh_1d(loc_mesh)
    mesh_1d <- fm_mesh_1d(loc_mesh)
  }

  if ((is.null(C) || is.null(G)) && (is.null(graph)) && (!is.null(loc_mesh) && d == 1)) {
    # fem <- rSPDE.fem1d(loc_mesh)
    # fem <- fmesher::fm_fem(mesh_1d)
    fem <- fm_fem(mesh_1d)
    C <- fem$c0
    G <- fem$g1
  }

  has_mesh <- FALSE
  has_graph <- FALSE

  if (!is.null(loc_mesh)) {
    if (!is.numeric(loc_mesh)) {
      stop("loc_mesh must be numerical.")
    }
    range_mesh <- abs(diff(range(loc_mesh)))
  }

  parameterization <- parameterization[[1]]

  if (!parameterization %in% c("matern", "spde")) {
    stop("parameterization should be either 'matern', 'spde' or 'graph'!")
  }

  if (parameterization == "spde") {
    if (!is.null(nu)) {
      stop("For 'spde' parameterization, you should NOT supply 'nu'. You need to provide 'alpha'!")
    }
    if (!is.null(sigma)) {
      stop("For 'spde' parameterization, you should NOT supply 'sigma'. You need to provide 'tau'!")
    }
    if (!is.null(range)) {
      stop("For 'spde' parameterization, you should NOT supply 'range'. You need to provide 'kappa'!")
    }
  } else {
    if (!is.null(alpha)) {
      stop("For 'matern' parameterization, you should NOT supply 'alpha'. You need to provide 'nu'!")
    }
    if (!is.null(tau)) {
      stop("For 'matern' parameterization, you should NOT supply 'tau'. You need to provide 'sigma'!")
    }
    if (!is.null(kappa)) {
      stop("For 'matern' parameterization, you should NOT supply 'kappa'. You need to provide 'range'!")
    }
  }

  if (!is.null(graph)) {
    if (!inherits(graph, "metric_graph")) {
      stop("graph should be a metric_graph object!")
    }
    d <- 1
    if (is.null(graph$mesh)) {
      warning("The graph object did not contain a mesh, one was created with h = 0.01. Use the build_mesh() method to replace this mesh with a new one.")
      graph$build_mesh(h = 0.01)
    }
    graph$compute_fem()
    C <- graph$mesh$C
    G <- graph$mesh$G
    has_graph <- TRUE
  }

  if (!is.null(mesh)) {
    d <- fmesher::fm_manifold_dim(mesh)
    # fem <- INLA::inla.mesh.fem(mesh)
    # fem <- fmesher::fm_fem(mesh)
    fem <- fm_fem(mesh)
    C <- fem$c0
    G <- fem$g1
    has_mesh <- TRUE
  }

  if (parameterization == "spde") {
    if (is.null(kappa) || is.null(tau)) {
      if (is.null(graph) && is.null(mesh) && is.null(range_mesh)) {
        stop("You should either provide all the parameters, or you should provide one of the following: mesh, range_mesh or graph.")
      }
      if (!is.null(mesh) || !is.null(range_mesh)) {
        param <- get.initial.values.rSPDE(mesh.range = range_mesh, dim = d, parameterization = parameterization, mesh = mesh, nu = nu)
      } else if (!is.null(graph)) {
        param <- get.initial.values.rSPDE(graph.obj = graph, parameterization = parameterization, nu = nu)
      }
    }

    if (is.null(kappa)) {
      kappa <- exp(param[2])
    } else {
      kappa <- rspde_check_user_input(kappa, "kappa", 0)
    }
    if (is.null(tau)) {
      tau <- exp(param[1])
    } else {
      tau <- rspde_check_user_input(tau, "tau", 0)
    }

    if (is.null(alpha)) {
      alpha <- 1 + d / 2
    } else {
      alpha <- rspde_check_user_input(alpha, "alpha", d / 2)
      alpha <- min(alpha, 10)
    }

    nu <- alpha - d / 2

    range <- sqrt(8 * nu) / kappa
    sigma <- sqrt(gamma(nu) / (tau^2 * kappa^(2 * nu) *
      (4 * pi)^(d / 2) * gamma(nu + d / 2)))
  } else if (parameterization == "matern") {
    if (is.null(sigma) || is.null(range)) {
      if (is.null(graph) && is.null(mesh) && is.null(range_mesh)) {
        stop("You should either provide all the parameters, or you should provide one of the following: mesh, range_mesh or graph.")
      }
      if (!is.null(mesh) || !is.null(range_mesh)) {
        param <- get.initial.values.rSPDE(mesh.range = range_mesh, dim = d, parameterization = parameterization, mesh = mesh, nu = nu)
      } else if (!is.null(graph)) {
        param <- get.initial.values.rSPDE(graph.obj = graph, parameterization = parameterization, nu = nu)
      }
    }

    if (is.null(range)) {
      range <- exp(param[2])
    } else {
      range <- rspde_check_user_input(range, "range", 0)
    }
    if (is.null(sigma)) {
      sigma <- exp(param[1])
    } else {
      sigma <- rspde_check_user_input(sigma, "sigma", 0)
    }

    if (is.null(nu)) {
      if(d == 1){
        nu <- 0.75
      } else {
        nu <- 1
      }
    } else {
      nu <- rspde_check_user_input(nu, "nu", 0)
    }
    kappa <- sqrt(8 * nu) / range
    tau <- sqrt(gamma(nu) / (sigma^2 * kappa^(2 * nu) *
      (4 * pi)^(d / 2) * gamma(nu + d / 2)))
    alpha <- nu + d / 2
  }

  # else if(parameterization == "graph"){
  #   if(is.null(kappa) || is.null(sigma)){
  #     if(is.null(graph) && is.null(mesh) && is.null(range_mesh)){
  #       stop("You should either provide all the parameters, or you should provide one of the following: mesh, range_mesh or graph.")
  #     }
  #     if(!is.null(mesh) || !is.null(range_mesh)){
  #        param <- get.initial.values.rSPDE(mesh.range = range_mesh, dim = d, parameterization = "spde", mesh = mesh, nu = nu)
  #     } else if(!is.null(graph)){
  #       param <- get.initial.values.rSPDE(graph.obj = graph, parameterization = "spde", nu = nu)
  #     }
  #   }

  #   if(is.null(kappa)){
  #     kappa <- exp(param[2])
  #   }
  #   if(is.null(sigma)){
  #     tau <- exp(param[1])
  #     sigma <- 1 / tau
  #   }
  #   range <- sqrt(8 * nu) / kappa
  #   tau <- 1 / sigma
  # }

  if (type == "operator") {
    beta <- (nu + d / 2) / 2
    operators <- fractional.operators(
      L = G + C * kappa^2,
      beta = beta,
      C = C,
      scale.factor = kappa^2,
      m = m,
      tau = tau
    )
    output <- operators
    output$kappa <- kappa
    output$range <- range
    output$tau <- tau
    output$sigma <- sigma
    output$alpha <- alpha
    output$nu <- nu
    output$d <- d
    output$G <- G
    output$type <- "Matern approximation"
    output$stationary <- TRUE
    output$parameterization <- parameterization
    output$has_mesh <- has_mesh
    output$has_graph <- has_graph
    output$alpha <- 2 * beta
    output$mesh <- mesh
    output$range_mesh <- range_mesh
    output$graph <- graph
    output$loc_mesh <- loc_mesh
    output$mesh_1d <- mesh_1d
    class(output) <- c("matern_operator", class(output))
    return(output)
  } else {
    C <- Matrix::Diagonal(dim(C)[1], rowSums(C))
    type_rational_approximation <- type_rational_approximation[[1]]
    out <- CBrSPDE.matern.operators(
      C = C, G = G, mesh = mesh, nu = nu, kappa = kappa, tau = tau,
      m = m, d = d, compute_higher_order = compute_higher_order,
      return_block_list = return_block_list,
      type_rational_approximation = type_rational_approximation,
      compute_logdet = compute_logdet
    )
    out$range <- range
    out$tau <- tau
    out$alpha <- alpha
    out$sigma <- sigma
    out$parameterization <- parameterization
    out$has_mesh <- has_mesh
    out$has_graph <- has_graph
    out$mesh <- mesh
    out$range_mesh <- range_mesh
    out$graph <- graph
    out$loc_mesh <- loc_mesh
    out$mesh_1d <- mesh_1d
    class(out) <- c("matern_operator", class(out))
    return(out)
  }
}



#' @name CBrSPDE.matern.operators
#' @title Covariance-based rational approximations of stationary Gaussian
#' Matern random fields
#' @description `CBrSPDE.matern.operators` is used for computing a
#' covariance-based rational SPDE approximation of stationary Gaussian random
#' fields on \eqn{R^d} with a Matern covariance function
#' \deqn{C(h) = \frac{\sigma^2}{2^{\nu-1}\Gamma(\nu)}(\kappa h)^\nu
#' K_\nu(\kappa h)}{C(h) =
#' (\sigma^2/(2^{\nu-1}\Gamma(\nu))(\kappa h)^\nu K_\nu(\kappa h)}
#' @param kappa Range parameter of the covariance function.
#  @param sigma Standard deviation of the covariance function.
#' @param tau Precision parameter.
#' @param nu Shape parameter of the covariance function.
#' @param G The stiffness matrix of a finite element discretization
#' of the domain of interest.
#' @param C The mass matrix of a finite element discretization of
#' the domain of interest.
#' @param mesh An inla mesh.
#' @param d The dimension of the domain.
#' @param m The order of the rational approximation, which needs
#' to be a positive integer.
#' The default value is 2.
#' @param compute_logdet Should log determinants be computed while building the model?
#' @return `CBrSPDE.matern.operators` returns an object of
#' class "CBrSPDEobj". This object is a list containing the
#' following quantities:
#' \item{C}{The mass lumped mass matrix.}
#' \item{Ci}{The inverse of `C`.}
#' \item{GCi}{The stiffness matrix G times `Ci`}
#' \item{Gk}{The stiffness matrix G along with the higher-order
#' FEM-related matrices G2, G3, etc.}
#' \item{fem_mesh_matrices}{A list containing the mass lumped mass
#' matrix, the stiffness matrix and
#' the higher-order FEM-related matrices.}
#' \item{m}{The order of the rational approximation.}
#' \item{alpha}{The fractional power of the precision operator.}
#' \item{type}{String indicating the type of approximation.}
#' \item{d}{The dimension of the domain.}
#' \item{nu}{Shape parameter of the covariance function.}
#' \item{kappa}{Range parameter of the covariance function}
#' \item{tau}{Scale parameter of the covariance function.}
#' \item{sigma}{Standard deviation of the covariance function.}
#' \item{type}{String indicating the type of approximation.}
#' @noRd
#' @seealso [matern.operators()], [spde.matern.operators()]
#' @details We use the covariance-based rational approximation of the
#' fractional operator. In the SPDE approach, we model \eqn{u} as the
#' solution of the following SPDE:
#' \deqn{L^{\alpha/2}(\tau u) = \mathcal{W},}
#' where
#' \eqn{L  = -\Delta +\kappa^2 I} and \eqn{\mathcal{W}} is
#' the standard Gaussian white noise. The covariance operator of \eqn{u}
#' is given by \eqn{L^{-\alpha}}. Now, let \eqn{L_h} be a finite-element
#' approximation of \eqn{L}. We can use a rational approximation of order
#' \eqn{m} on \eqn{L_h^{-\alpha}} to obtain the following approximation:
#' \deqn{L_{h,m}^{-\alpha} = L_h^{-m_\alpha}
#' p(L_h^{-1})q(L_h^{-1})^{-1},}
#' where \eqn{m_\alpha = \max\{1,\lfloor \alpha\rfloor\}},
#' \eqn{p} and \eqn{q} are polynomials arising from such rational approximation.
#' From this approximation we construct an approximate precision
#' matrix for \eqn{u}.
#'
#'
#' @examples
#' # Compute the covariance-based rational approximation of a
#' # Gaussian process with a Matern covariance function on R
#' kappa <- 10
#' sigma <- 1
#' nu <- 0.8
#'
#' # create mass and stiffness matrices for a FEM discretization
#' nobs <- 101
#' x <- seq(from = 0, to = 1, length.out = 101)
#' fem <- rSPDE.fem1d(x)
#'
#' # compute rational approximation of covariance function at 0.5
#' tau <- sqrt(gamma(nu) / (sigma^2 * kappa^(2 * nu) *
#'   (4 * pi)^(1 / 2) * gamma(nu + 1 / 2)))
#' op_cov <- CBrSPDE.matern.operators(
#'   C = fem$C, G = fem$G, nu = nu,
#'   kappa = kappa, tau = tau, d = 1, m = 2
#' )
#'
#' v <- t(rSPDE.A1d(x, 0.5))
#' # Compute the precision matrix
#' Q <- op_cov$Q
#' # A matrix here is the identity matrix
#' A <- Diagonal(nobs)
#' # We need to concatenate 3 A's since we are doing a covariance-based rational
#' # approximation of order 2
#' Abar <- cbind(A, A, A)
#' w <- rbind(v, v, v)
#' # The approximate covariance function:
#' c_cov.approx <- (Abar) %*% solve(Q, w)
#'
#' # plot the result and compare with the true Matern covariance
#' plot(x, matern.covariance(abs(x - 0.5), kappa, nu, sigma),
#'   type = "l", ylab = "C(h)",
#'   xlab = "h", main = "Matern covariance and rational approximations"
#' )
#' lines(x, c_cov.approx, col = 2)
CBrSPDE.matern.operators <- function(C,
                                     G,
                                     mesh,
                                     nu,
                                     kappa,
                                     tau,
                                     m = 2,
                                     d,
                                     compute_higher_order = FALSE,
                                     return_block_list = FALSE,
                                     type_rational_approximation = c(
                                       "chebfun",
                                       "brasil", "chebfunLB"
                                     ),
                                     fem_mesh_matrices = NULL,
                                     compute_logdet = FALSE) {
  type_rational_approximation <- type_rational_approximation[[1]]

  if (is.null(fem_mesh_matrices)) {
    if (!is.null(mesh)) {
      ## get alpha, m_alpha
      d <- fmesher::fm_manifold_dim(mesh)
      alpha <- nu + d / 2
      m_alpha <- floor(alpha)
      m_order <- m_alpha + 1
      # tau <- sqrt(gamma(nu) / (sigma^2 * kappa^(2 * nu) *
      # (4 * pi)^(d / 2) * gamma(nu + d / 2)))

      if (d > 1) {
        if (compute_higher_order) {
          # fem <- INLA::inla.mesh.fem(mesh, order = m_alpha + 1)
          # fem <- fmesher::fm_fem(mesh, order = m_alpha + 1)
          fem <- fm_fem(mesh, order = m_alpha + 1)
        } else {
          # fem <- INLA::inla.mesh.fem(mesh)
          # fem <- fmesher::fm_fem(mesh)
          fem <- fm_fem(mesh)
        }

        C <- fem$c0
        G <- fem$g1
        Ci <- Matrix::Diagonal(dim(C)[1], 1 / rowSums(C))
        GCi <- G %*% Ci

        Gk <- list()
        Gk[[1]] <- G
        for (i in 2:m_order) {
          Gk[[i]] <- fem[[paste0("g", i)]]
        }
      } else if (d == 1) {
        # fem <- INLA::inla.mesh.fem(mesh, order = 2)
        # fem <- fmesher::fm_fem(mesh)
        fem <- fm_fem(mesh)
        C <- fem$c0
        G <- fem$g1
        Ci <- Matrix::Diagonal(dim(C)[1], 1 / rowSums(C))
        GCi <- G %*% Ci

        Gk <- list()
        Gk[[1]] <- G
        if (compute_higher_order) {
          Gk[[2]] <- fem$g2
          if (m_order > 2) {
            for (i in 3:m_order) {
              Gk[[i]] <- GCi %*% Gk[[i - 1]]
            }
          }
        }
      }
    } else {
      ## get alpha, m_alpha
      alpha <- nu + d / 2
      m_alpha <- floor(alpha)
      m_order <- m_alpha + 1
      # tau <- sqrt(gamma(nu) / (sigma^2 * kappa^(2 * nu) *
      # (4 * pi)^(d / 2) * gamma(nu + d / 2)))

      ## get lumped mass matrix
      C <- Matrix::Diagonal(dim(C)[1], rowSums(C))

      ## get G_k matrix: k is up to m_alpha if alpha is integer,
      # k is up to m_alpha + 1 otherwise.
      # inverse lumped mass matrix
      Ci <- Matrix::Diagonal(dim(C)[1], 1 / rowSums(C))

      GCi <- G %*% Ci
      # create a list to store all the G_k matrix

      Gk <- list()

      Gk[[1]] <- G
      # determine how many G_k matrices we want to create
      if (compute_higher_order) {
        for (i in 2:m_order) {
          Gk[[i]] <- GCi %*% Gk[[i - 1]]
        }
      }
    }

    # create a list contains all the finite element related matrices
    fem_mesh_matrices <- list()
    fem_mesh_matrices[["c0"]] <- C
    fem_mesh_matrices[["g1"]] <- G

    if (compute_higher_order) {
      for (i in 1:m_order) {
        fem_mesh_matrices[[paste0("g", i)]] <- Gk[[i]]
      }
    }
  } else {
    C <- fem_mesh_matrices$c0
    G <- fem_mesh_matrices$g1
    Ci <- Matrix::Diagonal(dim(C)[1], 1 / rowSums(C))
    GCi <- G %*% Ci
    alpha <- nu + d / 2
    m_alpha <- floor(alpha)
    m_order <- m_alpha + 1
    if (!is.null(mesh)) {
      d <- fmesher::fm_manifold_dim(mesh)
    }
    # tau <- sqrt(gamma(nu) / (sigma^2 * kappa^(2 * nu) *
    # (4 * pi)^(d / 2) * gamma(nu + d / 2)))
    Gk <- list()
    Gk[[1]] <- G
    for (i in 2:m_order) {
      Gk[[i]] <- fem_mesh_matrices[[paste0("g", i)]]
    }
  }


  L <- (G + kappa^2 * C) / kappa^2

  if (compute_logdet) {
    Lchol <- Matrix::Cholesky(L)
    logdetL <- 2 * c(determinant(Lchol, logarithm = TRUE, sqrt = TRUE)$modulus)

    logdetC <- sum(log(diag(C)))
  } else {
    logdetL <- NULL
    logdetC <- NULL
  }

  CiL <- GCi / kappa^2 + Diagonal(dim(GCi)[1])

  if (m_alpha == 0) {
    aux_mat <- Diagonal(dim(L)[1])
  } else {
    aux_mat <- CiL
  }


  if (return_block_list) {
    Q.int <- aux_mat

    if (alpha %% 1 == 0) {
      Q.frac <- Matrix::Diagonal(dim(L)[1])
      Q <- L * kappa^2

      if (alpha > 1) {
        CinvL <- Matrix::Diagonal(dim(C)[1], 1 / diag(C)) %*% L * kappa^2

        for (k in 1:(alpha - 1)) {
          Q <- Q %*% CinvL
        }
      }

      Q <- tau^2 * Q
      Q.int <- Q
    } else {
      if (m == 0) {
        stop("Return block list does not work with m = 0, either increase m or set return_block_list to FALSE.")
      }
      Q.frac <- rspde.matern.precision(
        kappa = kappa, nu = nu, tau = tau,
        rspde.order = m, dim = d,
        fem_mesh_matrices = fem_mesh_matrices, only_fractional = TRUE,
        return_block_list = TRUE,
        type_rational_approx = type_rational_approximation
      )

      Q <- Q.frac

      if (m_alpha > 0) {
        for (j in seq_len(length(Q))) {
          for (i in 1:m_alpha) {
            Q[[j]] <- Q.int %*% Q[[j]]
          }
        }
      }
      Q.int <- list(Q.int = Q.int, order = m_alpha)
    }
  } else {
    Q.int <- list(Q.int = kronecker(Diagonal(m + 1), aux_mat), order = m_alpha)

    if (alpha %% 1 == 0) {
      Q.frac <- Matrix::Diagonal(dim(L)[1])
      Q <- L * kappa^2

      if (alpha > 1) {
        CinvL <- Matrix::Diagonal(dim(C)[1], 1 / diag(C)) %*% L * kappa^2

        for (k in 1:(alpha - 1)) {
          Q <- Q %*% CinvL
        }
      }

      Q <- tau^2 * Q
      Q.int <- list(Q.int = Q, order = m_alpha)
    } else if (m > 0) {
      Q.frac <- rspde.matern.precision(
        kappa = kappa, nu = nu, tau = tau,
        rspde.order = m, dim = d,
        fem_mesh_matrices = fem_mesh_matrices, only_fractional = TRUE,
        type_rational_approx = type_rational_approximation
      )

      Q <- Q.frac

      if (m_alpha > 0) {
        for (i in 1:m_alpha) {
          Q <- Q.int$Q.int %*% Q
        }
      }
    } else {
      Q <- markov.Q(alpha / 2, kappa, d, list(C = C, G = G))
      Q.frac <- Matrix::Diagonal(dim(L)[1])
    }
  }

  ## output
  output <- list(
    C = C, G = G, L = L, Ci = Ci, GCi = GCi, Gk = Gk,
    fem_mesh_matrices = fem_mesh_matrices,
    alpha = alpha, nu = nu, kappa = kappa, range = sqrt(8 * max(nu,0)) / kappa,
    tau = tau, m = m, d = d,
    # sigma = sigma,
    logdetL = logdetL, logdetC = logdetC,
    Q.frac = Q.frac, Q.int = Q.int,
    Q = Q, sizeC = dim(C)[1],
    higher_order = compute_higher_order,
    type_rational_approximation = type_rational_approximation,
    return_block_list = return_block_list,
    stationary = TRUE
  )
  output$type <- "Covariance-Based Matern SPDE Approximation"
  output$compute_logdet <- compute_logdet
  class(output) <- "CBrSPDEobj"
  return(output)
}


#' Rational approximations of non-stationary Gaussian SPDE Matern random fields
#'
#' `spde.matern.operators` is used for computing a rational SPDE
#' approximation of a Gaussian random
#' fields on \eqn{R^d} defined as a solution to the SPDE
#' \deqn{(\kappa(s) - \Delta)^\beta (\tau(s)u(s)) = W.}
#'
#' @param kappa Vector with the, possibly spatially varying, range parameter
#' evaluated at the locations of the mesh used for the finite element
#' discretization of the SPDE.
#' @param tau Vector with the, possibly spatially varying, precision
#' parameter evaluated at the locations
#' of the mesh used for the finite element discretization of the SPDE.
#' @param theta Theta parameter that connects B.tau and B.kappa to tau and kappa through a log-linear regression, in case the parameterization is `spde`,
#' and that connects B.sigma and B.range to tau and kappa in case the parameterization is `matern`. When both tau and kappa are constant after this
#' transformation, `spde.matern.operators()` delegates to `matern.operators()`.
#' @param B.sigma Matrix with specification of log-linear model for \eqn{\sigma}. Will be used if `parameterization = 'matern'`.
#' @param B.range Matrix with specification of log-linear model for \eqn{\rho}, which is a range-like parameter (it is exactly the range parameter in the stationary case). Will be used if `parameterization = 'matern'`.
#' @param parameterization Which parameterization to use? `matern` uses range, std. deviation and nu (smoothness). `spde` uses kappa, tau and nu (smoothness). The default is `matern`.
#' @param B.tau Matrix with specification of log-linear model for \eqn{\tau}. Will be used if `parameterization = 'spde'`.
#' @param B.kappa Matrix with specification of log-linear model for \eqn{\kappa}. Will be used if `parameterization = 'spde'`.
#' @param nu Shape parameter of the covariance function. Will be used if the parameterization is 'matern'.
#' @param alpha smoothness parameter. Will be used if the parameterization is 'spde'.
#' @param G The stiffness matrix of a finite element discretization of
#' the domain of interest.
#' @param C The mass matrix of a finite element discretization of the
#' domain of interest.
#' @param mesh An optional inla mesh. `d`, `C` and `G`
#' must be given if `mesh` is not given.
#' @param graph An optional `metric_graph` object. Replaces `d`, `C` and `G`.
#' @param range_mesh The range of the mesh. Will be used to provide starting values for the parameters. Will be used if `mesh` and `graph` are `NULL`, and if one of the parameters (kappa or tau for spde parameterization, or sigma or range for matern parameterization) are not provided.
#' @param loc_mesh The mesh locations used to construct the matrices C and G. This option should be provided if one wants to use the `rspde_lme()` function and will not provide neither graph nor mesh. Only works for 1d data. Does not work for metric graphs. For metric graphs you should supply the graph using the `graph` argument.
#' @param d The dimension of the domain. Does not need to be given if
#' `mesh` is used.
#' @param m The order of the rational approximation, which needs to be a
#' positive integer. The default value is 1.
#' @param type The type of the rational approximation. The options are
#' "covariance" and "operator". The default is "covariance".
#' @param type_rational_approximation Which type of rational
#' approximation should be used? The current types are
#' "brasil", "chebfun" or "chebfunLB".
#'
#'
#' @details The approximation is based on a rational approximation of the
#' fractional operator \eqn{(\kappa(s)^2 -\Delta)^\beta}, where
#' \eqn{\beta = (\nu + d/2)/2}. This results in an approximate model
#' on the form \deqn{P_l u(s) = P_r W,} where \eqn{P_j = p_j(L)} are
#' non-fractional operators defined in terms of polynomials \eqn{p_j} for
#' \eqn{j=l,r}. The order of \eqn{p_r} is given by `m` and the order
#' of \eqn{p_l} is \eqn{m + m_\beta} where \eqn{m_\beta} is the integer
#' part of \eqn{\beta} if \eqn{\beta>1} and \eqn{m_\beta = 1} otherwise.
#'
#' The discrete approximation can be written as \eqn{u = P_r x} where
#' \eqn{x \sim N(0,Q^{-1})}{x ~ N(0,Q^{-1})}
#' and \eqn{Q = P_l^T C^{-1} P_l}. Note that the matrices \eqn{P_r} and
#' \eqn{Q} may be be ill-conditioned for \eqn{m>1}.
#' In this case, the metehods in [operator.operations()]
#' should be used for operations involving the matrices, since
#' these methods are more numerically stable.
#'
#' @return `spde.matern.operators` returns an object of
#' class "rSPDEobj. This object contains the
#' quantities listed in the output of [fractional.operators()]
#' as well as the smoothness parameter \eqn{\nu}.
#' @export
#' @seealso [fractional.operators()],
#' [spde.matern.operators()],
#' [matern.operators()]
#'
#' @examples
#' # Sample non-stationary Matern field on R
#' tau <- 1
#' nu <- 0.8
#'
#' # create mass and stiffness matrices for a FEM discretization
#' x <- seq(from = 0, to = 1, length.out = 101)
#' fem <- rSPDE.fem1d(x)
#'
#' # define a non-stationary range parameter
#' kappa <- seq(from = 2, to = 20, length.out = length(x))
#' alpha <- nu + 1 / 2
#' # compute rational approximation
#' op <- spde.matern.operators(
#'   kappa = kappa, tau = tau, alpha = alpha,
#'   G = fem$G, C = fem$C, d = 1
#' )
#'
#' # sample the field
#' u <- simulate(op)
#'
#' # plot the sample
#' plot(x, u, type = "l", ylab = "u(s)", xlab = "s")
#'
spde.matern.operators <- function(kappa = NULL,
                                  tau = NULL,
                                  theta = NULL,
                                  B.tau = matrix(c(0, 1, 0), 1, 3),
                                  B.kappa = matrix(c(0, 0, 1), 1, 3),
                                  B.sigma = matrix(c(0, 1, 0), 1, 3),
                                  B.range = matrix(c(0, 0, 1), 1, 3),
                                  alpha = NULL,
                                  nu = NULL,
                                  parameterization = c("spde", "matern"),
                                  G = NULL,
                                  C = NULL,
                                  d = NULL,
                                  graph = NULL,
                                  mesh = NULL,
                                  range_mesh = NULL,
                                  loc_mesh = NULL,
                                  m = 1,
                                  type = c("covariance", "operator"),
                                  type_rational_approximation = c(
                                    "brasil",
                                    "chebfun",
                                    "chebfunLB"
                                  )) {
  is_constant_param <- function(x) {
    x <- as.numeric(x)
    length(unique(x)) == 1
  }

  type <- type[[1]]
  if (!type %in% c("covariance", "operator")) {
    stop("The type should be 'covariance' or 'operator'!")
  }

  has_mesh <- FALSE
  has_graph <- FALSE

  parameterization <- parameterization[[1]]

  if (!parameterization %in% c("matern", "spde")) {
    stop("parameterization should be either 'matern' or 'spde'!")
  }

  if (parameterization == "spde") {
    if (!is.null(nu)) {
      stop("For 'spde' parameterization, you should NOT supply 'nu'. You need to provide 'alpha'!")
    }
    if (!missing(B.sigma)) {
      stop("For 'spde' parameterization, you should NOT supply 'B.sigma'. You need to provide 'B.tau'!")
    }
    if (!missing(B.range)) {
      stop("For 'spde' parameterization, you should NOT supply 'B.range'. You need to provide 'B.kappa'!")
    }
  } else {
    if (!is.null(alpha)) {
      stop("For 'matern' parameterization, you should NOT supply 'alpha'. You need to provide 'nu'!")
    }
    if (!missing(B.tau)) {
      stop("For 'matern' parameterization, you should NOT supply 'B.tau'. You need to provide 'B.sigma'!")
    }
    if (!missing(B.kappa)) {
      stop("For 'matern' parameterization, you should NOT supply 'B.kappa'. You need to provide 'B.range'!")
    }
  }

  if (is.null(d) && is.null(mesh) && is.null(graph)) {
    stop("You should give either the dimension d, the mesh or graph!")
  }

  if ((is.null(C) || is.null(G)) && is.null(mesh) && is.null(graph) && (is.null(loc_mesh) || d != 1)) {
    stop("You should either provide mesh, graph, or provide both C *and* G!")
  }

  mesh_1d <- NULL

  if (!is.null(loc_mesh) && d != 1) {
    stop("loc_mesh only works with dimension 1.")
  } else if (!is.null(loc_mesh)) {
    # mesh_1d <- fmesher::fm_mesh_1d(loc_mesh)
    mesh_1d <- fm_mesh_1d(loc_mesh)
  }

  if ((is.null(C) || is.null(G)) && (is.null(graph)) && (!is.null(loc_mesh) && d == 1)) {
    # fem <- rSPDE.fem1d(loc_mesh)
    fem <- fm_fem(mesh_1d)
    C <- fem$c0
    G <- fem$g1
  }

  has_mesh <- FALSE
  has_graph <- FALSE

  if (!is.null(loc_mesh)) {
    if (!is.numeric(loc_mesh)) {
      stop("loc_mesh must be numerical.")
    }
    range_mesh <- abs(diff(range(loc_mesh)))
  }

  if (!is.null(mesh)) {
    d <- fmesher::fm_manifold_dim(mesh)
    # fem <- INLA::inla.mesh.fem(mesh)
    # fem <- fmesher::fm_fem(mesh)
    fem <- fm_fem(mesh)
    C <- fem$c0
    G <- fem$g1
    has_mesh <- TRUE
  }

  if (!is.null(graph)) {
    if (!inherits(graph, "metric_graph")) {
      stop("graph should be a metric_graph object!")
    }
    has_graph <- TRUE
    d <- 1
    if (is.null(graph$mesh)) {
      warning("The graph object did not contain a mesh, one was created with h = 0.01. Use the build_mesh() method to replace this mesh with a new one.")
      graph$build_mesh(h = 0.01)
    }
    graph$compute_fem()
    C <- graph$mesh$C
    G <- graph$mesh$G
  }

  if (!is.null(theta)) {
    if (parameterization == "matern") {
      if (any(dim(B.sigma) != dim(B.range))) {
        stop("B.sigma and B.range must have the same dimensions!")
      }

      B_matrices <- convert_B_matrices(B.sigma, B.range, ncol(C), nu, d)
      B.tau <- B_matrices[["B.tau"]]
      B.kappa <- B_matrices[["B.kappa"]]

      if (is.null(nu)) {
        if(d == 1){
          nu <- 0.75
        } else {
          nu <- 1
        }
      } else {
        nu <- rspde_check_user_input(nu, "nu", 0)
      }

      alpha <- nu + d / 2
    } else if (parameterization == "spde") {
      if (any(dim(B.tau) != dim(B.kappa))) {
        stop("B.tau and B.kappa must have the same dimensions!")
      }

      B.tau <- prepare_B_matrices(
        B.tau, ncol(C),
        ncol(B.tau) - 1
      )
      B.kappa <- prepare_B_matrices(
        B.kappa, ncol(C),
        ncol(B.kappa) - 1
      )

      if (is.null(alpha)) {
        alpha <- 1 + d / 2
      } else {
        alpha <- rspde_check_user_input(alpha, "alpha", d / 2)
        alpha <- min(alpha, 10)
      }

      nu <- alpha - d / 2
    }

    # else if(parameterization == "graph") {
    #   if(any(dim(B.sigma) != dim(B.kappa))){
    #     stop("B.sigma and B.kappa must have the same dimensions!")
    #   }

    #   B.sigma <- prepare_B_matrices(B.sigma, ncol(C),
    #       ncol(B.sigma)-1)
    #   B.kappa <- prepare_B_matrices(B.kappa, ncol(C),
    #       ncol(B.kappa)-1)
    #   B.tau <- -B.sigma
    # }

    if (is.null(tau) || is.null(kappa)) {
      new_theta <- c(1, theta)

      tau <- c(exp(B.tau %*% new_theta))
      kappa <- c(exp(B.kappa %*% new_theta))
    }
  } else if (is.null(kappa) || is.null(tau)) {
    if (parameterization == "matern") {
      if (any(dim(B.sigma) != dim(B.range))) {
        stop("B.sigma and B.range must have the same dimensions!")
      }
      if (is.null(nu)) {
        if(d == 1){
          nu <- 0.75
        } else {
          nu <- 1
        }
      } else {
        nu <- rspde_check_user_input(nu, "nu", 0)
      }
      B_matrices <- convert_B_matrices(B.sigma, B.range, ncol(C), nu, d)
      B.tau <- B_matrices[["B.tau"]]
      B.kappa <- B_matrices[["B.kappa"]]

      if (is.null(nu)) {
        if(d == 1){
          nu <- 0.75
        } else {
          nu <- 1
        }
      } else {
        nu <- rspde_check_user_input(nu, "nu", 0)
      }

      alpha <- nu + d / 2
    } else {
      if (any(dim(B.tau) != dim(B.kappa))) {
        stop("B.tau and B.kappa must have the same dimensions!")
      }
      B.tau <- prepare_B_matrices(
        B.tau, ncol(C),
        ncol(B.tau) - 1
      )
      B.kappa <- prepare_B_matrices(
        B.kappa, ncol(C),
        ncol(B.kappa) - 1
      )
      if (is.null(alpha)) {
        alpha <- 1 + d / 2
      } else {
        alpha <- rspde_check_user_input(alpha, "alpha", d / 2)
        alpha <- min(alpha, 10)
      }

      nu <- alpha - d / 2
    }

    if (is.null(graph) && is.null(mesh) && is.null(range_mesh)) {
      stop("You should either provide all the parameters, or you should provide one of the following: mesh, range_mesh or graph.")
    }
    if (!is.null(mesh) || !is.null(range_mesh)) {
      theta <- get.initial.values.rSPDE(
        mesh.range = range_mesh,
        dim = d, parameterization = parameterization,
        mesh = mesh, nu = nu, B.tau = B.tau,
        B.sigma = B.sigma, B.range = B.range,
        B.kappa = B.kappa, include.nu = FALSE
      )
    } else if (!is.null(graph)) {
      theta <- get.initial.values.rSPDE(
        graph.obj = graph,
        parameterization = parameterization, nu = nu,
        B.tau = B.tau, B.sigma = B.sigma, B.range = B.range,
        B.kappa = B.kappa, include.nu = FALSE
      )
    }

    new_theta <- c(1, theta)

    if (is.null(tau)) {
      tau <- c(exp(B.tau %*% new_theta))
    }
    if (is.null(kappa)) {
      kappa <- c(exp(B.kappa %*% new_theta))
    }
  } else {
    if (parameterization == "matern") {
      if (is.null(nu)) {
        if(d == 1){
          nu <- 0.75
        } else {
          nu <- 1
        }
      } else {
        nu <- rspde_check_user_input(nu, "nu", 0)
      }
      alpha <- nu + d / 2
    } else {
      if (is.null(alpha)) {
        alpha <- 1 + d / 2
      } else {
        alpha <- rspde_check_user_input(alpha, "alpha", d / 2)
        alpha <- min(alpha, 10)
      }

      nu <- alpha - d / 2
    }
  }

  if (!is.null(tau) && !is.null(kappa) &&
      is_constant_param(tau) && is_constant_param(kappa)) {
    if (parameterization == "spde") {
      return(matern.operators(
        kappa = as.numeric(kappa)[1],
        tau = as.numeric(tau)[1],
        alpha = alpha,
        G = G,
        C = C,
        d = d,
        mesh = mesh,
        graph = graph,
        range_mesh = range_mesh,
        loc_mesh = loc_mesh,
        m = m,
        type = type,
        parameterization = "spde",
        type_rational_approximation = type_rational_approximation
      ))
    }
    range <- sqrt(8 * nu) / as.numeric(kappa)[1]
    sigma <- sqrt(gamma(nu) / (as.numeric(tau)[1]^2 * as.numeric(kappa)[1]^(2 * nu) *
      (4 * pi)^(d / 2) * gamma(nu + d / 2)))
    return(matern.operators(
      range = range,
      sigma = sigma,
      nu = nu,
      G = G,
      C = C,
      d = d,
      mesh = mesh,
      graph = graph,
      range_mesh = range_mesh,
      loc_mesh = loc_mesh,
      m = m,
      type = type,
      parameterization = "matern",
      type_rational_approximation = type_rational_approximation
    ))
  }

  if (nu < 0) {
    stop("nu must be positive")
  }

  if (type == "operator") {
    beta <- (nu + d / 2) / 2
    kp <- Matrix::Diagonal(dim(C)[1], kappa^2)

    operators <- fractional.operators(
      L = G + C %*% kp,
      beta = beta,
      C = C,
      scale.factor = min(kappa)^2,
      m = m,
      tau = tau
    )
    output <- operators
    output$d <- d
    output$alpha <- alpha
    output$beta <- beta
    output$B.tau <- B.tau
    output$B.kappa <- B.kappa
    output$kappa <- kappa
    output$tau <- tau
    output$theta <- theta
    output$stationary <- FALSE
    output$type <- "Matern SPDE approximation"
    output$has_mesh <- has_mesh
    output$has_graph <- has_graph
    output$mesh <- mesh
    output$range_mesh <- range_mesh
    output$graph <- graph
    output$loc_mesh <- loc_mesh
    output$mesh_1d <- mesh_1d
    output$B.sigma <- B.sigma
    output$B.range <- B.range
    output$parameterization <- parameterization
    class(output) <- c("spde_matern_operator", class(output))
  } else {
    type_rational_approximation <- type_rational_approximation[[1]]
    m_alpha <- floor(alpha)
    m_order <- m_alpha + 1

    C <- Matrix::Diagonal(dim(C)[1], rowSums(C))
    Ci <- Matrix::Diagonal(dim(C)[1], 1 / rowSums(C))

    L <- Matrix::Diagonal(dim(C)[1], kappa^2 * diag(C))
    L <- L + G

    tau_matrix <- Matrix::Diagonal(dim(C)[1], tau)

    if (alpha %% 1 == 0) {
      Q <- L

      if (alpha > 1) {
        CinvL <- Matrix::Diagonal(dim(C)[1], 1 / diag(C)) %*% L

        for (k in 1:(alpha - 1)) {
          Q <- Q %*% CinvL
        }
      }

      Q <- tau_matrix %*% Q %*% tau_matrix
    } else {
      factor <- min(kappa^2)
      L <- L / factor

      if (m_alpha > 0) {
        CinvL <- Matrix::Diagonal(dim(C)[1], 1 / diag(C)) %*% L
      }

      mt <- get_rational_coefficients(m, type_rational_approximation)
      rspde.order <- m
      row_nu <- round(1000 * cut_decimals(alpha))
      r <- unlist(mt[row_nu, 2:(1 + rspde.order)])
      p <- unlist(mt[row_nu, (2 + rspde.order):(1 + 2 * rspde.order)])
      k_rat <- unlist(mt[row_nu, 2 + 2 * rspde.order])

      Q <- (L - p[1] * C) / r[1]
      if (m_alpha > 0) {
        for (count_m in 1:m_alpha) {
          Q <- Q %*% CinvL
        }
      }

      Q <- tau_matrix %*% Q %*% tau_matrix

      if (rspde.order > 1) {
        for (k_ind in 2:rspde.order) {
          Q_tmp <- (L - p[k_ind] * C) / r[k_ind]
          if (m_alpha > 0) {
            for (count_m in 1:m_alpha) {
              Q_tmp <- Q_tmp %*% CinvL
            }
          }
        }
        Q_tmp <- tau_matrix %*% Q_tmp %*% tau_matrix
        Q <- bdiag(Q, Q_tmp)
      }

      # K part

      if (m_alpha == 0) {
        Q_tmp <- Ci
      } else if (m_alpha == 1) {
        Q_tmp <- L
      } else {
        Q_tmp <- L
        for (count_m in 1:(m_alpha - 1)) {
          Q_tmp <- Q_tmp %*% CinvL
        }
      }

      Q_tmp <- tau_matrix %*% Q_tmp %*% tau_matrix
      Q_tmp <- Q_tmp / k_rat

      Q <- bdiag(Q, Q_tmp)

      Q <- factor^(alpha) * Q
    }


    ## output
    output <- list(
      C = C, G = G,
      alpha = alpha, nu = nu, kappa = kappa,
      B.tau = B.tau, B.kappa = B.kappa,
      tau = tau, theta = theta, m = m, d = d,
      Q = Q, sizeC = dim(C)[1],
      type_rational_approximation = type_rational_approximation,
      stationary = FALSE
    )
    output$type <- "Covariance-Based Matern SPDE Approximation"
    output$has_mesh <- has_mesh
    output$has_graph <- has_graph
    output$mesh <- mesh
    output$range_mesh <- range_mesh
    output$graph <- graph
    output$loc_mesh <- loc_mesh
    output$B.sigma <- B.sigma
    output$B.range <- B.range
    output$parameterization <- parameterization
    output$mesh_1d <- mesh_1d
    class(output) <- c("spde_matern_operator", "CBrSPDEobj")
  }

  return(output)
}




#' Rational approximations of stationary anisotropic Gaussian Matern random fields
#'
#' `matern2d.operators` is used for computing a rational SPDE approximation
#' of a stationary Gaussian random fields on \eqn{R^d} with a Matern covariance
#' function
#' \deqn{C(h) = \frac{\sigma^2}{2^{\nu-1}\Gamma(\nu)}(\sqrt{h^T H^{-1}h})^\nu K_\nu(\sqrt{h^T H^{-1}h})},
#' based on a SPDE representation of the form
#' \deqn{(I - \nabla\cdot(H\nabla))^{(\nu+1)/2}u = c\sigma W},
#' where $c>0$ is a constant. The matrix \eqn{H} is defined as
#' \deqn{\begin{bmatrix}
#' h_x^2 & h_xh_yh_{xy} \\
#' h_xh_yh_{xy} & h_y^2
#' \end{bmatrix}}
#'
#' @param hx Parameter in the H matrix.
#' @param hy Parameter in the H matrix.
#' @param hxy Parameter in the H matrix.
#' @param sigma standard deviation parameter.
#' @param nu Smoothness parameter.
#' @param mesh An `fmesher` mesh.
#' @param fem Optional precomputed FEM matrices.
#' @param m The order of the rational approximation, which needs to be a
#' positive integer. The default value is 1.
#' @param type_rational_approximation Which type of rational
#' approximation should be used? The current types are
#' "brasil", "chebfun" or "chebfunLB".
#' @param return_fem_matrices Should the FEM matrices be returned?
#' @return An object of type `CBrSPDEobj2d`
#'
#' @export
#' @seealso [fractional.operators()],
#' [spde.matern.operators()],
#' [matern.operators()]
#'
#' @examples
#' library(fmesher)
#' n_loc <- 2000
#' loc_2d_mesh <- matrix(runif(n_loc * 2), n_loc, 2)
#' mesh_2d <- fm_mesh_2d(loc = loc_2d_mesh, cutoff = 0.03, max.edge = c(0.1, 0.5))
#' op <- matern2d.operators(mesh = mesh_2d)
matern2d.operators <- function(hx = NULL,
                             hy = NULL,
                             hxy = NULL,
                             nu = NULL,
                             sigma = NULL,
                             mesh = NULL,
                             fem = NULL,
                             m = 1,
                             type_rational_approximation = c(
                                 "brasil",
                                 "chebfun",
                                 "chebfunLB"
                             ),
                             return_fem_matrices = FALSE) {

    if (is.null(mesh)) {
        stop("No mesh provided!")
    } else {
        if (!inherits(mesh, c("fm_mesh_2d"))) {
            stop("The mesh should be created using INLA or fmesher!")
        }

        d <- fmesher::fm_manifold_dim(mesh)
        if(d != 2) {
            stop("Only 2d domains supported")
        }
    }

    type_rational_approximation <- type_rational_approximation[[1]]

    if(is.null(nu)) {
        nu <- 1
    } else {
        nu <- rspde_check_user_input(nu, "nu", 0)
    }
    alpha <- nu + 1

    if(is.null(sigma)) {
        sigma <- 1
    } else {
        sigma <- rspde_check_user_input(sigma, "sigma", 0)
    }



    if (is.null(hx) || is.null(hy) || is.null(hxy)) {
        mesh.range <- max(c(diff(range(mesh$loc[,1])),
                            diff(range(mesh$loc[, 2])),
                            diff(range(mesh$loc[,3]))))

        hxy <- 0
        hx <- hy <- 0.2*mesh.range
    } else {
        hx <- rspde_check_user_input(hx, "hx", 0)
        hy <- rspde_check_user_input(hy, "hy", 0)
        hxy <- rspde_check_user_input(hxy, "hxy", -1)
        if(hxy>1) {
            stop("hxy must be in (-1,1)")
        }
    }

    tau <- sqrt(gamma(alpha-1) / (gamma(alpha) * 4 * pi * sigma^2))
    tau <- tau / sqrt(hx*hy*sqrt(1-hxy^2))

    m_alpha <- floor(alpha)
    m_order <- m_alpha + 1
    if(is.null(fem)) {
        fem <- rSPDE.fem2d(P = mesh$loc[,1:2], FV = mesh$graph$tv)
    }

    C <- Matrix::Diagonal(dim(fem$C)[1], rowSums(fem$C))
    Ci <- Matrix::Diagonal(dim(C)[1], 1 / rowSums(C))

    L <- C + hx^2*fem$Hxx + hy^2*fem$Hyy + hx*hy*hxy*(fem$Hxy + t(fem$Hxy))

    CiL <- Ci %*% L

    if (alpha %% 1 == 0) {
        Q <- L
        if (alpha > 1) {
            for (k in 1:(alpha - 1)) {
                Q <- Q %*% CiL
            }
        }

        Q <- tau^2 * Q
    } else if (m > 0) {
        Q <- CBrSPDE.L.precision(L = L, tau = tau,
                                 C = C,
                                 beta = alpha / 2,
                                 rspde.order = m,
                                 scale_factor = 1,
                                 only_fractional = FALSE,
                                 type_rational_approx = type_rational_approximation)

    } else {
        stop("m > 0 required")
    }


    ## output
    out <- list(hx = hx,
                hy = hy,
                hxy = hxy,
                nu = nu,
                alpha = alpha,
                sigma = sigma,
                m = m,
                Q = Q,
                type_rational_approximation = type_rational_approximation,
                stationary = TRUE,
                type = "Covariance-Based aniostropic Matern SPDE Approximation",
                mesh = mesh,
                fem = fem,
                d = 2,
                has_graph = FALSE)

    if(return_fem_matrices) {
        out$C <- C
        out$Ci <- Ci
        out$Hxx <- fem$Hxx
        out$Hyy <- fem$Hyy
        out$Hxy <- fem$Hxy
    }

    class(out) <- c("matern2d_operator", "CBrSPDEobj2d")

    return(out)
}

#' @noRd
CBrSPDE.L.precision <- function(L, tau = 1,
                        rspde.order,
                        C,
                        beta,
                        only_fractional = FALSE,
                        scale_factor = 1,
                        type_rational_approx = "chebfun") {

    Ci <- Matrix::Diagonal(dim(C)[1], 1 / rowSums(C))
    n_m <- rspde.order

    mt <- get_rational_coefficients(n_m, type_rational_approx)

    m_alpha <- floor(2 * beta)

    row_nu <- round(1000 * cut_decimals(2 * beta))
    r <- unlist(mt[row_nu, 2:(1 + rspde.order)])
    p <- unlist(mt[row_nu, (2 + rspde.order):(1 + 2 * rspde.order)])
    k <- unlist(mt[row_nu, 2 + 2 * rspde.order])

    CiL <- Ci %*% L
    M <- Matrix::Diagonal(dim(L)[1])
    if (!only_fractional && m_alpha > 0) {
      for(i in 1:m_alpha) {
          M <- M %*% CiL
      }
    }

    Q <- (L - p[1] * C) %*% M / r[1]
    if (length(r) > 1) {
        for (i in 2:length(r)) {
            Q <- bdiag(Q, (L - p[i] * C) %*% M / r[i])
        }
    }

    # add k_part into Q

    if(m_alpha > 0) {
        K <- C
        for(i in 1:m_alpha) {
            K <- K %*% CiL
        }
    } else {
        K <- Ci
    }

    K <- K / k

    Q <- bdiag(Q, K)

    Q <- Q * scale_factor^(2 * beta) * tau^2

    return(Q)
}

#' Projection matrix for model objects
#'
#' Generic for computing the projection matrix that links observation locations
#' to model mesh nodes.
#'
#' @param object A model object.
#' @param ... Additional arguments passed to methods.
#' @export
make_A <- function(object, ...) {
  UseMethod("make_A")
}

#' @export
#' @method make_A default
make_A.default <- function(object, ...) {
  if (is.function(object$make_A)) {
    return(object$make_A(...))
  }
  stop("make_A is not available for this object.")
}

#' @export
#' @method make_A matern_operator
make_A.matern_operator <- function(object, loc, ...) {
  if (!is.null(object$mesh)) {
    return(fm_basis(x = object$mesh, loc = loc))
  }
  if (!is.null(object$graph)) {
    return(object$graph$fem_basis(loc))
  }
  if (!is.null(object$loc_mesh) && object$d == 1) {
    mesh_1d <- object$mesh_1d
    if (is.null(mesh_1d)) {
      mesh_1d <- fm_mesh_1d(object$loc_mesh)
    }
    return(fm_basis(x = mesh_1d, loc = loc))
  }
  stop("Please, create the object using the mesh.")
}

#' @export
#' @method make_A spde_matern_operator
make_A.spde_matern_operator <- function(object, loc, ...) {
  if (!is.null(object$mesh)) {
    return(fm_basis(x = object$mesh, loc = loc))
  }
  if (!is.null(object$graph)) {
    return(object$graph$fem_basis(loc))
  }
  if (!is.null(object$loc_mesh) && object$d == 1) {
    return(rspde.make.A(mesh = object$loc_mesh, loc = loc))
  }
  stop("Please, create the object using the mesh.")
}

#' @export
#' @method make_A matern2d_operator
make_A.matern2d_operator <- function(object, loc, ...) {
  A <- fm_basis(x = object$mesh, loc = loc)
  if (object$alpha %% 1 == 0) {
    return(A)
  }
  return(kronecker(matrix(1, ncol = object$m + 1), A))
}

#' Covariance between mesh nodes and locations
#'
#' Generic for computing the covariance between mesh nodes and locations for
#' model objects.
#'
#' @param object A model object.
#' @param ... Additional arguments passed to methods.
#' @export
cov_function_mesh <- function(object, ...) {
  UseMethod("cov_function_mesh")
}

#' @export
#' @method cov_function_mesh default
cov_function_mesh.default <- function(object, ...) {
  if (is.function(object$cov_function_mesh)) {
    return(object$cov_function_mesh(...))
  }
  stop("cov_function_mesh is not available for this object.")
}

#' @export
#' @method cov_function_mesh matern_operator
cov_function_mesh.matern_operator <- function(object, p, direct = FALSE, ...) {
  if (inherits(object, "rSPDEobj")) {
    v <- make_A(object, loc = p)
    if (direct) {
      return(object$Pr %*% solve(object$Q, object$Pr %*% t(v)))
    }
    return(Sigma.mult(object, t(v)))
  }
  v <- t(make_A(object, loc = p))
  A <- Matrix::Diagonal(dim(object$C)[1])
  if (object$alpha %% 1 == 0) {
    return(A %*% solve(object$Q, v))
  }
  v_bar <- kronecker(matrix(1, nrow = object$m + 1), v)
  A_bar <- kronecker(matrix(1, ncol = object$m + 1), A)
  return(A_bar %*% solve(object$Q, v_bar))
}

#' @export
#' @method cov_function_mesh spde_matern_operator
cov_function_mesh.spde_matern_operator <- function(object, p, direct = FALSE, ...) {
  if (inherits(object, "rSPDEobj")) {
    v <- make_A(object, loc = p)
    if (direct) {
      return(object$Pr %*% solve(object$Q, object$Pr %*% t(v)))
    }
    return(Sigma.mult(object, t(v)))
  }
  v <- t(make_A(object, loc = p))
  A <- Matrix::Diagonal(dim(object$C)[1])
  if (object$alpha %% 1 == 0) {
    return(A %*% solve(object$Q, v))
  }
  v_bar <- kronecker(matrix(1, nrow = object$m + 1), v)
  A_bar <- kronecker(matrix(1, ncol = object$m + 1), A)
  return(A_bar %*% solve(object$Q, v_bar))
}

#' @export
#' @method cov_function_mesh matern2d_operator
cov_function_mesh.matern2d_operator <- function(object, p, ...) {
  v <- t(make_A(object, loc = p))
  A <- Matrix::Diagonal(dim(object$fem$C)[1])
  if (object$alpha %% 1 != 0) {
    A <- kronecker(matrix(1, ncol = object$m + 1), A)
  }
  return(A %*% solve(object$Q, v))
}

#' Covariance between mesh nodes
#'
#' Generic for computing the covariance between mesh nodes for model objects.
#'
#' @param object A model object.
#' @param ... Additional arguments passed to methods.
#' @export
covariance_mesh <- function(object, ...) {
  UseMethod("covariance_mesh")
}

#' @export
#' @method covariance_mesh default
covariance_mesh.default <- function(object, ...) {
  if (is.function(object$covariance_mesh)) {
    return(object$covariance_mesh(...))
  }
  stop("covariance_mesh is not available for this object.")
}

#' @export
#' @method covariance_mesh matern_operator
covariance_mesh.matern_operator <- function(object, ...) {
  if (inherits(object, "rSPDEobj")) {
    return(object$Pr %*% solve(object$Q, object$Pr))
  }
  A <- Matrix::Diagonal(dim(object$C)[1])
  if (object$alpha %% 1 == 0) {
    return(A %*% solve(object$Q, t(A)))
  }
  A_bar <- kronecker(matrix(1, ncol = object$m + 1), A)
  return(A_bar %*% solve(object$Q, t(A_bar)))
}

#' @export
#' @method covariance_mesh spde_matern_operator
covariance_mesh.spde_matern_operator <- function(object, ...) {
  if (inherits(object, "rSPDEobj")) {
    return(object$Pr %*% solve(object$Q, object$Pr))
  }
  A <- Matrix::Diagonal(dim(object$C)[1])
  if (object$alpha %% 1 == 0) {
    return(A %*% solve(object$Q, t(A)))
  }
  A_bar <- kronecker(matrix(1, ncol = object$m + 1), A)
  return(A_bar %*% solve(object$Q, t(A_bar)))
}

#' @export
#' @method covariance_mesh matern2d_operator
covariance_mesh.matern2d_operator <- function(object, ...) {
  C <- Matrix::Diagonal(dim(object$fem$C)[1], rowSums(object$fem$C))
  A <- Matrix::Diagonal(dim(C)[1])
  if (object$alpha %% 1 == 0) {
    return(A %*% solve(object$Q, t(A)))
  }
  A_bar <- kronecker(matrix(1, ncol = object$m + 1), A)
  return(A_bar %*% solve(object$Q, t(A_bar)))
}