New error: In sizeAssignAfterRecursing: 'rpen' is not available or its output type is unknown.

[Annie Shapiro]

Hi again, 

I was previously able to get my code to run, but I'm now getting the following error and I don't believe I've changed anything. I've attached sample data if that's helpful. I have the following code, and get an error that I need to define the random function for my user-defined distribution. After reading the user manual, it doesn't appear that I should need it. Furthermore, this code previously ran without issue. Any thoughts as to what is happening here? 

Thanks!
Annie

Here is my code:

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

f_vec <- nimbleFunction(
  run = function(x = double(0)){
   
    # #Phl:
    k <- 6.06
    r <- 2.94

    tmp_vec <- rep(NA, x)
   
    for (i in 1:x) {
      tmp_vec[i] <- 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)*lambda[1:i])
   
    #likelihood
    x[i] ~ dpois(mu[i])}
 
  #priors
  lambda[1:L] ~ dpen(theta = theta, L = L)
  # k ~ T(dnorm(6.06, .1), 0, )
  # r ~ T(dnorm(2.94, .05), 0, )
 
  #hyperprior
  theta ~ dunif(1, 100)
})



modelConsts <- list(L = nrow(dat))

modelData <- list(x = dat$Notifications)

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


model <- nimbleModel(modelCode, data = modelData,  constants = modelConsts) #inits = initsList,
#====================================================================
# Compile and run model
#====================================================================

# compile into nimble
Cmodel <- compileNimble(model)

And here is the error I get when compiling my model: 

Error: In sizeAssignAfterRecursing: 'rpen' is not available or its output type is unknown. This occurred for: eigenBlock(model_lambda,1:29) <<- rpen(=1,theta=model_theta[1],L=29) This was part of the call: { eigenBlock(model_lambda,1:29) <<- rpen(=1,theta=model_theta[1],L=29) }

[Annie Shapiro]

An update from my attempts to troubleshoot that. I only get this error if I have already run the compileNimble previously in the R session. If it is a new R session, my code runs just fine. Not sure if this helps explain anything, but I'd love to understand still what is going on here!

Thanks,

Annie

[Perry de Valpine]

Dear Annie,

Thanks for the question and follow up. Good timing, as I was literally sitting here looking at your example when your update email came in. I'm not able to reproduce the problem. I tried running compileNimble before your example code, and it still worked fine on my system. I am wondering if your previous call to compileNimble is creating some kind of specific clash later. Would it be easy to provide a fully reproducible example of the workflow that causes the error? Also can you recall more specifically any version information (nimble and/or R) when it worked previously, and if you are confident your workflow then was identical to what doesn't work now?

Perry

--
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/6b8c1b3f-36b6-4956-9e1f-8d6c865fe2c7n%40googlegroups.com.

[Shapiro, Annie]

Hi Perry, 

Thanks for the follow up. It occurs in the code below if I run it, clear the environment, and then run it again. If I don't clear the environment, I don't get the error. Is this helpful? I've also included my R session info. I haven't updated R or numble since it worked previously, and I'm confident that my workflow hasn't changed. 

R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 22000)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] basicMCMCplots_0.2.7 ggplot2_3.4.4        readxl_1.4.3         dplyr_1.1.3          nimble_1.0.1        

loaded via a namespace (and not attached):
 [1] vctrs_0.6.4         cli_3.6.1           rlang_1.1.1         generics_0.1.3      glue_1.6.2         
 [6] colorspace_2.1-0    pracma_2.4.2        scales_1.2.1        fansi_1.0.5         grid_4.3.1         
[11] cellranger_1.1.0    munsell_0.5.0       tibble_3.2.1        lifecycle_1.0.3     numDeriv_2016.8-1.1
[16] compiler_4.3.1      codetools_0.2-19    igraph_1.5.1        coda_0.19-4         pkgconfig_2.0.3    
[21] lattice_0.21-8      R6_2.5.1            tidyselect_1.2.0    utf8_1.2.3          pillar_1.9.0       
[26] parallel_4.3.1      magrittr_2.0.3      withr_2.5.1         tools_4.3.1         gtable_0.3.4       

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

## 10/5: implement penalized likelihood from Bacchetti 1993

dat <- read_excel("C:/Users/ashap/OneDrive/Desktop/Thesis/my countries notif data.xlsx") %>%
  filter(Country == "Philippines") %>%
  dplyr::select(-c(Country, Notes))

#====================================================================
# functions
#====================================================================
#Philippines:
k <- 6.06; r <- 1/.34 # shape; scale

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

     k <- 6.06 # shape
     r <- 1/.34 # scale = 1/rate

   
     tmp_vec <- rep(NA, x)
     
     for (i in 1:x) {
       tmp_vec[i] <- pgamma(x-(i-1), k, rate = r) - pgamma(x-i, k, rate = r)
     }
     
     return(tmp_vec)
     returnType(double(1))
})

# Cf_vec <- compileNimble(f_vec)

# 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))
  })

# registerDistributions(list(
#   dpen = list(
#     BUGSdist = "dpen(theta, L)",
#     Rdist = "dpen(theta, L)",
#     pqAvail = F
#   )))
 
#====================================================================
# set up MCMC1
#====================================================================
#
# # get parameters for lnorm
# m <- mean(dat$Notifications)
# s <- sd(dat$Notifications)
# location <- log(m^2 / sqrt(s^2 + m^2))
# shape <- sqrt(log(1 + (s^2 / m^2)))

modelCode <- nimbleCode({
 
  for (i in 1:L){
    #get mu
    mu[i] <- sum(f_vec(i)*lambda[1:i])
   
    #likelihood
    x[i] ~ dpois(mu[i])}
 

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

  theta ~ dunif(1, 100)
  })



modelConsts <- list(L = nrow(dat))

modelData <- list(x = dat$Notifications)

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


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

Cmodel <- compileNimble(model)

[Perry de Valpine]

Hi Annie,

I replaced your read_excel line with reading the file you provided previously (which looked already filtered and selected).

By "clear the environment", I tried "rm(list = ls())".

I was able to reproduce the problem. Thanks for this report. I will file a GitHub issue (with a minimal reproducible version of the problem) and it will take some investigation, especially if it is, as you said, a new problem.

If you want a workaround so you can use your workflow and not encounter this, you can do

deregisterDistributions("dpen")

when you are clearing the environment.

HTH

Perry