c++ - 如何在RcppParallel中调用用户定义的函数？

受艺术性http://gallery.rcpp.org/articles/parallel-distance-matrix/的启发，我尝试使用RcppParallel在高维参数空间中运行蛮力搜索，以使用多线程进行回测。我被困在如何在struct部分中调用自定义函数。这个想法是这样的:

首先，首先在R中创建一个参数矩阵NumericMatrix params_mat，并使用具有List, NumericVector, CharacterVector数据类型的回测数据，例如List Data_1, NumericVector Data_2, CharacterVector Data_3, ...，它们对于每个参数场景params_vec都是静态的(请注意，这是params_mat的行)。

接下来，定义回测函数，该函数输出一个包含3个关键变量的 vector ，以评估策略绩效。

这是我的params_mat和Backtesting_Fun的示例，可以分别在R和Rcpp中运行。

//[[Rcpp::export]]
NumericMatrix data_frame_rcpp(const Rcpp::List& list_params)
{
  NumericMatrix res = list_params[0];
  return res;
}

# R codes to generate params_mat
params <- expand.grid (x_1=seq(1,100,1), x_2=seq(3,100,2), ..., x_n=seq(4,200,1));
list_params = list(ts(params));
tmp_params_data = data_frame_rcpp(list_params);
params_mat = matrix(tmp_params_data, ncol = ncol(tmp_params_data), dimnames = NULL);
params_vec = params_mat[ii,];

# User-defined Rcpp codes for backtesting
NumericVector Backtesting_Fun (List Data_1, NumericVector Data_2, CharacterVector Data_3, ..., NumericVector params_vec)
{
  // Main function parts to run backtesting for each params_vec scenario.
  ... etc

  // save 3 key result variables together with each params_vec (just a simple illustration).
  NumericVector res = NumericVector::create(params_vec[0],...,params_vec[size-1],
                                            key_1, key_2, key_3);
  return res;
}

当然，我们需要使用RVector/RMatrix类型重写/修改原始的Rcpp Backtesting_Fun，然后使用以下RcppParallel代码在Backtesting_Fun中调用struct Backtest_parallel:

// [[Rcpp::depends(RcppParallel)]]
#include <RcppParallel.h>
using namespace RcppParallel;

RVector<double> Backtesting_Fun (const RVector<double> Data_1, const RVector<double> Data_2,
                                const RVector<string> Data_3,..., const RVector<double> params_vec)
{
   // Main function parts to run backtesting for each params_vec scenario.
   ... etc;

   // save 3 key result variables together with each params_vec
   ... etc;

   return res;
}

struct Backtest_parallel : public Worker
{
   // input matrix to read from
   const RVector<List> Data_1;
   const RVector<double> Data_2;
   const RVector<string> Data_3;
   ...
   const RMatrix<double> params_mat;

   // output matrix to write to
   RMatrix<double> rmat;

   // initialize from Rcpp input and output matrixes (the RMatrix class
   // can be automatically converted to from the Rcpp matrix type)
   Backtest_parallel(const List Data_1, const NumericVector Data_2,
   const CharacterVector Data_3, ..., const NumericMatrix params_mat)
      : Data_1(Data_1), Data_2(Data_2), Data_3(Data_3), ..., params_mat(params_mat) {}

   // function call operator that work for the specified range (begin/end)
   void operator()(std::size_t begin, std::size_t end)
   {
      for (std::size_t ii = begin; ii < end; i++)
      {
         // params rows that we will operate on
         RMatrix<double>::Row params_row = params_mat.row(ii);

         // Run the backtesting function defined above
         RVector<double> res = Backtesting_Fun(Data_1, Data_2, ..., params_row)
         for (std::size_t jj = 0; jj < res.length(); jj++)
         {
            // write to output matrix
            rmat(ii,jj) = res[jj];
         }
      }
   }
};

// [[Rcpp::export]]
NumericMatrix rcpp_parallel_backtest(List Data_1, NumericVector Data_2, CharacterVector Data_3,
                                     ..., NumericMatrix params_mat)
{
   // allocate the matrix we will return
   NumericMatrix rmat(params_mat.nrow(), params_mat.nrow()+3);

   // create the worker
   Backtest_parallel backtest_parallel(Data_1, Date_2, ..., params_mat);

   // call it with parallelFor
   parallelFor(0, rmat.nrow(), backtest_parallel);

   return rmat;
}

这是我的问题:

RVector可以包含List数据类型，还是在RcppParallel中有任何特定的容器包含List？

在Backtesting_Fun中，输入应为RVector/RMatrix类型，这是否意味着我们真的需要将原始的Rcpp主代码和NumericVector转换为RVector吗？

还是在RcppParallel中有更好的方式对我的案例进行并行计算？提前致谢。

编辑:

我看看http://gallery.rcpp.org/articles/parallel-matrix-transform/，http://gallery.rcpp.org/articles/parallel-inner-product/中有关RcppPararrel的其他示例，struct operator()的共同思想是使用指针来操作operator()的数据输入，那么在我的情况下，是否有任何方法可以使用指针输入来构建用户定义的函数？

如果上述方法无效，是否可以使用wrap将RVector/RMatrix转换回Rcpp数据类型，即NumericVector..中的operator()，以便用户定义函数Backtesting_Fun的输入类型保持不变。

最佳答案

我想我可能会找到解决此问题的另一种方法:关键是使用线程安全访问器在struct中包含变量，并在外部main函数中保留RVector/RMatrix，以便parallelFor可以正常工作，这是这个并行算法中最重要的部分。这是我的方法:

摆脱掉List数据类型:相反，我们可以使用List/NumericVector容器转换NumericMatrix变量并记录其对应的索引，以便subvector/submatrix指向与列表元素相同的元素。

将RVector/RMatrix转换为arma::vec/arma::mat :如RcppParallel Github中所述，C++ Armadillo在struct的运算符中是线程安全的。在这里，我使用这种想法修改了Parallel Distance Matrix Calculation with RcppParallel中给出的示例，该示例几乎保持相同的测试速度。

struct JsDistance : public Worker
{
  const RMatrix<double> tmp_MAT;  // input matrix to read from
  RMatrix<double> tmp_rmat;       // output matrix to write to
  std::size_t row_size, col_size;

  // Convert global input/output into RMatrix/RVector type
  JsDistance(const NumericMatrix& matrix_input, NumericMatrix& matrix_output,
             std::size_t row_size, std::size_t col_size)
    : tmp_MAT(matrix_input), tmp_rmat(matrix_output), row_size(row_size), col_size(col_size) {}

  // convert RVector/RMatrix into arma type for Rcpp function
  // and the follwing arma data will be shared in parallel computing
  arma::mat convert()
  {
    RMatrix<double> tmp_mat = tmp_MAT;
    arma::mat MAT(tmp_mat.begin(), row_size, col_size, false);
    return MAT;
  }


  void operator()(std::size_t begin, std::size_t end)
  {
    for (std::size_t i = begin; i < end; i++)
    {
      for (std::size_t j = 0; j < i; j++)
      {
        // rows we will operate on
        arma::mat MAT = convert();
        arma::rowvec row1 = MAT.row(i);          // get the row of arma matrix
        arma::rowvec row2 = MAT.row(j);

        // compute the average using std::tranform from the STL
        std::vector<double> avg(row1.n_elem);
        std::transform(row1.begin(), row1.end(), // input range 1
                       row2.begin(),             // input range 2
                       avg.begin(),              // output range
                       average);                 // function to apply

        // calculate divergences
        double d1 = kl_divergence(row1.begin(), row1.end(), avg.begin());
        double d2 = kl_divergence(row2.begin(), row2.end(), avg.begin());

        // write to output matrix
        tmp_rmat(i,j) = sqrt(.5 * (d1 + d2));
      }
    }
  }
};

// [[Rcpp::export]]
NumericMatrix rcpp_parallel_js_distance_modify(const Rcpp::NumericMatrix& matrix_input, int N_cores)
{
  // allocate the matrix we will return
  NumericMatrix matrix_output(matrix_input.nrow(), matrix_input.nrow());
  std::size_t row_size = matrix_input.nrow();
  std::size_t col_size = matrix_input.ncol();

  // create the worker
  JsDistance jsDistance(matrix_input, matrix_output, row_size, col_size);

  // call it with parallelFor
  parallelFor(0, matrix_input.nrow(), jsDistance, matrix_input.nrow()/N_cores);           // parallelFor with grain size setting

  return matrix_output;
 }

 // Example compare:
 n_row = 1E3;
 n_col = 1E2;
 m = matrix(runif(n_row*n_col), nrow = n_row, ncol = n_col);
 m = m/rowSums(m);

 res <- benchmark(rcpp_parallel_js_distance(m, 6),
         rcpp_parallel_js_distance_orignal(m, 6),
         order="relative")
 res[,1:4];

 #test                                    #elapsed   #relative
 rcpp_parallel_js_distance_orignal(m, 6)  128.069    1.000
 rcpp_parallel_js_distance(m, 6)          129.210    1.009

如我们所见，operator中的数据类型将为C++ arma，现在我们可以通过直接使用对象而不是仅使用指针来安全，快速地调用用户定义的函数，该对象可能不是通用的，也不容易设计。

现在，此parallelFor结构将共享相同的数据源，而无需在并行计算中进行额外的复制，然后，我们可以使用上述问题中提到的想法对回测做一些小的更改。