Многопоточность в C++ Windows

Question

Я использую модуль pystan в Windows, где многопоточность не поддерживается в 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.

Answer 1

Поскольку Stan - это собственный язык, вы можете достичь только того, для чего предназначен компилятор, предназначенный для анализа и генерации кода, что не включает поддержку произвольного кода C++.

Серверная часть Stan обеспечивает поддержку распараллеливания отдельных цепочек через MPI, но, как вы правильно заметили, в настоящее время это не распространяется на Windows. Кроме попыток компилировать ваши собственные версии бэкэндов, которые каким-то образом используют библиотеки MPI, на самом деле вы ничего не можете сделать на языке моделирования как таковом, чтобы исправить это.

Если вы работаете не на платформе Windows, вы можете начать использовать преимущества новой параллельной математической библиотеки с помощью операции уменьшения карты map_rect, Руководство пользователя (см. 22. Map-Reduce, стр. 237-244) содержит некоторые подробности использования этого метода с примерами.

Обратите внимание, что вы все еще можете использовать этот метод на любой платформе с установленным бэкэндом v2.18+ Stan, но MPI будет использоваться только в Linux и Mac OS X. Я не уверен, что планируется развернуть поддержку MPI в Windows, но, возможно, стоит присмотреться к этой вики-странице MPI Parallelism.

---

Обновление: возможна экспериментальная поддержка потоков

Похоже, кто-то играл с компиляцией бэкэнда Stan с предварительной версией RTools, которая включает в себя версию 8.x g++ и это может включить многопоточные математические библиотеки и MPI. Звучит как кроличья нора, но если вы хотите попробовать красную таблетку, вот она.