OpenMP 在 SEIR 模型的 Rcpp 代码中生成段错误
Posted
技术标签:
【中文标题】OpenMP 在 SEIR 模型的 Rcpp 代码中生成段错误【英文标题】：OpenMP generate segfault in Rcpp code for the SEIR model 【发布时间】：2013-10-16 21:39:44 【问题描述】：我编写了一个（可能效率低下，但无论如何..）Rcpp 代码，使用内联来模拟随机SEIR model。
串行版本可以完美地编译和工作，但是因为我需要从中模拟很多次，而且在我看来这似乎是一个令人尴尬的并行问题（只需要再次模拟其他参数值并返回一个带有结果的矩阵）我尝试添加#pragma omp parallel for 并使用-fopenmp -lgomp 进行编译，但是......繁荣！
即使是非常小的例子，我也会遇到段错误！
我尝试添加 setenv("OMP_STACKSIZE","24M",1); 并且值远远超过 24M，但仍然发生段错误。
我将简要解释一下代码，因为它有点长（我试图缩短它，但结果发生了变化，我无法重现它..）： 我有两个嵌套循环，内部循环执行给定参数集的模型，外部循环更改参数。
可能发生竞态条件的唯一原因是代码试图在 inner 循环内并行执行一组指令（由于模型结构，在迭代 @987654326 时无法完成@ 它取决于迭代 t-1) 而不是并行化 outer，但如果我没记错的话，如果将 parallel for 构造函数放在外部...
这基本上是我尝试运行的代码形式：
mat result(n_param,T_MAX);
#pragma omp parallel for
for(int i=0,i<n_param_set;i++)
t=0;
rowvec jnk(T_MAX);
while(t < T_MAX)
...
jnk(t) = something(jnk(t-1));
...
t++;
result.row(i)=jnk;
return wrap(result);
我的问题是：我如何告诉编译器我只想并行计算外部循环（甚至像每个线程的 n_loops/n_threads 一样静态地分配它们）而不是内部循环（实际上是不可并行化的） ?
真正的代码有点复杂，如果你真的愿意，我会在这里展示它以方便重现，但我只是询问 OpenMP 的行为。请注意，唯一的 OpenMP 指令出现在 122 行。
library(Rcpp);library(RcppArmadillo);library(inline)
misc='
#include <math.h>
#define _USE_MATH_DEFINES
#include <omp.h>
using namespace arma;
template <typename T> int sgn(T val)
return (T(0) < val) - (val < T(0));
uvec rmultinomial(int n,vec prob)
int K = prob.n_elem;
uvec rN = zeros<uvec>(K);
double p_tot = sum(prob);
double pp;
for(int k = 0; k < K-1; k++)
if(prob(k)>0)
pp = prob[k] / p_tot;
rN(k) = ((pp < 1.) ? (rbinom(1,(double) n, pp))(0) : n);
n -= rN[k];
else
rN[k] = 0;
if(n <= 0) /* we have all*/
return rN;
p_tot -= prob[k]; /* i.e. = sum(prob[(k+1):K]) */
rN[K-1] = n;
return rN;
'
model_and_summary='
mat SEIR_sim_plus_summaries()
vec alpha;
alpha << 0.002 << 0.0045;
vec beta;
beta << 0.01 << 0.01;
vec gamma;
gamma << 1.0/14.0 << 1.0/14.0;
vec sigma;
sigma << 1.0/(3.5) << 1.0/(3.5);
vec phi;
phi << 0.8 << 0.8;
int S_0 = 800;
int E_0 = 100;
int I_0 = 100;
int R_0 = 0;
int pop = 1000;
double tau = 0.01;
double t_0 = 0;
vec obs_time;
obs_time << 1 << 2 << 3 << 4 << 5 << 6 << 7 << 8 << 9 << 10 << 11 << 12 << 13 << 14 << 15 << 16 << 17 << 18 << 19 << 20 << 21 << 22 << 23 << 24;
const int n_obs = obs_time.n_elem;
const int n_part = alpha.n_elem;
mat stat(n_part,6);
//#pragma omp parallel for
for(int k=0;k<n_part;k++)
ivec INC_i(n_obs);
ivec INC_o(n_obs);
// Event variables
double alpha_t;
int nX; //current number of people moving
vec rates(8);
uvec trans(4); // current transitions, e.g. from S to E,I,R,Universe
vec r(4); // rates e.g. from S to E, I, R, Univ.
/*********************** Initialize **********************/
int S_curr = S_0;
int S_prev = S_0;
int E_curr = E_0;
int E_prev = E_0;
int I_curr = I_0;
int I_prev = I_0;
int R_curr = R_0;
int R_prev = R_0;
int IncI_curr = 0;
int IncI_prev = 0;
int IncO_curr = 0;
int IncO_prev = 0;
double t_curr = t_0;
int t_idx =0;
while( t_idx < n_obs )
// next time preparation
t_curr += tau;
S_prev = S_curr;
E_prev = E_curr;
I_prev = I_curr;
R_prev = R_curr;
IncI_prev = IncI_curr;
IncO_prev = IncO_curr;
/*********************** description (rates) of the events **********************/
alpha_t = alpha(k)*(1+phi(k)*sin(2*M_PI*(t_curr+0)/52)); //real contact rate, time expressed in weeks
rates(0) = (alpha_t * ((double)I_curr / (double)pop ) * ((double)S_curr)); //e+1, s-1, r,i one s get infected (goes in E, not yey infectous)
rates(1) = (sigma(k) * E_curr); //e-1, i+1, r,s one exposed become infectous (goes in I) INCIDENCE!!
rates(2) = (gamma(k) * I_curr); //i-1, s,e, r+1 one i recover
rates(3) = (beta(k) * I_curr); //i-1, s, r,e one i dies
rates(4) = (beta(k) * R_curr); //i,e, s, r-1 one r dies
rates(5) = (beta(k) * E_curr); //e-1, s, r,i one e dies
rates(6) = (beta(k) * S_curr); //s-1 e, i ,r one s dies
rates(7) = (beta(k) * pop); //s+1 one susc is born
// Let the events occour
/*********************** S compartement **********************/
if((rates(0)+rates(6))>0)
nX = rbinom(1,S_prev,1-exp(-(rates(0)+rates(6))*tau))(0);
r(0) = rates(0)/(rates(0)+rates(6)); r(1) = 0.0; r(2) = 0; r(3) = rates(6)/(rates(0)+rates(6));
trans = rmultinomial(nX, r);
S_curr -= nX;
E_curr += trans(0);
I_curr += trans(1);
R_curr += trans(2);
//trans(3) contains dead individual, who disappear...we could avoid this using sequential conditional binomial
/*********************** E compartement **********************/
if((rates(1)+rates(5))>0)
nX = rbinom(1,E_prev,1-exp(-(rates(1)+rates(5))*tau))(0);
r(0) = 0.0; r(1) = rates(1)/(rates(1)+rates(5)); r(2) = 0.0; r(3) = rates(5)/(rates(1)+rates(5));
trans = rmultinomial(nX, r);
S_curr += trans(0);
E_curr -= nX;
I_curr += trans(1);
R_curr += trans(2);
IncI_curr += trans(1);
/*********************** I compartement **********************/
if((rates(2)+rates(3))>0)
nX = rbinom(1,I_prev,1-exp(-(rates(2)+rates(3))*tau))(0);
r(0) = 0.0; r(1) = 0.0; r(2) = rates(2)/(rates(2)+rates(3)); r(3) = rates(3)/(rates(2)+rates(3));
trans = rmultinomial(nX, r);
S_curr += trans(0);
E_curr += trans(1);
I_curr -= nX;
R_curr += trans(2);
IncO_curr += trans(2);
/*********************** R compartement **********************/
if(rates(4)>0)
nX = rbinom(1,R_prev,1-exp(-rates(4)*tau))(0);
r(0) = 0.0; r(1) = 0.0; r(2) = 0.0; r(3) = rates(4)/rates(4);
trans = rmultinomial(nX, r);
S_curr += trans(0);
E_curr += trans(1);
I_curr += trans(2);
R_curr -= nX;
/*********************** Universe **********************/
S_curr += pop - (S_curr+E_curr+I_curr+R_curr); //it should be poisson, but since the pop is fixed...
/*********************** Save & Continue **********************/
// Check if the time is interesting for us
if(t_curr > obs_time[t_idx])
INC_i(t_idx) = IncI_curr;
INC_o(t_idx) = IncO_curr;
IncI_curr = IncI_prev = 0;
IncO_curr = IncO_prev = 0;
t_idx++;
//else just go on...
/*********************** Finished - Starting w/ stats **********************/
// INC_i is the useful variable, how can I change its reference withour copying it?
ivec incidence = INC_i; //just so if I want to use INC_o i have to change just this...
//Scan the epidemics to recover the summary stats (naively divide the data each 52 weeks)
double n_years = ceil((double)obs_time(n_obs-1)/52.0);
vec mu_attack(n_years);
vec ratio_attack(n_years-1);
vec peak(n_years);
vec atk(52);
peak(0)=0.0;
vec tmpExplo(52); //explosiveness
vec explo(n_years);
int year=0;
int week;
for(week=0 ; week<n_obs ; week++)
if(week - 52*year > 51)
mu_attack(year) = sum( atk )/(double)pop;
if(year>0)
ratio_attack(year-1) = mu_attack(year)/mu_attack(year-1);
for(int i=0;i<52;i++)
if(atk(i)>(peak(year)/2.0))
tmpExplo(i) = 1.0;
else
tmpExplo(i) = 0.0;
explo(year) = sum(tmpExplo);
year++;
peak(year)=0.0;
atk(week-52*year) = incidence(week);
if( peak(year) < incidence(week) )
peak(year)=incidence(week);
if(week - 52*year > 51)
mu_attack(year) = sum( atk )/(double)pop;
else
ivec idx(52);
for(int i=0;i<52;i++)
idx(i) = i; //take just the updated ones...
vec tmp = atk.elem(find(idx<(week - 52*year)));
mu_attack(year) = sum( tmp )/((double)pop * (tmp.n_elem/52.0));
ratio_attack(year-1) = mu_attack(year)/mu_attack(year-1);
for(int i=0;i<tmp.n_elem;i++)
if(tmp(i)>(peak(year)/2.0))
tmpExplo(i) = 1.0;
else
tmpExplo(i) = 0.0;
for(int i=tmp.n_elem;i<52;i++)
tmpExplo(i) = 0.0; //to reset the others
explo(year) = sum(tmpExplo);
double correlation2;
double correlation4;
vec autocorr = acf(peak);
/***** ACF *****/
if(n_years<3)
correlation2=0.0;
correlation4=0.0;
else
if(n_years<5)
correlation2 = autocorr(1);
correlation4 = 0.0;
else
correlation2 = autocorr(1);
correlation4 = autocorr(3);
rowvec jnk(6);
jnk << sum(mu_attack)/(year+1.0)
<< (sum( log(ratio_attack)%log(ratio_attack) )/(n_years-1)) - (pow(sum( log(ratio_attack) )/(n_years-1),2))
<< correlation2 << correlation4 << max(peak) << sum(explo)/n_years;
stat.row(k) = jnk;
return stat;
'
main='
std::cout << "max_num_threads " << omp_get_max_threads() << std::endl;
RNGScope scope;
mat summaries = SEIR_sim_plus_summaries();
return wrap(summaries);
'
plug = getPlugin("RcppArmadillo")
## modify the plugin for Rcpp to support OpenMP
plug$env$PKG_CXXFLAGS <- paste('-fopenmp', plug$env$PKG_CXXFLAGS)
plug$env$PKG_LIBS <- paste('-fopenmp -lgomp', plug$env$PKG_LIBS)
SEIR_sim_summary = cxxfunction(sig=signature(),main,settings=plug,inc = paste(misc,model_and_summary),verbose=TRUE)
SEIR_sim_summary()
感谢您的帮助！
注意：在你问之前，我稍微修改了 Rcpp 多项式采样函数，只是因为我更喜欢这种方式而不是使用指针的方式......没有任何其他特殊原因！ :)
【问题讨论】：
尝试在 我看到对rbinom和rmultinomial等函数的调用；请参阅***.com/questions/18967763/…，了解为什么这不能按您预期的那样工作。
如果您从 OpenMP 部分回调到 R，如果您的代码爆炸，请不要感到惊讶。 R 不是多线程或可重入的。
@KevinUshey、rbinom 和 rmultinomial 在串行版本中工作，在并行版本中我没有机会看到它们的行为是否奇怪.. @DirkEddelbuettel 这不是我的代码可以。在（并行）计算结果后，主线程应该返回 R 一个矩阵。
请注意，OMP_STACKSIZE 仅控制其他 OpenMP 线程的堆栈大小。它不影响主线程的堆栈大小。如果您怀疑堆栈是问题所在并且您在类 Unix 操作系统上运行，也许您应该使用 ulimit -s ...。
【参考方案1】：
R 中的核心伪随机数生成器 (PRNG) 并非设计用于多线程环境。也就是说，它们的状态存储在一个静态数组中（dummy 来自src/main/PRNG.c），因此在所有线程之间共享。此外，其他几个静态结构用于存储与核心 PRNG 的高级接口的状态。
一种可能的解决方案是，您将每次调用都放在rnorm() 或其他具有相同名称的命名关键部分内的采样函数，例如：
...
#pragma omp critical(random)
rN(k) = ((pp < 1.) ? (rbinom(1,(double) n, pp))(0) : n);
...
if((rates(0)+rates(6))>0)
#pragma omp critical(random)
nX = rbinom(1,S_prev,1-exp(-(rates(0)+rates(6))*tau))(0);
...
请注意，critical 构造在其后面的结构化块上运行，因此锁定了整个语句。如果在对耗时函数的调用中内联绘制随机数，例如
#pragma omp critical(random)
x = slow_computation(rbinom(...));
这最好转化为：
#pragma omp critical(random)
rb = rbinom(...);
x = slow_computation(rb);
这样只会保护rb = rbinom(...); 语句。
【讨论】：
再次感谢，我只想向未来的观众指出，即使这个解决方案最终有效，我也切换到 boost::random 并遵循这个解决方案使其有效！ ***.com/questions/15504264/boost-random-and-openmp干杯，M以上是关于OpenMP 在 SEIR 模型的 Rcpp 代码中生成段错误的主要内容，如果未能解决你的问题，请参考以下文章
R语言应用实战-基于R浅谈SEIR传染病模型以以及马尔萨斯，logistic模型(推导过程和源代码)
Python数学建模SEIR传染病模型模型延伸-SEIDR模型，加入疫苗接种政府管控病毒变异等因素的影响