2.13 リリースから、Stan プログラムで外部 C++ コードを使用することがより簡単になった。この vignette は、その方法を簡単に説明する。
以下のようなフィボナッチ数を扱うスタンプログラムが（一部）あるとする。
functions {
int fib(int n);
int fib(int n) {
if (n <= 0) reject("n must be positive");
return n <= 2 ? 1 : fib(n - 1) + fib(n - 2);
}
}
model {} // use the fib() function somehow2行目では、fib 関数を再帰的に呼び出すために、定義される前に宣言している。
再帰的でない関数については、定義する前に宣言する必要はないが、宣言した方が有利な場合もある。例えば、私はよく、#include "file" の仕組みを使って、邪魔なだけである複雑なユーティリティ関数の定義を隠したいことがある。
functions {
real complicated(real a, real b, real c, real d, real e, real f, real g);
#include "complicated.stan"
}
model {} // use the complicated() function somehowこのStanプログラムは、デフォルトの stanc 関数（ sampling と stan が内部で呼び出す）ではなく、rstan パッケージの stanc_builder 関数を使ってパースする必要がある。
フィボナッチの例に戻ると、fib 関数を Stan 言語で定義する必要はない。なぜなら、declared であるが defined ではない関数を持つ Stan プログラムは、C++ ツールチェインの標準機能を使って、C++ で関数定義を提供することができるためである。例えば、このプログラムはデフォルトでパーサエラーを発生させる。
mc <-
'
functions { int fib(int n); }
model {} // use the fib() function somehow
'
try(stan_model(model_code = mc, model_name = "parser_error"), silent = TRUE)しかし、stan_model 関数に allow_undefined と includes の引数を指定し、指定した C++ ヘッダーファイルで fib 関数を定義すると、パース・コンパイルされる。
stan_model(model_code = mc, model_name = "external", allow_undefined = TRUE,
includes = paste0('\n#include "',
file.path(getwd(), 'fib.hpp'), '"\n'))Stan プログラムの C++ 表現は、一時ディレクトリに書き込まれ、コンパイルされるため、includes の引数を指定するのは少し厄介である。したがって、includes 引数には fib.hpp ファイルへのフルパスを指定する必要があり、この場合は作業ディレクトリにある。また、パスは二重引用符で囲む必要がある。そのため、paste0 関数の個別の引数では一重引用符が使用され、二重引用符は文字通りに解釈される。最後に、includes の引数には改行文字を含める必要がある ( "\n " ) at the start and end. It is possible to specify multiple paths using additional newline characters or include a" 他の C++ ヘッダーファイルへの #include 指示を含む “meta-header” ファイル。
引数 includes の結果は、C++ ファイルの行末に直接挿入される（次の行の before に直接挿入される CmdStan とは対照的）。
#include <stan/model/model_header.hpp>
namespace some_namespace {
using std::istream;
using std::string;
using std::stringstream;
using std::vector;
using stan::io::dump;
using stan::math::lgamma;
using stan::model::prob_grad;
using namespace stan::math;
typedef Eigen::Matrix<double,Eigen::Dynamic,1> vector_d;
typedef Eigen::Matrix<double,1,Eigen::Dynamic> row_vector_d;
typedef Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> matrix_d;
static int current_statement_begin__;
stan::io::program_reader prog_reader__() {
stan::io::program_reader reader;
reader.add_event(0, 0, "start", "model181f4ea8e0c2_mc");
reader.add_event(3, 3, "end", "model181f4ea8e0c2_mc");
return reader;
}
// various function declarations and / or definitions
#include "/full/path/to/fib.hpp"したがって、fib.hpp ファイル内の fib 関数の定義は、特定の名前空間（デフォルトではランダムな文字列）で囲む必要はない。メタインクルード」stan/model/model_header.hpp ファイルは、次のように読み取る。
#ifndef STAN_MODEL_MODEL_HEADER_HPP
#define STAN_MODEL_MODEL_HEADER_HPP
#include <stan/math.hpp>
#include <stan/io/cmd_line.hpp>
#include <stan/io/dump.hpp>
#include <stan/io/reader.hpp>
#include <stan/io/writer.hpp>
#include <stan/lang/rethrow_located.hpp>
#include <stan/model/model_base.hpp>
#include <stan/model/model_base_crtp.hpp>
#include <stan/model/prob_grad.hpp>
#include <stan/model/indexing.hpp>
#include <stan/services/util/create_rng.hpp>
#include <boost/random/additive_combine.hpp>
#include <boost/random/linear_congruential.hpp>
#include <cmath>
#include <cstddef>
#include <fstream>
#include <iostream>
#include <sstream>
#include <stdexcept>
#include <utility>
#include <vector>
#endifしたがって、fib.hpp ファイル内の fib 関数の定義は、Stan Math Library の任意の関数（関数呼び出しの前に stan::math:: を付ける必要はない）、Eigen 行列代数ライブラリのクラスに対するいくつかの typedef、さらにストリーム、例外などを利用でき、 対応するヘッダーファイルについて気にかける必要はない。それでも、外部の C++ ファイルには、例えばクラス定義を持ち込む追加のインクルードディレクティブが含まれている場合がある。
次に、C++の行を含むfib.hppファイルを見てみよう。
int fib(const int&n, std::ostream* pstream__) {
if (n <= 0) {
stringstream errmsg;
errmsg << "n must be positive";
throw std::domain_error(errmsg.str());
}
return n <= 1 ? 1 : fib(n - 1, 0) + fib(n - 2, 0);
}このC++の関数は、本質的には、先の Stan 言語のユーザー定義関数が持つ
int fib(int n) {
if (n <= 0) reject("n must be positive");
return n <= 2 ? 1 : fib(n - 1) + fib(n - 2);
}にパースする。したがって、fib 関数を外部の fib.hpp ファイルで定義することには、speed の利点はない。しかし、関数の勾配を解析的に処理するために、Stan の自己微分機能を使用するのではなく、外部の C++ ファイルを使用することは可能である。この場合、fib 関数は整数しか扱わないので、微分を取るものはない。外部 C++ ファイルを使用する第一の利点は、Stan言語では直接できないことを柔軟に行えることである。また、__rstanarm__のような R パッケージが、パッケージのsrcディレクトリにC++関数を定義し、インストール時（またはインストール前）にコンパイルされる Stan プログラムで利用できるようにリンカに依存したい場合にも有用である。
C++版では、n が非正であるかどうかをチェックする。その場合、例外を投げる。MCMCサンプラーの設定により、もし例外が投げられ、それが std::domain_error であれば、その反復は拒絶されるが回復可能なエラーとして扱われる：その反復の対数確率値を負の無限大に設定し、次の反復に移行してみよう。std::invalid_argument が投げられた場合、MCMCは終了する。これらは回復不可能な、ユーザーのプログラミングエラーとして扱いる。 generated C++ ファイルに using std::stringstream; 行があるため、stringstream の前に std:: を付ける必要はない。しかし、対応する using std::domain_error; 行がないので、例外がスローされたときに適切に修飾する必要がある。
C++版の fib 関数の唯一の混乱する部分は、pstream__ という追加の引数（デフォルト値なし）があり、Stan によって fib 関数の declaration に追加されていることである。したがって、fib 関数の definition はこのシグネチャと一致する必要がある。この追加引数は std::ostream へのポインターで、関数が画面に何かをプリントする場合にのみ使用される。したがって、最後の行で fib 関数を再帰的に呼び出すときには、fib(n - 1, 0) + fib(n - 2, 0); を指定して、出力（もしあれば、この場合は何もない）がヌルポイン タに向けられるようにしている。
この vignette では、フィボナッチ関数をおもちゃのように使っているが、これは Stan プログラムではほとんど使い道がなく、もし有用であれば、冒頭で説明したように、functions ブロックのユーザー定義関数としてもっと簡単に実装できるはずである。外部 C++ コードの利用が可能になるのは、より複雑な C++ 関数の場合のみである。言うまでもなく、Stan は推定を行うために未知のパラメータに対する導関数を必要とするので、このメカニズムでは通常C、Fortran、R、その他の言語の関数を呼び出すことができない。この導関数はC++のカスタム型で処理されるため、double , float などのプリミティブ型しか扱えない他言語の関数では処理できない。
とはいえ、特に Stan Math Library を活用すれば、C++ で多くのことを実現することが可能である。詳しくは、 The Stan Math Library: Reverse-Mode Automatic Differentiation in C++ とその GitHub repositoryを参照。Stan プログラムの functions ブロックで_宣言する関数は、C++ にパースされるとき、通常、そのシグネチャにテンプレート化と型推進が含まれる（唯一の例外は、上記の fib 関数のように、引数が整数だけの関数）。例えば、引数が実数である（あるいは、少なくともどちらかの引数が実数である）関数を定義したいとする。たとえば
mc <-
'
functions { real sinc(real x); }
transformed data { real sinc_pi = sinc(pi()); }
'
stan_model(model_code = mc, model_name = "external", allow_undefined = TRUE,
includes = paste0('\n#include "',
file.path(getwd(), 'sinc.hpp'), '"\n'))## Trying to compile a simple C file## Running \
## /Volumes/HDD/opt/sw/Library/Frameworks/R.framework/Versions/4.2/Resources/bin/R \
## CMD SHLIB foo.c
## clang -I"/Volumes/HDD/opt/sw/Library/Frameworks/R.framework/Versions/4.2/Resources/include" -DNDEBUG -I"/Users/baba/Library/R/x86_64/4.2/library/Rcpp/include/" -I"/Users/baba/Library/R/x86_64/4.2/library/RcppEigen/include/" -I"/Users/baba/Library/R/x86_64/4.2/library/RcppEigen/include/unsupported" -I"/Users/baba/Library/R/x86_64/4.2/library/BH/include" -I"/Users/baba/Library/R/x86_64/4.2/library/StanHeaders/include/src/" -I"/Users/baba/Library/R/x86_64/4.2/library/StanHeaders/include/" -I"/Users/baba/Library/R/x86_64/4.2/library/RcppParallel/include/" -I"/Users/baba/Library/R/x86_64/4.2/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -DBOOST_NO_AUTO_PTR -include '/Users/baba/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/Volumes/HDD/opt/sw/include -fPIC -g -O3 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Users/baba/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Users/baba/Library/R/x86_64/4.2/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Users/baba/Library/R/x86_64/4.2/library/RcppEigen/include/Eigen/Core:88:
## /Users/baba/Library/R/x86_64/4.2/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Users/baba/Library/R/x86_64/4.2/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Users/baba/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Users/baba/Library/R/x86_64/4.2/library/RcppEigen/include/Eigen/Dense:1:
## /Users/baba/Library/R/x86_64/4.2/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1## Error in compileCode(f, code, language = language, verbose = verbose): ^/Users/baba/Library/R/x86_64/4.2/library/Rcpp/include/Rcpp/sugar/functions/var.h:89:33: note: candidate found by name lookup is 'Rcpp::var'inline sugar::Var<CPLXSXP,NA,T> var( const VectorBase<CPLXSXP,NA,T>& t){ ^18 warnings and 3 errors generated.make: *** [file126a87f90c5c4.o] Error 1## Error in sink(type = "output"): コネクションが不正ですsinc.hpp ファイルは次のように読む。
double
sinc(const double& x, std::ostream* pstream__) {
return x != 0.0 ? sin(x) / x : 1.0;
}
var
sinc(const var& x, std::ostream* pstream__) {
double x_ = x.val();
double f = x_ != 0.0 ? sin(x_) / x_ : 1.0;
double dfdx_ = x_ != 0.0 ? (cos(x_) - sin(x_)) / x_ : 0.0;
return var(new precomp_v_vari(f, x.vi_, dfdx_));
}最初の関数 sinc の本体は、単に三項演算子という形で数学的に定義されており、入力が double であればそれで十分である。
sinc.hpp の最後の行は、入力が未知の実数パラメータであるときのための特殊化で、 Stan Math Library では stan::math::var として表現されている。sinc 関数の導関数は解析的に計算しやすいので、入力された stan::math::var から基となる倍精度値を取り出し、それを使って関数の値とその1次導関数を計算する。そして、関数値( f )、他のパラメータに対する x の微分( x.vi_ )、そして f の一次微分( dfdx_ )を引数とする precomputed_gradients の結果を返す。後者の2つは実際には std::vector であるが、未知数が1つしかないため、それぞれ1つの要素しか持たない。
生成された関数シグネチャがどうなるかを見る簡単な方法は、rstan パッケージの stanc 関数を allow_undefined = TRUE で呼び出して、再定義された C++ コードを検査することである。この場合、私はまず
try(readLines(stanc(model_code = mc, allow_undefined = TRUE)$cppcode))## Warning in file(con, "r"): expanded path length 13273 would be too long for
## #ifndef USE_STANC3
## #define USE_STANC3
## #endif
## // Code generated by stanc v2.26.1-1-g67504470
## #include <stan/model/model_header.hpp>
## namespace model126a82d1baa9c_file126a82616f4b5_namespace {
## inline void validate_positive_index(const char* var_name, const char* expr,
## int val) {
## if (val < 1) {
## std::stringstream msg;
## msg << "Found dimension size less than one in simplex declaration"
## << "; variable=" << var_name << "; dimension size expression=" << expr
## << "; expression value=" << val;
## std::string msg_str(msg.str());
## throw std::invalid_argument(msg_str.c_str());
## }
## }
## inline void validate_unit_vector_index(const char* var_name, const char* expr,
## int val) {
## if (val <= 1) {
## std::stringstream msg;
## if (val == 1) {
## msg << "Found dimension size one in unit vector declaration."
## << " One-dimensional unit vector is dis [... 省略]
## Warning in file(con, "r"): expanded path length 13273 would be too long for
## #ifndef USE_STANC3
## #define USE_STANC3
## #endif
## // Code generated by stanc v2.26.1-1-g67504470
## #include <stan/model/model_header.hpp>
## namespace model126a82d1baa9c_file126a82616f4b5_namespace {
## inline void validate_positive_index(const char* var_name, const char* expr,
## int val) {
## if (val < 1) {
## std::stringstream msg;
## msg << "Found dimension size less than one in simplex declaration"
## << "; variable=" << var_name << "; dimension size expression=" << expr
## << "; expression value=" << val;
## std::string msg_str(msg.str());
## throw std::invalid_argument(msg_str.c_str());
## }
## }
## inline void validate_unit_vector_index(const char* var_name, const char* expr,
## int val) {
## if (val <= 1) {
## std::stringstream msg;
## if (val == 1) {
## msg << "Found dimension size one in unit vector declaration."
## << " One-dimensional unit vector is dis [... 省略]## Warning in file(con, "r"): ファイル '#ifndef USE_STANC3
## #define USE_STANC3
## #endif
## // Code generated by stanc v2.26.1-1-g67504470
## #include <stan/model/model_header.hpp>
## namespace model126a82d1baa9c_file126a82616f4b5_namespace {
## inline void validate_positive_index(const char* var_name, const char* expr,
## int val) {
## if (val < 1) {
## std::stringstream msg;
## msg << "Found dimension size less than one in simplex declaration"
## << "; variable=" << var_name << "; dimension size expression=" << expr
## << "; expression value=" << val;
## std::string msg_str(msg.str());
## throw std::invalid_argument(msg_str.c_str());
## }
## }
## inline void validate_unit_vector_index(const char* var_name, const char* expr,
## int val) {
## if (val <= 1) {
## std::stringstream msg;
## if (val == 1) {
## msg << "Found dimension size one in unit vector declaration."
## << " One-dimensional unit vector is discrete"
## << " but the targ [... 省略]## Error in file(con, "r") : コネクションを開くことができませんで、sinc.hpp にどのような関数シグネチャを書く必要があるのかを確認した。
数値的に安定な C++ 関数を書いた場合、そのC++関数を Stan Math Library に収録し、誰でも利用できるように GitHub にプルリクエストしていただきたい。
Stan Math Library は C++14 標準に準拠しているが、すべてのコンパイラが C++14 標準を完全にサポートしているわけではない。特に、RToolsに付属するコンパイラは、すべての C++14 のフィーチャをサポートしているわけではない。しかし、auto キーワードのような多くの C++14 のフィーチャを使用することができる。