Linker error when calling external C++ within Stan function

I have a rather complex model including an ODE component for which I would like to provide external C++ code that implements the analytical solution for the gradient calculation of the expensive (non-linear) bits while allowing automatic differentiation to handle the rest. Because of the complexity of the model I have separated the code into multiple functions to avoid cognitive overload when reading/modifying the code. All works fine when I use pure Stan. However, when I try to compile with external C++ code I keep getting errors like:

ok_wheat_pheno/reprex.o: In function boost::math::tools::promote_args<double, double, double, float, float, float>::type reprex_model_namespace::calling_function<double, double, double>(double const&, double const&, double const&, std::ostream*)': reprex.hpp:(.text._ZN22reprex_model_namespace16calling_functionIdddEEN5boost4math5tools12promote_argsIT_T0_T1_fffE4typeERKS5_RKS6_RKS7_PSo[_ZN22reprex_model_namespace16calling_functionIdddEEN5boost4math5tools12promote_argsIT_T0_T1_fffE4typeERKS5_RKS6_RKS7_PSo]+0x1d): undefined reference to boost::math::tools::promote_args<double, double, double, float, float, float>::type reprex_model_namespace::dof_fun<double, double, double>(double const&, double const&, double const&, std::ostream*)’

I have looked at this thread which seems to demonstrate the same/related problem, but I don’t see a real solution there apart from avoiding calling external C++ functions from within other Stan functions. Is there any other way to make this work? Any guidance here would be much appreciated. I provide a minimal reproducible example below.

I am using cmdstan (Stan code and continuous_functions.hpp code given below) and invoking like so:

STANCFLAGS=--allow_undefined USER_HEADER=/home/palderman/active/ok_wheat_pheno/stan/continuous_functions.hpp make ok_wheat_pheno/reprex

Stan code:

functions{
  
  real dof_fun(real To, real rn, real rd);
  
  real calling_function(real To, real rn, real rd){
    real mu = dof_fun(To, rn, rd);
    return mu;
  }
  
}

data {
  int<lower=0> N;
  vector[N] y;
}

parameters {
  real To;
  real rn;
  real rd;
  real<lower=0> sigma;
}

model {
  y ~ normal(calling_function(To,rn,rd), sigma);
}

External C++ code:

#include <stan/math/rev/core.hpp>

// Functions for dof and analytic partial gradients
namespace stan{
  namespace math{
  
   template<typename T1, typename T2, typename T3>
   typename boost::math::tools::promote_args<T1, T2, T3>::type
     dof_fun(const T1& To, const T2& rn, const T3& rd, std::ostream* pstream__) {
       throw std::logic_error("not implemented");  // this should never be called
     }
   
    double dof_fun(double To, double rn, double rd, std::ostream* pstream__){
      return (To*rd+log(rn)-log(rd-rn))/rd;
    }

    double ddof_drd(double dof, double To, double rn, double rd, std::ostream* pstream__){
      return (To - 1/(rd - rn) )/rd - dof/rd;
    }

    double ddof_drn(double rn, double rd, std::ostream* pstream__){
      return 1/(rn*(rd-rn));
    }

    // All var args
    var dof_fun(const var& To, const var& rn, const var& rd, std::ostream* pstream__){
  
      double dof = dof_fun(To.val(),rn.val(),rd.val(),pstream__);
  
      return var(new precomp_vvv_vari(
        dof,                                                  // value of output
        To.vi_,                                               // input gradient wrt To
        rn.vi_,                                               // input gradient wrt rn
        rd.vi_,                                               // input gradient wrt rd
        1,                                                    // partial wrt To
        ddof_drn(rn.val(), rd.val(),pstream__),               // partial wrt rn
        ddof_drd(dof, To.val(), rn.val(), rd.val(),pstream__) // partial wrt rd
      ));
    }
  }
}

This was an annoying error to work through (edit: I mean that like, it’s opaque what is going on, not saying I am annoyed at you for posting this lol). The model gets compiled to C++ with a fully-templated function definition right above it.

Like this:

template <typename T0__, typename T1__, typename T2__>
stan::promote_args_t<T0__, T1__,
T2__>
dof_fun(const T0__& To, const T1__& rn, const T2__& rd,
        std::ostream* pstream__) ;

So when it gets called in the model file, it has something to point at:

mu = dof_fun(To, rn, rd, pstream__);

So if it’s not finding other functions, the compiler will just point at those and you’ll get link time errors.

Anyhow, you can fix this by putting your code in the wheat_model_namespace namespace (not stan::math). This way your functions will get found and everything compiles as you expect.

Like this: dof.hpp (1.5 KB)

2 Likes

Thanks so much for the insight. I was at a loss for how to proceed. I will try to implement it as you’ve suggested and come back if I run into issues.

P.S Thanks for the clarification that you are not annoyed at me! :)

2 Likes