Merge pull request #241 from kenkellner/parboot_update

Add error handling to parboot

3 files changed, 101 insertions, 82 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 1bbd46b..a878727 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: unmarked
-Version: 1.2.5.9006
-Date: 2022-08-24
+Version: 1.2.5.9007
+Date: 2022-10-12
 Type: Package
 Title: Models for Data from Unmarked Animals
 Authors@R: c(
diff --git a/R/boot.R b/R/boot.R
index 37f0fd5..2f8cffc 100644
--- a/R/boot.R
+++ b/R/boot.R
@@ -30,81 +30,66 @@ setMethod("replaceY", "unmarkedFrameOccuMulti",
       object
 })
 
-setMethod("parboot", "unmarkedFit",
-    function(object, statistic=SSE, nsim=10, report, seed = NULL, parallel = TRUE, ncores, ...)
-{
-    dots <- list(...)
-    statistic <- match.fun(statistic)
-    call <- match.call(call = sys.call(-1))
-    formula <- object@formula
-    umf <- getData(object)
-    y <- getY(object)
-    ests <- as.numeric(coef(object))
-    starts <- ests
-    if(methods::.hasSlot(object, "TMB") && !is.null(object@TMB)) starts <- NULL
-    t0 <- statistic(object, ...)
-    lt0 <- length(t0)
-    t.star <- matrix(NA, nsim, lt0)
-    if(!missing(report))
-        cat("t0 =", t0, "\n")
-    simdata <- umf
-    if (!is.null(seed)) set.seed(seed)
-    simList <- simulate(object, nsim = nsim, na.rm = FALSE)
-    availcores <- detectCores()
-    if(missing(ncores)) ncores <- availcores - 1
-    if(ncores > availcores) ncores <- availcores
-
-    no_par <- ncores < 2 || nsim < 100 || !parallel
-
-    if (no_par) {
-      if (!missing(report)) {
-        for(i in 1:nsim) {
-          simdata <- replaceY(simdata, simList[[i]])
-          fit <- update(object, data=simdata, starts=starts, se=FALSE)
-          t.star[i,] <- statistic(fit, ...)
-          if(!missing(report)) {
-            if (nsim > report && i %in% seq(report, nsim, by=report))
-              cat("iter", i, ": ", t.star[i, ], "\n")
-          }
-        }
-      } else {
-        t.star <- pbsapply(1:nsim, function(i) {
-          simdata <- replaceY(simdata, simList[[i]])
-          fit <- update(object, data=simdata, starts=starts, se=FALSE)
-          t.star.tmp <- statistic(fit, ...)
-        })
-        if (lt0 > 1)
-          t.star <- t(t.star)
-        else
-          t.star <- matrix(t.star, ncol = lt0)
-      }
-    } else {
-      message("Running parametric bootstrap in parallel on ", ncores, " cores.")
-      if (!missing(report)) message("Bootstrapped statistics not reported during parallel processing.")
-      cl <- makeCluster(ncores)
-      if (!is.null(seed)) parallel::clusterSetRNGStream(cl, iseed = seed)
-      on.exit(stopCluster(cl))
-      varList <- c("simList", "y", "object", "simdata", "starts", "statistic", "dots")
-      # If call formula is an object, include it too
-      fm.nms <- all.names(object@call)
-      if (!any(grepl("~", fm.nms))) varList <- c(varList, fm.nms[2])
-      ## Hack to get piFun for unmarkedFitGMM and unmarkedFitMPois
-      if(.hasSlot(umf, "piFun")) varList <- c(varList, umf@piFun)
-      clusterExport(cl, varList, envir = environment())
-      clusterEvalQ(cl, library(unmarked))
-      clusterEvalQ(cl, list2env(dots))
-      t.star.parallel <- pblapply(1:nsim, function(i) {
-        simdata <- replaceY(simdata, simList[[i]])
-        fit <- update(object, data = simdata, starts = starts, se = FALSE)
-        t.star <- statistic(fit, ...)
-      }, cl = cl)
-      t.star <- matrix(unlist(t.star.parallel), nrow = length(t.star.parallel), byrow = TRUE)
-    }
-    if (!is.null(names(t0)))
-      colnames(t.star) <- names(t0)
-    else colnames(t.star) <- paste("t*", 1:lt0, sep="")
-    out <- new("parboot", call = call, t0 = t0, t.star = t.star)
-    return(out)
+
+setMethod("parboot", "unmarkedFit",  function(object, statistic=SSE, nsim=10,
+          report, seed = NULL, parallel = FALSE, ncores, ...){
+
+  if(!missing(report)){
+    warning("report argument is non-functional and will be deprecated in the next version", call.=FALSE)
+  }
+
+  dots <- list(...)
+  call <- match.call(call = sys.call(-1))
+  stopifnot(is.function(statistic))
+  starts <- as.numeric(coef(object))
+  # Get rid of starting values if model was fit with TMB
+  if(methods::.hasSlot(object, "TMB") && !is.null(object@TMB)) starts <- NULL
+
+  t0 <- statistic(object, ...)
+
+  simList <- simulate(object, nsim = nsim, na.rm = FALSE)
+
+  availcores <- parallel::detectCores() - 1
+  if(missing(ncores) || ncores > availcores) ncores <- availcores
+
+  cl <- NULL
+  if(parallel){
+    cl <- parallel::makeCluster(ncores)
+    on.exit(parallel::stopCluster(cl))
+    parallel::clusterEvalQ(cl, library(unmarked))
+    env_vars <- c("dots", "replaceY")
+    fm.nms <- all.names(object@call)
+    if (!any(grepl("~", fm.nms))) env_vars <- c(env_vars, fm.nms[2])
+    if(.hasSlot(object@data, "piFun")) env_vars <- c(env_vars, object@data@piFun)
+    parallel::clusterExport(cl, env_vars, envir = environment())
+    parallel::clusterEvalQ(cl, list2env(dots))
+  }
+
+  run_sim <- function(x, object, statistic, starts, t0, ...){
+    simdata <- replaceY(object@data, x)
+    tryCatch({
+      #if(runif(1,0,1) < 0.5) stop("fail") # for testing error trapping
+      fit <- update(object, data=simdata, starts=starts, se=FALSE)
+      statistic(fit, ...)
+    }, error=function(e){
+      t0[] <- NA
+      t0
+    })
+  }
+
+  t.star <- t(pbapply::pbsapply(simList, run_sim, object=object,
+                              statistic=statistic, starts=starts, t0=t0,
+                              cl=cl, ...))
+  if(length(t0) == 1) t.star <- matrix(t.star, ncol=1)
+
+  failed <- apply(t.star, 1, function(x) any(is.na(x)))
+  if(sum(failed) > 0){
+    warning(paste0("Model fitting failed in ",sum(failed), " sims."), call.=FALSE)
+    t.star <- t.star[!failed,,drop=FALSE]
+  }
+
+  new("parboot", call = call, t0 = t0, t.star = t.star)
+
 })
 
 
diff --git a/tests/testthat/test_parboot.R b/tests/testthat/test_parboot.R
index c3ece25..af62658 100644
--- a/tests/testthat/test_parboot.R
+++ b/tests/testthat/test_parboot.R
@@ -32,15 +32,49 @@ test_that("parboot works", {
   dev.off()
   expect_equal(pl, NULL)
 
-  # check that report works
-  rep_output <- capture.output(parboot(fm, fitstats, nsim=3, report=TRUE))
-  expect_equal(substr(rep_output[1], 1,2), "t0")
+  # check that report argument gives warning
+  expect_warning(parboot(fm, fitstats, nsim=3, report=TRUE))
 })
 
 test_that("parboot works in parallel",{
   skip_on_cran()
   skip_on_ci()
   # check parallel
-  pb <- parboot(fm, nsim=101, parallel=TRUE, ncores=2)
-  expect_equal(length(pb@t.star), 101)
+  pb <- parboot(fm, nsim=10, parallel=TRUE, ncores=2)
+  expect_equal(length(pb@t.star), 10)
+})
+
+test_that("parboot handles failing model fits", {
+
+  fail_func <- function(x){
+    rand <- rnorm(1)
+    if(rand > 0.5){
+      stop("fail")
+    }
+    return(rand)
+  }
+
+  set.seed(123)
+  expect_warning(pb <- parboot(fm, nsim=20, statistic=fail_func))
+  expect_equal(nrow(pb@t.star), 13)
+
+  expect_warning(pb <- parboot(fm, nsim=20, statistic=fail_func, parallel=TRUE))
+  expect_true(nrow(pb@t.star) < 20)
+
+})
+
+test_that("parboot handles statistic functions with additional arguments", {
+
+  opt_func <- function(x, y){
+    res <- mean(residuals(x), na.rm=TRUE)
+    c(res=res, y=y)
+  }
+
+  pb <- parboot(fm, nsim=10, statistic=opt_func, y=0.1)
+  expect_equal(colnames(pb@t.star), c("res", "y"))
+  expect_true(all(pb@t.star[,"y"]==0.1))
+
+  pb <- parboot(fm, nsim=10, statistic=opt_func, y=0.1, parallel=TRUE)
+  expect_equal(colnames(pb@t.star), c("res", "y"))
+  expect_true(all(pb@t.star[,"y"]==0.1))
 })