updated nimbleHMC to include userEnv argument

nimbleHMC/DESCRIPTION

@@ -1,13 +1,13 @@
 Package: nimbleHMC
 Title: Hamiltonian Monte Carlo and Other Gradient-Based MCMC Sampling Algorithms for 'nimble'
-Version: 0.2.3
-Date: 2024-12-18
+Version: 0.2.4
+Date: 2025-12-15
 Authors@R: c(person("Daniel", "Turek", role = c("aut", "cre"), email = "[email protected]"),
 	     person("Perry", "de Valpine", role = "aut"),
 	     person("Christopher", "Paciorek", role = "aut"))
 Maintainer: Daniel Turek <[email protected]>
 Description: Provides gradient-based MCMC sampling algorithms for use with the MCMC engine provided by the 'nimble' package.  This includes two versions of Hamiltonian Monte Carlo (HMC) No-U-Turn (NUTS) sampling, and (under development) Langevin samplers.  The `NUTS_classic` sampler implements the original HMC-NUTS algorithm as described in Hoffman and Gelman (2014) <doi:10.48550/arXiv.1111.4246>.  The `NUTS` sampler is a modern version of HMC-NUTS sampling matching the HMC sampler available in version 2.32.2 of Stan (Stan Development Team, 2023). In addition, convenience functions are provided for generating and modifying MCMC configuration objects which employ HMC sampling. Functionality of the 'nimbleHMC' package is described further in Turek, et al (2024) <doi: 10.21105/joss.06745>.
-Depends: R (>= 3.5.0), nimble (>= 1.0.0)
+Depends: R (>= 3.5.0), nimble (>= 1.4.0)
 Imports: methods
 Suggests: testthat
 License: BSD_3_clause + file LICENSE | GPL (>= 2)

nimbleHMC/R/HMC_configuration.R

@@ -256,9 +256,11 @@ buildHMC <- function(model, nodes = character(), type = 'NUTS', control = list()
 #' @param samplesAsCodaMCMC Logical argument.  If \code{TRUE}, then a \code{coda} \code{mcmc} object is returned instead of an R matrix of samples, or when \code{nchains > 1} a \code{coda} \code{mcmc.list} object is returned containing \code{nchains} \code{mcmc} objects.  This argument is only used when \code{samples} is \code{TRUE}.  Default value is \code{FALSE}.  See details.
 #' 
 #' @param summary Logical argument.  When \code{TRUE}, summary statistics for the posterior samples of each parameter are also returned, for each MCMC chain.  This may be returned in addition to the posterior samples themselves.  Default value is \code{FALSE}.  See details.
-#'z
+#' 
 #' @param WAIC Logical argument.  When \code{TRUE}, the WAIC (Watanabe, 2010) of the model is calculated and returned.  If multiple chains are run, then a single WAIC value is calculated using the posterior samples from all chains.  Default value is \code{FALSE}. Note that the version of WAIC used is the default WAIC conditional on random effects/latent states and without any grouping of data nodes. See \code{help(waic)} for more details.
 #' 
+#' @param userEnv Environment in which if-then-else statements in model code will be evaluated if needed information not found in \code{constants}; intended primarily for internal use only.
+#' 
 #' @return A list is returned with named elements depending on the arguments, unless only one among samples, summary, and WAIC are requested, in which case only that element is returned.  These elements may include \code{samples}, \code{summary}, and \code{WAIC}.  When \code{nchains = 1}, posterior samples are returned as a single matrix, and summary statistics as a single matrix.  When \code{nchains > 1}, posterior samples are returned as a list of matrices, one matrix for each chain, and summary statistics are returned as a list containing \code{nchains+1} matrices: one matrix corresponding to each chain, and the final element providing a summary of all chains, combined.  If \code{samplesAsCodaMCMC} is \code{TRUE}, then posterior samples are provided as \code{coda} \code{mcmc} and \code{mcmc.list} objects.  When \code{WAIC} is \code{TRUE}, a WAIC summary object is returned.
 #'
 #' @details
@@ -320,11 +322,12 @@ nimbleHMC <- function(code,
                       samples = TRUE,
                       samplesAsCodaMCMC = FALSE,
                       summary = FALSE,
-                      WAIC = FALSE) {
+                      WAIC = FALSE,
+                      userEnv = parent.frame()) {
     if(missing(code) && missing(model)) stop('must provide either code or model argument')
     if(!samples && !summary && !WAIC) stop('no output specified, use samples = TRUE, summary = TRUE, or WAIC = TRUE')
     if(!missing(code) && inherits(code, 'modelBaseClass')) model <- code   ## let's handle it, if model object is provided as un-named first argument
-    Rmodel <- mcmc_createModelObject(model, inits, nchains, setSeed, code, constants, data, dimensions, check, buildDerivs = TRUE)
+    Rmodel <- mcmc_createModelObject(model, inits, nchains, setSeed, code, constants, data, dimensions, check, buildDerivs = TRUE, userEnv = userEnv)
     Rmcmc <- buildHMC(Rmodel, type = type, monitors = monitors, thin = thin, enableWAIC = WAIC, print = FALSE)
     compiledList <- compileNimble(Rmodel, Rmcmc)    ## only one compileNimble() call
     Cmcmc <- compiledList$Rmcmc

nimbleHMC/R/HMC_samplers.R

@@ -21,7 +21,6 @@
 ## #' }
 ## #' 
 ## #' @import nimble
-## #' @import methods
 ## #' 
 ## #' @aliases Langevin langevin
 ## #' 

nimbleHMC/tests/testthat/test-HMC.R

@@ -316,7 +316,7 @@ test_that('HMC on conjugate Wishart', {
     OmegaSimTrueSDs <- apply(wishRV, c(1,2), sd)
     ##
     expect_equal(as.numeric(apply(samples, 2, mean)), as.numeric(OmegaTrueMean), tol = 0.1)
-    expect_equal(as.numeric(apply(samples, 2, sd)), as.numeric(OmegaSimTrueSDs), tol = 0.02)
+    expect_equal(as.numeric(apply(samples, 2, sd)), as.numeric(OmegaSimTrueSDs), tol = 0.1)
 })
 
 test_that('HMC on LKJ', {
@@ -564,8 +564,8 @@ test_that('HMC runs with various non-differentiable constructs', {
 })
 
 test_that('HMC results for CAR match non-HMC', {
-   set.seed(1)
-   code <- nimbleCode({
+    set.seed(1)
+    code <- nimbleCode({
         S[1:N] ~ dcar_normal(adj[1:L], weights[1:L], num[1:N], tau)
         tau ~ dunif(0, 5)
         for(i in 1:N)
@@ -575,32 +575,31 @@ test_that('HMC results for CAR match non-HMC', {
             Y[i] ~ dpois(lambda[i])
         }
     })
-
+    ##
     constants <- list(N = 6,
                       num = c(1,2,2,2,2,1),
                       adj = c(2, 1,3, 2,4, 3,5, 4,6, 5),
                       weights = rep(1, 10),
                       L = 10)
     data <- list(Y = c(1,0,2,1,4,3))
     inits <- list(tau = 1, S = c(0,0,0,0,0,0))
-
+    ##
     Rmodel <- nimbleModel(code, constants, data, inits, buildDerivs = TRUE)
     conf <- configureMCMC(Rmodel, monitors = c('tau','S'))
     Rmcmc <- buildMCMC(conf)
     Cmodel <- compileNimble(Rmodel)
     Cmcmc <- compileNimble(Rmcmc, project = Rmodel)
     out <- runMCMC(Cmcmc, niter = 505000, nburnin = 5000, thin = 500)
-
+    ##
     Rmodel <- nimbleModel(code, constants, data, inits, buildDerivs = TRUE)
     conf <- configureHMC(Rmodel, monitors = c('tau','S'))
     Rmcmc <- buildMCMC(conf)
     Cmodel <- compileNimble(Rmodel)
     Cmcmc <- compileNimble(Rmcmc, project = Rmodel)
     outHMC <- runMCMC(Cmcmc, niter = 22000, nburnin=2000, thin=20)
-
+    ##
     expect_equal(apply(out[,1:6],2,mean), apply(outHMC[,1:6],2,mean), tolerance = .06)
     expect_equal(mean(out[,7]),mean(outHMC[,7]), tolerance = .15)
-
     expect_equal(apply(out,2,quantile,c(.1,.9)), apply(outHMC,2,quantile,c(.1,.9)), tolerance = 0.15)
 })
 
@@ -624,35 +623,33 @@ test_that('HMC results for mixture model match non-HMC', {
     mu <- c(0,2,4)
     data <- list(y=sample(c(rnorm(50,mu[1],.35), rnorm(250,mu[2],.35), rnorm(200,mu[3],.35)), n, replace=FALSE))
     inits <- list(k=sample(1:K,n,replace=T),mu=c(-1,2,6),mu0=1)
-
+    ##
     m <- nimbleModel(code, constants = constants, data = data,
                      inits = inits, buildDerivs = TRUE)
     conf <- configureMCMC(m, monitors=c('mu','mu0','p'))
     mcmc <- buildMCMC(conf)
     cm <- compileNimble(m)
     cmcmc <- compileNimble(mcmc)
     out <- runMCMC(cmcmc, niter=50000, nburnin=10000, thin=40)
-
+    ##
     m <- nimbleModel(code, constants = constants, data = data,
                      inits = inits, buildDerivs = TRUE)
     conf <- configureMCMC(m, nodes=c('k'), monitors=c('mu','mu0','p'))
-    conf$addSampler(c('mu0','mu','p'),'NUTS')
+    conf$addSampler(c('mu0','mu','p'), 'NUTS')
     mcmc <- buildMCMC(conf)
     cm <- compileNimble(m)
     cmcmc <- compileNimble(mcmc)
     outHMC <- runMCMC(cmcmc, niter=22000, nburnin=2000, thin=20)
-
-    ## Deal with label switching.
+    ## deal with label switching
     sorter <- function(row) {
         ord <- order(row[1:3])
         return(c(row[1:3][ord], row[4], row[5:7][ord]))
     }
     out <- t(apply(out, 1, sorter))
     outHMC <- t(apply(outHMC, 1, sorter))
-
-    expect_equal(apply(out,2,mean), apply(outHMC,2,mean), tolerance = .1)
+    ##
+    expect_equal(apply(out,2,mean), apply(outHMC,2,mean), tolerance = 0.1)
     expect_equal(apply(out,2,quantile,c(.1,.9)), apply(outHMC,2,quantile,c(.1,.9)), tolerance = 0.1)
-
 })