fixes to sfd and vignettes

R/simulate_flocker_data.R

@@ -28,14 +28,53 @@
 #'   be simulated without any covariate influence on occupancy.
 #' @param rep_constant logical: create data with unit covariates only (TRUE) 
 #'   or data that includes event covariates (FALSE)
-#' @param params a named list containing of parameter values to use in simulation.
-#'   Any required parameters whose names are not in this list will be assigned their
-#'   default values. To see the parameter names and structures required, run with 
-#'   `params = NULL` (the default) and examine the `$params` element of the output.
-#' @param covariates a dataframe of covariate values to use in simulation, or
-#'   NULL to simulate values. To see the covariate names and structures required,
-#'   run with `covariates = NULL` (the default) and examine the `$covariates` element
-#'   of the output.
+#' @param params A named list controlling the coefficients used to generate the
+#'   latent occupancy state and detections. Supported names are:
+#'   \itemize{
+#'     \item \code{coefs}: A numeric data.frame of per-species coefficient
+#'       values. Must have \code{n_sp} rows and
+#'       one column for each coefficient required by the requested model type.
+#'       Column names must match the required set (order does not matter).
+#'       If supplied, \code{coef_means} and \code{Sigma} must not be supplied.
+#'     \item \code{coef_means}: Numeric vector of length \code{n_coef} giving the
+#'       mean of the multivariate normal distribution used to simulate per-species
+#'       coefficients when \code{coefs} is not supplied.
+#'     \item \code{Sigma}: Numeric \code{n_coef x n_coef} symmetric covariance
+#'       matrix for the multivariate normal distribution used to simulate per-species
+#'       coefficients when \code{coefs} is not supplied.
+#'   }
+#'
+#'   If \code{params = NULL} (the default), coefficients are simulated from a
+#'   multivariate normal prior with \code{coef_means = rep(0, n_coef)} and a
+#'   default \code{Sigma}. The dimension \code{n_coef} and the required coefficient
+#'   names depend on the requested model type:
+#'   \describe{
+#'     \item{Always}{\code{det_intercept}, \code{det_slope_unit}}
+#'     \item{If \code{rep_constant = FALSE}}{\code{det_slope_visit}}
+#'     \item{If \code{n_season == 1} or \code{multi_init == "explicit"}}{\code{occ_intercept}}
+#'     \item{If \code{n_season == 1} or \code{multi_init == "explicit"}, and \code{augmented = FALSE}}{\code{occ_slope_unit}}
+#'     \item{If \code{n_season > 1}}{\code{col_intercept}, \code{col_slope_unit}}
+#'     \item{If \code{multiseason == "colex"}}{\code{ex_intercept}, \code{ex_slope_unit}}
+#'     \item{If \code{multiseason == "autologistic"}}{\code{auto_intercept}, \code{auto_slope_unit}}
+#'   }
+#'
+#'   Notes:
+#'   \itemize{
+#'     \item When \code{augmented = TRUE}, \code{occ_slope_unit} is not used and
+#'       occupancy is simulated with no covariate effect.
+#'     \item Any names in \code{params} other than \code{coefs}, \code{coef_means},
+#'       and \code{Sigma} are rejected with an error.
+#'   }
+#' @param covariates A named list of realized covariate values to use in the 
+#'   simulation. Supported names are:
+#'   \itemize{
+#'     \item \code{uc1}: numeric vector of length \code{n_pt} giving the unit-level covariate.
+#'     \item \code{ec1}: numeric vector of length \code{n_pt * n_season * n_rep}
+#'       giving the event-level covariate (ignored if \code{rep_constant = TRUE}).
+#'   }
+#'   If \code{NULL} (the default), covariates are simulated internally from 
+#'   standard normal distributions.
+#'   \code{ec1} is assumed not to vary by species.
 #' @param seed random seed. NULL uses (and updates) the existing RNG state. Other values
 #'   do not update the global RNG state.
 #' @param ragged_rep logical: create data with variable (TRUE) or constant 
@@ -80,7 +119,7 @@ simulate_flocker_data <- function(
   )
   assertthat::assert_that(
     is.null(multi_init) | isTRUE(multi_init %in% c("explicit", "equilibrium")),
-    msg = "multiseason must be NULL, 'explicit', or 'equilibrium'"
+    msg = "multi_init must be NULL, 'explicit', or 'equilibrium'"
   )
   assertthat::assert_that(is_one_logical(augmented), msg = "augmented must be a single logical")
   assertthat::assert_that(is_one_logical(rep_constant), msg = "rep_constant must be a single logical")
@@ -160,6 +199,72 @@ sfd <- function(
   # If coefs not specified directly, simulate from the prior (either passed via
   # params or otherwise the default prior).
   n_coef <- length(coef_names)
+
+  np <- names(params)
+
+  # Disallow unknown params keys (they don't affect simulation)
+  allowed_keys <- c("coefs", "coef_means", "Sigma")
+  unknown_keys <- setdiff(np, allowed_keys)
+  if (length(unknown_keys) > 0) {
+    stop(
+      "simulate_flocker_data(): 'params' contains unsupported names that are ignored: ",
+      paste(unknown_keys, collapse = ", "),
+      call. = FALSE
+    )
+  }
+
+  if ("coefs" %in% np) {
+    if ("coef_means" %in% np || "Sigma" %in% np) {
+      stop(
+        "simulate_flocker_data(): if params$coefs is provided, do not also provide ",
+        "params$coef_means or params$Sigma (they would be ignored).",
+        call. = FALSE
+      )
+    }
+
+    assertthat::assert_that(
+      is.data.frame(params$coefs),
+      msg = "params$coefs, if provided, must be a data.frame"
+    )
+    assertthat::assert_that(
+      nrow(params$coefs) == n_sp,
+      msg = "params$coefs, if provided, must have n_sp rows (one per species)"
+    )
+    assertthat::assert_that(
+      all(vapply(params$coefs, is.numeric, logical(1))),
+      msg = "params$coefs must have all-numeric columns"
+    )
+    assertthat::assert_that(
+      setequal(colnames(params$coefs), coef_names),
+      msg = "params$coefs must contain exactly the required coefficient names"
+    )
+
+    params$coefs <- params$coefs[, coef_names, drop = FALSE]
+  } else {
+    # If coefs absent, validate coef_means/Sigma if present
+    if ("coef_means" %in% np) {
+      assertthat::assert_that(
+        length(params$coef_means) == length(coef_names),
+        msg = "params$coef_means must have length equal to the number of coefficients"
+      )
+    }
+    if ("Sigma" %in% np) {
+      assertthat::assert_that(
+        is.matrix(params$Sigma),
+        nrow(params$Sigma) == length(coef_names),
+        ncol(params$Sigma) == length(coef_names),
+        msg = "params$Sigma must be an n_coef x n_coef matrix"
+      )
+      assertthat::assert_that(
+        isTRUE(all.equal(params$Sigma, t(params$Sigma))),
+        msg = "params$Sigma must be symmetric"
+      )
+
+      # could add a check for PSD here, but for now we'll rely on errors from MASS::mvnorm
+    }
+  }
+
+
   np <- names(params)
   if (!("coefs" %in% np)) { # coef values not specified directly
     if (!("coef_means" %in% np)) {
@@ -194,9 +299,15 @@ sfd <- function(
 
     names(params$coefs) <- coef_names
   }
+
+  # Ensure all required coefficient names are present (order irrelevant)
   assertthat::assert_that(
-    identical(names(params$coefs), coef_names)
+    setequal(names(params$coefs), coef_names),
+    msg = "params$coefs must contain exactly the required coefficient names"
   )
+
+  # Reorder columns into canonical order
+  params$coefs <- params$coefs[, coef_names, drop = FALSE]
   coefs <- params$coefs
 
   # Create covariates (currently assume that unit covariates are the same
@@ -209,6 +320,13 @@ sfd <- function(
     covariates <- list(
       uc1 = stats::rnorm(n_pt)[unit_backbone$id_point]
     )
+  } else {
+    assertthat::assert_that(is.list(covariates), msg = "covariates must be a named list")
+    assertthat::assert_that("uc1" %in% names(covariates), msg = "covariates must include uc1")
+    assertthat::assert_that(
+      is.numeric(covariates$uc1) && length(covariates$uc1) == n_pt,
+      msg = "covariates$uc1 must be numeric of length n_pt"
+    )
   }
 
   unit_backbone$uc1 <- covariates$uc1
@@ -285,7 +403,12 @@ sfd <- function(
     coefs$det_slope_unit[as.integer(visit_full$species)] * visit_full$uc1
 
   if (!rep_constant) {
-    if (!("ec1" %in% names(covariates))) {
+    if ("ec1" %in% names(covariates)) {
+      assertthat::assert_that(
+        is.numeric(covariates$ec1) && length(covariates$ec1) == n_pt * n_season * n_rep,
+        msg = "covariates$ec1 must be numeric of length n_pt * n_season * n_rep"
+      )
+    } else {
       covariates$ec1 <- stats::rnorm(n_pt*n_season*n_rep)[ # we don't allow ec1 to vary by species;
         # this allows the covariate to play nicely with augmented models.
 
@@ -344,7 +467,7 @@ sfd <- function(
 
     out <- list(obs = remove_rownames(obs), 
                 unit_covs = remove_rownames(unit_covs),
-                event_covs = lapply(event_covs, remove_rownames) # this yields an empty list if event_covs is NULL
+                event_covs = if (is.null(event_covs)) NULL else lapply(event_covs, remove_rownames)
                 )
   } else if (augmented) {
     visit_full$id_unit <- interaction(visit_full$id_point, visit_full$species)
@@ -379,7 +502,7 @@ sfd <- function(
     out <- list(
       obs = remove_rownames(abind::abind(obs_temp, along = 3)),
       unit_covs = remove_rownames(unit_covs),
-      event_covs = lapply(event_covs, remove_rownames)
+      event_covs = if (is.null(event_covs)) NULL else lapply(event_covs, remove_rownames)
     )
 
   } else if (n_season > 1) {
@@ -405,7 +528,7 @@ sfd <- function(
     }
 
     if(rep_constant){
-      event_covs <- list()
+      event_covs <- NULL
     } else {
       event_covs <- list(ec1 = remove_rownames(abind::abind(events_temp, along = 3)))
     }

inst/cache_vignettes.R

@@ -1,37 +1,119 @@
-# This script creates .Rmd vignettes that include pre-computed R output so that
-# whent they get built on CI it's a lightweight operation.
+if (!file.exists("DESCRIPTION")) {
+  stop("Run this script from the package root (where DESCRIPTION lives).")
+}
 
-# This script should be run every time there is an update to the vignettes
-# or an update to the package for which the vignettes are a desirable
-# part of the test harness (but ideally the actual test harness should be 
-# sufficient for testing, and we shouldn't rely on vignettes).
+suppressPackageStartupMessages({
+  library(withr)
+  library(knitr)
+  library(rmarkdown)
+})
 
-# Edits that touch only the text of the vignettes, and not the R computation,
-# can be made with extreme care by changing both the .Rmd.orig files and the
-# .Rmd files with computation cached.
+find_pkg_root <- function() {
+  path <- normalizePath(getwd(), winslash = "/", mustWork = TRUE)
+  repeat {
+    if (file.exists(file.path(path, "DESCRIPTION"))) return(path)
+    parent <- normalizePath(file.path(path, ".."), winslash = "/", mustWork = TRUE)
+    if (identical(parent, path)) stop("Could not find package root (no DESCRIPTION found).")
+    path <- parent
+  }
+}
 
-old_wd <- getwd()
+pkg_root <- find_pkg_root()
 
-setwd(paste0(dirname(rstudioapi::getActiveDocumentContext()$path), "/.."))
-
-knitr::knit(
-  "vignettes/augmented_models.Rmd.orig", 
-  output = "vignettes/augmented_models.Rmd"
+vigs <- list(
+  c("vignettes/augmented_models.Rmd.orig", "vignettes/augmented_models.Rmd"),
+  c("vignettes/flocker_tutorial.Rmd.orig", "vignettes/flocker_tutorial.Rmd"),
+  c("vignettes/nonlinear_models.Rmd.orig", "vignettes/nonlinear_models.Rmd"),
+  c("vignettes/articles/sbc.Rmd.orig", "vignettes/articles/sbc.Rmd"),
+  c("vignettes/articles/sbc_multi.Rmd.orig", "vignettes/articles/sbc_multi.Rmd"),
+  c("vignettes/articles/sbc_aug.Rmd.orig", "vignettes/articles/sbc_aug.Rmd")
 )
 
-knitr::knit(
-  "vignettes/flocker_tutorial.Rmd.orig", 
-  output = "vignettes/flocker_tutorial.Rmd"
-)
+with_dir(pkg_root, {
+  # Ensure relative paths in chunks resolve from package root
+  old_root <- knitr::opts_knit$get("root.dir")
+  on.exit(knitr::opts_knit$set(root.dir = old_root), add = TRUE)
+  knitr::opts_knit$set(root.dir = pkg_root)
+
+  knit_one <- function(infile, outfile, fig_path = NULL) {
+    dir.create(dirname(outfile), recursive = TRUE, showWarnings = FALSE)
+
+    if (!is.null(fig_path)) {
+      dir.create(fig_path, recursive = TRUE, showWarnings = FALSE)
+      old_fig <- knitr::opts_chunk$get("fig.path")
+      on.exit(knitr::opts_chunk$set(fig.path = old_fig), add = TRUE)
+      knitr::opts_chunk$set(fig.path = paste0(fig_path, "/"))
+    }
+
+    message("Knitting ", infile, " -> ", outfile)
+
+    env <- new.env(parent = globalenv())
+    env$params <- rmarkdown::yaml_front_matter(infile)$params
+
+    knitr::knit(infile, output = outfile, envir = env)
+  }
+
+  for(i in 1:3){
+    knit_one(vigs[[i]][[1]], vigs[[i]][[2]])
+  }
+  for(i in 4:6){
+    # SBC: force figures into man/figures/sbc_vignette
+    knit_one(vigs[[i]][[1]], vigs[[i]][[2]], fig_path = "man/figures/sbc_vignette")
+  }
 
-knitr::knit(
-  "vignettes/nonlinear_models.Rmd.orig", 
-  output = "vignettes/nonlinear_models.Rmd"
-)
+})
 
-knitr::knit(
-  "vignettes/sbc.Rmd.orig", 
-  output = "vignettes/sbc.Rmd"
-)
+# ---- Post-process cached articles to fix figure paths (Pandoc resolves relative
+# paths from the input document directory, so vignettes/articles/* need ../man/...) ----
+
+fix_article_fig_paths <- function(rmd_path,
+                                  from = "man/figures/sbc_vignette/",
+                                  to   = "../../man/figures/sbc_vignette/") {
+  # Only makes changes if the file exists and actually contains the "from" string,
+  # and avoids double-prepending if it already has "../".
+  stopifnot(length(rmd_path) == 1)
+
+  if (!file.exists(rmd_path)) {
+    warning("Skipping missing file: ", rmd_path)
+    return(invisible(FALSE))
+  }
+
+  lines <- readLines(rmd_path, warn = FALSE)
+
+  # 1) HTML: <img src="man/figures/...">
+  # Replace src="man/figures/..." but NOT src="../man/figures/..."
+  lines2 <- gsub(
+    pattern = paste0('src="', from),
+    replacement = paste0('src="', to),
+    x = lines,
+    fixed = TRUE
+  )
+
+  # 2) Markdown: ![](man/figures/...)
+  # Replace (man/figures/...) but NOT (../man/figures/...)
+  lines2 <- gsub(
+    pattern = paste0("](", from),
+    replacement = paste0("](", to),
+    x = lines2,
+    fixed = TRUE
+  )
+
+  changed <- !identical(lines, lines2)
+  if (changed) {
+    writeLines(lines2, rmd_path)
+    message("Rewrote figure paths in ", rmd_path)
+  } else {
+    message("No figure paths to rewrite in ", rmd_path)
+  }
+
+  invisible(changed)
+}
+
+# Apply to the cached article outputs (the .Rmd files, not the .orig)
+article_outfiles <- vapply(vigs[4:6], `[[`, character(1), 2)
+
+with_dir(pkg_root, {
+  for (f in article_outfiles) fix_article_fig_paths(f)
+})
 
-setwd(old_wd)
+message("Done.")