Error when extracting matrix as array slice

[Jan-Ole Koslik]

Hello all, 

I recently started using RTMB and encountered an error where I need some help.

My minimal example is not so minimal, but the specific model details do not matter and the bottom line is: I build a 3D array from some parameters, and at the bottom of my likelihood function I extract a slice from it for each observation, which works just fine R. 

But when I want to build the AD function, I get the error:

Error in advector(y) :
  Invalid argument to 'advector' (lost class attribute?)

Minimal example:

library(RTMB)

# generating some data from Poisson HMM with covariate effect on transition probabilities
set.seed(123)
n = 100 # number of data points
lambda = c(1, 5) # state-dependent Poisson means
z = rnorm(n) # white noise covariate
beta = matrix(c(-2, -2, 1, -1), nrow = 2) # parameter matrix for t.p.m.
s = rep(NA, n) # simulating Markov chain of hidden states
s[1] = sample(1:2, 1) # initial state
expEta = exp(cbind(1, z) %*% t(beta))
for(t in 2:n){
  # computing t.p.m. from linear predictors via inverse multinomial logit link
  G = diag(2)
  G[!G] = expEta[t,] # covariate dependent transition probability matrix
  G = G / rowSums(G)
  s[t] = sample(1:2, 1, prob = G[s[t-1],])
}
x = rpois(n, lambda[s]) # simulating observations based on hidden states

# negative log-likelihood function for Poisson HMM with covariate effect on transition probabilities
nllk = function(par) {
  getAll(par, dat)
  lambda = exp(loglambda) # state-dependent Poisson mean
  beta = matrix(betavec, nrow = 2, ncol = 2) # parameter matrix for t.p.m.
  expEta = exp(cbind(1, z) %*% t(beta)) # exp linear predictors for t.p.m.
  Gamma = array(NA, dim = c(2, 2, nrow(expEta)))
  for(t in 1:nrow(expEta)) {
    # each tpm slice is computed via inverse multinomial logit link
    G = diag(2)
    G[!G] = expEta[t,]
    Gamma[,,t] = G / rowSums(G)
  }
  delta = c(0.5, 0.5) # initial state probabilities fixed
  allprobs = cbind(dpois(x, lambda[1]), dpois(x, lambda[2])) # state-dependent probabilities
  # forward algorithm
  foo = t(delta) * allprobs[1,]
  phi = foo / sum(foo)
  l = log(sum(foo))
  for(t in 2:nrow(allprobs)) {
    # I get an AD error here when extracting the matrix from the array
    foo = phi %*% Gamma[,,t] * allprobs[t,]
    phi = foo / sum(foo)
    l = l + log(sum(foo))
  }
  -l
}

par = list(loglambda = log(c(1,5)), betavec = c(-2, -2, 1, -1))
dat = list(x = x, z = z)

nllk(par) # works

obj = MakeADFun(nllk, par) # does not work
# Error in advector(y) :
#  Invalid argument to 'advector' (lost class attribute?)

Best regards and thanks in advance,

Ole

[Kasper Kristensen]

A remaining challenge with RTMB is the inplace sub-assignment operator ("[<-"). R doesn't support the case when assigning to a numeric LHS from AD RHS.
Your example uses the problematic assignment here:

G = diag(2) ## numeric
G[!G] = expEta[t,] ## sub-assign from AD

and here:

Gamma = array(NA, dim = c(2, 2, nrow(expEta))) ## numeric
Gamma[,,t] = G / rowSums(G)  ## sub-assign from AD

You'll currently need the following minimal changes for the code to work:

1. Add an explicit overload at the beginning of your function:

   nllk <- function(par) {
      "[<-" <- ADoverload("[<-")
      ...
   }


2. Replace NA by e.g. NaN in the line:

   Gamma = array(NA, dim = c(2, 2, nrow(expEta)))

   Reason: RTMB uses 'is.numeric' internally to decide if a cast is needed but surprisingly (to me) is.numeric(NA) returns FALSE...

Hopefully, in the future these rough edges can be ironed out in RTMB.

[Ben Bolker]

  For what it's worth I think you should be able to use NA_real_ instead of NaN if you prefer for the second part (class(NA) is 'logical', class(NA_real_) is 'numeric')

--
To post to this group, send email to us...@tmb-project.org. Before posting, please check the wiki and issuetracker at https://github.com/kaskr/adcomp/. Please try to create a simple repeatable example to go with your question (e.g issues 154, 134, 51). Use the issuetracker to report bugs.
---
You received this message because you are subscribed to the Google Groups "TMB Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to tmb-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/tmb-users/5a0720f4-aead-470c-a34f-2d274de9b642n%40googlegroups.com.

[Jan-Ole Koslik]

Thanks Kasper, that fixed it. 

However I'm confused that the same did not happen when I only built one matrix, but in the same way using NA.

But maybe I just don't get it :D

Anyway, thanks again!

[james...@gmail.com]

So as I understand it, there's three primitives that benefit from explicit overloading within functions:

```R

  "c" <- ADoverload("c")

  "[<-" <- ADoverload("[<-")

  "diag<-" <- ADoverload("diag<-")

```

Having tripped over the first two ... is there any reason not to just define these at the beginning of each function that we define (i.e., computational cost)?

[Kasper Kristensen]

> any reason not to just define these at the beginning of each function ?

The main reason for not defining them is that it's ugly. If a better solution comes up in the future these ADoverloads will likely become deprecated.

[Jan-Ole Koslik]

I mean this is probably naive, but can't you just let MakeADFun() create a wrapper first, so if f is the user-defined function do

g = function(par){

environment(f) = environment()

"c" <- ADoverload("c")

"[<-" <- ADoverload("[<-")

 "diag<-" <- ADoverload("diag<-")

f(par)

}

and then actually create the AD function from g?

Still kinda ugly, but not for the user :D

[Kasper Kristensen]

janole,

It's a reasonable suggestion but there are some issues with this construct:

1. While the function 'f' will see the ADoverloads, what about functions called by f?
2. If 'f' is created as part of a package it will no longer see the package namespace if we modify the environment.
3. In my experience, the R byte compiler might optimize away the ADoverloads if not visible from the function body of 'f'.
4. If the debugging flag is set on 'f', changing its environment will unset the debugging flag.

FWIW, RTMB already attempts a similar solution by attaching the methods temporarily at the beginning of the search path (see ?ADoverload). I made that choice mainly to address 1 and 4, but have later been tricked by 2 and 3.