0

我在 Windows 中使用 pystan 模块，其中模块中的 Windows 不支持多线程。pystan 模块部分是用 C++ 编写的，由于我试图减少模块的运行时间，我想知道是否有办法在模块的 C++ 部分中手动编写多线程代码以减少运行时间所以我可以增加迭代次数吗？下面是代码：

from __future__  import division
import pystan
import numpy as np
import os 

x=np.array([453.05,453.05,453.24,453.35,453.44,453.44,453.83,454.02,454.89])
y=np.array([3232.12,3231.45,3231.90,3231.67,3231.84,3231.95,3231.89,3231.67,3231.45])
x=np.array(zip(x,y))
c=np.array([0.01,0.07,0.001,0.1,0.05,0.001,0.001,0.05,0.001])
s = np.array([454.4062631951059,3230.808656891571])
st=np.array([12,12,12,12,12,12,12,12,12])

model='''
    data {
     int D; //number of dimensions
     int K; //number of gaussians
     int N; //number of data 

     vector[D] y[N]; // observation data
     real con[N]; //concentration
     vector[D] s;//oil spill location
     real st[N]; // sample time
    }

    parameters {
     simplex[K] theta; //mixing proportions
     vector[D] v[K];
     vector<lower=0>[D] Dif[K];
     cholesky_factor_corr[D] L[K]; //cholesky factor of correlation matrix
    }

    transformed parameters {
      cholesky_factor_cov[D,D] cov[K,N];
      vector<lower=0>[D] sigma[K,N]; // standard deviations  
      vector[D] mu[K,N];
      real ro[K];
      matrix[D,D] Omega[K];
      matrix[D,D] Sigma[K,N];  
      vector[N] lamba;  

      for (k in 1:K) {  
      Omega[k] = multiply_lower_tri_self_transpose(L[k]);
         for (n in 1:N){
          sigma[k,n] = 0.05 + sqrt(2*st[n]*Dif[k]);
          mu[k,n] = s+v[k]*st[n];
          cov[k,n] = diag_pre_multiply(sigma[k,n],L[k]);
      Sigma[k,n] = quad_form_diag(Omega[k], sigma[k,n]); 
         }
      ro[k]=Omega[k,2,1]; 
      }

      for (i in 1 : N) {lamba[i] = 1/(theta[1]*(1./2./3.1415926/sqrt (Sigma[1,i, 1, 1])/sqrt (Sigma[1,i, 2, 2])/sqrt (1 - ro[1]*ro[1]))*exp (-1./2./(1 - ro[1]*ro[1])*(-(y[i, 1] - mu[1,i, 1])*(y[i, 1] - mu[1,i, 1])/Sigma[1, i,1, 1] - (y[i, 2] - mu[1, i,2])*(y[i, 2] - mu[1, i,2])/Sigma[1,i, 2, 2] + 2.*ro[1]*(y[i, 1] - mu[1,i, 1])*(y[i, 2] - mu[1,i, 2])/sqrt (Sigma[1, i,1, 1])/sqrt (Sigma[1,i, 2, 2]))) + 
           theta[2]*(1./2./3.1415926/sqrt (Sigma[2, i,1, 1])/sqrt (Sigma[2,i, 2, 2])/sqrt (1 - ro[2]*ro[2]))*exp (-1./2./(1 - ro[2]*ro[2])*(-(y[i, 1] - mu[2, i,1])*(y[i, 1] - mu[2, i,1])/Sigma[2, i,1, 1] - (y[i, 2] - mu[2,i, 2])*(y[i, 2] - mu[2, i,2])/Sigma[2,i, 2, 2] + 2.*ro[2]*(y[i, 1] - mu[2, i,1])*(y[i, 2] - mu[2, i,2])/sqrt (Sigma[2, i,1, 1])/sqrt (Sigma[2, i,2, 2]))) +
           theta[3]*(1./2./3.1415926/sqrt (Sigma[3, i,1, 1])/sqrt (Sigma[3,i, 2, 2])/sqrt (1 - ro[3]*ro[3]))*exp (-1./2./(1 - ro[3]*ro[3])*(-(y[i, 1] - mu[3, i,1])*(y[i, 1] - mu[3, i,1])/Sigma[3, i,1, 1] - (y[i, 2] - mu[3,i, 2])*(y[i, 2] - mu[3, i,2])/Sigma[3,i, 2, 2] + 2.*ro[3]*(y[i, 1] - mu[3, i,1])*(y[i, 2] - mu[3, i,2])/sqrt (Sigma[3, i,1, 1])/sqrt (Sigma[3, i,2, 2]))) +
           theta[4]*(1./2./3.1415926/sqrt (Sigma[4, i,1, 1])/sqrt (Sigma[4,i, 2, 2])/sqrt (1 - ro[4]*ro[4]))*exp (-1./2./(1 - ro[4]*ro[4])*(-(y[i, 1] - mu[4, i,1])*(y[i, 1] - mu[4, i,1])/Sigma[4, i,1, 1] - (y[i, 2] - mu[4,i, 2])*(y[i, 2] - mu[4, i,2])/Sigma[4,i, 2, 2] + 2.*ro[4]*(y[i, 1] - mu[4, i,1])*(y[i, 2] - mu[4, i,2])/sqrt (Sigma[4, i,1, 1])/sqrt (Sigma[4, i,2, 2]))));}
    }

    model {
     real ps[K];
     theta ~ dirichlet(rep_vector(2.0, 4));
     for(k in 1:K){
     v[k,1] ~ normal(0.0,4.1);// uniform(340/100,380/100);//
     v[k,2] ~  normal(0.0,4.1);//uniform(3160/100,3190/100);//
     Dif[k] ~ normal(0.5,0.2);//exponential(0.05);//beta(2,5);
     L[k] ~ lkj_corr_cholesky(2);// contain rho 
     con ~ exponential(lamba);
     }

     for (n in 1:N){
     for (k in 1:K){
      ps[k] = log(theta[k])+multi_normal_cholesky_lpdf(y[n] | mu[k,N], cov[k,N]); //increment log probability of the gaussian
     }
     target += log_sum_exp(ps);
     }
       for(i in 1:N){
       target +=   - lamba[i]*con[i]+log(lamba[i]);
      }
    }
    '''

    dat={'D':2,'K':4,'N':9,'y':x,'con':c,'s':s,'st':st}
    fit = pystan.stan(model_code=model,data=dat,iter=1000,warmup=500, chains=1,init_r=0.5)
    print(fit)

我对 C++ 不是很精通，因为我一直在使用 python，而 pystan 模块需要用 C++ 编写代码。我希望有一种方法可以对我的 Windows 上不同内核的迭代次数进行多线程处理。

4

1 回答 1

0

由于 Stan 是它自己的语言，因此您只能实现编译器旨在解析和发出代码的功能，其中不包括对任意 C++ 代码的支持。

Stan 的后端确实通过 MPI 提供了对单链并行化的支持，但正如您正确指出的那样，不幸的是，这目前还没有扩展到 Windows。除了尝试编译您自己的以某种方式利用 MPI 库的后端版本之外，您实际上无法在建模语言本身中做任何事情来解决这个问题。

如果您在非 Windows 平台上，则可以通过 map-reduce 操作开始利用新的并行数学库map_rect。用户指南（参见22. Map-Reduce，第 237-244 页）提供了有关使用该方法的一些详细信息，并附有示例。

请注意，您仍然可以在安装了 v2.18+ Stan 后端的任何平台上使用此方法，但只有在 Linux 和 Mac OS X 上才能使用 MPI。我不确定在 Windows 上推出 MPI 支持的计划是什么，但可能值得关注这个 MPI Parallelism wiki 页面。

更新：可能的实验线程支持

似乎有些人一直在尝试使用 RTools 的预发布版本编译 Stan 后端，其中包括 8.x 版本，g++并且可以启用线程数学库和 MPI。听起来像兔子洞，但如果你想试试这里的红色药丸。

于 2018-12-04T18:29:30.617 回答