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
));
}
}
}
```