better error trap when constant equals 1

[paciorek]

This came up as a possible parallelization question, but the error seems like it is probably simply triggered by using T without defining it and coercion to the value 1. Try seeing if running it serially but with T not set also triggers and presuming so, trap the error.

gapCode <- nimbleCode({
  for (t in 1:mxpq){
    psi[t] <- omega
    g[t] ~ dgamma(shape=eta, scale=psi[t]/eta)
  }
  for( t in (mxpq+1):T){
    psi[t] <- omega+inprod(durpar[1:p],g[(t-p):(t-1)])+inprod(durpar[(p+1):(p+q)],psi[(t-q):(t-1)])
    g[t] ~ dgamma(shape=eta, scale=psi[t]/eta)
  }
    for (i in 1:k){
    gamdirch[i] ~ dgamma(shape=0.5,rate=1)
  }
  for (i in 1:k){
    durpar[i] <- gamdirch[i]/sum(gamdirch[1:k])
  }
  eta ~ dgamma(shape=0.001,scale=1/0.001)
  omega ~ dgamma(shape=0.001,scale=1/0.001)
})

# Gap between refresh times:
gamma_t=as.numeric(gamma_t)

T=length(gamma_t)-20
p=2
q=2
mxpq <- max(p,q)
k=p+q+1

useWAIC <- TRUE

run_MCMC_allcode <- function(seed, data, code, useWAIC = TRUE) {
  library(nimble)
 
  T=T
  p=2
  q=2
  mxpq <- max(p,q)
  k=p+q+1
 
  # Running nimble code:
  gapModel <- nimbleModel(code = code,
                               constants = list(T = T,p=p,q=q,k=k,mxpq=mxpq,alpha = rep(0.5, k))
                               ,data = data
                               ,inits = list(eta=1,omega=0.001))
 
  mcmcConf <- configureMCMC(gapModel, monitors = c("eta","omega","durpar"),
                            enableWAIC = TRUE,useConjugacy = FALSE)
  Rmcmc<-buildMCMC(mcmcConf)
  Cmodel<-compileNimble(gapModel)
  Cmcmc <- compileNimble(Rmcmc, project = gapModel)
  samplesList <- runMCMC(Cmcmc, niter = 5000, thin=1, nburnin = 0
                         ,samplesAsCodaMCMC = TRUE, WAIC = TRUE)
  return(samplesList)
 
}

chain_output <- parLapply(cl = this_cluster, X = 1:2,
                          fun = run_MCMC_allcode,
                          data = list(g=gamma_t[1:T]), code =  gapCode,
                          useWAIC = useWAIC)

stopCluster(this_cluster)

[paciorek]

Indeed this is unrelated to parallization and problem comes from having backward indexing of 3:1 because T is being set to 1. The error message does point to the problematic model code line but does not indicate the indexing issue.

I'm not sure why we don't warn users when we encounter backward indexing.

@danielturek @perrydv do either of you remember why, when we check allowNegativeIndexSequences and then if FALSE use nm_seq_noDecrease why we don't issue some sort of note/warning?

[danielturek]

@paciorek I think it was the case that some historical-types of model structures, backwards indexing can just "arise" somewhat naturally, and in those cases, it should be ignored... that was what BUGS did. So, we followed that model (of silently ignoring it), but (per usual) added an option for different behavior, if that's what a user wants (in this case, processing the backwards indexing).

Given all of that. Perhaps one path is to introduce a warning when we ignore such indexing, but perhaps that warning itself could be suppress-able. Or, maybe controlled by our verbose system option?

[perrydv]

A related issue came up recently on nimble-users here.

(The message subject was not really relevant.)

The issue was whether if nimbleModel sees "for(i in a:b)", where a:b is 1:0, does it result in no declarations? Evidently in BUGS/JAGS that is the case, and so in nested loops, e.g. "for(i in a:b[j])", one can skip some values of j for the inner loop. However we determined nimble is not supporting that. If you have your hands on these corner cases of for loop indices, is this is also relevant to decide about and handle?

My thoughts echo @danielturek 's. It does not seem like something is broken here. Given that it is a declarative language, index order should not matter. I can see some cases where it might be convenient to allow negative orders, but they seem rare. More often it seems it would be unintentional.

[paciorek]

Ok, I will add a warning when we notice backward indexing and it is being ignored (which is the default).

[paciorek]

Actually one other note here is that we actually get errors with backwards indexing. (And the exact error depends on other stuff in the model.) So I'm confused as to what use cases we might have checked when we were setting this up.

> code <- nimbleCode({
    for(i in 3:1)
        y[i] ~ dnorm(0,1)
    z ~ dnorm(0,1)
})
m <- nimbleModel(code)
Defining model
Error in genExpandedNodeAndParentNames3(debug = debug) : 
  confused assigning unrolledIndicesMatrixRow in case with no unrolledRows by numUnrolledNodes != 1 for code y[i] ~ dnorm(mean = 0, sd = 1, lower_ = -Inf, upper_ = Inf, .tau = 1, .var = 1)
> code <- nimbleCode({
    for(i in 3:1)
        y[i] ~ dnorm(0,1)
  #  z ~ dnorm(0,1)
})
m <- nimbleModel(code)
Defining model
Error: There are multiple definitions for node(s): NA.
In addition: Warning message:
In origVertexID_2_contigID[contigID_2_origVertexID] <- 1:length(contigID_2_origVertexID) :
  number of items to replace is not a multiple of replacement length