Incredibly Slow Model Runs

[heather...@gmail.com]

Hi all,

I'm running a complex model that seems to be taking ages to run, even if I were to run it without data to perform a prior predictive check. Right now the model runs about 1 iteration per minute, which is agonizing. 

I suspect the problem is in the way I'm dealing with some latent multinomials in the model, but I can't figure out a way to speed it up so I'm reaching out for help.

Essentially, there are 105 pixels and 12 life stages that an individual can be in. Individuals can move and they can age (change stages). The problem is that none of this is observed, so to get the likelihood, I combined two large multinomials into a single multinomial draw. This turns the problem into a 105*12 vector of probabilities. I think something in this process is making the model go very very slowly. 

Here's what my custom function looks like:

dmmm_psi_move_two <- nimbleFunction(
  run = function(x = double(2),
                 nn = double(2), #n_survivors[1:nState, t, 1:nPix]
                 psi = double(3), #psi[1:nState,1:nState,t,1:nPix]
                 nState = double(0), #12 states
                 nPix = double(0),
                 pixcount = double(1),
                 pixopts = double(2),
                 not0_psi = double(2),#binary indicating impossible transitions
                 probs_norm = double(3),  #movement probabilities
                 log = integer(0, default = 0)) {
    returnType(double(0))
   
    newsize <- sum(nn[1:nState, 1:nPix])
    if(newsize < 1)
      return(0) #unlikely but saves time if the pop dies out!
   
    x_size <- sum(x[1:nState,1:nPix])
    if(x_size != newsize)
      return(-Inf)
   
   
    all_prob = matrix(0, nr = nState*nPix, nc = nState*nPix)
    all_prob_norm = matrix(0, nr = nState*nPix, nc = nState*nPix)
   
    #Combine transitions and movement:
   
    for(q in 1:nPix){
      for(es in 1:nState){
        for(ss in 1:nState){
          if(not0_psi[ss,es] >0){ #skip all the known impossible transitions
            for(pp in 1:pixcount[q]){
              all_prob[ss+nState*(q-1), es+nState*(pixopts[q, pp]-1)] <- probs_norm[q, pixopts[q, pp], es]*psi[ss,es,q]
            } #end pp
          } #end if statement
        } #end ss
      } #end es
    } #end q
   
    for(i in 1:(nState*nPix)){
      if(sum(all_prob[i,]) > 0){
        all_prob_norm[i,] <- all_prob[i,]/sum(all_prob[i,])
      }
    }
    #eww, this will go very slowly
   
    #we want to do one multinomial, based on the outcomes of the size-weighted probabilities in each columns
    #but first we need to count up how many are in each state right now:
    colsizes = c(nn)#just rearranging nn; now goes Stage1pix1, Stage2pix1.... Stage12pix1, stage1pix2..... stage12pixN
    #nn <- n_survivors[1:nState, t, 1:nPix] 12*25
   
   
    newpsis = array(0, dim = nState*nPix)
    for(i in 1:(nState*nPix)){ #weight each psi value by size
      newpsis[i] = sum(all_prob_norm[1:(nState*nPix), i]*colsizes/newsize)
    }
    x2 <- c(x) #rearrange n vector from n[1:nState, t, 1:nPix] to a single vector
   
    #x = n[1:nState, t, 1:nPix], two dimensional argument. Hmm.
    loglike = dmulti(x2, size = newsize, prob = newpsis, log = TRUE)
   
    if(log) return(loglike)
    else return(exp(loglike))
  }
)

I realize it's a long shot, but I was hoping someone had some suggestions for how to increase the speed of this beast. I've attached an example that attempts to simplify most everything else in the model *except* for the slow down part.  

Best,

Heather 

[Perry de Valpine]

Hi Heather,

Thanks for posting this. I was not able to run your code because your zip file didn't include "bigmod.txt". I thought maybe "bigmod_test.txt" would be that so I changed the name but it didn't work and bigmod_test.txt looks like MCMC samples. (And I think files will be smaller saved via saveRDS, but no matter).

Let me make some general suggestions from reading your code.

help(buildMCMC) shows how you can use "time=TRUE" and then "getTimes()" from the compiled MCMC object to see the time spent in each sampler. This will give you an idea of where the computational costs are. Of course you also need to understand what is being computed when sampling a node, and this is usually given by model$getDependencies(node).

If you are confident the problem is your dmmm_psi_move_two, you can profile within that function using the function run.time. I see that help(run.time) doesn't give much, but it is in section 15.8 of the User Manual. You could just print the time of different sections of the code and see whether the huge for loops or the huge dmulti or both or something else is the slow part. There is a saying about pursuing code efficiency that it is better to profile than to guess what is slow, although you might have a good idea already.

Here are some ideas to try in changing how your model or nimbleFunctions are written.

You could remove psi from the model, pass eta_d to dmmm_psi_move_two, and manually use it there for the elements you need. psi looks like a large-sized expansion that duplicates the elements of eta_d only for purposes of using within dmmm_psi_move_two.

You could also remove the deterministic elements of n_post from the model and split the survivors[ss, t, q] into multiple declarations with some nested indexing to achieve what you have written. My guess would be that this is less of a win than removing psi from the model, but I can't be sure.

I'm not following all the details of dmmm_psi_move_two so I might miss the mark but here are some notes on that:

There is a cost to creating large objects such as all_prob and all_prob_norm. One thing you can do is create only one of them but re-use it, e.g.  all_prob[i,] <- all_prob[i,]/sum(all_prob[i,]), because after that you don't need the old all_prob. Another thing you can do is check out User Manual section 12.3 on user-defined distributions that include setup code and thus can retain their own internal variables (like a class). With this feature, you can write dmmm_psi_move_two and/or normalize_move so that they contain a *single* instance (each) of all_prob (created in setup code and assigned in run code using "<<-") that gets re-used each time without allocating+initializing and de-allocating. (The size of all_prob is (105*12)^2 * 8 / 2^20 = 12.1 MB, which is a lot to be doing many thousands or millions of times in MCMC.) Another thing you can do with this approach is hold gridDist in the user-defined function or distribution objects and thus keep it out of the model, where it only adds overhead.

You might be able to usefully combine normalize_move and dmmm_psi_move_two to form the large object(s) even less. You might be able to avoid repeating some of the work across t that does not depend on t. You might be able to compute the dmulti log probability manually as you get your hands on the relevant elements, to avoid forming large objects if their only purpose in existing as single objects is to pass the dmulti. In normalize_move, I think you could vectorize over k. Looking at the big nested loops in dmmm_psi_move_two, since there are (many?) impossible state transitions and so (guessing) many 0s in the large objects, I am wondering if there is a way to create smaller objects with a different indexing scheme.

HTH

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 visit https://groups.google.com/d/msgid/nimble-users/668a53be-dd94-4320-8b60-4e349bf87d3an%40googlegroups.com.

[Dr. Heather Gaya]

Thanks Perry! 

I was unaware of the getTimes() function, though I see it very clearly on the help page. I will try some of those suggestions you mentioned. Perhaps minimizing the re-creation of all_prob will save a lot of space. 

I appreciate the second set of eyes - hopefully some or all of your suggestions can get this working at a more reasonable speed. It doesn't really matter if it takes forever once I'm sure the model is working as I want, but in the meantime... Phew! 

Heather