* using log directory ‘/data/gannet/ripley/R/packages/tests-LENGTH1/MultiGHQuad.Rcheck’
* using R Under development (unstable) (2022-04-26 r82260)
* using platform: x86_64-pc-linux-gnu (64-bit)
* using session charset: UTF-8
* using option ‘--no-stop-on-test-error’
* checking for file ‘MultiGHQuad/DESCRIPTION’ ... OK
* checking extension type ... Package
* this is package ‘MultiGHQuad’ version ‘1.2.0’
* checking package namespace information ... OK
* checking package dependencies ... OK
* checking if this is a source package ... OK
* checking if there is a namespace ... OK
* checking for executable files ... OK
* checking for hidden files and directories ... OK
* checking for portable file names ... OK
* checking for sufficient/correct file permissions ... OK
* checking whether package ‘MultiGHQuad’ can be installed ... [10s/25s] OK
* checking package directory ... OK
* checking DESCRIPTION meta-information ... OK
* checking top-level files ... OK
* checking for left-over files ... OK
* checking index information ... OK
* checking package subdirectories ... OK
* checking R files for non-ASCII characters ... OK
* checking R files for syntax errors ... OK
* checking whether the package can be loaded ... OK
* checking whether the package can be loaded with stated dependencies ... OK
* checking whether the package can be unloaded cleanly ... OK
* checking whether the namespace can be loaded with stated dependencies ... OK
* checking whether the namespace can be unloaded cleanly ... OK
* checking loading without being on the library search path ... OK
* checking use of S3 registration ... OK
* checking dependencies in R code ... OK
* checking S3 generic/method consistency ... OK
* checking replacement functions ... OK
* checking foreign function calls ... OK
* checking R code for possible problems ... [13s/40s] OK
* checking Rd files ... OK
* checking Rd metadata ... OK
* checking Rd line widths ... OK
* checking Rd cross-references ... OK
* checking for missing documentation entries ... OK
* checking for code/documentation mismatches ... OK
* checking Rd \usage sections ... OK
* checking Rd contents ... OK
* checking for unstated dependencies in examples ... OK
* checking examples ... [4s/12s] ERROR
Running examples in ‘MultiGHQuad-Ex.R’ failed
The error most likely occurred in:

> ### Name: eval.quad
> ### Title: Evaluation of multivariate normal distributed expectations
> ### Aliases: eval.quad
> 
> ### ** Examples
> 
> ### Basic example; E(X), X ~ N(0,1)
> grid <- init.quad(Q = 1, prior = list(mu = 0, Sigma = diag(1)))
> eval.quad(X = grid)
[1] -2.178954e-16
attr(,"variance")
     [,1]
[1,]    1
> 
> ### Example; Rasch model person parameter
> # E(theta), theta ~ N(0,1) * P(X = 1 | theta, beta), P is simplified rasch model
> # set up rasch model with fixed beta, returns LL
> rasch <- function(theta, beta, responses){
+   p <- exp(theta - beta)/(1 + exp(theta - beta))
+   q <- 1 - p
+   return(log(p) * sum(responses == 1) + log(q) * sum(responses == 0))
+ }
> 
> # when theta == beta, P(X = 1) = .5, generate some bernoulli trials with p = .5
> responses <- rbinom(5, 1, .5)
> 
> # get EAP estimate for theta, prior N(0,1)
> eval.quad(rasch, grid, beta = 0, responses = responses)
[1] -0.2471169
attr(,"variance")
          [,1]
[1,] 0.4852284
> 
> # with more data, the estimate becomes more accurate, and variance decreases
> eval.quad(rasch, grid, beta = 0, responses = rbinom(20, 1, .5))
[1] 0.3393648
attr(,"variance")
          [,1]
[1,] 0.2671098
> eval.quad(rasch, grid, beta = 0, responses = rbinom(50, 1, .5))
[1] -0.5869554
attr(,"variance")
          [,1]
[1,] 0.0358114
> eval.quad(rasch, grid, beta = 0, responses = rbinom(100, 1, .5))
[1] -1.80146e-16
attr(,"variance")
         [,1]
[1,] 0.380327
> 
> ### problem; the result starts to 'snap' to the closest quadrature point when
> # the posterior distribution is too dissimilar to the prior.
> evals <- eval.quad(rasch, grid, beta = 0, responses = rbinom(100, 1, .5), debug = TRUE)
> evals.values <- attr(evals, "values")
> 
> # posterior density after 40 items
> p <- plot(function(x) exp(dnorm(x, log = TRUE) + 
+                             rasch(x, beta = 1, responses = rbinom(100, 1, .5))),
+           from = -3, to = 3)
> 
> # quadrature points used
> points(grid$X, exp(grid$W)*max(p$y), pch = 20)
> 
> # the evaluation relies almost completely on one quadrature point, 
> # which causes results to 'snap' to that point.
> # we could add more quadrature points...
> grid2 <- init.quad(Q = 1, ip = 20)
> points(grid2$X, exp(grid2$W)*max(p$y), pch = 20, col = "grey")
> 
> # but if the posterior is not centered on the prior, this quickly fails:
> p <- plot(function(x) exp(dnorm(x, log = TRUE) + 
+                             rasch(x, beta = 2, responses = rbinom(100, 1, .5))),
+           from = -3, to = 3)
> points(grid2$X, exp(grid2$W)*max(p$y), pch = 20, col = "grey")
> 
> # additionally, adding extra quadrature points in a multidimensional 
> # problem quickly grows out of control.
> 
> ### a better solution; adaptive quadrature grid.
> # say we have an idea of where our parameter is located, through another estimator, 
> # or a previous estimate.
> # we can then use this to adapt where our quadrature grid should be.
> # get an estimate;
> responses <- rbinom(10, 1, .5)
> est <- eval.quad(rasch, grid, beta = 2, responses = responses)
> print( est )
[1] 1.751766
attr(,"variance")
          [,1]
[1,] 0.2102369
> 
> # adapt the grid;
> grid3 <- init.quad(Q = 1, adapt = est)
 ----------- FAILURE REPORT -------------- 
 --- failure: length > 1 in coercion to logical ---
 --- srcref --- 
: 
 --- package (from environment) --- 
MultiGHQuad
 --- call from context --- 
init.quad(Q = 1, adapt = est)
 --- call from argument --- 
!is.list(adapt) || length(adapt$mu) != Q || dim(adapt$Sigma) != 
    c(Q, Q)
 --- R stacktrace ---
where 1: init.quad(Q = 1, adapt = est)

 --- value of length: 2 type: logical ---
[1] FALSE FALSE
 --- function from context --- 
function (Q = 2, prior = list(mu = rep(0, Q), Sigma = diag(Q)), 
    adapt = NULL, ip = 6, prune = FALSE, forcePD = FALSE, debug = FALSE) 
{
    if (!is.null(adapt) && !is.null(attr(adapt, "variance"))) {
        adapt <- list(mu = adapt, Sigma = attr(adapt, "variance"))
    }
    if (!is.null(adapt) && (!is.list(adapt) || length(adapt$mu) != 
        Q || dim(adapt$Sigma) != c(Q, Q))) 
        stop("Size or format of Adapt argument invalid.")
    x <- fastGHQuad::gaussHermiteData(ip)
    w <- x$w/sqrt(pi)
    x <- x$x * sqrt(2)
    X <- as.matrix(expand.grid(lapply(apply(replicate(Q, x), 
        2, list), unlist)))
    trans <- function(X, Sigma) {
        if (forcePD) 
            Sigma <- nearPD(Sigma)$mat
        lambda <- with(eigen(Sigma), {
            if (any(values < 0)) 
                warning("Matrix is not positive definite.")
            if (length(values) > 1) 
                vectors %*% diag(sqrt(values))
            else vectors * sqrt(values)
        })
        t(lambda %*% t(X))
    }
    if (debug) {
        noise <- rmvnorm(10000, prior$mu, prior$Sigma)
        plot(noise, col = "grey", pch = 20, cex = 0.2, xlim = range(noise[, 
            1]) * 1.1, ylim = range(noise[, 2]) * 1.1)
        points(X, pch = 20)
    }
    g <- as.matrix(expand.grid(lapply(apply(replicate(Q, w), 
        2, list), unlist)))
    W <- apply(g, 1, function(x) sum(log(x)))
    if (debug) {
        points(X, pch = 20, col = "green", cex = exp(W)/max(exp(W)) * 
            5)
        cat("\ngreen is transformed to prior")
    }
    if (is.null(adapt)) {
        X <- trans(X, prior$Sigma)
        X <- t(t(X) + prior$mu)
    }
    else {
        X <- trans(X, adapt$Sigma)
        X <- t(t(X) + adapt$mu)
        if (debug) {
            points(X, pch = 20, col = "blue", cex = exp(W)/max(exp(W)) * 
                5)
            cat("\nblue is transformed to posterior")
        }
        adapt$chol <- chol(adapt$Sigma)
        adapt$det <- sum(log(diag(adapt$chol)))
        adapt$aux <- colSums(backsolve(adapt$chol, t(X) - adapt$mu, 
            transpose = TRUE)^2)
        prior$chol <- chol(prior$Sigma)
        prior$det <- sum(log(diag(prior$chol)))
        prior$aux <- colSums(backsolve(prior$chol, t(X) - prior$mu, 
            transpose = TRUE)^2)
        fact <- (adapt$aux - prior$aux)/2 + adapt$det - prior$det
        W <- W + fact
    }
    if (debug) {
        points(X, pch = 20, col = "orange", cex = exp(W)/max(exp(W)) * 
            5)
        cat("\norange is mean shifted")
    }
    if (prune) {
        threshold <- log(min(w)^(Q - 1) * max(w))
        relevant <- W >= threshold
        W <- W[relevant]
        X <- X[relevant, , drop = FALSE]
        if (debug) {
            points(X, pch = 20, col = "pink", cex = exp(W)/max(exp(W)) * 
                5)
            cat("\npink is pruned.")
        }
    }
    return(invisible(list(X = X, W = W)))
}
<bytecode: 0xa0be698>
<environment: namespace:MultiGHQuad>
 --- function search by body ---
Function init.quad in namespace MultiGHQuad has this body.
 ----------- END OF FAILURE REPORT -------------- 
Fatal error: length > 1 in coercion to logical
* checking for unstated dependencies in ‘tests’ ... OK
* checking tests ...
  Running ‘testthat.R’ [25s/63s]
 [26s/64s] ERROR
Running the tests in ‘tests/testthat.R’ failed.
Complete output:
  > library(testthat)
  > library(MultiGHQuad)
  Loading required package: mvtnorm
  Loading required package: Matrix
  > 
  > test_check("MultiGHQuad")
   ----------- FAILURE REPORT -------------- 
   --- failure: length > 1 in coercion to logical ---
   --- srcref --- 
  : 
   --- package (from environment) --- 
  MultiGHQuad
   --- call from context --- 
  init.quad(2, adapt = estimate)
   --- call from argument --- 
  !is.list(adapt) || length(adapt$mu) != Q || dim(adapt$Sigma) != 
      c(Q, Q)
   --- R stacktrace ---
  where 1: init.quad(2, adapt = estimate)
  where 2: eval_bare(expr, quo_get_env(quo))
  where 3: quasi_label(enquo(object), label, arg = "object")
  where 4 at test_generate_quadpoints.R#22: expect_is(init.quad(2, adapt = estimate), "list")
  where 5: eval(code, test_env)
  where 6: eval(code, test_env)
  where 7: withCallingHandlers({
      eval(code, test_env)
      if (!handled && !is.null(test)) {
          skip_empty()
      }
  }, expectation = handle_expectation, skip = handle_skip, warning = handle_warning, 
      message = handle_message, error = handle_error)
  where 8: doTryCatch(return(expr), name, parentenv, handler)
  where 9: tryCatchOne(expr, names, parentenv, handlers[[1L]])
  where 10: tryCatchList(expr, names[-nh], parentenv, handlers[-nh])
  where 11: doTryCatch(return(expr), name, parentenv, handler)
  where 12: tryCatchOne(tryCatchList(expr, names[-nh], parentenv, handlers[-nh]), 
      names[nh], parentenv, handlers[[nh]])
  where 13: tryCatchList(expr, classes, parentenv, handlers)
  where 14: tryCatch(withCallingHandlers({
      eval(code, test_env)
      if (!handled && !is.null(test)) {
          skip_empty()
      }
  }, expectation = handle_expectation, skip = handle_skip, warning = handle_warning, 
      message = handle_message, error = handle_error), error = handle_fatal, 
      skip = function(e) {
      })
  where 15: test_code(desc, code, env = parent.frame(), reporter = reporter)
  where 16 at test_generate_quadpoints.R#20: test_that("Adapt accepts previous estimate", {
      expect_error(init.quad(1, adapt = estimate))
      expect_is(init.quad(2, adapt = estimate), "list")
      expect_error(init.quad(3, adapt = estimate))
      expect_is(init.quad(1, adapt = estimateb), "list")
      expect_error(init.quad(2, adapt = estimateb))
  })
  where 17: eval(code, test_env)
  where 18: eval(code, test_env)
  where 19: withCallingHandlers({
      eval(code, test_env)
      if (!handled && !is.null(test)) {
          skip_empty()
      }
  }, expectation = handle_expectation, skip = handle_skip, warning = handle_warning, 
      message = handle_message, error = handle_error)
  where 20: doTryCatch(return(expr), name, parentenv, handler)
  where 21: tryCatchOne(expr, names, parentenv, handlers[[1L]])
  where 22: tryCatchList(expr, names[-nh], parentenv, handlers[-nh])
  where 23: doTryCatch(return(expr), name, parentenv, handler)
  where 24: tryCatchOne(tryCatchList(expr, names[-nh], parentenv, handlers[-nh]), 
      names[nh], parentenv, handlers[[nh]])
  where 25: tryCatchList(expr, classes, parentenv, handlers)
  where 26: tryCatch(withCallingHandlers({
      eval(code, test_env)
      if (!handled && !is.null(test)) {
          skip_empty()
      }
  }, expectation = handle_expectation, skip = handle_skip, warning = handle_warning, 
      message = handle_message, error = handle_error), error = handle_fatal, 
      skip = function(e) {
      })
  where 27: test_code(NULL, exprs, env)
  where 28: source_file(path, child_env(env), wrap = wrap)
  where 29: FUN(X[[i]], ...)
  where 30: lapply(test_paths, test_one_file, env = env, wrap = wrap)
  where 31: doTryCatch(return(expr), name, parentenv, handler)
  where 32: tryCatchOne(expr, names, parentenv, handlers[[1L]])
  where 33: tryCatchList(expr, classes, parentenv, handlers)
  where 34: tryCatch(code, testthat_abort_reporter = function(cnd) {
      cat(conditionMessage(cnd), "\n")
      NULL
  })
  where 35: with_reporter(reporters$multi, lapply(test_paths, test_one_file, 
      env = env, wrap = wrap))
  where 36: test_files(test_dir = test_dir, test_package = test_package, 
      test_paths = test_paths, load_helpers = load_helpers, reporter = reporter, 
      env = env, stop_on_failure = stop_on_failure, stop_on_warning = stop_on_warning, 
      wrap = wrap, load_package = load_package)
  where 37: test_files(test_dir = path, test_paths = test_paths, test_package = package, 
      reporter = reporter, load_helpers = load_helpers, env = env, 
      stop_on_failure = stop_on_failure, stop_on_warning = stop_on_warning, 
      wrap = wrap, load_package = load_package, parallel = parallel)
  where 38: test_dir("testthat", package = package, reporter = reporter, 
      ..., load_package = "installed")
  where 39: test_check("MultiGHQuad")
  
   --- value of length: 2 type: logical ---
  [1] FALSE FALSE
   --- function from context --- 
  function (Q = 2, prior = list(mu = rep(0, Q), Sigma = diag(Q)), 
      adapt = NULL, ip = 6, prune = FALSE, forcePD = FALSE, debug = FALSE) 
  {
      if (!is.null(adapt) && !is.null(attr(adapt, "variance"))) {
          adapt <- list(mu = adapt, Sigma = attr(adapt, "variance"))
      }
      if (!is.null(adapt) && (!is.list(adapt) || length(adapt$mu) != 
          Q || dim(adapt$Sigma) != c(Q, Q))) 
          stop("Size or format of Adapt argument invalid.")
      x <- fastGHQuad::gaussHermiteData(ip)
      w <- x$w/sqrt(pi)
      x <- x$x * sqrt(2)
      X <- as.matrix(expand.grid(lapply(apply(replicate(Q, x), 
          2, list), unlist)))
      trans <- function(X, Sigma) {
          if (forcePD) 
              Sigma <- nearPD(Sigma)$mat
          lambda <- with(eigen(Sigma), {
              if (any(values < 0)) 
                  warning("Matrix is not positive definite.")
              if (length(values) > 1) 
                  vectors %*% diag(sqrt(values))
              else vectors * sqrt(values)
          })
          t(lambda %*% t(X))
      }
      if (debug) {
          noise <- rmvnorm(10000, prior$mu, prior$Sigma)
          plot(noise, col = "grey", pch = 20, cex = 0.2, xlim = range(noise[, 
              1]) * 1.1, ylim = range(noise[, 2]) * 1.1)
          points(X, pch = 20)
      }
      g <- as.matrix(expand.grid(lapply(apply(replicate(Q, w), 
          2, list), unlist)))
      W <- apply(g, 1, function(x) sum(log(x)))
      if (debug) {
          points(X, pch = 20, col = "green", cex = exp(W)/max(exp(W)) * 
              5)
          cat("\ngreen is transformed to prior")
      }
      if (is.null(adapt)) {
          X <- trans(X, prior$Sigma)
          X <- t(t(X) + prior$mu)
      }
      else {
          X <- trans(X, adapt$Sigma)
          X <- t(t(X) + adapt$mu)
          if (debug) {
              points(X, pch = 20, col = "blue", cex = exp(W)/max(exp(W)) * 
                  5)
              cat("\nblue is transformed to posterior")
          }
          adapt$chol <- chol(adapt$Sigma)
          adapt$det <- sum(log(diag(adapt$chol)))
          adapt$aux <- colSums(backsolve(adapt$chol, t(X) - adapt$mu, 
              transpose = TRUE)^2)
          prior$chol <- chol(prior$Sigma)
          prior$det <- sum(log(diag(prior$chol)))
          prior$aux <- colSums(backsolve(prior$chol, t(X) - prior$mu, 
              transpose = TRUE)^2)
          fact <- (adapt$aux - prior$aux)/2 + adapt$det - prior$det
          W <- W + fact
      }
      if (debug) {
          points(X, pch = 20, col = "orange", cex = exp(W)/max(exp(W)) * 
              5)
          cat("\norange is mean shifted")
      }
      if (prune) {
          threshold <- log(min(w)^(Q - 1) * max(w))
          relevant <- W >= threshold
          W <- W[relevant]
          X <- X[relevant, , drop = FALSE]
          if (debug) {
              points(X, pch = 20, col = "pink", cex = exp(W)/max(exp(W)) * 
                  5)
              cat("\npink is pruned.")
          }
      }
      return(invisible(list(X = X, W = W)))
  }
  <bytecode: 0xcaf3248>
  <environment: namespace:MultiGHQuad>
   --- function search by body ---
  Function init.quad in namespace MultiGHQuad has this body.
   ----------- END OF FAILURE REPORT -------------- 
  Fatal error: length > 1 in coercion to logical
* checking PDF version of manual ... OK
* checking for non-standard things in the check directory ... OK
* checking for detritus in the temp directory ... OK
* DONE

Status: 2 ERRORs
See
  ‘/data/gannet/ripley/R/packages/tests-LENGTH1/MultiGHQuad.Rcheck/00check.log’
for details.

Command exited with non-zero status 1
Time 4:55.94, 103.48 + 10.45