Refactor code to remove old int.quadrature function

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: inlabru
 Type: Package
 Title: Bayesian Latent Gaussian Modelling using INLA and Extensions
-Version: 2.11.1.9024
+Version: 2.11.1.9025
 Authors@R: c(
   person("Finn", "Lindgren", email = "[email protected]",
          role = c("aut", "cre", "cph"),
@@ -121,7 +121,6 @@ Collate:
     'ggplot.R'
     'inla.R'
     'inlabru-package.R'
-    'integration.R'
     'local_testthat.R'
     'mappers.R'
     'mapper_repeat.R'

R/bru.gof.R

@@ -18,15 +18,15 @@
 #'   [predict()].
 #' @param observations A vector of observed values
 #' @param breaks A vector of bin boundaries
-#' @param nint Number of integration points per bin. Increase this if the bins
-#'   are wide and
+#' @param nint Number of integration intervals per bin. Increase this if the
+#'   bins are wide and the LGCP is not smooth.
 #' @param probs numeric vector of probabilities with values in `[0,1]`
 #' @param \dots arguments passed on to [predict.bru()]
 #' @return An `data.frame` with a ggplot attribute `ggp`
 #' @importFrom rlang .data
 #' @examples
 #' \dontrun{
-#' if (require(ggplot2) && require(fmesher)) {
+#' if (require(ggplot2) && require(fmesher) && bru_safe_inla()) {
 #'   # Load a point pattern
 #'   data(Poisson2_1D)
 #'
@@ -46,8 +46,12 @@
 #'   # Fit an LGCP model
 #'   x <- seq(0, 55, length.out = 50)
 #'   mesh1D <- fm_mesh_1d(x, boundary = "free")
-#'   mdl <- x ~ spde1D(x, model = inla.spde2.matern(mesh1D)) + Intercept(1)
-#'   fit.spde <- lgcp(mdl, pts2, domain = list(x = c(0, 55)))
+#'   matern <- INLA::inla.spde2.pcmatern(mesh1D,
+#'                                       prior.range = c(1, 0.01),
+#'                                       prior.sigma = c(1, 0.01),
+#'                                       constr = TRUE)
+#'   mdl <- x ~ spde1D(x, model = matern) + Intercept(1)
+#'   fit.spde <- lgcp(mdl, pts2, domain = list(x = mesh1D))
 #'
 #'   # Calculate bin statistics
 #'   bc <- bincount(
@@ -57,7 +61,6 @@
 #'     predictor = x ~ exp(spde1D + Intercept)
 #'   )
 #'
-#'
 #'   # Plot them!
 #'   attributes(bc)$ggp
 #' }
@@ -80,43 +83,51 @@ bincount <- function(result, predictor, observations, breaks, nint = 20,
   mid <- breaks[1:(length(breaks) - 1)] + diff(breaks) / 2
 
   # Create integration points
-  ip <- int.quadrature(
-    breaks[1:(length(breaks) - 1)],
-    breaks[2:(length(breaks))],
-    scheme = "trapezoid",
-    n.points = nint
+  nsub <- nint
+  u <- rep(
+    (seq_len(nsub) - 1L) / (nsub + 1),
+    nbins
   )
-  points <- data.frame(tmp = as.vector(ip$ips))
-  colnames(points) <- all.vars(update.formula(predictor, ~.0))
-  points$bin <- rep(1:nbins, each = nint)
+  domain_breaks <-
+    c(
+      breaks[rep(seq_len(nbins), each = nsub)] * (1 - u) +
+        breaks[rep(seq_len(nbins) + 1, each = nsub)] * u,
+      breaks[length(breaks)]
+    )
+  points <- fmesher::fm_int(samplers = cbind(breaks[1:(length(breaks) - 1)],
+                                             breaks[2:(length(breaks))]),
+                            domain = fmesher::fm_mesh_1d(domain_breaks),
+                            int.args = list(method = "stable",
+                                            nsub1 = 1L),
+                            name = as.character(predictor)[2])
 
   # Sampler
-  smp <- generate(result, points, predictor, ..., format = "matrix")
+  smp <- generate(result,
+                  newdata = points,
+                  formula = predictor,
+                  ...,
+                  format = "matrix")
 
   # Integrate per bin
-  smp <- ip$wl * smp
-  qq <- aggregate(
-    smp,
-    by = list(rep(1:nbins, each = nint)),
-    FUN = sum,
-    drop = TRUE
-  )[, 2:(ncol(smp) + 1)]
-
-  # Normalize bin probabilities
-  for (s in seq_len(ncol(smp))) {
-    qq[, s] <- qq[, s] / sum(qq[, s])
+  qq <- matrix(0.0, nrow = nbins, ncol(smp))
+  for (k in seq_len(ncol(smp))) {
+    qq[, k] <- fmesher::fm_block_logsumexp_eval(points$.block,
+                                                weights = points$weight,
+                                                values = log(smp[, k]),
+                                                log = TRUE)
+    # Normalize bin probabilities
+    qq[, k] <- qq[, k] - max(qq[, k])
+    qq[, k] <- exp(qq[, k] - log(sum(exp(qq[, k]))))
   }
 
   # For each bin calculate predictive interval
   pint <- list()
   for (k in 1:nbins) {
-    xx <- 0:nobs
+    xx <- 0:(nobs + 1L)
     cdff <- function(p) pbinom(xx, size = nobs, prob = p)
     zz <- apply(qq[k, , drop = FALSE], MARGIN = 2, cdff)
     zz <- apply(zz, MARGIN = 1, mean)
-    # Suppress warnings about duplicate zz values; these may appear in the
-    # tails, with zz==0 or zz==1.
-    pint[[k]] <- suppressWarnings(approx(zz, xx, xout = probs, rule = 2)$y)
+    pint[[k]] <- vapply(probs, function(pr) xx[sum(zz < pr) + 1L], 0.0)
   }
   pint <- data.frame(do.call(rbind, pint))
   colnames(pint) <- paste0("q", probs)

inst/examples/spde.posterior.R

@@ -17,7 +17,7 @@ if (bru_safe_inla() && require(ggplot2, quietly = TRUE)) {
 
   x <- seq(0, 55, length.out = 20)
   mesh1D <- fm_mesh_1d(x, boundary = "free")
-  mdl <- x ~ spde1D(x, model = INLA::inla.spde2.matern(mesh1D)) + Intercept(1)
+  mdl <- x ~ spde1D(x, model = INLA::inla.spde2.matern(mesh1D, constr = TRUE)) + Intercept(1)
   fit <- lgcp(mdl, data = pts2, domain = list(x = mesh1D))
 
   # Calculate and plot the posterior range