error trapping unusual indexing in call to user-defined function

[paciorek]

A user (nimble-users list 2023-10-18) tried to use this model code:

for (i in 1:L){
    lambda[i] ~ dgamma(alpha,beta)
    mu[i] <- sum(f(i-1:i)*lambda[1:i])
    x[i] ~ dpois(mu[i]) - exp{-(theta/2)*sum((log(x[1:(L-2)]) - 2*log(x[2:(L-1)]) + log(x[3:L]))^2)}

Untrapped error in building the model was:

Error in if (currentVertexCounts[VID] != numWithinBlock) { : missing value where TRUE/FALSE needed

I think this is caused by f(i-1:i) but didn't investigate in detail. We should trap this and give an informative error message. I'm not remembering offhand what types of sequencing we're expecting to be valid.

[paciorek]

Just reported to the user that I can't reproduce this. the i-1:i indexing seems to work fine.

[paciorek]

Here's a reprex from the user (see 2024-03-22 email for data file):

Sorry for my delay in responding! I didn't see your request, but am now experiencing the error again and went back to this thread to remind myself of the solution. Below are the code and I've also attached data.

library(nimble)
library(dplyr)
library(readxl)
library(ggplot2)
library(basicMCMCplots)
set.seed(1)

dat <- read_excel("updated_country_dat.xlsx") %>%
  filter(Country == "Viet Nam",
         Year < 2020)

k <- 8.07
r <- 4.12
x0 <- 2012

#2007
p0 <- 1-(1-1/7.84)^(k/r)

#2017
p1 <- 1-(1-1/5.06)^(k/r)

#make bottom 0 because that is what the equation likes
l <- p1-p0

#slope
# s <- (p1-p0)/(2016-2007)

dat <- dat %>%
  mutate(p = l/(1+exp(.5*(Year-x0)))+ p0)

# tmp <- dat %>%
#   filter(Year < 2018)
#
# plot(tmp$Year, tmp$p, xlab = "Year", ylab = "p")
# title("(c)", adj = 0)

# remove Notifs_m that were incident before study period
dat <- dat %>%
  mutate(adj = pgamma(row_number(), k, rate = r),
         Notifs_m = round(Notifs_m*adj))

#====================================================================
# functions
#====================================================================

f_vec <- nimbleFunction(
  run = function(x = double(0),
                 p = double(0)){

    k <- 8.07
    r <- 4.12

    tmp_vec <- rep(NA, x)

    for (i in 1:x) {
      tmp_vec[i] <- p*(pgamma(x-(i-1), k, rate = r) - pgamma(x-i, k, rate = r))
    }

    return(tmp_vec)
    returnType(double(1))
  })

# code penalty as prior
dpen <- nimbleFunction(
  run = function(x = double(1),
                 theta = double(0), L = double(0),
                 log = integer(0, default = 0)) {
    returnType(double(0))
    logProb <- -(theta/2)*sum((log(x[1:(L-2)]) - 2*log(x[2:(L-1)]) + log(x[3:L]))^2)
    if(log) return(logProb)
    else return(exp(logProb))
  })

#====================================================================
# set up MCMC
#====================================================================
modelCode <- nimbleCode({

  for (i in 1:L){
    #get mu
    mu[i] <- sum(f_vec(i, p[i])*lambda[1:i])

    x[i] ~ dpois(mu[i])
    }

  #priors
  lambda[1:L] ~ dpen(theta = theta, L = L)
})

modelConsts <- list(L = nrow(dat),
                    p = dat$p,  
                    theta = mean(dat$Notifs_m))

modelData <- list(x = dat$Notifs_m)

initsList <- list(lambda = abs(rnorm(nrow(dat), mean(dat$Notifs_m), sd(dat$Notifs_m))))

model <- nimbleModel(modelCode, data = modelData,  inits = initsList, constants = modelConsts)

[paciorek]

This was fine on nimble 1.1.0 for user and for myself.