User-defined function causing Stan program to crash

General
1

Using the help provided here I have written user functions in C++ that are used in my likelihood calculation, and for calculating the derivatives of the likelihood.

I got an example working where the return type was univariate (double in my C++ source code).

I am now trying to implement another example where the return type is VectorXd.

I have managed to get my code to compile, but when I run the executable, there is some kind of run-time error:

 > ./fhn.exe sample data file=fhn.d
method = sample (Default)
  sample
    num_samples = 1000 (Default)
    num_warmup = 1000 (Default)
...
Aborted (core dumped)

Outside of Stan I have written a simple test program that just evaluates my user functions. This seems to work fine.

The next level up from my source code are wrappers that can take var arguments. Here are my wrapper definitions,

namespace stan{
  namespace math{

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

    // this is a template specialization for <double, double>
    template<>
    Eigen::Matrix<double, -1,1> spec_fun<double, double>(const Eigen::Matrix<double, -1,1>& 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<>
    Eigen::Matrix<var, -1, 1> spec_fun<double,var>(const Eigen::Matrix<double, -1,1>& p,
                                                   const Eigen::Matrix<var, -1, 1>& x,
                                                   std::ostream* pstream__) {
      vector<var> x_std(x.data(), x.data() + x.rows() * x.cols());
      Eigen::VectorXd a = value_of(x);
      Eigen::MatrixXd J = linear_sde::spec_fun_jacobian(p, a);
      Eigen::Matrix<var, -1, 1> result(a.rows());
      std::vector<double> grads(x.rows());
      for (int i = 0; i < J.rows(); ++i) {
        for (int j = 0; j < a.rows(); ++j)
          grads[j] = J(i, j);
        result(i) = precomputed_gradients(a(i), x_std, grads);
      }
      return result;
    }
  }
}

These wrappers were written using the guidance here under the heading of ‘Functions with more than one output’.

Note that in the guidance an Eigen::VectorXd is passed to the function precomputed_gradients(...). When I tried this I got a compiler error, so I now copy the data into an std::vector before passing it to precomputed_gradients(...).

My usual way of debugging issues like this is to use the Visual Studio debugger, and this typically gives me a call stack with line numbers. Is it worth trying to do this with a Stan program?

Do you have any other suggestions for how to debug Stan programs? Or what the specific issue with my program might be?

2 
result(i) = precomputed_gradients(a(i), x_std, grads);

In that, x_std and grads be the same length. Implementation is here if it helps: https://github.com/stan-dev/math/blob/develop/stan/math/rev/core/precomputed_gradients.hpp (they’re actually allocated the same space on lines 55-57)

The lingo is that a(i) is a function of all the x_std variables, and so the gradient of x_std is length x_std.size() (one for each var).

3

I have now got my code to compile and run without errors.

It turned out that there was a problem with my user function. But this problem was only an issue with the g++ compiler, which I hadn’t tested.

Here is my working source code wrapper. I think this fixes a couple of errors in the documentation here, and in my opinion it is more readable,

    Eigen::Matrix<var, -1, 1> spec_fun<double,var>(const Eigen::Matrix<double, -1,1>& pos,
                                                   const Eigen::Matrix<var,-1,1>& theta_eigen_var,
                                                   std::ostream* pstream__) {
      // extract parameter values from theta_eigen_var
      const var * theta_data = theta_eigen_var.data();
      int theta_size = theta_eigen_var.size();
      vector<var> theta_std_var(theta_data, theta_data + theta_size );
      VectorXd theta_val = value_of( theta_eigen_var );

      // evaluate f
      VectorXd f_val = linear_sde::spec_fun(pos, theta_val);
      // f_val_jacobian[i][j] is the partial derivative of f at pos[i] w.r.t. parameter j
      vector<vector<double>> f_val_jacobian = linear_sde::spec_fun_jacobian(pos, theta_val);

      int n_pos = f_val_jacobian.size(); // number of frequencies
      Eigen::Matrix<var, -1, 1> f_var_jacobian(n_pos);
      for (int i=0; i < n_pos; i++) {
        f_var_jacobian(i) = precomputed_gradients(f_val(i), theta_std_var, f_val_jacobian[i]);
      }
      return f_var_jacobian;
    }
4

Was it just the VectorXd thing or were there other things as well? If you could fix those errors while it’s still fresh on your mind that’d be super groovy. Wiki should be editable I think.

5

There was another thing - the example did not evaluate the function (only its derivatives). I tried editing the wiki, but when I did a git push I got a message saying Permission to stan-dev/math.wiki.git denied to pmaybank.

Editing in the browser didn’t work either for some reason.

6

Whoa. That looks like a GitHub problem that we can’t control / solve. The workaround is to edit the wiki through the web interface. Was it a major change to the wiki?

7

Yes, it’s fixed now, and it was some random GitHub error that one of their admins had to sort out. The Wiki has now been updated!

1 Like
8

Thanks!