rcpp — RcppParallel и C ++. Непоследовательные результаты

Я играл с RcppParallel и написал довольно простой пример, чтобы выяснить, как все работает. Код отображается ниже.

Функция float pdf (двойная х, двойная сигма) вычисляет масштабированную версию гауссовского распределения со средним значением 0 и сигмой стандартного отклонения.

Struct_1 — это структура, которая создает работника для выполнения некоторых вычислений. Я заполняю матрицу, чтобы выяснить, почему некоторые вещи не работают правильно.

void Struct_check () выполняет расчеты.

Функция, кажется, работает, но время от времени она не работает, как ожидалось. Я думаю, что это связано с типами, используемыми для выполнения вычислений в функции pdf!

Пример выполнения отображается под кодом.

Буду признателен за любую помощь помощь!

#include <RcppParallel.h>
#include <RcppArmadillo.h>
#include <RcppArmadilloExtensions/sample.h>
#include <math.h>

#define pi           3.14159265358979323846  /* pi */
using namespace arma;
using namespace Rcpp;
using namespace R;
using namespace sugar;
using namespace std;
using namespace RcppParallel;

// Enable C++11 via this plugin (Rcpp 0.10.3 or later)
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(RcppParallel)]]// Returns the probability of x, given the distribution described by mu and sigma.
float pdf(double x,  double sigma)
{
return exp( -1 * x * x / (2 * sigma * sigma)) / sigma;
}

struct Struct_1 : public Worker
{
arma::vec wr;
arma::vec sr;
NumericVector w2;

// source matrix
const RVector<double> input;

// destination matrix
RMatrix<double> output;

// initialize with source and destination
Struct_1(const NumericMatrix input, NumericMatrix output)
: input(input), output(output) {}

//what is done.
void operator()(std::size_t begin, std::size_t end) {

for (std::size_t i=begin; i<end; i++){ //the processor loop!

NumericVector w2(3);

for (int comp_j=0; comp_j<3; ++comp_j){
w2(comp_j) = wr(comp_j) * pdf( input[i], sr(comp_j) ) ;
}

double sw1 = sum(w2);

output(i,0) = w2(0);
output(i,1) = w2(1);
output(i,2) = w2(2);
output(i,3) = sw1;

w2 = w2/sw1;

output(i,4) = w2(0);
output(i,5) = w2(1);
output(i,6) = w2(2);

double sw2 = sum(w2);

output(i,7) = sw2;

}//end of i loop
}//end of operator
};// [[Rcpp::depends("RcppArmadillo")]]
// [[Rcpp::export]]
void Struct_check(){

//Some vecs defined
arma::vec wr = {0.2522, 0.58523, 0.16257};
arma::vec s2r = {1.2131, 2.9955, 7.5458};
arma::vec sr = sqrt(s2r);

//an arma mat that will be used in the struct
arma::mat arb_mat;
arb_mat.randn(20);
Rcout<<"Arb_mat=\n"<<arb_mat<<endl;

NumericMatrix r_i_x_NM = as<NumericMatrix>(wrap( arb_mat )); //convert to NumericMatrix
NumericMatrix output( r_i_x_NM.nrow() , 8 ); //define the output matrix

Struct_1 struct_1( r_i_x_NM , output);
struct_1.wr = wr;
struct_1.sr = sr;

Rcout<<"nrow output = "<<output.nrow()<<endl;
Rcout<<"ncol output = "<<output.ncol()<<endl;

parallelFor(0, r_i_x_NM.length(), struct_1);
Rcout<<"completed Parallell calculations"<<endl;
Rcout<<"output = \n"<<output<<endl;

}

Запустите изнутри Rstudio. Я использую OS X El Capitan, если это имеет значение.

Struct_check()
Arb_mat=
-0.4539
0.7915
0.2581
1.5917
0.3718
0.4452
0.1230
-1.4719
0.0024
2.6166
-0.4839
-1.2865
2.0492
-1.5980
-0.7531
-0.7312
-1.4482
0.0202
0.4434
-0.0224

nrow output = 20
ncol output = 8
completed Parallell calculations
output =
0.210336 0.326704 0.0583792 0.595419 0.353256 0.548696 0.0980473 1.00000
0.176872 0.304564 0.0567753 0.538211 0.328629 0.565882 0.105489 1.00000
0.222778 0.334398 0.0589211 0.616097 0.361596 0.542768 0.0956361 1.00000
0.0805904 0.221529 0.0500356 0.352155 0.228849 0.629067 0.142084 1.00000
0.216296 0.330423 0.0586421 0.605361 0.357301 0.545827 0.0968712 1.00000
0.211018 0.327133 0.0584096 0.596561 0.353724 0.548365 0.0979106 1.00000
0.227556 0.337284 0.0591224 0.623962 0.364695 0.540551 0.0947533 1.00000
0.0937487 0.235521 0.0512670 0.380537 0.246359 0.618918 0.134723 1.00000
0.228979 0.338136 0.0591817 0.626297 0.365608 0.539897 0.0944947 1.00000
0.0136216 0.107837 0.0375975 0.159056 0.0856401 0.677981 0.236379 1.00000
0.207911 0.325174 0.0582705 0.591355 0.351584 0.549879 0.0985372 1.00000
0.115751 0.256513 0.0530344 0.425298 0.272164 0.603137 0.124699 1.00000
0.0405607 0.167755 0.0448066 0.253123 0.160241 0.662743 0.177015 1.00000
0.0799309 0.220793 0.0499695 0.350694 0.227922 0.629590 0.142488 1.00000
0.181248 0.307594 0.0569989 0.545841 0.332053 0.563523 0.104424 1.00000
0.183689 0.309265 0.0571216 0.550075 0.333934 0.562222 0.103843 1.00000
**0.228941 0.338113 0.0591801 0.618557 0.591026 0.872861 0.152777 1.61666**
0.228941 0.338113 0.0591801 0.626234 0.365583 0.539915 0.0945016 1.61666
0.211153 0.327218 0.0584156 0.596786 0.353816 0.548300 0.0978837 1.00000
0.228932 0.338108 0.0591798 0.626220 0.365578 0.539919 0.0945032 1.00000

Ошибка возникает, когда -1,4482 оценивается для получения следующей строки 0,228941 0,338113 0,051801 0,618557 0,591026 0,872861 0,152777 1,61666

В R — проверке я получаю:

wr <- c(0.2522, 0.58523, 0.16257)
s2r <- c(1.2131, 2.9955, 7.5458)
sr <- sqrt(s2r)

w<-NULL
for (i in 1:3){
w[i] = wr[i]*exp( -0.5*((-1.4482/sr[i])^ 2))/(sr[i])
}

w
[1] 0.09646706 0.23826346 0.05150315

sum(w)
[1] 0.3862337

w = w/sum(w)

w
[1] 0.2497635 0.6168894 0.1333471