从 rcpp 返回 R 函数

Question

Rcpp 中有没有办法返回一个 R 函数，其中包含一些仅在第一次函数调用时计算的预计算值？考虑以下 R 代码:

1: func_generator<-function(X) {
2:  X_tot<-sum(X)
3:  function(b_vec) { (X_tot*b_vec) }
4: }
5: myfunc<-func_generator(c(3,4,5))
6: myfunc(1:2)
7: myfunc(5:6)
8: myfunc2<-func_generator(c(10,11,12,13))
...

这可以在 Rcpp 中编程吗？在实践中，假设在第 2 行的位置上做了一些计算量更大的事情。

为了添加上下文，给定向量 X 和标量 b，有一些似然函数 f(b|X)，对于一些足够的统计量 s(X)，它可以重新表示为 f(b,s(X))，它只是一个函数X 的，这涉及一些计算。这是在计算密集型计算机实验中，有许多向量 X(许多可能性)，以及对每个可能性对 f(bvec|X) 的许多单独调用，所以我宁愿计算 s(X) 一次(对于每个可能性)和以某种方式保存它，而不是多次重新计算。我已经开始通过简单地编程 f(bvec,X) 来计算 bvec=(b_1,...,b_n) 点处的 f(b|X)，但是这有额外的开销，因为我多次调用这个函数并且它在每次运行时计算 s(X)。我只想计算 s(X) 一次。

任何在 Rcpp 中有效完成此任务的建议将不胜感激(无论是通过返回函数；还是通过以其他方式存储中间计算)。

Answer 1

存储中间结果的一种简单方法是函数级别的静态变量:

// [[Rcpp::plugins(cpp11)]]
#include <thread>
#include <chrono>
#include <Rcpp.h>

// [[Rcpp::export]]
Rcpp::NumericVector foo(Rcpp::NumericVector X, Rcpp::NumericVector b, bool useCache = true) {
  static double cache;
  static bool initialized{false};
  if (!(useCache && initialized)) {
    // sleep to simulate actual work
    std::this_thread::sleep_for (std::chrono::seconds(1));
    cache = Rcpp::sum(X);
    initialized = true;
  }
  return cache * b;
}

/*** R
X <- 1:10
b <- 10:20

system.time(r1 <- foo(X, b))
system.time(r2 <- foo(X, b))
all.equal(r1, r2)
system.time(r3 <- foo(X, b, FALSE))
all.equal(r1, r3)
*/

输出:

> system.time(r1 <- foo(X, b))
   user  system elapsed 
      0       0       1 

> system.time(r2 <- foo(X, b))
   user  system elapsed 
  0.002   0.000   0.002 

> all.equal(r1, r2)
[1] TRUE

> system.time(r3 <- foo(X, b, FALSE))
   user  system elapsed 
      0       0       1 

> all.equal(r1, r3)
[1] TRUE

当在第二个函数调用中使用缓存时，结果几乎是即时计算的。

如果您可以遍历不同的 b，这种方法是有效的。在不同的 X 的循环中.如果此限制对您不起作用，那么您可以使用 memoise R 级别的包以有效地存储任意输入的昂贵函数的输出:

// [[Rcpp::plugins(cpp11)]]
#include <thread>
#include <chrono>
#include <Rcpp.h>

// [[Rcpp::export]]
Rcpp::NumericVector foo(double total, Rcpp::NumericVector b) {
  return total * b;
}

// [[Rcpp::export]]
double bar(Rcpp::NumericVector X) {
  // sleep to simulate actual work
  std::this_thread::sleep_for (std::chrono::seconds(1));
  return Rcpp::sum(X);
}


/*** R
X1 <- 1:10
b1 <- 10:20
X2 <- 10:1
b2 <- 20:10

library(memoise)
bar2 <- memoise(bar)

system.time(r11 <- foo(bar2(X1), b1))
system.time(r21 <- foo(bar2(X2), b2))
system.time(r12 <- foo(bar2(X1), b1))
system.time(r22 <- foo(bar2(X2), b2))
all.equal(r11, r12)
all.equal(r21, r22)
*/

输出:

> system.time(r11 <- foo(bar2(X1), b1))
   user  system elapsed 
  0.001   0.000   1.001 

> system.time(r21 <- foo(bar2(X2), b2))
   user  system elapsed 
  0.033   0.000   1.033 

> system.time(r12 <- foo(bar2(X1), b1))
   user  system elapsed 
      0       0       0 

> system.time(r22 <- foo(bar2(X2), b2))
   user  system elapsed 
      0       0       0 

> all.equal(r11, r12)
[1] TRUE

> all.equal(r21, r22)
[1] TRUE

作为替代方案，您还可以使用这两个函数作为函数生成器的构建块:

func_generator <- function(X) {
  X_tot <- bar(X)
  function(b_vec) { foo(X_tot, b_vec) }
}
myfunc <- func_generator(c(3,4,5))
myfunc2 <- func_generator(c(10,11,12,13))
myfunc(1:2)
myfunc(5:6)
myfunc2(1:2)
myfunc2(5:6)

因此，在 C++ 中保留数值昂贵的工作，但保持简单。然后可以使用 R 添加功能方面。