4

在 Rcpp 项目中，我希望能够调用 R 函数（包中的cobs函数进行cobs凹样条拟合）或调用它依赖的 fortran 代码cobs（函数使用rq.fit.sfnc用于拟合约束样条模型的函数，该模型又依赖于包中的 fortran 编码srqfnc函数）在pragma openmp 并行 for 循环中quantregquantreg（我的其余代码主要需要一些简单的线性代数，所以这没有问题，但遗憾的是，每个内部循环迭代也需要我进行凹样条拟合）。我想知道这是否会被允许或可能，因为我认为这样的调用不会是线程安全的？是否有一个简单的解决方法，比如用 a 包围这些电话#pragma omp critical？有人会有这方面的例子吗？或者在这种情况下，唯一的方法是首先使用线程安全的 Armadillo 类对&函数进行完整Rcpp的移植？cobsrq.fit.sfnc

4

2 回答 2

4

引用手册：

从线程代码调用任何 R API 是“仅供专家使用”，强烈建议不要这样做。R API 中的许多函数会修改内部 R 数据结构，如果从多个线程同时调用，可能会损坏这些数据结构。大多数 R API 函数都可以发出错误信号，这只能发生在 R 主线程上。此外，外部库（例如 LAPACK）可能不是线程安全的。

我一直将此解释为“不得从线程代码中调用 R API 函数”。无论内部使用什么，从 omp 并行区域内部调用 R 函数就是这样。使用#pragma omp critical 可能会起作用，但如果它坏了，你必须保留碎片......

重新实现相关代码或在 C++/C/Fortran 中查找现有实现并直接调用它会更安全。

于 2019-03-19T15:34:10.160 回答
3

所以我刚刚尝试过，似乎只有在#pragma openmp parallel for循环中调用 R 函数才有效#pragma omp critical（否则它会导致堆栈不平衡，并使 R 崩溃）。当然，这将导致那部分代码按顺序执行，但这在某些情况下可能仍然有用。

例子：

零件，Rcpp另存为文件"fitMbycol.cpp"：

// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// #define RCPP_ARMADILLO_RETURN_COLVEC_AS_VECTOR
using namespace Rcpp;
using namespace arma;
using namespace std;

#include <omp.h>
// [[Rcpp::plugins(openmp)]]

// [[Rcpp::export]]
arma::mat fitMbycol(arma::mat& M, Rcpp::Function f, const int nthreads) {

  // ARGUMENTS
  // M: matrix for which we want to fit given function f over each column
  // f: fitting function to use with one single argument (vector y) that returns the fitted values as a vector
  // nthreads: number of threads to use

  // we apply fitting function over columns
  int c = M.n_cols;
  int r = M.n_rows;
  arma::mat out(r,c);
  int i;
  omp_set_num_threads(nthreads);
#pragma omp parallel for shared(out)
  for (i = 0; i < c; i++) {
      arma::vec y = M.col(i); // ith column of M
#pragma omp critical
{
      out.col(i) = as<arma::colvec>(f(NumericVector(y.begin(),y.end())));
}
  }

  return out;

}

然后在 R 中：

首先是纯R版本：

（我们用泊松噪声模拟一些高斯峰形状，然后使用cobs函数对它们进行对数凹样条拟合）

x=1:100
n=length(x)
ncols=50
gauspeak=function(x, u, w, h=1) h*exp(((x-u)^2)/(-2*(w^2)))
Y_nonoise=do.call(cbind,lapply(seq(min(x), max(x), length.out=ncols), function (u) gauspeak(x, u=u, w=10, h=u*100)))
set.seed(123)
Y=apply(Y_nonoise, 2, function (col) rpois(n,col))

# log-concave spline fit on each column of matrix Y using cobs
require(cobs)
logconcobs = function(y, tau=0.5, nknots=length(y)/10) {
  x = 1:length(y)
  offs = max(y)*1E-6
  weights = y^(1/2)
  fit.y = suppressWarnings(cobs(x=x,y=log10(y+offs), 
              constraint = "concave", lambda=0, 
              knots = seq(min(x),max(x), length.out = nknots), 
              nknots=nknots, knots.add = FALSE, repeat.delete.add = FALSE,
              keep.data = FALSE, keep.x.ps = TRUE,
              w=weights, 
              tau=tau, print.warn = F, print.mesg = F, rq.tol = 0.1, maxiter = 100)$fitted)
  return(pmax(10^fit.y - offs, 0))
}
library(microbenchmark)
microbenchmark(Y.fitted <- apply(Y, 2, function(col) logconcobs(y=col, tau=0.5)),times=5L) # 363 ms, ie 363/50=7 ms per fit
matplot(Y,type="l",lty=1)
matplot(Y_nonoise,type="l",add=TRUE, lwd=3, col=adjustcolor("blue",alpha.f=0.2),lty=1)
matplot(Y.fitted,type="l",add=TRUE, lwd=3, col=adjustcolor("red",alpha.f=0.2),lty=1)

在此处输入图像描述

现在使用在 a中Rcpp调用我们的 R 拟合函数，用 a 括起来：logconcobs#pragma openmp parallel for loop#pragma omp critical

library(Rcpp)
library(RcppArmadillo)
Rcpp::sourceCpp('fitMbycol.cpp')
microbenchmark(Y.fitted <- fitMbycol(Y, function (y) logconcobs(y, tau=0.5, nknots=10), nthreads=8L ), times=5L) # 361 ms

在这种情况下，OpenMP 当然最终不会产生任何影响，因为这#pragma omp critical会导致所有内容都按顺序执行，但在更复杂的示例中，这仍然很有用。

于 2019-03-19T21:35:20.047 回答