Could not evaluate loop syntax

[Anthony S]

Hey everyone,

I am running a large N-mixture model. I am actually attempting to run two different models with the same coefficients in an attempt to leverage two different types of data. Therefore, I have a lot of nested loops. 

Here is the model:

P.lat.sv.mod <- nimble::nimbleCode({

  

        psi ~ dbeta(1, 1)

  

  for(k in 1:13){

  beta.lam[k] ~ dunif(-1,1)

  sigma.lam[k] ~ dbern(psi)}

  #var.beta.lam[k] ~ dunif(-1,1)}

  beta0.lam ~ dnorm(1,.01)

 # var.beta0.lam ~ dunif(1,2)

  #Detection Priors

  for(j in 1:5){

  beta.p[j] ~ dnorm(0,.01)}

  beta0.p ~ dnorm(-4,.01)

  tau.p <- pow(sd.p, -2)

  sd.p ~ dunif(1,4) 

  

  

 for(i in 1:Id.sv){

#loop through sites

for(y in 1:n_Year.sv[i]){

#loop through years

for(s in 1:n_Season.sv[i,y]){

#loop through seasons

  

  omega.sv[i,y,s] ~ dbeta(1,1)

#state model

N.sv[i,y,s] ~ dZIP(lambda = lambda.sv[i,y,s], zeroProb = omega.sv[i,y,s])

log(lambda.sv[i,y,s]) <- beta0.lam + sigma.lam[1]*beta.lam[1]*lam.var.sv[i,y,s,1] + sigma.lam[2]*beta.lam[2]*lam.var.sv[i,y,s,2] + sigma.lam[3]*beta.lam[3]*lam.var.sv[i,y,s,3] + sigma.lam[4]*beta.lam[4]*lam.var.sv[i,y,s,4] + sigma.lam[5]*beta.lam[5]*lam.var.sv[i,y,s,5] + sigma.lam[6]*beta.lam[6]*lam.var.sv[i,y,s,6] + sigma.lam[7]*beta.lam[7]*lam.var.sv[i,y,s,7] + sigma.lam[8]*beta.lam[8]*lam.var.sv[i,y,s,8] + sigma.lam[9]*beta.lam[9]*lam.var.sv[i,y,s,9] + sigma.lam[10]*beta.lam[10]*lam.var.sv[i,y,s,10] + sigma.lam[11]*beta.lam[11]*lam.var.sv[i,y,s,11] + sigma.lam[12]*beta.lam[12]*lam.var.sv[i,y,s,12] + sigma.lam[13]*beta.lam[13]*lam.var.sv[i,y,s,13] #+ sigma.eps.lam[i]*eps.lam[i]

  

  logit(p.sv[i,y,s])  ~ dnorm(mu.p.sv[i,y,s], tau.p)

  

  mu.p.sv[i,y,s] <- beta0.p + beta.p[1]*p.var.sv[i,y,s,1] + beta.p[2]*p.var.sv[i,y,s,2] + beta.p[3]*p.var.sv[i,y,s,3] + beta.p[4]*p.var.sv[i,y,s,4] + beta.p[5]*p.var.sv[i,y,s,5]

  

  Count.sv[i,y,s] ~ dbin(p.sv[i,y,s], N.sv[i,y,s])

  

  #expected count

ex.sv[i,y,s] <- N.sv[i,y,s] * p.sv[i,y,s]

        

#discrepancy

E.sv[i,y,s] <- pow((Count.sv[i,y,s] - ex.sv[i,y,s]), 2) / (ex.sv[i,y,s] + 0.1)

        

#simulated new count

rep.sv[i,y,s] ~ dbin(p.sv[i,y,s], N.sv[i,y,s])

        

#discrepancy

E.rep.sv[i,y,s] <- pow((rep.sv[i,y,s] - ex.sv[i,y,s]), 2) / (ex.sv[i,y,s] + 0.1)

  

}}}

  

 for(l in 1:Id.rv){

#loop through sites

for(t in 1:n_Year.rv[l]){

#loop through years

for(m in 1:n_Season.rv[l,t]){

#loop through seasons

  

  omega.rv[l,t,m] ~ dbeta(1,1)

#state model

N.rv[l,t,m] ~ dZIP(lambda = lambda.rv[l,t,m], zeroProb = omega.rv[l,t,m])

log(lambda.rv[l,t,m]) <- beta0.lam + sigma.lam[1]*beta.lam[1]*lam.var.rv[l,t,m,1] + sigma.lam[2]*beta.lam[2]*lam.var.rv[l,t,m,2] + sigma.lam[3]*beta.lam[3]*lam.var.rv[l,t,m,3] + sigma.lam[4]*beta.lam[4]*lam.var.rv[l,t,m,4] + sigma.lam[5]*beta.lam[5]*lam.var.rv[l,t,m,5] + sigma.lam[6]*beta.lam[6]*lam.var.rv[l,t,m,6] + sigma.lam[7]*beta.lam[7]*lam.var.rv[l,t,m,7] + sigma.lam[8]*beta.lam[8]*lam.var.rv[l,t,m,8] + sigma.lam[9]*beta.lam[9]*lam.var.rv[l,t,m,9] + sigma.lam[10]*beta.lam[10]*lam.var.rv[l,t,m,10] + sigma.lam[11]*beta.lam[11]*lam.var.rv[l,t,m,11] + sigma.lam[12]*beta.lam[12]*lam.var.rv[l,t,m,12] + sigma.lam[13]*beta.lam[13]*lam.var.rv[l,t,m,13] #+ sigma.eps.lam[i]*eps.lam[i]

for(v in 1:n_Survey.rv[l,t,m]){

  

  Count.rv[l,t,m,v] ~ dbin(p.rv[l,t,m,v], N.rv[l,t,m])

  

  logit(p.rv[l,t,m,v])  ~ dnorm(mu.p.rv[l,t,m,v], tau.p)

  

  mu.p.rv[l,t,m,v] <- beta0.p + beta.p[1]*p.var.rv[l,t,m,v,1] + beta.p[2]*p.var.rv[l,t,m,v,2] + beta.p[3]*p.var.rv[l,t,m,v,3] + beta.p[4]*p.var.rv[l,t,m,v,4] + beta.p[5]*p.var.rv[l,t,m,v,5]

  

  

#expected count

ex.rv[l,t,m,v] <- N.rv[l,t,m] * p.rv[l,t,m,v]

        

#discrepancy

E.rv[l,t,m,v] <- pow((Count.rv[l,t,m,v] - ex.rv[l,t,m,v]), 2) / (ex.rv[l,t,m,v] + 0.1)

        

#simulated new count

rep.rv[l,t,m,v] ~ dbin(p.rv[l,t,m,v], N.rv[l,t,m])

        

#discrepancy

E.rep.rv[l,t,m,v] <- pow((rep.rv[l,t,m,v] - ex.rv[l,t,m,v]), 2) / (ex.rv[l,t,m,v] + 0.1)

}}}}

     

     fit <- fit.rv + fit.sv

     fit.rep <- fit.rep.rv + fit.rep.sv

 fit.rv <- sum(fit.i.rv[1:Id.rv])

 fit.rep.rv <- sum(fit.rep.i.rv[1:Id.rv])

 

for(l in 1:Id.rv){

#loop through sites

fit.i.rv[l] <- sum(fit.y.rv[l,1:n_Year.rv[l]]) 

fit.rep.i.rv[l] <- sum(fit.rep.y.rv[l,1:n_Year.rv[l]])

for(t in 1:n_Year.rv[l]){

#loop through years

fit.y.rv[l,t] <- sum(fit.s.rv[l,t,1:n_Season.rv[l,t]]) 

fit.rep.y.rv[l,t] <- sum(fit.rep.s.rv[l,t,1:n_Season.rv[l,t]])

for(m in 1:n_Season.rv[l,t]){

#loop through seasons

fit.s.rv[l,t,m] <- sum(E.rv[l,t,m,1:n_Survey.rv[l,t,m]]) 

 fit.rep.s.rv[l,t,m] <- sum(E.rep.rv[l,t,m,1:n_Survey.rv[l,t,m]]) 

  }

 }

}

  

   fit.sv <- sum(fit.i.sv[1:Id.sv])

 fit.rep.sv  <- sum(fit.rep.i.sv[1:Id.sv])

 

for(i in 1:Id.sv){

#loop through sites

fit.i.sv[i] <- sum(fit.y.sv[i,1:n_Year.sv[i]]) 

fit.rep.i.sv[i] <- sum(fit.rep.y.sv[i,1:n_Year.sv[i]])

for(y in 1:n_Year.sv[i]){

#loop through years

fit.y.sv[i,y] <- sum(fit.s.sv[i,y,1:n_Season.sv[i,y]]) 

fit.rep.y.sv[i,y] <- sum(fit.rep.s.sv[i,y,1:n_Season.sv[i,y]])

 }

}

})

I keep getting : Error in checkForRemoteErrors(val) : 

  3 nodes produced errors; first error: Could not evaluate loop syntax: is indexing information provided via 'constants'?

However, I have all the indexing constants in the constants already. 

Constants <- list(Id.rv = n_Id.rv,

                      n_Year.rv = Year_Iterator.rv,

                      n_Season.rv = Season_Iterator.rv,

                      n_Survey.rv = Survey_Iterator.rv,

                      Id.sv = n_Id.sv,

                      n_Year.sv = Year_Iterator.sv,

                      n_Season.sv = Season_Iterator.sv)

I really am at a lose at this point. Any help would be greatly appreciated. 

[Chris Paciorek]

Hi Anthony,

If you're running this in parallel (which it looks like given the reference to checkForRemoteErrors and to 'nodes'), have you seen our guidance on using nimble in parallel?

If not, please see here:

https://r-nimble.org/nimbleExamples/parallelizing_NIMBLE.html

-chris

--
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/2feb496a-106a-4947-96e7-865b739f74bfn%40googlegroups.com.

[Anthony S]

Yes. I had a similar model already running in parallel. Here is the function:

run.nimble <- function(seed, Data, inits, parameters, Constants, model.code, iter, burn, thn, chains){

  library(Rlab)

  library(nimble)

  

dZIP <- nimbleFunction(

 run = function(x = integer(), lambda = double(), 

                zeroProb = double(), log = logical(0, default = 0)) {

   returnType(double())

   ## First handle non-zero data

   if (x != 0) {

       ## return the log probability if log = TRUE

       if (log) return(dpois(x, lambda, log = TRUE) + log(1 - zeroProb))

       ## or the probability if log = FALSE

      else return((1 - zeroProb) * dpois(x, lambda, log = FALSE))

   }

   ## From here down we know x is 0

   totalProbZero <- zeroProb + (1 - zeroProb) * dpois(0, lambda, log = FALSE)

   if (log) return(log(totalProbZero))

   return(totalProbZero)

 })

#here we create the function with the zero inflated poison distribution

rZIP <- nimbleFunction(

 run = function(n = integer(), lambda = double(), zeroProb = double()) {

   returnType(integer())

   imtructuralZero <- rbinom(1, prob = zeroProb, size = 1)

   if (imtructuralZero) return(0)

   return(rpois(1, lambda))

})

#random nummber genetor required by nimble

registerDistributions(list(

    dZIP = list(

        BUGSdist = "dZIP(lambda, zeroProb)",

        discrete = TRUE,

        range = c(0, Inf),

        types = c('value = integer()', 'lambda = double()', 'zeroProb = double()')

     )))

#registering the new distribution

  

assign('dZIP', dZIP, envir = .GlobalEnv)

assign('rZIP', rZIP, envir = .GlobalEnv)

#amign both function so they can be used. Must do because we are working in parallel 

  

  Model <- nimble::nimbleModel(code = model.code, constants = Constants, data = Data, inits = inits)

  #make nimble model

  

  modelconf <- configureMCMC(Model, monitors = parameters)

  #configure model

  

configureRJ(conf = modelconf,

            targetNodes = c("beta.lam[1]", "beta.lam[2]", "beta.lam[3]", "beta.lam[4]", "beta.lam[5]", "beta.lam[6]", "beta.lam[7]", "beta.lam[8]","beta.lam[9]", "beta.lam[10]", "beta.lam[11]", "beta.lam[12]", "beta.lam[13]"#, "eps.lam"

                            ),

            indicatorNodes = c("sigma.lam[1]", "sigma.lam[2]", "sigma.lam[3]", "sigma.lam[4]", "sigma.lam[5]", "sigma.lam[6]", "sigma.lam[7]", "sigma.lam[8]","sigma.lam[9]", "sigma.lam[10]", "sigma.lam[11]", "sigma.lam[12]", "sigma.lam[13]"#, "sigma.eps.lam"

                               ),

            control = list(mean = 0, scale = 2))

#configure reversible jump

MCMC <- nimble::buildMCMC(modelconf)

#build the MCMC

C.Model <- nimble::compileNimble(Model)

#compile the model

C.MCMC <- nimble::compileNimble(MCMC)

#compile the MCMC

results <- nimble::runMCMC(C.MCMC, niter = iter, nburnin = burn, thin = thn, nchains = chains, setSeed = seed, samplesAsCodaMCMC = FALSE, samples = TRUE, summary = FALSE)

#run the model

return(results)

#return the samples

  

} 

I have some sections commented out just for development. This code works for the simpler models without the two different nests.

[Chris Paciorek]

hi Anthony, could you please run this model not in parallel and see if you can see more of the error message? It should indicate which variable is not being found, as illustrated in this simple example, where it says 'n' cannot be found:

> code=nimbleCode({
+ for(i in 1:n)
+ for(j in 1:m)
+ y[i,j]~dnorm(0,1)
+ })
> m=nimbleModel(code)
defining model...
Error in nm_seq_noDecrease((1), (n)) : object 'n' not found
Error in determineContextSize(context, useContext, constantsEnvCopy) :

  Could not evaluate loop syntax: is indexing information provided via 'constants'?

If you're still stuck, please provide us with the actual values of the constants you are using so we can reproduce the error with the code you have already provided.

-chris

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/2cb323c3-99cf-4b91-97f0-26dd1722eadbn%40googlegroups.com.

[Anthony S]

I did it outside the function and I get this error:

defining model...

Error in if (a>b) { : missing value where True/False needed

Error in determineContextSize(context, useContext, constantsEnvCopy) : Could not evaluate loop synthax: is indexing information provided via 'constants'?

I do not know where a>b would be.

[Perry de Valpine]

Hi Anthony,

We might need to ask you to send your inputs (off-list if you prefer) so one of us can run the code and look.

But first let me suggest a strategy to debug this.  I think it is not happy with some for-loop extent, i.e whatever is on the right side of the colon.  A first step would be to double-check that the constants you are providing are correctly formed, e.g. have the correct number of indices.  A second step would be to replace some of those by a placeholder to rule out each one being the problem.  E.g. if you replace

for (m in 1:n_Season.rv[l,t])

with

for (m in 1:3) # I chose 3 arbitrarily

and the problem goes away, then you will know n_Season.rv[l, t] is part of the problem.  You can probe each for-loop in this way.  It's just for 

By the way, a possible problem would be if you have any 0s in n_Season.rv[l,t] and the like, I think.

A more adventurous approach would be to set options(error = recover).  That is an R flag that will through you into a browser (aka debugger) when the error occurs.  You can choose the layer of the call stack to enter and try to see the contents of context, useContext, and constantsEnvCopy as they are being used.  They might point you to the for-loop that is the problem.

I also want to point out you could simplify your model code somewhat by vectorizing and splitting.  E.g. you could replace

log(lambda.sv[i,y,s]) <- beta0.lam + sigma.lam[1]*beta.lam[1]*lam.var.sv[i,y,s,1] + sigma.lam[2]*beta.lam[2]*lam.var.sv[i,y,s,2] + sigma.lam[3]*beta.lam[3]*lam.var.sv[i,y,s,3] + sigma.lam[4]*beta.lam[4]*lam.var.sv[i,y,s,4] + sigma.lam[5]*beta.lam[5]*lam.var.sv[i,y,s,5] + sigma.lam[6]*beta.lam[6]*lam.var.sv[i,y,s,6] + sigma.lam[7]*beta.lam[7]*lam.var.sv[i,y,s,7] + sigma.lam[8]*beta.lam[8]*lam.var.sv[i,y,s,8] + sigma.lam[9]*beta.lam[9]*lam.var.sv[i,y,s,9] + sigma.lam[10]*beta.lam[10]*lam.var.sv[i,y,s,10] + sigma.lam[11]*beta.lam[11]*lam.var.sv[i,y,s,11] + sigma.lam[12]*beta.lam[12]*lam.var.sv[i,y,s,12] + sigma.lam[13]*beta.lam[13]*lam.var.sv[i,y,s,13] #+ sigma.eps.lam[i]*eps.lam[i]

with

log(lambda.sv[i,y,1:n_Season.sv[i,y] ]) <- beta0.lam + sigma.lam[1]*beta.lam[1]*lam.var.sv[i,y, 1:n_Season.sv[i,y] ,1] + sigma.lam[2]*beta.lam[2]*lam.var.sv[i,y, 1:n_Season.sv[i,y] ,2] + ... etc. 

and then with

part_1[i, y,  1:n_Season.sv[i,y] ] <- beta0.lam + sigma.lam[1]*beta.lam[1]*lam.var.sv[i,y, 1:n_Season.sv[i,y] ,1] + sigma.lam[2]*beta.lam[2]*lam.var.sv[i,y, 1:n_Season.sv[i,y] ,2]

part_2[i, y,  1:n_Season.sv[i,y] ] <- sigma.lam[3]*beta.lam[3]*lam.var.sv[i,y, 1:n_Season.sv[i,y],3] + ...

log(lambda.sv[i,y,1:n_Season.sv[i,y] ]) <-  part_1[i, y,  1:n_Season.sv[i,y] ] + part_2[i, y,  1:n_Season.sv[i,y] ]

You can vectorize in up to two dimensions, but in your problem the array looks doubly ragged so that won't easily work.  You could switch to a long format with nested indexing and possibly vectorize over more of it that way.

The reason for vectorizing is some speed gain but also a reduction in the number of nodes in the graphical model, which simplifies a lot of internal steps.  It only makes sense when all elements will be recalculated together anyway, which looks like the case.  E.g. if beta.lam[1] or one of the other coefficients is being updated by MCMC, all the lambda.sv[i, y, s] elements will be recalculated anyway, so might as well vectorize them.  

The reason for splitting is that you have a large number of terms and by splitting you will effectively cache half of the calculations at a time.  E.g. when beta.lam[1] is being updated, all of the part_1 calculations will be redone but the part_2 terms will not be redone, and their stored values will be used.  And vice-versa for some coefficient that enters part_2 but not part_1.

I hope it helps. 

-Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/5c4ddb41-9131-4948-91a0-fadf32d28f64n%40googlegroups.com.

[Anthony S]

Hey Perry,

Great tip on the loop troubleshooting. I was going through and actually every iterator was messed up. I was using reshape2 to create these arrays, but unfortunately they were not arranged in descending order. Therefore, the indexes were not aligning correctly. I had never thought to replace each index to diagnose the problem. I really appreciate the help. 

I also appreciate the tips to get rid of the lowest loop. With a model this size, I am all about efficiency. I do have a question regarding the advice to split the regression. If I am using reversible jump on those variables, would spitting the regression impact the indicator coefficients? Would splitting the regression have any impact other than speed?

Thanks!

Anthony

[Perry de Valpine]

Great to hear you got to the bottom of it!  Congrats on successful debugging.

In answer to your questions, I think the splitting approach should work with RJ and would not change the model's mathematical definition in any way, so it shouldn't change the posteriors or anything other than speed.  It's always hard to predict what efficiency ideas will make small vs. large impacts.

-Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/1e9ccaec-5b23-41af-b05f-640fe1ce57dbn%40googlegroups.com.

[Matthijs Hollanders]

Hi Perry,

Now that Anthony successfully troubleshot his model, perhaps I can swoop in and hijack the thread to ask you a question about vectorization. I'm still trying to wrap my head when it's a good idea to vectorize. Is vectorizing the matrix elements of a state transition matrix a good idea? It sounds like it could be because those parameters would get updated at the same time, presumably. And if vectorizing only goes in two dimensions, what would be the best course of action when there's more than two dimensions (e.g. which dimensions would you choose to vectorize over and which would you use a for loop for)?

Kind regards,

Matt

[Perry de Valpine]

Thanks for the follow-up question.

Here is an example from the workshop we started today ( https://github.com/nimble-training/nimble-virtual-2021, module 4)

A toy linear model with vectorization that makes sense:

code <- nimbleCode({

  intercept ~ dnorm(0, sd = 1000)
  slope ~ dnorm(0, sd = 1000)
  sigma ~ dunif(0, 100)
  predicted.y[1:4] <- intercept + slope * x[1:4] # vectorized node
  for(i in 1:4) {
    y[i] ~ dnorm(predicted.y[i], sd = sigma)
  }
})

predicted.y is vectorized.  This makes sense because if intercept or slope is sampled by MCMC, all of predicted.y will be calculated anyway.  By comparison, here is the unvectorized version:

code <- nimbleCode({
  intercept ~ dnorm(0, sd = 1000)
  slope ~ dnorm(0, sd = 1000)
  sigma ~ dunif(0, 100)
  for(i in 1:4) {
     predicted.y[i] <- intercept + slope * x[i]

    y[i] ~ dnorm(predicted.y[i], sd = sigma)
  }
})

But if the x[i] values are random effects, vectorizing would be inefficient.  This would look like:

code <- nimbleCode({
  intercept ~ dnorm(0, sd = 1000)
  slope ~ dnorm(0, sd = 1000)
  sigma ~ dunif(0, 100)
  predicted.y[1:4] <- intercept + slope * x[1:4] # vectorized node
  for(i in 1:4) {
    x[i] ~ dnorm(0, 1)   # scalar random effects
    y[i] ~ dnorm(predicted.y[i], sd = sigma)
  }
})

In this case, x[1], x[2], ... x[4] are subject to MCMC sampling.  If x[1] is being updated in MCMC, only predicted.y[i] and the log probability of y[i] should need to be calculated.  The non-vectorized version will allow NIMBLE to see that in the graphical model.  The vectorized version will tell nimble that predicted.y[1:4] always come together, which means that if x[1] is being updated, all of predicted.y[1:4] and y[1:4] will be calculated.  The MCMC won't be incorrect, but there will be wasteful computation.

-Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/5c5133ea-f2c3-4228-9773-9ce543374a55n%40googlegroups.com.