从 rcpp 返回一个 R 函数

Returning an R function from rcpp

在 Rcpp 中有没有一种方法可以 return 一个带有一些预计算值的 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 中有效完成此任务的建议都将不胜感激(无论是通过 return 函数;还是通过以其他方式存储中间计算)。

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

// [[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

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

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

// [[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.

添加功能方面