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) {
        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
#' "chebfun", "brasil" or "chebfunLB".
#' @param fem_mesh_matrices A list containing FEM-related matrices.
#' The list should contain elements c0, g1, g2, g3, etc.
#' @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("chebfun",
                             "brasil", "chebfunLB"),
                             fem_mesh_matrices = NULL,
                             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!")
  }

  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!")
  }

  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 <- get_inla_mesh_dimension(inla_mesh = 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)){
        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(!is.null(mesh)){
      make_A <- function(loc){
        # return(INLA::inla.spde.make.A(mesh = mesh, loc = loc))
        # return(fmesher::fm_basis(x = mesh, loc = loc))
        return(fm_basis(x = mesh, loc = loc))
      }
  } else if(!is.null(graph)){
      make_A <- function(loc){
        return(graph$fem_basis(loc))
      }
  } else if(!is.null(loc_mesh) && d == 1){
    make_A <- function(loc){
          # return(rSPDE::rSPDE.A1d(x = loc_mesh, loc = loc))
          # return(fmesher::fm_basis(x = mesh_1d, loc = loc))
          return(fm_basis(x = mesh_1d, loc = loc))
          }   
  } else {
    make_A <- NULL
  }

  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$make_A <- make_A
    output$alpha <- 2*beta
    output$mesh <- mesh
    output$range_mesh <- range_mesh
    output$graph <- graph
    output$loc_mesh <- loc_mesh
    output$cov_function_mesh <- function(p, direct = FALSE){
      if(is.null(output$make_A)){
        stop("Please, create the object using the mesh.")
      }
      v <- output$make_A(loc = p)
      if(direct){
        return(output$Pr %*% solve(output$Q, output$Pr %*% t(v)))
      } else{
        return(Sigma.mult(output, t(v)))
      } 
    }
    output$covariance_mesh <- function(){
      return(output$Pr %*% solve(output$Q, output$Pr))
  }    
    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$make_A <- make_A
    out$mesh <- mesh
    out$range_mesh <- range_mesh
    out$graph <- graph
    out$loc_mesh <- loc_mesh    
    out$cov_function_mesh <- function(p){
      if(is.null(out$make_A)){
        stop("Please, create the object using the mesh.")
      }
      v <- t(out$make_A(loc = p))
      A <- Matrix::Diagonal(dim(C)[1])
      v_bar <- kronecker(matrix(1, nrow = m + 1), v)
      A_bar <- kronecker(matrix(1, ncol = m+1), A)
      return((A_bar) %*% solve(out$Q, v_bar))
    }

  out$covariance_mesh <- function(){
      A <- Matrix::Diagonal(dim(C)[1])
      A_bar <- kronecker(matrix(1, ncol = m+1), A)
      return((A_bar) %*% solve(out$Q, t(A_bar)))
  }    
    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 <- get_inla_mesh_dimension(inla_mesh = 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 <- get_inla_mesh_dimension(inla_mesh = 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 * nu) / 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`.
#' @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
#' "chebfun", "brasil" 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("chebfun",
                             "brasil", "chebfunLB")) {
  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!")
  }

  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 <- get_inla_mesh_dimension(inla_mesh = 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)){
        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
    # }
    
      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(any(dim(B.tau) != dim(B.kappa))){
        stop("B.tau and B.kappa must have the same dimensions!")
      }
      if(any(dim(B.sigma) != dim(B.range))){
        stop("B.sigma and B.range must have the same dimensions!")
      }


    if(parameterization == "matern"){
      if(is.null(nu)){
        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)){
        nu <- 1
      } else{
      nu <- rspde_check_user_input(nu, "nu" , 0)
    }    

      alpha <- nu + d/2

    } else{
      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)){
        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(mesh)){
      make_A <- function(loc){
        # return(INLA::inla.spde.make.A(mesh = mesh, loc = loc))
        # return(fmesher::fm_basis(x = mesh, loc=loc))
        return(fm_basis(x = mesh, loc=loc))
      }
  } else if(!is.null(graph)){
      make_A <- function(loc){
        return(graph$fem_basis(loc))
      }
  } else if(!is.null(loc_mesh) && d == 1){
    make_A <- function(loc){
          return(rspde.make.A(mesh = loc_mesh, loc = loc))
          }   
  } else {
    make_A <- NULL
  }
  
  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$make_A <- make_A
      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$cov_function_mesh <- function(p, direct = FALSE){
      if(is.null(output$make_A)){
        stop("Please, create the object using the mesh.")
      }
      v <- output$make_A(loc = p)
      if(direct){
        return(output$Pr %*% solve(output$Q, output$Pr %*% t(v)))
      } else{
        return(Sigma.mult(output, t(v)))
      } 
    }
    output$covariance_mesh <- function(){
      return(output$Pr %*% solve(output$Q, output$Pr))
  }    
  } 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$make_A <- make_A
  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$cov_function_mesh <- function(p){
      if(is.null(output$make_A)){
        stop("Please, create the object using the mesh.")
      }
      v <- t(output$make_A(loc = p))
      A <- Matrix::Diagonal(dim(C)[1])
      v_bar <- kronecker(matrix(1, nrow = m + 1), v)
      A_bar <- kronecker(matrix(1, ncol = m+1), A)
      return((A_bar) %*% solve(output$Q, v_bar))
    }
  class(output) <- "CBrSPDEobj"
  }

  output$covariance_mesh <- function(){
      A <- Matrix::Diagonal(dim(C)[1])
      A_bar <- kronecker(matrix(1, ncol = m+1), A)
      return((A_bar) %*% solve(output$Q, t(A_bar)))
  }

  return(output)
}