nimbleHMC version 0.2.0 released, providing improved HMC performance

[Chris Paciorek]

All,

We’ve released version 0.2.0 of nimbleHMC, which includes a new default NUTS sampler inspired by Stan’s implementation of NUTS. It also provides an updated version of our previous NUTS sampler (which is based on the original Hoffman and Gelman paper, and is now called the ‘NUTS_classic’ sampler in NIMBLE) that fixes performance issues in version 0.1.1.

As we discussed when first releasing HMC for NIMBLE, nimbleHMC provides Hamiltonian Monte Carlo samplers for use with NIMBLE, in particular NUTS samplers. NIMBLE’s HMC samplers can be flexibly assigned to a subset of model parameters, allowing users to consider various sampling configurations. For example, this allows one to use HMC for the continuous parameters in a model with discrete parameters/latent states.

-Chris (for the NIMBLE development team)

[Adrian Monroe]

Thanks for this, Chris. I had a question about setting HMC samplers to a subset of model parameters. Do you have an example for how to do that? It seems that when I set buildDerivs = TRUE when running nimbleModel, it does so for the entire model and throws an error for distributions that do not support derivatives (like dconstraint) or for dynamic indexing. If I set buildDerivs = FALSE, can I later build derivatives only for the subset of model parameters? Or am I misunderstanding something?

Thank you for your time,

Adrian

--
You received this message because you are subscribed to the Google Groups "nimble-announce" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-announ...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-announce/CA%2B7Ts_qzLPhkcq8MHS5K04qa36EKB7WtvHHn8WP-7T7Lr_c%2BKw%40mail.gmail.com.

[Qing Zhao]

Dear Chris & team,

Congratulations for the new update!! One question, is there a block sampler for this NUTS sampler, or maybe NUTS is naturally a block sampler (e.g., updating intercept and slopes simultaneously)?

Best,

Qing

On Wed, Sep 20, 2023 at 10:42 AM Chris Paciorek <paci...@stat.berkeley.edu> wrote:

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CA%2B7Ts_qzLPhkcq8MHS5K04qa36EKB7WtvHHn8WP-7T7Lr_c%2BKw%40mail.gmail.com.

[Chris Paciorek]

Hi Qing,

NUTS is naturally a block sampler. In most cases, people use NUTS on the entire parameter vector (where the parameters are transformed to live on the whole real line).

-chris

[Qing Zhao]

Hi Chris,

Awesome. Thank you for the feedback!!

Best,

Qing

[heather...@gmail.com]

Relating back to Adrian's question - 

I was hoping I could just ignore the warnings it produces and it would figure itself out when it compiled. However, when I run my model using HMC, it produces the shared library error, which seems to be regarding the dimensions of objects in my custom functions. When I run the code with buildDerivs = FALSE, the code runs fine. 

I suspect the error occurs because the buildDerivs of my custom function fails and then later it tries to reference it, even though none of the nodes in the custom sampler are using the HMC part. 

Here's the basic setup:

 library(nimble)
  library(nimbleHMC)
  library(coda)
  nimbleOptions(MCMCusePredictiveDependenciesInCalculations = TRUE) #needed if using custom sampler
  nimbleOptions(MCMCsaveHistory = TRUE)
  nimbleOptions(buildInterfacesForCompiledNestedNimbleFunctions = TRUE)
  prepbirds <- nimbleModel(code = nimIBM, constants = nc,
                           data = nd, inits = ni, calculate = T, check = T, buildDerivs = FALSE) #goes fast w/ fake data :)
  #plenty of warnings about functions that we can ignore
  prepbirds$phi.p[,1] <- 1
  prepbirds$initializeInfo()
  #can ignore PrAtLeastOne and liftedPrAtLeastOne warnings as well
  prepbirds$calculate()
  mcmcbirds <- configureMCMC(prepbirds, monitors = pars, print = T, adaptive = TRUE)
  #mcmcbirds <- configureHMC(prepbirds, monitors = pars, print = T, adaptive = TRUE, nodes = c("beta0", "beta1", "delta0", "delta1", "age0", "age1", "omega", "p.s", "phi_HY", "alpha0", "alpha1", "age.prob", "g0", "sigma", "k", "nFledglings", "sex"), )
  mcmcbirds$removeSampler(c("s", "z", "parent", "age.pre", "k"))
  mcmcbirds$addSampler(target = c("s", "z", "parent", "age.pre"),
                       type = "myRWtrajectoryPotentialSampler", silent = F,
                       control = list(scale=SamplerScale, adaptScaleOnly = T, adaptive = TRUE,
                                      npix = nc$npix, pxw = nc$pxw, xmin = nc$xmin, ymin = nc$ymin,
                                      xmax = nc$xmax, ymax = nc$ymax,
                                      gcoords = nc$gcoords, pixels = nd$pixels, pixcol = nc$pixcol, pixrow = nc$pixrow,
                                      known.z = known.z,
                                      known.parent = known.parent,
                                      known.s = known.s, known.age = known.age,
                                      known.fledge = known.fledge, east = nd$east,
                                      north = nd$north, nocc = nc$nocc, M = nc$M, xoffset = nc$xoffset,
                                      yoffset = nc$yoffset, ntrap = nc$ntrap, nind = nc$nind,
                                      open = nc$open, trap = nc$trap, seen_t = seen_t, tune = Tune
                       ))
  mcmcbirds$addSampler(target = "k[1]", type = "RW", control = list(adapative = TRUE, scale = .05))
  mcmcbirds$addSampler(target = "k[2]", type = "RW", control = list(adapative = TRUE, scale = .025))
  mcmcbirds$getUnsampledNodes() #make sure this is 0
  birdsMCMC <- buildMCMC(mcmcbirds)
 
  Cmodel <- compileNimble(prepbirds) #takes a long time (sometimes)
  Compbirds <- compileNimble(birdsMCMC, project = prepbirds, showCompilerOutput = TRUE)

The error it produces is:

Error: Failed to create the shared library. Run 'printErrors()' to see the compilation errors.
> printErrors()
P_16_nimIBM_MID_7_nfCode.cpp:4952:74: error: no matching function for call to 'rcFun_R_GlobalEnv73'
(**ADproxyModel_logProb_s)((ARG1_INDEXEDNODEINFO__.info[0]) - 1, 0, 0) = rcFun_R_GlobalEnv73(Interm_4896, Interm_4897, 64, 0.125, 0, 1, 0, 1, Interm_4898, Interm_4899, 0, 0, 1);
                                                                         ^~~~~~~~~~~~~~~~~~~
P_16_nimIBM_MID_7_nfCode.cpp:4726:9: note: candidate function not viable: no known conversion from 'NimArr<1, TYPE_>' (aka 'NimArr<1, AD<double>>') to 'NimArr<1, double> &' for 1st argument
double  rcFun_R_GlobalEnv73 ( NimArr<1, double> & ARG1_x_, NimArr<1, double> & ARG2_pi_, double ARG3_npix_, double ARG4_pxw_, double ARG5_xmin_, double ARG6_xmax_, double ARG7_ymin_, double ARG8_ymax_, NimArr<2, double> & ARG9_pixMat_, NimArr<2, double> & ARG10_gcoords_, double ARG11_xoffset_, double ARG12_yoffset_, int ARG13_log_ )  {
        ^
P_16_nimIBM_MID_7_nfCode.cpp:5977:65: error: no matching function for call to 'rcFun_R_GlobalEnv71'
(**ADproxyModel_pix)((ARG1_INDEXEDNODEINFO__.info[0]) - 1, 0) = rcFun_R_GlobalEnv71(Interm_4931, 0.125, Interm_4932, 0, 0);
                                                                ^~~~~~~~~~~~~~~~~~~
P_16_nimIBM_MID_7_nfCode.cpp:5899:9: note: candidate function not viable: no known conversion from 'NimArr<1, TYPE_>' (aka 'NimArr<1, AD<double>>') to 'NimArr<1, double> &' for 1st argument
double  rcFun_R_GlobalEnv71 ( NimArr<1, double> & ARG1_s_, double ARG2_pxw_, NimArr<2, double> & ARG3_pixMat_, double ARG4_xoffset_, double ARG5_yoffset_ )  {
        ^
P_16_nimIBM_MID_7_nfCode.cpp:7386:108: error: no matching function for call to object of type 'NimArr<2, CppAD::AD<double>>'
 (**ADproxyModel_phi_dot_p)((ARG1_INDEXEDNODEINFO__.info[1]) - 1, (ARG1_INDEXEDNODEINFO__.info[0]) - 1) = ((**ADproxyModel_phi_dot_loc)(((**ADproxyModel_pix)((ARG1_INDEXEDNODEINFO__.info[1]) - 1, (ARG1_INDEXEDNODEINFO__.info[2]) - 1)) - 1, (ARG1_INDEXEDNODEINFO__.info[2]) - 1) * (**ADproxyModel_phi_dot_age)[((**ADproxyModel_age)((ARG1_INDEXEDNODEINFO__.info[1]) - 1, (ARG1_INDEXEDNODEINFO__.info[2]) - 1)) - 1]) * (**ADproxyModel_z)((ARG1_INDEXEDNODEINFO__.info[1]) - 1, (ARG1_INDEXEDNODEINFO__.info[2]) - 1) + (**ADproxyModel_a)((ARG1_INDEXEDNODEINFO__.info[1]) - 1, (ARG1_INDEXEDNODEINFO__.info[2]) - 1) * (**ADproxyModel_Ent_Prob)[(ARG1_INDEXEDNODEINFO__.info[2]) - 1];
                                                                                                           ^~~~~~~~~~~~~~~~~~~~~~~~~~~~
/Users/heather/Library/R/x86_64/4.1/library/nimble/include/nimble/NimArr.h:345:6: note: candidate function not viable: no known conversion from 'AD<double>' to 'int' for 1st argument
  T &operator()(int i, int j) const {
     ^
P_16_nimIBM_MID_7_nfCode.cpp:7386:309: error: no viable overloaded operator[] for type 'NimArr<1, CppAD::AD<double>>'
 (**ADproxyModel_phi_dot_p)((ARG1_INDEXEDNODEINFO__.info[1]) - 1, (ARG1_INDEXEDNODEINFO__.info[0]) - 1) = ((**ADproxyModel_phi_dot_loc)(((**ADproxyModel_pix)((ARG1_INDEXEDNODEINFO__.info[1]) - 1, (ARG1_INDEXEDNODEINFO__.info[2]) - 1)) - 1, (ARG1_INDEXEDNODEINFO__.info[2]) - 1) * (**ADproxyModel_phi_dot_age)[((**ADproxyModel_age)((ARG1_INDEXEDNODEINFO__.info[1]) - 1, (ARG1_INDEXEDNODEINFO__.info[2]) - 1)) - 1]) * (**ADproxyModel_z)((ARG1_INDEXEDNODEINFO__.info[1]) - 1, (ARG1_INDEXEDNODEINFO__.info[2]) - 1) + (**ADproxyModel_a)((ARG1_INDEXEDNODEINFO__.info[1]) - 1, (ARG1_INDEXEDNODEINFO__.info[2]) - 1) * (**ADproxyModel_Ent_Prob)[(ARG1_INDEXEDNODEINFO__.info[2]) - 1];
                                                                                                                                                                                                                                                                                        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/Users/heather/Library/R/x86_64/4.1/library/nimble/include/nimble/NimArrBase.h:127:6: note: candidate function not viable: no known conversion from 'AD<double>' to 'int' for 1st argument
  T &operator[](int i) const { return ((*vPtr)[offset + i * stride1]); }
     ^
P_16_nimIBM_MID_7_nfCode.cpp:8279:41: error: cannot convert 'AD<double>' to 'int' without a conversion operator
 Interm_5009.setMap((**ADproxyModel_s), static_cast<int>(((**ADproxyModel_parent)((ARG1_INDEXEDNODEINFO__.info[1]) - 1, (ARG1_INDEXEDNODEINFO__.info[0]) - 1) - 1) * (**ADproxyModel_s).strides()[0] + (ARG1_INDEXEDNODEINFO__.info[2] - 1) * (**ADproxyModel_s).strides()[2]), (**ADproxyModel_s).strides()[1], 2);
                                        ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
P_16_nimIBM_MID_7_nfCode.cpp:8285:110: error: no matching function for call to 'rcFun_R_GlobalEnv75'
 (**ADproxyModel_logProb_s)((ARG1_INDEXEDNODEINFO__.info[1]) - 1, 0, (ARG1_INDEXEDNODEINFO__.info[0]) - 1) = rcFun_R_GlobalEnv75(Interm_5002, Interm_5003, Interm_5004, Interm_5005, Interm_5006, Interm_5007, Interm_5008, Interm_5009, Interm_5010, Interm_5011, Interm_5012, Interm_5013, Interm_5014, 0, 1, 0, 1, 1);
                                                                                                             ^~~~~~~~~~~~~~~~~~~
P_16_nimIBM_MID_7_nfCode.cpp:7880:9: note: candidate function not viable: no known conversion from 'NimArr<1, TYPE_>' (aka 'NimArr<1, AD<double>>') to 'NimArr<1, double> &' for 1st argument
double  rcFun_R_GlobalEnv75 ( NimArr<1, double> & ARG1_x_, NimArr<1, double> & ARG2_kappa_, double ARG3_id_, double ARG4_t_, NimArr<1, double> & ARG5_s_dot_pre_, NimArr<1, double> & ARG6_omega_, NimArr<1, double> & ARG7_prev_dot_pix_, NimArr<1, double> & ARG8_s_dot_parent_, NimArr<1, double> & ARG9_a_, NimArr<1, double> & ARG10_z_, NimArr<2, double> & ARG11_east_, NimArr<2, double> & ARG12_north_, double ARG13_myparent_, double ARG14_xmin_, double ARG15_xmax_, double ARG16_ymin_, double ARG17_ymax_, int ARG18_log_ )  {
        ^
P_16_nimIBM_MID_7_nfCode.cpp:8663:100: error: no matching function for call to 'rcFun_R_GlobalEnv71'
(**ADproxyModel_pix)((ARG1_INDEXEDNODEINFO__.info[1]) - 1, (ARG1_INDEXEDNODEINFO__.info[0]) - 1) = rcFun_R_GlobalEnv71(Interm_5044, 0.125, Interm_5045, 0, 0);
                                                                                                   ^~~~~~~~~~~~~~~~~~~
P_16_nimIBM_MID_7_nfCode.cpp:5899:9: note: candidate function not viable: no known conversion from 'NimArr<1, TYPE_>' (aka 'NimArr<1, AD<double>>') to 'NimArr<1, double> &' for 1st argument
double  rcFun_R_GlobalEnv71 ( NimArr<1, double> & ARG1_s_, double ARG2_pxw_, NimArr<2, double> & ARG3_pixMat_, double ARG4_xoffset_, double ARG5_yoffset_ )  {
        ^
P_16_nimIBM_MID_7_nfCode.cpp:10068:15: error: no matching function for call to 'rcFun_R_GlobalEnv72'
Interm_5061 = rcFun_R_GlobalEnv72(Interm_5057, Interm_5058, Interm_5059, Interm_5060, 144);
              ^~~~~~~~~~~~~~~~~~~
P_16_nimIBM_MID_7_nfCode.cpp:9975:20: note: candidate function not viable: no known conversion from 'NimArr<1, TYPE_>' (aka 'NimArr<1, AD<double>>') to 'NimArr<1, double> &' for 1st argument
NimArr<1, double>  rcFun_R_GlobalEnv72 ( NimArr<1, double> & ARG1_s_, double ARG2_g0_, double ARG3_sigma_, NimArr<2, double> & ARG4_trap_, double ARG5_J_ )  {
                   ^
8 errors generated.
make: *** [P_16_nimIBM_MID_7_nfCode.o] Error 1

[Perry de Valpine]

It is definitely the case that error-trapping of problems with AD-related code is much more limited than of non-AD related code.

It is also correct that there is no mechanism to skip generation of AD code for one part of the model even when you know no algorithm will want to use AD for that part of the model.

I would need to see the user-defined functions or distributions and how they are called from the model.

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/d5d1b3b2-ee83-4121-ad2f-b804dce1f248n%40googlegroups.com.

[heather...@gmail.com]

Here's where I think the error (or one of them anyway) originally starts. It involves dcat, which I know is a sticking point for AD. 

dPointProcess <- nimbleFunction(
  run = function( x = double(1),pi=double(1), npix=double(0), pxw=double(0),
                  xmin=double(0), xmax = double(0), ymin = double(0), ymax = double(0),
                  pixMat=double(2), gcoords=double(2), xoffset= double(0), yoffset = double(0),
                  log = integer(0, default = 0)) {
    returnType(double(0))
   
    # assign continuous space location to grid cell
    sx <- nimRound((x[1]+ .5*pxw)/pxw-xoffset)
    sy <- nimRound((x[2]+ .5*pxw)/pxw-yoffset)
    pix <- pixMat[sx, sy]
   
    loglike <- dcat(pix, pi[1:npix], log=TRUE)
   
    if(log) return(loglike)
    else return(exp(loglike))
  }
)

#random value generator.
rPointProcess <- nimbleFunction(
  run = function(n = integer(0), pi=double(1), npix=double(0), pxw=double(0),
                 xmin=double(0), xmax = double(0), ymin = double(0), ymax = double(0),
                 pixMat=double(2), gcoords=double(2), xoffset= double(0), yoffset = double(0)) {
    returnType(double(1))
    # select grid cell
    pix <- rcat(1, pi[1:npix])
    # xy coordinates are uniform within grid cell
    x <- c(runif(1, max(xmin,gcoords[pix,1]-.5*pxw), min(xmax, gcoords[pix,1]+.5*pxw)), # x-location
           runif(1, max(ymin,gcoords[pix,2]-.5*pxw), min(ymax, gcoords[pix,2]+.5*pxw))) # y-location
   
    return(x[1:2])
  })

In the model, this gets called around here:

for(q in 1:npix){
    lambda[q] <- exp(beta0 + beta1*cov[q,1])*(1/npix) #relative probability of any particular location
    pi[q] <- lambda[q]/Lambda[1]
  } #end q

  Lambda[1] <- sum(lambda[1:npix])
 
 
  psi <- Lambda[1]/M #initial alive/dead state in time t = 1
 
    s[i, 1:2, 1] ~ dPointProcess(pi=pi[1:npix], ## expected distribution prob of N individuals at occasion 1
                                 npix=npix,
                                 xmin=xmin, xmax = xmax, ymin = ymin, ymax = ymax,
                                 pxw = pxw, xoffset = xoffset, yoffset = yoffset,
                                 pixMat = pixels[1:pixrow, 1:pixcol],
                                 gcoords=gcoords[1:npix, 1:2])

[Perry de Valpine]

Hi Heather,

Before I try diving in any further:

A nimbleFunction that is to be used with AD algorithms (in a model with buildDerivs=TRUE) needs to have buildDerivs=TRUE (this is also an argument to nimbleFunction).

And, dcat is not supported for derivatives, especially if the index will change. Somebody, I think, it will be supported, but it isn't currently.

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/3c2734de-2f62-43d5-aef9-fc55d8f0ebd0n%40googlegroups.com.

[Heather Gaya]

Absolutely, I understand that. However, I don't want this function to use AD algorithms. It's just also in the same model with other variables that I would like to use AD algorithms :) 

[Perry de Valpine]

Yes, that makes sense.

We realize it's not at all easy right now. We're going to try to come up with some workarounds where possible and find some better solutions in a future release.

Perry

[Heather Gaya]

Would it be accurate to say that, for the moment, Nimble allows for HMC only for models that:

•  are all continuous or 
•  are a mix of continuous and discrete, but don't involve any custom functions that include discrete variables? 

-Heather 

P.S. Thank you for working towards incorporating HMC in models! I imagine that will be a huge time-saver once all the kinks are worked out. Hope I haven't sounded ungracious in any of my previous emails - just excited to try out the new features :) 

[Perry de Valpine]

I'm not sure your second bullet point characterizes the problem.

The issue is discrete distributions for latent states, like dcat and dconstraint, for which there is no way around that the code won't compile (and/or won't give correct results) even if you don't want to use derivatives involving those nodes. I believe that is a problem whether the distribution "d" function is used directly in a model or in a user-defined function or distribution.

Building AD support was a huge change in nimble. So, we had our hands full. We definitely want to make it easier to do what you (and others) are asking about. It will take time. You definitely didn't sound ungracious. It's a bit harder than we realized it would be. Thanks for the thoughtful note.

Perry

[Qing Zhao]

Dear Team,

I'm trying to use the NUTS sampler, but got the following error. I must have missed something?

> mcmcConf$addSampler(target = c('beta_lambda0'), type = "NUTS")

Error in mcmcConf$addSampler(target = c("beta_lambda0"), type = "NUTS") : 

  cannot find sampler type 'NUTS'

Best,

Qing

On Wed, Sep 20, 2023 at 10:42 AM Chris Paciorek <paci...@stat.berkeley.edu> wrote:

--

[Daniel Turek]

Hmm, if you have the latest version of the package, that should work fine.  Are you certainly using version 0.2.0 of nimbleHMC?  If so, are you able to send a reproducible example giving this error?

Thanks,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CAE-HdUGvfypqt5vTftpsUgAJGDoNQJHjb%3DUORL%2BeDdbwZk-b8Q%40mail.gmail.com.

[Abby Keller]

Hi,

I'm trying to use HMC NUTS sampling on just one node, since one of my latent states is discrete. I'm a bit confused about how nimble and nimbleHMC interact if you are only adding a NUTS sampler to one node.

For example, I create the model with nimbleModel(), configure the MCMC with configureMCMC(), and add a NUTS block sampler for the 'trap_pmax_log' node:

#remove RW block sampler

mcmcConf_myModel$removeSamplers('trap_pmax_log')

#add NUTS block sampler
for(i in 1:nsite_sim){
  for(t in 1:ntime_sim){
    mcmcConf_myModel$addSampler(paste0('trap_pmax_log[',t,', 1:',nobs_sim,', ',i,']'),type='NUTS')
  }
}

When I build the MCMC with buildMCMC(mcmcConf_myModel), the line of code runs for ~20 min and then crashes my R session.

I'm wondering if I'm adding the NUTS sampler correctly? Should I be using functions from the nimbleHMC package?

I've attached my model -- the code for creating the model/building the MCMC starts at line 562.

Thanks so much!

Abby

[Qing Zhao]

Hi Daniel & All,

It has been a while. I finally got time to try to solve this problem, but didn't succeed. 

I first removed the "nimble" folder from the R library folder to uninstall it. I then installed nimble from R. Not sure if that's the correct way to reinstall nimble, but after that, it seems I now have nimble version 1.0.1. Under this version, when I ran the code below, I got the following error message after the line "mcmcConf$addSampler(target = c('beta'), type = "NUTS")": 

Error in mcmcConf$addSampler(target = c("beta"), type = "NUTS") : 

  cannot find sampler type 'NUTS'

Please let me know if I did something wrong. Thank you.

Best regards, 

Qing

library(boot)

library(nimble)

nyear <- 11 # number of years

nindv <- 100 # number of individuals

beta  <- c( 0.4,  0.8) # intercept and slopes for apparent survival probability

alpha <- c(-1.5, -0.5) # intercept and slopes for capture probability

x <- rnorm(nyear-1, 0, 1) # environmental covariates

# Simulate data

### Process

phi <- inv.logit(cbind(1,x) %*% beta) # apparent survival probability

f <- rep(1:(nyear-1), each=nindv/(nyear-1)) # year of first capture

z <- matrix(0, nindv, nyear) # individual survival history

for (v in 1:nindv) {

  z[v,f[v]] <- 1

  for (t in (f[v]+1):nyear) {

    z[v,t] <- rbinom(1, 1, z[v,t-1] * phi[t-1])

  } # t

} # v

### Capture

w <- rnorm(nyear-1, 0, 1) # capture covariates

pcap <- inv.logit(cbind(1,w) %*% alpha) # capture probability

c <- matrix(0, nindv, nyear) # individual capture history

for (v in 1:nindv) {

  c[v,f[v]] <- z[v,f[v]]

  for (t in (f[v]+1):nyear) {

    c[v,t] <- rbinom(1, 1, z[v,t] * pcap[t-1])

  } # t

} # v

#========================

# Define model in Nimble

#========================

code <- nimbleCode({

      

  # Priors

  for (k in 1:2) {

    beta[k] ~ dnorm(0, sd=10)

    alpha[k] ~ dnorm(0, sd=10)

  } # k

  # Process model

  for (t in 1:(nyear-1)) {

    phi[t] <- ilogit(beta[1] + beta[2] * x[t])

  } # t

  for (v in 1:nindv) {

    for (t in (f[v]+1):nyear) {

      z[v,t] ~ dbern(z[v,t-1] * phi[t-1])

    } # t

  } # v

  # Observation model

  for (t in 1:(nyear-1)) {

    pcap[t] <- ilogit(alpha[1] + alpha[2] * w[t])

  } # t

  for (v in 1:nindv) {

    for (t in (f[v]+1):nyear) {

      c[v,t] ~ dbern(z[v,t] * pcap[t-1])

    } # t

  } # v

}) # nimbleCode

#============

# Run Nimble

#============

zd <- c

for (v in 1:nindv) {

  if (sum(c[v,]) > 1) {

    zd[v, min(which(c[v,]==1)):max(which(c[v,]==1))] <- 1

  }

  if (max(which(c[v,]==1)) < nyear) {

    zd[v, (max(which(c[v,]==1))+1):nyear] <- NA

  }

} # v

zi <- c

for (v in 1:nindv) {

  zi[v,] <- rev(sort(sample(c(0,1), nyear, replace=T)))

} # v

zi[which(!(is.na(zd)))] <- NA

# Data

constants <- list(

  nindv=nindv, nyear=nyear, 

  f=f, x=x, w=w

)

data <- list(

  c=c, z=zd

)

# Initial values

inits <- list(

  beta=c(0,0), alpha=c(0,0), z=zi

)

model <- nimbleModel(code, constants=constants, data=data, inits=inits)

mcmcConf <- configureMCMC(model)

mcmcConf$printSamplers()

mcmcConf$removeSamplers(c('beta', 'alpha'))

mcmcConf$printSamplers()

mcmcConf$addSampler(target = c('beta'), type = "NUTS")

mcmcConf$addSampler(target = c('alpha'), type = "RW_block")

mcmcConf$printSamplers()

mcmc <- buildMCMC(mcmcConf)

compiled <- compileNimble(model, mcmc)

nmcmc <- 5000

fit <- runMCMC(compiled$mcmc, nchains = 3, niter = nmcmc, nburnin = 0)

[Daniel Turek]

Qing, I'm guessing that you just need to load the package nimbleHMC.  Note that the HMC ("NUTS") sampler is *not* contained in the core nimble package, but rather in the separate nimbleHMC package.  Try adding library(nimbleHMC) at the beginning of your script, when you load the packages.

Cheers,

Daniel

[Qing Zhao]

Hi Daniel,

Thank you for the response. I got the following error message when trying to install the nimbleHMC packages. Do you know how to solve this? Thank you.

Best,

Qing

> install.packages('nimbleHMC')

--- Please select a CRAN mirror for use in this session ---

installing the source package ‘nimbleHMC’

trying URL 'https://cloud.r-project.org/src/contrib/nimbleHMC_0.2.0.tar.gz'

Content type 'application/x-gzip' length 39112 bytes (38 KB)

downloaded 38 KB

* installing *source* package 'nimbleHMC' ...

** package 'nimbleHMC' successfully unpacked and MD5 sums checked

** using staged installation

** R

** inst

** byte-compile and prepare package for lazy loading

** help

*** installing help indices

  converting help for package 'nimbleHMC'

    finding HTML links ... done

    addHMC                                  html  

    finding level-2 HTML links ... done

    buildHMC                                html  

    configureHMC                            html  

    nimbleHMC                               html  

    sampler_NUTS                            html  

REDIRECT:topic   Previous alias or file overwritten by alias: C:/Users/qing.zhao/Documents/R/R-4.1.2/library/00LOCK-nimbleHMC/00new/nimbleHMC/help/nuts.html

REDIRECT:topic   Previous alias or file overwritten by alias: C:/Users/qing.zhao/Documents/R/R-4.1.2/library/00LOCK-nimbleHMC/00new/nimbleHMC/help/hmc.html

    sampler_NUTS_classic                    html  

REDIRECT:topic   Previous alias or file overwritten by alias: C:/Users/qing.zhao/Documents/R/R-4.1.2/library/00LOCK-nimbleHMC/00new/nimbleHMC/help/nuts-classic.html

REDIRECT:topic   Previous alias or file overwritten by alias: C:/Users/qing.zhao/Documents/R/R-4.1.2/library/00LOCK-nimbleHMC/00new/nimbleHMC/help/nuts_classic.html

    stateNL_NUTS                            html  

    treebranchNL_NUTS                       html  

** building package indices

** testing if installed package can be loaded from temporary location

*** arch - i386

Error : package 'nimble' is not installed for 'arch = i386'

Error: loading failed

Execution halted

*** arch - x64

ERROR: loading failed for 'i386'

* removing 'C:/Users/qing.zhao/Documents/R/R-4.1.2/library/nimbleHMC'

The downloaded source packages are in

        ‘C:\Users\qing.zhao\AppData\Local\Temp\Rtmp2N79FZ\downloaded_packages’

Warning message:

In install.packages("nimbleHMC") :

  installation of package ‘nimbleHMC’ had non-zero exit status

[Chris Paciorek]

Hi Qing,

I'm not sure exactly what is causing the installation problem (from the message it sounds like some inconsistency between `nimble` and `nimbleHMC` related to i386 vs x64). In any event, it looks like the installation is trying to install from the nimble source package (i.e., the original code files, which need  to be compiled), rather than from a binary package, in which CRAN has already compiled the files. Using the binary package is much less error-prone, so I think the solution here is to try to install the binary package.

I believe that the installation is using the source package because you are using a somewhat old version of R (4.1.2). I suggest you upgrade to the current version of R and try again and it should install the binary package.

-Chris

[Qing Zhao]

Hi Chris,

Thank you for the response. I will give it a try.

Best,

Qing

[Chris Paciorek]

Heather and others following this thread,

Version 1.1.0, just released last week, should substantially smooth out the process of using nimbleHMC for models with discrete distributions. In particular, automatic differentiation (AD) can be used in a wider range of  models, including models with stochastic indexing, discrete latent states, and CAR distributions. Support for AD for these models means that HMC sampling and Laplace approximation can be used, although of course HMC samplers can only be assigned to nodes following continuous distributions. In addition, version 1.1.0 now allows distributions and functions (whether user-defined or built-in) that lack AD support (such as dinterval, dconstraint, and truncated distributions) to be used and compiled in AD-enabled models. The added flexibility increases the range of models in which one can use AD methods (HMC or Laplace) on some parts of a model and other samplers or methods on other parts.

We've updated this section of the manual to try to be clear about what is and is not supported by AD and implications for derivative-based methods such as HMC.

-chris

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/CABeTmqgKWVKDiQrAgSRbNSi6rC%2B0Su6J8uib03V62448ZUA_cg%40mail.gmail.com.

[Max Rubino]

Hi Chris,

I'm very excited to hear that Nimble now supports AD for a wider range of models. I tried testing this on a multistate model that includes dcat, but am running into issues creating the shared library when compiling. The likelihood portion of the model is below

for(i in 1:marked){
   
    z[i,f[i]]~dcat(TLM[fs[i],1:10])
   
    for(t in (f[i]+1):l[i]){
      z[i,t] ~ dcat(TPM[t-1,z[i,t-1],1:10])
      y[i,t] ~ dcat(OPM[i,z[i,t],1:10])
    }
  }

The errors are:

> printErrors()
using C++ compiler: ‘Apple clang version 14.0.0 (clang-1400.0.29.202)’
using SDK: ‘MacOSX13.1.sdk’
In file included from P_1_MEcode_MID_1.cpp:7:
In file included from ./P_1_MEcode_MID_1.h:6:
In file included from /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/nimble/include/nimble/NimArr.h:26:
In file included from /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/nimble/include/nimble/NimArrBase.h:33:
In file included from /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/nimble/include/Eigen/Core:382:
/Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/nimble/include/Eigen/src/Core/util/ReenableStupidWarnings.h:14:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
    #pragma clang diagnostic pop
                             ^
In file included from P_1_MEcode_MID_1.cpp:7:
In file included from ./P_1_MEcode_MID_1.h:11:
In file included from /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/nimble/include/nimble/nimbleCppAD.h:38:
In file included from /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/nimble/include/nimble/EigenTypedefs.h:24:
In file included from /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/nimble/include/Eigen/Dense:2:
In file included from /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/nimble/include/Eigen/LU:45:
/Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/nimble/include/Eigen/src/Core/util/ReenableStupidWarnings.h:14:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
    #pragma clang diagnostic pop
                             ^
In file included from P_1_MEcode_MID_1.cpp:7:
In file included from ./P_1_MEcode_MID_1.h:11:
In file included from /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/nimble/include/nimble/nimbleCppAD.h:38:
In file included from /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/nimble/include/nimble/EigenTypedefs.h:24:
In file included from /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/nimble/include/Eigen/Dense:3:
In file included from /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/nimble/include/Eigen/Cholesky:12:
In file included from /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/nimble/include/Eigen/Jacobi:29:
/Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/library/nimble/include/Eigen/src/Core/util/ReenableStupidWarnings.h:14:30: warning: pragma diagnostic pop could not pop, no matching push [-Wunknown-pragmas]
    #pragma clang diagnostic pop

Any insight would be appreciated! Perhaps I am not understanding what is now supported by AD?

Thanks,

Max

[Perry de Valpine]

Hi Max,

I actually don't see any errors there. Those are warnings that turn out to be harmless.  It is voluminous and alarming at first, but there's not much we can do about it. They come from the Eigen linear algebra package that we use. Eigen actually has a file in its source code called "DisableStupidWarnings.h" that uses C++ #pragma (pre-processor) statements do disable these warnings, but CRAN policy disallows disabling of C++ warnings in this way and scans code to catch violations of its policy, so we have to let the warning fly at you (if you have showCompilerOutput=TRUE). In general it is difficult to manage C++ compiler warnings because users may be running many variants of C++ compilers and configurations. That's why by default the C++ compiler output is not shown by compileNimble. If you are getting output that shows errors (typically with the word "error"), please send a reproducible example. Thanks and HTH.

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/6290b8e7-700a-4181-91f1-d7a6ac9cbe67n%40googlegroups.com.

[Max Rubino]

Hi Perry,

Thanks for the quick response. I am running into an issue when compiling, Error: Failed to create the shared library. Run 'printErrors()' to see the compilation errors.

I have attached the model code, which should be reproducible. 

Thanks so very much!

Max

[Max Rubino]

Upon further inspection, this is the error preventing creation of the shared library.

P_15_MEcode_MID_20_nfCode.cpp:7156:42: error: cannot convert 'AD<double>' to 'int' without a conversion operator
Interm_2602.setMap((**ADproxyModel_TPM), static_cast<int>((ARG1_INDEXEDNODEINFO__.info[2] - 1) * (**ADproxyModel_TPM).strides()[0] + ((**ADproxyModel_z)((ARG1_INDEXEDNODEINFO__.info[0]) - 1, (ARG1_INDEXEDNODEINFO__.info[2]) - 1) - 1) * (**ADproxyModel_TPM).strides()[1]), (**ADproxyModel_TPM).strides()[2], 10);
                                         ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
P_15_MEcode_MID_20_nfCode.cpp:7424:42: error: cannot convert 'AD<double>' to 'int' without a conversion operator
Interm_2611.setMap((**ADproxyModel_OPM), static_cast<int>((ARG1_INDEXEDNODEINFO__.info[0] - 1) * (**ADproxyModel_OPM).strides()[0] + ((**ADproxyModel_z)((ARG1_INDEXEDNODEINFO__.info[0]) - 1, (ARG1_INDEXEDNODEINFO__.info[1]) - 1) - 1) * (**ADproxyModel_OPM).strides()[1]), (**ADproxyModel_OPM).strides()[2], 10);

Max

[Perry de Valpine]

Hi Max,

Thanks for this. Yes, I reproduced this as well. 

The problem is that the stochastic indexing in the dcat with derivatives enabled is not working for this case. Let me explain a little and give you a workaround. We'll try to extend AD to cover this in the future.

Here's the explanation. Doing AD correctly with stochastic indexing is actually surprisingly tricky. That's because the AD system (using CppAD) tracks (or records) all calculations in a model (or part of a model) and then reuses them forwards and backwards for derivatives, essentially for implementing the chain rule of calculus. When there is a stochastic index, the pathway of calculations may change. So this seemingly innocuous aspect of a model really complicates any system that manages calculations as a static graph or tree structure. One way or another, all elements rather than just a single element need to be included in the graph. (This is true in basic nimble models as well (but is more thoroughly implemented than for AD) and always introduces more computation.) Some cases now work, but the one if your code doesn't. 

Here's the workaround. Stochastic indexing does work one by one. So one can write a nimbleFunction to provide a special flavor of dcat that constructs the vector TPM[t-1,z[i,t-1],1:10] (for example) one by one and then calls the dcat distribution. That's easy enough and is dmycat in the attached code. Then one wants the categorical sampler to work, but that needs tweaking to determine how many categories are available and to disable an error-trap that gets in the way. I copied the sampler_categorical source code and made these modifications. And finally one needs a little code to assign this tweaked sampler to the dmycat nodes, which I've also done.

The modified model code looks like

      z[i,t] ~ dmycat(TPM[t-1, 1:10, 1:10], z[i,t-1]) # replaces dcat(TPM[t-1,z[i,t-1],1:10])
      y[i,t] ~ dmycat(OPM[i, 1:10, 1:10], z[i, t]).     # replaced dcat(OPM[i,z[i,t],1:10])

The other changes are a bit longer and are in the attached file. The model and MCMC build and compile on my system. I didn't investigate further whether they work well or not.

Let me know how it goes.

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/fd0f1e73-d060-4ae2-85dd-49a9e6891c90n%40googlegroups.com.

[Max Rubino]

That did the trick. I appreciate the help and the explanation of what was causing the issue. 

Best,

Max