Update AcceptReject.cpp

DESCRIPTION

@@ -24,6 +24,7 @@ Imports:
     numDeriv,
     pbmcapply,
     purrr,
+    Rcpp,
     rlang,
     scales
 Suggests: 
@@ -33,3 +34,6 @@ Suggests:
     tictoc,
     testthat (>= 3.0.0)
 Config/testthat/edition: 3
+LinkingTo: 
+    Rcpp,
+    RcppArmadillo

NAMESPACE

@@ -5,6 +5,7 @@ S3method(print,accept_reject)
 export(accept_reject)
 export(inspect)
 import(rlang)
+importFrom(Rcpp,evalCpp)
 importFrom(assertthat,assert_that)
 importFrom(cli,cli_alert_danger)
 importFrom(cli,cli_alert_info)
@@ -32,9 +33,7 @@ importFrom(graphics,hist)
 importFrom(lbfgs,lbfgs)
 importFrom(numDeriv,grad)
 importFrom(parallel,detectCores)
-importFrom(parallel,makeCluster)
-importFrom(parallel,stopCluster)
-importFrom(pbmcapply,pbmclapply)
+importFrom(parallel,mclapply)
 importFrom(purrr,map_dbl)
 importFrom(purrr,partial)
 importFrom(rlang,list2)
@@ -46,3 +45,4 @@ importFrom(stats,runif)
 importFrom(utils,capture.output)
 importFrom(utils,head)
 importFrom(utils,tail)
+useDynLib(AcceptReject, .registration = TRUE)

NEWS.md

@@ -4,7 +4,7 @@
 
 # AcceptReject 0.1.1
 
-* Improved performance in serial and parallel processing;
+* Improved performance in serial and parallel processing with Rcpp and RcppArmadillo;
 
 * Now it is possible to specify a different base density/probability mass function than the uniform one. If none is specified, the uniform density (either discrete or continuous) is assumed for the case of discrete or continuous random variables, respectively;
 

R/RcppExports.R

@@ -0,0 +1,7 @@
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+one_step <- function(n, f, f_base, random_base, c) {
+    .Call(`_AcceptReject_one_step`, n, f, f_base, random_base, c)
+}
+

R/accept_reject.r

@@ -106,13 +106,15 @@
 #' @importFrom lbfgs lbfgs
 #' @importFrom purrr partial map_dbl
 #' @importFrom numDeriv grad
-#' @importFrom parallel detectCores makeCluster stopCluster
+#' @importFrom parallel detectCores mclapply
 #' @importFrom stats dunif runif dweibull
 #' @importFrom utils capture.output
 #' @importFrom assertthat assert_that
 #' @importFrom cli cli_alert_danger cli_alert_warning
 #' @importFrom glue glue
-#' @importFrom pbmcapply pbmclapply
+#' @useDynLib AcceptReject, .registration = TRUE
+#' @importFrom Rcpp evalCpp
+#'
 #' @export
 accept_reject <-
   function(n = 1L,
@@ -218,17 +220,6 @@ accept_reject <-
         )$par[1L]
     }
 
-    one_step <- function(n) {
-      x <- numeric(0L)
-      while(length(x) < n) {
-        x_cand <- random_base(n = n)
-        u <- runif(n = n)
-        accepted <- u <= f(x = x_cand) / (c * f_base(x = x_cand))
-        x <- c(x, x_cand[accepted])
-      }
-      return(x[1L:n])
-    }
-
     if (parallel && .Platform$OS.type == "unix") {
       n_cores <- parallel::detectCores()
       n_per_core <- n %/% n_cores
@@ -238,14 +229,13 @@ accept_reject <-
           rep(n_per_core, n_cores - remainder),
           rep(n_per_core + 1L, remainder)
         )
-      capture.output(
-        r <- unlist(pbmcapply::pbmclapply(
-          X = n_each_core,
-          FUN = one_step,
-          mc.cores = n_cores
-          )))
+      r <- unlist(parallel::mclapply(
+        X = n_each_core,
+        FUN = function(n) one_step(n, f, f_base, random_base, c),
+        mc.cores = n_cores
+      ))
     } else {
-      r <- one_step(n)
+      r <- one_step(n, f, f_base, random_base, c)
     }
 
     class(r) <- "accept_reject"

README.md

@@ -289,7 +289,7 @@ case_1 <- accept_reject(
   xlim = c(0, 10)
 )
 toc()
-#> 0.491 sec elapsed
+#> 0.375 sec elapsed
 
 # Specifying the base probability density function
 tic()
@@ -305,7 +305,7 @@ case_2 <- accept_reject(
   c = 1.2
 )
 toc()
-#> 0.153 sec elapsed
+#> 0.145 sec elapsed
 
 # Visualizing the results
 p1 <- plot(case_1)

docs/pkgdown.yml

@@ -4,5 +4,5 @@ pkgdown_sha: ~
 articles:
   accept_reject: accept_reject.html
   inspect: inspect.html
-last_built: 2024-04-23T17:24Z
+last_built: 2024-04-24T13:46Z
 

src/.gitignore

@@ -0,0 +1,3 @@
+*.o
+*.so
+*.dll

src/AcceptReject.cpp

@@ -0,0 +1,52 @@
+// #include <RcppArmadillo.h>
+//
+// // [[Rcpp::depends(RcppArmadillo)]]
+// // [[Rcpp::plugins(cpp11)]]
+//
+// using namespace Rcpp;
+//
+// // [[Rcpp::export]]
+//  arma::vec one_step(int n, Function f, Function f_base, Function random_base, double c) {
+//    arma::vec x;
+//    while (x.size() < n) {
+//      arma::vec x_cand = as<arma::vec>(random_base(n));
+//      arma::vec u = as<arma::vec>(runif(n));
+//      arma::uvec accepted = find(u <= as<arma::vec>(f(wrap(x_cand))) / (c * as<arma::vec>(f_base(wrap(x_cand)))));
+//      x = join_cols(x, x_cand.elem(accepted));
+//    }
+//    return x.rows(0, n-1);
+//  }
+
+#include <RcppArmadillo.h>
+
+// [[Rcpp::depends(RcppArmadillo)]]
+// [[Rcpp::plugins(cpp11)]]
+
+using namespace Rcpp;
+
+// [[Rcpp::export]]
+arma::vec one_step(int n, Function f, Function f_base, Function random_base, double c) {
+  arma::vec x(n);
+  int filled = 0;
+
+  while (filled < n) {
+    arma::vec x_cand = as<arma::vec>(random_base(n));
+    arma::vec fx_cand = as<arma::vec>(f(wrap(x_cand)));
+    arma::vec u = as<arma::vec>(runif(n));
+    arma::uvec accepted = find(u <= fx_cand / (c * as<arma::vec>(f_base(wrap(x_cand)))));
+
+    int num_accepted = accepted.n_elem;
+    if (num_accepted > 0) {
+      int space_left = n - filled;
+      if (num_accepted > space_left) {
+        num_accepted = space_left;
+        accepted = accepted.subvec(0, space_left - 1);
+      }
+
+      x.subvec(filled, filled + num_accepted - 1) = x_cand.elem(accepted);
+      filled += num_accepted;
+    }
+  }
+
+  return x;
+}

src/RcppExports.cpp

@@ -0,0 +1,38 @@
+// Generated by using Rcpp::compileAttributes() -> do not edit by hand
+// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+#include <RcppArmadillo.h>
+#include <Rcpp.h>
+
+using namespace Rcpp;
+
+#ifdef RCPP_USE_GLOBAL_ROSTREAM
+Rcpp::Rostream<true>&  Rcpp::Rcout = Rcpp::Rcpp_cout_get();
+Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
+#endif
+
+// one_step
+arma::vec one_step(int n, Function f, Function f_base, Function random_base, double c);
+RcppExport SEXP _AcceptReject_one_step(SEXP nSEXP, SEXP fSEXP, SEXP f_baseSEXP, SEXP random_baseSEXP, SEXP cSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< int >::type n(nSEXP);
+    Rcpp::traits::input_parameter< Function >::type f(fSEXP);
+    Rcpp::traits::input_parameter< Function >::type f_base(f_baseSEXP);
+    Rcpp::traits::input_parameter< Function >::type random_base(random_baseSEXP);
+    Rcpp::traits::input_parameter< double >::type c(cSEXP);
+    rcpp_result_gen = Rcpp::wrap(one_step(n, f, f_base, random_base, c));
+    return rcpp_result_gen;
+END_RCPP
+}
+
+static const R_CallMethodDef CallEntries[] = {
+    {"_AcceptReject_one_step", (DL_FUNC) &_AcceptReject_one_step, 5},
+    {NULL, NULL, 0}
+};
+
+RcppExport void R_init_AcceptReject(DllInfo *dll) {
+    R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
+    R_useDynamicSymbols(dll, FALSE);
+}