inst/doc/RTMB-tips.R

In RTMB: 'R' Bindings for 'TMB'

## ----setup, include=FALSE, message=FALSE--------------------------------------
##library(knitr)
library(RTMB)
set.seed(1)
formals(MakeADFun)$silent <- TRUE

## ----benchmark, echo=FALSE----------------------------------------------------
library(microbenchmark)
set.seed(101)

## Loops *must* be fairly long for the caller overhead to vanish
x <- rnorm(1e5)
y <- rnorm(1e5)

## Benchmark functions *must* return something, or RTMB will optimize it all away!
f1_vect <- function(p) {
  sum(p*x+y)
}
f1_forloop <- function(p) {
  z <- 0
  for (i in seq_along(x)) {
    z = z + p*x[i] + y[i]
  }
  z
}

f1_vect_RTMB <- MakeADFun(f1_vect, c(p=1))
f1_forloop_RTMB <- MakeADFun(f1_forloop, c(p=1))

old <- TapeConfig()
TapeConfig(vectorize="enable")
f1_vect_RTMB2 <- MakeADFun(f1_vect, c(p=1))
invisible(TapeConfig(old))

## RTMB *must* be given a new argument for each eval, or it will just return the previous value!
if (requireNamespace("microbenchmark")) {
  mm <- microbenchmark(vect = f1_vect(rnorm(1)),
                       vect_RTMB = f1_vect_RTMB$fn(rnorm(1)),
                       vect_RTMB2 = f1_vect_RTMB2$fn(rnorm(1)),
                       forloop = f1_forloop(rnorm(1)),
                       forloop_RTMB = f1_forloop_RTMB$fn(rnorm(1)),
                       times = 100L)
}

## ----benchmark-plot, echo = FALSE, fig.width = 8, fig.height = 6, message=FALSE----
if (exists("mm") && requireNamespace("tinyplot")) {
  par(las=1)
  mm2 <- transform(as.data.frame(mm), ntime = time/1000,
                   expr = factor(expr,
                                 levels = c("vect", "vect_RTMB2", "vect_RTMB",
                                            "forloop_RTMB", "forloop"),
                                 labels = c("vectorized\nnative R",
                                            "TapeConfig\nRTMB",
                                            "vectorized\nRTMB",
                                            "for loop\nRTMB",
                                            "for loop\nnative R"))
                   )
  tinyplot::tinyplot(ntime ~ expr,
                     data = mm2, type = "violin", log = "y",
                     trim = TRUE, joint.bw = FALSE,
                     xlab = "",
                     ylab = "time (microseconds)")

}