1

我一直在玩 RcppParallel 并编写了一个相当简单的示例来弄清楚事情是如何工作的。代码如下所示。

函数float pdf(double x, double sigma)计算均值为 0 和标准差为 sigma 的高斯分布的缩放版本。

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.0591801 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
4

0 回答 0