From 74e2e0506a038684edbf8dc9fed9aa38bc6fc408 Mon Sep 17 00:00:00 2001 From: David Bolin Date: Wed, 18 Dec 2024 14:23:06 +0100 Subject: [PATCH] adding fractional intrinsic inla model --- NAMESPACE | 5 + R/inla_rspde.R | 887 +++++++++--------- R/inla_rspde_intrinsic.R | 150 +++ R/inla_rspde_spacetime.R | 41 +- R/intrinsic.R | 109 ++- R/spacetime.R | 70 +- R/util.R | 77 +- man/bru_get_mapper.inla_rspde.Rd | 8 +- ...bru_get_mapper.inla_rspde_anisotropic2d.Rd | 2 +- man/bru_get_mapper.inla_rspde_matern1d.Rd | 8 +- man/bru_get_mapper.inla_rspde_spacetime.Rd | 2 +- man/intrinsic.operators.Rd | 101 ++ man/rSPDE.Ast.Rd | 5 +- man/rspde.intrinsic.Rd | 84 ++ man/rspde.spacetime.Rd | 3 + man/spacetime.operators.Rd | 7 +- src/cgeneric_aux_anisotropic.cpp | 2 +- src/cgeneric_aux_intrinsic.cpp | 201 ++-- src/cgeneric_defs.h | 4 +- src/cgeneric_rspde_fintrinsic.c | 333 +++++++ vignettes/intrinsic.Rmd | 101 +- 21 files changed, 1604 insertions(+), 596 deletions(-) create mode 100644 man/intrinsic.operators.Rd create mode 100644 man/rspde.intrinsic.Rd create mode 100644 src/cgeneric_rspde_fintrinsic.c diff --git a/NAMESPACE b/NAMESPACE index 42a9404f..7667b2c9 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -74,6 +74,7 @@ export(gg_df) export(glance) export(graph_data_rspde) export(intrinsic.matern.operators) +export(intrinsic.operators) export(matern.covariance) export(matern.operators) export(matern.rational) @@ -137,6 +138,10 @@ if (getRversion() >= "3.6.0") { S3method(inlabru::ibm_n, bru_mapper_metric_graph) S3method(inlabru::ibm_values, bru_mapper_metric_graph) S3method(inlabru::ibm_jacobian, bru_mapper_metric_graph) + S3method(inlabru::bru_mapper, metric_graph_dirichlet) + S3method(inlabru::ibm_n, bru_mapper_metric_graph_dirichlet) + S3method(inlabru::ibm_values, bru_mapper_metric_graph_dirichlet) + S3method(inlabru::ibm_jacobian, bru_mapper_metric_graph_dirichlet) } import(Matrix) diff --git a/R/inla_rspde.R b/R/inla_rspde.R index e2a03d1b..af8b7dac 100644 --- a/R/inla_rspde.R +++ b/R/inla_rspde.R @@ -1742,465 +1742,472 @@ graph_index_rspde <- function(graph_spde, n.repl = 1, n.group = 1) { #' } #' # devel.tag #' } -rspde.result <- function(inla, name, rspde, compute.summary = TRUE, parameterization = "detect", n_samples = 5000, n_density = 1024) { +rspde.result <- function(inla, name, rspde, compute.summary = TRUE, + parameterization = "detect", + n_samples = 5000, n_density = 1024) { check_class_inla_rspde(rspde) - stationary <- rspde$stationary - - parameterization <- parameterization[[1]] - parameterization <- tolower(parameterization) - - if (!(parameterization %in% c("detect", "spde", "matern", "matern2"))) { - stop("The possible options for parameterization are 'detect', 'spde', 'matern' and 'matern2'.") - } - - nu.upper.bound <- rspde$nu.upper.bound - result <- list() - - par_model <- rspde$parameterization - - if (parameterization == "detect") { - parameterization <- rspde$parameterization - } - - if (stationary) { - if (!rspde$est_nu) { - if (parameterization == "spde") { - row_names <- c("tau", "kappa") - } else if (parameterization == "matern") { - row_names <- c("std.dev", "range") - } else if (parameterization == "matern2") { - row_names <- c("var", "r") - } - } else { - if (parameterization == "spde") { - row_names <- c("tau", "kappa", "nu") - } else if (parameterization == "matern") { - row_names <- c("std.dev", "range", "nu") - } else if (parameterization == "matern2") { - row_names <- c("var", "r", "nu") - } - } - - - result$summary.values <- inla$summary.random[[name]] - - if (!is.null(inla$marginals.random[[name]])) { - result$marginals.values <- inla$marginals.random[[name]] - } - - if (parameterization == "spde") { - name_theta1 <- "tau" - name_theta2 <- "kappa" - } else if (parameterization == "matern") { - name_theta1 <- "std.dev" - name_theta2 <- "range" - } else if (parameterization == "matern2") { - name_theta1 <- "var" - name_theta2 <- "r" - } - - if (par_model == "spde") { - name_theta1_model <- "tau" - name_theta2_model <- "kappa" - } else if (par_model == "matern") { - name_theta1_model <- "std.dev" - name_theta2_model <- "range" - } else if (par_model == "matern2") { - name_theta1_model <- "var" - name_theta2_model <- "r" - } - - result[[paste0("summary.log.", name_theta1_model)]] <- INLA::inla.extract.el( - inla$summary.hyperpar, - paste("Theta1 for ", name, "$", sep = "") - ) - rownames(result[[paste0("summary.log.", name_theta1_model)]]) <- paste0("log(", name_theta1_model, ")") - - result[[paste0("summary.log.", name_theta2_model)]] <- INLA::inla.extract.el( - inla$summary.hyperpar, - paste("Theta2 for ", name, "$", sep = "") - ) - rownames(result[[paste0("summary.log.", name_theta2_model)]]) <- paste0("log(", name_theta2_model, ")") - if (rspde$est_nu) { - result$summary.logit.nu <- INLA::inla.extract.el( - inla$summary.hyperpar, - paste("Theta3 for ", name, "$", sep = "") - ) - rownames(result$summary.logit.nu) <- "logit(nu)" - } - - if (!is.null(inla$marginals.hyperpar[[paste0("Theta1 for ", name)]])) { - result[[paste0("marginals.log.", name_theta1_model)]] <- INLA::inla.extract.el( - inla$marginals.hyperpar, - paste("Theta1 for ", name, "$", sep = "") - ) - names(result[[paste0("marginals.log.", name_theta1_model)]]) <- name_theta1_model - result[[paste0("marginals.log.", name_theta2_model)]] <- INLA::inla.extract.el( - inla$marginals.hyperpar, - paste("Theta2 for ", name, "$", sep = "") - ) - names(result[[paste0("marginals.log.", name_theta2_model)]]) <- name_theta2_model - - if (rspde$est_nu) { - result$marginals.logit.nu <- INLA::inla.extract.el( - inla$marginals.hyperpar, - paste("Theta3 for ", name, "$", sep = "") - ) - names(result$marginals.logit.nu) <- "nu" - } - - - if (par_model == parameterization) { - result[[paste0("marginals.", name_theta1)]] <- lapply( - result[[paste0("marginals.log.", name_theta1)]], - function(x) { - INLA::inla.tmarginal( - function(y) exp(y), - x - ) + if(inherits(rspde, "inla_rspde_fintrinsic")){ + return(rspde.intrinsic.result(inla = inla, + name = name, + rspde = rspde, + compute.summary = compute.summary, + n_samples = n_samples, + n_density = n_density)) + } else { + stationary <- rspde$stationary + + parameterization <- parameterization[[1]] + parameterization <- tolower(parameterization) + + if (!(parameterization %in% c("detect", "spde", "matern", "matern2"))) { + stop("The possible options for parameterization are 'detect', 'spde', 'matern' and 'matern2'.") + } + + nu.upper.bound <- rspde$nu.upper.bound + result <- list() + + par_model <- rspde$parameterization + + if (parameterization == "detect") { + parameterization <- rspde$parameterization + } + + if (stationary) { + if (!rspde$est_nu) { + if (parameterization == "spde") { + row_names <- c("tau", "kappa") + } else if (parameterization == "matern") { + row_names <- c("std.dev", "range") + } else if (parameterization == "matern2") { + row_names <- c("var", "r") + } + } else { + if (parameterization == "spde") { + row_names <- c("tau", "kappa", "nu") + } else if (parameterization == "matern") { + row_names <- c("std.dev", "range", "nu") + } else if (parameterization == "matern2") { + row_names <- c("var", "r", "nu") + } } - ) - result[[paste0("marginals.", name_theta2)]] <- lapply( - result[[paste0("marginals.log.", name_theta2)]], - function(x) { - INLA::inla.tmarginal( - function(y) exp(y), - x - ) + + + result$summary.values <- inla$summary.random[[name]] + + if (!is.null(inla$marginals.random[[name]])) { + result$marginals.values <- inla$marginals.random[[name]] } - ) - if (rspde$est_nu) { - result$marginals.nu <- lapply( - result$marginals.logit.nu, - function(x) { - INLA::inla.tmarginal( - function(y) { - nu.upper.bound * exp(y) / (1 + exp(y)) - }, - x - ) - } + + if (parameterization == "spde") { + name_theta1 <- "tau" + name_theta2 <- "kappa" + } else if (parameterization == "matern") { + name_theta1 <- "std.dev" + name_theta2 <- "range" + } else if (parameterization == "matern2") { + name_theta1 <- "var" + name_theta2 <- "r" + } + + if (par_model == "spde") { + name_theta1_model <- "tau" + name_theta2_model <- "kappa" + } else if (par_model == "matern") { + name_theta1_model <- "std.dev" + name_theta2_model <- "range" + } else if (par_model == "matern2") { + name_theta1_model <- "var" + name_theta2_model <- "r" + } + + result[[paste0("summary.log.", name_theta1_model)]] <- INLA::inla.extract.el( + inla$summary.hyperpar, + paste("Theta1 for ", name, "$", sep = "") ) - } - } else { - if (par_model == "spde") { - dim <- rspde$dim - hyperpar_sample <- INLA::inla.hyperpar.sample(n_samples, inla) + rownames(result[[paste0("summary.log.", name_theta1_model)]]) <- paste0("log(", name_theta1_model, ")") + + result[[paste0("summary.log.", name_theta2_model)]] <- INLA::inla.extract.el( + inla$summary.hyperpar, + paste("Theta2 for ", name, "$", sep = "") + ) + rownames(result[[paste0("summary.log.", name_theta2_model)]]) <- paste0("log(", name_theta2_model, ")") if (rspde$est_nu) { - nu_est <- rspde$nu.upper.bound * exp(hyperpar_sample[, paste0("Theta3 for ", name)]) / (1 + exp(hyperpar_sample[, paste0("Theta3 for ", name)])) - } else { - nu_est <- rspde[["nu"]] + result$summary.logit.nu <- INLA::inla.extract.el( + inla$summary.hyperpar, + paste("Theta3 for ", name, "$", sep = "") + ) + rownames(result$summary.logit.nu) <- "logit(nu)" } - tau_est <- exp(hyperpar_sample[, paste0("Theta1 for ", name)]) - kappa_est <- exp(hyperpar_sample[, paste0("Theta2 for ", name)]) - - sigma_est <- sqrt(gamma(0.5) / (tau_est^2 * kappa_est^(2 * nu_est) * - (4 * pi)^(dim / 2) * gamma(nu_est + dim / 2))) - - if (parameterization == "matern") { - range_est <- sqrt(8 * nu_est) / kappa_est - density_theta1 <- stats::density(sigma_est, n = n_density) - density_theta2 <- stats::density(range_est, n = n_density) - } else if (parameterization == "matern2") { - var_est <- sigma_est^2 - r_est <- 1 / kappa_est - density_theta1 <- stats::density(var_est, n = n_density) - density_theta2 <- stats::density(r_est, n = n_density) + + if (!is.null(inla$marginals.hyperpar[[paste0("Theta1 for ", name)]])) { + result[[paste0("marginals.log.", name_theta1_model)]] <- INLA::inla.extract.el( + inla$marginals.hyperpar, + paste("Theta1 for ", name, "$", sep = "") + ) + names(result[[paste0("marginals.log.", name_theta1_model)]]) <- name_theta1_model + result[[paste0("marginals.log.", name_theta2_model)]] <- INLA::inla.extract.el( + inla$marginals.hyperpar, + paste("Theta2 for ", name, "$", sep = "") + ) + names(result[[paste0("marginals.log.", name_theta2_model)]]) <- name_theta2_model + + if (rspde$est_nu) { + result$marginals.logit.nu <- INLA::inla.extract.el( + inla$marginals.hyperpar, + paste("Theta3 for ", name, "$", sep = "") + ) + names(result$marginals.logit.nu) <- "nu" + } + + + if (par_model == parameterization) { + result[[paste0("marginals.", name_theta1)]] <- lapply( + result[[paste0("marginals.log.", name_theta1)]], + function(x) { + INLA::inla.tmarginal( + function(y) exp(y), + x + ) + } + ) + result[[paste0("marginals.", name_theta2)]] <- lapply( + result[[paste0("marginals.log.", name_theta2)]], + function(x) { + INLA::inla.tmarginal( + function(y) exp(y), + x + ) + } + ) + if (rspde$est_nu) { + result$marginals.nu <- lapply( + result$marginals.logit.nu, + function(x) { + INLA::inla.tmarginal( + function(y) { + nu.upper.bound * exp(y) / (1 + exp(y)) + }, + x + ) + } + ) + } + } else { + if (par_model == "spde") { + dim <- rspde$dim + hyperpar_sample <- INLA::inla.hyperpar.sample(n_samples, inla) + if (rspde$est_nu) { + nu_est <- rspde$nu.upper.bound * exp(hyperpar_sample[, paste0("Theta3 for ", name)]) / (1 + exp(hyperpar_sample[, paste0("Theta3 for ", name)])) + } else { + nu_est <- rspde[["nu"]] + } + tau_est <- exp(hyperpar_sample[, paste0("Theta1 for ", name)]) + kappa_est <- exp(hyperpar_sample[, paste0("Theta2 for ", name)]) + + sigma_est <- sqrt(gamma(0.5) / (tau_est^2 * kappa_est^(2 * nu_est) * + (4 * pi)^(dim / 2) * gamma(nu_est + dim / 2))) + + if (parameterization == "matern") { + range_est <- sqrt(8 * nu_est) / kappa_est + density_theta1 <- stats::density(sigma_est, n = n_density) + density_theta2 <- stats::density(range_est, n = n_density) + } else if (parameterization == "matern2") { + var_est <- sigma_est^2 + r_est <- 1 / kappa_est + density_theta1 <- stats::density(var_est, n = n_density) + density_theta2 <- stats::density(r_est, n = n_density) + } + + result[[paste0("marginals.", name_theta1)]] <- list() + result[[paste0("marginals.", name_theta1)]][[name_theta1]] <- cbind(density_theta1$x, density_theta1$y) + colnames(result[[paste0("marginals.", name_theta1)]][[name_theta1]]) <- c("x", "y") + + result[[paste0("marginals.", name_theta2)]] <- list() + result[[paste0("marginals.", name_theta2)]][[name_theta2]] <- cbind(density_theta2$x, density_theta2$y) + colnames(result[[paste0("marginals.", name_theta2)]][[name_theta2]]) <- c("x", "y") + } else if (par_model == "matern") { + hyperpar_sample <- INLA::inla.hyperpar.sample(n_samples, inla) + dim <- rspde$dim + if (rspde$est_nu) { + nu_est <- rspde$nu.upper.bound * exp(hyperpar_sample[, paste0("Theta3 for ", name)]) / (1 + exp(hyperpar_sample[, paste0("Theta3 for ", name)])) + } else { + nu_est <- rspde[["nu"]] + } + sigma_est <- exp(hyperpar_sample[, paste0("Theta1 for ", name)]) + range_est <- exp(hyperpar_sample[, paste0("Theta2 for ", name)]) + + kappa_est <- sqrt(8 * nu_est) / range_est + + if (parameterization == "spde") { + tau_est <- sqrt(gamma(0.5) / (sigma_est^2 * kappa_est^(2 * nu_est) * + (4 * pi)^(dim / 2) * gamma(nu_est + dim / 2))) + density_theta1 <- stats::density(tau_est, n = n_density) + density_theta2 <- stats::density(kappa_est, n = n_density) + } else if (parameterization == "matern2") { + var_est <- sigma_est^2 + r_est <- 1 / kappa_est + density_theta1 <- stats::density(var_est, n = n_density) + density_theta2 <- stats::density(r_est, n = n_density) + } + + result[[paste0("marginals.", name_theta1)]] <- list() + result[[paste0("marginals.", name_theta1)]][[name_theta1]] <- cbind(density_theta1$x, density_theta1$y) + colnames(result[[paste0("marginals.", name_theta1)]][[name_theta1]]) <- c("x", "y") + + result[[paste0("marginals.", name_theta2)]] <- list() + result[[paste0("marginals.", name_theta2)]][[name_theta2]] <- cbind(density_theta2$x, density_theta2$y) + colnames(result[[paste0("marginals.", name_theta2)]][[name_theta2]]) <- c("x", "y") + } else if (par_model == "matern2") { + hyperpar_sample <- INLA::inla.hyperpar.sample(n_samples, inla) + dim <- rspde$dim + if (rspde$est_nu) { + nu_est <- rspde$nu.upper.bound * exp(hyperpar_sample[, paste0("Theta3 for ", name)]) / (1 + exp(hyperpar_sample[, paste0("Theta3 for ", name)])) + } else { + nu_est <- rspde[["nu"]] + } + var_est <- exp(hyperpar_sample[, paste0("Theta1 for ", name)]) + r_est <- exp(hyperpar_sample[, paste0("Theta2 for ", name)]) + + kappa_est <- 1 / r_est + sigma_est <- sqrt(var_est) + + if (parameterization == "spde") { + tau_est <- sqrt(gamma(0.5) / (sigma_est^2 * kappa_est^(2 * nu_est) * + (4 * pi)^(dim / 2) * gamma(nu_est + dim / 2))) + density_theta1 <- stats::density(tau_est, n = n_density) + density_theta2 <- stats::density(kappa_est, n = n_density) + } else if (parameterization == "matern") { + range_est <- sqrt(8 * nu_est) / kappa_est + density_theta1 <- stats::density(sigma_est, n = n_density) + density_theta2 <- stats::density(range_est, n = n_density) + } + + result[[paste0("marginals.", name_theta1)]] <- list() + result[[paste0("marginals.", name_theta1)]][[name_theta1]] <- cbind(density_theta1$x, density_theta1$y) + colnames(result[[paste0("marginals.", name_theta1)]][[name_theta1]]) <- c("x", "y") + + result[[paste0("marginals.", name_theta2)]] <- list() + result[[paste0("marginals.", name_theta2)]][[name_theta2]] <- cbind(density_theta2$x, density_theta2$y) + colnames(result[[paste0("marginals.", name_theta2)]][[name_theta2]]) <- c("x", "y") + } + + if (rspde$est_nu) { + result$marginals.nu <- lapply( + result$marginals.logit.nu, + function(x) { + INLA::inla.tmarginal( + function(y) { + nu.upper.bound * exp(y) / (1 + exp(y)) + }, + x + ) + } + ) + } + } } - - result[[paste0("marginals.", name_theta1)]] <- list() - result[[paste0("marginals.", name_theta1)]][[name_theta1]] <- cbind(density_theta1$x, density_theta1$y) - colnames(result[[paste0("marginals.", name_theta1)]][[name_theta1]]) <- c("x", "y") - - result[[paste0("marginals.", name_theta2)]] <- list() - result[[paste0("marginals.", name_theta2)]][[name_theta2]] <- cbind(density_theta2$x, density_theta2$y) - colnames(result[[paste0("marginals.", name_theta2)]][[name_theta2]]) <- c("x", "y") - } else if (par_model == "matern") { - hyperpar_sample <- INLA::inla.hyperpar.sample(n_samples, inla) - dim <- rspde$dim - if (rspde$est_nu) { - nu_est <- rspde$nu.upper.bound * exp(hyperpar_sample[, paste0("Theta3 for ", name)]) / (1 + exp(hyperpar_sample[, paste0("Theta3 for ", name)])) + + if (compute.summary) { + norm_const <- function(density_df) { + 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) + } + norm_const <- stats::integrate( + f = function(z) { + denstemp(z) + }, lower = min_x, upper = max_x, + subdivisions = nrow(density_df), + stop.on.error = FALSE + )$value + return(norm_const) + } + + norm_const_theta1 <- norm_const(result[[paste0("marginals.", name_theta1)]][[name_theta1]]) + result[[paste0("marginals.", name_theta1)]][[name_theta1]][, "y"] <- + result[[paste0("marginals.", name_theta1)]][[name_theta1]][, "y"] / norm_const_theta1 + + norm_const_theta2 <- norm_const(result[[paste0("marginals.", name_theta2)]][[name_theta2]]) + result[[paste0("marginals.", name_theta2)]][[name_theta2]][, "y"] <- + result[[paste0("marginals.", name_theta2)]][[name_theta2]][, "y"] / norm_const_theta2 + + result[[paste0("summary.", name_theta1)]] <- create_summary_from_density(result[[paste0("marginals.", name_theta1)]][[name_theta1]], + name = name_theta1 + ) + + result[[paste0("summary.", name_theta2)]] <- + create_summary_from_density(result[[paste0("marginals.", name_theta2)]][[name_theta2]], name = name_theta2) + + if (rspde$est_nu) { + norm_const_nu <- norm_const(result$marginals.nu$nu) + result$marginals.nu$nu[, "y"] <- + result$marginals.nu$nu[, "y"] / norm_const_nu + + result$summary.nu <- create_summary_from_density(result$marginals.nu$nu, + name = "nu" + ) + } + } + } else { + n_par <- length(rspde$start.theta) + + if (!rspde$est_nu) { + if (parameterization == "spde") { + row_names <- sapply(1:n_par, function(i) { + paste0("Theta", i, ".spde") + }) + } else if (parameterization == "matern") { + row_names <- sapply(1:n_par, function(i) { + paste0("Theta", i, ".matern") + }) + } else if (parameterization == "matern2") { + row_names <- sapply(1:n_par, function(i) { + paste0("Theta", i, ".matern2") + }) + } } else { - nu_est <- rspde[["nu"]] + if (parameterization == "spde") { + row_names <- sapply(1:n_par, function(i) { + paste0("Theta", i, ".spde") + }) + row_names <- c(row_names, "nu") + } else if (parameterization == "matern") { + row_names <- sapply(1:n_par, function(i) { + paste0("Theta", i, ".matern") + }) + row_names <- c(row_names, "nu") + } else if (parameterization == "matern2") { + row_names <- sapply(1:n_par, function(i) { + paste0("Theta", i, ".matern2") + }) + row_names <- c(row_names, "nu") + } } - sigma_est <- exp(hyperpar_sample[, paste0("Theta1 for ", name)]) - range_est <- exp(hyperpar_sample[, paste0("Theta2 for ", name)]) - - kappa_est <- sqrt(8 * nu_est) / range_est - - if (parameterization == "spde") { - tau_est <- sqrt(gamma(0.5) / (sigma_est^2 * kappa_est^(2 * nu_est) * - (4 * pi)^(dim / 2) * gamma(nu_est + dim / 2))) - density_theta1 <- stats::density(tau_est, n = n_density) - density_theta2 <- stats::density(kappa_est, n = n_density) - } else if (parameterization == "matern2") { - var_est <- sigma_est^2 - r_est <- 1 / kappa_est - density_theta1 <- stats::density(var_est, n = n_density) - density_theta2 <- stats::density(r_est, n = n_density) + + for (i in 1:n_par) { + result[[paste0("summary.", row_names[i])]] <- INLA::inla.extract.el( + inla$summary.hyperpar, + paste0("Theta", i, " for ", name, "$", sep = "") + ) + rownames(result[[paste0("summary.", row_names[i])]]) <- row_names[i] } - - result[[paste0("marginals.", name_theta1)]] <- list() - result[[paste0("marginals.", name_theta1)]][[name_theta1]] <- cbind(density_theta1$x, density_theta1$y) - colnames(result[[paste0("marginals.", name_theta1)]][[name_theta1]]) <- c("x", "y") - - result[[paste0("marginals.", name_theta2)]] <- list() - result[[paste0("marginals.", name_theta2)]][[name_theta2]] <- cbind(density_theta2$x, density_theta2$y) - colnames(result[[paste0("marginals.", name_theta2)]][[name_theta2]]) <- c("x", "y") - } else if (par_model == "matern2") { - hyperpar_sample <- INLA::inla.hyperpar.sample(n_samples, inla) - dim <- rspde$dim + if (rspde$est_nu) { - nu_est <- rspde$nu.upper.bound * exp(hyperpar_sample[, paste0("Theta3 for ", name)]) / (1 + exp(hyperpar_sample[, paste0("Theta3 for ", name)])) - } else { - nu_est <- rspde[["nu"]] + result$summary.logit.nu <- INLA::inla.extract.el( + inla$summary.hyperpar, + paste0("Theta", n_par + 1, " for ", name, "$", sep = "") + ) + rownames(result$summary.logit.nu) <- "logit(nu)" } - var_est <- exp(hyperpar_sample[, paste0("Theta1 for ", name)]) - r_est <- exp(hyperpar_sample[, paste0("Theta2 for ", name)]) - - kappa_est <- 1 / r_est - sigma_est <- sqrt(var_est) - - if (parameterization == "spde") { - tau_est <- sqrt(gamma(0.5) / (sigma_est^2 * kappa_est^(2 * nu_est) * - (4 * pi)^(dim / 2) * gamma(nu_est + dim / 2))) - density_theta1 <- stats::density(tau_est, n = n_density) - density_theta2 <- stats::density(kappa_est, n = n_density) - } else if (parameterization == "matern") { - range_est <- sqrt(8 * nu_est) / kappa_est - density_theta1 <- stats::density(sigma_est, n = n_density) - density_theta2 <- stats::density(range_est, n = n_density) + + + for (i in 1:n_par) { + if (!is.null(inla$marginals.hyperpar[[paste0("Theta", i, " for ", name)]])) { + result[[paste0("marginals.", row_names[i])]] <- INLA::inla.extract.el( + inla$marginals.hyperpar, + paste0("Theta", i, " for ", name, "$", sep = "") + ) + names(result[[paste0("marginals.", row_names[i])]]) <- row_names[i] + } } - - result[[paste0("marginals.", name_theta1)]] <- list() - result[[paste0("marginals.", name_theta1)]][[name_theta1]] <- cbind(density_theta1$x, density_theta1$y) - colnames(result[[paste0("marginals.", name_theta1)]][[name_theta1]]) <- c("x", "y") - - result[[paste0("marginals.", name_theta2)]] <- list() - result[[paste0("marginals.", name_theta2)]][[name_theta2]] <- cbind(density_theta2$x, density_theta2$y) - colnames(result[[paste0("marginals.", name_theta2)]][[name_theta2]]) <- c("x", "y") - } - - if (rspde$est_nu) { - result$marginals.nu <- lapply( - result$marginals.logit.nu, - function(x) { - INLA::inla.tmarginal( - function(y) { - nu.upper.bound * exp(y) / (1 + exp(y)) - }, - x + + + if (rspde$est_nu) { + result$marginals.logit.nu <- INLA::inla.extract.el( + inla$marginals.hyperpar, + paste0("Theta", n_par + 1, " for ", name, "$", sep = "") ) - } - ) - } - } - } - - - - - - if (compute.summary) { - norm_const <- function(density_df) { - 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) - } - norm_const <- stats::integrate( - f = function(z) { - denstemp(z) - }, lower = min_x, upper = max_x, - subdivisions = nrow(density_df), - stop.on.error = FALSE - )$value - return(norm_const) - } - - norm_const_theta1 <- norm_const(result[[paste0("marginals.", name_theta1)]][[name_theta1]]) - result[[paste0("marginals.", name_theta1)]][[name_theta1]][, "y"] <- - result[[paste0("marginals.", name_theta1)]][[name_theta1]][, "y"] / norm_const_theta1 - - norm_const_theta2 <- norm_const(result[[paste0("marginals.", name_theta2)]][[name_theta2]]) - result[[paste0("marginals.", name_theta2)]][[name_theta2]][, "y"] <- - result[[paste0("marginals.", name_theta2)]][[name_theta2]][, "y"] / norm_const_theta2 - - result[[paste0("summary.", name_theta1)]] <- create_summary_from_density(result[[paste0("marginals.", name_theta1)]][[name_theta1]], - name = name_theta1 - ) - - result[[paste0("summary.", name_theta2)]] <- - create_summary_from_density(result[[paste0("marginals.", name_theta2)]][[name_theta2]], name = name_theta2) - - if (rspde$est_nu) { - norm_const_nu <- norm_const(result$marginals.nu$nu) - result$marginals.nu$nu[, "y"] <- - result$marginals.nu$nu[, "y"] / norm_const_nu - - result$summary.nu <- create_summary_from_density(result$marginals.nu$nu, - name = "nu" - ) - } - } - } else { - n_par <- length(rspde$start.theta) - - if (!rspde$est_nu) { - if (parameterization == "spde") { - row_names <- sapply(1:n_par, function(i) { - paste0("Theta", i, ".spde") - }) - } else if (parameterization == "matern") { - row_names <- sapply(1:n_par, function(i) { - paste0("Theta", i, ".matern") - }) - } else if (parameterization == "matern2") { - row_names <- sapply(1:n_par, function(i) { - paste0("Theta", i, ".matern2") - }) - } - } else { - if (parameterization == "spde") { - row_names <- sapply(1:n_par, function(i) { - paste0("Theta", i, ".spde") - }) - row_names <- c(row_names, "nu") - } else if (parameterization == "matern") { - row_names <- sapply(1:n_par, function(i) { - paste0("Theta", i, ".matern") - }) - row_names <- c(row_names, "nu") - } else if (parameterization == "matern2") { - row_names <- sapply(1:n_par, function(i) { - paste0("Theta", i, ".matern2") - }) - row_names <- c(row_names, "nu") - } - } - - for (i in 1:n_par) { - result[[paste0("summary.", row_names[i])]] <- INLA::inla.extract.el( - inla$summary.hyperpar, - paste0("Theta", i, " for ", name, "$", sep = "") - ) - rownames(result[[paste0("summary.", row_names[i])]]) <- row_names[i] - } - - if (rspde$est_nu) { - result$summary.logit.nu <- INLA::inla.extract.el( - inla$summary.hyperpar, - paste0("Theta", n_par + 1, " for ", name, "$", sep = "") - ) - rownames(result$summary.logit.nu) <- "logit(nu)" - } - - - for (i in 1:n_par) { - if (!is.null(inla$marginals.hyperpar[[paste0("Theta", i, " for ", name)]])) { - result[[paste0("marginals.", row_names[i])]] <- INLA::inla.extract.el( - inla$marginals.hyperpar, - paste0("Theta", i, " for ", name, "$", sep = "") - ) - names(result[[paste0("marginals.", row_names[i])]]) <- row_names[i] - } - } - - - if (rspde$est_nu) { - result$marginals.logit.nu <- INLA::inla.extract.el( - inla$marginals.hyperpar, - paste0("Theta", n_par + 1, " for ", name, "$", sep = "") - ) - names(result$marginals.logit.nu) <- "nu" - - result$marginals.nu <- lapply( - result$marginals.logit.nu, - function(x) { - INLA::inla.tmarginal( - function(y) { - nu.upper.bound * exp(y) / (1 + exp(y)) - }, - x - ) - } - ) - } - - if (compute.summary) { - norm_const <- function(density_df) { - 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) - } - norm_const <- stats::integrate( - f = function(z) { - denstemp(z) - }, lower = min_x, upper = max_x, - stop.on.error = FALSE - )$value - return(norm_const) + names(result$marginals.logit.nu) <- "nu" + + result$marginals.nu <- lapply( + result$marginals.logit.nu, + function(x) { + INLA::inla.tmarginal( + function(y) { + nu.upper.bound * exp(y) / (1 + exp(y)) + }, + x + ) + } + ) + } + + if (compute.summary) { + norm_const <- function(density_df) { + 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) + } + norm_const <- stats::integrate( + f = function(z) { + denstemp(z) + }, lower = min_x, upper = max_x, + stop.on.error = FALSE + )$value + return(norm_const) + } + + if (rspde$est_nu) { + norm_const_nu <- norm_const(result$marginals.nu$nu) + result$marginals.nu$nu[, "y"] <- + result$marginals.nu$nu[, "y"] / norm_const_nu + + result$summary.nu <- create_summary_from_density(result$marginals.nu$nu, + name = "nu" + ) + } + } } - - if (rspde$est_nu) { - norm_const_nu <- norm_const(result$marginals.nu$nu) - result$marginals.nu$nu[, "y"] <- - result$marginals.nu$nu[, "y"] / norm_const_nu - - result$summary.nu <- create_summary_from_density(result$marginals.nu$nu, - name = "nu" - ) + + result$n_par <- length(rspde$start.theta) + + class(result) <- "rspde_result" + result$stationary <- stationary + result$parameterization <- parameterization + if (stationary) { + result$params <- c(name_theta1, name_theta2) + if (rspde$est_nu) { + result$params <- c(result$params, "nu") + } + } else { + result$params <- row_names } - } - } - - result$n_par <- length(rspde$start.theta) - - class(result) <- "rspde_result" - result$stationary <- stationary - result$parameterization <- parameterization - if (stationary) { - result$params <- c(name_theta1, name_theta2) - if (rspde$est_nu) { - result$params <- c(result$params, "nu") - } - } else { - result$params <- row_names - } - - if (!is.null(result$summary.nu)) { - if(nu.upper.bound - result$summary.nu$mean < 0.1 || nu.upper.bound - result$summary.nu$mode < 0.1){ - warning("the mean or mode of nu is very close to nu.upper.bound, please consider increasing nu.upper.bound, and refitting the model.") + + if (!is.null(result$summary.nu)) { + if(nu.upper.bound - result$summary.nu$mean < 0.1 || nu.upper.bound - result$summary.nu$mode < 0.1){ + warning("the mean or mode of nu is very close to nu.upper.bound, please consider increasing nu.upper.bound, and refitting the model.") + } } + + return(result) } - - return(result) } diff --git a/R/inla_rspde_intrinsic.R b/R/inla_rspde_intrinsic.R index 7612ceea..93256838 100644 --- a/R/inla_rspde_intrinsic.R +++ b/R/inla_rspde_intrinsic.R @@ -26,6 +26,7 @@ #' The current options are "beta" or "lognormal". The default is "lognormal". #' @param nu.prec.inc Amount to increase the precision in the beta prior #' distribution. Check details below. +#' @param diagonal Number added to diagonal of Q for increased stability. #' @param type.rational.approx Which type of rational approximation #' should be used? The current types are "chebfun", "brasil" or "chebfunLB". #' @param shared_lib String specifying which shared library to use for the Cgeneric @@ -54,6 +55,7 @@ rspde.intrinsic <- function(mesh, prior.nu = NULL, prior.nu.dist = "lognormal", nu.prec.inc = 0.01, + diagonal = 1e-5, type.rational.approx = "chebfun", shared_lib = "detect", debug = FALSE, @@ -119,6 +121,8 @@ rspde.intrinsic <- function(mesh, return(out) } Cmatrix <- to.inla.matrix(op$Q) + D <- Diagonal(dim(op$Q)[1], diagonal) + scaling <- RSpectra::eigs(as(G, "CsparseMatrix"), 2, which = "SM")$values[1] list_args <- list( @@ -138,11 +142,13 @@ rspde.intrinsic <- function(mesh, prior.nu.prec = prior.nu$prec, prior.nu.logscale = prior.nu$logscale, nu_upper_bound = nu.upper.bound, + scaling = scaling, prior.nu.dist = prior.nu.dist, Q = Cmatrix, C = C, Ci = Ci, G = G, + D = D, rational_table = rational_table ) @@ -444,3 +450,147 @@ rspde.intrinsic.matern <- function(mesh, model$rspde_version <- as.character(packageVersion("rSPDE")) return(model) } + + +#' result summary for intrinsic models +#' @noRd +rspde.intrinsic.result <- function(inla, name, rspde, + compute.summary = TRUE, + n_samples = 5000, + n_density = 1024) { + + nu.upper.bound <- rspde$nu_upper_bound + result <- list() + + if (!rspde$est_nu) { + row_names <- c("tau") + } else { + row_names <- c("tau", "nu") + } + + result$summary.values <- inla$summary.random[[name]] + + if (!is.null(inla$marginals.random[[name]])) { + result$marginals.values <- inla$marginals.random[[name]] + } + + name_theta1 <- "tau" + name_theta2 <- "nu" + + name_theta1_model <- "tau" + name_theta2_model <- "nu" + + result[[paste0("summary.log.", name_theta1_model)]] <- INLA::inla.extract.el( + inla$summary.hyperpar, + paste("Theta1 for ", name, "$", sep = "") + ) + rownames(result[[paste0("summary.log.", name_theta1_model)]]) <- paste0("log(", name_theta1_model, ")") + + if (rspde$est_nu) { + result$summary.logit.nu <- INLA::inla.extract.el(inla$summary.hyperpar, + paste("Theta2 for ", name, "$", sep = "")) + rownames(result$summary.logit.nu) <- "logit(nu)" + } + + if (!is.null(inla$marginals.hyperpar[[paste0("Theta1 for ", name)]])) { + result[[paste0("marginals.log.", name_theta1_model)]] <- INLA::inla.extract.el( + inla$marginals.hyperpar, + paste("Theta1 for ", name, "$", sep = "") + ) + names(result[[paste0("marginals.log.", name_theta1_model)]]) <- name_theta1_model + + if (rspde$est_nu) { + result$marginals.logit.nu <- INLA::inla.extract.el( + inla$marginals.hyperpar, + paste("Theta2 for ", name, "$", sep = "") + ) + names(result$marginals.logit.nu) <- "nu" + } + + + result[[paste0("marginals.", name_theta1)]] <- lapply( + result[[paste0("marginals.log.", name_theta1)]], + function(x) { + INLA::inla.tmarginal(function(y) exp(y), x) + } + ) + + if (rspde$est_nu) { + result$marginals.nu <- lapply(result$marginals.logit.nu, + function(x) { + INLA::inla.tmarginal( + function(y) { + nu.upper.bound * exp(y) / (1 + exp(y)) + }, x) + }) + } + + if (compute.summary) { + norm_const <- function(density_df) { + 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) + } + norm_const <- stats::integrate( + f = function(z) { + denstemp(z) + }, lower = min_x, upper = max_x, + subdivisions = nrow(density_df), + stop.on.error = FALSE + )$value + return(norm_const) + } + + norm_const_theta1 <- norm_const(result[[paste0("marginals.", name_theta1)]][[name_theta1]]) + result[[paste0("marginals.", name_theta1)]][[name_theta1]][, "y"] <- + result[[paste0("marginals.", name_theta1)]][[name_theta1]][, "y"] / norm_const_theta1 + + result[[paste0("summary.", name_theta1)]] <- create_summary_from_density(result[[paste0("marginals.", name_theta1)]][[name_theta1]], + name = name_theta1 + ) + + + if (rspde$est_nu) { + norm_const_nu <- norm_const(result$marginals.nu$nu) + result$marginals.nu$nu[, "y"] <- + result$marginals.nu$nu[, "y"] / norm_const_nu + + result$summary.nu <- create_summary_from_density(result$marginals.nu$nu, + name = "nu" + ) + } + } + } + + result$n_par <- length(rspde$start.theta) + + class(result) <- c("rspde_result", "rspde_intrinsic") + result$stationary <- TRUE + result$parameterization <- "spde" + + result$params <- c(name_theta1) + if (rspde$est_nu) { + result$params <- c(result$params, "nu") + } + + if (!is.null(result$summary.nu)) { + if(nu.upper.bound - result$summary.nu$mean < 0.1 || nu.upper.bound - result$summary.nu$mode < 0.1){ + warning("the mean or mode of nu is very close to nu.upper.bound, please consider increasing nu.upper.bound, and refitting the model.") + } + } + + return(result) +} \ No newline at end of file diff --git a/R/inla_rspde_spacetime.R b/R/inla_rspde_spacetime.R index 25daf320..2a4e1b82 100644 --- a/R/inla_rspde_spacetime.R +++ b/R/inla_rspde_spacetime.R @@ -39,6 +39,7 @@ #' default values will be used. The `mean` value is also used as starting value for gamma. #' @param prior.precision A precision matrix for \eqn{\log(\kappa), \log(\sigma), \log(\gamma), \rho}. This matrix replaces the precision #' element from `prior.kappa`, `prior.sigma`, `prior.gamma`, and `prior.rho` respectively. For dimension 1 `prior.precision` must be a 4x4 matrix. For dimension 2, \eqn{\rho} is a vector of length 2, so in this case `prior.precision` must be a 5x5 matrix. If `NULL`, a diagonal precision matrix with default values will be used. +#' @param graph_dirichlet For models on metric graphs, use Dirichlet vertex conditions at vertices of degree 1? #' @param bounded_rho Logical. Should `rho` be bounded to ensure the existence, uniqueness, and well-posedness of the solution? Defaults to `TRUE`. #' Note that this bounding is not a strict condition; there may exist values of rho beyond the upper bound that still satisfy these properties. #' When `bounded_rho = TRUE`, the `rspde_lme` models enforce bounded `rho` for consistency. @@ -82,6 +83,7 @@ rspde.spacetime <- function(mesh_space = NULL, prior.rho = NULL, prior.gamma = NULL, prior.precision = NULL, + graph_dirichlet = TRUE, bounded_rho = TRUE, shared_lib = "detect", debug = FALSE, @@ -119,7 +121,8 @@ rspde.spacetime <- function(mesh_space = NULL, rho = prior.rho$mean, alpha = alpha, beta = beta, - bounded_rho = bounded_rho + bounded_rho = bounded_rho, + graph_dirichlet = graph_dirichlet ) default_precision <- 0.1 @@ -234,6 +237,9 @@ rspde.spacetime <- function(mesh_space = NULL, model$mesh <- op$mesh_space } else{ model$mesh <- graph + if(graph_dirichlet) { + class(model$mesh) <- c("metric_graph_dirichlet", class(model$mesh)) + } } model$time_mesh <- op$mesh_time model$rspde_version <- as.character(packageVersion("rSPDE")) @@ -256,6 +262,11 @@ rspde.spacetime <- function(mesh_space = NULL, #' S3method(inlabru::ibm_n, bru_mapper_metric_graph) #' S3method(inlabru::ibm_values, bru_mapper_metric_graph) #' S3method(inlabru::ibm_jacobian, bru_mapper_metric_graph) +#' S3method(inlabru::bru_mapper, metric_graph_dirichlet) +#' S3method(inlabru::ibm_n, bru_mapper_metric_graph_dirichlet) +#' S3method(inlabru::ibm_values, bru_mapper_metric_graph_dirichlet) +#' S3method(inlabru::ibm_jacobian, bru_mapper_metric_graph_dirichlet) + #' } #' @@ -294,3 +305,31 @@ ibm_jacobian.bru_mapper_metric_graph <- function(mapper, input, ...) { } mapper[["mesh"]][["fem_basis"]](input) } + +## Dirichlet mappers +#' @noRd +bru_mapper.metric_graph_dirichlet <- function(mesh, ...) { + mapper <- list(mesh = mesh) + inlabru::bru_mapper_define(mapper, new_class = "bru_mapper_metric_graph_dirichlet") +} + +#' @noRd +ibm_n.bru_mapper_metric_graph_dirichlet <- function(mapper, ...) { + nrow(mapper[["mesh"]][["mesh"]][["VtE"]]) - sum(mapper[["mesh"]]$get_degrees()==1) +} + +#' @noRd +ibm_values.bru_mapper_metric_graph_dirichlet <- function(mapper, ...) { + seq_len(nrow(mapper[["mesh"]][["mesh"]][["VtE"]]) - sum(mapper[["mesh"]]$get_degrees()==1)) +} + +#' @noRd +ibm_jacobian.bru_mapper_metric_graph_dirichlet <- function(mapper, input, ...) { + if (is.null(input)) { + return(Matrix::Matrix(0, 0, inlabru::ibm_n(mapper))) + } + + ind.field <- setdiff(1:nrow(mapper[["mesh"]][["mesh"]][["VtE"]]), + which(mapper[["mesh"]]$get_degrees()==1)) + mapper[["mesh"]][["fem_basis"]](input)[,ind.field] +} diff --git a/R/intrinsic.R b/R/intrinsic.R index c03bd80f..a7ba8c13 100644 --- a/R/intrinsic.R +++ b/R/intrinsic.R @@ -1,6 +1,99 @@ #' @name intrinsic.operators #' @title Covariance-based approximations of intrinsic fields -#' @description `intrinsic.operators` is used for computing a +#' @description `intrinsic.matern.operators` is used for computing a +#' covariance-based rational SPDE approximation of intrinsic +#' fields on \eqn{R^d} defined through the SPDE +#' \deqn{(-\Delta)^{\beta/2} (\tau u) = \mathcal{W}}{(-\Delta)^{\beta/2} (\tau u) = \mathcal{W}} +#' @param tau precision parameter +#' @param beta Smoothness parameter +#' @param G The stiffness matrix of a finite element discretization +#' of the domain of interest. +#' @param C The mass matrix of a finite element discretization of +#' the domain of interest. +#' @param d The dimension of the domain. +#' @param mesh An inla mesh. +#' @param graph An optional `metric_graph` object. Replaces `d`, `C` and `G`. +#' @param loc_mesh locations for the mesh for `d=1`. +#' @param m The order of the rational approximation for the intrinsic part, +#' which needs to be a positive integer. The default value is 2. +#' @param compute_higher_order Logical. Should the higher order finite +#' element matrices be computed? +#' @param return_block_list Logical. For `type = "covariance"`, +#' should the block parts of the precision matrix be returned +#' separately as a list? +#' @param type_rational_approximation Which type of rational +#' approximation should be used? The current types are +#' "chebfun", "brasil" or "chebfunLB". +#' @param fem_mesh_matrices A list containing FEM-related matrices. +#' The list should contain elements c0, g1, g2, g3, etc. +#' @param scaling second lowest eigenvalue of g1 +#' @return `intrinsic.operators` returns an object of +#' class "intrinsicCBrSPDEobj". +#' @export +#' @details The covariance operator +#' \deqn{\tau^{-2}(-\Delta)^{\beta}}{\tau^{-2}(-\Delta)^{\beta}} +#' is approximated based on a rational approximation. The Laplacian is +#' equipped with homogeneous Neumann boundary +#' conditions and a zero-mean constraint is additionally imposed to obtained +#' a non-intrinsic model. +#' @examples +#' if (requireNamespace("RSpectra", quietly = TRUE)) { +#' x <- seq(from = 0, to = 10, length.out = 201) +#' beta <- 1 +#' alpha <- 1 +#' op <- intrinsic.operators(tau = 1, beta = beta, loc_mesh = x, d = 1) +#' # Compute and plot the variogram of the model +#' Sigma <- op$A[,-1] %*% solve(op$Q[-1,-1], t(op$A[,-1])) +#' One <- rep(1, times = ncol(Sigma)) +#' D <- diag(Sigma) +#' Gamma <- 0.5 * (One %*% t(D) + D %*% t(One) - 2 * Sigma) +#' k <- 100 +#' plot(x, Gamma[k, ], type = "l") +#' lines(x, +#' variogram.intrinsic.spde(x[k], x, kappa = 0, alpha = 0, +#' beta = beta, L = 10, d = 1), +#' col = 2, lty = 2 +#' ) +#' } +intrinsic.operators <- function(tau, + alpha, + beta = 1, + G = NULL, + C = NULL, + d = NULL, + mesh = NULL, + graph = NULL, + loc_mesh = NULL, + m_alpha = 2, + m_beta = 2, + compute_higher_order = FALSE, + return_block_list = FALSE, + type_rational_approximation = c( + "chebfun", + "brasil", "chebfunLB" + ), + fem_mesh_matrices = NULL, + scaling = NULL) { + return(intrinsic.matern.operators(tau = tau, + kappa = 0, + alpha = 0, + beta = beta, + G = G, + C = C, + d = d, + mesh = mesh, + graph = graph, + loc_mesh = loc_mesh, + m_beta = m, + compute_higher_order = compute_higher_order, + return_block_list = return_block_list, + type_rational_approximation = type_rational_approximation, + fem_mesh_matrices = fem_mesh_matrices, + scaling = scaling)) +} +#' @name intrinsic.operators.internal +#' @title Covariance-based approximations of intrinsic fields +#' @description `intrinsic.operators.internal` is used for computing a #' covariance-based rational SPDE approximation of intrinsic #' fields on \eqn{R^d} defined through the SPDE #' \deqn{(-\Delta)^{\alpha/2}u = \mathcal{W}}{(-\Delta)^{\alpha/2}u = \mathcal{W}} @@ -25,7 +118,7 @@ #' @param fem_mesh_matrices A list containing FEM-related matrices. #' The list should contain elements c0, g1, g2, g3, etc. #' @param scaling second lowest eigenvalue of g1 -#' @return `intrinsic.operators` returns an object of +#' @return `intrinsic.operators.internal` returns an object of #' class "CBrSPDEobj". This object is a list containing the #' following quantities: #' \item{C}{The mass lumped mass matrix.} @@ -46,7 +139,7 @@ #' so that the equation has a unique solution. This contraint needs to be #' imposed while working with the model later. #' @noRd -intrinsic.operators <- function(C, +intrinsic.operators.internal <- function(C, G, mesh, alpha, @@ -318,7 +411,7 @@ intrinsic.precision <- function(alpha, rspde.order, dim, fem_mesh_matrices, Q <- bdiag(Q, Kpart) - Q <- Q * scaling^alpha + Q <- Q * scaling^(alpha) return(Q) } else { @@ -575,7 +668,7 @@ intrinsic.matern.operators <- function(kappa, return_block_list = TRUE, type_rational_approximation = type_rational_approximation[[1]] ) - op2 <- intrinsic.operators( + op2 <- intrinsic.operators.internal( C = C, G = G, mesh = mesh, alpha = beta, m = m_beta, d = d, compute_higher_order = compute_higher_order, return_block_list = TRUE, fem_mesh_matrices = fem_mesh_matrices, @@ -656,14 +749,14 @@ intrinsic.matern.operators <- function(kappa, } else { alpha_beta <- beta } - op1 <- intrinsic.operators( + op1 <- intrinsic.operators.internal( C = C, G = G, mesh = mesh, alpha = alpha_beta, m = m_beta, d = d, compute_higher_order = compute_higher_order, return_block_list = TRUE, fem_mesh_matrices = fem_mesh_matrices, type_rational_approximation = type_rational_approximation[[1]] ) if (is.list(op1$Q)) { - Q.list1 <- op1$Q*tau^2 + Q.list1 <- lapply(op1$Q, function(x) { x * tau^2}) } else { Q.list1 <- list(op1$Q*tau^2) } @@ -746,7 +839,7 @@ simulate.intrinsicCBrSPDEobj <- function(object, nsim = 1, seed = NULL, X <- matrix(0,n,nsim) for(i in 1:m){ - ind <- (2+n*(i-1)) : n*i + ind <- (2+n*(i-1)) : (n*i) Q <- object$Q[ind,ind] Z <- rnorm((n - 1) * nsim) dim(Z) <- c(n-1, nsim) diff --git a/R/spacetime.R b/R/spacetime.R index 20c85a80..4e5a0f41 100644 --- a/R/spacetime.R +++ b/R/spacetime.R @@ -21,6 +21,7 @@ #' @param beta Integer smoothness parameter beta. #' @param bounded_rho Logical. Specifies whether `rho` should be bounded to ensure the existence, uniqueness, and well-posedness of the solution. Defaults to `TRUE`. #' Note that this bounding is not a strict condition; there may exist values of rho beyond the upper bound that still satisfy these properties. +#' @param graph_dirichlet For models on metric graphs, use Dirichlet vertex conditions at vertices of degree 1? #' When `bounded_rho = TRUE`, the `rspde_lme` models enforce bounded `rho` for consistency. #' If the estimated value of `rho` approaches the upper bound too closely, we recommend refitting the model with `bounded_rho = FALSE`. However, this should be done with caution, as it may lead to instability in some cases, though it can also result in a better model fit. #' The actual bound used for `rho` can be accessed from the `bound_rho` element of the returned object. @@ -52,6 +53,7 @@ spacetime.operators <- function(mesh_space = NULL, rho = NULL, alpha = NULL, beta = NULL, + graph_dirichlet = TRUE, bounded_rho = TRUE) { @@ -132,6 +134,17 @@ spacetime.operators <- function(mesh_space = NULL, G <- graph$mesh$G B <- graph$mesh$B has_graph <- TRUE + if(graph_dirichlet) { + ns <- dim(C)[1] + d1.ind <- which(graph$get_degrees()==1) + ind.field <- setdiff(1:dim(C)[1], d1.ind) + C <- C[ind.field,ind.field] + Ci <- Ci[ind.field,ind.field] + G <- G[ind.field,ind.field] + B <- B[ind.field,ind.field] + } else { + ns <- dim(C)[1] + } } else if (!is.null(mesh_space)) { mesh <- mesh_space d <- fmesher::fm_manifold_dim(mesh) @@ -152,6 +165,7 @@ spacetime.operators <- function(mesh_space = NULL, stop("Only supported for 1d and 2d meshes.") } has_mesh <- TRUE + ns <- dim(C)[1] } else { space_loc <- as.matrix(space_loc) if(min(dim(space_loc))>1) { @@ -164,6 +178,7 @@ spacetime.operators <- function(mesh_space = NULL, B <- fem$B d <- 1 mesh_space <- fm_mesh_1d(space_loc) + ns <- dim(C)[1] } if(alpha+beta= d, where d is the dimension of the spatial domain") @@ -303,9 +318,16 @@ spacetime.operators <- function(mesh_space = NULL, Q <- Q/sigma^2 if (!is.null(graph)) { - make_A <- function(loc, time) { - return(rSPDE.Ast(graph = graph, mesh_time = mesh_time, - obs.s = loc, obs.t = time)) + make_A <- function(loc, time, dirichlet = TRUE) { + if(graph_dirichlet && dirichlet) { + return(rSPDE.Ast(graph = graph, mesh_time = mesh_time, + obs.s = loc, obs.t = time, + ind.field = ind.field)) + } else { + return(rSPDE.Ast(graph = graph, mesh_time = mesh_time, + obs.s = loc, obs.t = time)) + } + } } else { make_A <- function(loc, time) { @@ -348,7 +370,14 @@ spacetime.operators <- function(mesh_space = NULL, plots <- list() for(j in 1:length(t.shift)) { ind <- ((t.ind-t.shift[j]-1)*n+1):((t.ind-t.shift[j])*n) - plots[[j]] <- graph$plot_function(as.vector(tmp[ind]), vertex_size = 0) + if(graph_dirichlet) { + f <- rep(0,dim(C)[1]) + f[ind.field] <- as.vector(tmp[ind]) + } else { + f <- as.vector(tmp[ind]) + } + + plots[[j]] <- graph$plot_function(f, vertex_size = 0) } pt <- ggplot2::ggplot(data.frame(t=mesh_time$loc, y=ct)) + @@ -484,6 +513,17 @@ spacetime.operators <- function(mesh_space = NULL, out$plot_covariances <- plot_covariances out$stationary <- TRUE out$bound_rho <- bound_rho + out$graph_dirichlet <- graph_dirichlet + if(has_graph && graph_dirichlet) { + out$ind.space <- ind.field + ind.st <- c() + for(i in 1:nt) { + ind.st <- c(ind.st, ind.field + (i-1)*ns) + } + out$ind.st <- ind.st + } + out$ns <- ns + out$nt <- nt out$is_bounded_rho <- bounded_rho class(out) <- "spacetimeobj" @@ -633,8 +673,13 @@ simulate.spacetimeobj <- function(object, nsim = 1, LQ <- chol(forceSymmetric(object$Q)) X <- solve(LQ, Z) - - return(as.matrix(X)) + if(object$has_graph && object$graph_dirichlet) { + Xout <- matrix(0,object$ns*object$nt, nsim) + Xout[object$ind.st,] <- as.matrix(X) + return(Xout) + } else { + return(as.matrix(X)) + } } @@ -1030,6 +1075,7 @@ aux2_lme_spacetime.loglike <- function(object, y, X_cov, repl, A_list, sigma_e, #' @param graph MetricGraph object for models on metric graphs #' @param obs.s spatial locations of observations #' @param obs.t time points for observations +#' @param ind.field indices for field (if Dirichlet conditions are used) #' #' @return Observation matrix linking observation locations to mesh nodes #' @export @@ -1048,7 +1094,8 @@ rSPDE.Ast <- function(mesh_space = NULL, time_loc = NULL, graph = NULL, obs.s = NULL, - obs.t = NULL) { + obs.t = NULL, + ind.field = NULL) { if ((!is.null(mesh_space) && !is.null(graph)) || (!is.null(mesh_space) && !is.null(space_loc)) || (!is.null(graph) && !is.null(space_loc))){ stop("You should provide only one of mesh_space, space_loc or graph.") @@ -1080,8 +1127,14 @@ rSPDE.Ast <- function(mesh_space = NULL, At <- rSPDE.A1d(time, obs.t) if(!is.null(graph)) { As <- graph$fem_basis(obs.s) + if(!is.null(ind.field)) { + As <- As[,ind.field] + } } else if (!is.null(mesh_space)){ As <- fm_basis(mesh_space, obs.s) + if(!is.null(ind.field)) { + As <- As[,ind.field] + } } else { space_loc <- as.matrix(space_loc) if(min(dim(space_loc))>1) { @@ -1089,6 +1142,9 @@ rSPDE.Ast <- function(mesh_space = NULL, } mesh_space <- fm_mesh_1d(space_loc) As <- fm_basis(mesh_space, obs.s) + if(!is.null(ind.field)) { + As <- As[,ind.field] + } } As.bar <- kronecker(matrix(1,nrow=1, ncol=length(time)),As) diff --git a/R/util.R b/R/util.R index 58bc4524..d499acca 100644 --- a/R/util.R +++ b/R/util.R @@ -507,7 +507,7 @@ cut_decimals <- function(nu) { #' @noRd check_class_inla_rspde <- function(model) { - if (!inherits(model, c("inla_rspde", "inla_rspde_matern1d"))) { + if (!inherits(model, c("inla_rspde", "inla_rspde_matern1d", "inla_rspde_fintrinsic"))) { stop("You should provide a rSPDE model!") } } @@ -2173,41 +2173,56 @@ transform_parameters_anisotropic <- function(theta, nu_upper_bound = NULL) { #' @noRd - rspde_check_cgeneric_symbol <- function(model) { - # Ensure the required fields exist in the model object - if (!"f" %in% names(model) || !"cgeneric" %in% names(model$f) || - !"shlib" %in% names(model$f$cgeneric) || !"model" %in% names(model$f$cgeneric)) { - stop("There was a problem with the model creation.") - } - - # Extract the shared library path and the symbol name - shlib <- model$f$cgeneric$shlib - symbol <- model$f$cgeneric$model - - # Check if the shared library exists - if (!file.exists(shlib)) { - stop(paste("The shared library", shlib, "does not exist.")) - } - - # Use the `dyn.load` and `is.loaded` functions to check for the symbol - tryCatch({ - dyn.load(shlib) # Load the shared library - if (is.loaded(symbol)) { - dyn.unload(shlib) # Unload if the symbol is available - return(invisible(TRUE)) # Return silently - } else { - warning(paste0("The symbol '", symbol, "' is not available in the shared library. Please install the latest testing version of INLA. + # Ensure the required fields exist in the model object + if (!"f" %in% names(model) || !"cgeneric" %in% names(model$f) || + !"shlib" %in% names(model$f$cgeneric) || !"model" %in% names(model$f$cgeneric)) { + stop("There was a problem with the model creation.") + } + + # Extract the shared library path and the symbol name + shlib <- model$f$cgeneric$shlib + symbol <- model$f$cgeneric$model + + # Check if the shared library exists + if (!file.exists(shlib)) { + stop(paste("The shared library", shlib, "does not exist.")) + } + + # Get R_HOME library path + r_lib_path <- file.path(R.home("lib")) + + # Add R_HOME to library path if it is not there + current_lib_path <- Sys.getenv("LD_LIBRARY_PATH") + if (!grepl(r_lib_path, current_lib_path)) { + if (current_lib_path == "") { + Sys.setenv(LD_LIBRARY_PATH = r_lib_path) + } else { + Sys.setenv(LD_LIBRARY_PATH = paste(current_lib_path, r_lib_path, sep = ":")) + } + } + + # Use the `dyn.load` and `is.loaded` functions to check for the symbol + tryCatch({ + dyn.load(shlib) # Load the shared library + if (is.loaded(symbol)) { + dyn.unload(shlib) # Unload if the symbol is available + return(invisible(TRUE)) # Return silently + } else { + warning(paste0("The symbol '", symbol, "' is not available in the shared library. Please install the latest testing version of INLA. If the problem persists after installing the latest testing version of INLA, please open an issue at https://github.com/davidbolin/rSPDE/issues, requesting that this model be added to INLA.")) - } - dyn.unload(shlib) # Ensure the library is unloaded - }, error = function(e) { - warning(paste0("Error while loading the shared library or checking the symbol: ", e$message, - ". Please install the latest testing version of INLA. If the problem persists after installing the + } + dyn.unload(shlib) # Ensure the library is unloaded + }, error = function(e) { + warning(paste0("Error while loading the shared library or checking the symbol: ", e$message, + ". Please install the latest testing version of INLA. If the problem persists after installing the latest testing version of INLA, please open an issue at https://github.com/davidbolin/rSPDE/issues, requesting that this model be added to INLA.")) - }) + }) + + # Restore original LD_LIBRARY_PATH + Sys.setenv(LD_LIBRARY_PATH = current_lib_path) } #' @noRd diff --git a/man/bru_get_mapper.inla_rspde.Rd b/man/bru_get_mapper.inla_rspde.Rd index 1ce70685..75728dc1 100644 --- a/man/bru_get_mapper.inla_rspde.Rd +++ b/man/bru_get_mapper.inla_rspde.Rd @@ -7,13 +7,13 @@ \alias{ibm_jacobian.bru_mapper_inla_rspde} \title{rSPDE inlabru mapper} \usage{ -\method{bru_get_mapper}{inla_rspde}(model, ...) +bru_get_mapper.inla_rspde(model, ...) -\method{ibm_n}{bru_mapper_inla_rspde}(mapper, ...) +ibm_n.bru_mapper_inla_rspde(mapper, ...) -\method{ibm_values}{bru_mapper_inla_rspde}(mapper, ...) +ibm_values.bru_mapper_inla_rspde(mapper, ...) -\method{ibm_jacobian}{bru_mapper_inla_rspde}(mapper, input, ...) +ibm_jacobian.bru_mapper_inla_rspde(mapper, input, ...) } \arguments{ \item{model}{An \code{inla_rspde} object for which to construct or extract a mapper} diff --git a/man/bru_get_mapper.inla_rspde_anisotropic2d.Rd b/man/bru_get_mapper.inla_rspde_anisotropic2d.Rd index f8c4d760..282769ed 100644 --- a/man/bru_get_mapper.inla_rspde_anisotropic2d.Rd +++ b/man/bru_get_mapper.inla_rspde_anisotropic2d.Rd @@ -4,7 +4,7 @@ \alias{bru_get_mapper.inla_rspde_anisotropic2d} \title{rSPDE anisotropic inlabru mapper} \usage{ -\method{bru_get_mapper}{inla_rspde_anisotropic2d}(model, ...) +bru_get_mapper.inla_rspde_anisotropic2d(model, ...) } \arguments{ \item{model}{An \code{inla_rspde_anisotropic2d} object for which to construct or extract a mapper} diff --git a/man/bru_get_mapper.inla_rspde_matern1d.Rd b/man/bru_get_mapper.inla_rspde_matern1d.Rd index 33d67587..cc2b8187 100644 --- a/man/bru_get_mapper.inla_rspde_matern1d.Rd +++ b/man/bru_get_mapper.inla_rspde_matern1d.Rd @@ -7,13 +7,13 @@ \alias{ibm_jacobian.bru_mapper_inla_rspde_matern1d} \title{rSPDE stationary inlabru mapper} \usage{ -\method{bru_get_mapper}{inla_rspde_matern1d}(model, ...) +bru_get_mapper.inla_rspde_matern1d(model, ...) -\method{ibm_n}{bru_mapper_inla_rspde_matern1d}(mapper, ...) +ibm_n.bru_mapper_inla_rspde_matern1d(mapper, ...) -\method{ibm_values}{bru_mapper_inla_rspde_matern1d}(mapper, ...) +ibm_values.bru_mapper_inla_rspde_matern1d(mapper, ...) -\method{ibm_jacobian}{bru_mapper_inla_rspde_matern1d}(mapper, input, ...) +ibm_jacobian.bru_mapper_inla_rspde_matern1d(mapper, input, ...) } \arguments{ \item{model}{An \code{inla_rspde_matern1d} object for which to construct or extract a mapper} diff --git a/man/bru_get_mapper.inla_rspde_spacetime.Rd b/man/bru_get_mapper.inla_rspde_spacetime.Rd index c71f532d..298a131a 100644 --- a/man/bru_get_mapper.inla_rspde_spacetime.Rd +++ b/man/bru_get_mapper.inla_rspde_spacetime.Rd @@ -4,7 +4,7 @@ \alias{bru_get_mapper.inla_rspde_spacetime} \title{rSPDE space time inlabru mapper} \usage{ -\method{bru_get_mapper}{inla_rspde_spacetime}(model, ...) +bru_get_mapper.inla_rspde_spacetime(model, ...) } \arguments{ \item{model}{An \code{inla_rspde_spacetime} object for which to construct or extract a mapper} diff --git a/man/intrinsic.operators.Rd b/man/intrinsic.operators.Rd new file mode 100644 index 00000000..fd11c729 --- /dev/null +++ b/man/intrinsic.operators.Rd @@ -0,0 +1,101 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/intrinsic.R +\name{intrinsic.operators} +\alias{intrinsic.operators} +\title{Covariance-based approximations of intrinsic fields} +\usage{ +intrinsic.operators( + tau, + alpha, + beta = 1, + G = NULL, + C = NULL, + d = NULL, + mesh = NULL, + graph = NULL, + loc_mesh = NULL, + m_alpha = 2, + m_beta = 2, + compute_higher_order = FALSE, + return_block_list = FALSE, + type_rational_approximation = c("chebfun", "brasil", "chebfunLB"), + fem_mesh_matrices = NULL, + scaling = NULL +) +} +\arguments{ +\item{tau}{precision parameter} + +\item{beta}{Smoothness parameter} + +\item{G}{The stiffness matrix of a finite element discretization +of the domain of interest.} + +\item{C}{The mass matrix of a finite element discretization of +the domain of interest.} + +\item{d}{The dimension of the domain.} + +\item{mesh}{An inla mesh.} + +\item{graph}{An optional \code{metric_graph} object. Replaces \code{d}, \code{C} and \code{G}.} + +\item{loc_mesh}{locations for the mesh for \code{d=1}.} + +\item{compute_higher_order}{Logical. Should the higher order finite +element matrices be computed?} + +\item{return_block_list}{Logical. For \code{type = "covariance"}, +should the block parts of the precision matrix be returned +separately as a list?} + +\item{type_rational_approximation}{Which type of rational +approximation should be used? The current types are +"chebfun", "brasil" or "chebfunLB".} + +\item{fem_mesh_matrices}{A list containing FEM-related matrices. +The list should contain elements c0, g1, g2, g3, etc.} + +\item{scaling}{second lowest eigenvalue of g1} + +\item{m}{The order of the rational approximation for the intrinsic part, +which needs to be a positive integer. The default value is 2.} +} +\value{ +\code{intrinsic.operators} returns an object of +class "intrinsicCBrSPDEobj". +} +\description{ +\code{intrinsic.matern.operators} is used for computing a +covariance-based rational SPDE approximation of intrinsic +fields on \eqn{R^d} defined through the SPDE +\deqn{(-\Delta)^{\beta/2} (\tau u) = \mathcal{W}}{(-\Delta)^{\beta/2} (\tau u) = \mathcal{W}} +} +\details{ +The covariance operator +\deqn{\tau^{-2}(-\Delta)^{\beta}}{\tau^{-2}(-\Delta)^{\beta}} +is approximated based on a rational approximation. The Laplacian is +equipped with homogeneous Neumann boundary +conditions and a zero-mean constraint is additionally imposed to obtained +a non-intrinsic model. +} +\examples{ +if (requireNamespace("RSpectra", quietly = TRUE)) { + x <- seq(from = 0, to = 10, length.out = 201) + beta <- 1 + alpha <- 1 + op <- intrinsic.operators(tau = 1, beta = beta, loc_mesh = x, d = 1) + # Compute and plot the variogram of the model + Sigma <- op$A[,-1] \%*\% solve(op$Q[-1,-1], t(op$A[,-1])) + One <- rep(1, times = ncol(Sigma)) + D <- diag(Sigma) + Gamma <- 0.5 * (One \%*\% t(D) + D \%*\% t(One) - 2 * Sigma) + k <- 100 + plot(x, Gamma[k, ], type = "l") + lines(x, + variogram.intrinsic.spde(x[k], x, kappa = 0, alpha = 0, + beta = beta, L = 10, d = 1), + col = 2, lty = 2 + ) +} +} diff --git a/man/rSPDE.Ast.Rd b/man/rSPDE.Ast.Rd index afdbbfe2..703d31f6 100644 --- a/man/rSPDE.Ast.Rd +++ b/man/rSPDE.Ast.Rd @@ -11,7 +11,8 @@ rSPDE.Ast( time_loc = NULL, graph = NULL, obs.s = NULL, - obs.t = NULL + obs.t = NULL, + ind.field = NULL ) } \arguments{ @@ -28,6 +29,8 @@ rSPDE.Ast( \item{obs.s}{spatial locations of observations} \item{obs.t}{time points for observations} + +\item{ind.field}{indices for field (if Dirichlet conditions are used)} } \value{ Observation matrix linking observation locations to mesh nodes diff --git a/man/rspde.intrinsic.Rd b/man/rspde.intrinsic.Rd new file mode 100644 index 00000000..9fcc5f2d --- /dev/null +++ b/man/rspde.intrinsic.Rd @@ -0,0 +1,84 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/inla_rspde_intrinsic.R +\name{rspde.intrinsic} +\alias{rspde.intrinsic} +\title{Rational approximations of fractional intrinsic fields} +\usage{ +rspde.intrinsic( + mesh, + nu = NULL, + nu.upper.bound = 2, + mean.correction = FALSE, + rspde.order = 1, + prior.tau = NULL, + prior.nu = NULL, + prior.nu.dist = "lognormal", + nu.prec.inc = 0.01, + diagonal = 0, + type.rational.approx = "chebfun", + shared_lib = "detect", + debug = FALSE, + ... +) +} +\arguments{ +\item{mesh}{Spatial mesh for the FEM approximation.} + +\item{nu}{If nu is set to a parameter, nu will be kept fixed and will not +be estimated. If nu is \code{NULL}, it will be estimated.} + +\item{nu.upper.bound}{Upper bound for the smoothness parameter \eqn{\nu}. If \code{NULL}, it will be set to 2.} + +\item{mean.correction}{Add mean correction for extreme value models?} + +\item{rspde.order}{The order of the covariance-based rational SPDE approach. The default order is 1.} + +\item{prior.tau}{A list specifying the prior for the variance parameter \eqn{\tau}. +This list may contain two elements: \code{mean} and/or \code{precision}, both of which must +be numeric scalars.} + +\item{prior.nu}{a list containing the elements \code{mean} and \code{prec} +for beta distribution, or \code{loglocation} and \code{logscale} for a +truncated lognormal distribution. \code{loglocation} stands for +the location parameter of the truncated lognormal distribution in the log +scale. \code{prec} stands for the precision of a beta distribution. +\code{logscale} stands for the scale of the truncated lognormal +distribution on the log scale. Check details below.} + +\item{prior.nu.dist}{The distribution of the smoothness parameter. +The current options are "beta" or "lognormal". The default is "lognormal".} + +\item{nu.prec.inc}{Amount to increase the precision in the beta prior +distribution. Check details below.} + +\item{diagonal}{Number added to diagonal of Q for increased stability.} + +\item{type.rational.approx}{Which type of rational approximation +should be used? The current types are "chebfun", "brasil" or "chebfunLB".} + +\item{shared_lib}{String specifying which shared library to use for the Cgeneric +implementation. Options are "detect", "INLA", or "rSPDE". You may also specify the +direct path to a .so (or .dll) file.} + +\item{debug}{Logical value indicating whether to enable INLA debug mode.} + +\item{...}{Additional arguments passed internally for configuration purposes.} +} +\value{ +An object of class \code{inla_rspde_intrinsic} representing the FEM approximation of +the intrinsic Gaussian random field. +} +\description{ +\code{rspde.intrinsic} computes a Finite Element Method (FEM) approximation of a +Gaussian random field defined as the solution to the stochastic partial +differential equation (SPDE): +\deqn{(-\Delta)^{(\nu+d/2)/2}\tau u = W}. +} +\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)) +model <- rspde.intrinsic(mesh = mesh_2d) + +} diff --git a/man/rspde.spacetime.Rd b/man/rspde.spacetime.Rd index 35248b18..f5292162 100644 --- a/man/rspde.spacetime.Rd +++ b/man/rspde.spacetime.Rd @@ -17,6 +17,7 @@ rspde.spacetime( prior.rho = NULL, prior.gamma = NULL, prior.precision = NULL, + graph_dirichlet = TRUE, bounded_rho = TRUE, shared_lib = "detect", debug = FALSE, @@ -66,6 +67,8 @@ default values will be used. The \code{mean} value is also used as starting valu \item{prior.precision}{A precision matrix for \eqn{\log(\kappa), \log(\sigma), \log(\gamma), \rho}. This matrix replaces the precision element from \code{prior.kappa}, \code{prior.sigma}, \code{prior.gamma}, and \code{prior.rho} respectively. For dimension 1 \code{prior.precision} must be a 4x4 matrix. For dimension 2, \eqn{\rho} is a vector of length 2, so in this case \code{prior.precision} must be a 5x5 matrix. If \code{NULL}, a diagonal precision matrix with default values will be used.} +\item{graph_dirichlet}{For models on metric graphs, use Dirichlet vertex conditions at vertices of degree 1?} + \item{bounded_rho}{Logical. Should \code{rho} be bounded to ensure the existence, uniqueness, and well-posedness of the solution? Defaults to \code{TRUE}. Note that this bounding is not a strict condition; there may exist values of rho beyond the upper bound that still satisfy these properties. When \code{bounded_rho = TRUE}, the \code{rspde_lme} models enforce bounded \code{rho} for consistency. diff --git a/man/spacetime.operators.Rd b/man/spacetime.operators.Rd index b2aa207d..e78ed97c 100644 --- a/man/spacetime.operators.Rd +++ b/man/spacetime.operators.Rd @@ -16,6 +16,7 @@ spacetime.operators( rho = NULL, alpha = NULL, beta = NULL, + graph_dirichlet = TRUE, bounded_rho = TRUE ) } @@ -44,11 +45,13 @@ one-dimensional spatial domains, a vector with two number on 2d domains.} \item{beta}{Integer smoothness parameter beta.} -\item{bounded_rho}{Logical. Specifies whether \code{rho} should be bounded to ensure the existence, uniqueness, and well-posedness of the solution. Defaults to \code{TRUE}. -Note that this bounding is not a strict condition; there may exist values of rho beyond the upper bound that still satisfy these properties. +\item{graph_dirichlet}{For models on metric graphs, use Dirichlet vertex conditions at vertices of degree 1? When \code{bounded_rho = TRUE}, the \code{rspde_lme} models enforce bounded \code{rho} for consistency. If the estimated value of \code{rho} approaches the upper bound too closely, we recommend refitting the model with \code{bounded_rho = FALSE}. However, this should be done with caution, as it may lead to instability in some cases, though it can also result in a better model fit. The actual bound used for \code{rho} can be accessed from the \code{bound_rho} element of the returned object.} + +\item{bounded_rho}{Logical. Specifies whether \code{rho} should be bounded to ensure the existence, uniqueness, and well-posedness of the solution. Defaults to \code{TRUE}. +Note that this bounding is not a strict condition; there may exist values of rho beyond the upper bound that still satisfy these properties.} } \value{ An object of type spacetimeobj. diff --git a/src/cgeneric_aux_anisotropic.cpp b/src/cgeneric_aux_anisotropic.cpp index 9a403000..346dca8c 100644 --- a/src/cgeneric_aux_anisotropic.cpp +++ b/src/cgeneric_aux_anisotropic.cpp @@ -91,7 +91,7 @@ Eigen::SparseMatrix anisotropic_precision( } // Apply scaling factor and tau - Q *= std::pow(scale_factor, 2 * alpha) * std::pow(tau, 2); + Q *= std::pow(scale_factor, alpha) * std::pow(tau, 2); return Q; } diff --git a/src/cgeneric_aux_intrinsic.cpp b/src/cgeneric_aux_intrinsic.cpp index 45c9ebd9..5c04077b 100644 --- a/src/cgeneric_aux_intrinsic.cpp +++ b/src/cgeneric_aux_intrinsic.cpp @@ -11,6 +11,45 @@ #include "cgeneric_cpp.h" +Eigen::VectorXd intrinsic_mean(Eigen::SparseMatrix Lq) { + + Eigen::SparseMatrix S = Lq.selfadjointView(); + + int n = Lq.rows(); + + for (int i = n - 1; i >= 0; --i) { + Eigen::SparseMatrix::ReverseInnerIterator Si(S, i); + for (Eigen::SparseMatrix::ReverseInnerIterator ij(Lq, i); ij; --ij) { + + Eigen::SparseMatrix::ReverseInnerIterator iL(Lq, i); + Eigen::SparseMatrix::ReverseInnerIterator iS(S, ij.row()); + Si.valueRef() = 0.0; + while (iL.row() > i) { + while (iS && (iL.row() < iS.row())) { + --iS; + } + if (iS && (iL.row() == iS.row())) { + Si.valueRef() -= iL.value() * iS.value(); + --iS; + } + --iL; + } + if (i == ij.row()) { + Si.valueRef() += 1 / iL.value(); + Si.valueRef() /= iL.value(); + } else { + Si.valueRef() /= iL.value(); + while (iS.row() > i) { + --iS; + } + iS.valueRef() = Si.value(); + } + --Si; + } + } + + return S.diagonal(); +} extern "C" void compute_Q_fintrinsic( double tau, double nu, @@ -22,7 +61,9 @@ extern "C" void compute_Q_fintrinsic( const inla_cgeneric_mat_tp *rational_table, int est_nu, int d, int compute_Q, int compute_mean, int compute_logdet, - double*ld_out, double *mean_out); + double*ld_out, double *mean_out, + double scaling, + const inla_cgeneric_smat_tp *D); void compute_Q_fintrinsic( double tau, double nu, @@ -34,28 +75,31 @@ void compute_Q_fintrinsic( double tau, double nu, const inla_cgeneric_mat_tp *rational_table, int est_nu, int d, int compute_Q, int compute_mean, int compute_logdet, - double*ld_out, double *mean_out) { + double*ld_out, double *mean_out, + double scaling, + const inla_cgeneric_smat_tp *D) { double alpha = nu + d/2; - - std::cout << "set matrices" << std::endl; + //std::cout << "tau = " << tau << ", nu = " << nu << std::endl; + // Assembling the FEM matrices - SparseMatrixColMajor C_eigen = convertInlaToEigen(C); - SparseMatrixColMajor Ci_eigen = convertInlaToEigen(Ci); - SparseMatrixColMajor G_eigen = convertInlaToEigen(G); + Eigen::SparseMatrix C_eigen = convertInlaToEigen(C); + Eigen::SparseMatrix Ci_eigen = convertInlaToEigen(Ci); + Eigen::SparseMatrix G_eigen = convertInlaToEigen(G); + Eigen::SparseMatrix D_eigen = convertInlaToEigen(D); - std::cout << "set table" << std::endl; // Assembling the rational table Eigen::MatrixXd rational_table_eigen = inlaToEigenMatrix(rational_table); int m_alpha = static_cast(std::floor(alpha)); - SparseMatrixColMajor L = G_eigen; - SparseMatrixColMajor CiL = Ci_eigen * L; + Eigen::SparseMatrix L = G_eigen / scaling; + Eigen::SparseMatrix CiL = Ci_eigen * L; - SparseMatrixColMajor Q; - std::cout << "alpha = " << alpha << std::endl; - if (std::floor(alpha) == alpha && est_nu == 0) { // Check if alpha is an integer + Eigen::SparseMatrix Q; + + bool int_alpha = std::floor(alpha) == alpha && est_nu == 0; // Check if alpha is an integer + if (int_alpha) { Q = L; if (alpha > 1) { for (int k = 1; k < alpha; ++k) { @@ -64,33 +108,23 @@ void compute_Q_fintrinsic( double tau, double nu, } Q *= tau * tau; } else if (rspde_order > 0) { - SparseMatrixColMajor Q_graph_eigen = convertInlaToEigen(Q_graph); - SparseMatrixColMajor Q_graph_transpose = Q_graph_eigen.transpose(); + Eigen::SparseMatrix Q_graph_eigen = convertInlaToEigen(Q_graph); + Eigen::SparseMatrix Q_graph_transpose = Q_graph_eigen.transpose(); Q_graph_eigen = Q_graph_eigen + Q_graph_transpose; - std::cout << "create precision" << std::endl; - std::cout << L.rows() << " " << L.cols() << std::endl; - std::cout << C_eigen.rows() << " " << C_eigen.cols() << std::endl; - std::cout << Ci_eigen.rows() << " " << Ci_eigen.cols() << std::endl; - std::cout << CiL.rows() << " " << CiL.cols() << std::endl; - std::cout << alpha << " " << m_alpha << " " << rspde_order << std::endl; Q = anisotropic_precision(L, tau, C_eigen, Ci_eigen, CiL, alpha, m_alpha, - rspde_order, 1.0, rational_table_eigen); - std::cout << "add graph " << Q_graph_eigen.rows() << " " << Q_graph_eigen.cols() << std::endl; - std::cout << "add graph " << Q.rows() << " " << Q.cols() << std::endl; + rspde_order, scaling, rational_table_eigen); Q = Q + 0 * Q_graph_eigen; } else { throw std::invalid_argument("rspde_order > 0 required"); } - - if(compute_Q== 1) { - std::cout << "save Q" << std::endl; + if(compute_Q == 1) { + Eigen::SparseMatrix Qadj = Q + tau * tau * D_eigen; // add diagonal correction // Extract the values from Q into the result array, using only lower triangular part - Eigen::SparseMatrix Q_triang = Q.triangularView(); - + Eigen::SparseMatrix Q_triang = Qadj.triangularView(); int count = 0; for (int k = 0; k < Q_triang.outerSize(); ++k) { - for (Eigen::SparseMatrix::InnerIterator it(Q_triang, k); it; ++it) { + for (Eigen::SparseMatrix::InnerIterator it(Q_triang, k); it; ++it) { ret[count] = it.value(); count++; } @@ -98,66 +132,53 @@ void compute_Q_fintrinsic( double tau, double nu, } if(compute_mean || compute_logdet){ - std::cout << "compute mean" << std::endl; - int size = Q.rows(); - // Get submatrix Q[-n,-,n] - Eigen::SparseMatrix Qsub(size-1,size-1); - Qsub = Q.topLeftCorner(size-1,size-1); - - Eigen::SimplicialLLT > R; - R.analyzePattern(Qsub); - R.factorize(Qsub); - Eigen::SparseMatrix L = R.matrixL(); - if(compute_logdet){ - // constant = log(|Q|*/(2pi)^((d-1)/2), |Q|* = d|Qsub| - //log const = log(|Q|*) - (d-1)log(2pi)/2 - // = log(2d) + log(|R|*) - (d-1)log(2pi)/2 - double ldet = L.diagonal().array().log().sum(); - ld_out[0] = ldet + log(2*size) - (size- 1) * log(2 * M_PI) / 2; - } - - - if(compute_mean){ - Eigen::SparseMatrix S = L.selfadjointView(); - - int n = L.rows(); + Eigen::SimplicialLLT > R; + Eigen::SparseMatrix Lq; + Eigen::VectorXd Qidiag_perm; + Eigen::VectorXd Qidiag; + int n = C_eigen.rows(); + if(int_alpha) { + // Get submatrix Q[-n,-,n] + Eigen::SparseMatrix Qsub = Q.topLeftCorner(n-1,n-1); + R.analyzePattern(Qsub); + R.factorize(Qsub); - for (int i = n - 1; i >= 0; --i) { - Eigen::SparseMatrix::ReverseInnerIterator Si(S, i); - for (Eigen::SparseMatrix::ReverseInnerIterator ij(L, i); ij; --ij) { - - Eigen::SparseMatrix::ReverseInnerIterator iL(L, i); - Eigen::SparseMatrix::ReverseInnerIterator iS(S, ij.row()); - Si.valueRef() = 0.0; - while (iL.row() > i) { - while (iS && (iL.row() < iS.row())) { - --iS; - } - if (iS && (iL.row() == iS.row())) { - Si.valueRef() -= iL.value() * iS.value(); - --iS; - } - --iL; - } - if (i == ij.row()) { - Si.valueRef() += 1 / iL.value(); - Si.valueRef() /= iL.value(); - } else { - Si.valueRef() /= iL.value(); - while (iS.row() > i) { - --iS; - } - iS.valueRef() = Si.value(); - } - --Si; - } + Lq = R.matrixL(); + if(compute_logdet){ + // constant = log(|Q|*/(2pi)^((d-1)/2), |Q|* = d|Qsub| + //log const = log(|Q|*) - (d-1)log(2pi)/2 + // = log(2d) + log(|R|*) - (d-1)log(2pi)/2 + double ldet = Lq.diagonal().array().log().sum(); + ld_out[0] = ldet + log(2*n) - (n - 1) * log(2 * M_PI) / 2; } - - Eigen::VectorXd Qidiag_perm = S.diagonal(); - Eigen::VectorXd Qidiag = R.permutationPinv() * Qidiag_perm; - for(int i = 0; i < n; i++){ - mean_out[i] = Qidiag[i]; + if(compute_mean){ + Qidiag_perm = intrinsic_mean(Lq); + Qidiag = R.permutationPinv() * Qidiag_perm; + for(int i = 0; i < n; i++){ + mean_out[i] = Qidiag[i]; + } } + } else { + int start; + ld_out[0] = (rspde_order + 1)*log(2*n) - (rspde_order + 1)* (n - rspde_order) * log(2 * M_PI) / 2; + for(int k = 0; k < rspde_order + 1; k ++) { + start = k*n; + Eigen::SparseMatrix Qsub = Q.block(start,start,n-1,n-1); + R.analyzePattern(Qsub); + R.factorize(Qsub); + Lq = R.matrixL(); + if(compute_logdet){ + ld_out[0] += Lq.diagonal().array().log().sum(); + } + + if(compute_mean){ + Qidiag_perm = intrinsic_mean(Lq); + Qidiag = R.permutationPinv() * Qidiag_perm; + for(int i = 0; i < n; i++){ + mean_out[start + i] = Qidiag[i]; + } + } + } } } } @@ -244,10 +265,10 @@ void compute_Q_intrinsic(int size, Eigen::SparseMatrix Qsub(size-1,size-1); Qsub = Q.topLeftCorner(size-1,size-1); - Eigen::SimplicialLLT > R; + Eigen::SimplicialLLT > R; R.analyzePattern(Qsub); R.factorize(Qsub); - Eigen::SparseMatrix L = R.matrixL(); + Eigen::SparseMatrix L = R.matrixL(); if(compute_logdet){ // constant = log(|Q|*/(2pi)^((d-1)/2), |Q|* = d|Qsub| //log const = log(|Q|*) - (d-1)log(2pi)/2 @@ -258,7 +279,7 @@ void compute_Q_intrinsic(int size, if(compute_mean){ - Eigen::SparseMatrix S = L.selfadjointView(); + Eigen::SparseMatrix S = L.selfadjointView(); int n = L.rows(); diff --git a/src/cgeneric_defs.h b/src/cgeneric_defs.h index 497c8f86..fadadab5 100644 --- a/src/cgeneric_defs.h +++ b/src/cgeneric_defs.h @@ -137,7 +137,9 @@ void compute_Q_fintrinsic(double tau, double nu, const inla_cgeneric_mat_tp *rational_table, int est_nu, int d, int compute_Q, int compute_mean, int compute_logdet, - double*ld_out, double *mean_out); + double*ld_out, double *mean_out, + double scaling, + const inla_cgeneric_smat_tp *D); #ifdef __cplusplus } diff --git a/src/cgeneric_rspde_fintrinsic.c b/src/cgeneric_rspde_fintrinsic.c new file mode 100644 index 00000000..277b1eda --- /dev/null +++ b/src/cgeneric_rspde_fintrinsic.c @@ -0,0 +1,333 @@ +#include "cgeneric_defs.h" +#include "stdio.h" + + + +typedef struct +{ + double *Q; + double *mu; + double *lconst; + double *theta; +} +my_cache_tp; + +double *inla_cgeneric_rspde_fintrinsic_model(inla_cgeneric_cmd_tp cmd, double *theta, + inla_cgeneric_data_tp *data) { + + double *ret = NULL; + double *mu_store = NULL; + double *const_store = NULL; + double*Q_store = NULL; + int k, i; + double lnu, ltau, tau; + char *prior_nu_dist; + int est_nu = 0; + + // Retrieve parameter values from `data` + assert(!strcasecmp(data->ints[0]->name, "n")); + int N = data->ints[0]->ints[0]; + assert(N > 0); + + // Basic parameter assertions and retrievals + assert(!strcasecmp(data->ints[2]->name, "est_nu")); + est_nu = data->ints[2]->ints[0]; + + assert(!strcasecmp(data->ints[3]->name, "rspde_order")); + int rspde_order = data->ints[3]->ints[0]; + + assert(!strcasecmp(data->ints[4]->name, "dim")); + int d = data->ints[4]->ints[0]; + + assert(!strcasecmp(data->ints[5]->name, "mean_correction")); + int mean_correction = data->ints[5]->ints[0]; + + // Retrieve prior means for each parameter + assert(!strcasecmp(data->doubles[0]->name, "prior.tau.mean")); + double prior_tau_mean = data->doubles[0]->doubles[0]; + + assert(!strcasecmp(data->doubles[1]->name, "prior.tau.precision")); + double prior_tau_prec = data->doubles[1]->doubles[0]; + + assert(!strcasecmp(data->doubles[2]->name, "start_nu")); + double start_nu = data->doubles[2]->doubles[0]; + + assert(!strcasecmp(data->doubles[3]->name, "nu")); + double nu = data->doubles[3]->doubles[0]; + + assert(!strcasecmp(data->doubles[4]->name, "prior.nu.loglocation")); + double prior_nu_loglocation = data->doubles[4]->doubles[0]; + + assert(!strcasecmp(data->doubles[5]->name, "prior.nu.mean")); + double prior_nu_mean = data->doubles[5]->doubles[0]; + + assert(!strcasecmp(data->doubles[6]->name, "prior.nu.prec")); + double prior_nu_prec = data->doubles[6]->doubles[0]; + + assert(!strcasecmp(data->doubles[7]->name, "prior.nu.logscale")); + double prior_nu_logscale = data->doubles[7]->doubles[0]; + + assert(!strcasecmp(data->doubles[8]->name, "nu_upper_bound")); + double nu_upper_bound = data->doubles[8]->doubles[0]; + + assert(!strcasecmp(data->doubles[9]->name, "scaling")); + double scaling = data->doubles[9]->doubles[0]; + + assert(!strcasecmp(data->chars[2]->name, "prior.nu.dist")); + prior_nu_dist = &data->chars[2]->chars[0]; + + // Retrieve sparse matrix Q + assert(!strcasecmp(data->smats[0]->name, "Q")); + inla_cgeneric_smat_tp *Q = data->smats[0]; + int M = Q->n; + + assert(!strcasecmp(data->smats[1]->name, "C")); + inla_cgeneric_smat_tp *C = data->smats[1]; + + assert(!strcasecmp(data->smats[2]->name, "Ci")); + inla_cgeneric_smat_tp *Ci = data->smats[2]; + + assert(!strcasecmp(data->smats[3]->name, "G")); + inla_cgeneric_smat_tp *G = data->smats[3]; + + assert(!strcasecmp(data->smats[4]->name, "D")); + inla_cgeneric_smat_tp *D = data->smats[4]; + + assert(!strcasecmp(data->mats[0]->name, "rational_table")); + inla_cgeneric_mat_tp *rational_table = data->mats[0]; + + int n_par = 1; + if(nu < 0){ + est_nu = 1; + n_par = 2; + } + + if (theta) { + ltau = theta[0]; + tau = exp(ltau); + if(est_nu == 1){ + lnu = theta[1]; + nu = forward_nu(lnu, nu_upper_bound); + } + } else { + ltau = tau = NAN; + if(est_nu == 1){ + nu = lnu = NAN; + } + } + if (!(data->cache)) { + #pragma omp critical (Name_7c3b4712ebb2dda8def3a5273e2a7e6cf1794b5d) + if (!(data->cache)) { + data->cache = (void **) Calloc(CGENERIC_CACHE_LEN(data), my_cache_tp *); + } + } + + int cache_idx; + CGENERIC_CACHE_ASSIGN_IDX(cache_idx, data); + my_cache_tp *cache = ((my_cache_tp **) data->cache)[cache_idx]; + + switch (cmd) { + + case INLA_CGENERIC_VOID: + { + assert(!(cmd == INLA_CGENERIC_VOID)); + break; + } + + case INLA_CGENERIC_GRAPH: + { + k = 2; + ret = Calloc(k + 2 * M, double); + ret[0] = N; + ret[1] = M; + + for (i = 0; i < M; i++) { + ret[k++] = Q->i[i]; // Row indices + } + for (i = 0; i < M; i++) { + ret[k++] = Q->j[i]; // Column indices + } + break; + } + + case INLA_CGENERIC_Q: + { + k = 2; + ret = Calloc(k + M, double); // Adjust based on Q matrix size + ret[0] = -1; // Required value + ret[1] = M; + if(cache) { + if( (est_nu && cache->theta[0] == theta[0] && cache->theta[1] == theta[1]) || + cache->theta[0] == theta[0]) { + memcpy(ret + 2, cache->Q, M * sizeof(double)); + } else { + compute_Q_fintrinsic(tau, nu, C, Ci, G, Q, &ret[k], rspde_order, + rational_table, est_nu, d, 1, 0, 0, + const_store, mu_store, scaling, D); + } + } else { + compute_Q_fintrinsic(tau, nu, C, Ci, G, Q, &ret[k], rspde_order, + rational_table, est_nu, d, 1, 0, 0, + const_store, mu_store, scaling, D); + } + break; + } + + case INLA_CGENERIC_MU: + { + if (mean_correction) { + ret = Calloc(1 + N, double); + ret[0] = N; /* REQUIRED */ + if(cache) { + if((est_nu && cache->theta[0] == theta[0] && cache->theta[1] == theta[1]) || + cache->theta[0] == theta[0]) { + memcpy(ret + 1, cache->mu, N * sizeof(double)); + } else { + Q_store = Calloc(M, double); + const_store = Calloc(1, double); + compute_Q_fintrinsic(tau, nu, C, Ci, G, Q, Q_store, rspde_order, + rational_table, est_nu, d, 1, 1, 1, + const_store, &ret[1], scaling, D); + memcpy(cache->Q, Q_store, M * sizeof(double)); + memcpy(cache->mu, ret + 1, N * sizeof(double)); + memcpy(cache->lconst, const_store, sizeof(double)); + memcpy(cache->theta, theta, n_par * sizeof(double)); + } + } else { + Q_store = Calloc(M, double); + const_store = Calloc(1, double); + compute_Q_fintrinsic(tau, nu, C, Ci, G, Q, Q_store, rspde_order, + rational_table, est_nu, d, 1, 1, 1, + const_store, &ret[1], scaling, D); + ((my_cache_tp **) data->cache)[cache_idx] = cache = Calloc(1, my_cache_tp); + cache->Q = Calloc(M, double); + cache->mu = Calloc(N, double); + cache->lconst = Calloc(1, double); + cache->theta = Calloc(n_par, double); + memcpy(cache->Q, Q_store, M * sizeof(double)); + memcpy(cache->mu, ret + 1, N * sizeof(double)); + memcpy(cache->lconst, const_store, sizeof(double)); + memcpy(cache->theta, theta, n_par * sizeof(double)); + } + } else { + ret = Calloc(1, double); + ret[0] = 0.0; + } + break; + } + + case INLA_CGENERIC_INITIAL: + { + if(est_nu == 1){ + ret = Calloc(3, double); + ret[0] = 2; // Number of hyperparameters + } else{ + ret = Calloc(2, double); + ret[0] = 1; + } + + ret[1] = prior_tau_mean; + if(est_nu == 1){ + ret[2] = inverse_lnu(start_nu, nu_upper_bound); + } + break; + } + + case INLA_CGENERIC_LOG_NORM_CONST: + { + ret = Calloc(1, double); + if(cache) { + if((est_nu && cache->theta[0] == theta[0] && cache->theta[1] == theta[1]) || + cache->theta[0] == theta[0]) { + memcpy(ret, cache->lconst, sizeof(double)); + } else { + if(mean_correction == 1) { + Q_store = Calloc(M, double); + mu_store = Calloc(N, double); + compute_Q_fintrinsic(tau, nu, C, Ci, G, Q, Q_store, rspde_order, + rational_table, est_nu, d, 1, 1, 1, + &ret[0], mu_store, scaling, D); + + //then cache them and update theta + memcpy(cache->Q, Q_store, M * sizeof(double)); + memcpy(cache->mu, mu_store, N * sizeof(double)); + memcpy(cache->lconst, ret, sizeof(double)); + memcpy(cache->theta, theta, n_par * sizeof(double)); + } else { + Q_store = Calloc(M, double); + compute_Q_fintrinsic(tau, nu, C, Ci, G, Q, Q_store, rspde_order, + rational_table, est_nu, d, 1, 0, 1, + &ret[0], mu_store, scaling, D); + + //then cache them and update theta + memcpy(cache->Q, Q_store, M * sizeof(double)); + memcpy(cache->lconst, ret, sizeof(double)); + memcpy(cache->theta, theta, n_par * sizeof(double)); + } + } + } else { + if(mean_correction == 1) { + Q_store = Calloc(M, double); + mu_store = Calloc(N, double); + compute_Q_fintrinsic(tau, nu, C, Ci, G, Q, Q_store, rspde_order, + rational_table, est_nu, d, 1, 1, 1, + &ret[0], mu_store, scaling, D); + + //then cache them and update theta + ((my_cache_tp **) data->cache)[cache_idx] = cache = Calloc(1, my_cache_tp); + cache->Q = Calloc(M, double); + cache->mu = Calloc(N, double); + cache->lconst = Calloc(1, double); + cache->theta = Calloc(n_par, double); + memcpy(cache->Q, Q_store, M * sizeof(double)); + memcpy(cache->mu, mu_store, N * sizeof(double)); + memcpy(cache->lconst, ret, sizeof(double)); + memcpy(cache->theta, theta, n_par * sizeof(double)); + } else { + Q_store = Calloc(M, double); + compute_Q_fintrinsic(tau, nu, C, Ci, G, Q, Q_store, rspde_order, + rational_table, est_nu, d, 1, 0, 1, + &ret[0], mu_store, scaling, D); + + //then cache them and update theta + ((my_cache_tp **) data->cache)[cache_idx] = cache = Calloc(1, my_cache_tp); + cache->Q = Calloc(M, double); + cache->lconst = Calloc(1, double); + cache->theta = Calloc(n_par, double); + memcpy(cache->Q, Q_store, M * sizeof(double)); + memcpy(cache->lconst, ret, sizeof(double)); + memcpy(cache->theta, theta, n_par * sizeof(double)); + } + } + break; + } + + case INLA_CGENERIC_LOG_PRIOR: + { + ret = Calloc(1, double); + + ret[0] = 0.5 * log(prior_tau_prec) - 0.5 * log(2.0*M_PI); + ret[0] += 0.5 * SQR(prior_tau_prec * (prior_tau_mean - ltau)); + + if(est_nu == 1){ + if(!strcasecmp(prior_nu_dist, "lognormal")){ + ret[0] += -0.5 * SQR(lnu - prior_nu_loglocation)/(SQR(prior_nu_logscale)); + ret[0] += -log(prior_nu_logscale) - 0.5 * log(2.0*M_PI); + ret[0] -= log(pnorm(log(nu_upper_bound), prior_nu_loglocation, prior_nu_logscale)); + } else { // if(!strcasecmp(prior_nu_dist, "beta")){ + double s_1 = (prior_nu_mean / nu_upper_bound) * prior_nu_prec; + double s_2 = (1 - prior_nu_mean / nu_upper_bound) * prior_nu_prec; + ret[0] += logdbeta(nu / nu_upper_bound, s_1, s_2) - log(nu_upper_bound); + } + } + + break; + } + + case INLA_CGENERIC_QUIT: + default: + break; + } + + return ret; +} \ No newline at end of file diff --git a/vignettes/intrinsic.Rmd b/vignettes/intrinsic.Rmd index 283cda59..b5cef038 100644 --- a/vignettes/intrinsic.Rmd +++ b/vignettes/intrinsic.Rmd @@ -32,12 +32,15 @@ package. ## Model specification and simulation -The intrinsic models are defined as +The most general intrinsic model is defined as $$ (-\Delta)^{\beta/2}(\kappa^2-\Delta)^{\alpha/2}(\tau u) = \mathcal{W}, $$ where $\alpha + \beta > d/2$ and $d$ is the dimension of the spatial domain. These models are handled -by performing two rational approximations, one for each fractional operator. +by performing two rational approximations, one for each fractional operator. The package also contains a separate implementation for the special case +$$ +(-\Delta)^{\beta/2}(\tau u) = \mathcal{W}. +$$ To illustrate these models, we begin by defining a mesh over $[0,2]\times [0, 2]$: @@ -54,7 +57,8 @@ plot(mesh_2d, main = "") ``` We now use the `intrinsic.matern.operators()` function to construct the `rSPDE` representation -of the model. +of the general model. One can similarly use `intrinsic.operators()` to construct the `rSPDE` representation +of the special case with $\alpha = 0$. ```{r, message=FALSE} library(rSPDE) @@ -320,4 +324,93 @@ posterior_df_fit <- gg_df(result_fit) ggplot(posterior_df_fit) + geom_line(aes(x = x, y = y)) + facet_wrap(~parameter, scales = "free") + labs(y = "Density") -``` \ No newline at end of file +``` + +## Fitting fractional models with `R-INLA` + +Although the general intrinsic model currently only can be fitted to data in `R-INLA` having fixed values of +$alpha$ and $\beta$, the special case with $\alpha=0$ can be estimated with general smoothness parameter. We now illustrate this. We start by simulating some data using the `intrinsic.operators()` specification + +```{r, message=FALSE} +tau <- 1 +beta <- 1.4 +op <- intrinsic.operators(tau = tau, beta = beta, mesh = mesh_2d) +u <- simulate(op,nsim = 1) +field <- fm_evaluate(proj, field = as.vector(u)) +field.df <- data.frame(x1 = proj$lattice$loc[,1], + x2 = proj$lattice$loc[,2], + y = as.vector(field)) + +ggplot(field.df, aes(x = x1, y = x2, fill = y)) + + geom_raster() + + scale_fill_viridis() + +``` +We now generate some data for this model as above. + +```{r} +n_loc <- 2000 +loc_2d_mesh <- matrix(2*runif(n_loc * 2), n_loc, 2) +A <- spde.make.A( + mesh = mesh_2d, + loc = loc_2d_mesh +) +sigma.e <- 0.1 +y <- A %*% u + rnorm(n_loc) * sigma.e +``` + +The generated data can be seen in the following image. + +```{r,fig.align = "center", echo=TRUE} +df <- data.frame(x1 = as.double(loc_2d_mesh[, 1]), + x2 = as.double(loc_2d_mesh[, 2]), y = as.double(y)) +ggplot(df, aes(x = x1, y = x2, col = y)) + + geom_point() + + scale_color_viridis() +``` + + +We will now fit the model using our `r inla_link()` implementation of the +rational SPDE approach. As we will estimate the smoothness, we have to use the `rspde.make.index` and `rspde.make.A` functions when creating the stack. + +```{r} +mesh.index <- rspde.make.index(name = "field", mesh = mesh_2d) +Abar <- rspde.make.A(mesh = mesh_2d, loc = loc_2d_mesh) +st.dat <- inla.stack(data = list(y = as.vector(y)), A = Abar, effects = mesh.index) +``` + +We now create the model object. + +```{r} +rspde_model <- rspde.intrinsic(mesh = mesh_2d) + +``` + +Finally, we create the formula and fit the model to the data: + +```{r message=FALSE, warning=FALSE} +f <- y ~ -1 + f(field, model = rspde_model) +rspde_fit <- inla(f, + data = inla.stack.data(st.dat), + family = "gaussian", + control.predictor = list(A = inla.stack.A(st.dat))) + +``` + + +To compare the estimated parameters to the true parameters, we can do the following: +```{r} +result_fit <- rspde.result(rspde_fit, "field", rspde_model) +summary(result_fit) +tau <- op$tau +nu <- op$beta - 1 #beta = nu + d/2 +result_df <- data.frame( + parameter = c("tau", "nu", "sigma.e"), + true = c(tau, nu, sigma.e), + mean = c(result_fit$summary.tau$mean,result_fit$summary.nu$mean, + sqrt(1/rspde_fit$summary.hyperpar[1,1])), + mode = c(result_fit$summary.tau$mode, result_fit$summary.nu$mode, + sqrt(1/rspde_fit$summary.hyperpar[1,6])) +) +print(result_df) +```