Including autodifferentiation c++ code in Rstan

[sirmurphalot]

In this (https://mc-stan.org/rstan/articles/external.html) article, it is suggested that the rstan “includes” option can be used to allow a stan model to calculate a derivative in c++. I’m having trouble finding examples of this, and my recent attempts at implementing this have been fruitless. The code for the derivative is coming from this (https://arxiv.org/pdf/1509.07164.pdf) paper; for this instance I’m trying to use a gamma. The R call is:

library(rstan)
library(rstanarm)
library(ggplot2)
library(StanHeaders)
fitted_model_fid1 = stan_model( "gamma.stan", allow_undefined = TRUE,verbose = TRUE,
           includes = paste0('\n#include "', 
                             file.path(getwd(), 'lgammaJacLinf.hpp'), '"\n')) 

where my stan model is

functions {
        real lgammaJacLinf(vector x, real alpha, real lambda);
}
data {
        int<lower=1> N;
        vector[N] y;
}
parameters {
        real<lower=0>   alpha;
        real<lower=0>   lambda;
}
model {
        y ~ gamma_lpdf(alpha,1/lambda);
        target+=lgammaJacLinf(y,alpha,lambda);
}

and the two header files I want to use are attached. Any advice would be quite appreciated. I asked something similar on this (Trying to run autodifferentiation example) thread, and have since confirmed that I am not running Catalina. I decided to start a new thread to make my issue clear and to lay out completely my progress.

gammaCDFderiv.hpp (226 Bytes) lgammaJacLinf.hpp (737 Bytes)

[bgoodri]

In your lgammaJacLinf.hpp, the first lines is

#includes <vector>

which should be

#include <vector>

After I fixed that, the next error was

/tmp/gammaCDFderiv.hpp:3:1: error: unknown type name ‘real’
real gammaCDFderiv(real alpha, real x){

because you wrote Stan syntax in a C++ header file. Those should be templated types.

After that, I got bored. Just go through the error messages one by one until it is fixed.

[sirmurphalot]

Hey Ben – how are you getting that error output? I’d be happy to follow things error by error, but the output from R has been less informative than what you’ve posted above. I tried adding the verbose = TRUE argument, but this just gives me a deluge of meaningless output.

For example, attached are my updated files (using C++ templates), which give the error :

running command '/Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB fileb73a60e25b0f.cpp 2> fileb73a60e25b0f.cpp.err.txt' had status 1Error in compileCode(f, code, language = language, verbose = verbose) : 
  Compilation ERROR, function(s)/method(s) not created! In file included from fileb73a60e25b0f.cpp:8:
In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/src/stan/model/model_header.hpp:4:
In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math.hpp:4:
In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/mat.hpp:4:
In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core.hpp:5:
In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/rev/core/build_vari_array.hpp:4:
In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:4:
In file included from /Library/Frameworks/R.framework/Versions/3.5/Resources/l
Error in sink(type = "output") : invalid connection

from my R console. Existing documentation on Error in sink(type = "output") : invalid connection suggest a change in my CFLAGS for my R’s Markvars file. I’ve done this, but the error remains.

lgammaJacLinf.hpp (735 Bytes) gamma.stan (253 Bytes) gammaCDFderiv.hpp (249 Bytes)

[bgoodri]

verbose = TRUE and then search for error: (with colon).

[sirmurphalot]

I’m beginning to feel a little sheepish. I’ve done what you’ve suggested, but none of the errors are coming from files that I have written. This makes me think (perhaps erroneously) that I’m doing something worse than a syntax error. Attached is the output.
output.txt (65.7 KB)

[bgoodri]

It looks like your C++ toolchain is also messed up. Does

library(rstan)
example(stan_mode, run.dontrun = TRUE)

work for you?

[sirmurphalot]

Hmmm…

> library(rstan)
> example(stan_mode, run.dontrun = TRUE)
no help found for ‘stan_mode’

that would be a no.

[bgoodri]

I meant to type stan_model.

[sirmurphalot]

It looks like a similar group of 20 errors (see attached).
output.txt (62.6 KB)

[bgoodri]

OK, your C++ toolchain is messed up. Did you use

rmacoslib/r-macos-rtools

[martinmodrak]

Hi,
following the thread, it seems you are a novice in C++ development (this can be hard to distinguish from people not putting much work into their questions, which is probably why @bgoodri 's tone was less than friendly). While learning C++ can be hugely empowering and lets you do cool stuff you can’t do otherwise, it is also very challenging. And Stan C++ code is currently optimized for maintainability and speed, not for friendliness to beginners. If you plan to undertake this journey, keep in mind that we are mostly unable to provide support to C++ learners here (there are other sites better suited for that). As well as basic syntax, you would need at least a basic knowledge of templates and template specializations.

An example where I do some custom gradient computations that might get you started faster can be found at

BesselK with the order (v) as real (parameter)

@lionel68 The PR was closed for now, but there is hope :-) It is currently out of my short-term plan to complete this, but if there is need and will to cooperate on a final push to make it happen, then a lot of the work has been done. If you read the comments on the PR (start at the end) you’ll notice that there is a remaining part of the parameter space where the code has large error. This needs to be resolved if the code is to be included in the main codebase and I ran out of ideas. The comp…

It is way more complex than you need but I don’t have a better example at hand.

But I would suggest you think twice whether you really need to delve into C++. You seem to be trying to invoke the autodiff library. But this library is already used when you write your functions in Stan, so unless you have a better formula for the derivative, no improvement will happen.

Are you sure the derivatives are what is slowing down your model? I also have the experience that in some cases custom derivatives were slower than the autodiff derivatives. Or do you have other concerns than speed that motivate you to write your own derivatives?

[sirmurphalot]

It’s quite alright. My impression of @bgoodri is that they have a sort of gruff caring I would expect from a grandfather, or a very old physical therapist.

My overarching aim is to compute a gradient of an arbitrarily defined function inside the STAN function block – my concerns were not actually runtimes (although they may yet become runtimes). The only examples close to this I’ve found online do it in C++, and I did not find a general gradient function here. The .grad() function isn’t available in STAN, but in C++, I believe. Perhaps allow me to rephrase now that I’ve realized I’m out of my depth: is it possible to compute the gradient, in the STAN function block, of another function that has also been defined in the STAN function block?

Thanks so much for all this help. While I wait on your response I’ll start reading through what each of you have linked.

[martinmodrak]

My overarching aim is to compute a gradient of an arbitrarily defined function inside the STAN function block

Ah! Well that is different. There are two possiblities:

1) The gradient will not be part of computation for the log density of the model.
i.e. you only need the gradient in transformed data or generated quantities block. This is easy and can be done with nested autodiff using something like:

Stan side

functions {
  real my_gradient_da(real a, real b);
}

C++ side:

double
  my_derivative_da(const double& a, const double& b, std::ostream *msgs) {
    start_nested();

    var a_var(a);
    var res = the_function_i_want_to_differentiate(a_var, b);
    res.grad();
    double result = a_var.adj();
    recover_memory_nested();
    return result;
  }

It can also be modified to return all gradients as vector (would you need that?).
(But I’d like to check with @bgoodri whether using nested autodiff is really necessary/the best option here)

The above will result in nasty compilation errors when used in transformed parameters or model blocks but will handle both built-in and user-defined functions well.

2) The computed gradient will eventually influence the target density.
This is more tricky because then you need the gradient of the gradient itself wrt. parameters of your model (i.e. second-order derivatives). While Stan’s autodiff should be capable to do this for you, it is AFAIK a) still an experimental feautre and b) advanced stuff with little documentation. In this case I would suggest you find formula for your derivative directly (e.g. with Wolfram Alpha) and implement it as a Stan function. Everything else is IMHO likely to cause you headaches.

Obviously, first make sure your toolchain is setup and you can compile Stan models without external C++.

Hope that helps.