R/util.R

Defines functions select_indexes idx_not_any_NA idx_not_all_NA return_same_input_type_matrix_vector change_parameterization_lme convert_B_matrices likelihood_process_parameters likelihood_process_inputs_matern likelihood_process_inputs_spde rspde_check_user_input rational.order `rational.order<-` rational.type `rational.type<-` get_rational_coefficients print.CBrSPDEobj print.summary.CBrSPDEobj summary.CBrSPDEobj create_summary_from_density symmetric_part_matrix analyze_sparsity_rspde build_sparse_matrix_rspde get.sparsity.graph.rspde generic_fem_mesh_order fem_mesh_order_1d get_inla_mesh_dimension check_class_inla_rspde cut_decimals get.initial.values.rSPDE require.nowarnings rSPDE.fem2d rSPDE.fem1d rSPDE.A1d print.rSPDEobj print.summary.rSPDEobj summary.rSPDEobj folded.matern.covariance.2d folded.matern.covariance.1d matern.covariance get.roots

Documented in folded.matern.covariance.1d folded.matern.covariance.2d get.initial.values.rSPDE matern.covariance print.CBrSPDEobj print.rSPDEobj print.summary.CBrSPDEobj print.summary.rSPDEobj rational.order rational.type require.nowarnings rSPDE.A1d rSPDE.fem1d rSPDE.fem2d summary.CBrSPDEobj summary.rSPDEobj

## util.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/>.


# Internal function to get the roots of the polynomials used
# in the rational approximation.
get.roots <- function(m, beta) {
  if (beta > 2) {
    beta <- beta - floor(beta - 1)
  }

  rb <- rep(0, m + 1)
  rc <- rep(0, m)
  if (m == 1) {
    rc <- approx(m1table$beta, m1table$rc, beta)$y
    rb[1] <- approx(m1table$beta, m1table$rb.1, beta)$y
    rb[2] <- approx(m1table$beta, m1table$rb.2, beta)$y
    factor <- approx(m1table$beta, m1table$factor, beta)$y
  } else if (m == 2) {
    rc[1] <- approx(m2table$beta, m2table$rc.1, beta)$y
    rc[2] <- approx(m2table$beta, m2table$rc.2, beta)$y
    rb[1] <- approx(m2table$beta, m2table$rb.1, beta)$y
    rb[2] <- approx(m2table$beta, m2table$rb.2, beta)$y
    rb[3] <- approx(m2table$beta, m2table$rb.3, beta)$y
    factor <- approx(m2table$beta, m2table$factor, beta)$y
  } else if (m == 3) {
    rc[1] <- approx(m3table$beta, m3table$rc.1, beta)$y
    rc[2] <- approx(m3table$beta, m3table$rc.2, beta)$y
    rc[3] <- approx(m3table$beta, m3table$rc.3, beta)$y
    rb[1] <- approx(m3table$beta, m3table$rb.1, beta)$y
    rb[2] <- approx(m3table$beta, m3table$rb.2, beta)$y
    rb[3] <- approx(m3table$beta, m3table$rb.3, beta)$y
    rb[4] <- approx(m3table$beta, m3table$rb.4, beta)$y
    factor <- approx(m3table$beta, m3table$factor, beta)$y
  } else if (m == 4) {
    rc[1] <- approx(m4table$beta, m4table$rc.1, beta)$y
    rc[2] <- approx(m4table$beta, m4table$rc.2, beta)$y
    rc[3] <- approx(m4table$beta, m4table$rc.3, beta)$y
    rc[4] <- approx(m4table$beta, m4table$rc.4, beta)$y
    rb[1] <- approx(m4table$beta, m4table$rb.1, beta)$y
    rb[2] <- approx(m4table$beta, m4table$rb.2, beta)$y
    rb[3] <- approx(m4table$beta, m4table$rb.3, beta)$y
    rb[4] <- approx(m4table$beta, m4table$rb.4, beta)$y
    rb[5] <- approx(m4table$beta, m4table$rb.5, beta)$y
    factor <- approx(m4table$beta, m4table$factor, beta)$y
  } else {
    stop("m must be one of the values 1,2,3,4.")
  }

  return(list(rb = rb, rc = rc, factor = factor))
}

#' The Matern covariance function
#'
#' `matern.covariance` evaluates the Matern covariance function
#' \deqn{C(h) = \frac{\sigma^2}{2^{\nu-1}\Gamma(\nu)}(\kappa h)^\nu
#' K_\nu(\kappa h).}
#'
#' @param h Distances to evaluate the covariance function at.
#' @param kappa Range parameter.
#' @param nu Shape parameter.
#' @param sigma Standard deviation.
#'
#' @return A vector with the values C(h).
#' @export
#'
#' @examples
#' x <- seq(from = 0, to = 1, length.out = 101)
#' plot(x, matern.covariance(abs(x - 0.5), kappa = 10, nu = 1 / 5, sigma = 1),
#'   type = "l", ylab = "C(h)", xlab = "h"
#' )
#'
matern.covariance <- function(h,
                              kappa,
                              nu,
                              sigma) {
  if (nu == 1 / 2) {
    C <- sigma^2 * exp(-kappa * abs(h))
  } else {
    C <- (sigma^2 / (2^(nu - 1) * gamma(nu))) *
    ((kappa * abs(h))^nu) * besselK(kappa * abs(h), nu)
  }
  C[h == 0] <- sigma^2
  return(as.matrix(C))
}

#' The 1d folded Matern covariance function
#'
#' @description
#' `folded.matern.covariance.1d` evaluates the 1d
#' folded Matern covariance function over an interval \eqn{[0,L]}.
#'
#' @details
#' `folded.matern.covariance.1d` evaluates the 1d folded Matern
#' covariance function over an interval \eqn{[0,L]} under different
#' boundary conditions. For periodic boundary conditions
#' \deqn{C_{\mathcal{P}}(h,m) = \sum_{k=-\infty}^{\infty} (C(h-m+2kL),}
#' for Neumann boundary conditions
#' \deqn{C_{\mathcal{N}}(h,m) = \sum_{k=-\infty}^{\infty}
#' (C(h-m+2kL)+C(h+m+2kL)),}
#' and for Dirichlet boundary conditions:
#' \deqn{C_{\mathcal{D}}(h,m) = \sum_{k=-\infty}^{\infty}
#' (C(h-m+2kL)-C(h+m+2kL)),}
#' where \eqn{C(\cdot)} is the Matern covariance function:
#' \deqn{C(h) = \frac{\sigma^2}{2^{\nu-1}\Gamma(\nu)}(\kappa h)^\nu
#' K_\nu(\kappa h).}
#'
#' We consider the truncation:
#' \deqn{C_{{\mathcal{P}},N}(h,m) = \sum_{k=-N}^{N} C(h-m+2kL),
#' C_{\mathcal{N},N}(h,m) = \sum_{k=-\infty}^{\infty}
#' (C(h-m+2kL)+C(h+m+2kL)),}
#' and
#' \deqn{C_{\mathcal{D},N}(h,m) = \sum_{k=-N}^{N}
#' (C(h-m+2kL)-C(h+m+2kL)).}
#'
#' @param h,m Vectors of arguments of the covariance function.
#' @param kappa Range parameter.
#' @param nu Shape parameter.
#' @param sigma Standard deviation.
#' @param L The upper bound of the interval \eqn{[0,L]}. By default, `L=1`.
#' @param N The truncation parameter.
#' @param boundary The boundary condition. The possible conditions
#' are `"neumann"` (default), `"dirichlet"` or
#' `"periodic"`.
#'
#' @return A matrix with the corresponding covariance values.
#' @export
#'
#' @examples
#' x <- seq(from = 0, to = 1, length.out = 101)
#' plot(x, folded.matern.covariance.1d(rep(0.5, length(x)), x,
#' kappa = 10, nu = 1 / 5, sigma = 1),
#'   type = "l", ylab = "C(h)", xlab = "h"
#' )
#'
folded.matern.covariance.1d <- function(h, m, kappa, nu, sigma,
                                        L = 1, N = 10,
                                        boundary = c("neumann",
                                        "dirichlet", "periodic")) {
  boundary <- tolower(boundary[1])
  if (!(boundary %in% c("neumann", "dirichlet", "periodic"))) {
    stop("The possible boundary conditions are 'neumann',
    'dirichlet' or 'periodic'!")
  }
  if (length(h) != length(m)) {
    stop("h and m should have the same length!")
  }

  s1 <- sapply(-N:N, function(j) {
    h - m + 2 * j * L
  })
  s2 <- sapply(-N:N, function(j) {
    h + m + 2 * j * L
  })
  if (boundary == "neumann") {
    C <- rowSums(matern.covariance(h = s1, kappa = kappa,
    nu = nu, sigma = sigma) +
      matern.covariance(h = s2, kappa = kappa,
      nu = nu, sigma = sigma))
  } else if (boundary == "dirichlet") {
    C <- rowSums(matern.covariance(h = s1, kappa = kappa,
    nu = nu, sigma = sigma) -
      matern.covariance(h = s2, kappa = kappa,
      nu = nu, sigma = sigma))
  } else {
    C <- rowSums(matern.covariance(h = s1,
    kappa = kappa, nu = nu, sigma = sigma))
  }

  if (length(h) == 1) {
    return(sum(C))
  }
  return(as.matrix(C))
}

#' The 2d folded Matern covariance function
#'
#' @description
#' `folded.matern.covariance.2d` evaluates the 2d
#' folded Matern covariance function over an interval
#' \eqn{[0,L]\times [0,L]}.
#'
#' @details
#' `folded.matern.covariance.2d` evaluates the 1d folded
#' Matern covariance function over an interval
#' \eqn{[0,L]\times [0,L]} under different boundary conditions.
#' For periodic boundary conditions
#' \deqn{C_{\mathcal{P}}((h_1,h_2),(m_1,m_2)) =
#' \sum_{k_2=-\infty}^\infty \sum_{k_1=-\infty}^{\infty}
#' (C(\|(h_1-m_1+2k_1L,h_2-m_2+2k_2L)\|),}
#' for Neumann boundary conditions
#' \deqn{C_{\mathcal{N}}((h_1,h_2),(m_1,m_2)) =
#' \sum_{k_2=-\infty}^\infty \sum_{k_1=-\infty}^{\infty}
#' (C(\|(h_1-m_1+2k_1L,h_2-m_2+2k_2L)\|)+C(\|(h_1-m_1+2k_1L,
#' h_2+m_2+2k_2L)\|)+C(\|(h_1+m_1+2k_1L,h_2-m_2+2k_2L)\|)+
#' C(\|(h_1+m_1+2k_1L,h_2+m_2+2k_2L)\|)),}
#' and for Dirichlet boundary conditions:
#' \deqn{C_{\mathcal{D}}((h_1,h_2),(m_1,m_2)) = \sum_{k_2=-\infty}^\infty
#' \sum_{k_1=-\infty}^{\infty} (C(\|(h_1-m_1+2k_1L,h_2-m_2+2k_2L)\|)-
#' C(\|(h_1-m_1+2k_1L,h_2+m_2+2k_2L)\|)-C(\|(h_1+m_1+2k_1L,
#' h_2-m_2+2k_2L)\|)+C(\|(h_1+m_1+2k_1L,h_2+m_2+2k_2L)\|)),}
#' where \eqn{C(\cdot)} is the Matern covariance function:
#' \deqn{C(h) = \frac{\sigma^2}{2^{\nu-1}\Gamma(\nu)}
#' (\kappa h)^\nu K_\nu(\kappa h).}
#'
#' We consider the truncation for \eqn{k_1,k_2} from \eqn{-N} to \eqn{N}.
#'
#' @param h,m Vectors with two coordinates.
#' @param kappa Range parameter.
#' @param nu Shape parameter.
#' @param sigma Standard deviation.
#' @param L The upper bound of the square \eqn{[0,L]\times [0,L]}.
#' By default, `L=1`.
#' @param N The truncation parameter.
#' @param boundary The boundary condition. The possible conditions
#' are `"neumann"` (default), `"dirichlet"`,
#' `"periodic"` or `"R2"`.
#'
#' @return The correspoding covariance.
#' @export
#'
#' @examples
#' h <- c(0.5, 0.5)
#' m <- c(0.5, 0.5)
#' folded.matern.covariance.2d(h, m, kappa = 10, nu = 1 / 5, sigma = 1)
#'
folded.matern.covariance.2d <- function(h, m, kappa, nu, sigma,
                                        L = 1, N = 10,
                                        boundary = c("neumann",
                                        "dirichlet",
                                        "periodic",
                                        "R2")) {
  boundary <- tolower(boundary[1])
  if (!(boundary %in% c("neumann", "dirichlet",
  "periodic", "r2"))) {
    stop("The possible boundary conditions are
    'neumann', 'dirichlet', 'periodic' or 'R2'!")
  }

  if (is.vector(h)) {
    if (!is.vector(m)) {
      stop("since 'h' is a vector, 'm' should be a vector!")
    }

    if ((length(h) != 2) || (length(m) != 2)) {
      stop("The vectors h and m should have length 2!")
    }
  } else if (is.matrix(h)) {
    if (!is.matrix(m)) {
      stop("since 'h' is a matrix, 'm' should be a matrix!")
    }
    if (ncol(h) != 2) {
      stop("h must have two columns!")
    }
    if (!all(dim(h) == dim(m))) {
      stop("h and m must have the same dimensions!")
    }
  }

  if (!is.vector(h) && !is.matrix(h)) {
    stop("h should be either a vector or a matrix!")
  }

  list.comb <- expand.grid(-N:N, -N:N)

  if (is.matrix(h)) {
    h_matrix_1 <- matrix(rep(h[, 1], length(list.comb[, 1])), nrow = nrow(h))
    h_matrix_2 <- matrix(rep(h[, 2], length(list.comb[, 1])), nrow = nrow(h))
    m_matrix_1 <- matrix(rep(m[, 1], length(list.comb[, 1])), nrow = nrow(m))
    m_matrix_2 <- matrix(rep(m[, 2], length(list.comb[, 1])), nrow = nrow(m))
    list_comb_matrix_1 <- t(matrix(rep(list.comb[, 1],
    nrow(h)), ncol = nrow(h)))
    list_comb_matrix_2 <- t(matrix(rep(list.comb[, 2],
    nrow(h)), ncol = nrow(h)))

    s11 <- sqrt((h_matrix_1 - m_matrix_1 + 2 *
    list_comb_matrix_1 * L)^2 + (h_matrix_2 -
    m_matrix_2 + 2 * list_comb_matrix_2 * L)^2)
    s12 <- sqrt((h_matrix_1 - m_matrix_1 + 2 *
    list_comb_matrix_1 * L)^2 + (h_matrix_2 +
    m_matrix_2 + 2 * list_comb_matrix_2 * L)^2)
    s21 <- sqrt((h_matrix_1 + m_matrix_1 + 2 *
    list_comb_matrix_1 * L)^2 + (h_matrix_2 -
    m_matrix_2 + 2 * list_comb_matrix_2 * L)^2)
    s22 <- sqrt((h_matrix_1 + m_matrix_1 + 2 *
    list_comb_matrix_1 * L)^2 + (h_matrix_2 +
    m_matrix_2 + 2 * list_comb_matrix_2 * L)^2)

    if (boundary == "neumann") {
      C <- rowSums(matern.covariance(h = s11,
      kappa = kappa, nu = nu, sigma = sigma) +
        matern.covariance(h = s12, kappa = kappa,
        nu = nu, sigma = sigma) +
        matern.covariance(h = s21, kappa = kappa,
        nu = nu, sigma = sigma) +
        matern.covariance(h = s22, kappa = kappa,
        nu = nu, sigma = sigma))
    } else if (boundary == "dirichlet") {
      C <- rowSums(matern.covariance(h = s11,
      kappa = kappa, nu = nu, sigma = sigma) -
        matern.covariance(h = s12, kappa = kappa,
        nu = nu, sigma = sigma) -
        matern.covariance(h = s21, kappa = kappa,
        nu = nu, sigma = sigma) +
        matern.covariance(h = s22, kappa = kappa,
        nu = nu, sigma = sigma))
    } else if (boundary == "r2") {
      C <- matern.covariance(h = sqrt((h[1] - m[1])^2 +
      (h[2] - m[2])^2), kappa = kappa, sigma = sigma,
      nu = nu)
    } else {
      C <- rowSums(matern.covariance(h = s11,
      kappa = kappa, nu = nu, sigma = sigma))
    }
    return(as.double(C))
  } else {
    s11 <- sqrt((h[1] - m[1] + 2 * list.comb[, 1] * L)^2 +
    (h[2] - m[2] + 2 * list.comb[, 2] * L)^2)
    s12 <- sqrt((h[1] - m[1] + 2 * list.comb[, 1] * L)^2 +
    (h[2] + m[2] + 2 * list.comb[, 2] * L)^2)
    s21 <- sqrt((h[1] + m[1] + 2 * list.comb[, 1] * L)^2 +
    (h[2] - m[2] + 2 * list.comb[, 2] * L)^2)
    s22 <- sqrt((h[1] + m[1] + 2 * list.comb[, 1] * L)^2 +
    (h[2] + m[2] + 2 * list.comb[, 2] * L)^2)

    if (boundary == "neumann") {
      C <- sum(matern.covariance(h = s11, kappa = kappa,
      nu = nu, sigma = sigma) +
        matern.covariance(h = s12, kappa = kappa,
        nu = nu, sigma = sigma) +
        matern.covariance(h = s21, kappa = kappa,
        nu = nu, sigma = sigma) +
        matern.covariance(h = s22, kappa = kappa,
        nu = nu, sigma = sigma))
    } else if (boundary == "dirichlet") {
      C <- sum(matern.covariance(h = s11,
      kappa = kappa, nu = nu, sigma = sigma) -
        matern.covariance(h = s12, kappa = kappa,
        nu = nu, sigma = sigma) -
        matern.covariance(h = s21, kappa = kappa,
        nu = nu, sigma = sigma) +
        matern.covariance(h = s22, kappa = kappa,
        nu = nu, sigma = sigma))
    } else if (boundary == "r2") {
      C <- matern.covariance(h = sqrt((h[1] - m[1])^2 +
      (h[2] - m[2])^2), kappa = kappa, sigma = sigma,
      nu = nu)
    } else {
      C <- sum(matern.covariance(h = s11,
      kappa = kappa, nu = nu, sigma = sigma))
    }
    return(as.double(C))
  }
}




#' Summarise rSPDE objects
#'
#' Summary method for class "rSPDEobj"
#'
#' @param object an object of class "rSPDEobj", usually, a result of a call
#'   to [fractional.operators()], [matern.operators()], or
#'   [spde.matern.operators()].
#' @param ... further arguments passed to or from other methods.
#' @export
#' @method summary rSPDEobj
summary.rSPDEobj <- function(object, ...) {
  out <- list()
  class(out) <- "summary.rSPDEobj"
  out$type <- object$type
  if (out$type == "Matern approximation") {
    out$kappa <- object$kappa
    out$sigma <- object$sigma
    out$tau <- object$tau
    out[["range"]] <- object[["range"]]
    out$nu <- object$nu
  }
  out$m <- object$m
  out$stationary <- object$stationary
  out$parameterization <- object$parameterization
  out$n <- dim(object$L)[1]
  return(out)
}

#' @param x an object of class "summary.rSPDEobj", usually, a result of a call
#'   to [summary.rSPDEobj()].
#' @export
#' @method print summary.rSPDEobj
#' @rdname summary.rSPDEobj
print.summary.rSPDEobj <- function(x, ...) {
  cat("Type of approximation: ", x$type, "\n")
  cat("Parameterization: ", x$parameterization, "\n")  
  if (x$type == "Matern approximation") {
    if(x$parameterization == "spde"){
    cat(
      "Parameters of covariance function: kappa = ",
      x$kappa, ", tau = ", x$tau, ", nu = ", x$nu, "\n"
    )
    } else{
      cat(
      "Parameters of covariance function: range = ",
      x[["range"]], ", sigma = ", x$sigma, ", nu = ", x$nu, "\n"
    )
    }
  }
  cat("Order or rational approximation: ", x$m, "\n")
  cat("Size of discrete operators: ", x$n, " x ", x$n, "\n")
  if(x$stationary){
    cat("Stationary Model\n")
  } else{
    cat("Non-Stationary Model")
  }
}

#' @export
#' @method print rSPDEobj
#' @rdname summary.rSPDEobj
print.rSPDEobj <- function(x, ...) {
  print.summary.rSPDEobj(summary(x))
}

#' Observation matrix for finite element discretization on R
#'
#' A finite element discretization on R can be written as
#' \eqn{u(s) = \sum_i^n u_i \varphi_i(s)}{u(s) = \sum_i^n u_i \varphi_i(s)}
#' where \eqn{\varphi_i(s)} is a piecewise linear
#' "hat function" centered at location
#' \eqn{x_i}{x_i}. This function computes an
#' \eqn{m\times n}{m x n} matrix \eqn{A}{A}
#' that links the basis function in the expansion to specified locations
#' \eqn{s = (s_1,\ldots, s_m)} in the domain through
#' \eqn{A_ij = \varphi_j(s_i)}{A_ij = \varphi_j(s_i)}.
#'
#' @param x The locations of the nodes in the FEM discretization.
#' @param loc The locations \eqn{(s_1,\ldots, s_m)}
#'
#' @return The sparse matrix `A`.
#' @export
#' @author David Bolin \email{davidbolin@@gmail.com}
#' @seealso [rSPDE.fem1d()]
#'
#' @examples
#' # create mass and stiffness matrices for a FEM discretization on [0,1]
#' x <- seq(from = 0, to = 1, length.out = 101)
#' fem <- rSPDE.fem1d(x)
#'
#' # create the observation matrix for some locations in the domain
#' obs.loc <- runif(n = 10, min = 0, max = 1)
#' A <- rSPDE.A1d(x, obs.loc)
rSPDE.A1d <- function(x, loc) {
  if (min(loc) < min(x) || max(loc) > max(x)) {
    stop("locations outside support of basis")
  }

  n.x <- length(x)
  n.loc <- length(loc)
  i <- as.vector(cbind(1:n.loc, 1:n.loc))
  j <- matrix(0, n.loc, 2)
  vals <- matrix(1, n.loc, 2)
  for (ii in seq_len(n.loc)) {
    j[ii, 1] <- sum(sum((loc[ii] - x) >= 0))
    vals[ii, 1] <- loc[ii] - x[j[ii, 1]]
    j[ii, 2] <- j[ii, 1] + 1
    if (j[ii, 2] <= n.x) {
      vals[ii, 2] <- x[j[ii, 2]] - loc[ii]
    } else {
      j[ii, 2] <- j[ii, 2] - 2
    }
  }
  j <- as.vector(j)
  vals <- as.vector(matrix(1 - vals / rowSums(vals)))

  A <- sparseMatrix(i = i, j = j, x = vals, dims = c(n.loc, n.x))
  return(A)
}


#' Finite element calculations for problems on R
#'
#' This function computes mass and stiffness matrices
#' for a FEM approximation on R, assuming
#' Neumann boundary conditions.
#' These matrices are needed when discretizing the
#' operators in rational approximations.
#'
#' @param x Locations of the nodes in the FEM approximation.
#'
#' @return The function returns a list with the following elements
#' \item{G }{The stiffness matrix.}
#' \item{C }{The mass matrix.}
#' @export
#' @author David Bolin \email{davidbolin@@gmail.com}
#' @seealso [rSPDE.A1d()]
#' @examples
#' # create mass and stiffness matrices for a FEM discretization on [0,1]
#' x <- seq(from = 0, to = 1, length.out = 101)
#' fem <- rSPDE.fem1d(x)
rSPDE.fem1d <- function(x) {
  n <- length(x)
  d <- c(Inf, diff(x))
  dm1 <- c(d[2:n], Inf)
  G <- -bandSparse(
    n = n, m = n, k = c(-1, 0, 1),
    diagonals = cbind(1 / dm1, -(1 / dm1 + 1 / d), 1 / dm1)
  )
  C <- bandSparse(
    n = n, m = n, k = c(-1, 0, 1),
    diagonals = cbind(dm1 / 6, (dm1 + d) / 3, c(d[2:n],Inf) / 6)
  )
  C[1, 1:2] <- c(d[2], d[2] / 2) / 3
  C[n, (n - 1):n] <- c(d[n] / 2, d[n]) / 3

  return(list(G = G, C = C))
}

#' Finite element calculations for problems in 2D
#'
#' This function computes mass and stiffness matrices for a mesh in 2D, assuming
#' Neumann boundary conditions.
#'
#' @param FV Matrix where each row defines a triangle
#' @param P Locations of the nodes in the mesh.
#'
#' @return The function returns a list with the following elements
#' \item{G }{The stiffness matrix with elements \eqn{(\nabla \phi_i, \nabla \phi_j)}.}
#' \item{C }{The mass matrix with elements \eqn{(\phi_i, \phi_j)}.}
#' \item{Cd }{The mass lumped matrix with diagonal elements \eqn{(\phi_i, 1)}.}
#' \item{Hxx }{Matrix with elements \eqn{(\partial_x \phi_i, \partial_x \phi_j)}.}
#' \item{Hyy }{Matrix with elements \eqn{(\partial_y \phi_i, \partial_y \phi_j)}.}
#' \item{Hxy }{Matrix with elements \eqn{(\partial_x \phi_i, \partial_y \phi_j)}.}
#' \item{Hyx }{Matrix with elements \eqn{(\partial_y \phi_i, \partial_x \phi_j)}.}
#' @export
#' @author David Bolin \email{davidbolin@@gmail.com}
#' @seealso [rSPDE.fem1d()]
#' @examples
#' P <- rbind(c(0,0), c(1,0), c(1,1), c(0,1))
#' FV <- rbind(c(1,2,3), c(2,3,4))
#' fem <- rSPDE.fem2d(FV,P)
rSPDE.fem2d <- function(FV, P) {
  
  d <- ncol(FV) - 1
  if(d != 2){stop("Only 2d supported")}
  if(ncol(P) != d){P <- t(P)}
  if(ncol(P) != d){stop("Wrong dimension of P")}
  
  nV <- nrow(P)
  nF <- nrow(FV)
  Gi <- matrix(0, nrow = nF*3, ncol = 3)
  Gj <- Gz <- Ci <- Cj <- Cz <- Gxx <- Gxy <- Gyx <- Gyy <- Gi
  
  Mxx <- matrix(c(1, -1, 0,-1, 1, 0, 0, 0, 0), 3, 3)   
  Myy <- matrix(c(1, 0, -1, 0, 0, 0, -1, 0, 1), 3, 3)
  Mxy <- matrix(c(1, -1, 0, 0, 0, 0, -1, 1, 0), 3, 3)
  Myx <- matrix(c(1, 0, -1, -1, 0, 1, 0, 0, 0), 3, 3)
  for(f in 1:nF) {
    dd <- 3*(f-1)+(1:3)
    Gi[dd, ] <- Ci[dd, ] <- FV[f, ] %*% t(rep(1,3))
    Gj[dd, ] <- Cj[dd, ] <- t(Gi[dd, ])
    
    xy <- t(P[FV[f, ], ])
    m1 <- rbind(rep(1, 3), xy)
    m2 <- rbind(rep(0, 2), diag(1, 2))
    m <- solve(m1, m2)
    ddet <- abs(det(m1))
    Gz[dd,] <- ddet * (m %*% t(m)) / 2
    Cz[dd,] <- ddet * (rep(1, 3) + diag(3)) / 24
    
    Bk <- matrix(c(xy[1, 2] - xy[1, 1], 
                   xy[2, 2] - xy[2, 1],
                   xy[1, 3] - xy[1, 1], 
                   xy[2, 3] - xy[2, 1]), 2, 2)
    
    Bki <- solve(Bk)
    Cxx <- Bki%*%matrix(c(1, 0, 0, 0), 2, 2)%*%t(Bki)
    Cyy <- Bki%*%matrix(c(0, 0, 0, 1), 2, 2)%*%t(Bki)
    Cxy <- Bki%*%matrix(c(0, 0, 1, 0), 2, 2)%*%t(Bki)
    Cyx <- Bki%*%matrix(c(0, 1, 0, 0), 2, 2)%*%t(Bki)
    
    Gxx[dd, ] <- ddet*(Cxx[1, 1]*Mxx + Cxx[1, 2]*Mxy + Cxx[2, 1]*Myx + Cxx[2, 2]*Myy) / 2
    Gyy[dd, ] <- ddet*(Cyy[1, 1]*Mxx + Cyy[1, 2]*Mxy + Cyy[2, 1]*Myx + Cyy[2, 2]*Myy) / 2
    Gxy[dd, ] <- ddet*(Cxy[1, 1]*Mxx + Cxy[1, 2]*Mxy + Cxy[2, 1]*Myx + Cxy[2, 2]*Myy) / 2
    Gyx[dd, ] <- ddet*(Cyx[1, 1]*Mxx + Cyx[1, 2]*Mxy + Cyx[2, 1]*Myx + Cyx[2, 2]*Myy) / 2
  }
  
  G <- Matrix::sparseMatrix(i = as.vector(Gi), j = as.vector(Gj), 
                            x = as.vector(Gz), dims = c(nV,nV))
  Hxx <- Matrix::sparseMatrix(i = as.vector(Gi), j = as.vector(Gj), 
                              x = as.vector(Gxx), dims = c(nV,nV))
  Hyy <- Matrix::sparseMatrix(i = as.vector(Gi), j = as.vector(Gj), 
                              x = as.vector(Gyy), dims = c(nV,nV))
  Hxy <- Matrix::sparseMatrix(i = as.vector(Gi), j = as.vector(Gj), 
                              x = as.vector(Gxy), dims = c(nV,nV))
  Hyx <- Matrix::sparseMatrix(i = as.vector(Gi), j = as.vector(Gj), 
                              x = as.vector(Gyx), dims = c(nV,nV))
  Ce <- Matrix::sparseMatrix(i = as.vector(Ci), j = as.vector(Cj), 
                             x = as.vector(Cz), dims = c(nV,nV))
  C <- Matrix::Diagonal(n = nV, x = Matrix::colSums(Ce))
  return(list(G = G, C = Ce, Cd = C, Hxx = Hxx, Hyy = Hyy, Hxy = Hxy, Hyx = Hyx))
}
#' Warnings free loading of add-on packages
#'
#' Turn off all warnings for require(), to allow clean completion
#' of examples that require unavailable Suggested packages.
#'
#' @param package The name of a package, given as a character string.
#' @param lib.loc a character vector describing the location of R library trees
#' to search through, or `NULL`.  The default value of `NULL`
#' corresponds to all libraries currently known to `.libPaths()`.
#' Non-existent library trees are silently ignored.
#' @param character.only a logical indicating whether `package` can be
#' assumed to be a character string.
#'
#' @return `require.nowarnings` returns (invisibly)
#' `TRUE` if it succeeds, otherwise `FALSE`
#' @details `require(package)` acts the same as
#' `require(package, quietly = TRUE)` but with warnings turned off.
#' In particular, no warning or error is given if the package is unavailable.
#' Most cases should use `requireNamespace(package,
#' quietly = TRUE)` instead,
#' which doesn't produce warnings.
#' @seealso [require()]
#' @export
#' @examples
#' ## This should produce no output:
#' if (require.nowarnings(nonexistent)) {
#'   message("Package loaded successfully")
#' }
#'
require.nowarnings <- function(package, lib.loc = NULL,
character.only = FALSE) {
  if (!character.only) {
    package <- as.character(substitute(package))
  }
  suppressWarnings(
    require(package,
      lib.loc = lib.loc,
      quietly = TRUE,
      character.only = TRUE
    )
  )
}

#' @name get.initial.values.rSPDE
#' @title Initial values for log-likelihood optimization in rSPDE models
#' with a latent stationary Gaussian Matern model
#' @description Auxiliar function to obtain domain-based initial values for
#' log-likelihood optimization in rSPDE models
#' with a latent stationary Gaussian Matern model
#' @param mesh An in INLA mesh
#' @param mesh.range The range of the mesh.
#' @param graph.obj A `metric_graph` object. To be used in case both `mesh` and `mesh.range` are `NULL`.
#' @param dim The dimension of the domain.
#' @param B.sigma Matrix with specification of log-linear model for \eqn{\sigma}. Will be used if `parameterization = 'matern'`.
#' @param B.range Matrix with specification of log-linear model for \eqn{\rho}, which is a range-like parameter (it is exactly the range parameter in the stationary case). Will be used if `parameterization = 'matern'`.
#' @param parameterization Which parameterization to use? `matern` uses range, std. deviation and nu (smoothness). `spde` uses kappa, tau and nu (smoothness). The default is `matern`.
#' @param B.tau Matrix with specification of log-linear model for \eqn{\tau}. Will be used if `parameterization = 'spde'`.
#' @param B.kappa Matrix with specification of log-linear model for \eqn{\kappa}. Will be used if `parameterization = 'spde'`.
#' @param nu The smoothness parameter.
#' @param include.nu Should we also provide an initial guess for nu?
#' @param n.spde The number of basis functions in the mesh model.
#' @param log.scale Should the results be provided in log scale?
#' @param nu.upper.bound Should an upper bound for nu be considered?
#' @return A vector of the form (theta_1,theta_2,theta_3) or where
#' theta_1 is the initial guess for tau, theta_2 is the initial guess for kappa
#' and theta_3 is the initial guess for nu.
#' @export
#'

get.initial.values.rSPDE <- function(mesh = NULL, mesh.range = NULL,
                                     graph.obj = NULL,
                                    n.spde = 1,
                                    dim = NULL, B.tau = NULL, B.kappa = NULL,
                                     B.sigma = NULL, B.range = NULL, nu = NULL,
                                    parameterization = c("matern", "spde"),
                                    include.nu = TRUE, log.scale = TRUE,
                                    nu.upper.bound = NULL) {
  if (is.null(mesh) && is.null(mesh.range) && is.null(graph.obj)) {
    stop("You should either provide mesh, mesh.range or graph_obj!")
  }

    parameterization <- parameterization[[1]]

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

  if (is.null(mesh) && is.null(graph.obj) && is.null(dim)) {
    stop("If you don't provide mesh, you have to provide dim!")
  }

  if(!is.null(mesh)){
    if(!inherits(mesh, c("inla.mesh", "inla.mesh.1d"))){
      stop("The mesh should be created using INLA!")
    }

    dim <- ifelse(inherits(mesh, "inla.mesh"), 2, 1)
  } 

  if(!is.null(graph.obj)){
    if(!inherits(graph.obj, "metric_graph")){
      stop("graph_obj should be a metric_graph object.")
    }
    dim <- 1
  }

  if (include.nu) {
    if (!is.null(nu.upper.bound)) {
      nu <- min(1, nu.upper.bound / 2)
    } else {
      nu <- 1
    }
  } else{
    if(is.null(nu)){
      stop("If include.nu is FALSE, then nu must be provided!")
    }
  }

  if(parameterization == "matern"){
    if(is.null(B.sigma)){
      B.sigma = matrix(c(0, 1, 0), 1, 3)
    }
    if(is.null(B.range)){
      B.range = matrix(c(0, 0, 1), 1, 3)
    }

    if(is.null(graph.obj)){
      param <- get_parameters_rSPDE(mesh = mesh,
                                  alpha = nu + dim/2,
                                  B.tau = B.tau,
                                  B.kappa = B.kappa,
                                  B.sigma = B.sigma,
                                  B.range = B.range,
                                  nu.nominal = nu,
                                  alpha.nominal = nu + dim/2,
                                  parameterization = parameterization,
                                  prior.std.dev.nominal = 1,
                                  prior.range.nominal = NULL,
                                  prior.tau = NULL,
                                  prior.kappa = NULL,
                                  theta.prior.mean = NULL,
                                  theta.prior.prec = 0.1,
                                  mesh.range = mesh.range,
                                  d = dim,
                                  n.spde = n.spde
                                  )
    } else{
      param <- get_parameters_rSPDE_graph(graph_obj = graph.obj,
                                  alpha = nu + 1/2,
                                  B.tau = B.tau,
                                  B.kappa = B.kappa,
                                  B.sigma = B.sigma,
                                  B.range = B.range,
                                  nu.nominal = nu,
                                  alpha.nominal = nu + 1/2,
                                  parameterization = parameterization,
                                  prior.std.dev.nominal = 1,
                                  prior.range.nominal = NULL,
                                  prior.tau = NULL,
                                  prior.kappa = NULL,
                                  theta.prior.mean = NULL,
                                  theta.prior.prec = 0.1
                                  )
    }

    initial <- param$theta.prior.mean
  } else{
    if(is.null(B.tau)){
      B.tau = matrix(c(0, 1, 0), 1, 3)
    }
    if(is.null(B.kappa)){
      B.kappa = matrix(c(0, 0, 1), 1, 3)
    }
    if(is.null(graph.obj)){
      param <- get_parameters_rSPDE(mesh = mesh,
                                  alpha = nu + dim/2,
                                  B.tau = B.tau,
                                  B.kappa = B.kappa,
                                  B.sigma = B.sigma,
                                  B.range = B.range,
                                  nu.nominal = nu,
                                  alpha.nominal = nu + dim/2,
                                  parameterization = parameterization,
                                  prior.std.dev.nominal = 1,
                                  prior.range.nominal = NULL,
                                  prior.tau = NULL,
                                  prior.kappa = NULL,
                                  theta.prior.mean = NULL,
                                  theta.prior.prec = 0.1,
                                  mesh.range = mesh.range,
                                  d = dim,
                                  n.spde = n.spde
                                  )
    } else {
      param <- get_parameters_rSPDE_graph(graph_obj = graph.obj,
                                  alpha = nu + 1/2,
                                  B.tau = B.tau,
                                  B.kappa = B.kappa,
                                  B.sigma = B.sigma,
                                  B.range = B.range,
                                  nu.nominal = nu,
                                  alpha.nominal = nu + 1/2,
                                  parameterization = parameterization,
                                  prior.std.dev.nominal = 1,
                                  prior.range.nominal = NULL,
                                  prior.tau = NULL,
                                  prior.kappa = NULL,
                                  theta.prior.mean = NULL,
                                  theta.prior.prec = 0.1
                                  )
    }

  initial <- param$theta.prior.mean
  }

  if(include.nu){
    initial <- c(initial, log(nu))
  }

  if (log.scale) {
    return(initial)
  } else {
    return(exp(initial))
  }
}




#' @name cut_decimals
#' @title Approximation function for covariance-based rSPDE models
#' @description Approximation function to be used to compute the
#' precision matrix for covariance-based rSPDE models
#' @param nu A real number
#' @return An approximation
#' @noRd

cut_decimals <- function(nu) {
  temp <- nu - floor(nu)
  if (temp < 10^(-3)) {
    temp <- 10^(-3)
  }
  if (temp > 0.999) {
    temp <- 0.999
  }
  return(temp)
}

#' @name check_class_inla_rspde
#' @title Check if the object inherits from inla_rspde class
#' @description Check if the object inherits from inla_rspde class
#' @param model A model to test if it inherits from inla_rspde
#' @return Gives an error if the object does not inherit from inla_rspde
#' @noRd

check_class_inla_rspde <- function(model) {
  if (!inherits(model, "inla_rspde")) {
    stop("You should provide a rSPDE model!")
  }
}

#' @name get_inla_mesh_dimension
#' @title Get the dimension of an INLA mesh
#' @description Get the dimension of an INLA mesh
#' @param inla_mesh An INLA mesh
#' @return The dimension of an INLA mesh.
#' @noRd


get_inla_mesh_dimension <- function(inla_mesh) {
  cond1 <- inherits(inla_mesh, "inla.mesh.1d")
  cond2 <- inherits(inla_mesh, "inla.mesh")
  stopifnot(cond1 || cond2)
  if (inla_mesh$manifold == "R1") {
    d <- 1
  } else if (inla_mesh$manifold %in% c("R2","S2")) {
    d <- 2
  } else {
    stop("The mesh should be from a flat manifold.")
  }
  return(d)
}

#' @name fem_mesh_order_1d
#' @title Get fem_mesh_matrices for 1d inla.mesh objects
#' @description Get fem_mesh_matrices for 1d inla.mesh objects
#' @param inla_mesh An INLA mesh
#' @param m_order the order of the FEM matrices
#' @return A list with fem_mesh_matrices
#' @noRd


fem_mesh_order_1d <- function(inla_mesh, m_order) {
  # fem_mesh <- rSPDE.fem1d(inla_mesh[["loc"]])
  # mesh_1d <- fmesher::fm_mesh_1d(inla_mesh[["loc"]])
  # fem_mesh <- fmesher::fm_fem(mesh_1d)
  mesh_1d <- fm_mesh_1d(inla_mesh[["loc"]])
  fem_mesh <- fm_fem(mesh_1d)
  C <- fem_mesh$c0
  C <- Matrix::Diagonal(dim(C)[1], rowSums(C))
  C <- as(C,"TsparseMatrix")
  G <- fem_mesh$g1
  Gk <- list()
  Ci <- C
  Ci@x <- 1 / (C@x)

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

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

  for (i in 1:m_order) {
    fem_mesh_matrices[[paste0("g", i)]] <- Gk[[i]]
  }
  return(fem_mesh_matrices)
}

#' @name generic_fem_mesh_order
#' @title Get fem_mesh_matrices from C and G matrices
#' @description Get fem_mesh_matrices from C and G matrices
#' @param fem_matrices A list with objects C and G
#' @param m_order the order of the FEM matrices
#' @return A list with fem_mesh_matrices
#' @noRd


generic_fem_mesh_order <- function(fem_matrices, m_order) {
  C <- fem_matrices$C
  C <- Matrix::Diagonal(dim(C)[1], rowSums(C))
  C <- INLA::inla.as.sparse(C)
  # C <- as(C,"TsparseMatrix")
  G <- fem_matrices$G
  Gk <- list()
  Ci <- C
  Ci@x <- 1 / (C@x)

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

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

  for (i in 1:m_order) {
    fem_mesh_matrices[[paste0("g", i)]] <- Gk[[i]]
  }
  return(fem_mesh_matrices)
}


#' @name get.sparsity.graph.rspde
#' @title Sparsity graph for rSPDE models
#' @description Creates the sparsity graph for rSPDE models
#' @param mesh An INLA mesh, optional
#' @param fem_mesh_matrices A list containing the FEM-related matrices.
#' The list should contain elements C, G, G_2, G_3, etc. Optional,
#' should be provided if mesh is not provided.
#' @param dim The dimension, optional. Should be provided if mesh
#' is not provided.
#' @param nu The smoothness parameter
#' @param force_non_integer Should nu be treated as non_integer?
#' @param rspde.order The order of the covariance-based rational SPDE approach.
#' @return The sparsity graph for rSPDE models to be used in R-INLA interface.
#' @noRd

get.sparsity.graph.rspde <- function(mesh = NULL,
                                     fem_mesh_matrices = NULL,
                                     nu,
                                     force_non_integer = FALSE,
                                     rspde.order = 2,
                                     dim = NULL) {
  if (!is.null(mesh)) {
    stopifnot(inherits(mesh, "inla.mesh"))
    if (mesh$manifold == "R1") {
      dim <- 1
    } else if (mesh$manifold == "R2") {
      dim <- 2
    } else {
      stop("The mesh should be from a flat manifold.")
    }
  } else if (is.null(dim)) {
    stop("If an INLA mesh is not provided, you should provide the dimension!")
  }
  sharp = TRUE
  alpha <- nu + dim / 2

  m_alpha <- floor(alpha)

  integer_alpha <- (alpha %% 1 == 0)

  if (force_non_integer) {
    integer_alpha <- FALSE
  }

  if (!is.null(fem_mesh_matrices)) {
    if (integer_alpha) {
      return(fem_mesh_matrices[[paste0("g", m_alpha)]])
    } else {
      if (sharp) {
        if (m_alpha > 0) {
          return(bdiag(
            kronecker(
              diag(rep(1, rspde.order)),
              fem_mesh_matrices[[paste0("g", m_alpha + 1)]]
            ),
            fem_mesh_matrices[[paste0("g", m_alpha)]]
          ))
        } else {
          return(bdiag(
            kronecker(
              diag(rep(1, rspde.order)),
              fem_mesh_matrices[["g1"]]
            ),
            fem_mesh_matrices[["c0"]]
          ))
        }
      } else {
        return(kronecker(
          diag(rep(1, rspde.order + 1)),
          fem_mesh_matrices[[paste0("g", m_alpha + 1)]]
        ))
      }
    }
  } else if (!is.null(mesh)) {
    if (integer_alpha) {
      # fem_mesh_matrices <- INLA::inla.mesh.fem(mesh, order = m_alpha)
      # fem_mesh_matrices <- fmesher::fm_fem(mesh, order = m_alpha)
      fem_mesh_matrices <- fm_fem(mesh, order = m_alpha)
      return(fem_mesh_matrices[[paste0("g", m_alpha)]])
    } else {
      if (dim == 2) {
        # fem_mesh_matrices <- INLA::inla.mesh.fem(mesh, order = m_alpha + 1)
        # fem_mesh_matrices <- fmesher::fm_fem(mesh, order = m_alpha + 1)
        fem_mesh_matrices <- fm_fem(mesh, order = m_alpha + 1)
      } else {
        fem_mesh_matrices <- fem_mesh_order_1d(mesh, m_order = m_alpha + 1)
      }


      if (sharp) {
        if (m_alpha > 0) {
          return(bdiag(
            kronecker(
              diag(rep(1, rspde.order)),
              fem_mesh_matrices[[paste0("g", m_alpha + 1)]]
            ),
            fem_mesh_matrices[[paste0("g", m_alpha)]]
          ))
        } else {
          return(bdiag(
            kronecker(
              diag(rep(1, rspde.order)),
              fem_mesh_matrices[["g1"]]
            ),
            fem_mesh_matrices[["c0"]]
          ))
        }
      } else {
        return(kronecker(
          diag(rep(1, rspde.order + 1)),
          fem_mesh_matrices[[paste0("g", m_alpha + 1)]]
        ))
      }
    }
  } else {
    stop("You should provide either mesh or fem_mesh_matrices!")
  }
}


#' @name build_sparse_matrix_rspde
#' @title Create sparse matrix from entries and graph
#' @description Create sparse matrix from entries and graph
#' @param entries The entries of the precision matrix
#' @param graph The sparsity graph of the precision matrix
#' @return index for rSPDE models.
#' @noRd

build_sparse_matrix_rspde <- function(entries, graph) {
  if (!is.null(graph)) {
    # graph <- as(graph, "dgTMatrix")
    graph <- as(graph,"TsparseMatrix")
    idx <- which(graph@i <= graph@j)
    Q <- Matrix::sparseMatrix(
      i = graph@i[idx], j = graph@j[idx], x = entries,
      symmetric = TRUE, index1 = FALSE
    )
  }
  return(Q)
}


#' @name analyze_sparsity_rspde
#' @title Analyze sparsity of matrices in the rSPDE approach
#' @description Auxiliar function to analyze sparsity of matrices
#' in the rSPDE approach
#' @param nu.upper.bound Upper bound for the smoothness parameter
#' @param dim The dimension of the domain
#' @param rspde.order The order of the rational approximation
#' @param fem_mesh_matrices A list containing FEM-related matrices.
#' The list should contain elements c0, g1, g2, g3, etc.
#' @param include_lower_order Logical. Should the lower-order terms
#' be included? They are needed for the cases
#' when alpha = nu + d/2 is integer or for when sharp is set to TRUE.
#' @param include_higher_order Logical. Should be included for when nu
#' is estimated or for when alpha = nu + d/2 is not an integer.
#' @return A list containing informations on sparsity of the precision matrices
#' @noRd

analyze_sparsity_rspde <- function(nu.upper.bound, dim, rspde.order,
                                   fem_mesh_matrices,
                                   include_lower_order = TRUE,
                                   include_higher_order = TRUE) {
  beta <- nu.upper.bound / 2 + dim / 4

  m_alpha <- floor(2 * beta)

  positions_matrices <- list()

  C_list <- symmetric_part_matrix(fem_mesh_matrices$c0)
  G_1_list <- symmetric_part_matrix(fem_mesh_matrices$g1)
  if (m_alpha < 2) {
    G_2_list <- symmetric_part_matrix(fem_mesh_matrices[["g2"]])
  }
  if (m_alpha > 1) {
    for (j in 2:(m_alpha)) {
      assign(paste0("G_", j, "_list"),
      symmetric_part_matrix(fem_mesh_matrices[[paste0("g", j)]]))
    }
  }

    if (include_higher_order) {
      assign(paste0("G_", m_alpha + 1, "_list"),
      symmetric_part_matrix(fem_mesh_matrices[[paste0("g",
      m_alpha + 1)]]))

      positions_matrices[[1]] <- match(C_list$M,
      get(paste0("G_", m_alpha + 1, "_list"))[["M"]])
    }

  idx_matrices <- list()

  idx_matrices[[1]] <- C_list$idx

  if (m_alpha > 0) {
    for (i in 1:m_alpha) {
      if (include_higher_order) {
        positions_matrices[[i + 1]] <- match(get(paste0("G_", i,
        "_list"))[["M"]], get(paste0("G_", m_alpha + 1,
        "_list"))[["M"]])
      }
      idx_matrices[[i + 1]] <- get(paste0("G_", i, "_list"))[["idx"]]
    }
  }

  if (include_higher_order) {
    idx_matrices[[m_alpha + 2]] <- get(paste0("G_", m_alpha + 1,
    "_list"))[["idx"]]
  }

  if (include_lower_order) {
    positions_matrices_less <- list()
    if (m_alpha > 0) {
      positions_matrices_less[[1]] <- match(C_list$M, get(paste0("G_",
      m_alpha, "_list"))[["M"]])
    } else {
      positions_matrices_less[[1]] <- match(C_list$M, get(paste0("G_",
      1, "_list"))[["M"]])
    }

    if (m_alpha > 1) {
      for (i in 1:(m_alpha - 1)) {
        positions_matrices_less[[i + 1]] <- match(get(paste0("G_", i,
        "_list"))[["M"]], get(paste0("G_", m_alpha, "_list"))[["M"]])
      }
    } else if (m_alpha == 1) {
      positions_matrices_less[[2]] <- seq_len(length(get(paste0("G_",
      m_alpha, "_list"))[["M"]]))
    }
  } else {
    positions_matrices_less <- NULL
  }

  return(list(
    positions_matrices = positions_matrices,
    idx_matrices = idx_matrices,
    positions_matrices_less = positions_matrices_less
  ))
}

#' @name symmetric_part_matrix
#' @title Gets the upper triangular part of a matrix
#' @description Gets the upper triangular part of a matrix
#' @param M A matrix or a sparse matrix
#' @return A sparse matrix formed by the upper triangular part of `M`.
#' @noRd

symmetric_part_matrix <- function(M) {
  # M <- as(M, "dgTMatrix")
  M <- as(M,"TsparseMatrix")
  idx <- which(M@i <= M@j)
  sM <- cbind(M@i[idx], M@j[idx])
  colnames(sM) <- NULL
  return(list(M = split(sM, seq(nrow(sM))), idx = idx))
}


#' @name create_summary_from_density
#' @title Creates a summary from a density data frame
#' @description Auxiliar function to create summaries from density data drames
#' @param density_df A density data frame
#' @param name Name of the parameter
#' @return A data frame containing a basic summary
#' @noRd

create_summary_from_density <- function(density_df, name) {
  min_x <- min(density_df[, "x"])
  max_x <- max(density_df[, "x"])
  denstemp <- function(x) {
    dens <- sapply(x, function(z) {
      if (z < min_x) {
        return(0)
      } else if (z > max_x) {
        return(0)
      } else {
        return(approx(x = density_df[, "x"], y = density_df[, "y"], xout = z)$y)
      }
    })
    return(dens)
  }

  ptemp <- function(q) {
    prob_temp <- sapply(q, function(v) {
      if (v <= min_x) {
        return(0)
      } else if (v >= max_x) {
        return(1)
      } else {
        stats::integrate(
          f = denstemp, lower = min_x, upper = v,
                  stop.on.error = FALSE
        )$value
      }
    })
    return(prob_temp)
  }

  mean_temp <- stats::integrate(
    f = function(z) {
      denstemp(z) * z
    }, lower = min_x, upper = max_x,
    subdivisions = nrow(density_df),
                  stop.on.error = FALSE
  )$value

  sd_temp <- sqrt(stats::integrate(
    f = function(z) {
      denstemp(z) * (z - mean_temp)^2
    }, lower = min_x, upper = max_x,
                  stop.on.error = FALSE
  )$value)

  mode_temp <- density_df[which.max(density_df[, "y"]), "x"]

  qtemp <- function(p) {
    quant_temp <- sapply(p, function(x) {
      if (x < 0 | x > 1) {
        return(NaN)
      } else {
        return(stats::uniroot(function(y) {
          ptemp(y) - x
        }, lower = min_x, upper = max_x)$root)
      }
    })
    return(quant_temp)
  }

  out <- data.frame(
    mean = mean_temp, sd = sd_temp, `0.025quant` = qtemp(0.025),
    `0.5quant` = qtemp(0.5), `0.975quant` = qtemp(0.975), mode = mode_temp
  )
  rownames(out) <- name
  colnames(out) <- c("mean", "sd", "0.025quant",
  "0.5quant", "0.975quant", "mode")
  return(out)
}






#' Summarise CBrSPDE objects
#'
#' Summary method for class "CBrSPDEobj"
#'
#' @param object an object of class "CBrSPDEobj", usually, a result of a call
#'   to [matern.operators()].
#' @param ... further arguments passed to or from other methods.
#' @export
#' @method summary CBrSPDEobj
#' @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
#' 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"
#' )
#'
#' op_cov
summary.CBrSPDEobj <- function(object, ...) {
  out <- list()
  class(out) <- "summary.CBrSPDEobj"
  out$type <- object$type
  out$kappa <- object$kappa
  out$sigma <- object$sigma
  out$theta <- object$theta
  out$tau <- object$tau
  out[["range"]] <- object[["range"]]  
  out$nu <- object$nu
  out$m <- object$m
  out$stationary <- object$stationary
  out$parameterization <- object$parameterization  
  out$n <- dim(object$C)[1]
  out[["type_rational_approximation"]] <-
  object[["type_rational_approximation"]]
  return(out)
}

#' @param x an object of class "summary.CBrSPDEobj", usually, a result of a call
#'   to [summary.CBrSPDEobj()].
#' @export
#' @method print summary.CBrSPDEobj
#' @rdname summary.CBrSPDEobj
print.summary.CBrSPDEobj <- function(x, ...) {
  cat("Type of approximation: ", x$type, "\n")
  cat("Parameterization: ", x$parameterization, "\n")
  cat("Type of rational approximation: ",
  x[["type_rational_approximation"]], "\n")
  if(x$stationary){
    if(x$parameterization == "spde"){
      cat(
      "Parameters of covariance function: kappa = ",
      x$kappa, ", tau = ", x$tau, ", nu = ", x$nu, "\n"
    )} else{
        cat(
      "Parameters of covariance function: range = ",
      x[["range"]], ", sigma = ", x$sigma, ", nu = ", x$nu, "\n"
    )
    }
  } else if (!is.null(x$theta)){
        cat(
      "Parameters of covariance function: theta = ",
      x$theta, ", nu = ", x$nu, "\n"
    )
  }

  cat("Order or rational approximation: ", x$m, "\n")
  cat("Size of discrete operators: ", x$n, " x ", x$n, "\n")
  if(x$stationary){
    cat("Stationary Model\n")
  } else{
    cat("Non-Stationary Model")
  }
}

#' @export
#' @method print CBrSPDEobj
#' @rdname summary.CBrSPDEobj
print.CBrSPDEobj <- function(x, ...) {
  print.summary.CBrSPDEobj(summary(x))
}



#' @name get_rational_coefficients
#' @title Get matrix with rational coefficients
#' @description Get matrix with rational coefficients
#' @param order order of the rational approximation
#' @param type_rational_approx Type of the rational
#' approximation. Options are "chebfun", "brasil"
#' and "chebfunLB"
#' @return A matrix with rational approximations.
#' @noRd

get_rational_coefficients <- function(order, type_rational_approx) {
  if (type_rational_approx == "chebfun") {
    mt <- get(paste0("m", order, "t"))
  } else if (type_rational_approx == "brasil") {
    mt <- get(paste0("m_brasil", order, "t"))
  } else if (type_rational_approx == "chebfunLB") {
    mt <- get(paste0("m_chebfun", order, "t"))
  } else {
    stop("The options are 'chebfun', 'brasil' and 'chebfunLB'!")
  }
  return(mt)
}




#' Changing the type of the rational approximation
#'
#' @param x A `CBrSPDE` or an `rpsde.inla` object
#' @param value The type of rational approximation.
#' The current options are "chebfun", "brasil" and "chebfunLB"
#'
#' @return An object of the same class with the new rational approximation.
#' @export
#'
`rational.type<-` <- function(x, value) {
  object <- x

  type_rational_approximation <- value
  type_rational_approximation <- type_rational_approximation[[1]]
  if (!(type_rational_approximation %in% c("chebfun", "brasil", "chebfunLB"))) {
    stop('The possible types are "chebfun", "brasil" and "chebfunLB"!')
  }
  if (inherits(object, "CBrSPDEobj")) {
    model <- update(x, type_rational_approximation = value)
  } else if (inherits(object, "inla_rspde")) {
    nu.upper.bound <- object$nu.upper.bound
    prior.nu.dist <- object$prior.nu.dist
    mesh <- object$mesh
    nu <- object[["nu"]]
    rspde.order <- object$rspde.order
    parameterization <- object$parameterization
    theta.prior.prec <- object$theta.prior.prec
    theta.prior.mean <- object$theta.prior.mean
    start.theta <- object$start.theta
    prior.nu <- object$prior.nu
    start.nu <- object$start.nu
    debug <- object$debug


    model <- rspde.matern(mesh,
      nu.upper.bound = nu.upper.bound,
      rspde.order = rspde.order,
      nu = nu,
      debug = debug,
      parameterization = parameterization,
      theta.prior.mean = theta.prior.mean,
      theta.prior.prec = theta.prior.prec,
      start.theta = start.theta,
      prior.nu = prior.nu,
      start.nu = start.nu,
      prior.nu.dist = prior.nu.dist,
      type.rational.approx = type_rational_approximation
    )
    
  } else {
    stop("The object must be of class 'CBrSPDE' or 'inla_rspde'!")
  }
  return(model)
}



#' Get type of rational approximation.
#'
#' @param object A `CBrSPDEobj` object or an `inla_rspde` object.
#'
#' @return The type of rational approximation.
#' @export
#'
rational.type <- function(object) {
  if (inherits(object, "CBrSPDEobj")) {
    return(object$type_rational_approximation)
  } else if (inherits(object, "inla_rspde")) {
    return(object$type.rational.approx)
  } else if (inherits(object, "rSPDEobj")) {
    return("chebfun")
  } else {
    stop("Not a valid rSPDE object!")
  }
}


#' Changing the order of the rational approximation
#'
#' @param x A `CBrSPDE` or an `rpsde.inla` object
#' @param value The order of rational approximation.
#'
#' @return An object of the same class with the new order
#' of rational approximation.
#' @export
#'
`rational.order<-` <- function(x, value) {
  object <- x

  rspde.order <- value
  rspde.order <- rspde.order[[1]]

  if (inherits(object, "CBrSPDEobj") || inherits(object, "rSPDEobj")) {
    model <- update(object, user_m = rspde.order)
  } else if (inherits(object, "inla_rspde")) {
    if (rspde.order > 0 && object$integer.nu) {
      warning("The order was not changed since there is no
      rational approximation (an integer model was
      considered).")
      return(object)
    }
    nu.upper.bound <- object$nu.upper.bound
    prior.nu.dist <- object$prior.nu.dist
    mesh <- object$mesh
    nu <- object[["nu"]]
    parameterization <- object$parameterization
    theta.prior.prec <- object$theta.prior.prec
    theta.prior.mean <- object$theta.prior.mean
    start.theta <- object$start.theta
    prior.nu <- object$prior.nu
    start.nu <- object$start.nu
    type_rational_approximation <- object$type.rational.approx
    debug <- object$debug


    model <- rspde.matern(mesh,
      nu.upper.bound = nu.upper.bound,
      rspde.order = rspde.order,
      nu = nu,
      debug = debug,
      parameterization = parameterization,
      theta.prior.mean = theta.prior.mean,
      theta.prior.prec = theta.prior.prec,
      start.theta = start.theta,
      prior.nu = prior.nu,
      start.nu = start.nu,
      prior.nu.dist = prior.nu.dist,
      type.rational.approx = type_rational_approximation
    )
  } else if (!is.null(attr(object, "inla_rspde_Amatrix"))) {
    n_temp <- ncol(object)
    old_rspde.order <- attr(object, "rspde.order")
    orig_dim <- n_temp / (old_rspde.order + 1)
    A <- object[, 1:orig_dim]
    Abar <- kronecker(matrix(1, 1, rspde.order + 1), A)
    attr(Abar, "inla_rspde_Amatrix") <- TRUE
    attr(Abar, "rspde.order") <- rspde.order
    integer_nu <- attr(object, "integer_nu")
    if (integer_nu && rspde.order > 0) {
      warning("The order was not changed since there is
      no rational approximation (an integer model was
      considered).")
      return(object)
    }
    attr(Abar, "integer_nu") <- integer_nu
    return(Abar)
  } else if (inherits(object, "inla_rspde_index")) {
    integer_nu <- attr(object, "integer_nu")

    if (integer_nu && rspde.order > 0) {
      warning("The order was not changed since there is
      no rational approximation (an integer model was
      considered).")
      return(object)
    }

    n_mesh <- attr(object, "n.mesh")
    name <- attr(object, "name")
    n.group <- attr(object, "n.group")
    n.repl <- attr(object, "n.repl")

    factor_rspde <- rspde.order + 1

    name.group <- paste(name, ".group", sep = "")
    name.repl <- paste(name, ".repl", sep = "")

    out <- list()
    out[[name]] <- as.vector(sapply(1:factor_rspde, function(i) {
      rep(rep(((i - 1) * n_mesh + 1):(i * n_mesh),
      times = n.group), times = n.repl)
    }))
    out[[name.group]] <- rep(rep(rep(1:n.group, each = n_mesh),
    times = n.repl), times = factor_rspde)
    out[[name.repl]] <- rep(rep(1:n.repl, each = n_mesh * n.group),
    times = factor_rspde)
    class(out) <- c("inla_rspde_index", class(out))
    attr(out, "rspde.order") <- rspde.order
    attr(out, "integer_nu") <- integer_nu
    attr(out, "n.mesh") <- n_mesh
    attr(out, "name") <- name
    attr(out, "n.group") <- n.group
    attr(out, "n.repl") <- n.repl
    return(out)
  } else {
    stop("The object must be of class 'CBrSPDE' or 'inla_rspde'!")
  }
  return(model)
}


#' Get the order of rational approximation.
#'
#' @param object A `CBrSPDEobj` object or an `inla_rspde` object.
#'
#' @return The order of rational approximation.
#' @export
#'
rational.order <- function(object) {
  if (inherits(object, "CBrSPDEobj") || inherits(object, "rSPDEobj")) {
    return(object$m)
  } else if (inherits(object, "inla_rspde")) {
    return(object$rspde.order)
  } else if (!is.null(attr(object, "inla_rspde_Amatrix"))) {
    return(attr(object, "rspde.order"))
  } else if (inherits(object, "inla_rspde_index")) {
    return(attr(object, "rspde.order"))
  } else {
    stop("Not a valid rSPDE object!")
  }
}


#' Check user input.
#'
#' @param parameter A parameter.
#' @param label Label for the parameter
#' @param check_null Check if parameter is null.
#'
#' @return Check the parameter.
#' @noRd
#'

rspde_check_user_input <- function(param, label, lower_bound = NULL){
  if(is.null(lower_bound)){
    if (!is.numeric(param)) {
      stop(paste(param,"should be a number!"))
    }
    if (length(param) > 1) {
      stop(paste(param,"should be a number!"))
    }
    return(param)
  } else{
    if (!is.numeric(param)) {
      stop(paste(param,"should be a number greater or equal to",lower_bound))
    }
    if (length(param) > 1) {
      stop(paste(param,"should be a number greater or equal to",lower_bound))
    }
    if (param < lower_bound) {
      stop(paste(param,"should be a number greater or equal to",lower_bound))
    }
    return(param)
  }
  }


  #' Process inputs likelihood
  #'
  #' @param user_kappa kappa
  #' @param user_tau tau
  #' @param user_nu nu
  #' @param sigma.e sigma.e
  #'
  #' @return List with the positions
  #' @noRd

likelihood_process_inputs_spde <- function(user_kappa, user_tau, user_nu, sigma.e){
  param_vector <- c("tau", "kappa", "nu", "sigma.e")
  if(!is.null(user_tau)){
    param_vector <- setdiff(param_vector, "tau")
  } 
  if(!is.null(user_kappa)){
    param_vector <- setdiff(param_vector, "kappa")
  }
  if(!is.null(user_nu)){
    param_vector <- setdiff(param_vector, "nu")
  }
  if(!is.null(sigma.e)){
    param_vector <- setdiff(param_vector, "sigma.e")
  }
  if(length(param_vector)==0){
    stop("You should leave at least one parameter free.")
  }
  return(param_vector)
}

  #' Process inputs likelihood
  #'
  #' @param user_kappa kappa
  #' @param user_tau tau
  #' @param user_nu nu
  #' @param sigma.e sigma.e
  #'
  #' @return List with the positions
  #' @noRd

likelihood_process_inputs_matern <- function(user_range, user_sigma, user_nu, sigma.e){
  param_vector <- c("sigma", "range", "nu", "sigma.e")
  if(!is.null(user_sigma)){
    param_vector <- setdiff(param_vector, "sigma")
  } 
  if(!is.null(user_range)){
    param_vector <- setdiff(param_vector, "range")
  }
  if(!is.null(user_nu)){
    param_vector <- setdiff(param_vector, "nu")
  }
  if(!is.null(sigma.e)){
    param_vector <- setdiff(param_vector, "sigma.e")
  }
  if(length(param_vector)==0){
    stop("You should leave at least one parameter free.")
  }
  return(param_vector)
}

#' Process parameters likelihood
#'
#' @param theta vector of parameters
#' @param param_vector vector of parameters to be used
#' @param which_par which parameter to consider
#' @param logscale log scale?
#'
#' @return The value in the correct scale
#' @noRd

likelihood_process_parameters <- function(theta, param_vector, which_par, logscale){
  coord_par <- which(which_par == param_vector)
  if(logscale){
    param_value <- exp(theta[[coord_par]])
  } else{
    param_value <- theta[[coord_par]]
  }
  return(param_value)
}


#' @noRd 
# Get priors and starting values
# Based on INLA::param2.matern.orig()

get_parameters_rSPDE <- function (mesh, alpha, 
    B.tau, 
    B.kappa, 
    B.sigma,
    B.range, 
    nu.nominal,
    alpha.nominal,
    parameterization,
    prior.std.dev.nominal, 
    prior.range.nominal, 
    prior.tau, 
    prior.kappa, 
    theta.prior.mean, 
    theta.prior.prec,
    mesh.range = NULL,
    d = NULL,
    n.spde = NULL) 
{
  if(!is.null(mesh)){
    if(!inherits(mesh, c("inla.mesh", "inla.mesh.1d"))){
      stop("The mesh should be created using INLA!")
    }

    d <- ifelse(inherits(mesh, "inla.mesh"), 2, 1)
    n.spde <- ifelse(d == 2, mesh$n, mesh$m)
  } else{
    if(is.null(d)){
      stop("If you do not provide the mesh, you must provide the dimension!")
    }
    if(is.null(n.spde)){
      stop("If you do not provide the mesh, you must provide n.spde!")
    }
  } 

    if (is.null(B.tau) && is.null(B.sigma)) 
        stop("One of B.tau or B.sigma must not be NULL.")
    if (is.null(B.kappa) && is.null(B.range)) 
        stop("One of B.kappa or B.range must not be NULL.")



    if(parameterization == "spde"){
      n.theta <- ncol(B.kappa) - 1L

      B.kappa <- prepare_B_matrices(B.kappa, n.spde, 
          n.theta)
      B.tau <- prepare_B_matrices(B.tau, n.spde, n.theta)
    } else if(parameterization == "matern"){
      n.theta <- ncol(B.sigma) - 1L
      
      B.sigma <- prepare_B_matrices(B.sigma, n.spde, 
          n.theta)
      B.range <- prepare_B_matrices(B.range, n.spde, 
          n.theta)

      B.kappa <- cbind(0.5 * log(8 * nu.nominal) - B.range[, 1], 
        -B.range[, -1, drop = FALSE])

      B.tau <- cbind(0.5 * (lgamma(nu.nominal) - lgamma(alpha.nominal) - 
                d/2 * log(4 * pi)) - nu.nominal * B.kappa[, 1] - 
                B.sigma[,1], 
                - nu.nominal * B.kappa[, -1, drop = FALSE] -
                B.sigma[, -1, drop = FALSE])
  } else if(parameterization == "matern2"){
      n.theta <- ncol(B.sigma) - 1L
      
      B.sigma <- prepare_B_matrices(B.sigma, n.spde, 
          n.theta)
      B.range <- prepare_B_matrices(B.range, n.spde, 
          n.theta)

      B.kappa <- -B.range

      B.tau <- cbind(0.5 * (lgamma(nu.nominal) - lgamma(alpha.nominal) - 
                d/2 * log(4 * pi)) - nu.nominal * B.kappa[, 1] - 
                0.5 * B.sigma[,1], 
                - nu.nominal * B.kappa[, -1, drop = FALSE] -
                0.5 * B.sigma[, -1, drop = FALSE])
    }


    if (is.null(theta.prior.prec)) {
        theta.prior.prec = diag(0.1, n.theta, n.theta)
    }
    else {
        theta.prior.prec = as.matrix(theta.prior.prec)
        if (ncol(theta.prior.prec) == 1) {
            theta.prior.prec = diag(as.vector(theta.prior.prec), 
                n.theta, n.theta)
        }
        if ((nrow(theta.prior.prec) != n.theta) || (ncol(theta.prior.prec) != 
            n.theta)) {
            stop(paste("Size of theta.prior.prec is (", paste(dim(theta.prior.prec), 
                collapse = ",", sep = ""), ") but should be (", 
                paste(c(n.theta, n.theta), collapse = ",", sep = ""), 
                ")."))
        }
    }

    
    if (is.null(theta.prior.mean)) {
        if (is.null(prior.range.nominal)) {
          if(is.null(mesh.range)){
            mesh.range = ifelse(d == 2, (max(c(diff(range(mesh$loc[, 
                1])), diff(range(mesh$loc[, 2])), diff(range(mesh$loc[, 
                3]))))), diff(mesh$interval))
          }
            prior.range.nominal = mesh.range * 0.2
        }
        if (is.null(prior.kappa)) {
            prior.kappa = sqrt(8 * nu.nominal)/prior.range.nominal
        }
        if (is.null(prior.tau)) {
            prior.tau = sqrt(gamma(nu.nominal)/gamma(alpha.nominal)/((4 * 
                pi)^(d/2) * prior.kappa^(2 * nu.nominal) * prior.std.dev.nominal^2))
        }
        if (n.theta > 0) {
          if(parameterization == "spde"){
              theta.prior.mean = qr.solve(rbind(B.tau[, -1, drop = FALSE], 
                  B.kappa[, -1, drop = FALSE]), c(log(prior.tau) - 
                  B.tau[, 1], log(prior.kappa) - B.kappa[, 1]))
          } else if(parameterization == "matern"){
              theta.prior.mean = qr.solve(rbind(B.sigma[, -1, drop = FALSE], 
                  B.range[, -1, drop = FALSE]), c(log(prior.std.dev.nominal) - 
                  B.sigma[, 1], log(prior.range.nominal) - B.range[, 1]))
          } else if(parameterization == "matern2"){
              theta.prior.mean = qr.solve(rbind(B.sigma[, -1, drop = FALSE], 
                  B.range[, -1, drop = FALSE]), c(2*log(prior.std.dev.nominal) - 
                  B.sigma[, 1], -log(prior.kappa) - B.range[, 1]))
          }
        }
        else {
            theta.prior.mean = rep(0, n.theta)
        }
    }
    param = list(B.tau = B.tau, 
        B.kappa = B.kappa, theta.prior.mean = theta.prior.mean, 
        theta.prior.prec = theta.prior.prec)
    return(param)
}

#' @noRd 
# Check B matrices and adjust the number of lines
# Based on INLA:::inla.spde.homogenise_B_matrix()

prepare_B_matrices <- function (B, n.spde, n.theta) 
{
    if (!is.numeric(B)) {
        stop("B matrix must be numeric.")
    }
    if (is.matrix(B)) {
        if ((nrow(B) != 1) && (nrow(B) != n.spde)) {
            stop(paste("B matrix must have either 1 or", as.character(n.spde),"rows."))
        }
        if ((ncol(B) != 1) && (ncol(B) != 1 + n.theta)) {
            stop(paste("B matrix must have 1 or", as.character(1 + 
                  n.theta),"columns."))
        }
        if (ncol(B) == 1) {
            return(cbind(as.vector(B), matrix(0, n.spde, n.theta)))
        }
        else if (ncol(B) == 1 + n.theta) {
            if (nrow(B) == 1) {
                return(matrix(as.vector(B), n.spde, 1 + n.theta, 
                  byrow = TRUE))
            }
            else if (nrow(B) == n.spde) {
                return(B)
            }
        }
    }
    else {
        if ((length(B) == 1) || (length(B) == n.spde)) {
            return(cbind(B, matrix(0, n.spde, n.theta)))
        }
        else if (length(B) == 1 + n.theta) {
            return(matrix(B, n.spde, 1 + n.theta, byrow = TRUE))
        }
        else {
            stop(paste("Length of B must be 1,", as.character(1 + n.theta), 
                "or", as.character(n.spde)))
        }
    }
    stop("Unrecognised structure for B matrix")
}



#' @noRd 
# Get priors and starting values
# Based on INLA::param2.matern.orig()

get_parameters_rSPDE_graph <- function (graph_obj, alpha, 
    B.tau, 
    B.kappa, 
    B.sigma,
    B.range, 
    nu.nominal,
    alpha.nominal,
    parameterization,
    prior.std.dev.nominal, 
    prior.range.nominal, 
    prior.tau, 
    prior.kappa, 
    theta.prior.mean, 
    theta.prior.prec) 
{
    if(!inherits(graph_obj, "metric_graph")){
      stop("The graph object should be of class metric_graph!")
    }
    if (is.null(B.tau) && is.null(B.sigma)) 
        stop("One of B.tau or B.sigma must not be NULL.")
    if (is.null(B.kappa) && is.null(B.range)) 
        stop("One of B.kappa or B.range must not be NULL.")

    d <- 1
    n.spde <- nrow(graph_obj$mesh$C)

    if(parameterization == "spde"){
      n.theta <- ncol(B.kappa) - 1L
      B.kappa <- prepare_B_matrices(B.kappa, n.spde, 
          n.theta)
      B.tau <- prepare_B_matrices(B.tau, n.spde, n.theta)
    } else if(parameterization == "matern"){
      n.theta <- ncol(B.sigma) - 1L
      B.sigma <- prepare_B_matrices(B.sigma, n.spde, 
          n.theta)
      B.range <- prepare_B_matrices(B.range, n.spde, 
          n.theta)

      B.kappa <- cbind(0.5 * log(8 * nu.nominal) - B.range[, 1], 
        -B.range[, -1, drop = FALSE])

      B.tau <- cbind(0.5 * (lgamma(nu.nominal) - lgamma(alpha.nominal) - 
                d/2 * log(4 * pi)) - nu.nominal * B.kappa[, 1] - 
                B.sigma[,1], - nu.nominal * B.kappa[, -1, drop = FALSE] -
                B.sigma[, -1, drop = FALSE])
    }


    if (is.null(theta.prior.prec)) {
        theta.prior.prec = diag(0.1, n.theta, n.theta)
    }
    else {
        theta.prior.prec = as.matrix(theta.prior.prec)
        if (ncol(theta.prior.prec) == 1) {
            theta.prior.prec = diag(as.vector(theta.prior.prec), 
                n.theta, n.theta)
        }
        if ((nrow(theta.prior.prec) != n.theta) || (ncol(theta.prior.prec) != 
            n.theta)) {
            stop(paste("Size of theta.prior.prec is (", paste(dim(theta.prior.prec), 
                collapse = ",", sep = ""), ") but should be (", 
                paste(c(n.theta, n.theta), collapse = ",", sep = ""), 
                ")."))
        }
    }

    
    if (is.null(theta.prior.mean)) {
        if (is.null(prior.range.nominal)) {
            if(is.null(graph_obj$geo_dist)){
              graph_obj$compute_geodist(obs=FALSE)
            } else if(is.null(graph_obj$geo_dist[[".vertices"]])){
              graph_obj$compute_geodist(obs=FALSE)
            }
            finite_geodist <- is.finite(graph_obj$geo_dist[[".vertices"]])
            finite_geodist <- graph_obj$geo_dist[[".vertices"]][finite_geodist]
            prior.range.nominal <- max(finite_geodist) * 0.2
        }
        if (is.null(prior.kappa)) {
            prior.kappa = sqrt(8 * nu.nominal)/prior.range.nominal
        }
        if (is.null(prior.tau)) {
            prior.tau = sqrt(gamma(nu.nominal)/gamma(alpha.nominal)/(4 * 
                pi * prior.kappa^(2 * nu.nominal) * prior.std.dev.nominal^2))
        }
        if (n.theta > 0) {
          if(parameterization == "spde"){
              theta.prior.mean = qr.solve(rbind(B.tau[, -1, drop = FALSE], 
                  B.kappa[, -1, drop = FALSE]), c(log(prior.tau) - 
                  B.tau[, 1], log(prior.kappa) - B.kappa[, 1]))
          } else if(parameterization == "matern"){
              theta.prior.mean = qr.solve(rbind(B.sigma[, -1, drop = FALSE], 
                  B.range[, -1, drop = FALSE]), c(log(prior.std.dev.nominal) - 
                  B.sigma[, 1], log(prior.range.nominal) - B.range[, 1]))
          }
        }
        else {
            theta.prior.mean = rep(0, n.theta)
        }
    }
    param = list(B.tau = B.tau, 
        B.kappa = B.kappa, theta.prior.mean = theta.prior.mean, 
        theta.prior.prec = theta.prior.prec)
    return(param)
}



#' @noRd 

# Function to convert B.sigma and B.range to B.tau and B.kappa

convert_B_matrices <- function(B.sigma, B.range, n.spde, nu.nominal, d){
      n.theta <- ncol(B.sigma) - 1L

      alpha.nominal <- nu.nominal + d / 2
      
      B.sigma <- prepare_B_matrices(B.sigma, n.spde, 
          n.theta)
      B.range <- prepare_B_matrices(B.range, n.spde, 
          n.theta)

      B.kappa <- cbind(0.5 * log(8 * nu.nominal) - B.range[, 1], 
        -B.range[, -1, drop = FALSE])

      B.tau <- cbind(0.5 * (lgamma(nu.nominal) - lgamma(alpha.nominal) - 
                d/2 * log(4 * pi)) - nu.nominal * B.kappa[, 1] - 
                B.sigma[,1], 
                - nu.nominal * B.kappa[, -1, drop = FALSE] -
                B.sigma[, -1, drop = FALSE])
    
    return(list(B.tau = B.tau, B.kappa = B.kappa))
}

#' @noRd 
# Change parameterization in rspde_lme to matern

change_parameterization_lme <- function(likelihood, d, nu, par, hessian
#, improve_gradient, gradient_args
){
  tau <- par[1]
  kappa <- par[2]

  C1 <- sqrt(8*nu)
  C2 <- sqrt(gamma(nu) / ((4 * pi)^(d / 2) * gamma(nu + d / 2)))

  sigma <- C2/(tau * kappa^nu)
  range <- C1/kappa

  grad_par <- matrix(c(-C2/(kappa^nu * sigma^2),0,
                    nu * range^(nu-1) * C2/(sigma * C1^nu), 
                    -C1/range^2), nrow = 2, ncol=2)
  

  new_observed_fisher <- t(grad_par) %*% hessian %*% (grad_par)

  # No need to include the additional term as the gradient is approximately zero.
  # from some numerical experiments, the approximation without the additional term
  # seems to be better in general.

  # hess_par <- matrix(c(2*C2/(kappa^nu * sigma^3), 0,
  #                     -nu * C2/((sigma^2) * (C1^nu)) * range^(nu-1),
  #                     2*C1/range^3) , ncol=2, nrow=2)
  
  # if(!improve_gradient){
  #   grad_lik <- numDeriv::grad(likelihood, log(par), method = "simple", method.args = gradient_args)
  # } else{
  #   grad_lik <- numDeriv::grad(likelihood, log(par), method = "Richardson", method.args = gradient_args)
  # }

  # grad_lik <- c(1/tau, 1/kappa) * grad_lik

  # add_mat <- diag(grad_lik) %*% hess_par

  # add_mat <- 0.5 * (add_mat + t(add_mat))

  # new_observed_fisher <- new_observed_fisher + add_mat

  inv_fisher <- tryCatch(solve(new_observed_fisher), error = function(e) matrix(NA, nrow(new_observed_fisher), ncol(new_observed_fisher)))
  
  std_err <- sqrt(diag(inv_fisher))

  # new_lik <- function(theta){
  #       sigma <- exp(theta[1])
  #       range <- exp(theta[2])

  #       kappa <- C1/range
  #       tau <- C2/(sigma * kappa^nu)
  #       return(likelihood(log(c(tau,kappa))))
  # }

  # hess_tmp <- numDeriv::hessian(new_lik, log(c(sigma,range)))

  # hess_tmp <- diag(c(1/sigma, 1/range)) %*% hess_tmp %*% diag(c(1/sigma, 1/range))

  return(list(coeff = c(sigma, range), std_random = std_err))
}



#' @noRd 
#' 

return_same_input_type_matrix_vector <- function(v, orig_v){
  if(isS4(orig_v)){
    return(v)
  } else{
    v_out <- as.matrix(v)
    dim(v_out) <- dim(orig_v)
    return(v_out)
  }
}



#' find indices of the rows with all NA's in lists
#' @noRd
#'
idx_not_all_NA <- function(data_list){
     data_list[[".edge_number"]] <- NULL
     data_list[[".distance_on_edge"]] <- NULL
     data_list[[".coord_x"]] <- NULL
     data_list[[".coord_y"]] <- NULL
     data_list[[".group"]] <- NULL
     data_names <- names(data_list)
     n_data <- length(data_list[[data_names[1]]])
     idx_non_na <- logical(n_data)
     for(i in 1:n_data){
        na_idx <- lapply(data_list, function(dat){
          return(is.na(dat[i]))
        })
        idx_non_na[i] <- !all(unlist(na_idx))
     }
     return(idx_non_na)
}

#' find indices of the rows with at least one NA's in lists
#' @noRd
#'
idx_not_any_NA <- function(data_list){
     data_list[[".edge_number"]] <- NULL
     data_list[[".distance_on_edge"]] <- NULL
     data_list[[".coord_x"]] <- NULL
     data_list[[".coord_y"]] <- NULL
     data_list[[".group"]] <- NULL
     data_names <- names(data_list)
     n_data <- length(data_list[[data_names[1]]])
     idx_non_na <- logical(n_data)
     for(i in 1:n_data){
        na_idx <- lapply(data_list, function(dat){
          return(is.na(dat[i]))
        })
        idx_non_na[i] <- !any(unlist(na_idx))
     }
     return(idx_non_na)
}


#' @noRd 
#' 

select_indexes <- function(data, idx){
    if(inherits(data,"SpatialPointsDataFrame")){
      data <- data[idx,]
  } else{
      data <- lapply(data, function(dat){
        if(is.null(nrow(dat))){
          return(dat[idx])
        } else{
          return(dat[idx,])
        }
        })
  }
  return(data)
}

Try the rSPDE package in your browser

Any scripts or data that you put into this service are public.

rSPDE documentation built on Nov. 6, 2023, 1:06 a.m.