Rcpp wrapper for algorithm to generate r x c contingency tables

Question

I'm trying to use this C++ implementation of an algorithm for efficiently generating r x c contingency tables with fixed margins: https://people.sc.fsu.edu/%7Ejburkardt/cpp_src/asa159/asa159.html

I'm working in Rcpp and am trying to create a function that I can utilize in R. I'm aware that there's a C implementation of this algorithm in R. That doesn't work for me because I need to be able to utilize this function within Rcpp. As an intermediate step, I'm just trying to wrap this algorithm and export it to R. Seems simple but I've struggled mighty with this.

The original function is defined as follows:

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

This is a link directly to the algorithm: https://people.sc.fsu.edu/%7Ejburkardt/cpp_src/asa159/asa159.cpp

I don't need the error handling and couldn't figure out how to deal with the void type so I altered the algorithm slightly and ended up with this:

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

I've made several attempts. Here's one example. This will compile, but when I try to run the function in R, it crashes.

#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);

  }

I assume I'm making some elementary mistakes here. Thank you in advance to anyone who has time to take a look.

Answer 1

This is actually an excellent question, and I too have had good luck with the nice and very comprehensive site full of routines by John Burkardt at FSU.

But nothing is life is entirely free, and one needs to know a little bit about how C and C++ work to integrate such routines into R via Rcpp. My version of the final file is much shorter: we just ship rowsums and colsums to it, and get the desired matrix back. We also add the R call for the example at the bottom (a nice trick!)

Code

#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
*/

Output

This contains three example calls.

> 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
> 

Diff to original asa159.cpp

The change to the original asa159.cpp is fairly small: we use R's own RNG, and we simplify the calling interface. Propagation of error condition via Rcpp::stop() is left as an exercise, as they say :)

--- 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;
 
       }