Any simple working example of external C++ lpdf function?

I am trying to write a custom distribution function (lpdf) in C++ that is called from R with rstan. However, when I compile and fit it, I am getting some errors. Here the contents of my files.

Rscript:

library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

model_code = "
functions {
  real custom_function_lpdf(real y, real A, real B);
}
data {
    int N;
    real y[N];
}
parameters {
    real A;
    real B; 
}
model {
  for(i in 1:N){
    y[i] ~ custom_function(A, B);
  }
}
"

data <- list(N = 1, A = 1, B = 1)

model <- stan_model(model_code = model_code,
                    allow_undefined = TRUE,
                    includes = paste0('\n#include "', file.path(getwd(), 'custom_function_lpdf.hpp'), '"\n'))

fit <- sampling(model, data = data, iter=1000, chains=1)

print(fit)

custom_function_lpdf.hpp:

#include <boost/math/tools/promotion.hpp>

using namespace boost::math;
using namespace std;

template <typename T_y, typename T_A, typename T_B>
typename boost::math::tools::promote_args<T_y, T_A, T_B>::type
parabola_function_lpdf(const T_y& y,
                       const T_A& A,
                       const T_B& B,
                       std::ostream* pstream__) { 
  return A*B*y;
}

I get the following error:

Error in dyn.load(libLFile) : 
  unable to load shared object '/tmp/RtmpKiWiMx/file2dee16490e35f8.so':
  /tmp/RtmpKiWiMx/file2dee16490e35f8.so: undefined symbol: _ZN62model2dee1627232579_76aac3ff274b7f875f40dd19472b96f5_namespace20custom_function_lpdfILb0EdddEEN5boost4math5tools12promote_argsIT0_T1_T2_fffE4typeERKS5_RKS6_RKS7_PSo

The C++ code seems to compile but the error is related to linking. What am I doing wrong? Is there any working example of external custom C++ distribution function?

For example, if I want to implement the following std_normal_lpdf function from Stan as a custom_function_lpdf, how can I do it?

#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/err.hpp>
#include <stan/math/prim/fun/constants.hpp>
#include <stan/math/prim/fun/size.hpp>
#include <stan/math/prim/fun/size_zero.hpp>
#include <stan/math/prim/fun/value_of.hpp>
#include <stan/math/prim/functor/operands_and_partials.hpp>

namespace stan {
namespace math {

template <bool propto, typename T_y>
return_type_t<T_y> std_normal_lpdf(const T_y& y) {
  static const char* function = "std_normal_lpdf";
  using T_partials_return = partials_return_t<T_y>;
  check_not_nan(function, "Random variable", y);

  if (size_zero(y)) {
    return 0.0;
  }
  if (!include_summand<propto, T_y>::value) {
    return 0.0;
  }

  T_partials_return logp(0.0);
  operands_and_partials<T_y> ops_partials(y);

  scalar_seq_view<T_y> y_vec(y);
  size_t N = stan::math::size(y);

  for (size_t n = 0; n < N; n++) {
    const T_partials_return y_val = value_of(y_vec[n]);
    logp += y_val * y_val;
    if (!is_constant_all<T_y>::value) {
      ops_partials.edge1_.partials_[n] -= y_val;
    }
  }
  logp *= -0.5;
  if (include_summand<propto>::value) {
    logp += NEG_LOG_SQRT_TWO_PI * N;
  }

  return ops_partials.build(logp);
}

template <typename T_y>
inline return_type_t<T_y> std_normal_lpdf(const T_y& y) {
  return std_normal_lpdf<false>(y);
}

}  // namespace math
}  // namespace stan
3 Likes

@bgoodri Do you know what’s going wrong here?

Good question, I’m sure there are examples of implementing a distribution in C++ but we should have a place to find them more easily.

Here’s a thread I found where @maxbiostat was trying to implement a distribution in external C++ but that’s a long thread and I’m not sure how much is related to the linking issue you’re getting; Problem implementing Zipf (power law) distribution -- external C++

1 Like

I think you’re missing a template argument for propto.

You can look at the code generated by rstan for this model by doing:

model_source = stanc(model_code = model_code, allow_undefined = TRUE)
cat(model_source$cppcode)

Scrolling through this there is a function signature:

template <bool propto, typename T0__, typename T1__, typename T2__>
typename boost::math::tools::promote_args<T0__, T1__, T2__>::type
custom_function_lpdf(const T0__& y,
                     const T1__& A,
                     const T2__& B, std::ostream* pstream__);

You need to match that definition exactly with your code (so I think that means adding the propto template argument). If you don’t match that function exactly, the compiler will just think this function is defined elsewhere, point to it instead of your custom code, and then fail at link.

This stuff is hard to debug. I would like a more formal way to interfaces with other languages so we could have better errors.

4 Likes

2 Likes

@jonah and @bgoodri thank you for the links. Before posting here, I did quite a bit of searching and came across these links. I also went through the complete thread therein but they only discuss calling C++ functions not lpdf. There was one with lpdf but it is implemented in Stan code, not C++.

As pointed out by @bbbales2, C++ lpdf requires the template argument for propto. I have not found any example of it. After adding it, it works! Thank you for the tip to use model_source$cppcode to view the generated C++ source.

I just have another issue. What is the best way to debug the code in the external C++ function? Right now I am not getting the full error message from the compiler. For example, when I purposely added a undefined variable, ys in the code below:

#include <boost/math/tools/promotion.hpp>

using namespace boost::math;
using namespace std;

template <bool propto, typename T0__, typename T1__, typename T2__>
typename boost::math::tools::promote_args<T0__, T1__, T2__>::type
custom_function_lpdf(const T0__& y,
                         const T1__& A,
                         const T2__& B, std::ostream* pstream__) {
  return A*B*ys;
}

The error message that I am getting is:

Error in compileCode(f, code, language = language, verbose = verbose) : 
  /home/satya/R/x86_64-pc-linux-gnu-library/4.0/StanHeaders/include/stan/math/rev/mat/functor/adj_jac_apply.hpp:535:10:   required from ‘void stan::math::adj_jac_vari<F, Targs>::chain() [with F = stan::math::internal::softmax_op; Targs = {Eigen::Matrix<stan::math::var, -1, 1, 0, -1, 1>}]’/home/satya/R/x86_64-pc-linux-gnu-library/4.0/StanHeaders/include/stan/math/rev/mat/functor/adj_jac_apply.hpp:531:8:   required from here/home/satya/R/x86_64-pc-linux-gnu-library/4.0/RcppEigen/include/Eigen/src/Core/DenseCoeffsBase.h:650:34: warning: ignoring attributes on template argument ‘Eigen::internal::packet_traits<double>::type’ {aka ‘__vector(4) double’} [-Wignored-attributes]  650 |   return internal::first_aligned<int(unpacket_traits<DefaultPacketType>::alignment),Derived>(m);      |                                  ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~make: *** [/usr/lib/R/etc/Makeconf:176: file17a2171b56b63b.o] Error 1

It only provides the warning part of the compiler message. How can I view the entire error message to help me debug?

1 Like

@bgoodri would know for sure, but try adding verbose = TRUE to your stan or stan_model call.

(Edit: I’m not really sure, just saw a verbose = verbose thing in the output and it looks like verbose is an argument for more output in the docs for stan)

Thanks but adding verbose = TRUE to both stanc() and stan_model() functions gives the same limited warning message. Any other ideas?

I found the source of the truncated error messages in inline’s compileCode function (https://github.com/eddelbuettel/inline/blob/master/R/cfunction.R).

It truncates the error messages if they are more than the warning.length size (default 1000). Setting the size to the max warning.length size of 8100 doesn’t help because it is not enough. The stop() function also truncates the error messages if the length is more than the warning.length.

One workaround I found is to change the following lines near the end of compileCode function:

    if ( verbose ) {
      writeLines(errmsg)
    }
    else if ( sum(nchar(errmsg)) > getOption("warning.length") ) {
      stop(tail(errmsg))
    }
    else stop(errmsg)

Another is to comment out unlink( errfile ) and set the tempdir to the current working directory using https://rdrr.io/github/s-u/unixtools/man/set.tempdir.html. The complete error messages can then be read from the error file.

1 Like

Now that I can see the complete error messages, I have managed to build a working lpdf function.
For future reference to anyone wanting to build a custom C++ lpdf function, here is a working example of custom_von_mises_lpdf with minimal modification of the original von_mises_lpdf C++ code:

Rscript:

library(rstan)
library(circular)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

mu <- 2*atan(-1.0) 
kappa <- exp(0.4)
N <- 10000
K <- 1
df <- data.frame(x=rep(1, N))
X <- model.matrix(~0+., data=df)
y <- rvonmises(N, mu, kappa)

y[y > pi] <- -(2*pi - y[y > pi])

model_code <- "
functions{
  real custom_von_mises_lpdf(real y, real mu, row_vector kappa);
}
data {
  int<lower=0> N; //the number of observations
  int<lower=0> K; //the number of columns in the model matrix
  matrix[N,K] X; //the model matrix
  vector[N] y; //the response
}
parameters {
  real kappa; //the scale parameter
  real mu; //the loc parameter
}
model {
  for(i in 1:N)
    y[i] ~ custom_von_mises(2*atan(mu), exp(X[i,]*kappa));
}"

init_mu = mean(y)
init_kappa = 1

model_source = stanc(model_code = model_code, allow_undefined = TRUE,
                     verbose = TRUE)

model <- stan_model(model_code = model_code,
        allow_undefined = TRUE,
        verbose = TRUE, 
        includes = paste0('\n#include "',
                          file.path(getwd(), 'custom_von_mises_lpdf.hpp'), '"\n'))

m <- sampling(model, data = list(X=X, K=K, N=N, y=y), chains = 1, cores=16,
          init=list(list(mu=init_mu, kappa=init_kappa))) 

summary(m)

custom_von_mises_lpdf.hpp:

#include <stan/model/model_base.hpp>
#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/scal/err/check_consistent_sizes.hpp>
#include <stan/math/prim/scal/err/check_finite.hpp>
#include <stan/math/prim/scal/err/check_positive_finite.hpp>
#include <stan/math/prim/scal/fun/log_modified_bessel_first_kind.hpp>
#include <stan/math/prim/scal/fun/modified_bessel_first_kind.hpp>
#include <stan/math/prim/scal/fun/constants.hpp>
#include <stan/math/prim/scal/fun/value_of.hpp>
#include <stan/math/prim/scal/fun/size_zero.hpp>
#include <cmath>

#include <boost/math/tools/promotion.hpp>

using namespace boost::math;
using namespace std;
using namespace stan;

template <bool propto, typename T_y, typename T_loc, typename T_scale>
typename boost::math::tools::promote_args<T_y, T_loc, T_scale>::type
custom_von_mises_lpdf(T_y const& y,
                      T_loc const& mu,
                      Eigen::Matrix<T_scale, 1, Eigen::Dynamic> const& kappa,
                      std::ostream* pstream__) {

  using T_partials_return = partials_return_t<T_y, T_loc, T_scale>;

  if (size_zero(y, mu, kappa)) {
    return 0.0;
  }

  using std::log;

  T_partials_return logp = 0.0;

  if (!include_summand<propto, T_y, T_loc, T_scale>::value) {
    return logp;
  }

  const bool y_const = is_constant_all<T_y>::value;
  const bool mu_const = is_constant_all<T_loc>::value;
  const bool kappa_const = is_constant_all<T_scale>::value;

  const bool compute_bessel0 = include_summand<propto, T_scale>::value;
  const bool compute_bessel1 = !kappa_const;
  const double TWO_PI = 2.0 * pi();

  scalar_seq_view<T_y> y_vec(y);
  scalar_seq_view<T_loc> mu_vec(mu);
  scalar_seq_view<Eigen::Matrix<T_scale, 1, Eigen::Dynamic>> kappa_vec(kappa);

  VectorBuilder<true, T_partials_return, T_scale> kappa_dbl(length(kappa));
  VectorBuilder<include_summand<propto, T_scale>::value, T_partials_return,
                T_scale>
      log_bessel0(length(kappa));
  for (size_t i = 0; i < length(kappa); i++) {
    kappa_dbl[i] = value_of(kappa_vec[i]);
    if (include_summand<propto, T_scale>::value) {
      log_bessel0[i]
          = log_modified_bessel_first_kind(0, value_of(kappa_vec[i]));
    }
  }

  operands_and_partials<T_y, T_loc, Eigen::Matrix<T_scale, 1, Eigen::Dynamic>>
    ops_partials(y, mu, kappa);

  size_t N = max_size(y, mu, kappa);

  for (size_t n = 0; n < N; n++) {
    const T_partials_return y_ = value_of(y_vec[n]);
    const T_partials_return y_dbl = y_ - floor(y_ / TWO_PI) * TWO_PI;
    const T_partials_return mu_dbl = value_of(mu_vec[n]);

    T_partials_return bessel0 = 0;
    if (compute_bessel0) {
      bessel0 = modified_bessel_first_kind(0, kappa_dbl[n]);
    }
    T_partials_return bessel1 = 0;
    if (compute_bessel1) {
      bessel1 = modified_bessel_first_kind(-1, kappa_dbl[n]);
    }
    const T_partials_return kappa_sin = kappa_dbl[n] * sin(mu_dbl - y_dbl);
    const T_partials_return kappa_cos = kappa_dbl[n] * cos(mu_dbl - y_dbl);

    if (include_summand<propto>::value) {
      logp -= LOG_TWO_PI;
    }
    if (include_summand<propto, T_scale>::value) {
      logp -= log_bessel0[n];
    }
    if (include_summand<propto, T_y, T_loc, T_scale>::value) {
      logp += kappa_cos;
    }

    if (!y_const) {
      ops_partials.edge1_.partials_[n] += kappa_sin;
    }
    if (!mu_const) {
      ops_partials.edge2_.partials_[n] -= kappa_sin;
    }
    if (!kappa_const) {
      ops_partials.edge3_.partials_[n]
          += kappa_cos / kappa_dbl[n] - bessel1 / bessel0;
    }
  }
  return ops_partials.build(logp);
}

If you find any issues with the code above, please let me know. I have only started using Stan/rstan a few weeks ago. Cheers!

3 Likes

This was partially fixed by me in the latest inline on CRAN.

Ben, I am having this issue of truncated error message in the latest inline. I think it may be caused by your pull request: https://github.com/eddelbuettel/inline/pull/12.

I can only see the warning message at the end, not the full error message when there is an error in the external C++ function, even when I set verbose = TRUE and warning.length = 8170.

It should print the tail of the messages when it exceeds the line length.

Yes it prints the tail of the messages, but in my case the error is somewhere in the middle followed by warning messages. I end up only seeing the warning messages. You can reproduce this issue with my model above.

That compiles for me, and I assume it would run if y were defined.

Thanks for figuring this out and posting what worked!

Sorry Ben, if I was not clear. Below are the steps to reproduce the problem of truncated error messages.

  1. Rscript:
library(rstan)
library(circular)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
options(warning.length = 8170)

mu <- 2*atan(-1.0) 
kappa <- exp(0.4)
N <- 500
K <- 1
df <- data.frame(x=rep(1, N))
X <- model.matrix(~0+., data=df)
y <- rvonmises(N, mu, kappa)

y[y > pi] <- -(2*pi - y[y > pi])

model_code <- "
functions{
  real custom_von_mises_lpdf(real y, real mu, row_vector kappa);
}
data {
  int<lower=0> N; //the number of observations
  int<lower=0> K; //the number of columns in the model matrix
  matrix[N,K] X; //the model matrix
  vector[N] y; //the response
}
parameters {
  real kappa; //the scale parameter
  real mu; //the loc parameter
}
model {
  for(i in 1:N)
    y[i] ~ custom_von_mises(2*atan(mu), exp(X[i,]*kappa));
}"

init_mu = mean(y)
init_kappa = 1

model_source = stanc(model_code = model_code, allow_undefined = TRUE,
                     verbose = TRUE)

model <- stan_model(model_code = model_code,
        allow_undefined = TRUE,
        verbose = TRUE, 
        includes = paste0('\n#include "',
                          file.path(getwd(), 'custom_von_mises_lpdf.hpp'), '"\n'))

m <- sampling(model, data = list(X=X, K=K, N=N, y=y), chains = 1, cores=1,
          init=list(list(mu=init_mu, kappa=init_kappa))) 

summary(m)
  1. custom_von_mises_lpdf.hpp:
#include <boost/math/tools/promotion.hpp>

using namespace boost::math;
using namespace std;
using namespace stan;

template <bool propto, typename T_y, typename T_loc, typename T_scale>
typename boost::math::tools::promote_args<T_y, T_loc, T_scale>::type
custom_von_mises_lpdf(T_y const& y,
                      T_loc const& mu,
                      Eigen::Matrix<T_scale, 1, Eigen::Dynamic> const& kappa,
                      std::ostream* pstream__) {
  return ys;
}

The output with latest inline from cran, which does not show the actual error message:output.txt (43.9 KB)

The output with my workaround below of the inline’s compileCode function:output_with_workaround.txt (1.3 MB)

    if ( verbose ) {
      writeLines(errmsg)
    }
    else if ( sum(nchar(errmsg)) > getOption("warning.length") ) {
      stop(tail(errmsg))
    }
    else stop(errmsg)

As you can see, the output with the current version of inline truncates the actual error message. There is no way to know what the error is. With my workaround, the error message is shown entirely.

1 Like

I get

error: ys was not declared in this scope

That’s weird, the error doesn’t show up on my system. I am using Ubuntu 20.04 with gcc 9.3.0.

You can try it with CXX14FLAGS = -O3 -g0 -Wno-ignored-attributes -Wno-deprecated-declarations -Wno-unknown-pragmas in ~/.R/Makevars .

Yes suppressing the warnings with those flags does work. Thanks.