Разрешено ли / возможно ли вызвать функцию R или код Fortran в прагме openmp, параллельной циклу for в Rcpp?

Question

В проекте Rcpp я хотел бы иметь возможность вызвать функцию R (функция cobs из пакета cobs для выполнения вогнутой подгонки сплайна) или вызовите код Fortran , на который он опирается (функция cobs использует функцию quantreg rq.fit.sfnc для соответствия модели с ограниченным сплайном, которая, в свою очередь, опирается на кодированная * fortune функция srqfnc в пакете quantreg) в рамках прагмы openmp параллельной для цикла (остальная часть моего кода в основном требует некоторой простой линейной алгебры, так что это не будет проблемой, но, к сожалению, каждая итерация внутреннего цикла также требует от меня вогнутой подгонки сплайна). Мне было интересно, если это будет разрешено или возможно, хотя, как я полагаю, такие вызовы не будут потокобезопасными? Будет ли это легко исправить, например, окружить эти вызовы #pragma omp critical? Есть ли у кого-нибудь примеры этого? Или единственный способ в этом случае заключается в том, чтобы сначала выполнить полный Rcpp порт функций cobs & rq.fit.sfnc с использованием поточно-ориентированных классов Armadillo?

Answer 1

Цитирование руководство :

Вызов любого из API R из многопоточного кода "только для экспертов" и настоятельно не рекомендуется. Многие функции в R API модифицируют внутренние структуры данных R и могут повредить эти структуры данных, если их вызвать одновременно из нескольких потоков. Большинство функций R API могут сигнализировать об ошибках, которые должны происходить только в главном потоке R. Кроме того, внешние библиотеки (например, LAPACK) могут быть не поточно-ориентированными.

Я всегда интерпретировал это как «один не должен вызывать функцию R API из многопоточного кода». Вызов функции R, независимо от того, что используется внутри, внутри параллельной области omp, был бы именно этим. Использование #pragma omp critical может сработать, но если оно сломается, вы должны сохранить куски ...

Было бы безопаснее либо заново реализовать рассматриваемый код, либо искать существующую реализацию в C ++ / C / Fortran и вызывать ее напрямую.

Answer 2

Итак, я только что попробовал, и кажется, что вызов R-функций в цикле #pragma openmp parallel for работает, только если ему предшествует #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)

А теперь с помощью Rcpp вызывается наша функция подгонки R 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 заставляет все выполняться последовательно, но в более сложномпримеры это еще может быть полезно.