I have written a function in C++ to evaluate a spectral density that I would then like to use in a likelihood calculation. I have also written a function that evaluates the derivatives of the spectral density with respect to model parameters.

I am aware of 3 different ways of including a user function in a Stan program,

- Write the function in the Stan modelling language and include it in the functions block of the .stan file.
- Declare the function in the functions block of the .stan file, define the function in a separate C++ header file and add appropriate options to the make command, e.g.,
`STANCFLAGS=--allow_undefined USER_HEADER=spec_fun.hpp`

in CmdStan. See http://mc-stan.org/rstan/articles/external.html for more details. - Add a new function to the Stan Math library using the guidance here and here.

I have ruled out Option 1 because my calculation involves complex numbers, which I believe the Stan programming language does not support.

The guidance in the link for Option 2 states *“it is possible to use an external C++ file to handle the gradient of a function analytically as opposed to using Stan’s autodifferentiation capabilities”*. However this statement is not illustrated by any examples, as far as I can see.

I have tried Option 2, but have had difficulty getting my code to compile. I get the following (redacted) error message,

```
g++ -Wall [...] -O3 -o [...]/stan-spec/fhn.exe src/cmdstan/main.cpp -include [...]/stan-spec/fhn.hpp -include [...]/stan-spec/USER_HEADER.hpp [...]
/tmp/ccC7WbsZ.o:main.cpp:(.text$_ZNK19fhn_model_namespace9fhn_model8log_prob[...]: undefined reference to `boost::math::tools::promote_args<[...]>::type fhn_model_namespace::spec_fun<[...]>
/tmp/ccC7WbsZ.o:main.cpp:(.text$_ZNK19fhn_model_namespace9fhn_model8log_prob[...]: relocation truncated to fit: R_X86_64_PC32 against undefined symbol `boost::math::tools::promote_args<[...]>::type fhn_model_namespace::spec_fun<[...]>(...)
/tmp/ccC7WbsZ.o:main.cpp:(.text$_ZNK19fhn_model_namespace9fhn_model8log_prob[...]: undefined reference to `boost::math::tools::promote_args<[...]>::type fhn_model_namespace::spec_fun<[...](...)
```

Option 3 seems a bit involved, but I am happy to give it a try if it is the best way of doing it?

Here is the function I would like to include,

```
real spec_fun(const real& p, const real& omega0, const real& zeta, const real& sd_in, std::ostream* pstream__) {
real spec_dens;
// define Jacobian for simple harmonic oscillator
MatrixXd J(2,2);
J(0,0) = 0; J(0,1) = 1;
J(1,0) = -10E3 * omega0; J(1,1) = -2 * zeta * omega0;
// calculate eigendecomposition
EigenSolver<MatrixXd> es(J);
MatrixXcd R = es.eigenvectors();
VectorXcd lambda_vec = es.eigenvalues();
MatrixXcd L = R.inverse();
// define variables for components of eigenvectors that contribute to spectral density
VectorXcd r = R.row(0), l = L.col(1);
complex<real> c = 2*M_PI * p * 1i;
VectorXcd diag_vec = (VectorXcd::Constant(2,c) - lambda_vec).cwiseInverse();
complex<real> trans_func = (r.transpose() * diag_vec.asDiagonal() * l)(0);
// evaluate spectral density
real spec_dens = pow(abs(trans_func),2) * pow(sd_in,2.0);
return spec_dens;
}
```

Apologies if this overlaps with material on some other posts. I think it is useful to bring together information on different ways of defining a user function, so that best practices in this area are more obvious to people.