What is the best way to add a user function that is used in likelihood calculation?

General
1

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,

  1. Write the function in the Stan modelling language and include it in the functions block of the .stan file.
  2. 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.
  3. 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.

2

With CmdStan, you have to define your C++ function in the stan::math namespace and include the appropriate headers. To tie in the derivatives, see

3

I have now defined my C++ function and its derivative inside the stan::math namespace, as suggested by @bgoodri . However I still get an error from the C++ linker complaining about an undefined reference,

/tmp/cczMg06L.o:main.cpp:... 
undefined reference to `fhn_model_namespace::spec_fun<double, stan::math::var>(...)'

What is the best way of making my function (spec_fun) visible to my Stan program?

So far the only way I can see is Option 3 in my original post.

In case it is helpful, here is the definition of my function inside the stan::math namespace,

namespace stan {
  namespace math {
    double spec_fun(const double& p, const Eigen::Matrix<double,-1,1>& x, std::ostream* pstream__) {
      return linear_sde::spec_fun(p, x);
    }

    var spec_fun(const double& p, const Eigen::Matrix<var,-1,1>& x, std::ostream* pstream__) {
      Eigen::Matrix<double,-1,1> a = value_of(x);
      double fa = linear_sde::spec_fun(p, a);
      vector<double> grad_fa = linear_sde::grad_spec_fun(p, a);
      vector<var> x_std(x.data(), x.data() + x.rows() * x.cols());
      return precomputed_gradients(fa, x_std, grad_fa);
    }
  }
}
4

This indicates that it’s looking for a templated function with two template parameters.

The code you posted has two functions with overloading that have 0 parameters. I think one way to get this to work is to do something like (not tested, may need to adjust some syntax):

template<typename T1, typename T2>
T1 spec_fun(const T1& p, const Eigen::Matrix<T2,-1,1>& x, std::ostream* pstream__) {
  throw std::logic_error("not implemented");  // this should never be called
}

// this is a template specialization for <double, double>
template<>
double spec_fun<double, double>(const double& p, const Eigen::Matrix<double,-1,1>& x, std::ostream* pstream__) {
  return linear_sde::spec_fun(p, x);
}

// this is a template specialization for <double, var>
template<>
var spec_fun<double, var>(const double& p, const Eigen::Matrix<var,-1,1>& x, std::ostream* pstream__) {
      Eigen::Matrix<double,-1,1> a = value_of(x);
      double fa = linear_sde::spec_fun(p, a);
      vector<double> grad_fa = linear_sde::grad_spec_fun(p, a);
      vector<var> x_std(x.data(), x.data() + x.rows() * x.cols());
      return precomputed_gradients(fa, x_std, grad_fa);
}

Hopefully that helps a bit. If you need more info, search for “function template specialization” and it should explain some more of the syntax.

1 Like
5

Like @syclik said, the functions you write need to be fully templated. If you have a function:

template<typename T1, typename T2>
T1 foo(const T2& arg) {
   ...
}

It needs to compile with T2 being a var but also a plain double. In the case T2 is a var, you’ll probably want to return a var, and likewise if T2 is a double you’ll probly want to return a double. To handle this you’ll want to do some fancy type checking return type stuff. I think the reason both versions are needed is sometimes the log density is evaluated without autodiff enabled, so that’s doubles.

The easiest way to get the right function signature (w/ return type calculations) is to define the function in a Stan model and compile that to C++ with the allow_undefined = TRUE flag on. From the C++ file that gets generated, you can see exactly the function signature that the compiler will be looking for.

This walks through custom gradients using the external C++ stuff for Rstan: https://cdn.rawgit.com/bbbales2/stancon_2018/718654d2/rus.html

But it looks like you understand how everything works, so maybe just skip to the “Eigenvalues in templated C++ (with Eigen)” part.

3 Likes
6

I have now managed to get my code to compile. As @syclik and @bbbales2 pointed out I needed to template my functions.

@syclik suggested using template specialization (a C++ feature that I was not previously aware of). This is nice because I don’t have to modify my source code.

The link that @bbbales2 posted was an example where all the templating is done in the source code. The text in the relevant section starts “if we have some fancy templated C++ code that solves part of our problem”.

I don’t think my code needs to be fancy or templated so I prefer @syclik’s solution of creating a fancy templated wrapper around my basic code.

In @syclik’s post there is no fancy type-checking of return arguments, which I think is an error (at least my compiler said it was). For completeness, here are my (working) function definitions,

namespace fhn_model_namespace{

  template<typename T1, typename T2>
  typename boost::math::tools::promote_args<T1, T2>::type
  spec_fun(const T1& p, const Eigen::Matrix<T2,-1,1>& x, std::ostream* pstream__) {
    throw std::logic_error("not implemented");  // this should never be called
  }

  // this is a template specialization for <double, double>
  template<>
  double spec_fun<double, double>(const double& p, const Eigen::Matrix<double,-1,1>& x, std::ostream* pstream__) {
    return linear_sde::spec_fun(p, x);
  }

  // this is a template specialization for <double, var>
  template<>
  var spec_fun<double, var>(const double& p, const Eigen::Matrix<var,-1,1>& x, std::ostream* pstream__) {
    Eigen::Matrix<double,-1,1> a = value_of(x);
    double fa = linear_sde::spec_fun(p, a);
    vector<double> grad_fa = linear_sde::grad_spec_fun(p, a);
    vector<var> x_std(x.data(), x.data() + x.rows() * x.cols());
    return precomputed_gradients(fa, x_std, grad_fa);
  }

}
3 Likes