How to get incidence matrix between inputs and outputs?

[Tomas T.]

Hello,

I am writing a package that will use TMB interface. User will supply the RTMB template like for example this one (this is a simplified one just to illustrate the question):

nll <- function(par, data)
{
getAll(data, par, warn=FALSE)
s <- plogis(beta)
p <- numeric() 
p[1] <- s[1]
p[2] <- s[1]*s[2]
p[3] <- s[3]*s[4]*s[5]
p[4] <- s[3]*s[5]
p[5] <- s[4]
p[6] <- s[1]*s[5]
p[7] <- s[3]
y <- OBS(y)
REPORT(y)
REPORT(p)
-sum(dbinom(x = y, size = 1, prob = p, log = TRUE))
}
par <- list(beta = numeric(5))
data <- list(y = c(0, 1, 1, 0, 1, 1, 1))
cmb <- function(f, d) function(p) f(p, d)
obj <- MakeADFun(cmb(nll, data), par)

Now, I need an incidence matrix between the input parameter vector (here, beta) and the probabilities vector p. I have seen the "Details behind the AD tape" vignette, so I know there is some internal graph representation, so it should be possible.

I did two attempts, but none of them are quite satisfactory:

#### Attempt #1 - use tape with nested MakeADFun, so that I can use the user-supplied template:

F <- MakeTape(cmb(function (par, data) {
obj <- MakeADFun(cmb(nll, data), par)
obj$report()$p
}, data), par)

# Fails with error: Problem with these parameter entries: beta: 1
# Error in (function (data, parameters, map = list(), type = c("ADFun",  :
#  Only numeric matrices, vectors and arrays can be interfaced


#### Attempt #2 - user has to supply another template that returns the vector `p` instead of the likelihood. I wanted to avoid that to simplify my package interface

nll <- function(par, data)
{
getAll(data, par, warn=FALSE)
s <- plogis(beta)
p <- numeric()
p[1] <- s[1]
p[2] <- s[1]*s[2]
p[3] <- s[3]*s[4]*s[5]
p[4] <- s[3]*s[5]
p[5] <- s[4]
p[6] <- s[1]*s[5]
p[7] <- s[3]
y <- OBS(y)
#REPORT(y)
#REPORT(p)
#-sum(dbinom(x = y, size = 1, prob = p, log = TRUE))
return (p) # return the probability vector instead of the likelihood
}

F <- MakeTape(cmb(nll, data), par)

F2 <- F$jacfun(sparse = TRUE)

incidence_matrix <- as(F2(par), "nsparseMatrix")

# incidence_matrix
7 x 5 sparse Matrix of class "ngCMatrix"
             
[1,] | . . . .
[2,] | | . . .
[3,] . . | | |
[4,] . . | . |
[5,] . . . | .
[6,] | . . . |
[7,] . . | . .

The second attempt actually seems to work, but feels like a dirty hack for two reasons:

1) Getting the incidence matrix from structural non-zeros of sparse jacobian feels like a dirty hack. I am not quite sure that's matematically equivalent in all situations. Would feel safer to have it directly somehow from the graph.

2) Requires the user to supply also the modified template, which doesn't return likelihood but the vector of interest. This makes my package interface more complicated.

Thank you very much for tips!

Best regards

Tomas