generic "^" for class "simref"

[James-Thorson-NOAA]

First off, RTMB is awesome! I think it'll be useful for real-world problems involving models that require a user-specified nonlinear function as input (e.g., allowing the user to define a function for nonlinear dynamics, and then using RTMB to estimate a state-space path that results from fitting to data). I'm happy to talk about the user-case I'm exploring if anyone's interested.

Second, I'm trying to use obj$simulate() based on docs at ?Simulation, and I think the docs claim that all standard operations should be supported for Class 'simref'. However, I get an error in my use-case:

Error in Ops.simref(u[1], 1.01) : 
  Class 'simref' does not know generic '^'

I can replicate the problem by modifying the example below. Am I missing something?

s <- simref(4)
s2 <- 2 * s[1:2] + 1
s2[] <- 7
s ## 3 3 NA NA
## Random walk
func <- function(p) {
  u <- p$u
  ans <- -dnorm(u[1]^1.01, log=TRUE) ## u[1] ~ N(0,1)
  ans <- ans - sum(dnorm(diff(u), log=TRUE)) ## u[i]-u[i-1] ~ N(0,1)
}
obj <- MakeADFun(func, list(u=numeric(20)), random="u")
obj$simulate()

[kaskr]

Your example is interpreted as the indirect assignment

u[1]^1.01 ~ N(0,1)

'^' is not generally invertible and as the documentation says:

Indirect assignment works for a limited set of easily invertible functions

'^' can only be used for direct assignment, e.g.

func <- function(p) {
  u <- p$u
  ans <- -dnorm(u[1], log=TRUE)
  ans <- ans - dnorm(u[2] - u[1]^1.01, log=TRUE)
  ans
}
obj <- MakeADFun(func, list(u=numeric(2)), random="u")
obj$simulate()

This works because u[1]^1.01 is used after u[1] has been simulated.

There are some glitches though. For example the following should work but doesn't:

func <- function(p) {
  u <- p$u
  ans <- -dnorm(u[1], log=TRUE)
  ans <- ans - sum(dnorm(tail(u,-1) - head(u,-1)^1.01, log=TRUE))
  ans
}
obj <- MakeADFun(func, list(u=numeric(20)), random="u")
obj$simulate()