Custom update works when run in R, but not after compiled

[baug...@vt.edu]

Hello,

I'm having an issue where nimble behaves correctly when running a sampler with a custom update in R, but not after it is compiled. I do not get any errors from Nimble, so I don't know what I've done in the custom update that does not work after compilation.

The model is a spatial capture-recapture model for individuals identified by their left and right flanks, but in the case where some or all of the flank matches are not known. The custom sampler involves updating the latent left-flank and right-flank capture histories. When I run the sampler in R through nimble (Rmcmc$run()), all the flanks are being updated, but after it compiles (Cmcmc$run()), only some of the flanks are being updated. This is the model code:

NimModel <- nimbleCode({
  #detection function priors
  p0B~dunif(0,1)
  #if you do not want to share p0L and p0R
  # p0L~dunif(0,1)
  # p0R~dunif(0,1)
  # #if you want to share p0L and p0R
  p0S ~ dunif(0,1)
  p0L <- p0S
  p0R <- p0S
  #Make sure upper bound on sigma uniform prior is appropriate
  sigma~dunif(0,20)
  #data augmentation prior
  psi~dunif(0,1)
  #likelihoods (except for s priors)
  for(i in 1:M) {
    z[i] ~ dbern(psi)
    s[i,1] ~ dunif(xlim[1],xlim[2])
    s[i,2] ~ dunif(ylim[1],ylim[2])
    d2[i,1:J] <- GetD2(s = s[i,1:2], X = X[1:J,1:2],z=z[i])
    pd.B[i,1:J] <- GetDetectionProb(p0=p0B,sigma=sigma,d2=d2[i,1:J],z=z[i])
    pd.L[i,1:J] <- GetDetectionProb(p0=p0L,sigma=sigma,d2=d2[i,1:J],z=z[i])
    pd.R[i,1:J] <- GetDetectionProb(p0=p0R,sigma=sigma,d2=d2[i,1:J],z=z[i])
    y.B.true[i,1:J,1:K] ~ dBernoulliVectorBoth(pd.B[i,1:J],K2D=K2D[1:J,1:K],z=z[i])
    y.L.true[i,1:J,1:K] ~ dBernoulliVectorSingle(pd.L[i,1:J],K2D=K2D[1:J,1:K],z=z[i])
    y.R.true[i,1:J,1:K] ~ dBernoulliVectorSingle(pd.R[i,1:J],K2D=K2D[1:J,1:K],z=z[i])
  }
  #have to trick Nimble to know ID is part of model by using it somewhere in model statement
  IDdummy <- IDdummyfun(ID.L=ID.L[1:n.L],ID.R=ID.R[1:n.R])
  N <- sum(z[1:M])
})# end model

and this is the left flank update. The right flank update is exactly the same, except it considers the other flank. I replace the default nimble update for y.L.true with this. With no known flank linkages, "n.fixed" will be 0, so I should be looping over 1:n.L.

#Left flank update
IDLSampler <- nimbleFunction(
  contains = sampler_BASE,
  setup = function(model, mvSaved, target, control) {
    K2D <- control$K2D
    n.L <- control$n.L
    n.fixed <- control$n.fixed
    prop.scale <- control$prop.scale #can scale the distance proposal, but a value of 1 should be a good choice.
    calcNodes <- model$getDependencies(target)
  },
  run = function() {
    if(n.L>0){ #skip if no lefts
      z <- model$z
      s <- model$s
      y.L.true <- model$y.L.true
      ID.L<- model$ID.L
      pd.L <- model$pd.L
      M <- nimDim(y.L.true)[1]
      J <- nimDim(y.L.true)[2]
      K <- nimDim(y.L.true)[3]
      #precalculate log likelihoods. Can pull from nimble, but structure depends on nimble model structure (vectorized or not)
      ll.y.L <- array(0,dim=c(M,J,K))
      for(i in 1:M){
        if(z[i]==1){
          for(j in 1:J){
            for(k in 1:K){
              if(K2D[j,k]==1){
                ll.y.L[i,j,k]=dbinom(y.L.true[i,j,k],prob=pd.L[i,j],size=1,log=TRUE)
              }else if(K2D[j,k]==2){
                ll.y.L[i,j,k]=dbinom(y.L.true[i,j,k],prob=2*pd.L[i,j]-pd.L[i,j]^2,size=1,log=TRUE)
              }
            }
          }
        }
      }
      ll.y.L.cand <- ll.y.L
      for(l in (n.fixed+1):n.L){
        y.L.true.cand <- y.L.true
        this.i=ID.L[l]
        #select an individual to swap it to
        d2=(s[this.i,1]-s[,1])^2+(s[this.i,2]-s[,2])^2
        propprobs=exp(-d2/(2*(prop.scale*model$sigma)^2)) #distance-based proposal
        propprobs[this.i]=0 #don't choose focal
        propprobs[z==0]=0 #must select an individual with z==1
        if(n.fixed>0){
          propprobs[1:n.fixed]=0 #cannot select an individual with fixed flanks
        }
        propprobs=propprobs/sum(propprobs)
        cand.i=rcat(1,propprobs)
        cand.ID.idx=which(ID.L==cand.i)#which ID.L index is this individual. Will be empty if not in ID.L
        y.L.true.cand[this.i,,] <- y.L.true[cand.i,,]
        y.L.true.cand[cand.i,,] <- y.L.true[this.i,,]
        #update ll.y.L
        for(j in 1:J){
          for(k in 1:K){
            if(K2D[j,k]==1){
              ll.y.L.cand[this.i,j,k]=dbinom(y.L.true.cand[this.i,j,k],prob=pd.L[this.i,j],size=1,log=TRUE)
              ll.y.L.cand[cand.i,j,k]=dbinom(y.L.true.cand[cand.i,j,k],prob=pd.L[cand.i,j],size=1,log=TRUE)
            }else if(K2D[j,k]==2){
              ll.y.L.cand[this.i,j,k]=dbinom(y.L.true.cand[this.i,j,k],prob=2*pd.L[this.i,j]-pd.L[this.i,j]^2,size=1,log=TRUE)
              ll.y.L.cand[cand.i,j,k]=dbinom(y.L.true.cand[cand.i,j,k],prob=2*pd.L[cand.i,j]-pd.L[cand.i,j]^2,size=1,log=TRUE)
            }
          }
        }
        #calculate backwards proposal probs
        d2=(s[cand.i,1]-s[,1])^2+(s[cand.i,2]-s[,2])^2
        backprobs=exp(-d2/(2*(prop.scale*model$sigma)^2)) #distance-based proposal
        backprobs[cand.i]=0 #don't choose focal
        backprobs[z==0]=0 #must select an individual with z==1
        if(n.fixed>0){
          backprobs[1:n.fixed]=0 #cannot select an individual with fixed flanks
        }
        backprobs=backprobs/sum(backprobs)
        lp_initial <- sum(ll.y.L[this.i,,])+sum(ll.y.L[cand.i,,])
        lp_proposed <- sum(ll.y.L.cand[this.i,,])+sum(ll.y.L.cand[cand.i,,])
        log_MH_ratio <- (lp_proposed+log(backprobs[this.i])) - (lp_initial+log(propprobs[cand.i]))
        accept <- decide(log_MH_ratio)
        if(accept){
          y.L.true[this.i,,]=y.L.true.cand[this.i,,]
          y.L.true[cand.i,,]=y.L.true.cand[cand.i,,]
          ll.y.L[this.i,,]=ll.y.L.cand[this.i,,]
          ll.y.L[cand.i,,]=ll.y.L.cand[cand.i,,]
          #swap flank indices
          ID.L[l]=cand.i
          #if cand.i was in ID.L, update this ID index
          tmp=nimDim(cand.ID.idx)[1]
          if(tmp>0){
            ID.L[cand.ID.idx]=this.i
          }
        }
      }
      #put everything back into model$stuff
      model$y.L.true <<- y.L.true
      model$ID.L <<- ID.L
      model$calculate(calcNodes) #update logprob
      copy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE)
    }
  },
  methods = list( reset = function () {} )
)

I can provide a reproducible example, but perhaps someone can easily spot what I've done here that does not work after compiling (and escapes nimble's error checks).

Any help would be much appreciated.

Ben Augustine

[Perry de Valpine]

Hi Ben,

Nothing jumped out at me.  (You'll want to be sure n.fixed + 1 <= n.L to avoid any decreasing sequences, but I bet that is the case.)

I can point you to some deeper debugging tricks.  If you do:

nimbleOptions(buildInterfacesForCompiledNestedNimbleFunctions = TRUE)

then you can more fully access the compiled samplers.  E.g. you can do:

Cmcmc$samplerFunctions[[1]]$run(). ## run one sampler at a time from the list of samplers.

You will need to find the index of the samplerFunctions list that corresponds to your custom sampler.

E.g. here is a sketched debugging workflow to compare uncompiled and compiled behavior for a single sampler.  They should give identical behavior if you carefully set up identical conditions.

R prefixes uncompiled stuff. C prefixes compiled stuff.

for(v in Rmodel$getVarNames()) {

   Cmodel[[v]] <- Rmodel[[v]] # ensure all variables are the same

}

Rmodel$calculate()

Cmodel$calculate() # ensure model log probs all match

Rmcmc$run(niter = 0)

Cmcmc$run(niter = 0) # Do MCMC initialization steps for both uncompiled and compiled, including initialization of the mvSaved object, but no iterations

## Find the sampler of interest

Rmcmc$samplerFunctions # Look at list of samplers

Cmcmc$samplerFunctions

# Dig into these to find yours if it is not obvious.  Using the name argument to nimbleFunction when creating IDLsampler might make it easier to see which is yours.

# Check by target node. Say it is sampler 23

Cmcmc$samplerFunctions[[23]]$target

Rmcmc$samplerFunctions[[23]]$target

# So let's say you have determined sampler 23 is the one you want to dig into.

## Before doing all this, modify IDLsampler to handle as member data any objects you want to inspect later for debugging purposes.  

## E.g. say you want to check whether backprobs is correct.  

## Add this to setup code

backprobs <- numeric(2) ## create it to be recognized as a vector

## And in run code change "backprobs <-" to "backprobs <<-"

# Run the matching samplers uncompiled and compiled with matching RNG seeds.

set.seed(1)

Cmcmc$samplerFunctions[[23]]$run()

Cmcmc$samplerFunctions[[23]]$backprobs # Look at the value of what you want to check, as it was modified by running the sampler.

Cmodel[[target]] # check the compiled model's value of the target node of sampler 23

# Repeat with uncompiled

set.seed(1)

Rmcmc$samplerFunctions[[23]]$run()

Rmcmc$samplerFunctions[[23]]$backprobs # Check if it matches Cmcmc's value

Rmodel[[target]] # check if the uncompiled model's value of the target node was updated identically to that of the compiled model's target node

etc.  If you find a discrepancy between uncompiled and compiled behavior, that is an issue.

I hope this helps!  The nimbleOption setting above is what lets you go more than one level via these R interface objects (R objects that let you access C++ objects).  E.g you can always look into the first level of a compiled nimbleFunction, say Cmcmc$samplerFunctions (samplerFunctions is an object one layer inside of the object that was compiled, Cmcmc).  However to get to the next layer, you need the option set.  It does add overhead in compilation time and memory use, which is why it's off by default.

-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/2b8c900b-6c6b-43df-b1cb-aac9e584540bn%40googlegroups.com.

[baug...@vt.edu]

Much appreciated, Perry. This is exactly what I needed. I feel like I have a new super power! (I deleted my first response a couple hours ago because I failed to follow directions).

The problem was where I was scaling the distance-based proposal by sigma:

propprobs=exp(-d2/(2*(prop.scale*model$sigma)^2))

In R, this is a scalar numeric. When compiled, I'm not sure what it is, but it is not a scalar, despite having only 1 value being updated. I guess it is treated as a vector. Dividing by model$sigma zeroed out every other value of propprobs. I've noticed before in a custom nimble function that I needed to reference the first element of an object that I though should be a scalar. I'm sure there is a good reason for this that I'm not aware of. Anyways, I fixed it by pulling out sigma from the model object and then referencing sigma[1].

propprobs=exp(-d2/(2*(prop.scale*sigma[1])^2))

Thanks for the help!

Ben

[Perry de Valpine]

Hi Ben,

Great to hear that this method worked and got you deeper into the system and you could find your way around.  It has us thinking about documenting these kinds of steps more easily for others in the future.

I'll try to clarify a bit of the run-time logic that can make the types (in this case, vector vs. scalar) confusing, especially for all of us used to the way R works.  On the C++ side, all types are static.  That means it won't decide at run-time that since the length of sigma is 1, it will treat it as vector * scalar multiplication.  The compiler sees component-wise vector*vector and generates C++ code for that. It also does not implement a general recycling rule (although we do do that for distribution arguments).  So it will attempt to do component-wise multiplication, and if the lengths don't match, you won't get what you want.  An added layer that can be confusing is that your sigma comes from a model, and all scalars in models are set up as length-1 vectors, so when you use that model variable directly in nimbleFunction code, you need to explicitly pull out the first element.  We know this can be confusing and hope to make it less so in the future.

-Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/dedb0075-3a78-4678-9e63-3cbeca82233fn%40googlegroups.com.