Возврат функции R из rcpp

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), которая может быть выражена как f (b, s (X)) для некоторой достаточной статистики s (X), которая функция только X, и которая включает в себя некоторые вычисления. Это компьютерный эксперимент с интенсивными вычислениями, со многими векторами X (много вероятностей) и многими отдельными вызовами f (bvec | X) для каждой вероятности, поэтому я бы предпочел вычислить s (X) один раз (для каждой вероятности) и сохранить его каким-либо образом, а не пересчитывать его много раз. Я начал с простого программирования f (bvec, X) для оценки f (b | X) в точках bvec = (b_1, ..., b_n), но это требует дополнительных затрат, так как я вызываю эту функцию несколько раз и вычисляет 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.