version 0.2.2

DESCRIPTION

@@ -1,13 +1,13 @@
 Package: nimbleHMC
 Title: Hamiltonian Monte Carlo and Other Gradient-Based MCMC Sampling
         Algorithms for 'nimble'
-Version: 0.2.1
-Date: 2024-02-08
+Version: 0.2.2
+Date: 2024-06-27
 Authors@R: c(person("Daniel", "Turek", role = c("aut", "cre"), email = "danielturek@gmail.com"),
 	     person("Perry", "de Valpine", role = "aut"),
 	     person("Christopher", "Paciorek", role = "aut"))
 Maintainer: Daniel Turek <danielturek@gmail.com>
-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) <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.
+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.
 Depends: R (>= 3.5.0), nimble (>= 1.0.0)
 Imports: methods
 Suggests: testthat
@@ -16,9 +16,9 @@ Copyright: See COPYRIGHTS file.
 Encoding: UTF-8
 RoxygenNote: 7.3.1
 NeedsCompilation: no
-Packaged: 2024-02-08 22:54:18 UTC; dturek
+Packaged: 2024-06-28 11:35:15 UTC; dturek
 Author: Daniel Turek [aut, cre],
   Perry de Valpine [aut],
   Christopher Paciorek [aut]
 Repository: CRAN
-Date/Publication: 2024-02-08 23:10:02 UTC
+Date/Publication: 2024-06-28 12:10:08 UTC

MD5

@@ -1,18 +1,18 @@
-2ad30bd9b808122eb1890a9d4d6d099d *DESCRIPTION
+e20f9b73538a14f64577d744afab0c0b *DESCRIPTION
 a3a1c853e2e3aa67ba36342bb4218bc7 *LICENSE
 b10b0438ac32674a8e766059023cf015 *NAMESPACE
-6dfc735f7cc26d3c79f3cb4ee3595f8b *R/HMC_configuration.R
-83a169b7c3b9b65922d98d6b29d48afb *R/HMC_samplers.R
+530f613033795b743ca9bda7180b9762 *R/HMC_configuration.R
+0108781398172a169b4be7d80ca84780 *R/HMC_samplers.R
 d2f76ad6e582ee78786dfa3208897e55 *inst/CITATION
 f539080f0f68c70212583907268dc4c6 *inst/COPYRIGHTS
 b636c1edf65a7c8d36c5c17f5f476938 *inst/NEWS.md
-e6766867432569db4172c7e2efe4b31c *man/addHMC.Rd
-ac8620f90e34294b481b72a2a8792123 *man/buildHMC.Rd
-f1297f9513d88001da4fc1a58c66a6b1 *man/configureHMC.Rd
-eb25cef9be5e4241cd8d5bf31fff635a *man/nimbleHMC.Rd
-d9155c8bdcfaaad3e6a269c3dce6f54c *man/sampler_NUTS.Rd
-f1a5501842e998b35d78f13866f5a641 *man/sampler_NUTS_classic.Rd
+9bcead8b994a33da89e7428896af79e5 *man/addHMC.Rd
+95e9c40be328c6a364f7502e76cf92a4 *man/buildHMC.Rd
+68750411d30e483766aa32982080ba1c *man/configureHMC.Rd
+657e22eefc3dfc364254618a5642114b *man/nimbleHMC.Rd
+6ca884544b8e71a66b492ef7003e43b3 *man/sampler_NUTS.Rd
+56b03e104ba8fc825203ff66bd0005bc *man/sampler_NUTS_classic.Rd
 3d275abbc4531b94441fa77447e01009 *man/stateNL_NUTS.Rd
 b204d777bea658e6c72cdb10d61e8eb2 *man/treebranchNL_NUTS.Rd
 215a55be36e7570466f0ac1499c5e86d *tests/testthat/test-HMC.R
-d6833c1389d479453ac91a6cece9e4df *tests/testthat/test-tuning.R
+5e99a8aee95cd92ed82e3279bfae07ce *tests/testthat/test-tuning.R

R/HMC_configuration.R

@@ -5,7 +5,7 @@
 #' Add a No-U-Turn (NUTS) Hamiltonian Monte Carlo (HMC) sampler to an existing nimble MCMC configuration object
 #'
 #' @param conf A nimble MCMC configuration object, as returned by `configureMCMC`.
-#' @param target A character vector of continuous-valued stochastic node names to sample.  If this argument contains any discrete-valued nodes, an error is produced and no sampler is added.
+#' @param target A character vector of continuous-valued stochastic node names to sample.  If this argument contains any discrete-valued nodes, an error is produced and no HMC sampler is added.  If this argument is omitted, then no HMC sampler is added.
 #' @param type A character string specifying the type of HMC sampler to add, either "NUTS" or "NUTS_classic".  See `help(NUTS)` or `help(NUTS_classic)` for details of each sampler.  The default sampler type is "NUTS".
 #' @param control Optional named list of control parameters to be passed as the `control` argument to the HMC sampler.  See `help(NUTS)` or `help(NUTS_classic)` for details of the control list elements accepted by each sampler.
 #' @param replace Logical argument.  If `TRUE`, any existing samplers operating on the specified nodes will be removed, prior to adding the HMC sampler.  Default value is `FALSE`.
@@ -26,7 +26,7 @@
 #'
 #' @return Invisibly returns an object of class `MCMCconf`, but this function is primary called for its side effect.
 #'
-#' @seealso \code{\link{configureHMC}} \code{\link{buildHMC}} \code{\link{configureMCMC}} \code{\link{addSampler}} \code{\link{sampler_NUTS}} \code{\link{sampler_NUTS_classic}}
+#' @seealso \code{\link{configureHMC}} \code{\link{buildHMC}} \code{\link[nimble]{configureMCMC}} \code{\link[nimble]{addSampler}} \code{\link{sampler_NUTS}} \code{\link{sampler_NUTS_classic}}
 #' 
 #' @examples
 #' code <- nimbleCode({
@@ -46,11 +46,14 @@
 #'
 #' Rmodel <- nimbleModel(code, constants, data, inits, buildDerivs = TRUE)
 #'
-#' ## create MCMC configuration object
-#' conf <- configureMCMC(Rmodel, nodes = NULL)
+#' ## create default MCMC configuration object
+#' conf <- configureMCMC(Rmodel)
 #'
-#' ## add HMC sampler operating on all stochastic model nodes
-#' addHMC(conf)
+#' ## remove default samplers operating on b0 and b1
+#' conf$removeSamplers(c("b0", "b1"))
+#'
+#' ## add an HMC sampler operating on b0 and b1
+#' addHMC(conf, target = c("b0", "b1"))
 #'
 #' Rmcmc <- buildMCMC(conf)
 #' 
@@ -98,7 +101,7 @@ addHMC <- function(conf, target = character(), type = 'NUTS', control = list(),
 #'
 #' @return An object of class `MCMCconf`.
 #'
-#' @seealso \code{\link{addHMC}} \code{\link{buildHMC}} \code{\link{configureMCMC}} \code{\link{addSampler}} \code{\link{sampler_NUTS}} \code{\link{sampler_NUTS_classic}}
+#' @seealso \code{\link{addHMC}} \code{\link{buildHMC}} \code{\link[nimble]{configureMCMC}} \code{\link[nimble]{addSampler}} \code{\link{sampler_NUTS}} \code{\link{sampler_NUTS_classic}}
 #' 
 #' @examples
 #' code <- nimbleCode({
@@ -171,7 +174,7 @@ configureHMC <- function(model, nodes = character(), type = 'NUTS', control = li
 #'
 #' @author Daniel Turek
 #' 
-#' @seealso \code{\link{addHMC}} \code{\link{configureHMC}} \code{\link{configureMCMC}} \code{\link{addSampler}} \code{\link{sampler_NUTS}} \code{\link{sampler_NUTS_classic}}
+#' @seealso \code{\link{addHMC}} \code{\link{configureHMC}} \code{\link[nimble]{configureMCMC}} \code{\link[nimble]{addSampler}} \code{\link{sampler_NUTS}} \code{\link{sampler_NUTS_classic}}
 #' 
 #' @examples
 #' code <- nimbleCode({
@@ -286,7 +289,7 @@ buildHMC <- function(model, nodes = character(), type = 'NUTS', control = list()
 #'                          summary = TRUE, WAIC = TRUE)
 #' }
 #'
-#' @seealso \code{\link{configureHMC}} \code{\link{buildHMC}} \code{\link{configureMCMC}} \code{\link{buildMCMC}} \code{\link{runMCMC}}
+#' @seealso \code{\link{configureHMC}} \code{\link{buildHMC}} \code{\link[nimble]{configureMCMC}} \code{\link[nimble]{buildMCMC}} \code{\link[nimble]{runMCMC}}
 #' 
 #' @author Daniel Turek
 #'

R/HMC_samplers.R

@@ -38,9 +38,12 @@
 ## #'     }
 ## #' })
 ## #' 
-## #' N <- 10
-## #' constants <- list(N = N, x = 1:N)
-## #' data <- list(y = 1:N)
+## #' set.seed(0)
+## #' N <- 100
+## #' x <- rnorm(N)
+## #' y <- 1 + 0.3*x + rnorm(N)
+## #' constants <- list(N = N, x = x)
+## #' data <- list(y = y)
 ## #' inits <- list(b0 = 1, b1 = 0.1, sigma = 1)
 ## #' 
 ## #' Rmodel <- nimbleModel(code, constants, data, inits, buildDerivs = TRUE)
@@ -127,31 +130,62 @@ hmc_checkTarget <- function(model, targetNodes, hmcType) {
     ## checks for:
     ## - target with discrete or truncated distribution
     ## - dependencies with truncated, dinterval, or dconstraint distribution
-    calcNodes <- model$getDependencies(targetNodes, stochOnly = TRUE)
-    if(any(model$isDiscrete(targetNodes)))
-        stop(paste0(hmcType, ' sampler cannot operate on discrete-valued nodes: ', paste0(targetNodes[model$isDiscrete(targetNodes)], collapse = ', ')))
-    if(any(model$isTruncated(targetNodes)))
-        stop(paste0(hmcType, ' sampler cannot operate on nodes with truncated prior distributions: ', paste0(targetNodes[model$isTruncated(targetNodes)], collapse = ', ')))
-    if(any(model$isTruncated(calcNodes)))
-        stop(paste0(hmcType, ' sampler cannot operate since these dependent nodes have truncated distributions, which do not support AD calculations: ', paste0(calcNodes[model$isTruncated(calcNodes)], collapse = ', ')))
-    if(any(model$getDistribution(calcNodes) == 'dinterval'))
-        stop(paste0(hmcType, ' sampler cannot operate since these dependent nodes have dinterval distributions, which do not support AD calculations: ', paste0(calcNodes[which(model$getDistribution(calcNodes) == 'dinterval')], collapse = ', ')))
-    if(any(model$getDistribution(calcNodes) == 'dconstraint'))
-        stop(paste0(hmcType, ' sampler cannot operate since these dependent nodes have dconstraint distributions, which do not support AD calculations: ', paste0(calcNodes[which(model$getDistribution(calcNodes) == 'dconstraint')], collapse = ', ')))
+    ##
+    targetDeclIDs_unique <- unique(model$getDeclID(targetNodes))
+    targetDeclInfo_unique <- model$getModelDef()$declInfo[targetDeclIDs_unique]
+    ##
+    ##if(any(model$isDiscrete(targetNodes)))
+    ##    stop(paste0(hmcType, ' sampler cannot operate on discrete-valued nodes: ', paste0(targetNodes[model$isDiscrete(targetNodes)], collapse = ', ')))
+    targetDists_unique <- unique(sapply(targetDeclInfo_unique, function(x) x$getDistributionName()))
+    targetDiscreteBool <- sapply(targetDists_unique, isDiscrete)
+    if(any(targetDiscreteBool)) {
+        stop(paste0(hmcType, ' sampler cannot operate on nodes with discrete-valued distributions: ', paste0(targetDists_unique[targetDiscreteBool], collapse = ', ')))
+    }
+    ##
+    ##if(any(model$isTruncated(targetNodes)))
+    ##    stop(paste0(hmcType, ' sampler cannot operate on nodes with truncated prior distributions: ', paste0(targetNodes[model$isTruncated(targetNodes)], collapse = ', ')))
+    targetTruncatedBool <- any(sapply(targetDeclInfo_unique, function(x) x$isTruncated()))
+    if(any(targetTruncatedBool)) {
+        targetExpanded <- model$expandNodeNames(targetNodes)
+        stop(paste0(hmcType, ' sampler cannot operate on nodes with truncated prior distributions: ', paste0(targetExpanded[model$isTruncated(targetExpanded)], collapse = ', ')))
+    }
+    ##
+    ##calcNodes <- model$getDependencies(targetNodes, stochOnly = TRUE)
+    depNodes <- model$getDependencies(targetNodes, self = FALSE, stochOnly = TRUE)
+    depDeclIDs_unique <- unique(model$getDeclID(depNodes))
+    depDeclInfo_unique <- model$getModelDef()$declInfo[depDeclIDs_unique]
+    ##
+    ##if(any(model$isTruncated(calcNodes)))
+    ##    stop(paste0(hmcType, ' sampler cannot operate since these dependent nodes have truncated distributions, which do not support AD calculations: ', paste0(calcNodes[model$isTruncated(calcNodes)], collapse = ', ')))
+    depTruncatedBool <- any(sapply(depDeclInfo_unique, function(x) x$isTruncated()))
+    if(any(depTruncatedBool)) {
+        stop(paste0(hmcType, ' sampler cannot operate since these dependent nodes have truncated distributions, which do not support AD calculations: ', paste0(depNodes[model$isTruncated(depNodes)], collapse = ', ')))
+    }
+    ##
+    depDists_unique <- sapply(depDeclInfo_unique, function(x) x$getDistributionName())
+    ##
+    ##if(any(model$getDistribution(calcNodes) == 'dinterval'))
+    ##    stop(paste0(hmcType, ' sampler cannot operate since these dependent nodes have dinterval distributions, which do not support AD calculations: ', paste0(calcNodes[which(model$getDistribution(calcNodes) == 'dinterval')], collapse = ', ')))
+    depIntervalBool <- (depDists_unique == 'dinterval')
+    if(any(depIntervalBool)) {
+        stop(paste0(hmcType, ' sampler cannot operate since these dependent nodes have dinterval distributions, which do not support AD calculations: ', paste0(depNodes[which(model$getDistribution(depNodes) == 'dinterval')], collapse = ', ')))
+    }
+    ##
+    ##if(any(model$getDistribution(calcNodes) == 'dconstraint'))
+    ##    stop(paste0(hmcType, ' sampler cannot operate since these dependent nodes have dconstraint distributions, which do not support AD calculations: ', paste0(calcNodes[which(model$getDistribution(calcNodes) == 'dconstraint')], collapse = ', ')))
+    depConstraintBool <- (depDists_unique == 'dconstraint')
+    if(any(depConstraintBool)) {
+        stop(paste0(hmcType, ' sampler cannot operate since these dependent nodes have dconstraint distributions, which do not support AD calculations: ', paste0(depNodes[which(model$getDistribution(depNodes) == 'dconstraint')], collapse = ', ')))
+    }
+    ##
     ## next, check for:
     ## - target with user-defined distribution (without AD support)
-    dists <- model$getDistribution(targetNodes)
+    ##
+    ##dists <- model$getDistribution(targetNodes)
+    dists <- targetDists_unique
     ADok <- rep(TRUE, length(dists))
     for(i in seq_along(dists)) {
-        ## these distributions get re-named to a nimble-version, and won't be found:
-        if(dists[i] %in% c('dweib', 'dmnorm', 'dmvt', 'dwish', 'dinvwish'))   next
-        ## find the function or this distribution:
-        nfObj <- get(dists[i], envir = parent.frame(4))    ## this took a bit of an investigation to make work
-        ## is a user-defined distribution:
-        if(!is.null(environment(nfObj)$nfMethodRCobject)) {
-            ## check for AD support:
-            ADok[i] <- !isFALSE(environment(nfObj)$nfMethodRCobject[['buildDerivs']])
-        }
+        ADok[i] <- model$getModelDef()$checkADsupportForDistribution(dists[i])
     }
     if(!all(ADok))
         stop(paste0(hmcType, ' sampler cannot operate on user-defined distributions which do not support AD calculations.  Try using buildDerivs = TRUE in the definition the distributions: ', paste0(dists[!ADok], collapse = ', ')))
@@ -276,9 +310,12 @@ hmc_setWarmup <- nimbleFunction(
 #'     }
 #' })
 #' 
-#' N <- 10
-#' constants <- list(N = N, x = 1:N)
-#' data <- list(y = 1:N)
+#' set.seed(0)
+#' N <- 100
+#' x <- rnorm(N)
+#' y <- 1 + 0.3*x + rnorm(N)
+#' constants <- list(N = N, x = x)
+#' data <- list(y = y)
 #' inits <- list(b0 = 1, b1 = 0.1, sigma = 1)
 #' 
 #' Rmodel <- nimbleModel(code, constants, data, inits, buildDerivs = TRUE)
@@ -741,9 +778,12 @@ treebranchNL_NUTS <- nimbleList(p_beg = double(1), p_end = double(1), rho = doub
 #'     }
 #' })
 #' 
-#' N <- 10
-#' constants <- list(N = N, x = 1:N)
-#' data <- list(y = 1:N)
+#' set.seed(0)
+#' N <- 100
+#' x <- rnorm(N)
+#' y <- 1 + 0.3*x + rnorm(N)
+#' constants <- list(N = N, x = x)
+#' data <- list(y = y)
 #' inits <- list(b0 = 1, b1 = 0.1, sigma = 1)
 #' 
 #' Rmodel <- nimbleModel(code, constants, data, inits, buildDerivs = TRUE)
@@ -890,7 +930,6 @@ sampler_NUTS <- nimbleFunction(
         divergent <<- FALSE
         branch <- treebranchNL$new()
         done <- FALSE
-        old_stepsize <- epsilon
         old_M <- M
         while((depth < maxTreeDepth) & (!done)) {
             checkInterrupt()
@@ -1140,7 +1179,6 @@ sampler_NUTS <- nimbleFunction(
         adapt_stepsize = function(adapt_stat = double()) {
             ## following Stan code, this is the same as what we have from Hoffman and Gelman, but with adapt_stat instead of a/na
             if(adapt_stat > 1)   adapt_stat <- 1
-            old_epsilon <- epsilon
             stepsizeCounter <<- stepsizeCounter + 1
             eta <- 1/(stepsizeCounter + t0)
             Hbar <<- (1-eta) * Hbar + eta * (delta - adapt_stat)          ## s_bar in Stan