R/fractional.computations.R
In rSPDE: Rational Approximations of Fractional Stochastic Partial Differential Equations

Documented in construct.spde.matern.loglike precision precision.CBrSPDEobj precision.CBrSPDEobj2d predict.CBrSPDEobj predict.CBrSPDEobj2d predict.rSPDEobj rSPDE.construct.matern.loglike rSPDE.loglike rSPDE.matern.loglike simulate.CBrSPDEobj simulate.CBrSPDEobj2d simulate.rSPDEobj spde.matern.loglike update.CBrSPDEobj update.CBrSPDEobj2d update.rSPDEobj

###########################################
## fractional.computations.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/>.

#' @title Simulation of a fractional SPDE using a rational SPDE approximation
#'
#' @description The function samples a Gaussian random field based on a
#' pre-computed rational SPDE approximation.
#'
#' @param object The rational SPDE approximation, computed
#' using [fractional.operators()],
#' [matern.operators()], or [spde.matern.operators()].
#' @param nsim The number of simulations.
#' @param seed an object specifying if and how the random number generator should be initialized (‘seeded’).
#' @param ... Currently not used.
#'
#' @return A matrix with the `n` samples as columns.
#' @seealso [simulate.CBrSPDEobj()]
#' @export
#' @method simulate rSPDEobj
#'
#' @examples
#' # Sample a Gaussian Matern process on R using a rational approximation
#' 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
#' op <- matern.operators(
#'   range = range, sigma = sigma,
#'   nu = nu, loc_mesh = x, d = 1,
#'   parameterization = "matern"
#' )
#'
#' # Sample the model and plot the result
#' Y <- simulate(op)
#' plot(x, Y, type = "l", ylab = "u(x)", xlab = "x")
#'
simulate.rSPDEobj <- function(object,
                              nsim = 1,
                              seed = NULL,
                              ...) {
  if (!is.null(seed)) {
    set.seed(seed)
  }

  if (!inherits(object, "rSPDEobj")) {
    stop("input object is not of class rSPDEobj")
  }
  m <- dim(object$Q)[1]
  z <- rnorm(nsim * m)
  dim(z) <- c(m, nsim)
  x <- Qsqrt.solve(object, z)
  x <- Pr.mult(object, x)

  return(x)
}


#' @name update.CBrSPDEobj
#' @title Update parameters of CBrSPDEobj objects
#' @description Function to change the parameters of a CBrSPDEobj object
#' @param object The covariance-based rational SPDE approximation,
#' computed using [matern.operators()]
#' @param kappa If non-null, update the parameter kappa of the SPDE. Will be used if parameterization is 'spde'.
#' @param tau If non-null, update the parameter tau of the SPDE. Will be used if parameterization is 'spde'.
#' @param sigma If non-null, update the standard deviation of
#' the covariance function. Will be used if parameterization is 'matern'.
#' @param range If non-null, update the range parameter
#' of the covariance function. Will be used if parameterization is 'matern'.
#' @param theta For non-stationary models. If non-null, update the vector of parameters.
#' @param nu If non-null, update the shape parameter of the
#' covariance function. Will be used if parameterization is 'matern'.
#' @param alpha If non-null, update the fractional SPDE order parameter. Will be used if parameterization is 'spde'.
#' @param m If non-null, update the order of the rational
#' approximation, which needs to be a positive integer.
#' @param mesh An optional inla 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 parameterization If non-null, update the parameterization. Only works for stationary models.
#' @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 ... Currently not used.
#' @return It returns an object of class "CBrSPDEobj. This object contains the
#' same quantities listed in the output of [matern.operators()].
#' @method update CBrSPDEobj
#' @seealso [simulate.CBrSPDEobj()], [matern.operators()]
#' @export
#' @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
#' 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"
#' )
#' op_cov
#'
#' # Update the range parameter of the model:
#' op_cov <- update(op_cov, kappa = 20)
#' op_cov
#'
update.CBrSPDEobj <- function(object, nu = NULL, alpha = NULL,
                              kappa = NULL,
                              tau = NULL,
                              sigma = NULL,
                              range = NULL,
                              theta = NULL,
                              m = NULL,
                              mesh = NULL,
                              loc_mesh = NULL,
                              graph = NULL,
                              range_mesh = NULL,
                              compute_higher_order = object$higher_order,
                              parameterization = NULL,
                              type_rational_approximation =
                                object$type_rational_approximation,
                              return_block_list = object$return_block_list,
                              ...) {
  new_object <- object
  d <- object$d

  if (object$stationary) {
    fem_mesh_matrices <- object$fem_mesh_matrices

    if (is.null(nu) && !(object$higher_order) && compute_higher_order) {
      nu <- object$nu
    }

    new_object[["fem_mesh_matrices"]] <- fem_mesh_matrices

    if (is.null(parameterization)) {
      parameterization <- new_object$parameterization
    } else {
      parameterization <- parameterization[[1]]
      if (!parameterization %in% c("matern", "spde")) {
        stop("parameterization should be either 'matern' or 'spde'!")
      }
    }

    if (parameterization == "spde") {
      if (!is.null(kappa)) {
        new_object$kappa <- rspde_check_user_input(kappa, "kappa", 0)
        new_object$range <- NULL
        new_object$sigma <- NULL
      }

      if (!is.null(tau)) {
        new_object$tau <- rspde_check_user_input(tau, "tau", 0)
        new_object$sigma <- NULL
      }

      if (!is.null(alpha)) {
        alpha <- rspde_check_user_input(alpha, "alpha", d / 2)
        nu <- alpha - d / 2
        new_object$nu <- nu
        new_object$alpha <- alpha
      }
    } else if (parameterization == "matern") {
      if (!is.null(range)) {
        new_object$range <- rspde_check_user_input(range, "range", 0)
        new_object$kappa <- NULL
        new_object$tau <- NULL
      }

      if (!is.null(sigma)) {
        new_object$sigma <- rspde_check_user_input(sigma, "sigma", 0)
        new_object$tau <- NULL
      }
      if (!is.null(nu)) {
        new_object$nu <- rspde_check_user_input(nu, "nu")
      }
      alpha <- new_object$nu + d / 2
      new_object$alpha <- alpha
    }
    # else if(parameterization == "graph"){
    #   if (!is.null(kappa)) {
    #     new_object$kappa <- rspde_check_user_input(kappa, "kappa")
    #     new_object$range <- NULL
    #     new_object$tau <- NULL
    #   }

    #   if (!is.null(sigma)) {
    #     new_object$sigma <- rspde_check_user_input(sigma, "sigma")
    #     new_object$tau <- NULL
    #   }
    # }

    ## get parameters
    alpha <- new_object$alpha

    m_alpha <- floor(alpha)
    m_order <- m_alpha + 1

    if (compute_higher_order) {
      if (m_order + 1 > length(object$fem_mesh_matrices)) {
        old_m_order <- length(object$fem_mesh_matrices) - 1
        GCi <- object$GCi
        for (i in (old_m_order + 1):m_order) {
          fem_mesh_matrices[[paste0("g", i)]] <- GCi %*%
            fem_mesh_matrices[[paste0("g", i - 1)]]
        }
      }
    }


    if (!is.null(m)) {
      new_object$m <- as.integer(rspde_check_user_input(m, "m", 0))
    }

    if (is.null(mesh)) {
      mesh <- new_object[["mesh"]]
    }
    if (is.null(range_mesh)) {
      range_mesh <- new_object[["range_mesh"]]
    }
    if (is.null(loc_mesh)) {
      loc_mesh <- new_object[["loc_mesh"]]
    }
    if (is.null(graph)) {
      graph <- new_object$graph
    }

    if (parameterization == "spde") {
      new_object <- matern.operators(
        kappa = new_object$kappa,
        tau = new_object$tau,
        alpha = new_object$alpha,
        G = new_object$G,
        C = new_object$C,
        d = new_object$d,
        m = new_object$m,
        mesh = mesh,
        loc_mesh = loc_mesh,
        range_mesh = range_mesh,
        graph = graph,
        parameterization = parameterization,
        type = "covariance",
        return_block_list = return_block_list,
        type_rational_approximation = type_rational_approximation,
        # fem_mesh_matrices = new_object$fem_mesh_matrices,
        compute_logdet = new_object$compute_logdet
      )
    } else {
      new_object <- matern.operators(
        sigma = new_object$sigma,
        range = new_object$range,
        nu = new_object$nu,
        G = new_object$G,
        C = new_object$C,
        d = new_object$d,
        m = new_object$m,
        mesh = mesh,
        loc_mesh = loc_mesh,
        range_mesh = range_mesh,
        graph = graph,
        parameterization = parameterization,
        type = "covariance",
        return_block_list = return_block_list,
        type_rational_approximation = type_rational_approximation,
        # fem_mesh_matrices = new_object$fem_mesh_matrices,
        compute_logdet = new_object$compute_logdet
      )
    }
  } else {
    ## get parameters

    if (!is.null(tau)) {
      new_object$tau <- rspde_check_user_input(tau, "tau", 0)
    }

    if (!is.null(kappa)) {
      new_object$kappa <- rspde_check_user_input(kappa, "kappa", 0)
    }

    if (!is.null(theta)) {
      if (!is.numeric(theta)) {
        stop("theta must be numeric!")
      }
      new_object$theta <- theta
    }

    if (!is.null(m)) {
      new_object$m <- as.integer(rspde_check_user_input(m, "m", 0))
    }

    if (is.null(mesh)) {
      mesh <- new_object[["mesh"]]
    }
    if (is.null(range_mesh)) {
      range_mesh <- new_object[["range_mesh"]]
    }
    if (is.null(loc_mesh)) {
      loc_mesh <- new_object[["loc_mesh"]]
    }
    if (is.null(graph)) {
      graph <- new_object$graph
    }

    if (is.null(parameterization)) {
      parameterization <- new_object$parameterization
    } else {
      parameterization <- parameterization[[1]]
      if (!parameterization %in% c("matern", "spde")) {
        stop("parameterization should be either 'matern', 'spde' or 'graph'!")
      }
    }

    if (parameterization == "spde") {
      if (!is.null(alpha)) {
        alpha <- rspde_check_user_input(alpha, "alpha", d / 2)
        nu <- alpha - d / 2
        new_object$nu <- nu
        new_object$alpha <- alpha
      }
    } else if (parameterization == "matern") {
      if (!is.null(nu)) {
        new_object$nu <- rspde_check_user_input(nu, "nu")
      }
      alpha <- new_object$nu + d / 2
      new_object$alpha <- alpha
    }
    if (parameterization == "spde") {
      new_object <- spde.matern.operators(
        kappa = new_object$kappa,
        tau = new_object$tau,
        theta = new_object$theta,
        alpha = new_object$alpha,
        B.tau = new_object$B.tau,
        B.kappa = new_object$B.kappa,
        G = new_object$G,
        C = new_object$C,
        d = new_object$d,
        m = new_object$m,
        mesh = mesh,
        loc_mesh = loc_mesh,
        range_mesh = range_mesh,
        graph = graph,
        parameterization = parameterization,
        type = "covariance",
        type_rational_approximation = new_object$type_rational_approximation
      )
    } else {
      new_object <- spde.matern.operators(
        kappa = new_object$kappa,
        tau = new_object$tau,
        theta = new_object$theta,
        nu = new_object$nu,
        G = new_object$G,
        C = new_object$C,
        d = new_object$d,
        m = new_object$m,
        mesh = mesh,
        loc_mesh = loc_mesh,
        range_mesh = range_mesh,
        graph = graph,
        parameterization = parameterization,
        B.sigma = new_object$B.sigma,
        B.range = new_object$B.range,
        type = "covariance",
        type_rational_approximation = new_object$type_rational_approximation
      )
    }
  }
  return(new_object)
}

#' @name update.CBrSPDEobj2d
#' @title Update parameters of CBrSPDEobj2d objects
#' @description Function to change the parameters of a CBrSPDEobj object
#' @param object The covariance-based rational SPDE approximation,
#' computed using [matern2d.operators()]
#' @param hx If non-null, update the hx parameter.
#' @param hy If non-null, update the hy parameter.
#' @param hxy If non-null, update the hxy parameter.
#' @param sigma If non-null, update the standard deviation of
#' the covariance function. 
#' @param nu If non-null, update the shape parameter of the
#' covariance function. Will be used if parameterization is 'matern'.
#' @param m If non-null, update the order of the rational
#' approximation, which needs to be a positive integer.
#' @param ... Currently not used.
#' @return It returns an object of class "CBrSPDEobj2d. 
#' @method update CBrSPDEobj2d
#' @seealso [simulate.CBrSPDEobj2d()], [matern2d.operators()]
#' @export
#' @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)
#' op <- update(op, nu = 0.5)
update.CBrSPDEobj2d <- function(object, 
                                hx = NULL,
                                hy = NULL,
                                hxy = NULL,
                                sigma = NULL,
                                nu = NULL, 
                                m = NULL,
                                ...) {
    new_object <- object
    
    
    if (!is.null(nu)) {
        new_object$nu <- rspde_check_user_input(nu, "nu", 0)
    }
     
    if (!is.null(hx)) {
        new_object$hx <- rspde_check_user_input(hx, "hx", 0)
    }
    if (!is.null(hy)) {
        new_object$hy <- rspde_check_user_input(hy, "hy", 0)
    }
    if (!is.null(hxy)) {
        new_object$hxy <- rspde_check_user_input(hxy, "hxy", -1)
        if(new_object$hxy > 1) {
            stop("hxy must be in (-1,1)")
        }
    }
            
    if (!is.null(sigma)) {
        new_object$sigma <- rspde_check_user_input(sigma, "sigma", 0)
    }
            

    if (!is.null(m)) {
        new_object$m <- as.integer(rspde_check_user_input(m, "m", 0))
    }
        
    new_object <- matern2d.operators(hx = new_object$hx,
                                     hy = new_object$hy,
                                     hxy = new_object$hxy,
                                     nu = new_object$nu,
                                     sigma = new_object$sigma,
                                     fem = new_object$fem,
                                     m = new_object$m,
                                     mesh = new_object$mesh)
       
    return(new_object)
}


#' @name update.rSPDEobj
#' @title Update parameters of rSPDEobj objects
#' @description Function to change the parameters of a rSPDEobj object
#' @param object The operator-based rational SPDE approximation,
#' computed using [matern.operators()] with `type="operator"`
#' @param kappa If non-null, update the range parameter
#' of the covariance function.
#' @param tau If non-null, update the parameter tau.
#' @param sigma If non-null, update the standard deviation of
#' the covariance function.
#' @param range If non-null, update the range parameter
#' of the covariance function.
#' @param theta If non-null, update the parameter theta, that connects
#' tau and kappa to the model matrices.
#' @param nu If non-null, update the shape parameter
#' of the covariance function.
#' @param alpha If non-null, update the fractional order.
#' @param m If non-null, update the order of the rational
#' approximation, which needs to be a positive integer.
#' @param mesh An optional inla 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 parameterization If non-null, update the parameterization. Only works for stationary models.
#' @param ... Currently not used.
#' @return It returns an object of class "rSPDEobj. This object contains the
#' same quantities listed in the output of [matern.operators()].
#' @method update rSPDEobj
#' @seealso [simulate.rSPDEobj()], [matern.operators()]
#' @export
#' @examples
#' # 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(
#'   loc_mesh = x, nu = nu,
#'   range = range, sigma = sigma, d = 1, m = 2, type = "operator",
#'   parameterization = "matern"
#' )
#' op
#'
#' # Update the range parameter of the model:
#' op <- update(op, kappa = 20)
#' op
#'
update.rSPDEobj <- function(object, nu = NULL,
                            alpha = NULL,
                            kappa = NULL,
                            sigma = NULL,
                            range = NULL,
                            tau = NULL,
                            theta = NULL,
                            m = NULL,
                            mesh = NULL,
                            loc_mesh = NULL,
                            graph = NULL,
                            range_mesh = NULL,
                            parameterization = NULL, ...) {
  new_object <- object

  ## get parameters

  d <- object$d

  if (!is.null(m)) {
    new_object$m <- as.integer(rspde_check_user_input(m, "m", 1))
  }

  if (!is.null(theta)) {
    new_object$theta <- rspde_check_user_input(theta, "theta")
  }

  if (new_object$stationary) {
    if (is.null(parameterization)) {
      parameterization <- new_object$parameterization
    } else {
      parameterization <- parameterization[[1]]
      if (!parameterization %in% c("matern", "spde", "graph")) {
        stop("parameterization should be either 'matern', 'spde' or 'graph'!")
      }
    }

    if (parameterization == "spde") {
      if (!is.null(kappa)) {
        new_object$kappa <- rspde_check_user_input(kappa, "kappa", 0)
        new_object$range <- NULL
        new_object$sigma <- NULL
      }

      if (!is.null(tau)) {
        new_object$tau <- rspde_check_user_input(tau, "tau", 0)
        new_object$sigma <- NULL
      }

      if (!is.null(alpha)) {
        alpha <- rspde_check_user_input(alpha, "alpha", d / 2)
        nu <- alpha - d / 2
        new_object$nu <- nu
        new_object$alpha <- alpha
      }
    } else if (parameterization == "matern") {
      if (!is.null(range)) {
        new_object$range <- rspde_check_user_input(range, "range", 0)
        new_object$kappa <- NULL
        new_object$tau <- NULL
      }

      if (!is.null(sigma)) {
        new_object$sigma <- rspde_check_user_input(sigma, "sigma", 0)
        new_object$tau <- NULL
      }
      if (!is.null(nu)) {
        new_object$nu <- rspde_check_user_input(nu, "nu")
      }
      alpha <- new_object$nu + d / 2
      new_object$alpha <- alpha
    }

    if (!is.null(m)) {
      new_object$m <- as.integer(rspde_check_user_input(m, "m", 0))
    }

    if (is.null(mesh)) {
      mesh <- new_object[["mesh"]]
    }
    if (is.null(range_mesh)) {
      range_mesh <- new_object[["range_mesh"]]
    }
    if (is.null(loc_mesh)) {
      loc_mesh <- new_object[["loc_mesh"]]
    }
    if (is.null(graph)) {
      graph <- new_object$graph
    }

    if (parameterization == "spde") {
      new_object <- matern.operators(
        kappa = new_object$kappa,
        tau = new_object$tau,
        alpha = new_object$alpha,
        G = new_object$G,
        C = new_object$C,
        d = new_object$d,
        m = new_object$m,
        mesh = mesh,
        loc_mesh = loc_mesh,
        range_mesh = range_mesh,
        graph = graph,
        parameterization = parameterization,
        type = "operator"
      )
    } else {
      new_object <- matern.operators(
        sigma = new_object$sigma,
        range = new_object$range,
        nu = new_object$nu,
        G = new_object$G,
        C = new_object$C,
        d = new_object$d,
        m = new_object$m,
        mesh = mesh,
        loc_mesh = loc_mesh,
        range_mesh = range_mesh,
        graph = graph,
        parameterization = parameterization,
        type = "operator"
      )
    }
  } else {
    if (!is.null(tau)) {
      new_object$tau <- rspde_check_user_input(tau, "tau", 0)
    }

    if (!is.null(kappa)) {
      new_object$kappa <- rspde_check_user_input(kappa, "kappa", 0)
    }

    if (!is.null(theta)) {
      if (!is.numeric(theta)) {
        stop("theta must be numeric!")
      }
      new_object$theta <- theta
    }

    if (!is.null(m)) {
      new_object$m <- as.integer(rspde_check_user_input(m, "m", 0))
    }

    if (is.null(mesh)) {
      mesh <- new_object[["mesh"]]
    }
    if (is.null(range_mesh)) {
      range_mesh <- new_object[["range_mesh"]]
    }
    if (is.null(loc_mesh)) {
      loc_mesh <- new_object[["loc_mesh"]]
    }
    if (is.null(graph)) {
      graph <- new_object$graph
    }

    if (is.null(parameterization)) {
      parameterization <- new_object$parameterization
    } else {
      parameterization <- parameterization[[1]]
      if (!parameterization %in% c("matern", "spde", "graph")) {
        stop("parameterization should be either 'matern', 'spde' or 'graph'!")
      }
    }

    if (parameterization == "spde") {
      if (!is.null(alpha)) {
        alpha <- rspde_check_user_input(alpha, "alpha", d / 2)
        nu <- alpha - d / 2
        new_object$nu <- nu
        new_object$alpha <- alpha
      }
    } else if (parameterization == "matern") {
      if (!is.null(nu)) {
        new_object$nu <- rspde_check_user_input(nu, "nu")
      }
      alpha <- new_object$nu + d / 2
      new_object$alpha <- alpha
    }

    if (parameterization == "spde") {
      new_object <- spde.matern.operators(
        kappa = new_object$kappa,
        tau = new_object$tau,
        theta = new_object$theta,
        alpha = new_object$alpha,
        B.tau = new_object$B.tau,
        B.kappa = new_object$B.kappa,
        G = new_object$G,
        C = new_object$C,
        d = new_object$d,
        m = new_object$m,
        mesh = mesh,
        loc_mesh = loc_mesh,
        range_mesh = range_mesh,
        parameterization = parameterization,
        graph = graph,
        type = "operator"
      )
    } else {
      new_object <- spde.matern.operators(
        kappa = new_object$kappa,
        tau = new_object$tau,
        theta = new_object$theta,
        nu = new_object$nu,
        B.range = new_object$B.range,
        B.sigma = new_object$B.sigma,
        G = new_object$G,
        C = new_object$C,
        d = new_object$d,
        m = new_object$m,
        mesh = mesh,
        loc_mesh = loc_mesh,
        range_mesh = range_mesh,
        parameterization = parameterization,
        graph = graph,
        type = "operator"
      )
    }
  }



  return(new_object)
}


#' @name simulate.CBrSPDEobj
#' @title Simulation of a fractional SPDE using the
#' covariance-based rational SPDE approximation
#' @description The function samples a Gaussian random field based using the
#' covariance-based rational SPDE approximation.
#' @param object The covariance-based rational SPDE approximation,
#' computed using [matern.operators()]
#' @param nsim The number of simulations.
#' @param seed An object specifying if and how the random number generator should be initialized (‘seeded’).
#' @param kappa If non-null, update the range parameter
#' of the covariance function.
#' @param tau If non-null, update the parameter tau.
#' @param sigma If non-null, update the standard deviation of
#' the covariance function.
#' @param range If non-null, update the range parameter
#' of the covariance function.
#' @param theta For non-stationary models. If non-null, update the vector of parameters.
#' @param nu If non-null, update the shape parameter of the
#' covariance function.
#' @param m If non-null, update the order of the rational
#' approximation, which needs to be a positive integer.
#' @param ... Currently not used.
#' @return A matrix with the `n` samples as columns.
#' @method simulate CBrSPDEobj
#' @export
#' @examples
#' # Sample a Gaussian Matern process on R using a rational approximation
#' 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
#' tau <- sqrt(gamma(nu) / (sigma^2 * kappa^(2 * nu) *
#'   (4 * pi)^(1 / 2) * gamma(nu + 1 / 2)))
#' op_cov <- matern.operators(
#'   loc_mesh = x, nu = nu,
#'   range = range, sigma = sigma, d = 1, m = 2,
#'   parameterization = "matern"
#' )
#'
#' # Sample the model and plot the result
#' Y <- simulate(op_cov)
#' plot(x, Y, type = "l", ylab = "u(x)", xlab = "x")
#'
simulate.CBrSPDEobj <- function(object, nsim = 1,
                                seed = NULL,
                                nu = NULL,
                                kappa = NULL,
                                sigma = NULL,
                                range = NULL,
                                tau = NULL,
                                theta = NULL,
                                m = NULL,
                                ...) {
  if (!is.null(seed)) {
    set.seed(seed)
  }
  d <- object$d
  nu_temp <- ifelse(is.null(nu), object$nu, nu)
  alpha <- nu_temp + d / 2


  ## simulation
  if ((alpha %% 1 == 0) && object$stationary) { # simulation in integer case
    object <- update.CBrSPDEobj(
      object = object,
      nu = nu,
      kappa = kappa,
      sigma = sigma,
      tau = tau,
      range = range,
      m = m,
      parameterization = object$parameterization,
      compute_higher_order = TRUE
    )


    kappa <- object$kappa
    nu <- object$nu
    tau <- object$tau
    m <- object$m

    alpha <- nu + d / 2

    fem_mesh_matrices <- object$fem_mesh_matrices

    L <- fem_mesh_matrices[["c0"]] + fem_mesh_matrices[["g1"]] / kappa^2
    sizeL <- dim(L)[1]


    # Q <- rspde.matern.precision.integer(
    #   kappa = kappa, nu = nu, tau = tau,
    #   dim = d,
    #   fem_mesh_matrices = fem_mesh_matrices
    # )

    Q <- object$Q

    Z <- rnorm(sizeL * nsim)
    dim(Z) <- c(sizeL, nsim)


    # LQ <- Matrix::Cholesky(forceSymmetric(Q), LDL = FALSE)
    # X <- solve(LQ, Z, system = "Lt")
    # X <- solve(LQ, X, system = "Pt")
    LQ <- chol(forceSymmetric(Q))
    X <- solve(LQ, Z)
  } else {
    object <- update.CBrSPDEobj(
      object = object,
      nu = nu,
      kappa = kappa,
      sigma = sigma,
      m = m,
      range = range,
      tau = tau,
      theta = theta,
      parameterization = object$parameterization,
    )

    m <- object$m
    Q <- object$Q

    Z <- rnorm(dim(Q)[1] * nsim)
    dim(Z) <- c(dim(Q)[1], nsim)

    # LQ <- Matrix::Cholesky(forceSymmetric(Q), LDL = FALSE)
    # X <- solve(LQ, Z, system = "Lt")
    # X <- solve(LQ, X, system = "Pt")
    LQ <- chol(forceSymmetric(Q))
    X <- solve(LQ, Z)
    if (alpha %% 1 == 0) {
      A <- Diagonal(dim(Q)[1])
      Abar <- A
    } else {
      A <- Diagonal(dim(Q)[1] / (m + 1))
      Abar <- kronecker(matrix(1, 1, m + 1), A)
    }

    X <- Abar %*% X
  }
  return(X)
}

#' @name simulate.CBrSPDEobj2d
#' @title Simulation of a fractional SPDE using the
#' covariance-based rational SPDE approximation
#' @description The function samples a Gaussian random field based using the
#' covariance-based rational SPDE approximation.
#' @param object The covariance-based rational SPDE approximation,
#' computed using [matern2d.operators()]
#' @param nsim The number of simulations.
#' @param seed An object specifying if and how the random number generator should be initialized (‘seeded’).
#' @param hx If non-null, update the hx parameter.
#' @param hy If non-null, update the hy parameter.
#' @param hxy If non-null, update the hxy parameter.
#' @param sigma If non-null, update the standard deviation of
#' the covariance function.
#' @param nu If non-null, update the shape parameter of the
#' covariance function.
#' @param m If non-null, update the order of the rational
#' approximation, which needs to be a positive integer.
#' @param ... Currently not used.
#' @return A matrix with the `n` samples as columns.
#' @method simulate CBrSPDEobj2d
#' @export
#' @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, sigma = 1, nu = 1, hx = 0.1, hy = 0.1, hxy = 0)
#' u <- simulate(op)
simulate.CBrSPDEobj2d <- function(object, 
                                  nsim = 1,
                                  seed = NULL,
                                  nu = NULL,
                                  hx = NULL,
                                  hy = NULL,
                                  hxy = NULL,
                                  sigma = NULL,
                                  m = NULL,
                                  ...) {
    if (!is.null(seed)) {
        set.seed(seed)
    }
    
    nu <- ifelse(is.null(nu), object$nu, nu)
    alpha <- nu + 1
    
    
    object <- update.CBrSPDEobj2d(object = object,
                                  nu = nu,
                                  hx = hx,
                                  hy = hy,
                                  hxy = hxy,
                                  sigma = sigma,
                                  m = m)
    
    Q <- object$Q
    sizeQ <- dim(Q)[1]
    
    Z <- rnorm(sizeQ * nsim)
    dim(Z) <- c(sizeQ, nsim)
    
    LQ <- chol(forceSymmetric(Q))
    X <- solve(LQ, Z)
    
    if(object$alpha %% 1 != 0) {
        A <- Diagonal(dim(Q)[1] / (object$m + 1))
        Abar <- kronecker(matrix(1, 1, object$m + 1), A)
        X <- Abar %*% X
    }
    
    return(X)
}


#' Prediction of a fractional SPDE using a rational SPDE approximation
#'
#' The function is used for computing kriging predictions based on data
#' \eqn{Y_i = u(s_i) + \epsilon_i},
#' where \eqn{\epsilon}{\epsilon} is mean-zero Gaussian measurement noise
#' and \eqn{u(s)}{u(s)} is defined by
#' a fractional SPDE \eqn{L^\beta u(s) = W}{L^\beta u(s) = W}, where
#' \eqn{W}{W} is Gaussian white noise.
#'
#' @param object The rational SPDE approximation, computed using
#' [fractional.operators()],
#' [matern.operators()], or [spde.matern.operators()].
#' @param A A matrix linking the measurement locations to the basis of the
#' FEM approximation of the latent model.
#' @param Aprd A matrix linking the prediction locations to the basis of the
#' FEM approximation of the latent model.
#' @param Y A vector with the observed data, can also be a matrix where the
#' columns are observations
#' of independent replicates of \eqn{u}.
#' @param sigma.e The standard deviation of the Gaussian measurement noise.
#' Put to zero if the model
#' does not have measurement noise.
#' @param compute.variances Set to also TRUE to compute the kriging variances.
#' @param posterior_samples If `TRUE`, posterior samples will be returned.
#' @param n_samples Number of samples to be returned. Will only be used if `sampling` is `TRUE`.
#' @param only_latent Should the posterior samples be only given to the latent model?
#' @param ... further arguments passed to or from other methods.
#'
#' @return A list with elements
#' \item{mean }{The kriging predictor (the posterior mean of u|Y).}
#' \item{variance }{The posterior variances (if computed).}
#' \item{samples }{A matrix containing the samples if `sampling` is `TRUE`.}
#' @export
#' @method predict rSPDEobj
#'
#' @examples
#' # Sample a Gaussian Matern process on R using a rational approximation
#' kappa <- 10
#' sigma <- 1
#' nu <- 0.8
#' sigma.e <- 0.3
#' 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
#' op <- matern.operators(
#'   range = range, sigma = sigma,
#'   nu = nu, loc_mesh = x, d = 1,
#'   parameterization = "matern"
#' )
#'
#' # Sample the model
#' u <- simulate(op)
#'
#' # Create some data
#' obs.loc <- runif(n = 10, min = 0, max = 1)
#' A <- rSPDE.A1d(x, obs.loc)
#' Y <- as.vector(A %*% u + sigma.e * rnorm(10))
#'
#' # compute kriging predictions at the FEM grid
#' A.krig <- rSPDE.A1d(x, x)
#' u.krig <- predict(op,
#'   A = A, Aprd = A.krig, Y = Y, sigma.e = sigma.e,
#'   compute.variances = TRUE
#' )
#'
#' plot(obs.loc, Y,
#'   ylab = "u(x)", xlab = "x", main = "Data and prediction",
#'   ylim = c(
#'     min(u.krig$mean - 2 * sqrt(u.krig$variance)),
#'     max(u.krig$mean + 2 * sqrt(u.krig$variance))
#'   )
#' )
#' lines(x, u.krig$mean)
#' lines(x, u.krig$mean + 2 * sqrt(u.krig$variance), col = 2)
#' lines(x, u.krig$mean - 2 * sqrt(u.krig$variance), col = 2)
predict.rSPDEobj <- function(object,
                             A,
                             Aprd,
                             Y,
                             sigma.e,
                             compute.variances = FALSE,
                             posterior_samples = FALSE,
                             n_samples = 100,
                             only_latent = FALSE,
                             ...) {
  Y <- as.matrix(Y)
  if (dim(Y)[1] != dim(A)[1]) {
    stop("the dimensions of A does not match the number of observations")
  }

  out <- list()
  if (length(sigma.e) == 1) {
    if (sigma.e < 0) {
      stop("sigma.e must be non-negative")
    } else if (sigma.e > 0) {
      A <- A %*% object$Pr
      AA <- Aprd %*% object$Pr
      Qhat <- object$Q + (t(A) %*% A) / sigma.e^2
      out$mean <- as.matrix(AA %*% solve(Qhat, t(A) %*% Y / sigma.e^2))
      if (compute.variances) {
        out$variance <- diag(AA %*% solve(Qhat, t(AA)))
      }
    } else { # no nugget
      Ahat <- A %*% object$Pr
      QiAt <- solve(object$Q, t(Ahat))
      AQiA <- Ahat %*% QiAt
      xhat <- solve(object$Q, t(Ahat) %*% solve(AQiA, Y))
      out$mean <- as.vector(Aprd %*% xhat)
      if (compute.variances) {
        AA <- Aprd %*% object$Pr
        M <- object$Q - QiAt %*% solve(AQiA, t(QiAt))
        out$variance <- diag(AA %*% M %*% t(AA))
      }
    }
  } else if (dim(Y)[1] == length(sigma.e)) {
    Q.e <- Diagonal(length(sigma.e), 1 / sigma.e^2)
    A <- A %*% object$Pr
    AA <- Aprd %*% object$Pr
    Qhat <- object$Q + t(A) %*% Q.e %*% A
    out$mean <- as.matrix(AA %*% solve(Qhat, t(A) %*% Q.e %*% Y))
    if (compute.variances) {
      out$variance <- diag(AA %*% solve(Qhat, t(AA)))
    }
  }

  if (posterior_samples) {
    if (sigma.e > 0) {
      post_cov <- AA %*% solve(Qhat, t(AA))
    } else {
      AA <- Aprd %*% object$Pr
      M <- object$Q - QiAt %*% solve(AQiA, t(QiAt))
      post_cov <- AA %*% M %*% t(AA)
    }
    Y_tmp <- as.matrix(Y)
    mean_tmp <- as.matrix(out$mean)
    out$samples <- lapply(1:ncol(Y_tmp), function(i) {
      Z <- rnorm(dim(post_cov)[1] * n_samples)
      dim(Z) <- c(dim(post_cov)[1], n_samples)
      LQ <- chol(forceSymmetric(post_cov))
      X <- LQ %*% Z
      X <- X + mean_tmp[, i]
      if (!only_latent) {
        X <- X + matrix(rnorm(n_samples * dim(Aprd)[1], sd = sigma.e), nrow = dim(Aprd)[1])
      }
      return(X)
    })
  }

  return(out)
}

#' Object-based log-likelihood function for latent Gaussian
#' fractional SPDE model
#'
#' This function evaluates the log-likelihood function for a
#' fractional SPDE model
#' \eqn{L^\beta u(s) = W}{L^\beta u(s) = W} that is observed under
#' Gaussian measurement noise:
#' \eqn{Y_i = u(s_i) + \epsilon_i}{Y_i = x(s_i) + \epsilon_i},
#' where \eqn{\epsilon_i}{\epsilon_i}
#' are iid mean-zero Gaussian variables and \eqn{x(s) =
#' \mu(s) + u(s)}{x(s) = \mu(s) + u(s)}, where
#' \eqn{\mu(s)}{\mu(s)} is the expectation vector of the latent field.
#'
#' @param obj The rational SPDE approximation, computed using
#' [fractional.operators()],
#' [matern.operators()], or [spde.matern.operators()].
#' @param Y The observations, either a vector or a matrix where
#' the columns correspond to independent replicates of observations.
#' @param A An observation matrix that links the measurement location
#' to the finite element basis.
#' @param sigma.e The standard deviation of the measurement noise.
#' @param mu Expectation vector of the latent field (default = 0).
#' @return The log-likelihood value.
#' @export
#' @note This example below shows how the function can be used to evaluate the likelihood of a latent Matern model.
#' @seealso [spde.matern.loglike()]
#'
#' @examples
#' # Sample a Gaussian Matern process on R using a rational approximation
#' kappa <- 10
#' sigma <- 1
#' nu <- 0.8
#' sigma.e <- 0.3
#' range <- 0.2
#'
#' # 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
#' op <- matern.operators(
#'   range = range, sigma = sigma, nu = nu,
#'   loc_mesh = x, d = 1,
#'   type = "operator", parameterization = "matern"
#' )
#'
#' # Sample the model
#' u <- simulate(op)
#'
#' # Create some data
#' obs.loc <- runif(n = 10, min = 0, max = 1)
#' A <- rSPDE.A1d(x, obs.loc)
#' Y <- as.vector(A %*% u + sigma.e * rnorm(10))
#'
#' # compute log-likelihood of the data
#' lik1 <- rSPDE.loglike(op, Y, A, sigma.e)
#' cat(lik1)
rSPDE.loglike <- function(obj,
                          Y,
                          A,
                          sigma.e,
                          mu = 0) {
  Y <- as.matrix(Y)
  if (length(dim(Y)) == 2) {
    n.rep <- dim(Y)[2]
    n <- dim(Y)[1]
  } else {
    n.rep <- 1
    if (length(dim(Y)) == 1) {
      n <- dim(Y)[1]
    } else {
      n <- length(Y)
    }
  }
  if (length(sigma.e) == 1) {
    Q.e <- Diagonal(n) / sigma.e^2
    nugget <- rep(sigma.e^2, n)
  } else {
    if (length(sigma.e) != n) {
      stop("the length of sigma.e does not match the number of observations")
    }
    Q.e <- Diagonal(length(sigma.e), 1 / sigma.e^2)
    nugget <- sigma.e^2
  }



  R <- Matrix::Cholesky(obj$Pl)

  prior.ld <- 4 * c(determinant(R, logarithm = TRUE, sqrt = TRUE)$modulus) -
    sum(log(diag(obj$C)))


  A <- A %*% obj$Pr
  Q.post <- obj$Q + t(A) %*% Q.e %*% A


  R.post <- Matrix::Cholesky(Q.post)

  posterior.ld <- 2 * c(determinant(R.post, logarithm = TRUE, sqrt = TRUE)$modulus)

  AtY <- t(A) %*% Q.e %*% Y


  mu.post <- mu + solve(R.post, AtY, system = "A")


  lik <- n.rep * (prior.ld - posterior.ld - dim(A)[1] *
    log(2 * pi) - sum(log(nugget))) / 2

  if (n.rep > 1) {
    lik <- lik - 0.5 * sum(colSums((mu.post - mu) * (obj$Q %*% (mu.post - mu))))
    v <- Q.e %*% (Y - A %*% mu.post)
    lik <- lik - 0.5 * sum(colSums((Y - A %*% mu.post) * v))
  } else {
    lik <- lik - 0.5 * (t(mu.post - mu) %*% obj$Q %*% (mu.post - mu) +
      t(Y - A %*% mu.post) %*% Q.e %*% (Y - A %*% mu.post))
  }
  return(as.double(lik))
}


#' @name rSPDE.matern.loglike
#' @title Object-based log-likelihood function for latent Gaussian fractional
#' SPDE model using the rational approximations
#' @description This function evaluates the log-likelihood function for a
#' Gaussian process with a Matern covariance
#' function, that is observed under Gaussian measurement noise:
#' \eqn{Y_i = u(s_i) + \epsilon_i}{Y_i = u(s_i) + \epsilon_i}, where
#' \eqn{\epsilon_i}{\epsilon_i} are
#' iid mean-zero Gaussian variables. The latent model is approximated using
#' the a rational approximation
#' of the fractional SPDE model corresponding to the Gaussian process.
#' @param object The rational SPDE approximation,
#' computed using [matern.operators()]
#' @param Y The observations, either a vector or a matrix where
#' the columns correspond to independent replicates of observations.
#' @param A An observation matrix that links the measurement location to the
#' finite element basis.
#' @param sigma.e The standard deviation of the measurement noise.
#' @param mu Expectation vector of the latent field (default = 0).
#' @param kappa If non-null, update the range parameter of the covariance
#' function.
#' @param tau If non-null, update the parameter tau.
#' @param sigma If non-null, update the standard deviation of
#' the covariance function.
#' @param range If non-null, update the range parameter
#' of the covariance function.
#' @param nu If non-null, update the shape parameter of the covariance
#' function.
#' @param m If non-null, update the order of the rational approximation,
#' which needs to be a positive integer.
#' @return The log-likelihood value.
#' @export
#' @seealso [matern.operators()], [predict.CBrSPDEobj()]
#' @examples
#' # this example illustrates how the function can be used for maximum likelihood estimation
#'
#' set.seed(123)
#' # Sample a Gaussian Matern process on R using a rational approximation
#' nu <- 0.8
#' kappa <- 5
#' sigma <- 1
#' sigma.e <- 0.1
#' n.rep <- 10
#' n.obs <- 100
#' n.x <- 51
#' range <- 0.2
#'
#' # create mass and stiffness matrices for a FEM discretization
#' x <- seq(from = 0, to = 1, length.out = n.x)
#' fem <- rSPDE.fem1d(x)
#'
#' tau <- sqrt(gamma(nu) / (sigma^2 * kappa^(2 * nu) *
#'   (4 * pi)^(1 / 2) * gamma(nu + 1 / 2)))
#'
#' # Compute the covariance-based rational approximation
#' op_cov <- matern.operators(
#'   loc_mesh = x, nu = nu,
#'   range = range, sigma = sigma, d = 1, m = 2,
#'   parameterization = "matern"
#' )
#'
#' # Sample the model
#' u <- simulate(op_cov, n.rep)
#'
#' # Create some data
#' obs.loc <- runif(n = n.obs, min = 0, max = 1)
#' A <- rSPDE.A1d(x, obs.loc)
#' noise <- rnorm(n.obs * n.rep)
#' dim(noise) <- c(n.obs, n.rep)
#' Y <- as.matrix(A %*% u + sigma.e * noise)
#'
#' # Define the negative likelihood function for optimization
#' # using CBrSPDE.matern.loglike
#'
#' # Notice that we are also using sigma instead of tau, so it can be compared
#' # to matern.loglike()
#' mlik_cov <- function(theta, Y, A, op_cov) {
#'   kappa <- exp(theta[1])
#'   sigma <- exp(theta[2])
#'   nu <- exp(theta[3])
#'   return(-rSPDE.matern.loglike(
#'     object = op_cov, Y = Y,
#'     A = A, kappa = kappa, sigma = sigma,
#'     nu = nu, sigma.e = exp(theta[4])
#'   ))
#' }
#'
#' # The parameters can now be estimated by minimizing mlik with optim
#' \donttest{
#' # Choose some reasonable starting values depending on the size of the domain
#' theta0 <- log(c(sqrt(8), 1 / sqrt(var(c(Y))), 0.9, 0.01))
#'
#' # run estimation and display the results
#' theta <- optim(theta0, mlik_cov,
#'   Y = Y, A = A, op_cov = op_cov,
#'   method = "L-BFGS-B"
#' )
#'
#' print(data.frame(
#'   range = c(range, exp(theta$par[1])), sigma = c(sigma, exp(theta$par[2])),
#'   nu = c(nu, exp(theta$par[3])), sigma.e = c(sigma.e, exp(theta$par[4])),
#'   row.names = c("Truth", "Estimates")
#' ))
#' }
#'
rSPDE.matern.loglike <- function(object, Y, A, sigma.e, mu = 0,
                                 nu = NULL,
                                 kappa = NULL,
                                 sigma = NULL,
                                 range = NULL,
                                 tau = NULL,
                                 m = NULL) {
  if (inherits(object, "CBrSPDEobj")) {
    return(CBrSPDE.matern.loglike(
      object = object,
      Y = Y, A = A,
      sigma.e = sigma.e,
      mu = mu,
      nu = nu,
      kappa = kappa,
      sigma = sigma,
      tau = tau,
      range = range,
      m = m
    ))
  } else {
    if (inherits(object, "rSPDEobj")) {
      if (object$type == "Matern approximation") {
        object <- update.rSPDEobj(object,
          nu = nu,
          kappa = kappa,
          sigma = sigma,
          tau = tau,
          range = range,
          m = m
        )
        return(rSPDE.loglike(
          obj = object, Y = Y, A = A,
          sigma.e = sigma.e, mu = mu
        ))
      } else {
        stop("The fractional operator should be of type
        'Matern approximation'!")
      }
    } else {
      stop("The object should be of class 'CBrSPDEobj' or 'rSPDEobj'!")
    }
  }
}

#' @name CBrSPDE.matern.loglike
#' @title Object-based log-likelihood function for latent Gaussian fractional
#' SPDE model using the covariance-based rational approximations
#' @description This function evaluates the log-likelihood function for a
#' Gaussian process with a Matern covariance function, that is observed under
#' Gaussian measurement noise:
#' \eqn{Y_i = u(s_i) + \epsilon_i}{Y_i = u(s_i) + \epsilon_i}, where
#' \eqn{\epsilon_i}{\epsilon_i} are iid mean-zero Gaussian variables.
#' The latent model is approximated using
#' the covariance-based rational approximation
#' of the fractional SPDE model corresponding to the Gaussian process.
#' @param object The covariance-based rational SPDE approximation,
#' computed using [matern.operators()]
#' @param Y The observations, either a vector or a matrix where
#' the columns correspond to independent replicates of observations.
#' @param A An observation matrix that links the measurement location
#' to the finite element basis.
#' @param sigma.e The standard deviation of the measurement noise.
#' @param mu Expectation vector of the latent field (default = 0).
#' @param kappa If non-null, update the range parameter of the
#' covariance function.
#' @param tau If non-null, update the parameter tau.
#' @param sigma If non-null, update the standard deviation of
#' the covariance function.
#' @param range If non-null, update the range parameter
#' of the covariance function.
#' @param nu If non-null, update the shape parameter of the
#' covariance function.
#' @param m If non-null, update the order of the rational approximation,
#' which needs to be a positive integer.
#' @return The log-likelihood value.
#' @noRd
#' @seealso [matern.operators()], [predict.CBrSPDEobj()]
#' @examples
#' # this example illustrates how the function can be used for maximum likelihood estimation
#' set.seed(123)
#' # Sample a Gaussian Matern process on R using a rational approximation
#' nu <- 0.8
#' kappa <- 5
#' sigma <- 1
#' sigma.e <- 0.1
#' n.rep <- 10
#' n.obs <- 100
#' n.x <- 51
#' range <- 0.2
#'
#' # create mass and stiffness matrices for a FEM discretization
#' x <- seq(from = 0, to = 1, length.out = n.x)
#' fem <- rSPDE.fem1d(x)
#'
#' tau <- sqrt(gamma(nu) / (sigma^2 * kappa^(2 * nu) * (4 * pi)^(1 / 2) *
#'   gamma(nu + 1 / 2)))
#'
#' # Compute the covariance-based rational approximation
#' op_cov <- matern.operators(
#'   loc_mesh = x, nu = nu,
#'   range = range, sigma = sigma, d = 1, m = 2,
#'   parameterization = "matern"
#' )
#'
#' # Sample the model
#' u <- simulate(op_cov, n.rep)
#'
#' # Create some data
#' obs.loc <- runif(n = n.obs, min = 0, max = 1)
#' A <- rSPDE.A1d(x, obs.loc)
#' noise <- rnorm(n.obs * n.rep)
#' dim(noise) <- c(n.obs, n.rep)
#' Y <- as.matrix(A %*% u + sigma.e * noise)
#'
#' # Define the negative likelihood function for optimization using
#' # CBrSPDE.matern.loglike
#' # Notice that we are also using sigma instead of tau, so it can be compared
#' # to matern.loglike()
#' mlik_cov <- function(theta, Y, A, op_cov) {
#'   kappa <- exp(theta[1])
#'   sigma <- exp(theta[2])
#'   nu <- exp(theta[3])
#'   return(-rSPDE.matern.loglike(
#'     object = op_cov, Y = Y,
#'     A = A, kappa = kappa, sigma = sigma,
#'     nu = nu, sigma.e = exp(theta[4])
#'   ))
#' }
#'
#' # The parameters can now be estimated by minimizing mlik with optim
#' \donttest{
#' # Choose some reasonable starting values depending on the size of the domain
#' theta0 <- log(c(sqrt(8), 1 / sqrt(var(c(Y))), 0.9, 0.01))
#'
#' # run estimation and display the results
#' theta <- optim(theta0, mlik_cov,
#'   Y = Y, A = A, op_cov = op_cov,
#'   method = "L-BFGS-B"
#' )
#'
#' print(data.frame(
#'   kappa = c(kappa, exp(theta$par[1])), sigma = c(sigma, exp(theta$par[2])),
#'   nu = c(nu, exp(theta$par[3])), sigma.e = c(sigma.e, exp(theta$par[4])),
#'   row.names = c("Truth", "Estimates")
#' ))
#' }
CBrSPDE.matern.loglike <- function(object, Y, A, sigma.e, mu = 0,
                                   nu = NULL,
                                   kappa = NULL,
                                   sigma = NULL,
                                   range = NULL,
                                   tau = NULL,
                                   m = NULL) {
  Y <- as.matrix(Y)
  if (length(dim(Y)) == 2) {
    n.rep <- dim(Y)[2]
    n <- dim(Y)[1]
  } else {
    n.rep <- 1
    if (length(dim(Y)) == 1) {
      n <- dim(Y)[1]
    } else {
      n <- length(Y)
    }
  }

  ## get relevant parameters

  object <- update.CBrSPDEobj(
    object = object,
    nu = nu,
    kappa = kappa,
    sigma = sigma,
    range = range,
    tau = tau,
    m = m,
    parameterization = object$parameterization
  )

  m <- object$m

  if (length(sigma.e) == 1) {
    Q.e <- Diagonal(n) / sigma.e^2
    nugget <- rep(sigma.e^2, n)
  } else {
    if (length(sigma.e) != n) {
      stop("the length of sigma.e does not match the number of observations")
    }
    Q.e <- Diagonal(length(sigma.e), 1 / sigma.e^2)
    nugget <- sigma.e^2
  }

  if (object$compute_logdet) {
    Q.frac <- object$Q.frac

    Q.fracR <- Matrix::Cholesky(Q.frac)

    logdetL <- object$logdetL
    logdetC <- object$logdetC
    Q.int.order <- object$Q.int$order

    if (Q.int.order > 0) {
      # logQ <- 2 * sum(log(diag(Q.fracR))) + (Q.int.order) *
      # (m + 1) * (logdetL - logdetC)

      logQ <- 2 * c(determinant(Q.fracR, logarithm = TRUE, sqrt = TRUE)$modulus) + (Q.int.order) *
        (m + 1) * (logdetL - logdetC)
    } else {
      logQ <- 2 * c(determinant(Q.fracR, logarithm = TRUE, sqrt = TRUE)$modulus)
    }
  } else {
    Q <- object$Q
    Q.R <- Matrix::Cholesky(Q)

    logQ <- 2 * c(determinant(Q.R, logarithm = TRUE, sqrt = TRUE)$modulus)
  }
  ## compute Q_x|y
  Q <- object$Q
  if (object$alpha %% 1 == 0) {
    Abar <- A
  } else {
    Abar <- kronecker(matrix(1, 1, m + 1), A)
  }
  Q_xgiveny <- t(Abar) %*% Q.e %*% Abar + Q
  ## construct mu_x|y

  mu_xgiveny <- t(Abar) %*% Q.e %*% Y
  # upper triangle with reordering


  R <- Matrix::Cholesky(Q_xgiveny)

  mu_xgiveny <- solve(R, mu_xgiveny, system = "A")

  mu_xgiveny <- mu + mu_xgiveny

  ## compute log|Q_xgiveny|
  log_Q_xgiveny <- 2 * determinant(R, logarithm = TRUE, sqrt = TRUE)$modulus
  ## compute mu_x|y*Q*mu_x|y
  if (n.rep > 1) {
    mu_part <- sum(colSums((mu_xgiveny - mu) * (Q %*% (mu_xgiveny - mu))))
  } else {
    mu_part <- t(mu_xgiveny - mu) %*% Q %*% (mu_xgiveny - mu)
  }
  ## compute central part
  if (n.rep > 1) {
    central_part <- sum(colSums((Y - Abar %*% mu_xgiveny) * (Q.e %*%
      (Y - Abar %*% mu_xgiveny))))
  } else {
    central_part <- t(Y - Abar %*% mu_xgiveny) %*% Q.e %*% (Y -
      Abar %*% mu_xgiveny)
  }
  ## compute log|Q_epsilon|
  log_Q_epsilon <- -sum(log(nugget))
  ## wrap up
  log_likelihood <- n.rep * (logQ + log_Q_epsilon - log_Q_xgiveny) -
    mu_part - central_part
  if (n.rep > 1) {
    log_likelihood <- log_likelihood - dim(A)[1] * n.rep * log(2 * pi)
  } else {
    log_likelihood <- log_likelihood - length(Y) * log(2 * pi)
  }
  log_likelihood <- log_likelihood / 2

  return(as.double(log_likelihood))
}

#' @noRd

aux_CBrSPDE.matern.loglike <- function(object, Y, A, sigma.e, mu = 0,
                                       nu = NULL,
                                       kappa = NULL,
                                       tau = NULL,
                                       theta = NULL,
                                       m = NULL) {
  Y <- as.matrix(Y)
  if (length(dim(Y)) == 2) {
    n.rep <- dim(Y)[2]
    n <- dim(Y)[1]
  } else {
    n.rep <- 1
    if (length(dim(Y)) == 1) {
      n <- dim(Y)[1]
    } else {
      n <- length(Y)
    }
  }


  ## get relevant parameters

  object <- update.CBrSPDEobj(
    object = object,
    nu = nu,
    kappa = kappa,
    theta = theta,
    tau = tau,
    m = m
  )

  m <- object$m

  if (length(sigma.e) == 1) {
    Q.e <- Diagonal(n) / sigma.e^2
    nugget <- rep(sigma.e^2, n)
  } else {
    if (length(sigma.e) != n) {
      stop("the length of sigma.e does not match the number of observations")
    }
    Q.e <- Diagonal(length(sigma.e), 1 / sigma.e^2)
    nugget <- sigma.e^2
  }

  Q <- object$Q

  if (object$stationary && object$compute_logdet) {
    Q.frac <- object$Q.frac

    Q.fracR <- Matrix::Cholesky(Q.frac)

    logdetL <- object$logdetL
    logdetC <- object$logdetC
    Q.int.order <- object$Q.int$order

    if (Q.int.order > 0) {
      # logQ <- 2 * sum(log(diag(Q.fracR))) + (Q.int.order) *
      # (m + 1) * (logdetL - logdetC)

      logQ <- 2 * c(determinant(Q.fracR, logarithm = TRUE, sqrt = TRUE)$modulus) + (Q.int.order) *
        (m + 1) * (logdetL - logdetC)
    } else {
      logQ <- 2 * c(determinant(Q.fracR, logarithm = TRUE, sqrt = TRUE)$modulus)
    }
  } else {
    Q.R <- Matrix::Cholesky(Q)

    logQ <- 2 * c(determinant(Q.R, logarithm = TRUE, sqrt = TRUE)$modulus)
  }


  ## compute Q_x|y
  if (object$alpha %% 1 == 0) {
    Abar <- A
  } else {
    Abar <- kronecker(matrix(1, 1, m + 1), A)
  }

  Q_xgiveny <- t(Abar) %*% Q.e %*% Abar + Q
  ## construct mu_x|y

  mu_xgiveny <- t(Abar) %*% Q.e %*% Y
  # upper triangle with reordering


  R <- Matrix::Cholesky(Q_xgiveny)

  mu_xgiveny <- solve(R, mu_xgiveny, system = "A")

  mu_xgiveny <- mu + mu_xgiveny

  ## compute log|Q_xgiveny|
  log_Q_xgiveny <- 2 * determinant(R, logarithm = TRUE, sqrt = TRUE)$modulus
  ## compute mu_x|y*Q*mu_x|y
  if (n.rep > 1) {
    mu_part <- sum(colSums((mu_xgiveny - mu) * (Q %*% (mu_xgiveny - mu))))
  } else {
    mu_part <- t(mu_xgiveny - mu) %*% Q %*% (mu_xgiveny - mu)
  }
  ## compute central part
  if (n.rep > 1) {
    central_part <- sum(colSums((Y - Abar %*% mu_xgiveny) * (Q.e %*%
      (Y - Abar %*% mu_xgiveny))))
  } else {
    central_part <- t(Y - Abar %*% mu_xgiveny) %*% Q.e %*% (Y -
      Abar %*% mu_xgiveny)
  }
  ## compute log|Q_epsilon|
  log_Q_epsilon <- -sum(log(nugget))
  ## wrap up
  log_likelihood <- n.rep * (logQ + log_Q_epsilon - log_Q_xgiveny) -
    mu_part - central_part
  if (n.rep > 1) {
    log_likelihood <- log_likelihood - dim(A)[1] * n.rep * log(2 * pi)
  } else {
    log_likelihood <- log_likelihood - length(Y) * log(2 * pi)
  }
  log_likelihood <- log_likelihood / 2

  return(as.double(log_likelihood))
}


#' Parameter-based log-likelihood for a latent Gaussian Matern SPDE model
#' using a rational SPDE approximation
#'
#' This function evaluates the log-likelihood function for observations
#' of a Gaussian process defined as the solution to the SPDE
#' \deqn{(\kappa(s) - \Delta)^\beta (\tau(s)u(s)) = W.}
#'
#' The observations are assumed to be generated as
#' \eqn{Y_i = u(s_i) + \epsilon_i}{Y_i = u(s_i) + \epsilon_i}, where
#' \eqn{\epsilon_i}{\epsilon_i} are
#' iid mean-zero Gaussian variables. The latent model is approximated using a
#' rational approximation of the fractional SPDE model.
#'
#' @param object The rational SPDE approximation,
#' computed using [spde.matern.operators()]
#' @param Y The observations, either a vector or a matrix where
#' the columns correspond to independent replicates of observations.
#' @param A An observation matrix that links the measurement location to the
#' finite element basis.
#' @param sigma.e IF non-null, the standard deviation of the measurement noise will be kept fixed in
#' the returned likelihood.
#' @param mu Expectation vector of the latent field (default = 0).
#' @param kappa If non-null, updates the range parameter.
#' @param tau If non-null, updates the parameter tau.
#' @param theta If non-null, updates the parameter theta (that connects tau and kappa to the model matrices in `object`).
#' @param nu If non-null, the shape parameter will be kept fixed in the returned likelihood.
#' @param m If non-null, update the order of the rational approximation,
#' which needs to be a positive integer.
#'
#' @return The log-likelihood value.
#' @export
#' @seealso [rSPDE.loglike()].
#'
#' @examples
#' # this example illustrates how the function can be used for maximum
#' # likelihood estimation
#' # Sample a Gaussian Matern process on R using a rational approximation
#' sigma.e <- 0.1
#' n.rep <- 10
#' n.obs <- 100
#' n.x <- 51
#' # create mass and stiffness matrices for a FEM discretization
#' x <- seq(from = 0, to = 1, length.out = n.x)
#' fem <- rSPDE.fem1d(x)
#' tau <- rep(0.5, n.x)
#' nu <- 0.8
#' alpha <- nu + 1 / 2
#' kappa <- rep(1, n.x)
#' # compute rational approximation
#' op <- spde.matern.operators(
#'   kappa = kappa, tau = tau, alpha = alpha,
#'   parameterization = "spde", d = 1,
#'   loc_mesh = x
#' )
#' # Sample the model
#' u <- simulate(op, n.rep)
#' # Create some data
#' obs.loc <- runif(n = n.obs, min = 0, max = 1)
#' A <- rSPDE.A1d(x, obs.loc)
#' noise <- rnorm(n.obs * n.rep)
#' dim(noise) <- c(n.obs, n.rep)
#' Y <- as.matrix(A %*% u + sigma.e * noise)
#' # define negative likelihood function for optimization using matern.loglike
#' mlik <- function(theta) {
#'   return(-spde.matern.loglike(op, Y, A,
#'     sigma.e = exp(theta[4]),
#'     nu = exp(theta[3]),
#'     kappa = exp(theta[2]),
#'     tau = exp(theta[1])
#'   ))
#' }
#' #' #The parameters can now be estimated by minimizing mlik with optim
#' \donttest{
#' # Choose some reasonable starting values depending on the size of the domain
#' theta0 <- log(c(1 / sqrt(var(c(Y))), sqrt(8), 0.9, 0.01))
#' # run estimation and display the results
#' theta <- optim(theta0, mlik)
#' print(data.frame(
#'   tau = c(tau[1], exp(theta$par[1])), kappa = c(kappa[1], exp(theta$par[2])),
#'   nu = c(nu, exp(theta$par[3])), sigma.e = c(sigma.e, exp(theta$par[4])),
#'   row.names = c("Truth", "Estimates")
#' ))
#' }
spde.matern.loglike <- function(object, Y, A, sigma.e, mu = 0,
                                nu = NULL,
                                kappa = NULL,
                                tau = NULL,
                                theta = NULL,
                                m = NULL) {
  if (inherits(object, "CBrSPDEobj")) {
    object <- update.CBrSPDEobj(object,
      nu = nu,
      kappa = kappa,
      tau = tau,
      theta = theta,
      m = m
    )

    return(aux_CBrSPDE.matern.loglike(
      object = object, Y = Y, A = A, sigma.e = sigma.e, mu = mu,
      nu = nu,
      kappa = kappa,
      tau = tau,
      theta = theta,
      m = m
    ))
  } else {
    if (inherits(object, "rSPDEobj")) {
      if (object$type == "Matern approximation") {
        object <- update.rSPDEobj(object,
          nu = nu,
          kappa = kappa,
          tau = tau,
          theta = theta,
          m = m
        )
        return(rSPDE.loglike(
          obj = object, Y = Y, A = A,
          sigma.e = sigma.e, mu = mu
        ))
      } else {
        stop("The fractional operator should be of type
        'Matern approximation'!")
      }
    } else {
      stop("The object should be of class 'CBrSPDEobj' or 'rSPDEobj'!")
    }
  }
}



#' @name predict.CBrSPDEobj
#' @title Prediction of a fractional SPDE using the covariance-based
#' rational SPDE approximation
#' @description The function is used for computing kriging predictions based
#' on data \eqn{Y_i = u(s_i) + \epsilon_i}, where \eqn{\epsilon}{\epsilon}
#' is mean-zero Gaussian measurement noise and \eqn{u(s)}{u(s)} is defined by
#' a fractional SPDE \eqn{(\kappa^2 I - \Delta)^{\alpha/2} (\tau u(s)) = W},
#' where \eqn{W}{W} is Gaussian white noise and \eqn{\alpha = \nu + d/2},
#' where \eqn{d} is the dimension of the domain.
#' @param object The covariance-based rational SPDE approximation,
#' computed using [matern.operators()]
#' @param A A matrix linking the measurement locations to the basis of the FEM
#' approximation of the latent model.
#' @param Aprd A matrix linking the prediction locations to the basis of the
#' FEM approximation of the latent model.
#' @param Y A vector with the observed data, can also be a matrix where the
#' columns are observations
#' of independent replicates of \eqn{u}.
#' @param sigma.e The standard deviation of the Gaussian measurement noise.
#' Put to zero if the model does not have measurement noise.
#' @param mu Expectation vector of the latent field (default = 0).
#' @param compute.variances Set to also TRUE to compute the kriging variances.
#' @param posterior_samples If `TRUE`, posterior samples will be returned.
#' @param n_samples Number of samples to be returned. Will only be used if `sampling` is `TRUE`.
#' @param compute_var_method Method to compute the variances. Options are:
#' "direct", "loop", and "selected_inv". The "direct" method computes the variance directly,
#' which is faster but can be memory intensive. The "loop" method computes the variance by looping
#' over the elements, which is more memory efficient but slower.
#' @param only_latent Should the posterior samples be only given to the laten model?
#' @param ... further arguments passed to or from other methods.
#' @return A list with elements
#' \item{mean }{The kriging predictor (the posterior mean of u|Y).}
#' \item{variance }{The posterior variances (if computed).}
#' @export
#' @method predict CBrSPDEobj
#' @examples
#' set.seed(123)
#' # Sample a Gaussian Matern process on R using a rational approximation
#' kappa <- 10
#' sigma <- 1
#' nu <- 0.8
#' sigma.e <- 0.3
#' range <- 0.2
#'
#' # create mass and stiffness matrices for a FEM discretization
#' x <- seq(from = 0, to = 1, length.out = 101)
#' fem <- rSPDE.fem1d(x)
#'
#' tau <- sqrt(gamma(nu) / (sigma^2 * kappa^(2 * nu) *
#'   (4 * pi)^(1 / 2) * gamma(nu + 1 / 2)))
#'
#' # Compute the covariance-based rational approximation
#' op_cov <- matern.operators(
#'   loc_mesh = x, nu = nu,
#'   range = range, sigma = sigma, d = 1, m = 2,
#'   parameterization = "matern"
#' )
#'
#' # Sample the model
#' u <- simulate(op_cov)
#'
#' # Create some data
#' obs.loc <- runif(n = 10, min = 0, max = 1)
#' A <- rSPDE.A1d(x, obs.loc)
#' Y <- as.vector(A %*% u + sigma.e * rnorm(10))
#'
#' # compute kriging predictions at the FEM grid
#' A.krig <- rSPDE.A1d(x, x)
#' u.krig <- predict(op_cov,
#'   A = A, Aprd = A.krig, Y = Y, sigma.e = sigma.e,
#'   compute.variances = TRUE
#' )
#'
#' plot(obs.loc, Y,
#'   ylab = "u(x)", xlab = "x", main = "Data and prediction",
#'   ylim = c(
#'     min(u.krig$mean - 2 * sqrt(u.krig$variance)),
#'     max(u.krig$mean + 2 * sqrt(u.krig$variance))
#'   )
#' )
#' lines(x, u.krig$mean)
#' lines(x, u.krig$mean + 2 * sqrt(u.krig$variance), col = 2)
#' lines(x, u.krig$mean - 2 * sqrt(u.krig$variance), col = 2)
predict.CBrSPDEobj <- function(object, A, Aprd, Y, sigma.e, mu = 0,
                               compute.variances = FALSE, posterior_samples = FALSE,
                               n_samples = 100, only_latent = FALSE, compute_var_method = c("direct", "loop"),
                               ...) {
  compute_var_method <- match.arg(compute_var_method)
  Y <- as.matrix(Y)
  if (dim(Y)[1] != dim(A)[1]) {
    stop("the dimensions of A does not match the number of observations")
  }

  n <- dim(Y)[1]
  out <- list()

  d <- object$d
  kappa <- object$kappa
  nu <- object$nu
  tau <- object$tau
  m <- object$m

  fem_mesh_matrices <- object$fem_mesh_matrices

  no_nugget <- FALSE

  if (length(sigma.e) == 1) {
    if (sigma.e == 0) {
      no_nugget <- TRUE
    } else {
      Q.e <- Diagonal(n) / sigma.e^2
    }
  } else {
    if (length(sigma.e) != n) {
      stop("the length of sigma.e does not match the number of observations")
    }
    Q.e <- Diagonal(length(sigma.e), 1 / sigma.e^2)
  }

  alpha <- nu + d / 2
  
  if(length(mu) == 1) {
      mu <- rep(mu, dim(object$Q)[1])
  } else {
      if(length(mu) == dim(object$C)[1] && dim(object$C)[1] < dim(object$Q)[1]) {
          mu <- rep(mu,object$m+1) / (object$m+1)
      } else if (length(mu) != dim(object$Q)[1]) {
          stop("the length of mu is wrong.")
      }
  }

  if (!no_nugget) {
    if (alpha %% 1 == 0) { # loglikelihood in integer case
      ## construct Q
      Q <- object$Q
      ## compute Q_x|y
      Q_xgiveny <- (t(A) %*% Q.e %*% A) + Q
      ## construct mu_x|y
      mu_xgiveny <- t(A) %*% Q.e %*% (Y - A%*%mu)

      R <- Matrix::Cholesky(forceSymmetric(Q_xgiveny))
      mu_xgiveny <- solve(R, mu_xgiveny, system = "A")

      mu_xgiveny <- mu + mu_xgiveny
      out$mean <- Aprd %*% mu_xgiveny

      if (compute.variances) {
        out$variance <- diag(Aprd %*% solve(Q_xgiveny, t(Aprd)))
      }
    } else { # loglikelihood in non-integer case

      Q <- object$Q

      ## compute Q_x|y
      Q_xgiveny <- kronecker(matrix(1, m + 1, m + 1), t(A) %*% Q.e %*% A) + Q
      ## construct mu_x|y
      Abar <- kronecker(matrix(1, 1, m + 1), A)
      
      mu_xgiveny <- t(Abar) %*% Q.e %*% (Y - Abar%*%mu)
      # upper triangle with reordering
      R <- Matrix::Cholesky(forceSymmetric(Q_xgiveny))

      mu_xgiveny <- solve(R, mu_xgiveny, system = "A")

      mu_xgiveny <- mu + mu_xgiveny

      Aprd_bar <- kronecker(matrix(1, 1, m + 1), Aprd)

      out$mean <- Aprd_bar %*% mu_xgiveny

      if (compute.variances) {
        if(compute_var_method == "direct"){
          out$variance <- diag(Aprd_bar %*% solve(Q_xgiveny, t(Aprd_bar)))
        } else if (compute_var_method == "loop"){
          out$variance <- rep(0, nrow(Aprd_bar))
          for (i in 1:nrow(Aprd_bar)) {
            out$variance[i] <- Aprd_bar[i, ] %*% solve(Q_xgiveny, Aprd_bar[i, ])
          }
        } 
        # else{
        #   sel_inv_Qxgiveny <- MetricGraph::selected_inv(Q_xgiveny)
        #   out$variance <- diag(Aprd_bar %*% sel_inv_Qxgiveny %*% t(Aprd_bar))
        # }
      }
    }
  } else {
    integer_alpha <- (alpha %% 1 == 0)
    if (integer_alpha) {
      Abar <- A
      Aprd_bar <- Aprd
      Q <- object$Q
    } else {
      Abar <- kronecker(matrix(1, 1, m + 1), A)
      Aprd_bar <- kronecker(matrix(1, 1, m + 1), Aprd)
      Q <- object$Q
    }

    QiAt <- solve(Q, t(Abar))
    AQiA <- Abar %*% QiAt
    xhat <- solve(Q, t(Abar) %*% solve(AQiA, Y))

    out$mean <- as.vector(Aprd_bar %*% xhat)
    if (compute.variances) {
      M <- Q - QiAt %*% solve(AQiA, t(QiAt))
      out$variance <- diag(Aprd_bar %*% M %*% t(Aprd_bar))
    }
  }


  if (posterior_samples) {
      if (!no_nugget) {
          Z <- rnorm(dim(object$Q)[1] * n_samples)
          dim(Z) <- c(dim(object$Q)[1], n_samples)
          LQ <-  chol(forceSymmetric(Q_xgiveny))
          X <- as.matrix(solve(LQ, Z)) + kronecker(as.matrix(mu_xgiveny), 
                                                   matrix(rep(1,n_samples),1,n_samples))
          X <- Aprd %*% X
          if (!only_latent) {
              X <- X + matrix(rnorm(n_samples * dim(Aprd)[1], sd = sigma.e), nrow = dim(Aprd)[1])
          }
          return(X)
      } else {
          M <- Q - QiAt %*% solve(AQiA, t(QiAt))
          post_cov <- Aprd %*% M %*% t(Aprd)
          Y_tmp <- as.matrix(Y)
          mean_tmp <- as.matrix(out$mean)
          out$samples <- lapply(1:ncol(Y_tmp), function(i) {
              Z <- rnorm(dim(post_cov)[1] * n_samples)
              dim(Z) <- c(dim(post_cov)[1], n_samples)
              LQ <-  Matrix::Cholesky(forceSymmetric(post_cov))
              X <- LQ %*% Z
              X <- X + mean_tmp[, i]
              if (!only_latent) {
                  X <- X + matrix(rnorm(n_samples * dim(Aprd)[1], sd = sigma.e), nrow = dim(Aprd)[1])
              }
              return(X)
          })
      }
  }


  return(out)
}


#' @name predict.CBrSPDEobj2d
#' @title Prediction of an anisotropic Whittle-Matern field
#' @description The function is used for computing kriging predictions based
#' on data \eqn{Y_i = u(s_i) + \epsilon_i}, where \eqn{\epsilon}{\epsilon}
#' is mean-zero Gaussian measurement noise and \eqn{u(s)}{u(s)} is defined by
#' a SPDE as described in [matern2d.operators()].
#' @param object The covariance-based rational SPDE approximation,
#' computed using [matern2d.operators()]
#' @param A A matrix linking the measurement locations to the basis of the FEM
#' approximation of the latent model.
#' @param Aprd A matrix linking the prediction locations to the basis of the
#' FEM approximation of the latent model.
#' @param Y A vector with the observed data, can also be a matrix where the
#' columns are observations
#' of independent replicates of \eqn{u}.
#' @param sigma.e The standard deviation of the Gaussian measurement noise.
#' Put to zero if the model does not have measurement noise.
#' @param mu Expectation vector of the latent field (default = 0).
#' @param compute.variances Set to also TRUE to compute the kriging variances.
#' @param posterior_samples If `TRUE`, posterior samples will be returned.
#' @param n_samples Number of samples to be returned. Will only be used if `sampling` is `TRUE`.
#' @param only_latent Should the posterior samples be only given to the laten model?
#' @param ... further arguments passed to or from other methods.
#' @return A list with elements
#' \item{mean }{The kriging predictor (the posterior mean of u|Y).}
#' \item{variance }{The posterior variances (if computed).}
#' @export
#' @method predict CBrSPDEobj2d
#' @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.01, max.edge = c(0.1, 0.5))
#'  op <- matern2d.operators(hx = 0.08, hy = 0.08, hxy = 0.5, nu = 0.5, 
#'  sigma = 1, mesh = mesh_2d)
#'  u <- simulate(op)
#'  n.obs <- 2000
#'  obs.loc <- cbind(runif(n.obs),runif(n.obs))
#'  A <- fm_basis(mesh_2d,obs.loc)
#'  sigma.e <- 0.1
#'  Y <- as.vector(A%*%u + sigma.e*rnorm(n.obs))
#'  A <- make_A(op, obs.loc)
#'  proj <- fm_evaluator(mesh_2d, dims = c(100, 100),
#'              xlim = c(0,1), ylim = c(0,1))
#'  Aprd <- make_A(op, proj$lattice$loc)
#'  u.krig <- predict(op, A = A, Aprd = Aprd, Y = Y, sigma.e = sigma.e)
predict.CBrSPDEobj2d <- function(object, A, Aprd, Y, sigma.e, mu = 0,
                                 compute.variances = FALSE, posterior_samples = FALSE,
                                 n_samples = 100, only_latent = FALSE,
                                 ...) {
    Y <- as.matrix(Y)
    if (dim(Y)[1] != dim(A)[1]) {
        stop("the dimensions of A does not match the number of observations")
    }
    
    n <- dim(Y)[1]
    out <- list()
    
    no_nugget <- FALSE
    
    if (length(sigma.e) == 1) {
        if (sigma.e == 0) {
            no_nugget <- TRUE
        } else {
            Q.e <- Diagonal(n) / sigma.e^2
        }
    } else {
        if (length(sigma.e) != n) {
            stop("the length of sigma.e does not match the number of observations")
        }
        Q.e <- Diagonal(length(sigma.e), 1 / sigma.e^2)
    }
    
    
    if (!no_nugget) {
        ## construct Q
        Q <- object$Q
        ## compute Q_x|y
        Q_xgiveny <- (t(A) %*% Q.e %*% A) + Q
        ## construct mu_x|y
        mu_xgiveny <- t(A) %*% Q.e %*% Y
        
        R <- Matrix::Cholesky(forceSymmetric(Q_xgiveny))
        mu_xgiveny <- solve(R, mu_xgiveny, system = "A")
        
        mu_xgiveny <- mu + mu_xgiveny
        out$mean <- Aprd %*% mu_xgiveny
        
        if (compute.variances) {
            out$variance <- diag(Aprd %*% solve(Q_xgiveny, t(Aprd)))
        }
    } else {
        Q <- object$Q
        
        QiAt <- solve(Q, t(A))
        AQiA <- A %*% QiAt
        xhat <- solve(Q, t(A) %*% solve(AQiA, Y))
        
        out$mean <- as.vector(Aprd %*% xhat)
        if (compute.variances) {
            M <- Q - QiAt %*% solve(AQiA, t(QiAt))
            out$variance <- diag(Aprd %*% M %*% t(Aprd))
        }
    }
    
    
    if (posterior_samples) {
        if (!no_nugget) {
            post_cov <- Aprd %*% solve(Q_xgiveny, t(Aprd))
        } else {
            M <- Q - QiAt %*% solve(AQiA, t(QiAt))
            post_cov <- Aprd %*% M %*% t(Aprd)
        }
        Y_tmp <- as.matrix(Y)
        mean_tmp <- as.matrix(out$mean)
        out$samples <- lapply(1:ncol(Y_tmp), function(i) {
            Z <- rnorm(dim(post_cov)[1] * n_samples)
            dim(Z) <- c(dim(post_cov)[1], n_samples)
            LQ <- chol(forceSymmetric(post_cov))
            X <- LQ %*% Z
            X <- X + mean_tmp[, i]
            if (!only_latent) {
                X <- X + matrix(rnorm(n_samples * dim(Aprd)[1], sd = sigma.e), nrow = dim(Aprd)[1])
            }
            return(X)
        })
    }
    return(out)
}

#' @rdname precision.CBrSPDEobj
#' @export
precision <- function(object, ...) {
  UseMethod("precision", object)
}

#' @name precision.CBrSPDEobj
#' @title Get the precision matrix of CBrSPDEobj objects
#' @description Function to get the precision matrix of a CBrSPDEobj object
#' @param object The covariance-based rational SPDE approximation,
#' computed using [matern.operators()]
#' @param kappa If non-null, update the range parameter of
#' the covariance function.
#' @param tau If non-null, update the parameter tau.
#' @param sigma If non-null, update the standard deviation of
#' the covariance function.
#' @param range If non-null, update the range parameter
#' of the covariance function.
#' @param nu If non-null, update the shape parameter of the
#' covariance function.
#' @param m If non-null, update the order of the rational approximation,
#' which needs to be a positive integer.
#' @param ... Currently not used.
#' @return The precision matrix.
#' @method precision CBrSPDEobj
#' @seealso [simulate.CBrSPDEobj()], [matern.operators()]
#' @export
#' @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 <- 0.2
#'
#' # 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_cov <- matern.operators(
#'   loc_mesh = x, nu = nu,
#'   range = range, sigma = sigma, d = 1, m = 2,
#'   parameterization = "matern"
#' )
#'
#' # Get the precision matrix:
#' prec_matrix <- precision(op_cov)
#'
precision.CBrSPDEobj <- function(object,
                                 nu = NULL,
                                 kappa = NULL,
                                 sigma = NULL,
                                 range = NULL,
                                 tau = NULL,
                                 m = NULL,
                                 ...) {
  object <- update.CBrSPDEobj(
    object = object,
    nu = nu,
    kappa = kappa,
    sigma = sigma,
    range = range,
    tau = tau,
    m = m
  )

  Q <- object$Q
  return(Q)
}

#' @name precision.CBrSPDEobj2d
#' @title Get the precision matrix of CBrSPDEobj2d objects
#' @description Function to get the precision matrix of a CBrSPDEobj2d object
#' @param object The covariance-based rational SPDE approximation,
#' computed using [matern2d.operators()]
#' @param nu If non-null, update the shape parameter of the
#' covariance function.
#' @param hx If non-null, update the hx parameter. 
#' @param hy If non-null, update the hy parameter. 
#' @param hxy If non-null, update the hxy parameter. 
#' @param sigma If non-null, update the standard deviation of
#' the covariance function.
#' @param m If non-null, update the order of the rational approximation,
#' which needs to be a positive integer.
#' @param ... Currently not used.
#' @return The precision matrix.
#' @method precision CBrSPDEobj2d
#' @seealso [simulate.CBrSPDEobj2d()], [matern2d.operators()]
#' @export
#' @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)
#' Q <- precision(op)
precision.CBrSPDEobj2d <- function(object,
                                 nu = NULL,
                                 hx = NULL,
                                 hy = NULL,
                                 hxy = NULL,
                                 sigma = NULL,
                                 m = NULL,
                                 ...) {
    object <- update.CBrSPDEobj2d(
        object = object,
        nu = nu,
        hx = hx,
        hy = hy,
        hxy = hxy,
        sigma = sigma,
        m = m
    )
    
    Q <- object$Q
    return(Q)
}


#' @name rSPDE.construct.matern.loglike
#' @title Constructor of Matern loglikelihood functions.
#' @description This function returns a log-likelihood function for a
#' Gaussian process with a Matern covariance
#' function, that is observed under Gaussian measurement noise:
#' \eqn{Y_i = u(s_i) + \epsilon_i}{Y_i = u(s_i) + \epsilon_i}, where
#' \eqn{\epsilon_i}{\epsilon_i} are
#' iid mean-zero Gaussian variables. The latent model is approximated using
#' the a rational approximation
#' of the fractional SPDE model corresponding to the Gaussian process.
#' @param object The rational SPDE approximation,
#' computed using [matern.operators()]
#' @param Y The observations, either a vector or a matrix where
#' the columns correspond to independent replicates of observations.
#' @param A An observation matrix that links the measurement location to the
#' finite element basis.
#' @param sigma.e IF non-null, the standard deviation of the measurement noise will be kept fixed in
#' the returned likelihood.
#' @param mu Expectation vector of the latent field (default = 0).
#' @param kappa If non-null, the range parameter will be kept fixed in the returned likelihood.
#' @param range If non-null, the range parameter will be kept fixed in the returned likelihood. (Replaces kappa)
#' @param sigma If non-null, the standard deviation will be kept fixed in the returned likelihood.
#' @param tau If non-null, the tau parameter will be kept fixed in the returned likelihood. (Replaces sigma)
#' @param nu If non-null, the shape parameter will be kept fixed in the returned likelihood.
#' @param parameterization If `spde`, then one will use the parameters `tau` and `kappa`. If `matern`, then one will use the parameters `sigma` and `range`.
#' @param m If non-null, update the order of the rational approximation,
#' which needs to be a positive integer.
#' @param log_scale Should the parameters be evaluated in log-scale?
#' @param return_negative_likelihood Return minus the likelihood to turn the maximization into a minimization?
#' @return The log-likelihood function. The parameters of the returned function
#' are given in the order sigma, kappa, nu, sigma.e, whenever they are available.
#' @export
#' @seealso [matern.operators()], [predict.CBrSPDEobj()]
#' @examples
#' # this example illustrates how the function can be used for maximum
#' # likelihood estimation
#'
#' set.seed(123)
#' # Sample a Gaussian Matern process on R using a rational approximation
#' nu <- 0.8
#' sigma <- 1
#' sigma.e <- 0.1
#' n.rep <- 10
#' n.obs <- 200
#' n.x <- 51
#' range <- 0.2
#' # create mass and stiffness matrices for a FEM discretization
#' x <- seq(from = 0, to = 1, length.out = n.x)
#' # Compute the covariance-based rational approximation
#' op_cov <- matern.operators(
#'   loc_mesh = x, nu = nu,
#'   range = range, sigma = sigma, d = 1, m = 2,
#'   parameterization = "matern"
#' )
#' # Sample the model
#' u <- simulate(op_cov, n.rep)
#' # Create some data
#' obs.loc <- runif(n = n.obs, min = 0, max = 1)
#' A <- rSPDE.A1d(x, obs.loc)
#' noise <- rnorm(n.obs * n.rep)
#' dim(noise) <- c(n.obs, n.rep)
#' Y <- as.matrix(A %*% u + sigma.e * noise)
#' \donttest{
#' # Define the negative likelihood function for optimization
#' # using CBrSPDE.matern.loglike
#' # Matern parameterization
#' loglike <- rSPDE.construct.matern.loglike(op_cov, Y, A, parameterization = "matern")
#'
#' # The parameters can now be estimated by minimizing mlik with optim
#'
#' # Choose some reasonable starting values depending on the size of the domain
#' theta0 <- c(
#'   get.initial.values.rSPDE(mesh.range = 1, dim = 1),
#'   log(0.1 * sd(as.vector(Y)))
#' )
#' # run estimation and display the results
#' theta <- optim(theta0, loglike,
#'   method = "L-BFGS-B"
#' )
#' print(data.frame(
#'   sigma = c(sigma, exp(theta$par[1])), range = c(range, exp(theta$par[2])),
#'   nu = c(nu, exp(theta$par[3])), sigma.e = c(sigma.e, exp(theta$par[4])),
#'   row.names = c("Truth", "Estimates")
#' ))
#' }
rSPDE.construct.matern.loglike <- function(object, Y, A,
                                           sigma.e = NULL, mu = 0,
                                           nu = NULL,
                                           tau = NULL,
                                           kappa = NULL,
                                           sigma = NULL,
                                           range = NULL,
                                           parameterization = c("spde", "matern"),
                                           m = NULL,
                                           log_scale = TRUE,
                                           return_negative_likelihood = TRUE) {
  parameterization <- parameterization[[1]]

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

  if (parameterization == "spde") {
    param_vector <- likelihood_process_inputs_spde(kappa, tau, nu, sigma.e)
  } else {
    param_vector <- likelihood_process_inputs_matern(range, sigma, nu, sigma.e)
  }


  if (parameterization == "spde") {
    loglik <- function(theta) {
      if (is.null(tau)) {
        tau <- likelihood_process_parameters(
          theta = theta,
          param_vector = param_vector,
          which_par = "tau",
          logscale = log_scale
        )
      } else {
        tau <- tau
      }
      if (is.null(kappa)) {
        kappa <- likelihood_process_parameters(
          theta = theta,
          param_vector = param_vector,
          which_par = "kappa",
          logscale = log_scale
        )
      } else {
        kappa <- kappa
      }
      if (is.null(nu)) {
        nu <- likelihood_process_parameters(
          theta = theta,
          param_vector = param_vector,
          which_par = "nu",
          logscale = log_scale
        )
      } else {
        nu <- nu
      }
      if (is.null(sigma.e)) {
        sigma.e <- likelihood_process_parameters(
          theta = theta,
          param_vector = param_vector,
          which_par = "sigma.e",
          logscale = log_scale
        )
      }
      nu <- min(nu, 10)
      if (nu %% 1 == 0) {
        nu <- nu + 1e-10
      }
      loglike <- rSPDE.matern.loglike(
        object = object, Y = Y, A = A,
        sigma.e = sigma.e,
        mu = mu,
        kappa = kappa,
        nu = nu,
        tau = tau,
        m = m
      )
      if (return_negative_likelihood) {
        return(-loglike)
      } else {
        return(loglike)
      }
    }
  } else {
    loglik <- function(theta) {
      if (is.null(sigma)) {
        sigma <- likelihood_process_parameters(
          theta = theta,
          param_vector = param_vector,
          which_par = "sigma",
          logscale = log_scale
        )
      } else {
        sigma <- sigma
      }
      if (is.null(range)) {
        range <- likelihood_process_parameters(
          theta = theta,
          param_vector = param_vector,
          which_par = "range",
          logscale = log_scale
        )
      } else {
        range <- range
      }
      if (is.null(nu)) {
        nu <- likelihood_process_parameters(
          theta = theta,
          param_vector = param_vector,
          which_par = "nu",
          logscale = log_scale
        )
      } else {
        nu <- nu
      }
      if (is.null(sigma.e)) {
        sigma.e <- likelihood_process_parameters(
          theta = theta,
          param_vector = param_vector,
          which_par = "sigma.e",
          logscale = log_scale
        )
      }
      nu <- min(nu, 10)
      if (nu %% 1 == 0) {
        nu <- nu + 1e-10
      }
      loglike <- rSPDE.matern.loglike(
        object = object, Y = Y, A = A,
        sigma.e = sigma.e,
        mu = mu,
        range = range,
        nu = nu,
        sigma = sigma,
        m = m
      )
      if (return_negative_likelihood) {
        return(-loglike)
      } else {
        return(loglike)
      }
    }
  }


  return(loglik)
}

#' @name construct.spde.matern.loglike
#' @title Constructor of Matern loglikelihood functions for non-stationary models.
#' @description This function evaluates the log-likelihood function for observations
#' of a non-stationary Gaussian process defined as the solution to the SPDE
#' \deqn{(\kappa(s) - \Delta)^\beta (\tau(s)u(s)) = W.}
#'
#' The observations are assumed to be generated as
#' \eqn{Y_i = u(s_i) + \epsilon_i}{Y_i = u(s_i) + \epsilon_i}, where
#' \eqn{\epsilon_i}{\epsilon_i} are
#' iid mean-zero Gaussian variables. The latent model is approximated using a
#' rational approximation of the fractional SPDE model.
#'
#' @param object The rational SPDE approximation,
#' computed using [matern.operators()]
#' @param Y The observations, either a vector or a matrix where
#' the columns correspond to independent replicates of observations.
#' @param A An observation matrix that links the measurement location to the
#' finite element basis.
#' @param sigma.e IF non-null, the standard deviation of the measurement noise will be kept fixed in
#' the returned likelihood.
#' @param mu Expectation vector of the latent field (default = 0).
#' @param nu If non-null, the shape parameter will be kept fixed in the returned likelihood.
#' @param m If non-null, update the order of the rational approximation,
#' which needs to be a positive integer.
#' @param log_scale Should the parameters be evaluated in log-scale?
#' @param return_negative_likelihood Return minus the likelihood to turn the maximization into a minimization?
#' @return The log-likelihood function. The parameters of the returned function
#' are given in the order theta, nu, sigma.e, whenever they are available.
#' @export
#' @seealso [matern.operators()], [predict.CBrSPDEobj()]
#' @examples
#' # this example illustrates how the function can be used for maximum
#' # likelihood estimation
#' # Sample a Gaussian Matern process on R using a rational approximation
#' set.seed(123)
#' sigma.e <- 0.1
#' n.rep <- 10
#' n.obs <- 100
#' n.x <- 51
#' # create mass and stiffness matrices for a FEM discretization
#' x <- seq(from = 0, to = 1, length.out = n.x)
#' fem <- rSPDE.fem1d(x)
#' tau <- rep(0.5, n.x)
#' nu <- 0.8
#' alpha <- nu + 0.5
#' kappa <- rep(1, n.x)
#' # Matern parameterization
#' # compute rational approximation
#' op <- spde.matern.operators(
#'   loc_mesh = x,
#'   kappa = kappa, tau = tau, alpha = alpha,
#'   parameterization = "spde", d = 1
#' )
#' # Sample the model
#' u <- simulate(op, n.rep)
#' # Create some data
#' obs.loc <- runif(n = n.obs, min = 0, max = 1)
#' A <- rSPDE.A1d(x, obs.loc)
#' noise <- rnorm(n.obs * n.rep)
#' dim(noise) <- c(n.obs, n.rep)
#' Y <- as.matrix(A %*% u + sigma.e * noise)
#' # define negative likelihood function for optimization using matern.loglike
#' mlik <- construct.spde.matern.loglike(op, Y, A)
#' #' #The parameters can now be estimated by minimizing mlik with optim
#' \donttest{
#' # Choose some reasonable starting values depending on the size of the domain
#' theta0 <- log(c(1 / sqrt(var(c(Y))), sqrt(8), 0.9, 0.01))
#' # run estimation and display the results
#' theta <- optim(theta0, mlik)
#' print(data.frame(
#'   tau = c(tau[1], exp(theta$par[1])), kappa = c(kappa[1], exp(theta$par[2])),
#'   nu = c(nu, exp(theta$par[3])), sigma.e = c(sigma.e, exp(theta$par[4])),
#'   row.names = c("Truth", "Estimates")
#' ))
#' }
#'
#' # SPDE parameterization
#' # compute rational approximation
#' op <- spde.matern.operators(
#'   kappa = kappa, tau = tau, alpha = alpha,
#'   loc_mesh = x, d = 1,
#'   parameterization = "spde"
#' )
#' # Sample the model
#' u <- simulate(op, n.rep)
#' # Create some data
#' obs.loc <- runif(n = n.obs, min = 0, max = 1)
#' A <- rSPDE.A1d(x, obs.loc)
#' noise <- rnorm(n.obs * n.rep)
#' dim(noise) <- c(n.obs, n.rep)
#' Y <- as.matrix(A %*% u + sigma.e * noise)
#' # define negative likelihood function for optimization using matern.loglike
#' mlik <- construct.spde.matern.loglike(op, Y, A)
#' #' #The parameters can now be estimated by minimizing mlik with optim
#' \donttest{
#' # Choose some reasonable starting values depending on the size of the domain
#' theta0 <- log(c(1 / sqrt(var(c(Y))), sqrt(8), 0.9, 0.01))
#' # run estimation and display the results
#' theta <- optim(theta0, mlik)
#' print(data.frame(
#'   tau = c(tau[1], exp(theta$par[1])), kappa = c(kappa[1], exp(theta$par[2])),
#'   nu = c(nu, exp(theta$par[3])), sigma.e = c(sigma.e, exp(theta$par[4])),
#'   row.names = c("Truth", "Estimates")
#' ))
#' }
construct.spde.matern.loglike <- function(object, Y, A,
                                          sigma.e = NULL, mu = 0,
                                          nu = NULL,
                                          m = NULL,
                                          log_scale = TRUE,
                                          return_negative_likelihood = TRUE) {
  loglik <- function(theta) {
    n_tmp <- length(theta)
    if (is.null(nu)) {
      if (log_scale) {
        nu <- exp(theta[n_tmp - 1])
      } else {
        nu <- theta[n_tmp - 1]
      }
      n_tmp <- n_tmp - 1
    } else {
      nu <- nu
    }

    if (is.null(sigma.e)) {
      if (log_scale) {
        sigma.e <- exp(theta[length(theta)])
      } else {
        sigma.e <- theta[length(theta)]
      }
      n_tmp <- n_tmp - 1
    }

    if (nu %% 1 == 0) {
      nu <- nu + 1e-10
    }

    loglike <- spde.matern.loglike(
      object = object, Y = Y, A = A,
      sigma.e = sigma.e,
      mu = mu,
      theta = theta[1:n_tmp],
      nu = nu,
      m = m
    )
    if (return_negative_likelihood) {
      return(-loglike)
    } else {
      return(loglike)
    }
  }


  return(loglik)
}





#' @noRd

aux2_lme_CBrSPDE.matern.loglike <- function(object, y, X_cov, repl, A_list, sigma_e, beta_cov) {
  m <- object$m

  Q <- object$Q

  # R <- tryCatch(Matrix::chol(Matrix::forceSymmetric(Q)), error=function(e){return(NULL)})
  R <- Matrix::Cholesky(Q)

  prior.ld <- c(determinant(R, logarithm = TRUE, sqrt = TRUE)$modulus)

  repl_val <- unique(repl)

  l <- 0

  for (i in repl_val) {
    ind_tmp <- (repl %in% i)
    y_tmp <- y[ind_tmp]

    if (ncol(X_cov) == 0) {
      X_cov_tmp <- 0
    } else {
      X_cov_tmp <- X_cov[ind_tmp, , drop = FALSE]
    }

    na_obs <- is.na(y_tmp)

    y_ <- y_tmp[!na_obs]

    # y_ <- y_list[[as.character(i)]]
    n.o <- length(y_)
    A_tmp <- A_list[[as.character(i)]]
    Q.p <- Q + t(A_tmp) %*% A_tmp / sigma_e^2
    # R.p <- tryCatch(Matrix::chol(Q.p), error=function(e){return(NULL)})
    # if(is.null(R.p)){
    #   return(-10^100)
    # }
    R.p <- Matrix::Cholesky(Q.p)

    posterior.ld <- c(determinant(R.p, logarithm = TRUE, sqrt = TRUE)$modulus)

    # l <- l + sum(log(diag(R))) - sum(log(diag(R.p))) - n.o*log(sigma_e)

    l <- l + prior.ld - posterior.ld - n.o * log(sigma_e)

    v <- y_

    # if(has_cov){
    if (ncol(X_cov) > 0) {
      X_cov_tmp <- X_cov_tmp[!na_obs, , drop = FALSE]
      # X_cov_tmp <- X_cov_list[[as.character(i)]]
      v <- v - X_cov_tmp %*% beta_cov
    }

    # mu.p <- solve(Q.p,as.vector(t(A_tmp) %*% v / sigma_e^2))
    mu.p <- solve(R.p, as.vector(t(A_tmp) %*% v / sigma_e^2), system = "A")

    v <- v - A_tmp %*% mu.p

    l <- l - 0.5 * (t(mu.p) %*% Q %*% mu.p + t(v) %*% v / sigma_e^2) -
      0.5 * n.o * log(2 * pi)
  }

  return(as.double(l))
}

#' @noRd

aux_lme_CBrSPDE.matern.loglike <- function(object, y, X_cov, repl, A_list, sigma_e, beta_cov) {
  l_tmp <- tryCatch(
    aux2_lme_CBrSPDE.matern.loglike(
      object = object,
      y = y, X_cov = X_cov, repl = repl, A_list = A_list,
      sigma_e = sigma_e, beta_cov = beta_cov
    ),
    error = function(e) {
      return(NULL)
    }
  )
  if (is.null(l_tmp)) {
    return(-10^100)
  }
  return(l_tmp)
}



#' @noRd
aux2_lme_rSPDE.matern.loglike <- function(object, y, X_cov, repl, A_list, sigma_e, beta_cov) {
  m <- object$m
  Q <- object$Q

  R <- tryCatch(Matrix::Cholesky(object$Pl), error = function(e) {
    return(NULL)
  })
  if (is.null(R)) {
    return(-10^100)
  }

  prior.ld <- 2 * c(determinant(R, logarithm = TRUE, sqrt = TRUE)$modulus) -
    sum(log(diag(object$C))) / 2

  repl_val <- unique(repl)

  l <- 0

  for (i in repl_val) {
    ind_tmp <- (repl %in% i)
    y_tmp <- y[ind_tmp]

    if (ncol(X_cov) == 0) {
      X_cov_tmp <- 0
    } else {
      X_cov_tmp <- X_cov[ind_tmp, , drop = FALSE]
    }
    na_obs <- is.na(y_tmp)

    y_ <- y_tmp[!na_obs]
    n.o <- length(y_)
    A_tmp <- A_list[[as.character(i)]]
    A_tmp <- A_tmp %*% object$Pr
    Q.p <- Q + t(A_tmp) %*% A_tmp / sigma_e^2
    # R.p <- tryCatch(Matrix::chol(Q.p), error = function(e){return(NULL)})
    # if(is.null(R.p)){
    #   return(-10^100)
    # }

    v <- y_

    if (ncol(X_cov) != 0) {
      X_cov_tmp <- X_cov_tmp[!na_obs, , drop = FALSE]
      v <- v - X_cov_tmp %*% beta_cov
    }

    R.post <- tryCatch(Matrix::Cholesky(Q.p), error = function(e) {
      return(NULL)
    })
    if (is.null(R.post)) {
      return(-10^100)
    }

    AtY <- t(A_tmp) %*% v / sigma_e^2

    mu.p <- solve(R.post, AtY, system = "A")

    v <- v - A_tmp %*% mu.p

    posterior.ld <- c(determinant(R.post, logarithm = TRUE, sqrt = TRUE)$modulus)

    l <- l + prior.ld - posterior.ld - n.o * log(sigma_e)

    l <- l - 0.5 * (t(mu.p) %*% Q %*% mu.p + t(v) %*% v / sigma_e^2) -
      0.5 * n.o * log(2 * pi)
  }
  return(as.double(l))
}


#' @noRd
aux_lme_CBrSPDE.matern2d.loglike <- function(object, y, X_cov, repl, A_list, sigma_e, beta_cov) {
    l_tmp <- tryCatch(
        aux2_lme_rSPDE.matern2d.loglike(
            object = object,
            y = y, X_cov = X_cov, repl = repl, A_list = A_list,
            sigma_e = sigma_e, beta_cov = beta_cov
        ),
        error = function(e) {
            return(NULL)
        }
    )
    if (is.null(l_tmp)) {
        return(-10^100)
    }
    return(l_tmp)
}



#' @noRd
aux2_lme_rSPDE.matern2d.loglike <- function(object, y, X_cov, repl, A_list, sigma_e, beta_cov) {
        m <- object$m
        
        Q <- object$Q
        
        # R <- tryCatch(Matrix::chol(Matrix::forceSymmetric(Q)), error=function(e){return(NULL)})
        R <- Matrix::Cholesky(Q)
        
        prior.ld <- c(determinant(R, logarithm = TRUE, sqrt = TRUE)$modulus)
        
        repl_val <- unique(repl)
        
        l <- 0
        
        for (i in repl_val) {
            ind_tmp <- (repl %in% i)
            y_tmp <- y[ind_tmp]
            
            if (ncol(X_cov) == 0) {
                X_cov_tmp <- 0
            } else {
                X_cov_tmp <- X_cov[ind_tmp, , drop = FALSE]
            }
            
            na_obs <- is.na(y_tmp)
            
            y_ <- y_tmp[!na_obs]
            
            # y_ <- y_list[[as.character(i)]]
            n.o <- length(y_)
            A_tmp <- A_list[[as.character(i)]]
            Q.p <- Q + t(A_tmp) %*% A_tmp / sigma_e^2
            # R.p <- tryCatch(Matrix::chol(Q.p), error=function(e){return(NULL)})
            # if(is.null(R.p)){
            #   return(-10^100)
            # }
            R.p <- Matrix::Cholesky(Q.p)
            
            posterior.ld <- c(determinant(R.p, logarithm = TRUE, sqrt = TRUE)$modulus)
            
            # l <- l + sum(log(diag(R))) - sum(log(diag(R.p))) - n.o*log(sigma_e)
            
            l <- l + prior.ld - posterior.ld - n.o * log(sigma_e)
            
            v <- y_
            
            # if(has_cov){
            if (ncol(X_cov) > 0) {
                X_cov_tmp <- X_cov_tmp[!na_obs, , drop = FALSE]
                # X_cov_tmp <- X_cov_list[[as.character(i)]]
                v <- v - X_cov_tmp %*% beta_cov
            }
            
            # mu.p <- solve(Q.p,as.vector(t(A_tmp) %*% v / sigma_e^2))
            mu.p <- solve(R.p, as.vector(t(A_tmp) %*% v / sigma_e^2), system = "A")
            
            v <- v - A_tmp %*% mu.p
            
            l <- l - 0.5 * (t(mu.p) %*% Q %*% mu.p + t(v) %*% v / sigma_e^2) -
                0.5 * n.o * log(2 * pi)
        }
        
        return(as.double(l))
    }

#' @noRd

aux_lme_rSPDE.matern.loglike <- function(object, y, X_cov, repl, A_list, sigma_e, beta_cov) {
  l_tmp <- tryCatch(
    aux2_lme_rSPDE.matern.loglike(
      object = object,
      y = y, X_cov = X_cov, repl = repl, A_list = A_list,
      sigma_e = sigma_e, beta_cov = beta_cov
    ),
    error = function(e) {
      return(NULL)
    }
  )
  if (is.null(l_tmp)) {
    return(-10^100)
  }
  return(l_tmp)
}
Any scripts or data that you put into this service are public.
rSPDE documentation built on Jan. 26, 2026, 9:06 a.m.
rdrr.io home R language documentation Run R code online
CRAN packages Bioconductor packages R-Forge packages GitHub packages
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
rSPDE
Rational Approximations of Fractional Stochastic Partial Differential Equations

R/fractional.computations.R
In rSPDE: Rational Approximations of Fractional Stochastic Partial Differential Equations

Try the rSPDE package in your browser

R Package Documentation

Browse R Packages

We want your feedback!

rSPDE Rational Approximations of Fractional Stochastic Partial Differential Equations

R/fractional.computations.R In rSPDE: Rational Approximations of Fractional Stochastic Partial Differential Equations

Try the rSPDE package in your browser

R Package Documentation

Browse R Packages

We want your feedback!

rSPDE
Rational Approximations of Fractional Stochastic Partial Differential Equations

R/fractional.computations.R
In rSPDE: Rational Approximations of Fractional Stochastic Partial Differential Equations
