用于生成 r x c 列联表的算法的 Rcpp 包装器

Question

我正在尝试使用算法的这个 C++ 实现来有效地生成具有固定边距的 r x c 列联表：https://people.sc.fsu.edu/%7Ejburkardt/cpp_src/asa159/asa159.html

我正在使用 Rcpp 并尝试创建一个可以在 R 中使用的函数。我知道在 R 中有一个 C 实现的这个算法。这对我不起作用，因为我需要能够在 Rcpp 中使用此功能。作为一个中间步骤，我只是想包装这个算法并将其导出到 R。看起来很简单，但我一直在努力解决这个问题。

原函数定义如下：

void rcont2 ( int nrow, int ncol, int nrowt[], int ncolt[], bool *key, 
  int *seed, int matrix[],  int *ierror )

这是一个link直接算法：https://people.sc.fsu.edu/%7Ejburkardt/cpp_src/asa159/asa159.cpp

我不需要错误处理，也不知道如何处理 void 类型，所以我稍微修改了算法，结果是这样的：

int* rcont2 ( int nrow, int ncol, int nrowt[], int ncolt[], bool *key, 
  int *seed, int matrix[])

我做了几次尝试。这是一个例子。这将编译，但是当我尝试运行 R 中的函数时，它崩溃了。

#include <Rcpp.h>
#include "asa159.h"

using namespace Rcpp;

//[[Rcpp::export]]

IntegerVector rcont2_cpp(int nrow,
                         int ncol,
                         IntegerVector nrowt_r,
                         IntegerVector ncolt_r,
                         bool key_r,
                         int seed_r) {

  std::vector<int>  nrowt_rr = as< std::vector<int> >(nrowt_r);
  std::vector<int>  ncolt_rr = as< std::vector<int> >(ncolt_r);

  int* nrowt = &nrowt_rr[0];
  int* ncolt = &ncolt_rr[0];

  int matrix[nrow * ncol];

  bool *key = &key_r;
  int *seed = &seed_r;
 
  int* out = rcont2(nrow, ncol, nrowt, ncolt, key, seed, matrix);

  int len_out = *(&out + 1) - out;

  std::vector<int> value_vec(out, out + len_out); 

  return wrap(value_vec);

  }

我想我在这里犯了一些基本错误。提前谢谢有时间看的人。

Answer 1

这实际上是一个很好的问题，我也很幸运地遇到了来自 FSU 的 John Burkardt 的 very comprehensive site 例程。

但没有什么是完全免费的，需要了解 一点 C 和 C++ 的工作原理 才能通过 Rcpp 将此类例程集成到 R 中。我的最终文件版本比 短得多 ：我们只是将行总和和列总和发送给它，然后得到所需的矩阵。我们还在底部添加了示例的 R 调用（一个不错的技巧！）

代码

#include <Rcpp.h>
#include <asa159.cpp>

// [[Rcpp::export]]
Rcpp::IntegerMatrix rcont2_cpp(Rcpp::IntegerVector rowsums, Rcpp::IntegerVector colsums) {
    int nrow = rowsums.length();
    int ncol = colsums.length();
    Rcpp::IntegerMatrix matrix(nrow, ncol);
    rcont2(nrow, ncol, rowsums.begin(), colsums.begin(), matrix.begin());
    return matrix;
}


/*** R
rs <- c(6, 5)
cs <- c(3, 4, 4)
rcont2_cpp(rs, cs)
rcont2_cpp(rs, cs)              // different answer as randomized algo
rcont2_cpp(rs, cs)              // different answer as randomized algo
*/

输出

这包含三个示例调用。

> sourceCpp("rcont2.cpp")
> rs <- c(6, 5)
> cs <- c(3, 4, 4)
> rcont2_cpp(rs, cs)
     [,1] [,2] [,3]
[1,]    3    1    2
[2,]    0    3    2
> rcont2_cpp(rs, cs)              
     [,1] [,2] [,3]
[1,]    1    3    2
[2,]    2    1    2
> rcont2_cpp(rs, cs)              
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    2    2    1
>

与原始版本不同 asa159.cpp

与原来asa159.cpp的变化很小：我们使用了R自带的RNG，我们简化了调用接口。正如他们所说，通过 Rcpp::stop() 传播错误条件留作练习:)

--- asa159.orig.cpp     2020-01-30 08:20:42.000000000 -0600
+++ asa159.cpp  2022-02-06 19:52:20.999042887 -0600
@@ -428,8 +428,8 @@
 }
 //****************************************************************************80
 
-void rcont2 ( int nrow, int ncol, int nrowt[], int ncolt[], bool *key, 
-  int *seed, int matrix[],  int *ierror )
+void rcont2 ( int nrow, int ncol, int nrowt[], int ncolt[], /*bool *key, */
+              /*int *seed,*/ int matrix[] /*,  int *ierror*/ )
 
 //****************************************************************************80
 //
@@ -517,6 +517,10 @@
   double x;
   double y;
 
+  bool keyflag = false;         // added
+  bool *key = &keyflag;         // added
+  int ierr = 0;                 // added
+  int *ierror = &ierr;          // added
   *ierror = 0;
 //
 //  On user's signal, set up the factorial table.
@@ -638,7 +642,7 @@
 //
 //  Generate a pseudo-random number.
 //
-      r = r8_uniform_01 ( seed );
+      r = R::runif(0,1); //r8_uniform_01 ( seed );
 //
 //  Compute the conditional expected value of MATRIX(L,M).
 //
@@ -740,7 +744,7 @@
           break;
         }
 
-        r = r8_uniform_01 ( seed );
+        r = R::runif(0,1); // r8_uniform_01 ( seed );
         r = sumprb * r;
 
       }

用于生成 r x c 列联表的算法的 Rcpp 包装器

Rcpp wrapper for algorithm to generate r x c contingency tables

c++

contingency

rcpp

代码

输出

与原始版本不同 asa159.cpp