diff --git a/nCompiler/R/Rcpp_nCompiler_plugin.R b/nCompiler/R/Rcpp_nCompiler_plugin.R index 3fc4f749..7db2a93b 100644 --- a/nCompiler/R/Rcpp_nCompiler_plugin.R +++ b/nCompiler/R/Rcpp_nCompiler_plugin.R @@ -45,7 +45,8 @@ make_nCompiler_plugin <- function(nCompiler_pluginEnv) { else "") result$env$PKG_CXXFLAGS <- "-std=c++17" - result$env$PKG_LIBS <- get_nCompLocal_PKG_LIBS_entry() + result$env$PKG_LIBS <- paste(get_nCompLocal_PKG_LIBS_entry(), + "$(LAPACK_LIBS) $(BLAS_LIBS)") if(isTRUE(nCompiler_pluginEnv$uses_TBB)) result$env <- setEnvTBB(result$env) ## Makevars doesn't work for sourceCpp plugins @@ -73,7 +74,8 @@ make_nCompiler_Eigen_plugin <- function(nCompiler_pluginEnv) { else "") result$env$PKG_CXXFLAGS <- "-std=c++17" - result$env$PKG_LIBS <- get_nCompLocal_PKG_LIBS_entry() + result$env$PKG_LIBS <- paste(get_nCompLocal_PKG_LIBS_entry(), + "$(LAPACK_LIBS) $(BLAS_LIBS)") if(isTRUE(get_nOption('compilerOptions')$throwEigenErrors)) { # replace include directives to enable Eigen errors #preamble = system.file(file.path('include', 'nCompiler', diff --git a/nCompiler/R/compile_aaa_operatorLists.R b/nCompiler/R/compile_aaa_operatorLists.R index ec9e67bb..ac9152ce 100644 --- a/nCompiler/R/compile_aaa_operatorLists.R +++ b/nCompiler/R/compile_aaa_operatorLists.R @@ -1136,6 +1136,111 @@ updateOperatorDef('rt_nonstandard', 'matchDef', val = function(n, df = 1, mu = 0 updateOperatorDef('runif', 'matchDef', val = function(n, min = 0, max = 1) {}) updateOperatorDef('rweibull', 'matchDef', val = function(n, shape, scale = 1) {}) +assignOperatorDef( + c("dmnorm_chol"), + list( + matchDef = function(x, mean, cholesky, prec_param = TRUE, log = FALSE) {}, + labelAbstractTypes = list( + handler = 'MVDist' + ) + ) +) + +assignOperatorDef( + c("rmnorm_chol"), + list( + matchDef = function(n = 1, mean, cholesky, prec_param = TRUE) {}, + labelAbstractTypes = list( + handler = 'MVDist', + return_nDim = 1 + ) + ) +) + +assignOperatorDef( + c("dmvt_chol"), + list( + matchDef = function(x, mu, cholesky, df, prec_param = TRUE, log = FALSE) {}, + labelAbstractTypes = list( + handler = 'MVDist' + ) + ) +) + +assignOperatorDef( + c("rmvt_chol"), + list( + matchDef = function(n = 1, mu, cholesky, df, prec_param = TRUE) {}, + labelAbstractTypes = list( + handler = 'MVDist', + return_nDim = 1 + ) + ) +) + +assignOperatorDef( + c("dwish_chol"), + list( + matchDef = function(x, cholesky, df, scale_param = TRUE, log = FALSE) {}, + labelAbstractTypes = list( + handler = 'MVDist' + ) + ) +) + +assignOperatorDef( + c("rwish_chol"), + list( + matchDef = function(n = 1, cholesky, df, scale_param = TRUE) {}, + labelAbstractTypes = list( + handler = 'MVDist', + return_nDim = 2 + ) + ) +) + +assignOperatorDef( + c("ddirch"), + list( + matchDef = function(x, alpha, log = FALSE) {}, + labelAbstractTypes = list( + handler = 'MVDist' + ) + ) +) + +assignOperatorDef( + c("rdirch"), + list( + matchDef = function(n = 1, alpha) {}, + labelAbstractTypes = list( + handler = 'MVDist', + return_nDim = 1 + ) + ) +) + +assignOperatorDef( + c("dmulti"), + list( + matchDef = function(x, size, prob, log = FALSE) {}, + labelAbstractTypes = list( + handler = 'MVDist' + ) + ) +) + +assignOperatorDef( + c("rmulti"), + list( + matchDef = function(n = 1, size, prob) {}, + labelAbstractTypes = list( + handler = 'MVDist', + return_nDim = 1 + ) + ) +) + assignOperatorDef( c('nSwitch'), list( diff --git a/nCompiler/R/compile_labelAbstractTypes.R b/nCompiler/R/compile_labelAbstractTypes.R index 232010ec..2e6fb106 100644 --- a/nCompiler/R/compile_labelAbstractTypes.R +++ b/nCompiler/R/compile_labelAbstractTypes.R @@ -568,6 +568,20 @@ inLabelAbstractTypesEnv( } ) +inLabelAbstractTypesEnv( + # Used for both multivariate density functions (e.g. dmnorm_chol), which + # return a scalar (nDim = 0, the default via setReturn_nDim), and their + # random-generation counterparts (e.g. rmnorm_chol), which return a + # vector -- set via `labelAbstractTypes = list(handler = 'MVDist', + # return_nDim = 1)` on the latter's operatorDef. + MVDist <- function(code, symTab, auxEnv, handlingInfo) { + inserts <- recurse_labelAbstractTypes(code, symTab, auxEnv, handlingInfo) + code$type <- symbolBasic$new(type = 'double', + nDim = setReturn_nDim(handlingInfo, 0)) + invisible(inserts) + } +) + inLabelAbstractTypesEnv( nChol <- function(code, symTab, auxEnv, handlingInfo) { inserts <- recurse_labelAbstractTypes(code, symTab, auxEnv, handlingInfo) diff --git a/nCompiler/R/cppDefs_nClass.R b/nCompiler/R/cppDefs_nClass.R index 68d542a1..bcffeb7e 100644 --- a/nCompiler/R/cppDefs_nClass.R +++ b/nCompiler/R/cppDefs_nClass.R @@ -12,11 +12,13 @@ nClassBaseClass_init_impl <- function(cppDef) { cppDef$Hpreamble <- c(cppDef$Hpreamble, "#define NCOMPILER_USES_EIGEN", "#define NCOMPILER_USES_NCPPVEC", + "#define NCOMPILER_USES_MVDISTS", "#define USES_NCOMPILER") cppDef$CPPpreamble <- pluginIncludes cppDef$CPPpreamble <- c(cppDef$CPPpreamble, "#define NCOMPILER_USES_EIGEN", "#define NCOMPILER_USES_NCPPVEC", + "#define NCOMPILER_USES_MVDISTS", "#define USES_NCOMPILER") cppDef$Hincludes <- c(cppDef$Hincludes, diff --git a/nCompiler/R/cppDefs_nFunction.R b/nCompiler/R/cppDefs_nFunction.R index df355e2f..0b055dd6 100644 --- a/nCompiler/R/cppDefs_nFunction.R +++ b/nCompiler/R/cppDefs_nFunction.R @@ -9,16 +9,21 @@ cpp_nFunctionClass_init_impl <- function(cppDef) { cppDef$Hpreamble <- c(cppDef$Hpreamble, "#define NCOMPILER_USES_EIGEN", "#define NCOMPILER_USES_NCPPVEC", + "#define NCOMPILER_USES_MVDISTS", "#define USES_NCOMPILER") ## handler nCppVec in labelAbstractTypes does record in auxEnv if an ## explicit call to nCppVec() was uses. That is the beginning of a smarter ## system for determining what #include (via #define) components are ## really needed. But I punted on further extension for now and ## simply tacked NCPPVEC onto the universal set of includes for now. + ## NCOMPILER_USES_MVDISTS (multivariate distributions, e.g. dmnorm_chol) + ## is likewise just tacked on universally for now; a future pass could + ## set it only when generated code actually calls one. cppDef$CPPpreamble <- pluginIncludes cppDef$CPPpreamble <- c(cppDef$CPPpreamble, "#define NCOMPILER_USES_EIGEN", "#define NCOMPILER_USES_NCPPVEC", + "#define NCOMPILER_USES_MVDISTS", "#define USES_NCOMPILER") cppDef$Hincludes <- c(cppDef$Hincludes)#, ## nCompilerIncludeFile("nCompiler_omnibus_first_h.h")) diff --git a/nCompiler/inst/include/nCompiler/dists/ddirch.h b/nCompiler/inst/include/nCompiler/dists/ddirch.h new file mode 100644 index 00000000..77e90cc7 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/ddirch.h @@ -0,0 +1,54 @@ +// The ddirch (Dirichlet) density kernel (raw double* arguments), ported +// from nimble's dists.cpp. +// +// Header-only and self-contained for the same reason as dmnorm_chol.h: the +// package DLL and code nCompile generates and compiles on a user's machine +// are separate shared objects that cannot reliably link against each other, +// so each side compiles its own copy of this source instead. +// +// Unlike the other dist kernels ported so far, ddirch never mutates x or +// alpha (no BLAS/LAPACK, no in-place overwrite), so there's no +// overwrite_inputs parameter here, and the Eigen::Tensor wrapper +// (ddirch_tensor.h) can use the no-copy reference path for both arguments. +// +// See also ddirch_tensor.h for an Eigen::Tensor-argument wrapper around +// this kernel, for use from nCompile-generated code. + +#ifndef _NCOMPILER_DISTS_DDIRCH +#define _NCOMPILER_DISTS_DDIRCH + +#include // DBL_EPSILON +#include +#include "dists_utils.h" + +inline double ddirch(double* x, double* alpha, int K, int give_log) { + double sumAlpha(0.0); + double sumX(0.0); + double dens(0.0); + + if (R_IsNA_ANY(x, K) || R_IsNA_ANY(alpha, K)) + return NA_REAL; + if (R_IsNaN_ANY(x, K) || R_IsNaN_ANY(alpha, K)) + return R_NaN; + + for(int i = 0; i < K; i++) { + if(alpha[i] <= 0.0) ML_ERR_return_NAN; + if(x[i] < 0.0 || x[i] > 1.0) return give_log ? R_NegInf : 0.0; + // Rf_lgammafn (not the "lgammafn" convenience macro) is used here + // deliberately: Rcpp's sugar layer #undefs that macro (Rcpp/sugar/ + // undoRmath.h, to avoid polluting the global namespace with names like + // beta/choose/gamma), so this header needs to work whether or not it's + // included after Rcpp.h has already run. + dens += (alpha[i]-1) * log(x[i]) - Rf_lgammafn(alpha[i]); + sumAlpha += alpha[i]; + sumX += x[i]; + } + if(sumX > 1.0 + 10*DBL_EPSILON || sumX < 1.0 - 10*DBL_EPSILON) { + return give_log ? R_NegInf : 0.0; + } + + dens += Rf_lgammafn(sumAlpha); + return give_log ? dens : exp(dens); +} + +#endif // _NCOMPILER_DISTS_DDIRCH diff --git a/nCompiler/inst/include/nCompiler/dists/ddirch_tensor.h b/nCompiler/inst/include/nCompiler/dists/ddirch_tensor.h new file mode 100644 index 00000000..749b1358 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/ddirch_tensor.h @@ -0,0 +1,36 @@ +// Eigen::Tensor-argument wrapper around the ddirch density kernel (see +// ddirch.h). Templated code like this can only ever be header-only -- the +// compiler must see the definition wherever it's instantiated -- so this +// lives in inst/include regardless of the package-DLL/on-the-fly- +// compilation split that motivates ddirch.h itself. + +#ifndef _NCOMPILER_DISTS_DDIRCH_TENSOR +#define _NCOMPILER_DISTS_DDIRCH_TENSOR + +#include +#include "dists_tensor_utils.h" +#include "ddirch.h" + +// Eigen::Tensor wrapper for ddirch. +// +// Unlike dmnorm_chol/dmvt_chol/dwish_chol, ddirch never mutates x or alpha +// (no BLAS/LAPACK, no in-place overwrite), so both arguments can use +// asDenseTensor's no-copy reference path when they're already concrete +// tensors, and only get materialized when the call site passed a lazy +// expression. +template +double ddirch(const TensorExprX &x, const TensorExprAlpha &alpha, int give_log) { + const auto &xEval = asDenseTensor<1>(x); + const auto &alphaEval = asDenseTensor<1>(alpha); + + int K = static_cast(alphaEval.dimension(0)); + + // xEval/alphaEval may alias the caller's own tensors (when no copy was + // needed above); ddirch's double* signature is non-const only because it + // predates const-correctness conventions here, but it never writes + // through x or alpha, so this cast is safe. + return ddirch(const_cast(xEval.data()), + const_cast(alphaEval.data()), K, give_log); +} + +#endif // _NCOMPILER_DISTS_DDIRCH_TENSOR diff --git a/nCompiler/inst/include/nCompiler/dists/dists_post_Rcpp.h b/nCompiler/inst/include/nCompiler/dists/dists_post_Rcpp.h new file mode 100644 index 00000000..d2078828 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/dists_post_Rcpp.h @@ -0,0 +1,33 @@ +// Umbrella header for nCompiler's distribution kernels, for use from +// nCompile-generated code. Pulled into nCompiler_omnibus.h under +// #ifdef NCOMPILER_USES_MVDISTS. +// +// Individual dist headers are self-contained (each brings its own Eigen +// includes, etc.), so this file's only job is to gather them in one place. +// As more dist kernels are added, they just get an #include line here. + +#ifndef _NCOMPILER_DISTS_POST_RCPP +#define _NCOMPILER_DISTS_POST_RCPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#endif // _NCOMPILER_DISTS_POST_RCPP diff --git a/nCompiler/inst/include/nCompiler/dists/dists_tensor_utils.h b/nCompiler/inst/include/nCompiler/dists/dists_tensor_utils.h new file mode 100644 index 00000000..965f3fe0 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/dists_tensor_utils.h @@ -0,0 +1,62 @@ +// Shared Eigen::Tensor materialization helpers for nCompiler's dist +// wrappers (e.g. dmnorm_chol_tensor.h, dmvt_chol_tensor.h). Split out into +// its own header, rather than being duplicated in each *_tensor.h, because +// multiple such wrappers get included together (via dists_post_Rcpp.h) and +// would otherwise redefine the same templates in one translation unit. + +#ifndef _NCOMPILER_DISTS_TENSOR_UTILS +#define _NCOMPILER_DISTS_TENSOR_UTILS + +#include +#include + +// True only when TensorExpr is exactly Eigen::Tensor, i.e. an +// already-materialized, densely-packed tensor rather than a lazy expression +// (a block, a cast, an arithmetic chain, a TensorMap, etc). +template +struct is_dense_double_tensor + : std::is_same> {}; + +// Returns x itself (no copy) when it is already a concrete +// Eigen::Tensor; otherwise assigns it into a fresh +// Eigen::Tensor, which forces Eigen to evaluate the expression +// into a contiguous buffer we can call .data() on. The .cast() is +// required here, not just a safety net: Eigen::Tensor's assignment/ +// construction from another tensor expression requires matching Scalar +// types, so a non-double input (e.g. Eigen::Tensor, or a bool +// expression from a comparison op) would otherwise fail to compile. +template +decltype(auto) asDenseTensor(const TensorExpr &x) { + if constexpr (is_dense_double_tensor::value) { + return (x); + } else { + return Eigen::Tensor(x.template cast()); + } +} + +// Applies nimble's/R's recycling rule for a "mean"-like argument that may +// be shorter than the target length n: if v has fewer than n elements, its +// elements are reused cyclically to fill out length n (matching the +// full_mean/full_mu construction in dists.cpp's C_dmnorm_chol, +// C_rmnorm_chol, C_dmvt_chol, and C_rmvt_chol). If v already has at least n +// elements, it's returned unchanged (extra elements beyond n are simply +// never read by the kernels that consume it, same as the original code). +// +// Unlike asDenseTensor, this always returns by value: whether recycling is +// needed is a runtime property (v's length vs. n), not something knowable +// from the type system, so there's no way to give the "no recycling +// needed" case a reference-typed fast path here without splitting the +// return type per branch. +inline Eigen::Tensor recycleToLength(const Eigen::Tensor &v, int n) { + int m = static_cast(v.dimension(0)); + if (m >= n) return v; + Eigen::Tensor out(n); + int j = 0; + for (int i = 0; i < n; i++) { + out(i) = v(j++); + if (j == m) j = 0; + } + return out; +} + +#endif // _NCOMPILER_DISTS_TENSOR_UTILS diff --git a/nCompiler/inst/include/nCompiler/dists/dists_utils.h b/nCompiler/inst/include/nCompiler/dists/dists_utils.h new file mode 100644 index 00000000..2a0d4e11 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/dists_utils.h @@ -0,0 +1,113 @@ +// Small helpers used by nCompiler's distribution kernels (ported from +// nimble's dists.cpp). Kept header-only and inline, rather than living only +// in src/dists.cpp, so the same implementation can be compiled both into +// the package DLL and into code nCompile generates and compiles on a +// user's machine -- see dmnorm_chol.h for why that split is necessary. + +#ifndef _NCOMPILER_DISTS_UTILS +#define _NCOMPILER_DISTS_UTILS + +// R_IsNA, R_IsNaN, ISNAN, and R_FINITE are declared in R_ext/Arith.h, not +// Rmath.h; NA_REAL, R_NaN, and R_NegInf (used by dmnorm_chol.h) live there +// too. +#include +#include // Rf_warning, used by MATHLIB_WARNING below +#include + +// Detects NA +inline bool R_IsNA_ANY(double* P, int s) { + for(int i = 0; i < s; ++i) if(R_IsNA(P[i])) return(true); + return(false); +} + +// Detects NaN +inline bool R_IsNaN_ANY(double* P, int s) { + for(int i = 0; i < s; ++i) if(R_IsNaN(P[i])) return(true); + return(false); +} + +// Detects NA or NaN +inline bool ISNAN_ANY(double* P, int s) { + for(int i = 0; i < s; ++i) if(ISNAN(P[i])) return(true); + return(false); +} + +inline bool R_FINITE_ANY(double* P, int s) { + for(int i = 0; i < s; ++i) if(!R_FINITE(P[i])) return(false); + return(true); +} + +// Mirrors src/Utils.h's ML_ERROR/ML_ERR_return_NAN, guarded so that if +// Utils.h has already been included in this translation unit (as it is via +// dists.h when building the package DLL), its definitions win and these are +// skipped. This is currently a hollow shell -- the ME_DOMAIN case (the only +// one dist kernels use today) never reaches MATHLIB_WARNING -- but it's +// kept in this form, not collapsed to "return ML_NAN;", so that filling in +// real warnings for other error classes later only requires editing here +// (and in Utils.h), not re-auditing every dist kernel that calls it. +#ifndef MATHLIB_WARNING +#define MATHLIB_WARNING(fmt,x) Rf_warning(fmt,x) +#endif + +// These are originally defined in Utils.h +// For on-the-fly compilation these will be read? + +#ifndef ME_NONE +#define ME_NONE 0 +#define ME_DOMAIN 1 +#define ME_RANGE 2 +#define ME_NOCONV 4 +#define ME_PRECISION 8 +#define ME_UNDERFLOW 16 +#endif + +#ifndef ML_ERROR +#define ML_ERROR(x, s) { \ + if(x > ME_DOMAIN) { \ + const char *msg = ""; \ + switch(x) { \ + case ME_DOMAIN: \ + msg = "argument out of domain in '%s'\n"; \ + break; \ + case ME_RANGE: \ + msg = "value out of range in '%s'\n"; \ + break; \ + case ME_NOCONV: \ + msg = "convergence failed in '%s'\n"; \ + break; \ + case ME_PRECISION: \ + msg = "full precision may not have been achieved in '%s'\n"; \ + break; \ + case ME_UNDERFLOW: \ + msg = "underflow occurred in '%s'\n"; \ + break; \ + } \ + MATHLIB_WARNING(msg, s); \ + } \ +} +#endif + +#ifndef ML_ERR_return_NAN +#define ML_ERR_return_NAN { ML_ERROR(ME_DOMAIN, ""); return R_NaN; } +#endif + +// Mirrors src/Utils.h's R_D_forceint/R_D_nonint/R_D_negInonint (from the +// "dpq.h" section), as plain inline functions rather than macros since +// none of the three reference give_log/log_p. Guarded the same way as the +// ML_ERROR family above, so Utils.h's macro versions win when this header +// is included alongside it (they'd otherwise clash: Utils.h defines these +// as object-like function macros, which would mangle these very +// definitions via textual substitution if both were active). +#ifndef R_D_forceint +inline double R_D_forceint(double x) { return floor(x + 0.5); } +#endif + +#ifndef R_D_nonint +inline bool R_D_nonint(double x) { return fabs(x - floor(x + 0.5)) > 1e-7; } +#endif + +#ifndef R_D_negInonint +inline bool R_D_negInonint(double x) { return x < 0. || R_D_nonint(x); } +#endif + +#endif // _NCOMPILER_DISTS_UTILS diff --git a/nCompiler/inst/include/nCompiler/dists/dmnorm_chol.h b/nCompiler/inst/include/nCompiler/dists/dmnorm_chol.h new file mode 100644 index 00000000..5ec4b659 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/dmnorm_chol.h @@ -0,0 +1,94 @@ +// The dmnorm_chol density kernel (raw double* arguments), ported from +// nimble's dists.cpp. +// +// This is deliberately header-only and self-contained (rather than living +// only in src/dists.cpp) so that the identical implementation can be +// compiled both into the package DLL (for the R-facing C_dmnorm_chol / +// .Call path) and into code nCompile generates and compiles on a user's +// machine. Those are two separate shared objects, and a package DLL (often +// built by CRAN) cannot be relied on to link against code compiled locally +// on a user's computer, across all supported OSes/compilers. Rather than +// linking across that boundary, each side just compiles its own copy of +// this same source, the way header-only libraries like Eigen already do. +// +// See also dmnorm_chol_tensor.h for an Eigen::Tensor-argument wrapper +// around this kernel, for use from nCompile-generated code. + +#ifndef _NCOMPILER_DISTS_DMNORM_CHOL +#define _NCOMPILER_DISTS_DMNORM_CHOL + +// Based on Writing R Extensions 6.6.1 and mgcv's mgcv.h +#define USE_FC_LEN_T +#include +#include +/* If we are compiling with a version of R before FCONE and the explicit + supplying of extra arguments was introduced, then FCONE has to be + defined */ +#ifndef FCONE +#define FCONE +#endif + +#include +#include "dists_utils.h" + +inline double dmnorm_chol(double* x, double* mean, double* chol, int n, + double prec_param, int give_log, + int overwrite_inputs) { + char uplo('U'); + char transPrec('N'); + char transCov('T'); + char diag('N'); + int lda(n); + int incx(1); + double* xCopy; + + double dens = -n * M_LN_SQRT_2PI; + int i; + // add diagonals of Cholesky + + if (R_IsNA_ANY(x, n) || R_IsNA_ANY(mean, n) || R_IsNA_ANY(chol, n*n) || R_IsNA(prec_param)) + return NA_REAL; + if (R_IsNaN_ANY(x, n) || R_IsNaN_ANY(mean, n) || R_IsNaN_ANY(chol, n*n) || R_IsNaN(prec_param)) + return R_NaN; + + if(!R_FINITE_ANY(x, n) || !R_FINITE_ANY(mean, n) || !R_FINITE_ANY(chol, n*n)) + return give_log ? R_NegInf : 0.0; + + if(prec_param) { + for(i = 0; i < n*n; i += n + 1) + dens += log(chol[i]); + } else { + for(i = 0; i < n*n; i += n + 1) + dens -= log(chol[i]); + } + + if(overwrite_inputs) { + xCopy = x; + for(i = 0; i < n; i++) + xCopy[i] -= mean[i]; + } else { + xCopy = new double[n]; + for(i = 0; i < n; i++) + xCopy[i] = x[i] - mean[i]; + } + + // do matrix-vector multiply with upper-triangular matrix stored column-wise as full n x n matrix (prec parameterization) + // or upper-triangular (transpose) solve (cov parameterization) + // dtr{m,s}v is a BLAS level-2 function + if(prec_param) F77_CALL(dtrmv)(&uplo, &transPrec, &diag, &n, chol, &lda, xCopy, &incx FCONE FCONE FCONE); + else F77_CALL(dtrsv)(&uplo, &transCov, &diag, &n, chol, &lda, xCopy, &incx FCONE FCONE FCONE); + + // sum of squares to calculate quadratic form + double tmp = 0.0; + for(i = 0; i < n; i++) + tmp += xCopy[i] * xCopy[i]; + + dens += -0.5 * tmp; + + if(!overwrite_inputs) + delete [] xCopy; + + return give_log ? dens : exp(dens); +} + +#endif // _NCOMPILER_DISTS_DMNORM_CHOL diff --git a/nCompiler/inst/include/nCompiler/dists/dmnorm_chol_tensor.h b/nCompiler/inst/include/nCompiler/dists/dmnorm_chol_tensor.h new file mode 100644 index 00000000..5ba40141 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/dmnorm_chol_tensor.h @@ -0,0 +1,50 @@ +// Eigen::Tensor-argument wrapper around the dmnorm_chol density kernel +// (see dmnorm_chol.h). Templated code like this can only ever be header- +// only -- the compiler must see the definition wherever it's instantiated +// -- so this lives in inst/include regardless of the package-DLL/ +// on-the-fly-compilation split that motivates dmnorm_chol.h itself. + +#ifndef _NCOMPILER_DISTS_DMNORM_CHOL_TENSOR +#define _NCOMPILER_DISTS_DMNORM_CHOL_TENSOR + +#include +#include "dists_tensor_utils.h" +#include "dmnorm_chol.h" + +// Eigen::Tensor wrapper for dmnorm_chol. +// +// x, mean, and chol are accepted as arbitrary Eigen::Tensor expressions +// built up at the call site. chol is read-only, so asDenseTensor reuses it +// directly (no copy) when it's already a concrete tensor, and only +// materializes a copy when the call site passed a lazy expression. x is +// always copied into its own buffer: dmnorm_chol overwrites its working copy +// of x in place, and we must not mutate the caller's tensor as a side +// effect. mean follows nimble's recycling rule (see recycleToLength in +// dists_tensor_utils.h): if it's shorter than x, its elements are reused +// cyclically to fill out x's length, matching the original C_dmnorm_chol's +// full_mean construction -- so mean is always materialized to x's length +// here rather than taking asDenseTensor's no-copy path. +template +double dmnorm_chol(const TensorExprX &x, const TensorExprMean &mean, + const TensorExprChol &chol, double prec_param, int give_log) { + const auto &cholEval = asDenseTensor<2>(chol); + // x is materialized unconditionally (never via asDenseTensor's reference + // path), so cast() is needed here too in case TensorExprX has a + // non-double Scalar. + Eigen::Tensor xEval = x.template cast(); + + int n = static_cast(xEval.dimension(0)); + + Eigen::Tensor meanEval = + recycleToLength(asDenseTensor<1>(mean), n); + + // cholEval may alias the caller's own tensor (when no copy was needed + // above); dmnorm_chol's double* signature is non-const only because it + // predates const-correctness conventions here, but it never writes + // through mean or chol, so this cast is safe. + return dmnorm_chol(xEval.data(), meanEval.data(), + const_cast(cholEval.data()), n, + prec_param, give_log, /*overwrite_inputs=*/1); +} + +#endif // _NCOMPILER_DISTS_DMNORM_CHOL_TENSOR diff --git a/nCompiler/inst/include/nCompiler/dists/dmulti.h b/nCompiler/inst/include/nCompiler/dists/dmulti.h new file mode 100644 index 00000000..e7cec198 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/dmulti.h @@ -0,0 +1,72 @@ +// The dmulti (multinomial) density kernel (raw double* arguments), ported +// from nimble's dists.cpp. +// +// Header-only and self-contained for the same reason as dmnorm_chol.h: the +// package DLL and code nCompile generates and compiles on a user's machine +// are separate shared objects that cannot reliably link against each other, +// so each side compiles its own copy of this source instead. +// +// dmulti always mutates x in place (rounding each entry to the nearest +// integer via R_D_forceint), unconditionally -- there's no +// overwrite_inputs parameter here, unlike dmnorm_chol/dmvt_chol/dwish_chol +// -- so the Eigen::Tensor wrapper (dmulti_tensor.h) must always materialize +// its own private copy of x. +// +// See also dmulti_tensor.h for an Eigen::Tensor-argument wrapper around +// this kernel, for use from nCompile-generated code. + +#ifndef _NCOMPILER_DISTS_DMULTI +#define _NCOMPILER_DISTS_DMULTI + +#include // DBL_EPSILON +#include +#include "dists_utils.h" + +// Calling functions need to copy first arg to int if needed. +inline double dmulti(double* x, double size, double* prob, int K, int give_log) { + double sumProb(0.0); + double sumX(0.0); + double logSumProb; + + if (R_IsNA_ANY(x, K) || R_IsNA_ANY(prob, K) || R_IsNA(size)) + return NA_REAL; + if (R_IsNaN_ANY(x, K) || R_IsNaN_ANY(prob, K) || R_IsNaN(size)) + return R_NaN; + + if(R_D_negInonint(size)) + ML_ERR_return_NAN; + size = R_D_forceint(size); + + // Rf_lgammafn (not the "lgammafn" convenience macro) is used here + // deliberately: Rcpp's sugar layer #undefs that macro (Rcpp/sugar/ + // undoRmath.h, to avoid polluting the global namespace with names like + // beta/choose/gamma), so this header needs to work whether or not it's + // included after Rcpp.h has already run. + double dens = Rf_lgammafn(size + 1); + for(int i = 0; i < K; i++) { + if (prob[i] < 0) ML_ERR_return_NAN; + if (R_D_nonint(x[i])) { + MATHLIB_WARNING("non-integer x = %f", x[i]); + return give_log ? R_NegInf : 0.0; + } + if (x[i] < 0 || !R_FINITE(x[i])) return give_log ? R_NegInf : 0.0; + + x[i] = R_D_forceint(x[i]); + sumProb += prob[i]; + sumX += x[i]; + } + logSumProb = log(sumProb); + + for(int i = 0; i < K; i++) { + if(!(x[i] == 0.0 && prob[i] == 0.0)) + dens += x[i]*(log(prob[i]) - logSumProb) - Rf_lgammafn(x[i] + 1); + } + + if(sumX > size + 10*DBL_EPSILON || sumX < size - 10*DBL_EPSILON) { + return give_log ? R_NegInf : 0.0; + } + + return give_log ? dens : exp(dens); +} + +#endif // _NCOMPILER_DISTS_DMULTI diff --git a/nCompiler/inst/include/nCompiler/dists/dmulti_tensor.h b/nCompiler/inst/include/nCompiler/dists/dmulti_tensor.h new file mode 100644 index 00000000..1a467c5a --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/dmulti_tensor.h @@ -0,0 +1,38 @@ +// Eigen::Tensor-argument wrapper around the dmulti density kernel (see +// dmulti.h). Templated code like this can only ever be header-only -- the +// compiler must see the definition wherever it's instantiated -- so this +// lives in inst/include regardless of the package-DLL/on-the-fly- +// compilation split that motivates dmulti.h itself. + +#ifndef _NCOMPILER_DISTS_DMULTI_TENSOR +#define _NCOMPILER_DISTS_DMULTI_TENSOR + +#include +#include "dists_tensor_utils.h" +#include "dmulti.h" + +// Eigen::Tensor wrapper for dmulti. +// +// prob is read-only, so asDenseTensor reuses it directly (no copy) when +// it's already a concrete tensor, and only materializes a copy when the +// call site passed a lazy expression. x is always copied into its own +// buffer: dmulti unconditionally rounds each entry of its working copy of +// x to the nearest integer in place, and we must not mutate the caller's +// tensor as a side effect. +template +double dmulti(const TensorExprX &x, double size, const TensorExprProb &prob, + int give_log) { + const auto &probEval = asDenseTensor<1>(prob); + Eigen::Tensor xEval = x.template cast(); + + int K = static_cast(probEval.dimension(0)); + + // probEval may alias the caller's own tensor (when no copy was needed + // above); dmulti's double* signature is non-const only because it + // predates const-correctness conventions here, but it never writes + // through prob, so this cast is safe. + return dmulti(xEval.data(), size, const_cast(probEval.data()), K, + give_log); +} + +#endif // _NCOMPILER_DISTS_DMULTI_TENSOR diff --git a/nCompiler/inst/include/nCompiler/dists/dmvt_chol.h b/nCompiler/inst/include/nCompiler/dists/dmvt_chol.h new file mode 100644 index 00000000..8b856688 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/dmvt_chol.h @@ -0,0 +1,94 @@ +// The dmvt_chol (multivariate t) density kernel (raw double* arguments), +// ported from nimble's dists.cpp. +// +// Header-only and self-contained for the same reason as dmnorm_chol.h: the +// package DLL and code nCompile generates and compiles on a user's machine +// are separate shared objects that cannot reliably link against each other, +// so each side compiles its own copy of this source instead. +// +// See also dmvt_chol_tensor.h for an Eigen::Tensor-argument wrapper around +// this kernel, for use from nCompile-generated code. + +#ifndef _NCOMPILER_DISTS_DMVT_CHOL +#define _NCOMPILER_DISTS_DMVT_CHOL + +// Based on Writing R Extensions 6.6.1 and mgcv's mgcv.h +#define USE_FC_LEN_T +#include +#include +/* If we are compiling with a version of R before FCONE and the explicit + supplying of extra arguments was introduced, then FCONE has to be + defined */ +#ifndef FCONE +#define FCONE +#endif + +#include +#include "dists_utils.h" + +inline double dmvt_chol(double* x, double* mu, double* chol, double df, + int n, double prec_param, int give_log, + int overwrite_inputs) { + char uplo('U'); + char transPrec('N'); + char transCov('T'); + char diag('N'); + int lda(n); + int incx(1); + double* xCopy; + + // Rf_lgammafn (not the "lgammafn" convenience macro) is used here + // deliberately: Rcpp's sugar layer #undefs that macro (Rcpp/sugar/ + // undoRmath.h, to avoid polluting the global namespace with names like + // beta/choose/gamma), so this header needs to work whether or not it's + // included after Rcpp.h has already run. + double dens = Rf_lgammafn((df + n) / 2) - Rf_lgammafn(df / 2) - n * M_LN_SQRT_PI - n * log(df) / 2; + int i; + + if (R_IsNA_ANY(x, n) || R_IsNA_ANY(mu, n) || R_IsNA_ANY(chol, n*n) || R_IsNA(df) || R_IsNA(prec_param)) + return NA_REAL; + if (R_IsNaN_ANY(x, n) || R_IsNaN_ANY(mu, n) || R_IsNaN_ANY(chol, n*n) || R_IsNA(df) || R_IsNaN(prec_param)) + return R_NaN; + + if(!R_FINITE_ANY(x, n) || !R_FINITE_ANY(mu, n) || !R_FINITE_ANY(chol, n*n)) + return give_log ? R_NegInf : 0.0; + + // add diagonals of Cholesky + if(prec_param) { + for(i = 0; i < n*n; i += n + 1) + dens += log(chol[i]); + } else { + for(i = 0; i < n*n; i += n + 1) + dens -= log(chol[i]); + } + + if(overwrite_inputs) { + xCopy = x; + for(i = 0; i < n; i++) + xCopy[i] -= mu[i]; + } else { + xCopy = new double[n]; + for(i = 0; i < n; i++) + xCopy[i] = x[i] - mu[i]; + } + + // do matrix-vector multiply with upper-triangular matrix stored column-wise as full n x n matrix (prec parameterization) + // or upper-triangular (transpose) solve (cov parameterization) + // dtr{m,s}v is a BLAS level-2 function + if(prec_param) F77_CALL(dtrmv)(&uplo, &transPrec, &diag, &n, chol, &lda, xCopy, &incx FCONE FCONE FCONE); + else F77_CALL(dtrsv)(&uplo, &transCov, &diag, &n, chol, &lda, xCopy, &incx FCONE FCONE FCONE); + + // sum of squares to calculate quadratic form + double tmp = 0.0; + for(i = 0; i < n; i++) + tmp += xCopy[i] * xCopy[i]; + + dens += -0.5 * (df + n) * log(1 + tmp / df); + + if(!overwrite_inputs) + delete [] xCopy; + + return give_log ? dens : exp(dens); +} + +#endif // _NCOMPILER_DISTS_DMVT_CHOL diff --git a/nCompiler/inst/include/nCompiler/dists/dmvt_chol_tensor.h b/nCompiler/inst/include/nCompiler/dists/dmvt_chol_tensor.h new file mode 100644 index 00000000..54518601 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/dmvt_chol_tensor.h @@ -0,0 +1,51 @@ +// Eigen::Tensor-argument wrapper around the dmvt_chol density kernel +// (see dmvt_chol.h). Templated code like this can only ever be header- +// only -- the compiler must see the definition wherever it's instantiated +// -- so this lives in inst/include regardless of the package-DLL/ +// on-the-fly-compilation split that motivates dmvt_chol.h itself. + +#ifndef _NCOMPILER_DISTS_DMVT_CHOL_TENSOR +#define _NCOMPILER_DISTS_DMVT_CHOL_TENSOR + +#include +#include "dists_tensor_utils.h" +#include "dmvt_chol.h" + +// Eigen::Tensor wrapper for dmvt_chol. +// +// x, mu, and chol are accepted as arbitrary Eigen::Tensor expressions built +// up at the call site. chol is read-only, so asDenseTensor reuses it +// directly (no copy) when it's already a concrete tensor, and only +// materializes a copy when the call site passed a lazy expression. x is +// always copied into its own buffer: dmvt_chol overwrites its working copy +// of x in place, and we must not mutate the caller's tensor as a side +// effect. mu follows nimble's recycling rule (see recycleToLength in +// dists_tensor_utils.h): if it's shorter than x, its elements are reused +// cyclically to fill out x's length, matching the original C_dmvt_chol's +// full_mu construction -- so mu is always materialized to x's length here +// rather than taking asDenseTensor's no-copy path. +template +double dmvt_chol(const TensorExprX &x, const TensorExprMu &mu, + const TensorExprChol &chol, double df, double prec_param, + int give_log) { + const auto &cholEval = asDenseTensor<2>(chol); + // x is materialized unconditionally (never via asDenseTensor's reference + // path), so cast() is needed here too in case TensorExprX has a + // non-double Scalar. + Eigen::Tensor xEval = x.template cast(); + + int n = static_cast(xEval.dimension(0)); + + Eigen::Tensor muEval = + recycleToLength(asDenseTensor<1>(mu), n); + + // cholEval may alias the caller's own tensor (when no copy was needed + // above); dmvt_chol's double* signature is non-const only because it + // predates const-correctness conventions here, but it never writes + // through mu or chol, so this cast is safe. + return dmvt_chol(xEval.data(), muEval.data(), + const_cast(cholEval.data()), df, n, + prec_param, give_log, /*overwrite_inputs=*/1); +} + +#endif // _NCOMPILER_DISTS_DMVT_CHOL_TENSOR diff --git a/nCompiler/inst/include/nCompiler/dists/dwish_chol.h b/nCompiler/inst/include/nCompiler/dists/dwish_chol.h new file mode 100644 index 00000000..4f01cd22 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/dwish_chol.h @@ -0,0 +1,151 @@ +// The dwish_chol (Wishart) density kernel (raw double* arguments), ported +// from nimble's dists.cpp. Unlike dmnorm_chol/dmvt_chol, the random +// variable x here is itself a p x p matrix, not a vector. +// +// Header-only and self-contained for the same reason as dmnorm_chol.h: the +// package DLL and code nCompile generates and compiles on a user's machine +// are separate shared objects that cannot reliably link against each other, +// so each side compiles its own copy of this source instead. +// +// See also dwish_chol_tensor.h for an Eigen::Tensor-argument wrapper around +// this kernel, for use from nCompile-generated code. + +#ifndef _NCOMPILER_DISTS_DWISH_CHOL +#define _NCOMPILER_DISTS_DWISH_CHOL + +// Based on Writing R Extensions 6.6.1 and mgcv's mgcv.h +#define USE_FC_LEN_T +#include +#include +#include +/* If we are compiling with a version of R before FCONE and the explicit + supplying of extra arguments was introduced, then FCONE has to be + defined */ +#ifndef FCONE +#define FCONE +#endif + +#include +#include "dists_utils.h" + +inline double dwish_chol(double* x, double* chol, double df, int p, + double scale_param, int give_log, + int overwrite_inputs) { + char uplo('U'); + char sideL('L'); + char sideR('R'); + char diag('N'); + char transN('N'); + int info(0); + double alpha(1.0); + double* xChol; + + int i, j; + + if (R_IsNA_ANY(x, p*p) || R_IsNA_ANY(chol, p*p) || R_IsNA(df) || R_IsNA(scale_param)) + return NA_REAL; + if (R_IsNaN_ANY(x, p*p) || R_IsNaN_ANY(chol, p*p) || R_IsNaN(df) || R_IsNaN(scale_param)) + return R_NaN; + + // also covers df < 0 + if(df < (double) p) ML_ERR_return_NAN; + + if(!R_FINITE_ANY(x, p*p) || !R_FINITE_ANY(chol, p*p)) + return give_log ? R_NegInf : 0.0; + + // Rf_lgammafn (not the "lgammafn" convenience macro) is used here + // deliberately: Rcpp's sugar layer #undefs that macro (Rcpp/sugar/ + // undoRmath.h, to avoid polluting the global namespace with names like + // beta/choose/gamma), so this header needs to work whether or not it's + // included after Rcpp.h has already run. + double dens = -(df*p/2 * M_LN2 + p*(p-1)*M_LN_SQRT_PI/2); + for(i = 0; i < p; i++) + dens -= Rf_lgammafn((df - i) / 2); + + if(scale_param) { + for(i = 0; i < p*p; i += p + 1) + dens -= df * log(chol[i]); + } else { + for(i = 0; i < p*p; i += p + 1) + dens += df * log(chol[i]); + } + + // determinant of x using Cholesky: + if(overwrite_inputs && (int) scale_param) // if !scale_param we need x below + xChol = x; + else { + xChol = new double[p*p]; + // only need upper triangle for dpotrf chol calculation + for(j = 0; j < p; j++) + for(i = 0; i <= j; i++) + xChol[j*p+i] = x[j*p+i]; + } + F77_CALL(dpotrf)(&uplo, &p, xChol, &p, &info FCONE); + for(i = 0; i < p*p; i += p + 1) + dens += (df - p - 1) * log(xChol[i]); + + // R %*% x = t(chol) %*% chol %*% x (could also do with chol(x) but no more efficient + // solve(S, x) = crossproduct( chol(x) %*% inverse(chol) ) + + // dtr{m,s}m is a BLAS level-3 function + double tmp_dens = 0.0; + if(scale_param) { + // chol(x) %*% inverse(chol) + // need lower triangle of xChol to have zeros as dtrsm assumes it is full matrix + for(j = 0; j < p-1; j++) + for(i = j+1; i < p; i++) + xChol[j*p+i] = 0.0; + F77_CALL(dtrsm)(&sideR, &uplo, &transN, &diag, &p, &p, &alpha, + chol, &p, xChol, &p FCONE FCONE FCONE FCONE); + // trace of crossproduct of result is sum of squares of elements + for(j = 0; j < p; j++) + for(i = 0; i <= j; i++) + tmp_dens += xChol[j*p+i]*xChol[j*p+i]; + } else { + double* xCopy; + if(overwrite_inputs) + xCopy = x; + else { + xCopy = new double[p*p]; + for(i = 0; i < p*p; i++) + xCopy[i] = x[i]; + } + // chol %*% x + F77_CALL(dtrmm)(&sideL, &uplo, &transN, &diag, &p, &p, &alpha, + chol, &p, xCopy, &p FCONE FCONE FCONE FCONE); + // trace crossproduct of t(chol) with result is sum of product of upper-triangular elements + for(j = 0; j < p; j++) { + for(i = 0; i <= j; i++) { + tmp_dens += xCopy[j*p+i] * chol[j*p+i]; + } + } + if(!overwrite_inputs) + delete [] xCopy; + } + + if(!(overwrite_inputs && (int) scale_param)) + delete [] xChol; + + // attempt to improve above calcs by doing efficient U^T U multiply followed by direct product multiply, however this would not make use of threading provided by BLAS and even with one thread seems to be no faster + // U^T*U directly followed by direct product with x + /* + double tmp_summand; + int minij; + for(j = 0; j < p; j++) + for(i = 0; i < p; i++) { + tmp_summand = 0.0; + minij = i <= j ? i : j; + for(int k = 0; k < minij; k++) + tmp_summand += chol[j*p+k]*chol[i*p+k]; // U^T U + // double if not on diagonal to account for direct product of lower triangle too + if(i != j) tmp_summand *= 2; + tmp_dens += xCopy[j*p+i] * tmp_summand; + } + */ + + dens += -0.5 * tmp_dens; + + return give_log ? dens : exp(dens); +} + +#endif // _NCOMPILER_DISTS_DWISH_CHOL diff --git a/nCompiler/inst/include/nCompiler/dists/dwish_chol_tensor.h b/nCompiler/inst/include/nCompiler/dists/dwish_chol_tensor.h new file mode 100644 index 00000000..24dac713 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/dwish_chol_tensor.h @@ -0,0 +1,45 @@ +// Eigen::Tensor-argument wrapper around the dwish_chol density kernel +// (see dwish_chol.h). Templated code like this can only ever be header- +// only -- the compiler must see the definition wherever it's instantiated +// -- so this lives in inst/include regardless of the package-DLL/ +// on-the-fly-compilation split that motivates dwish_chol.h itself. + +#ifndef _NCOMPILER_DISTS_DWISH_CHOL_TENSOR +#define _NCOMPILER_DISTS_DWISH_CHOL_TENSOR + +#include +#include "dists_tensor_utils.h" +#include "dwish_chol.h" + +// Eigen::Tensor wrapper for dwish_chol. +// +// Unlike dmnorm_chol/dmvt_chol, the random variable x here is itself a +// p x p matrix (rank-2 tensor), not a vector. +// +// x and chol are accepted as arbitrary Eigen::Tensor expressions built up +// at the call site. chol is read-only, so asDenseTensor reuses it directly +// (no copy) when it's already a concrete tensor, and only materializes a +// copy when the call site passed a lazy expression. x is always copied into +// its own buffer: dwish_chol overwrites its working copy of x in place +// (regardless of scale_param), and we must not mutate the caller's tensor +// as a side effect. +template +double dwish_chol(const TensorExprX &x, const TensorExprChol &chol, + double df, double scale_param, int give_log) { + const auto &cholEval = asDenseTensor<2>(chol); + // x is materialized unconditionally (never via asDenseTensor's reference + // path), so cast() is needed here too in case TensorExprX has a + // non-double Scalar. + Eigen::Tensor xEval = x.template cast(); + + int p = static_cast(xEval.dimension(0)); + + // cholEval may alias the caller's own tensor (when no copy was needed + // above); dwish_chol's double* signature is non-const only because it + // predates const-correctness conventions here, but it never writes + // through chol, so this cast is safe. + return dwish_chol(xEval.data(), const_cast(cholEval.data()), df, p, + scale_param, give_log, /*overwrite_inputs=*/1); +} + +#endif // _NCOMPILER_DISTS_DWISH_CHOL_TENSOR diff --git a/nCompiler/inst/include/nCompiler/dists/rdirch.h b/nCompiler/inst/include/nCompiler/dists/rdirch.h new file mode 100644 index 00000000..bdead611 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/rdirch.h @@ -0,0 +1,50 @@ +// The rdirch (Dirichlet random generation) kernel (raw double* arguments), +// ported from nimble's dists.cpp. +// +// Header-only and self-contained for the same reason as dmnorm_chol.h: the +// package DLL and code nCompile generates and compiles on a user's machine +// are separate shared objects that cannot reliably link against each other, +// so each side compiles its own copy of this source instead. +// +// Like ddirch, alpha is never mutated (no overwrite_inputs concept), and +// ans is always a fresh output buffer. +// +// See also rdirch_tensor.h for an Eigen::Tensor-argument wrapper around +// this kernel, for use from nCompile-generated code. + +#ifndef _NCOMPILER_DISTS_RDIRCH +#define _NCOMPILER_DISTS_RDIRCH + +#include // rgamma is declared here (via Rmath.h too) +#include +#include "dists_utils.h" + +inline void rdirch(double* ans, double* alpha, int K) { + int i, j; + + if (ISNAN_ANY(alpha, K)) { + for(j = 0; j < K; j++) + ans[j] = R_NaN; + return; + } + + double sum(0.0); + // Rf_rgamma (not the "rgamma" convenience macro) is used here + // deliberately: Rcpp's sugar layer #undefs that macro (Rcpp/sugar/ + // undoRmath.h), so this header needs to work whether or not it's + // included after Rcpp.h has already run. + for(i = 0; i < K; i++) { + if(alpha[i] <= 0.0) { + for(j = 0; j < K; j++) + ans[j] = R_NaN; + return; + } + ans[i] = Rf_rgamma(alpha[i], 1); + sum += ans[i]; + } + for(i = 0; i < K; i++) { + ans[i] /= sum; + } +} + +#endif // _NCOMPILER_DISTS_RDIRCH diff --git a/nCompiler/inst/include/nCompiler/dists/rdirch_tensor.h b/nCompiler/inst/include/nCompiler/dists/rdirch_tensor.h new file mode 100644 index 00000000..e4879657 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/rdirch_tensor.h @@ -0,0 +1,47 @@ +// Eigen::Tensor-argument wrapper around the rdirch random-generation +// kernel (see rdirch.h). Templated code like this can only ever be +// header-only -- the compiler must see the definition wherever it's +// instantiated -- so this lives in inst/include regardless of the +// package-DLL/on-the-fly-compilation split that motivates rdirch.h itself. + +#ifndef _NCOMPILER_DISTS_RDIRCH_TENSOR +#define _NCOMPILER_DISTS_RDIRCH_TENSOR + +#include +#include "dists_tensor_utils.h" +#include "rdirch.h" + +// Eigen::Tensor wrapper for rdirch. +// +// alpha is accepted as an arbitrary Eigen::Tensor expression built up at +// the call site. It's read-only (like in ddirch), so asDenseTensor reuses +// it directly (no copy) when it's already a concrete tensor, and only +// materializes a copy when the call site passed a lazy expression. The +// output is always a fresh buffer. +// +// No GetRNGstate()/PutRNGstate() bracketing here, matching the other +// r-function wrappers (rmnorm_chol_tensor.h, rmvt_chol_tensor.h, +// rwish_chol_tensor.h): code that reaches this wrapper is invoked from R +// through Rcpp, which already brackets RNG state around the outer .Call +// boundary (RNGScope). The plain-.Call DLL path (dists.cpp's C_rdirch) +// isn't routed through Rcpp, so it keeps its own bracketing -- that's +// unchanged and still correct. +template +Eigen::Tensor rdirch(int n, // ignored for now, to be used in the future. + const TensorExprAlpha &alpha) { + const auto &alphaEval = asDenseTensor<1>(alpha); + + int K = static_cast(alphaEval.dimension(0)); + + Eigen::Tensor ans(K); + + // alphaEval may alias the caller's own tensor (when no copy was needed + // above); rdirch's double* signature is non-const only because it + // predates const-correctness conventions here, but it never writes + // through alpha, so this cast is safe. + rdirch(ans.data(), const_cast(alphaEval.data()), K); + + return ans; +} + +#endif // _NCOMPILER_DISTS_RDIRCH_TENSOR diff --git a/nCompiler/inst/include/nCompiler/dists/rmnorm_chol.h b/nCompiler/inst/include/nCompiler/dists/rmnorm_chol.h new file mode 100644 index 00000000..d9f6f541 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/rmnorm_chol.h @@ -0,0 +1,73 @@ +// The rmnorm_chol (multivariate normal random generation) kernel (raw +// double* arguments), ported from nimble's dists.cpp. +// +// Header-only and self-contained for the same reason as dmnorm_chol.h: the +// package DLL and code nCompile generates and compiles on a user's machine +// are separate shared objects that cannot reliably link against each other, +// so each side compiles its own copy of this source instead. +// +// Unlike dmnorm_chol, there's no overwrite_inputs concept here: ans is +// always a fresh output buffer, never aliased to an input. +// +// See also rmnorm_chol_tensor.h for an Eigen::Tensor-argument wrapper +// around this kernel, for use from nCompile-generated code. That wrapper +// (like the original C_rmnorm_chol R-facing function) also brackets the +// call with GetRNGstate()/PutRNGstate() and applies the recycling rule for +// mean -- neither of those live in this raw kernel, matching how the +// original C_rmnorm_chol handled them at the R-call boundary rather than +// inside rmnorm_chol itself. + +#ifndef _NCOMPILER_DISTS_RMNORM_CHOL +#define _NCOMPILER_DISTS_RMNORM_CHOL + +// Based on Writing R Extensions 6.6.1 and mgcv's mgcv.h +#define USE_FC_LEN_T +#include +#include +/* If we are compiling with a version of R before FCONE and the explicit + supplying of extra arguments was introduced, then FCONE has to be + defined */ +#ifndef FCONE +#define FCONE +#endif + +#include // GetRNGstate, PutRNGstate, norm_rand +#include +#include "dists_utils.h" + +inline void rmnorm_chol(double *ans, double* mean, double* chol, int n, + double prec_param) { + char uplo('U'); + char transPrec('N'); + char transCov('T'); + char diag('N'); + int lda(n); + int incx(1); + + int i, j; + + if (ISNAN_ANY(mean, n) || ISNAN_ANY(chol, n*n) || ISNAN(prec_param)) { + for(j = 0; j < n; j++) + ans[j] = R_NaN; + return; + } + + if(!R_FINITE_ANY(chol, n*n)) { + for(j = 0; j < n; j++) + ans[j] = R_NaN; + return; + } + + for(i = 0; i < n; i++) + ans[i] = norm_rand(); + + // do upper-triangular solve or (transpose) multiply + // dtr{s,m}v is a BLAS level-2 function + if(prec_param) F77_CALL(dtrsv)(&uplo, &transPrec, &diag, &n, chol, &lda, ans, &incx FCONE FCONE FCONE); + else F77_CALL(dtrmv)(&uplo, &transCov, &diag, &n, chol, &lda, ans, &incx FCONE FCONE FCONE); + + for(i = 0; i < n; i++) + ans[i] += mean[i]; +} + +#endif // _NCOMPILER_DISTS_RMNORM_CHOL diff --git a/nCompiler/inst/include/nCompiler/dists/rmnorm_chol_tensor.h b/nCompiler/inst/include/nCompiler/dists/rmnorm_chol_tensor.h new file mode 100644 index 00000000..76ac099e --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/rmnorm_chol_tensor.h @@ -0,0 +1,60 @@ +// Eigen::Tensor-argument wrapper around the rmnorm_chol random-generation +// kernel (see rmnorm_chol.h). Templated code like this can only ever be +// header-only -- the compiler must see the definition wherever it's +// instantiated -- so this lives in inst/include regardless of the +// package-DLL/on-the-fly-compilation split that motivates rmnorm_chol.h +// itself. + +#ifndef _NCOMPILER_DISTS_RMNORM_CHOL_TENSOR +#define _NCOMPILER_DISTS_RMNORM_CHOL_TENSOR + +#include +#include "dists_tensor_utils.h" +#include "rmnorm_chol.h" + +// Eigen::Tensor wrapper for rmnorm_chol. +// +// mean and chol are accepted as arbitrary Eigen::Tensor expressions built +// up at the call site. chol is read-only, so asDenseTensor reuses it +// directly (no copy) when it's already a concrete tensor, and only +// materializes a copy when the call site passed a lazy expression. mean +// follows nimble's recycling rule (see recycleToLength in +// dists_tensor_utils.h): if it's shorter than the output length n +// (determined by chol's dimension), its elements are reused cyclically to +// fill out length n, matching the original C_rmnorm_chol's full_mean +// construction. +// +// Unlike dmnorm_chol_tensor.h, there's no x/overwrite_inputs concern here: +// the output is always a fresh buffer. +// +// No GetRNGstate()/PutRNGstate() bracketing here, unlike the original +// C_rmnorm_chol: code that reaches this wrapper is invoked from R through +// Rcpp, which brackets RNG state around the outer .Call boundary already +// (RNGScope). Bracketing individual inner draws like this one would be +// both redundant and the wrong altitude for it. The plain-.Call DLL path +// (dists.cpp's C_rmnorm_chol) isn't routed through Rcpp, so it keeps its +// own bracketing -- that's unchanged and still correct. +template +Eigen::Tensor rmnorm_chol(int n, // ignored for now, to be used in the future. + const TensorExprMean &mean, + const TensorExprChol &chol, + double prec_param) { + const auto &cholEval = asDenseTensor<2>(chol); + int len = static_cast(cholEval.dimension(0)); + + Eigen::Tensor meanEval = + recycleToLength(asDenseTensor<1>(mean), len); + + Eigen::Tensor ans(len); + + // cholEval may alias the caller's own tensor (when no copy was needed + // above); rmnorm_chol's double* signature is non-const only because it + // predates const-correctness conventions here, but it never writes + // through chol, so this cast is safe. + rmnorm_chol(ans.data(), meanEval.data(), + const_cast(cholEval.data()), len, prec_param); + + return ans; +} + +#endif // _NCOMPILER_DISTS_RMNORM_CHOL_TENSOR diff --git a/nCompiler/inst/include/nCompiler/dists/rmulti.h b/nCompiler/inst/include/nCompiler/dists/rmulti.h new file mode 100644 index 00000000..6ada7e97 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/rmulti.h @@ -0,0 +1,66 @@ +// The rmulti (multinomial random generation) kernel (raw arguments), ported +// from nimble's dists.cpp. +// +// Header-only and self-contained for the same reason as dmnorm_chol.h: the +// package DLL and code nCompile generates and compiles on a user's machine +// are separate shared objects that cannot reliably link against each other, +// so each side compiles its own copy of this source instead. +// +// Unlike every other dist kernel ported so far, the output here is int*, +// not double*: Rf_rmultinom (the underlying Rmath routine) requires an int +// output buffer. As the original comment below notes, calling functions +// need to convert to/from double themselves -- see rmulti_tensor.h, which +// does exactly that to match the double-based convention of every other +// wrapper in this directory (matching how nimble's own nimArr_rmulti also +// converts). + +#ifndef _NCOMPILER_DISTS_RMULTI +#define _NCOMPILER_DISTS_RMULTI + +#include // Rf_rmultinom +#include +#include "dists_utils.h" + +// Calling functions need to copy first arg back and forth to double if +// needed. +// IMPORTANT: ans and size are int when sent to Rf_rmultinom, since Rmath's +// rmultinom has these types. +inline void rmulti(int* ans, double size, double* prob, int K) { + double sumProb = 0.0; + int i; + + // The original nimble code assigned R_NaN (a double) directly into this + // int* array here, relying on the resulting int-truncation-of-NaN to + // land on NA_INTEGER -- which it reliably does on the platforms nimble + // targets, but is technically undefined behavior in C++. NA_INTEGER is + // used explicitly instead, which is well-defined and produces the exact + // same value (NA_INTEGER is INT_MIN, same as the original cast's result). + if (ISNAN_ANY(prob, K) || ISNAN(size)) { + for(i = 0; i < K; i++) + ans[i] = NA_INTEGER; + return; + } + + for(i = 0; i < K; i++) { + if(prob[i] < 0) { + for(i = 0; i < K; i++) + ans[i] = NA_INTEGER; + return; + } + sumProb += prob[i]; + } + if (sumProb <= 0.0) { // given above check for neg probs, this now will only catch '== 0' cases + for(i = 0; i < K; i++) + ans[i] = NA_INTEGER; + return; + } + for(i = 0; i < K; i++) + prob[i] /= sumProb; + // Rf_rmultinom (not the "rmultinom" convenience macro) is used here + // deliberately: Rcpp's sugar layer #undefs that macro (Rcpp/sugar/ + // undoRmath.h), so this header needs to work whether or not it's + // included after Rcpp.h has already run. + Rf_rmultinom((int) size, prob, K, ans); +} + +#endif // _NCOMPILER_DISTS_RMULTI diff --git a/nCompiler/inst/include/nCompiler/dists/rmulti_tensor.h b/nCompiler/inst/include/nCompiler/dists/rmulti_tensor.h new file mode 100644 index 00000000..16fdeef0 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/rmulti_tensor.h @@ -0,0 +1,55 @@ +// Eigen::Tensor-argument wrapper around the rmulti random-generation +// kernel (see rmulti.h). Templated code like this can only ever be +// header-only -- the compiler must see the definition wherever it's +// instantiated -- so this lives in inst/include regardless of the +// package-DLL/on-the-fly-compilation split that motivates rmulti.h itself. + +#ifndef _NCOMPILER_DISTS_RMULTI_TENSOR +#define _NCOMPILER_DISTS_RMULTI_TENSOR + +#include +#include "dists_tensor_utils.h" +#include "rmulti.h" + +// Eigen::Tensor wrapper for rmulti. +// +// prob is accepted as an arbitrary Eigen::Tensor expression built up at the +// call site, but -- unlike the read-only prob arguments elsewhere in this +// directory -- rmulti's kernel normalizes its working copy of prob in +// place (prob[i] /= sumProb) before calling Rf_rmultinom. So prob is +// always materialized into a private buffer here (never asDenseTensor's +// no-copy reference path), the same treatment given to arguments that get +// mutated elsewhere (e.g. chol in rwish_chol_tensor.h). +// +// rmulti's raw kernel writes into an int* buffer (Rf_rmultinom requires +// one), but every other wrapper in this directory works in double, so this +// wrapper draws into a temporary int buffer and converts to +// Eigen::Tensor for the result, mapping NA_INTEGER to R_NaN +// (matching the R_NaN used by every other dist kernel here to signal +// invalid input, rather than R's separate NA_REAL) -- matching what +// nimble's own nimArr_rmulti does at its call site. +// +// No GetRNGstate()/PutRNGstate() bracketing here, matching the other +// r-function wrappers: code that reaches this wrapper is invoked from R +// through Rcpp, which already brackets RNG state around the outer .Call +// boundary (RNGScope). The plain-.Call DLL path (dists.cpp's C_rmulti) +// isn't routed through Rcpp, so it keeps its own bracketing -- that's +// unchanged and still correct. +template +Eigen::Tensor rmulti(int n, // ignored for now, to be used in the future. + double size, const TensorExprProb &prob) { + Eigen::Tensor probEval = prob.template cast(); + + int K = static_cast(probEval.dimension(0)); + + Eigen::Tensor ansInt(K); + rmulti(ansInt.data(), size, probEval.data(), K); + + Eigen::Tensor ans(K); + for (int i = 0; i < K; i++) + ans(i) = (ansInt(i) == NA_INTEGER) ? R_NaN : static_cast(ansInt(i)); + + return ans; +} + +#endif // _NCOMPILER_DISTS_RMULTI_TENSOR diff --git a/nCompiler/inst/include/nCompiler/dists/rmvt_chol.h b/nCompiler/inst/include/nCompiler/dists/rmvt_chol.h new file mode 100644 index 00000000..1497f9b6 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/rmvt_chol.h @@ -0,0 +1,79 @@ +// The rmvt_chol (multivariate t random generation) kernel (raw double* +// arguments), ported from nimble's dists.cpp. +// +// Header-only and self-contained for the same reason as dmnorm_chol.h: the +// package DLL and code nCompile generates and compiles on a user's machine +// are separate shared objects that cannot reliably link against each other, +// so each side compiles its own copy of this source instead. +// +// Like rmnorm_chol, there's no overwrite_inputs concept here: ans is +// always a fresh output buffer, never aliased to an input. +// +// See also rmvt_chol_tensor.h for an Eigen::Tensor-argument wrapper around +// this kernel, for use from nCompile-generated code. That wrapper (like the +// original C_rmvt_chol R-facing function) applies the recycling rule for +// mu, which doesn't live in this raw kernel, matching how the original +// C_rmvt_chol handled it at the R-call boundary rather than inside +// rmvt_chol itself. + +#ifndef _NCOMPILER_DISTS_RMVT_CHOL +#define _NCOMPILER_DISTS_RMVT_CHOL + +// Based on Writing R Extensions 6.6.1 and mgcv's mgcv.h +#define USE_FC_LEN_T +#include +#include +/* If we are compiling with a version of R before FCONE and the explicit + supplying of extra arguments was introduced, then FCONE has to be + defined */ +#ifndef FCONE +#define FCONE +#endif + +#include // norm_rand +#include +#include "dists_utils.h" + +inline void rmvt_chol(double *ans, double* mu, double* chol, double df, + int n, double prec_param) { + char uplo('U'); + char transPrec('N'); + char transCov('T'); + char diag('N'); + int lda(n); + int incx(1); + + int i, j; + + if (ISNAN_ANY(mu, n) || ISNAN_ANY(chol, n*n) || ISNAN(df) || ISNAN(prec_param)) { + for(j = 0; j < n; j++) + ans[j] = R_NaN; + return; + } + + if(!R_FINITE_ANY(chol, n*n)) { + for(j = 0; j < n; j++) + ans[j] = R_NaN; + return; + } + + for(i = 0; i < n; i++) + ans[i] = norm_rand(); + + // sample from chi-squared and calculate scaling factor. + // Rf_rchisq (not the "rchisq" convenience macro) is used here + // deliberately: Rcpp's sugar layer #undefs that macro (Rcpp/sugar/ + // undoRmath.h), so this header needs to work whether or not it's + // included after Rcpp.h has already run. + double scaling = sqrt(df / Rf_rchisq(df)); + + // do upper-triangular solve or (transpose) multiply + // dtr{s,m}v is a BLAS level-2 function + if(prec_param) F77_CALL(dtrsv)(&uplo, &transPrec, &diag, &n, chol, &lda, ans, &incx FCONE FCONE FCONE); + else F77_CALL(dtrmv)(&uplo, &transCov, &diag, &n, chol, &lda, ans, &incx FCONE FCONE FCONE); + + for(i = 0; i < n; i++) + ans[i] = mu[i] + ans[i] * scaling; +} + +#endif // _NCOMPILER_DISTS_RMVT_CHOL diff --git a/nCompiler/inst/include/nCompiler/dists/rmvt_chol_tensor.h b/nCompiler/inst/include/nCompiler/dists/rmvt_chol_tensor.h new file mode 100644 index 00000000..d15bb062 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/rmvt_chol_tensor.h @@ -0,0 +1,58 @@ +// Eigen::Tensor-argument wrapper around the rmvt_chol random-generation +// kernel (see rmvt_chol.h). Templated code like this can only ever be +// header-only -- the compiler must see the definition wherever it's +// instantiated -- so this lives in inst/include regardless of the +// package-DLL/on-the-fly-compilation split that motivates rmvt_chol.h +// itself. + +#ifndef _NCOMPILER_DISTS_RMVT_CHOL_TENSOR +#define _NCOMPILER_DISTS_RMVT_CHOL_TENSOR + +#include +#include "dists_tensor_utils.h" +#include "rmvt_chol.h" + +// Eigen::Tensor wrapper for rmvt_chol. +// +// mu and chol are accepted as arbitrary Eigen::Tensor expressions built up +// at the call site. chol is read-only, so asDenseTensor reuses it directly +// (no copy) when it's already a concrete tensor, and only materializes a +// copy when the call site passed a lazy expression. mu follows nimble's +// recycling rule (see recycleToLength in dists_tensor_utils.h): if it's +// shorter than the output length n (determined by chol's dimension), its +// elements are reused cyclically to fill out length n, matching the +// original C_rmvt_chol's full_mu construction. +// +// Unlike dmvt_chol_tensor.h, there's no x/overwrite_inputs concern here: +// the output is always a fresh buffer. +// +// No GetRNGstate()/PutRNGstate() bracketing here, matching +// rmnorm_chol_tensor.h: code that reaches this wrapper is invoked from R +// through Rcpp, which already brackets RNG state around the outer .Call +// boundary (RNGScope). The plain-.Call DLL path (dists.cpp's C_rmvt_chol) +// isn't routed through Rcpp, so it keeps its own bracketing -- that's +// unchanged and still correct. +template +Eigen::Tensor rmvt_chol(int n, // ignored for now, to be used in the future. +const TensorExprMu &mu, + const TensorExprChol &chol, double df, + double prec_param) { + const auto &cholEval = asDenseTensor<2>(chol); + int len = static_cast(cholEval.dimension(0)); + + Eigen::Tensor muEval = + recycleToLength(asDenseTensor<1>(mu), len); + + Eigen::Tensor ans(len); + + // cholEval may alias the caller's own tensor (when no copy was needed + // above); rmvt_chol's double* signature is non-const only because it + // predates const-correctness conventions here, but it never writes + // through chol, so this cast is safe. + rmvt_chol(ans.data(), muEval.data(), + const_cast(cholEval.data()), df, len, prec_param); + + return ans; +} + +#endif // _NCOMPILER_DISTS_RMVT_CHOL_TENSOR diff --git a/nCompiler/inst/include/nCompiler/dists/rwish_chol.h b/nCompiler/inst/include/nCompiler/dists/rwish_chol.h new file mode 100644 index 00000000..2f68a41c --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/rwish_chol.h @@ -0,0 +1,110 @@ +// The rwish_chol (Wishart random generation) kernel (raw double* +// arguments), ported from nimble's dists.cpp. Like dwish_chol, the output +// Z here is itself a p x p matrix, not a vector. +// +// Header-only and self-contained for the same reason as dmnorm_chol.h: the +// package DLL and code nCompile generates and compiles on a user's machine +// are separate shared objects that cannot reliably link against each other, +// so each side compiles its own copy of this source instead. +// +// Unlike rmnorm_chol/rmvt_chol, this kernel does have an overwrite_inputs +// parameter: it controls whether chol itself is used as scratch space +// during the computation (overwrite_inputs = 1) or an internal copy is +// allocated (overwrite_inputs = 0). The original C_rwish_chol always +// passed 0, since chol there is the user's own R object and must not be +// mutated. See rwish_chol_tensor.h for how the wrapper handles this. +// +// See also rwish_chol_tensor.h for an Eigen::Tensor-argument wrapper +// around this kernel, for use from nCompile-generated code. + +#ifndef _NCOMPILER_DISTS_RWISH_CHOL +#define _NCOMPILER_DISTS_RWISH_CHOL + +// Based on Writing R Extensions 6.6.1 and mgcv's mgcv.h +#define USE_FC_LEN_T +#include +#include +/* If we are compiling with a version of R before FCONE and the explicit + supplying of extra arguments was introduced, then FCONE has to be + defined */ +#ifndef FCONE +#define FCONE +#endif + +#include // norm_rand +#include +#include "dists_utils.h" + +inline void rwish_chol(double *Z, double* chol, double df, int p, + double scale_param, int overwrite_inputs) { + char uplo('U'); + char sideL('L'); + char diag('N'); + char transT('T'); + char transN('N'); + double alpha(1.0); + double beta(0.0); + + double* cholCopy; + int i, j, uind, lind; + + if (ISNAN_ANY(chol, p*p) || ISNAN(df) || ISNAN(scale_param)) { + for(j = 0; j < p*p; j++) + Z[j] = R_NaN; + return; + } + + // also covers df < 0 + if(df < (double) p) { + for(j = 0; j < p*p; j++) + Z[j] = R_NaN; + return; + } + + // fill diags with sqrts of chi-squares and upper triangle (for scale_param) with std normals - crossproduct of result is standardized Wishart; based on rWishart in stats package + // + // Rf_rchisq (not the "rchisq" convenience macro) is used here + // deliberately: Rcpp's sugar layer #undefs that macro (Rcpp/sugar/ + // undoRmath.h), so this header needs to work whether or not it's + // included after Rcpp.h has already run. + for(j = 0; j < p; j++) { + Z[j*p + j] = sqrt(Rf_rchisq(df - (double) j)); + for(i = 0; i < j; i++) { + uind = i + j * p; /* upper triangle index */ + lind = j + i * p; /* lower triangle index */ + Z[(scale_param ? uind : lind)] = norm_rand(); + Z[(scale_param ? lind : uind)] = 0; + } + } + + // multiply Z*chol, both upper triangular or solve(chol, Z^T) + // would be more efficient if make use of fact that right-most matrix is triangular, but no available BLAS routine and hand-coding would eliminate use of threading and might well not be faster + if(overwrite_inputs) + cholCopy = chol; + else { + cholCopy = new double[p*p]; + if(scale_param) + for(i = 0; i < p*p; i++) + cholCopy[i] = chol[i]; + } + if(scale_param) F77_CALL(dtrmm)(&sideL, &uplo, &transN, &diag, &p, &p, &alpha, Z, &p, cholCopy, &p FCONE FCONE FCONE FCONE); + else F77_CALL(dtrsm)(&sideL, &uplo, &transN, &diag, &p, &p, &alpha, chol, &p, Z, &p FCONE FCONE FCONE FCONE); + + // cp result to Z or chol so can be used as matrix to multiply against and overwrite + if(scale_param) { + for(j = 0; j < p*p; j++) + Z[j] = cholCopy[j]; + } else { + for(j = 0; j < p*p; j++) + cholCopy[j] = Z[j]; + } + + // do crossprod of result + // for dtrmm call, again this would be more efficient if use fact that RHS upper triangular, but no available BLAS routine and hand-coding would eliminate use of threading and might well not be faster + if(scale_param) F77_CALL(dtrmm)(&sideL, &uplo, &transT, &diag, &p, &p, &alpha, cholCopy, &p, Z, &p FCONE FCONE FCONE FCONE); + else F77_CALL(dgemm)(&transN, &transT, &p, &p, &p, &alpha, cholCopy, &p, cholCopy, &p, &beta, Z, &p FCONE FCONE); + if(!overwrite_inputs) + delete [] cholCopy; +} + +#endif // _NCOMPILER_DISTS_RWISH_CHOL diff --git a/nCompiler/inst/include/nCompiler/dists/rwish_chol_tensor.h b/nCompiler/inst/include/nCompiler/dists/rwish_chol_tensor.h new file mode 100644 index 00000000..b983f459 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/dists/rwish_chol_tensor.h @@ -0,0 +1,51 @@ +// Eigen::Tensor-argument wrapper around the rwish_chol random-generation +// kernel (see rwish_chol.h). Templated code like this can only ever be +// header-only -- the compiler must see the definition wherever it's +// instantiated -- so this lives in inst/include regardless of the +// package-DLL/on-the-fly-compilation split that motivates rwish_chol.h +// itself. + +#ifndef _NCOMPILER_DISTS_RWISH_CHOL_TENSOR +#define _NCOMPILER_DISTS_RWISH_CHOL_TENSOR + +#include +#include "dists_tensor_utils.h" +#include "rwish_chol.h" + +// Eigen::Tensor wrapper for rwish_chol. +// +// chol is accepted as an arbitrary Eigen::Tensor expression built up at the +// call site. Unlike the read-only chol arguments in the density-function +// wrappers, rwish_chol can use chol itself as scratch space when +// overwrite_inputs = 1 (see rwish_chol.h), so chol is always materialized +// into a private buffer here (never asDenseTensor's no-copy reference +// path) before being handed to the kernel with overwrite_inputs = 1 -- +// matching how the original C_rwish_chol instead passed overwrite_inputs = +// 0 to force the kernel to make its own internal copy, since there chol +// aliased the user's own R object and could not be touched. +// +// The output p x p matrix Z is always a fresh buffer. +// +// No GetRNGstate()/PutRNGstate() bracketing here, matching +// rmnorm_chol_tensor.h/rmvt_chol_tensor.h: code that reaches this wrapper +// is invoked from R through Rcpp, which already brackets RNG state around +// the outer .Call boundary (RNGScope). The plain-.Call DLL path +// (dists.cpp's C_rwish_chol) isn't routed through Rcpp, so it keeps its +// own bracketing -- that's unchanged and still correct. +template +Eigen::Tensor rwish_chol(int n, // ignored for now, to be used in the future. +const TensorExprChol &chol, double df, + double scale_param) { + Eigen::Tensor cholEval = chol.template cast(); + + int p = static_cast(cholEval.dimension(0)); + + Eigen::Tensor Z(p, p); + + rwish_chol(Z.data(), cholEval.data(), df, p, scale_param, + /*overwrite_inputs=*/1); + + return Z; +} + +#endif // _NCOMPILER_DISTS_RWISH_CHOL_TENSOR diff --git a/nCompiler/inst/include/nCompiler/nCompiler_omnibus.h b/nCompiler/inst/include/nCompiler/nCompiler_omnibus.h index a565860c..cbbbf341 100644 --- a/nCompiler/inst/include/nCompiler/nCompiler_omnibus.h +++ b/nCompiler/inst/include/nCompiler/nCompiler_omnibus.h @@ -52,6 +52,14 @@ using namespace Rcpp; #include #endif +// dists (e.g. dmnorm_chol and its Eigen::Tensor wrapper) are template-only, +// so an unused one costs nothing beyond being parsed. NCOMPILER_USES_MVDISTS +// is unconditionally defined in cppDefs_nFunction.R/cppDefs_nClass.R for +// now; a future pass can set it only when generated code actually needs it. +#ifdef NCOMPILER_USES_MVDISTS +#include +#endif + #ifdef NCOMPILER_USES_TBB #include diff --git a/nCompiler/src/dists.cpp b/nCompiler/src/dists.cpp index d30aef3f..2d5afc74 100644 --- a/nCompiler/src/dists.cpp +++ b/nCompiler/src/dists.cpp @@ -29,140 +29,13 @@ #include "dists.h" #include -// Detects NA -bool R_IsNA_ANY(double* P, int s) { - for(int i = 0; i < s; ++i) if(R_IsNA(P[i])) return(true); - return(false); -} - -// Detects NaN -bool R_IsNaN_ANY(double* P, int s) { - for(int i = 0; i < s; ++i) if(R_IsNaN(P[i])) return(true); - return(false); -} - -// Detects NA or NaN -bool ISNAN_ANY(double* P, int s) { - for(int i = 0; i < s; ++i) if(ISNAN(P[i])) return(true); - return(false); -} - -bool R_FINITE_ANY(double* P, int s) { - for(int i = 0; i < s; ++i) if(!R_FINITE(P[i])) return(false); - return(true); -} - -double dwish_chol(double* x, double* chol, double df, int p, double scale_param, int give_log, int overwrite_inputs) { - char uplo('U'); - char sideL('L'); - char sideR('R'); - char diag('N'); - char transN('N'); - int info(0); - double alpha(1.0); - double* xChol; - - int i, j; - - if (R_IsNA_ANY(x, p*p) || R_IsNA_ANY(chol, p*p) || R_IsNA(df) || R_IsNA(scale_param)) - return NA_REAL; - if (R_IsNaN_ANY(x, p*p) || R_IsNaN_ANY(chol, p*p) || R_IsNaN(df) || R_IsNaN(scale_param)) - return R_NaN; - - // also covers df < 0 - if(df < (double) p) ML_ERR_return_NAN; - - if(!R_FINITE_ANY(x, p*p) || !R_FINITE_ANY(chol, p*p)) return R_D__0; - - double dens = -(df*p/2 * M_LN2 + p*(p-1)*M_LN_SQRT_PI/2); - for(i = 0; i < p; i++) - dens -= lgammafn((df - i) / 2); - - if(scale_param) { - for(i = 0; i < p*p; i += p + 1) - dens -= df * log(chol[i]); - } else { - for(i = 0; i < p*p; i += p + 1) - dens += df * log(chol[i]); - } - - // determinant of x using Cholesky: - if(overwrite_inputs && (int) scale_param) // if !scale_param we need x below - xChol = x; - else { - xChol = new double[p*p]; - // only need upper triangle for dpotrf chol calculation - for(j = 0; j < p; j++) - for(i = 0; i <= j; i++) - xChol[j*p+i] = x[j*p+i]; - } - F77_CALL(dpotrf)(&uplo, &p, xChol, &p, &info FCONE); - for(i = 0; i < p*p; i += p + 1) - dens += (df - p - 1) * log(xChol[i]); - - // R %*% x = t(chol) %*% chol %*% x (could also do with chol(x) but no more efficient - // solve(S, x) = crossproduct( chol(x) %*% inverse(chol) ) - - // dtr{m,s}m is a BLAS level-3 function - double tmp_dens = 0.0; - if(scale_param) { - // chol(x) %*% inverse(chol) - // need lower triangle of xChol to have zeros as dtrsm assumes it is full matrix - for(j = 0; j < p-1; j++) - for(i = j+1; i < p; i++) - xChol[j*p+i] = 0.0; - F77_CALL(dtrsm)(&sideR, &uplo, &transN, &diag, &p, &p, &alpha, - chol, &p, xChol, &p FCONE FCONE FCONE FCONE); - // trace of crossproduct of result is sum of squares of elements - for(j = 0; j < p; j++) - for(i = 0; i <= j; i++) - tmp_dens += xChol[j*p+i]*xChol[j*p+i]; - } else { - double* xCopy; - if(overwrite_inputs) - xCopy = x; - else { - xCopy = new double[p*p]; - for(i = 0; i < p*p; i++) - xCopy[i] = x[i]; - } - // chol %*% x - F77_CALL(dtrmm)(&sideL, &uplo, &transN, &diag, &p, &p, &alpha, - chol, &p, xCopy, &p FCONE FCONE FCONE FCONE); - // trace crossproduct of t(chol) with result is sum of product of upper-triangular elements - for(j = 0; j < p; j++) { - for(i = 0; i <= j; i++) { - tmp_dens += xCopy[j*p+i] * chol[j*p+i]; - } - } - if(!overwrite_inputs) - delete [] xCopy; - } - - if(!(overwrite_inputs && (int) scale_param)) - delete [] xChol; +// R_IsNA_ANY, R_IsNaN_ANY, ISNAN_ANY, and R_FINITE_ANY now come from +// nCompiler/dists/dists_utils.h (included via dists.h). - // attempt to improve above calcs by doing efficient U^T U multiply followed by direct product multiply, however this would not make use of threading provided by BLAS and even with one thread seems to be no faster - // U^T*U directly followed by direct product with x - /* - double tmp_summand; - int minij; - for(j = 0; j < p; j++) - for(i = 0; i < p; i++) { - tmp_summand = 0.0; - minij = i <= j ? i : j; - for(int k = 0; k < minij; k++) - tmp_summand += chol[j*p+k]*chol[i*p+k]; // U^T U - // double if not on diagonal to account for direct product of lower triangle too - if(i != j) tmp_summand *= 2; - tmp_dens += xCopy[j*p+i] * tmp_summand; - } - */ - - dens += -0.5 * tmp_dens; - - return give_log ? dens : exp(dens); -} +// dwish_chol (raw double* kernel) now comes from nCompiler/dists/dwish_chol.h +// (included via dists.h). The Eigen::Tensor-argument wrapper around it +// lives in nCompiler/dists/dwish_chol_tensor.h, for use from +// nCompile-generated code; it's not needed here. SEXP C_dwish_chol(SEXP x, SEXP chol, SEXP df, SEXP scale_param, SEXP return_log) @@ -200,79 +73,11 @@ SEXP C_dwish_chol(SEXP x, SEXP chol, SEXP df, SEXP scale_param, SEXP return_log) } -void rwish_chol(double *Z, double* chol, double df, int p, double scale_param, int overwrite_inputs) { - char uplo('U'); - char sideL('L'); - char diag('N'); - char transT('T'); - char transN('N'); - double alpha(1.0); - double beta(0.0); - - double* cholCopy; - int i, j, uind, lind; - - if (ISNAN_ANY(chol, p*p) || ISNAN(df) || ISNAN(scale_param)) { - for(j = 0; j < p*p; j++) - Z[j] = R_NaN; - return; - } - - // also covers df < 0 - if(df < (double) p) { - for(j = 0; j < p*p; j++) - Z[j] = R_NaN; - return; - } - - // fill diags with sqrts of chi-squares and upper triangle (for scale_param) with std normals - crossproduct of result is standardized Wishart; based on rWishart in stats package - for(j = 0; j < p; j++) { - // double *Z_j = &Z[j*p]; - //Z_j[j] = sqrt(rchisq(df - (double) j)); - Z[j*p + j] = sqrt(rchisq(df - (double) j)); - for(i = 0; i < j; i++) { - uind = i + j * p, /* upper triangle index */ - lind = j + i * p; /* lower triangle index */ - Z[(scale_param ? uind : lind)] = norm_rand(); - Z[(scale_param ? lind : uind)] = 0; - } - /* - for(i = 0; i < j; i++) - Z_j[i] = norm_rand(); - for (i = j + 1; i < p; i++) - Z_j[i] = 0; - */ - } - - // multiply Z*chol, both upper triangular or solve(chol, Z^T) - // would be more efficient if make use of fact that right-most matrix is triangular, but no available BLAS routine and hand-coding would eliminate use of threading and might well not be faster - if(overwrite_inputs) - cholCopy = chol; - else { - cholCopy = new double[p*p]; - if(scale_param) - for(i = 0; i < p*p; i++) - cholCopy[i] = chol[i]; - } - if(scale_param) F77_CALL(dtrmm)(&sideL, &uplo, &transN, &diag, &p, &p, &alpha, Z, &p, cholCopy, &p FCONE FCONE FCONE FCONE); - else F77_CALL(dtrsm)(&sideL, &uplo, &transN, &diag, &p, &p, &alpha, chol, &p, Z, &p FCONE FCONE FCONE FCONE); - - // cp result to Z or chol so can be used as matrix to multiply against and overwrite - if(scale_param) { - for(j = 0; j < p*p; j++) - Z[j] = cholCopy[j]; - } else { - for(j = 0; j < p*p; j++) - cholCopy[j] = Z[j]; - } - - // do crossprod of result - // for dtrmm call, again this would be more efficient if use fact that RHS upper triangular, but no available BLAS routine and hand-coding would eliminate use of threading and might well not be faster - if(scale_param) F77_CALL(dtrmm)(&sideL, &uplo, &transT, &diag, &p, &p, &alpha, cholCopy, &p, Z, &p FCONE FCONE FCONE FCONE); - else F77_CALL(dgemm)(&transN, &transT, &p, &p, &p, &alpha, cholCopy, &p, cholCopy, &p, &beta, Z, &p FCONE FCONE); - if(!overwrite_inputs) - delete [] cholCopy; -} +// rwish_chol (raw double* kernel) now comes from +// nCompiler/dists/rwish_chol.h (included via dists.h). The +// Eigen::Tensor-argument wrapper around it lives in +// nCompiler/dists/rwish_chol_tensor.h, for use from nCompile-generated +// code; it's not needed here. SEXP C_rwish_chol(SEXP chol, SEXP df, SEXP scale_param) // generates single Wishart draw given Cholesky of scale or rate matrix @@ -665,58 +470,15 @@ SEXP C_rlkj_corr_cholesky(SEXP eta, SEXP p) } -double ddirch(double* x, double* alpha, int K, int give_log) -// scalar function that can be called directly by NIMBLE with same name as in R -{ - double sumAlpha(0.0); - double sumX(0.0); - double dens(0.0); - - if (R_IsNA_ANY(x, K) || R_IsNA_ANY(alpha, K)) - return NA_REAL; - if (R_IsNaN_ANY(x, K) || R_IsNaN_ANY(alpha, K)) - return R_NaN; - - for(int i = 0; i < K; i++) { - if(alpha[i] <= 0.0) ML_ERR_return_NAN; - if(x[i] < 0.0 || x[i] > 1.0) return R_D__0; - dens += (alpha[i]-1) * log(x[i]) - lgammafn(alpha[i]) ; - sumAlpha += alpha[i]; - sumX += x[i]; - } - if(sumX > 1.0 + 10*DBL_EPSILON || sumX < 1.0 - 10*DBL_EPSILON) { - return R_D__0; - } - - dens += lgammafn(sumAlpha); - return give_log ? dens : exp(dens); -} - -void rdirch(double* ans, double* alpha, int K) -// scalar function that can be called directly by NIMBLE with same name as in R -{ - int i, j; - - if (ISNAN_ANY(alpha, K)) { - for(j = 0; j < K; j++) - ans[j] = R_NaN; - return; - } +// ddirch (raw double* kernel) now comes from nCompiler/dists/ddirch.h +// (included via dists.h). The Eigen::Tensor-argument wrapper around it +// lives in nCompiler/dists/ddirch_tensor.h, for use from nCompile-generated +// code; it's not needed here. - double sum(0.0); - for(i = 0; i < K; i++) { - if(alpha[i] <= 0.0) { - for(j = 0; j < K; j++) - ans[j] = R_NaN; - return; - } - ans[i] = rgamma(alpha[i], 1); - sum += ans[i]; - } - for(i = 0; i < K; i++) { - ans[i] /= sum; - } -} +// rdirch (raw double* kernel) now comes from nCompiler/dists/rdirch.h +// (included via dists.h). The Eigen::Tensor-argument wrapper around it +// lives in nCompiler/dists/rdirch_tensor.h, for use from nCompile-generated +// code; it's not needed here. SEXP C_ddirch(SEXP x, SEXP alpha, SEXP return_log) { @@ -769,84 +531,16 @@ SEXP C_rdirch(SEXP alpha) { -double dmulti(double* x, double size, double* prob, int K, int give_log) // Calling functions need to copy first arg to int if needed -// scalar function that can be called directly by NIMBLE with same name as in R -{ - double sumProb(0.0); - double sumX(0.0); - double logSumProb; - - if (R_IsNA_ANY(x, K) || R_IsNA_ANY(prob, K) || R_IsNA(size)) - return NA_REAL; - if (R_IsNaN_ANY(x, K) || R_IsNaN_ANY(prob, K) || R_IsNaN(size)) - return R_NaN; - - if(R_D_negInonint(size)) - ML_ERR_return_NAN; - size = R_D_forceint(size); - - double dens = lgammafn(size + 1); - for(int i = 0; i < K; i++) { - if (prob[i] < 0) ML_ERR_return_NAN; - R_D_nonint_check(x[i]); - if (x[i] < 0 || !R_FINITE(x[i])) return R_D__0; - - x[i] = R_D_forceint(x[i]); - sumProb += prob[i]; - sumX += x[i]; - } - logSumProb = log(sumProb); - - for(int i = 0; i < K; i++) { - if(!(x[i] == 0.0 && prob[i] == 0.0)) - dens += x[i]*(log(prob[i]) - logSumProb) - lgammafn(x[i] + 1); - } - - if(sumX > size + 10*DBL_EPSILON || sumX < size - 10*DBL_EPSILON) { - return R_D__0; - } - - return give_log ? dens : exp(dens); -} - - - +// dmulti (raw double* kernel) now comes from nCompiler/dists/dmulti.h +// (included via dists.h). The Eigen::Tensor-argument wrapper around it +// lives in nCompiler/dists/dmulti_tensor.h, for use from nCompile-generated +// code; it's not needed here. -void rmulti(int* ans, double size, double* prob, int K) // Calling functions need to copy first arg back and forth to double if needed -// scalar function that can be called directly by NIMBLE with same name as in R -// just call Rmath's rmultinom, which passes result by pointer -// IMPORTANT: have ans and size as int when sent to rmultinom as Rmath rmultinom has these types -// Nimble does a copy in nimArr_rmulti -{ - /* rmultinom requires normalized probs (in R this is done in Rmultinom interface - function via FixupProb before passing to rmultinom) */ - double sumProb = 0.0; - int i; - - if (ISNAN_ANY(prob, K) || ISNAN(size) ) { - for(i = 0; i < K; i++) - ans[i] = R_NaN; // but casting to integer in C_rmulti gives NA - return; - } - - for(i = 0; i < K; i++) { - if(prob[i] < 0) { - for(i = 0; i < K; i++) - ans[i] = R_NaN; - return; - } - sumProb += prob[i]; - } - if (sumProb <= 0.0) { // given above check for neg probs, this now will only catch '== 0' cases - for(i = 0; i < K; i++) - ans[i] = R_NaN; - return; - } - for(i = 0; i < K; i++) - prob[i] /= sumProb; - rmultinom((int) size, prob, K, ans); -} +// rmulti (raw kernel) now comes from nCompiler/dists/rmulti.h (included +// via dists.h). The Eigen::Tensor-argument wrapper around it lives in +// nCompiler/dists/rmulti_tensor.h, for use from nCompile-generated code; +// it's not needed here. SEXP C_dmulti(SEXP x, SEXP size, SEXP prob, SEXP return_log) { @@ -903,7 +597,7 @@ SEXP C_rmulti(SEXP size, SEXP prob) { GetRNGstate(); PROTECT(ans = Rf_allocVector(INTSXP, K)); - // note that if NaN set in rmulti, the INTEGER() casts it to NA + // note that rmulti sets NA_INTEGER directly for invalid input rmulti(INTEGER(ans), c_size, c_prob, K); PutRNGstate(); UNPROTECT(1); @@ -1023,64 +717,10 @@ SEXP C_rcat(SEXP n, SEXP prob) { -double dmnorm_chol(double* x, double* mean, double* chol, int n, double prec_param, int give_log, int overwrite_inputs) { - char uplo('U'); - char transPrec('N'); - char transCov('T'); - char diag('N'); - int lda(n); - int incx(1); - double* xCopy; - - double dens = -n * M_LN_SQRT_2PI; - int i; - // add diagonals of Cholesky - - if (R_IsNA_ANY(x, n) || R_IsNA_ANY(mean, n) || R_IsNA_ANY(chol, n*n) || R_IsNA(prec_param)) - return NA_REAL; - if (R_IsNaN_ANY(x, n) || R_IsNaN_ANY(mean, n) || R_IsNaN_ANY(chol, n*n) || R_IsNaN(prec_param)) - return R_NaN; - - if(!R_FINITE_ANY(x, n) || !R_FINITE_ANY(mean, n) || !R_FINITE_ANY(chol, n*n)) return R_D__0; - - - - if(prec_param) { - for(i = 0; i < n*n; i += n + 1) - dens += log(chol[i]); - } else { - for(i = 0; i < n*n; i += n + 1) - dens -= log(chol[i]); - } - - if(overwrite_inputs) { - xCopy = x; - for(i = 0; i < n; i++) - xCopy[i] -= mean[i]; - } else { - xCopy = new double[n]; - for(i = 0; i < n; i++) - xCopy[i] = x[i] - mean[i]; - } - - // do matrix-vector multiply with upper-triangular matrix stored column-wise as full n x n matrix (prec parameterization) - // or upper-triangular (transpose) solve (cov parameterization) - // dtr{m,s}v is a BLAS level-2 function - if(prec_param) F77_CALL(dtrmv)(&uplo, &transPrec, &diag, &n, chol, &lda, xCopy, &incx FCONE FCONE FCONE); - else F77_CALL(dtrsv)(&uplo, &transCov, &diag, &n, chol, &lda, xCopy, &incx FCONE FCONE FCONE); - - // sum of squares to calculate quadratic form - double tmp = 0.0; - for(i = 0; i < n; i++) - tmp += xCopy[i] * xCopy[i]; - - dens += -0.5 * tmp; - - if(!overwrite_inputs) - delete [] xCopy; - - return give_log ? dens : exp(dens); -} +// dmnorm_chol (raw double* kernel) now comes from +// nCompiler/dists/dmnorm_chol.h (included via dists.h). The Eigen::Tensor- +// argument wrapper around it lives in nCompiler/dists/dmnorm_chol_tensor.h, +// for use from nCompile-generated code; it's not needed here. SEXP C_dmnorm_chol(SEXP x, SEXP mean, SEXP chol, SEXP prec_param, SEXP return_log) // calculates mv normal density given Cholesky of precision matrix or covariance matrix @@ -1129,39 +769,11 @@ SEXP C_dmnorm_chol(SEXP x, SEXP mean, SEXP chol, SEXP prec_param, SEXP return_lo return ans; } -void rmnorm_chol(double *ans, double* mean, double* chol, int n, double prec_param) { - char uplo('U'); - char transPrec('N'); - char transCov('T'); - char diag('N'); - int lda(n); - int incx(1); - - int i, j; - - if (ISNAN_ANY(mean, n) || ISNAN_ANY(chol, n*n) || ISNAN(prec_param)) { - for(j = 0; j < n; j++) - ans[j] = R_NaN; - return; - } - - if(!R_FINITE_ANY(chol, n*n)) { - for(j = 0; j < n; j++) - ans[j] = R_NaN; - return; - } - - for(i = 0; i < n; i++) - ans[i] = norm_rand(); - - // do upper-triangular solve or (transpose) multiply - // dtr{s,m}v is a BLAS level-2 function - if(prec_param) F77_CALL(dtrsv)(&uplo, &transPrec, &diag, &n, chol, &lda, ans, &incx FCONE FCONE FCONE); - else F77_CALL(dtrmv)(&uplo, &transCov, &diag, &n, chol, &lda, ans, &incx FCONE FCONE FCONE); - - for(i = 0; i < n; i++) - ans[i] += mean[i]; -} +// rmnorm_chol (raw double* kernel) now comes from +// nCompiler/dists/rmnorm_chol.h (included via dists.h). The +// Eigen::Tensor-argument wrapper around it lives in +// nCompiler/dists/rmnorm_chol_tensor.h, for use from nCompile-generated +// code; it's not needed here. SEXP C_rmnorm_chol(SEXP mean, SEXP chol, SEXP prec_param) // generates single mv normal draw given Cholesky of precision matrix or covariance matrix @@ -1480,62 +1092,10 @@ if(!Rf_isReal(mean)) // ...existing code... // Begin multivariate t -double dmvt_chol(double* x, double* mu, double* chol, double df, int n, double prec_param, int give_log, int overwrite_inputs) { - char uplo('U'); - char transPrec('N'); - char transCov('T'); - char diag('N'); - int lda(n); - int incx(1); - double* xCopy; - - double dens = lgammafn((df + n) / 2) - lgammafn(df / 2) - n * M_LN_SQRT_PI - n * log(df) / 2; - int i; - - if (R_IsNA_ANY(x, n) || R_IsNA_ANY(mu, n) || R_IsNA_ANY(chol, n*n) || R_IsNA(df) || R_IsNA(prec_param)) - return NA_REAL; - if (R_IsNaN_ANY(x, n) || R_IsNaN_ANY(mu, n) || R_IsNaN_ANY(chol, n*n) || R_IsNA(df) || R_IsNaN(prec_param)) - return R_NaN; - - if(!R_FINITE_ANY(x, n) || !R_FINITE_ANY(mu, n) || !R_FINITE_ANY(chol, n*n)) return R_D__0; - - // add diagonals of Cholesky - if(prec_param) { - for(i = 0; i < n*n; i += n + 1) - dens += log(chol[i]); - } else { - for(i = 0; i < n*n; i += n + 1) - dens -= log(chol[i]); - } - - if(overwrite_inputs) { - xCopy = x; - for(i = 0; i < n; i++) - xCopy[i] -= mu[i]; - } else { - xCopy = new double[n]; - for(i = 0; i < n; i++) - xCopy[i] = x[i] - mu[i]; - } - - // do matrix-vector multiply with upper-triangular matrix stored column-wise as full n x n matrix (prec parameterization) - // or upper-triangular (transpose) solve (cov parameterization) - // dtr{m,s}v is a BLAS level-2 function - if(prec_param) F77_CALL(dtrmv)(&uplo, &transPrec, &diag, &n, chol, &lda, xCopy, &incx FCONE FCONE FCONE); - else F77_CALL(dtrsv)(&uplo, &transCov, &diag, &n, chol, &lda, xCopy, &incx FCONE FCONE FCONE); - - // sum of squares to calculate quadratic form - double tmp = 0.0; - for(i = 0; i < n; i++) - tmp += xCopy[i] * xCopy[i]; - - dens += -0.5 * (df + n) * log(1 + tmp / df); - - if(!overwrite_inputs) - delete [] xCopy; - - return give_log ? dens : exp(dens); -} +// dmvt_chol (raw double* kernel) now comes from nCompiler/dists/dmvt_chol.h +// (included via dists.h). The Eigen::Tensor-argument wrapper around it +// lives in nCompiler/dists/dmvt_chol_tensor.h, for use from +// nCompile-generated code; it's not needed here. SEXP C_dmvt_chol(SEXP x, SEXP mu, SEXP chol, SEXP df, SEXP prec_param, SEXP return_log) // calculates mv normal density given Cholesky of precision matrix or covariance matrix @@ -1585,42 +1145,10 @@ SEXP C_dmvt_chol(SEXP x, SEXP mu, SEXP chol, SEXP df, SEXP prec_param, SEXP retu return ans; } -void rmvt_chol(double *ans, double* mu, double* chol, double df, int n, double prec_param) { - char uplo('U'); - char transPrec('N'); - char transCov('T'); - char diag('N'); - int lda(n); - int incx(1); - - int i, j; - - if (ISNAN_ANY(mu, n) || ISNAN_ANY(chol, n*n) || ISNAN(df) || ISNAN(prec_param)) { - for(j = 0; j < n; j++) - ans[j] = R_NaN; - return; - } - - if(!R_FINITE_ANY(chol, n*n)) { - for(j = 0; j < n; j++) - ans[j] = R_NaN; - return; - } - - for(i = 0; i < n; i++) - ans[i] = norm_rand(); - - // sample from chi-squared and calculate scaling factor - double scaling = sqrt(df / rchisq(df)); - - // do upper-triangular solve or (transpose) multiply - // dtr{s,m}v is a BLAS level-2 function - if(prec_param) F77_CALL(dtrsv)(&uplo, &transPrec, &diag, &n, chol, &lda, ans, &incx FCONE FCONE FCONE); - else F77_CALL(dtrmv)(&uplo, &transCov, &diag, &n, chol, &lda, ans, &incx FCONE FCONE FCONE); - - for(i = 0; i < n; i++) - ans[i] = mu[i] + ans[i] * scaling; -} +// rmvt_chol (raw double* kernel) now comes from nCompiler/dists/rmvt_chol.h +// (included via dists.h). The Eigen::Tensor-argument wrapper around it +// lives in nCompiler/dists/rmvt_chol_tensor.h, for use from +// nCompile-generated code; it's not needed here. SEXP C_rmvt_chol(SEXP mu, SEXP chol, SEXP df, SEXP prec_param) // generates single mv normal draw given Cholesky of precision matrix or covariance matrix diff --git a/nCompiler/src/dists.h b/nCompiler/src/dists.h index f5a54ae8..8e98a907 100644 --- a/nCompiler/src/dists.h +++ b/nCompiler/src/dists.h @@ -34,14 +34,10 @@ #define __DISTS #include "Utils.h" +#include using std::max; -bool R_IsNA_ANY(double*, int); -bool R_IsNaN_ANY(double*, int); -bool ISNAN_ANY(double*, int); -bool R_FINITE_ANY(double*, int); - extern "C" { // BLAS/LAPACK routines // prototypes for dgemm_, dtrmv_, dtrsv_, dtrsm_, dtrmm_, dpotr_ not needed when #include "R_ext/lapack.h" @@ -91,20 +87,20 @@ extern "C" { // core scalar d/r functions provided by NIMBLE to extend R double dcat(double, double*, int, int); double rcat(double*, int); -double dmulti(double*, double, double*, int, int); -void rmulti(int*, double, double*, int); -double ddirch(double*, double*, int, int); -void rdirch(double*, double*, int); +#include +#include +#include +#include -double dmnorm_chol(double*, double*, double*, int, double, int, int); -void rmnorm_chol(double *, double*, double*, int, double); +#include +#include void PDinverse_logdet_internal(double *matPtr, double *ans, int n); double dmnorm_inv_ld(double*, double*, double*, double*, int, int, int, int); void rmnorm_inv_ld(double*, double*, double*, double*, int, int); -double dmvt_chol(double*, double*, double*, double, int, double, int, int); -void rmvt_chol(double *, double*, double*, double, int, double); -double dwish_chol(double*, double*, double, int, double, int, int); -void rwish_chol(double*, double*, double, int, double, int); +#include +#include +#include +#include double dinvwish_chol(double*, double*, double, int, double, int, int); void rinvwish_chol(double*, double*, double, int, double, int); double dlkj_corr_cholesky(double*, double, int, int); diff --git a/nCompiler/tests/testthat/math_tests/test-math.R b/nCompiler/tests/testthat/math_tests/test-math.R index 62159cf1..02e4ba0e 100644 --- a/nCompiler/tests/testthat/math_tests/test-math.R +++ b/nCompiler/tests/testthat/math_tests/test-math.R @@ -17,6 +17,7 @@ utils <- system.file( c('testing_utils.R', 'testing_operatorLists.R', if(respect_known_failures) 'known_failures.R' else character(), 'math_utils.R') + ), package = 'nCompiler' ) diff --git a/nCompiler/tests/testthat/specificOp_tests/test-mv-dists.R b/nCompiler/tests/testthat/specificOp_tests/test-mv-dists.R new file mode 100644 index 00000000..628ecdb8 --- /dev/null +++ b/nCompiler/tests/testthat/specificOp_tests/test-mv-dists.R @@ -0,0 +1,180 @@ +# Tests of multivariate distributions + +library(nCompiler) +library(testthat) + +test_that("multivariate dists work", { + nc <- nClass( + Cpublic = list( + dmnc = nFunction( + function(x = double(1), mean = double(1), cholesky = double(2), + prec_param = logical(), log = logical()) { + ans <- dmnorm_chol(x, mean, cholesky, prec_param, log) + return(ans); returnType(double()) + } + ), + rmnc = nFunction( + function(mean = double(1), cholesky = double(2), + prec_param = logical()) { + ans <- rmnorm_chol(n = 1, mean, cholesky, prec_param) + return(ans); returnType(double(1)) + } + ), + dmvtc = nFunction( + function(x = double(1), mu = double(1), cholesky = double(2), df = integer(), + prec_param = logical(), log = logical()) { + ans <- dmvt_chol(x, mu, cholesky, df, prec_param, log) + return(ans); returnType(double()) + } + ), + rmvtc = nFunction( + function(mu = double(1), cholesky = double(2), df = integer(), + prec_param = logical()) { + ans <- rmvt_chol(n = 1, mu, cholesky, df, prec_param) + return(ans); returnType(double(1)) + } + ), + dwishc = nFunction( + function(x = double(2), cholesky = double(2), df = integer(), + scale_param = logical(), log = logical()) { + ans <- dwish_chol(x, cholesky, df, scale_param, log) + return(ans); returnType(double()) + } + ), + rwishc = nFunction( + function(cholesky = double(2), df = integer(), + scale_param = logical()) { + ans <- rwish_chol(n = 1, cholesky, df, scale_param) + return(ans); returnType(double(2)) + } + ), + ddirch_ = nFunction( + function(x = double(1), alpha = double(1), + log = logical()){ + ans <- ddirch(x, alpha, log) + return(ans); returnType(double()) + } + ), + rdirch_ = nFunction( + function(alpha = double(1)){ + ans <- rdirch(n = 1, alpha) + return(ans); returnType(double(1)) + } + ), + dmulti_ = nFunction( + function(x = double(1), prob = double(1), + log = logical()){ + ans <- dmulti(x, size = sum(x), prob = prob, log = log) + return(ans); returnType(double()) + } + ), + rmulti_ = nFunction( + function(size = integer(), prob = double(1)){ + ans <- rmulti(n = 1, size = size, prob = prob) + return(ans); returnType(double(1)) + } + ) + ) + ) + + comp <- nCompile(nc) + set.seed(1) + x <- rnorm(3) + mean <- 1:3 + cov <- matrix(c(1, 0.5, 0.25, 0.5, 1, 0.5, 0.25, 0.5, 1), byrow = TRUE, nrow = 3) + chol_cov <- chol(cov) + prec <- solve(cov) + chol_prec <- chol(prec) + obj <- nc$new() + Cobj <- comp$new() + # dmnorm + expect_equal(obj$dmnc(x, mean, chol_cov, FALSE, FALSE), + Cobj$dmnc(x, mean, chol_cov, FALSE, FALSE)) + expect_equal(obj$dmnc(x, mean, chol_cov, FALSE, TRUE), + Cobj$dmnc(x, mean, chol_cov, FALSE, TRUE)) + expect_equal(obj$dmnc(x, mean, chol_prec, TRUE, FALSE), + Cobj$dmnc(x, mean, chol_prec, TRUE, FALSE)) + expect_equal(obj$dmnc(x, mean, chol_prec, TRUE, TRUE), + Cobj$dmnc(x, mean, chol_prec, TRUE, TRUE)) + expect_equal(Cobj$dmnc(x, mean, chol_cov, FALSE, TRUE), + Cobj$dmnc(x, mean, chol_prec, TRUE, TRUE)) + + # rmnorm + expect_equal( + {set.seed(10); obj$rmnc(mean, chol_cov, FALSE)}, + {set.seed(10); Cobj$rmnc(mean, chol_cov, FALSE)}) + expect_equal( + {set.seed(10); obj$rmnc(mean, chol_prec, TRUE)}, + {set.seed(10); Cobj$rmnc(mean, chol_prec, TRUE)}) + + # dmvt + expect_equal(obj$dmvtc(x, mean, chol_cov, 5, FALSE, FALSE), + Cobj$dmvtc(x, mean, chol_cov, 5, FALSE, FALSE)) + expect_equal(obj$dmvtc(x, mean, chol_cov, 5, FALSE, TRUE), + Cobj$dmvtc(x, mean, chol_cov, 5, FALSE, TRUE)) + expect_equal(obj$dmvtc(x, mean, chol_prec, 5, TRUE, FALSE), + Cobj$dmvtc(x, mean, chol_prec, 5, TRUE, FALSE)) + expect_equal(obj$dmvtc(x, mean, chol_prec, 5, TRUE, TRUE), + Cobj$dmvtc(x, mean, chol_prec, 5, TRUE, TRUE)) + expect_equal(Cobj$dmvtc(x, mean, chol_cov, 5, FALSE, TRUE), + Cobj$dmvtc(x, mean, chol_prec, 5, TRUE, TRUE)) + + # rmvt + expect_equal( + {set.seed(10); obj$rmvtc(mean, chol_cov, 5, FALSE)}, + {set.seed(10); Cobj$rmvtc(mean, chol_cov, 5, FALSE)}) + expect_equal( + {set.seed(10); obj$rmnc(mean, chol_prec, TRUE)}, + {set.seed(10); Cobj$rmnc(mean, chol_prec, TRUE)}) + + # dwish + x <- matrix(c(0.9, 0.45, 0.2, 0.45, 0.9, 0.45, 0.2, 0.45, 0.9), nrow = 3, byrow =TRUE) + expect_equal(obj$dwishc(x, chol_cov, 5, FALSE, FALSE), + Cobj$dwishc(x, chol_cov, 5, FALSE, FALSE)) + expect_equal(obj$dwishc(x, chol_cov, 5, FALSE, TRUE), + Cobj$dwishc(x, chol_cov, 5, FALSE, TRUE)) + expect_equal(obj$dwishc(x, chol_prec, 5, TRUE, FALSE), + Cobj$dwishc(x, chol_prec, 5, TRUE, FALSE)) + expect_equal(obj$dwishc(x, chol_prec, 5, TRUE, TRUE), + Cobj$dwishc(x, chol_prec, 5, TRUE, TRUE)) + expect_equal(Cobj$dwishc(x, chol_cov, 5, FALSE, TRUE), + Cobj$dwishc(x, chol_prec, 5, TRUE, TRUE)) + + #rwish + expect_equal( + {set.seed(10); obj$rwishc(chol_cov, 5, FALSE)}, + {set.seed(10); Cobj$rwishc(chol_cov, 5, FALSE)}) + expect_equal( + {set.seed(10); obj$rwishc(chol_prec, 5, TRUE)}, + {set.seed(10); Cobj$rwishc(chol_prec, 5, TRUE)}) + + # ddirch + alpha <- c(3, 4, 2, 6) + x <- c(0.4, 0.1, 0.2, 0.3) + expect_equal(obj$ddirch_(x, alpha, FALSE), + Cobj$ddirch_(x, alpha, FALSE)) + expect_equal(obj$ddirch_(x, alpha, TRUE), + Cobj$ddirch_(x, alpha, TRUE)) + + # rdirch + expect_equal( + {set.seed(10); obj$rdirch_(alpha)}, + {set.seed(10); Cobj$rdirch_(alpha)}) + + # dmulti + x <- c(3, 4, 2, 6) + prob <- c(0.4, 0.1, 0.2, 0.3) + expect_equal(obj$dmulti_(x, prob, FALSE), + Cobj$dmulti_(x, prob, FALSE)) + expect_equal(obj$dmulti_(x, prob, TRUE), + Cobj$dmulti_(x, prob, TRUE)) + + + # rmulti + expect_equal( + {set.seed(10); obj$rmulti_(50, prob)}, + {set.seed(10); Cobj$rmulti_(50, prob)}) + + rm(obj, Cobj); gc() + +})