strange slowdown for compiled model relative to uncompiled

[John Clare]

Hi all,

I've been working on trying to find a tractable way to fit a model, and have been perplexed by some strange slowdowns I'm seeing after the model is compiled. (I'm seeing some understandable warnings, but no errors during the building/compiling phases). For example, summing the gettimes() call for 1, 10, and 25 iterations using the uncompiled model returns ~10s, ~105s, and 253s time sampling; for the compiled version, these values are ~3000s, ~19750s, and ~40000s (so, about 11 hours vs. 4+ minutes to sample 25 iterations when compiled or not compiled). In some cases, the particular samplers seem maybe 1-2 x faster compiled, and in others, they seem to run about 200x more slowly. 

Testing the processing time and values returned by calculate(model) vs. calculate(CMCMC$model) for different nodes, by and large, things are relatively similar (CMCMC$model a little faster for some things). There are two multivariate nodes (one deterministic, one stochastic) where the uncompiled model calculates much, much faster (.1 s vs. 13 seconds, .06 seconds vs 2.5 minutes). For the latter, calculate also returns different logProbs, which was surprising in that they share initial values/data and  model$initializeInfo() does not seem to return any missing inits. 

I'll share the nimbleFunctions that seem to exhibit the discrepancies in returned values and time using calculate at the bottom in the hopes that maybe there's just an issue with how the functions are defined that I'm not seeing. Also happy to share full code, data, and a table where I've been keeping track of some of the getTimes discrepancies if nothing immediately pops up and it would be helpful. Thanks in advance for any resolution here.

John

###Function one

gpp<-nimbleFunction(
  run = function(D_star = double(2),
                 phi=double(0),
                 m=integer(0), n=integer(0),
                 w_z =double(1), D_site_star=double(2), e_z=double(1)){ 
    returnType(double(1))  # return a vector of  effects
    w<-numeric(length = n, init = FALSE)
    sigma_e_tilde<-rep(NA, n)
    Cstar<-exp(-D_star * phi)
    inv_Cstar <- inverse(Cstar)
    w_star <- t(chol(Cstar)) %*% w_z ##w_z~N(0, 1)
    C_site_star <-  exp(-D_site_star * phi)
    w <- c(C_site_star %*% inv_Cstar %*% w_star)
   
    C_ss_inv_Cstar <- C_site_star %*% inv_Cstar
    for (b in 1:n) {
      sigma_e_tilde[b] <- 1 - inprod(C_ss_inv_Cstar[b, ], C_site_star[b, ])
      if (sigma_e_tilde[b]<0) {sigma_e_tilde[b]<-0} ###shouldn't be necessary, but...
      w[b]<- w[b] + e_z[b] * sqrt(sigma_e_tilde[b]) 
    }
    return(w)
  }) 

###Function 2

dberngllvm<-nimbleFunction(
  run=function(x=double(1), alpha=double(0), beta=double(1), bigX=double(2),
               lambda=double(1), eta=double(2), indexing=double(1),
               log = integer(0, default = 0)) {
    eta2<-matrix(NA, length(x), dim(eta)[2])
    for (col in 1:dim(eta)[2]){
      eta2[,col]<-eta[indexing, col]
      ###I wonder if "indexing" is causing the problem in compiled code?
      ###Maybe another loop (i) helps? e.g., eta2[i, col]<-eta[indexing[i], col]
    }
   
    p<-c(iprobit((eta2 %*% lambda)+(bigX %*% beta)+alpha))
    ans<-sum(dbinom(x, 1, p, log=TRUE))
    returnType(double(0))
    if(log) return(ans)
    else return(exp(ans))
  })