Includes in USER_HEADER

Developers
1

I’m having some trouble figuring out how includes work in the USER_HEADER case: I want to implement the Stan-C bridge in an external header but variations of including Eigen result in different errors. To make the question clearer, here’s a gist with all the files laid out,

where the errors are of form

stan/lib/stan_math/stan/math/rev/fun/gp_exp_quad_cov.hpp:94:33: error: ‘const class stan::math::arena_matrix<Eigen::Matrix<stan::math::var_value<double>, -1, 1> >’ has no member named ‘val’; did you mean ‘eval’?
   94 |           adjsigma += (cov_diag.val().array() * cov_diag.adj().array()).sum();

Thanks in advance.

2

The errors:

has no member named ‘val’; did you mean ‘eval’?

Are occurring because your user header is including Eigen headers:

#include <Eigen/Eigen>
#include "support.h"

double
usersum(const Eigen::MatrixXd m, std::ostream* pstream__) {
        const double *x = m.data();
        int n = m.rows() * m.cols();
        return support_sum(n, x);
}

Stan uses Eigen’s plugin system for handling matrices/vectors of autodiff types. However, this requires that the plugins are defined before any Eigen headers are included, otherwise they don’t get loaded. Instead of including the Eigen headers directly, try including the Stan Math headers instead, as these make sure that the plugins are appropriately defined before the Eigen headers get loaded:

#include <stan/math.hpp>
#include "support.h"

double usersum(const Eigen::MatrixXd m, std::ostream* pstream__) {
        const double *x = m.data();
        int n = m.rows() * m.cols();
        return support_sum(n, x);
}

Also, I should mention that your current usersum definition will only work in the transformed data and generated quantities blocks, as it is only defined for double inputs. If you’re only working with basic arithmetic operations, then you can instead define the headers in a templated fashion, so that the operations can be handled through autodiff:

template <typename ScalarT>
ScalarT support_sum(int n, const ScalarT* els) {
        ScalarT acc = 0.0;
        for (int i=0; i<n; i++) acc += els[i];
        return acc;
}

#include <stan/math.hpp>
#include "support.h"

template <typename MatrixT, typename ScalarT = scalar_type_t<MatrixT>>
ScalarT usersum(const MatrixT m, std::ostream* pstream__) {
        const ScalarT *x = m.data();
        int n = m.rows() * m.cols();
        return support_sum(n, x);
}
1 Like
3

If your external headers are going to be using functions from other libraries (e.g., Boost), then the above templating approach will likely fail, as Stan can’t auto-differentiate through functions in external libraries (in most cases). In that case, you would also need to define the gradients wrt each input for the function.

If this is something that you need, I put together an example of how to do this for a Boost function in this post (and rest of thread): Gamma distribution quantile and inverse quantile function - #28 by andrjohns

1 Like
4

Thanks, this is what I was missing and couldn’t find.

Thanks, I have the gradients defined in another C file which is why I modeled this example like that. I do wish it were easier to scaffold the FFI stuff like this.

I will try to use reverse_pass_callback for this once I’ve got the gradients from the C function, as outlined here.

5

Trying to follow the examples linked to so far, I’ve tried to implement the usersum, as follows

#include <stan/math.hpp>

using namespace stan;
using namespace std;

template <typename M, require_matrix_t<M>* = nullptr>
inline auto usersum(const M& m, ostream *pstream__) {
        using ret_t = return_var_matrix_t<M>;
        const double *x = m.data();
        int n = m.rows() * m.cols();
        // auto ret = support_sum(n, x);
        arena_t<ret_t> one(m.rows(), m.cols());
        one = 1.0;
        arena_t<ret_t> ret = m.sum();
        arena_t<ret_t> arena_m = m;
        reverse_pass_callback([ret, arena_m, one]() mutable {
                arena_m.adj() += one;
        });
        return ret_t(ret);
}

but there is a pile of link errors of the sort

 undefined reference to `boost::math::tools::promote_args<stan::value_type<Eigen::Matrix<stan::math::var_value<double, void>, -1, -1, 0, -1, -1>, void>::type, float, float, float, float, float>::type model_model_namespace::usersum<Eigen::Matrix<stan::math::var_value<double, void>, -1, -1, 0, -1, -1> >(Eigen::Matrix<stan::math::var_value<double, void>, -1, -1, 0, -1, -1> const&, std::ostream*)'

I’m not sure if this is a namespace issue or that the types are wrong or missing an implementation without var? I’m also trying to understand this from the generalized inverse PR but haven’t got it yet.

edit 1 if I namespace the above function in the model_model_namespace, it seems to resolve the link error, but there’s an ambiguity in the overload which I also don’t understand.

/home/duke/src/gist-stan-user-func-ext-c/model.hpp:96:31: error: call of overloaded ‘usersum(Eigen::Matrix<stan::math::var_value<double>, -1, -1>&, std::ostream*&)’ is ambiguous
   96 |         lp_accum__.add(usersum(m, pstream__));
      |                        ~~~~~~~^~~~~~~~~~~~~~
In file included from <command-line>:
/home/duke/src/gist-stan-user-func-ext-c/support.hpp:9:1: note: candidate: ‘stan::promote_args_t<typename stan::value_type<T, void>::type> model_model_namespace::usersum(const M&, std::ostream*) [with M = Eigen::Matrix<stan::math::var_value<double>, -1, -1>; stan::require_matrix_t<T>* <anonymous> = 0; stan::promote_args_t<typename stan::value_type<T, void>::type> = stan::math::var_value<double>; std::ostream = std::basic_ostream<char>]’
    9 | usersum(const M& m, ostream *pstream__) {
      | ^~~~~~~
/home/duke/src/gist-stan-user-func-ext-c/model.hpp:28:1: note: candidate: ‘stan::promote_args_t<typename stan::value_type<T, void>::type> model_model_namespace::usersum(const T0__&, std::ostream*) [with T0__ = Eigen::Matrix<stan::math::var_value<double>, -1, -1>; stan::promote_args_t<typename stan::value_type<T, void>::type> = stan::math::var_value<double>; std::ostream = std::basic_ostream<char>]’
   28 | usersum(const T0__& m_arg__, std::ostream* pstream__) ;
      | ^~~~~~~

edit 2 if I remove the require_matrix_t<M>* = nullptr then I get errors relevant to the content of the function, I’ll try to sort those out before another question…

edit 3 narrowing it down it requires two implementations for double and … something else, so now my user header looks like this

#include <stan/math.hpp>
#include <stan/math/rev/core/reverse_pass_callback.hpp>

namespace model_model_namespace {

using namespace stan;
using namespace std;

double usersum(const Eigen::MatrixXd& m, ostream *pstream__) {
        return m.sum();
}

template <typename M>
promote_args_t<value_type_t<M>>
usersum(const M& m, ostream *pstream__) {
        using value_t = value_type_t<M>;
        using ret_type = promote_var_matrix_t<M,M>;
        value_t ret(m.sum());
        arena_t<Eigen::MatrixXd> one(m.rows(), m.cols());
        one.fill(1.0);
        arena_t<M> arena_m(m);
        math::reverse_pass_callback([ret, arena_m, one]() mutable {
                arena_m.adj() += one;
        });
        return ret;
}
}

(one is just a placeholder for an externally computed gradient) This compiles successfully.

edit 5 it seems value_t ret(m.sum()) is doing AD w/o being obvious and a value_of is required,

        value_t ret(value_of(m).sum());
        arena_t<M> arena_m(m);
        math::reverse_pass_callback([ret, arena_m]() mutable {
                arena_m.adj().array() += ret.adj();
        });

This compiles and passes ./model diagnose.

6

Just to follow up on the example since the original goal was to interface with external C code, instead of computing with Eigen types, one needs to extract it to something C-ish, compute with it, and return something Stan wants,

        Eigen::MatrixXd mv = value_of(m);
        const double *md = mv.data();
        value_t ret(support_sum(m.size(), md));

then, in the backwards pass, a suitable adjoint Eigen type should be constructed and then filled by the C code before incrementing the AD adjoint,

        math::reverse_pass_callback([aret, arena_m]() mutable {
                Eigen::MatrixXd adj(arena_m.rows(), arena_m.cols());
                support_sum_adj(adj.size(), aret.adj(), adj.data());
                arena_m.adj() += adj;
        });

(or at least, I couldn’t figure out how to update arena_m.adj() in-place w/o a segfault). The whole working example is updated in the gist linked above.