Error: "Atomic order not implemented"

[awstringer1]

Description:

I am using autodiff to calculate 1st, 2nd, and higher order derivatives. I am getting the error "Error in getParameterOrder(data, parameters, new.env(), DLL = DLL) : Atomic 'XXX' order not implemented."

where XXX refers to some atomic function. I've seen this happen for more than one function; I include a specific example below.

I think the problem is related to the manner by which CppAD chooses to execute forward or reverse mode AD. Lines 455 – 458 here describe this error as possibly occurring when CppAD chooses forward mode AD automatically, but the TMB only has reverse mode implemented, which I gather is the case for atomic functions (I have replicated the error with a user-defined function which I know only has reverse mode implemented, since I implemented it). Lines 199 – 237 here show the logic that CppAD uses to choose, and I can see the part that was (I assume) added specifically for TMB.

I am wondering if (a) this is the problem, and if so then (b) if there is a way to force CppAD to always use reverse mode, OR (c) if there is a way to implement the relevant forward mode calculations for TMB's atomic functions.

Or perhaps the problem is something else entirely! Either way, thank you for your time.

Reproducible Steps:

The following is the shortest reproducible example I could think up (my real example is very long). I call the atomic function dbinom_robust, and then pass this into autodiff::jacobian. On the R side, I first create a template using data and parameters that cause CppAD to use reverse mode, and the error does not occur. I then create a template using data and parameters that causes CppAD to choose forward mode, and the error is reproduced.

C++:

#include <TMB.hpp>

// Create a function for passing to autodiff::jacobian()
template<class Type>
struct do_dbinom {
  vector<Type> y;
  do_dbinom(vector<Type> y_) : y(y_) {}
  template<typename T>
  vector<T> operator()(vector<T> x) {
    vector<T> out(y.size());
    vector<T> yt = y.template cast<T>();
    for (int i=0;i<yt.size();i++)
      out(i) = -dbinom_robust(yt(i),T(1),x(0),1); // This is an example of an atomic function that yields the error; it happens for others too
    return out; // Return vector dimension equals that of y
  }
};

template<class Type>
Type objective_function<Type>::operator() () {
  DATA_VECTOR(y);
  PARAMETER_VECTOR(x);
  do_dbinom<Type> lsa(y);
  matrix<Type> g = autodiff::jacobian(lsa,x); // This calls CppAD's Jacobian, causing the error depending on whether it chooses forward or reverse mode.
  return g(0);
}

In R:

library(TMB)
tmbpath <- ".../00_atomic_error_example_grad"
compile(paste0(tmbpath,".cpp"))
dyn.load(dynlib(tmbpath))

x <- 0
# Call it with dim(y) = dim(x), so that CppAD chooses reverse mode
template <- MakeADFun(
  data = list(y=1),
  parameters = list(x=x),
  DLL = '00_atomic_error_example_grad'
)
# No error

# Now call it with dim(y) = 10 > dim(x), so that CppAD chooses forward mode
yy <- rbinom(10,1,.5)
template2 <- MakeADFun(
  data = list(y=yy),
  parameters = list(x=x),
  DLL = '00_atomic_error_example_grad'
)
# This throws an error:
# Error in getParameterOrder(data, parameters, new.env(), DLL = DLL) : 
#   Atomic 'log_dbinom_robust' order not implemented.

TMB Version:

> packageVersion("TMB")
[1] ‘1.9.0’

R Version:

> R.version.string
[1] "R version 4.1.2 (2021-11-01)"

Operating System:

macOS Monterey Version 12.5 (M1 processor)

Thank you!

[kaskr]

From version 1.9.2 there's an argument 'max.order' to the compile() function. Please see the help page.

[awstringer1]

Thank you! I updated my R and TMB, and re-compiled trying different max.order up to max.order = 20. The error persists. From ?compile:

The argument max.order controls the maximum derivative order of special functions (e.g. pbeta) generated by the compiler. By default the value is set to 3

My understanding is that my example above should require at most 3 derivatives of log_dbinom_robust; one for the autodiff::jacobian and then 2 more for the Hessian function created by MakeADFun. Therefore, max.order=3 (default) ought to have worked, and if the error had been caused by an insufficiently large max.order, max.order=20 ought to be more than enough.

I created a smaller example in which I don't use autodiff:

#include <TMB.hpp>

template<class Type>
Type objective_function<Type>::operator() () {
  DATA_VECTOR(y);
  PARAMETER_VECTOR(x);
  vector<Type> out(y.size());
  for (int i=0;i<y.size();i++)
    out(i) = dbinom_robust(y(i),Type(1),x(0),1);
  return out.sum();
}

library(TMB)
tmbpath <- ".../00_atomic_error_example"
compile(paste0(tmbpath,".cpp"),max.order = 1)
dyn.load(dynlib(tmbpath))

y <- rbinom(10,1,.5)
template <- MakeADFun(
  data = list(y = y[1]),
  parameters = list(x=0),
  DLL = '00_atomic_error_example'
) # Reverse mode

template2 <- MakeADFun(
  data = list(y = y),
  parameters = list(x=0),
  DLL = '00_atomic_error_example'
) # Forward mode

However, neither of these throw an error, even compiled with max.order = 1, and both give correct derivatives up to order 2 using template2$gr and template2$he.

I appreciate your help!

[kaskr]

Forgot to say that framework ="TMBad" is required.

[awstringer1]

Thank you! This does appear to have solved it. I will use TMBad going forward.