Hi all,
I am trying to write a custom lpdf in C++. I started by trying to run the example from this thread, but I get an error:
Here is the C++ code:
#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);
}
And here is the R script that does the sampling:
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=4,
init=list(list(mu=init_mu, kappa=init_kappa)))
summary(m)
I get this error when I try to run the sampling function:
SAMPLING FOR MODEL 'b785e1b459e1a1bedfb490735b98118b' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.039457 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 394.57 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
[1] "Error in sampler$call_sampler(args_list[[i]]) : "
[2] " c++ exception (unknown reason)"
error occurred during calling the sampler; sampling not done
I don’t get any errors when running Stan code that doesn’t involve an external lpdf. I have also tried some of the Catalina-related suggestions discussed here to no avail:
- Operating System: Catalina 10.15.6
- RStan Version:
rstan_2.21.1 - Output of
writeLines(readLines(file.path(Sys.getenv("HOME"), ".R/Makevars"))):
# The following statements are required to use the clang4 binary
LDFLAGS= -L/usr/local/clang4/lib
# End clang4 inclusion statements
CC=/usr/local/clang4/bin/clang
CXX=/usr/local/clang4/bin/clang++
CXX1X=/usr/local/clang4/bin/clang++
CXX98=/usr/local/clang4/bin/clang++
CXX11=/usr/local/clang4/bin/clang++
CXX14=/usr/local/clang4/bin/clang++
CXX17=/usr/local/clang4/bin/clang++
- Output of
devtools::session_info("rstan"):
─ Session info ──────────────────────────────────────────────────────────────
setting value
version R version 4.0.3 (2020-10-10)
os macOS Catalina 10.15.6
system x86_64, darwin17.0
ui RStudio
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz America/New_York
date 2020-10-15
─ Packages ──────────────────────────────────────────────────────────────────
package * version date lib source
assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.0.2)
backports 1.1.10 2020-09-15 [1] CRAN (R 4.0.2)
BH 1.72.0-3 2020-01-08 [1] CRAN (R 4.0.2)
callr 3.5.1 2020-10-13 [1] CRAN (R 4.0.3)
checkmate 2.0.0 2020-02-06 [1] CRAN (R 4.0.2)
cli 2.1.0 2020-10-12 [1] CRAN (R 4.0.3)
colorspace 1.4-1 2019-03-18 [1] CRAN (R 4.0.2)
crayon 1.3.4 2017-09-16 [1] CRAN (R 4.0.2)
curl 4.3 2019-12-02 [1] CRAN (R 4.0.1)
desc 1.2.0 2018-05-01 [1] CRAN (R 4.0.2)
digest 0.6.25 2020-02-23 [1] CRAN (R 4.0.2)
ellipsis 0.3.1 2020-05-15 [1] CRAN (R 4.0.2)
evaluate 0.14 2019-05-28 [1] CRAN (R 4.0.1)
fansi 0.4.1 2020-01-08 [1] CRAN (R 4.0.2)
farver 2.0.3 2020-01-16 [1] CRAN (R 4.0.2)
ggplot2 * 3.3.2 2020-06-19 [1] CRAN (R 4.0.2)
glue 1.4.2 2020-08-27 [1] CRAN (R 4.0.2)
gridExtra 2.3 2017-09-09 [1] CRAN (R 4.0.2)
gtable 0.3.0 2019-03-25 [1] CRAN (R 4.0.2)
inline 0.3.16 2020-09-06 [1] CRAN (R 4.0.2)
isoband 0.2.2 2020-06-20 [1] CRAN (R 4.0.2)
jsonlite 1.7.1 2020-09-07 [1] CRAN (R 4.0.2)
labeling 0.3 2014-08-23 [1] CRAN (R 4.0.2)
lattice 0.20-41 2020-04-02 [1] CRAN (R 4.0.3)
lifecycle 0.2.0 2020-03-06 [1] CRAN (R 4.0.2)
loo 2.3.1 2020-07-14 [1] CRAN (R 4.0.2)
magrittr 1.5 2014-11-22 [1] CRAN (R 4.0.2)
MASS 7.3-53 2020-09-09 [1] CRAN (R 4.0.3)
Matrix 1.2-18 2019-11-27 [1] CRAN (R 4.0.3)
matrixStats 0.57.0 2020-09-25 [1] CRAN (R 4.0.2)
mgcv 1.8-33 2020-08-27 [1] CRAN (R 4.0.3)
munsell 0.5.0 2018-06-12 [1] CRAN (R 4.0.2)
nlme 3.1-149 2020-08-23 [1] CRAN (R 4.0.3)
pillar 1.4.6 2020-07-10 [1] CRAN (R 4.0.2)
pkgbuild 1.1.0 2020-07-13 [1] CRAN (R 4.0.2)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.0.2)
pkgload 1.1.0 2020-05-29 [1] CRAN (R 4.0.2)
praise 1.0.0 2015-08-11 [1] CRAN (R 4.0.2)
prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.0.2)
processx 3.4.4 2020-09-03 [1] CRAN (R 4.0.2)
ps 1.4.0 2020-10-07 [1] CRAN (R 4.0.2)
R6 2.4.1 2019-11-12 [1] CRAN (R 4.0.2)
RColorBrewer 1.1-2 2014-12-07 [1] CRAN (R 4.0.2)
Rcpp 1.0.5 2020-07-06 [1] CRAN (R 4.0.2)
RcppEigen 0.3.3.7.0 2019-11-16 [1] CRAN (R 4.0.2)
RcppParallel 5.0.2 2020-06-24 [1] CRAN (R 4.0.2)
rlang 0.4.8 2020-10-08 [1] CRAN (R 4.0.2)
rprojroot 1.3-2 2018-01-03 [1] CRAN (R 4.0.2)
rstan * 2.21.1 2020-07-08 [1] CRAN (R 4.0.2)
rstudioapi 0.11 2020-02-07 [1] CRAN (R 4.0.2)
scales 1.1.1 2020-05-11 [1] CRAN (R 4.0.2)
StanHeaders * 2.21.0-6 2020-08-16 [1] CRAN (R 4.0.2)
testthat 2.3.2 2020-03-02 [1] CRAN (R 4.0.2)
tibble 3.0.4 2020-10-12 [1] CRAN (R 4.0.3)
utf8 1.1.4 2018-05-24 [1] CRAN (R 4.0.2)
V8 3.2.0 2020-06-19 [1] CRAN (R 4.0.2)
vctrs 0.3.4 2020-08-29 [1] CRAN (R 4.0.2)
viridisLite 0.3.0 2018-02-01 [1] CRAN (R 4.0.1)
withr 2.3.0 2020-09-22 [1] CRAN (R 4.0.2)
[1] /Library/Frameworks/R.framework/Versions/4.0/Resources/library
Thanks in advance.