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

受到文章的启发http://gallery.rcpp.org/articles/parallel-distance-matrix/ http://gallery.rcpp.org/articles/parallel-distance-matrix/，我尝试使用 RcppParallel 在高维参数空间中运行强力搜索，以便使用多线程进行回溯测试。我陷入了如何调用自定义函数的困境struct部分。这个想法是这样的：

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

接下来，定义回测函数，输出由 3 个关键变量组成的向量来评估策略性能。

这是我的一个例子params_mat and 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;
}

当然我们需要重写/修改原来的RcppBacktesting_Fun与 RVector/RMatrix 类型，然后使用以下RcppParallel调用代码Backtesting_Fun in 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;
}

这是我的问题：

1. Can RVector包含List数据类型，或者是否有任何特定的容器RcppParallel包含装有List;

2. In the Backtesting_Fun，输入应该是RVector/RMatrix类型，这是否意味着我们真的需要将原始 Rcpp 主代码转换为NumericVector into RVector?

或者有没有更好的方法可以在 RcppParallel 中为我的情况进行并行计算？提前致谢。

EDIT:

1. 我查看了有关 RcppParallel 的其他示例http://gallery.rcpp.org/articles/parallel-matrix-transform/ http://gallery.rcpp.org/articles/parallel-matrix-transform/, http://gallery.rcpp.org/articles/parallel-inner-product/ http://gallery.rcpp.org/articles/parallel-inner-product/，共同的想法是struct operator()是使用指针来操作数据输入operator()，那么在我的例子中，有没有办法使用指针输入构建用户定义的函数？

2. 如果以上方法不行的话是否可以使用wrap转换RVector/RMatrix返回 Rcpp 数据类型，即NumericVector.. in operator()这样用户定义函数的输入类型Backtesting_Fun可以保持不变。

我想我可能会找到一种替代方法来解决这个问题：关键是使用线程安全访问器来包含变量struct，并保持RVector / RMatrix在外部主函数中，以便parallelFor可以很好地工作，这是这个并行算法中最重要的部分。以下是我的方法：

1. 摆脱List数据类型: 相反，我们可以将List变量通过使用NumericVector / NumericMatrix容器并记录其相应的索引，以便子向量/子矩阵将指向与列表元素相同的元素。

2. Convert RVector / RMatrix into arma::vec / arma::mat: 如中提到的Rcpp并行Github http://rcppcore.github.io/RcppParallel/, C++ Armadillo在结构体的运算符中是线程安全的。这里我修改一下给出的例子使用 RcppParallel 进行并行距离矩阵计算 http://gallery.rcpp.org/articles/parallel-distance-matrix/通过使用这个想法，几乎保持相同的测试速度。

   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结构体在并行计算中会共享相同的数据源，无需额外的副本，然后我们可以利用上个问题中提到的想法做一些轻微的改变来进行回测。