# Using external C++ code

The `--allow-undefined`

flag can be passed to the call to stanc, which will allow undefined functions in the Stan language to be parsed without an error. We can then include a definition of the function in a C++ header file.

This requires specifying two makefile variables:

`STANCFLAGS=--allow-undefined`

`USER_HEADER=<header_file.hpp>`

, where`<header_file.hpp>`

is the name of a header file that defines a function with the same name and a compatible signature. This function can appear in the global namespace or in the model namespace, which is defined as the name of the model (either the file name, or the`--name`

argument to stanc) followed by`_namespace`

.

This is an advanced feature which is only recommended to users familiar with the internals of Stan’s Math library. Most existing C++ code will need to be modified to work with Stan, to varying degrees.

As an example, consider the following variant of the Bernoulli example

```
functions {
real make_odds(data real theta);
}data {
int<lower=0> N;
array[N] int<lower=0, upper=1> y;
}parameters {
real<lower=0, upper=1> theta;
}model {
1, 1); // uniform prior on interval 0, 1
theta ~ beta(
y ~ bernoulli(theta);
}generated quantities {
real odds;
odds = make_odds(theta); }
```

Here the `make_odds`

function is declared but not defined, which would ordinarily result in a parser error. However, if you put `STANCFLAGS = --allow-undefined`

into the `make/local`

file or into the `stanc`

call, then the stanc compiler will translate this program to C++, but the generated C++ code will not compile unless you write a file such as `examples/bernoulli/make_odds.hpp`

with the following lines

```
#include <ostream>
double make_odds(const double& theta, std::ostream *pstream__) {
return theta / (1 - theta);
}
```

The signature for this function needs to fulfill all the usages in the C++ class emitted by stanc. The `pstream__`

argument is mandatory in the signature but need not be used if your function does not print any output. Because `make_odds`

was declared with a `data`

argument and only used in `generated quantites`

, a signature which accepts and returns `double`

is acceptable. Functions which will have parameters passed as input in the `transformed parameters`

or `model`

blocks will require the ability to accept Stan’s autodiff types. If you wish to autodiff through this function, the simplest option is to make it a template, like

```
template <typename T>
(const T &theta, std::ostream *pstream__)
T make_odds{
return theta / (1 - theta);
}
```

Given the above, the following `make`

invocation should work

`> make STANCFLAGS=--allow-undefined USER_HEADER=examples/bernoulli/make_odds.hpp examples/bernoulli/bernoulli # on Windows add .exe`

Alternatively, you could put `STANCFLAGS`

and `USER_HEADER`

into the `make/local`

file instead of specifying them on the command-line.

If the function were more complicated and involved functions in the Stan Math Library, then you would need to add `#include <stan/model/model_header.hpp>`

and prefix the function calls with `stan::math::`

.

## Derivative specializations

External C++ functions are currently the only way to encode a function with a known analytic gradient outside the Stan Math Library. This is done very similarly to how a function would be added to the Math library with a reverse-mode specialization. The following code is adapted from the Stan Math documentation.

Suppose you have the following (nonsensical) model which relies on a function called `my_dot_self`

. We will implement this as a copy of the built-in `dot_self`

function.

```
functions {
// both overloads end up using the same C++ template
real my_dot_self(vector theta);
real my_dot_self(row_vector theta);
}data {
int<lower=0> N;
vector[N] input_data;
}transformed data {
// no autodiff for data - will call using doubles
real ds = my_dot_self(input_data);
}parameters {
row_vector[N] thetas;
}model {
0,1);
thetas ~ normal(// autodiff - will call using stan::math::var types
input_data ~ normal(thetas, my_dot_self(thetas)); }
```

If you wanted to autodiff through this function, the following header would suffice^{1}:

```
#include <stan/model/model_header.hpp>
#include <ostream>
template <typename EigVec, stan::require_eigen_vector_t<EigVec> * = nullptr>
inline stan::value_type_t<EigVec> my_dot_self(const EigVec &x, std::ostream *pstream__)
{
const auto &x_ref = stan::math::to_ref(x);
::value_type_t<EigVec> sum_x = 0.0;
stanfor (int i = 0; i < x.size(); ++i)
{
+= x_ref.coeff(i) * x_ref.coeff(i);
sum_x }
return sum_x;
}
```

However, we know the derivative of this function directly. To leverage this, we could use a more complicated form which has two function templates that differentiate themselves based on whether or not derivatives are required:

```
#include <stan/model/model_header.hpp>
#include <ostream>
template <typename EigVec, stan::require_eigen_vector_t<EigVec> * = nullptr,
::require_not_st_var<EigVec> * = nullptr>
staninline double my_dot_self(const EigVec &x, std::ostream *pstream__)
{
auto x_ref = stan::math::to_ref(x);
double sum = 0.0;
for (int i = 0; i < x.size(); ++i)
{
+= x_ref.coeff(i) * x_ref.coeff(i);
sum }
return sum;
}
template <typename EigVec, stan::require_eigen_vt<stan::is_var, EigVec> * = nullptr>
inline stan::math::var my_dot_self(const EigVec &v, std::ostream *pstream__)
{
// (1) put v into our memory arena
::arena_t<EigVec> arena_v(v);
stan// (2) calculate forward pass using
// (3) the .val() method for matrices of var types
::math::var res = my_dot_self(arena_v.val(), pstream__);
stan// (4) Place a callback for the reverse pass on the callback stack.
::math::reverse_pass_callback(
stan[res, arena_v]() mutable
{ arena_v.adj() += 2.0 * res.adj() * arena_v.val(); });
return res;
}
```

For more details about how to write C++ code using the Stan Math Library, see the Math library documentation at https://mc-stan.org/math/ or the paper at https://arxiv.org/abs/1509.07164.

## Special functions: RNGs, distributions, editing `target`

Some functions have special meanings in Stan and place additional requirements on their signatures if used in external C++.

- RNGs must end with
`_rng`

. They will be passed a “base RNG object” as the second to last argument, before the pointer to the ostream. We recommend making this a template, since it may change. This is currently a`stan::rng_t`

object (a type alias to`boost::rng::mixmax`

). - Functions which edit the
`target`

directly must end with`_lp`

and will be passed a reference to`lp__`

and a reference to a`stan::math::accumulator`

object as the final parameters before the ostream pointer. They are also expected to have a boolean template parameter`propto__`

which controls whether or not constant terms can be dropped. - Probability distributions must end with
`_lpdf`

or`_lpmf`

and will be passed a boolean template parameter`propto__`

which controls whether or not constant terms can be dropped.

## Footnotes

Details of programming in the Stan Math style are omitted from this section, it is presented only as an example↩︎