A pull request for a new entry?

[rsparapa]

Sorry, I am too pitiful with git to do a pull request. Here's the file that I'm trying to pull. 2020-07-18-accessing-3d-array.cpp

/**
 * @title Accessing entries of a three-dimensional array
 * @author Rodney Sparapani
 * @license GPL (>= 2)
 * @tags array
 * @summary An example of how to access the entries of an array 
 *   since generic (i, j, k)-like operators don't exist.
 */

#include <Rcpp.h>
// [[Rcpp::export]]
using namespace Rcpp;
IntegerVector get3d(IntegerVector x, IntegerVector args) {
  IntegerVector rc(1), dim=x.attr("dim");
  rc[0]=R_NaN;
  const size_t K=args.size();
  if(K!=dim.size()) return rc;
  size_t i;
  if(K==1) {
    i=args[0]-1;
    if(i>=0 && i<dim[0]) rc[0]=x[i];
  }
  else if(K==2) {
    i=(args[1]-1)*dim[0]+args[0]-1;
    if(i>=0 && i<(dim[0]*dim[1])) rc[0]=x[i];
  }
  else if(K==3) {
    i=((args[2]-1)*dim[1]+(args[1]-1))*dim[0]+args[0]-1;
    if(i>=0 && i<(dim[0]*dim[1]*dim[2])) rc[0]=x[i];
  }
  return rc;
}

/**
 * This function returns the entry corresponding to "args" from "x" which 
 * is either a 1, 2 or 3-dimensional arrays.  Since the `(i, j, k)` operator
 * doesn't exist, we resort to "args" which is an integer vector.
 */

/*** R
library(Rcpp)

b = array(1:8, dim=8)
c = array(1:8, dim=c(2, 4))
a <- array(1:24, dim=c(2, 3, 4))

get3d(b, 3)
b[3]
for(k in 1:8) print(b[k]-get3d(b, k))

get3d(c, c(2, 1))
c[2, 1]
for(i in 1:2)
        for(k in 1:4) print(c[i, k]-get3d(c, c(i, k)))

get3d(a, c(2, 1, 1))
a[2, 1, 1]

for(i in 1:2)
    for(j in 1:3)
        for(k in 1:4) 
            print(a[i, j, k]-get3d(a, c(i, j, k)))
*/

[eddelbuettel]

Providing it as a file or gist or .. is actually also easier for me too than committing forcefully but in the wrong directory :) as happened recently. It would be more in line with the style of the gallery if it has more than 30 some words.

Intro, problem definition, discussion of approach, possible gotchas, ... there are well over a 100 example to guide you.

[rsparapa]

Oh, I didn't realize that we could attach files. Here you go... 2020-07-18-accessing-3d-array.txt

[eddelbuettel]

Can you maybe add an opening paragraph motivating this? (As an aide, explaining why arma::cube did not work you may help there, but we leave that for another day,)

Rudimentary error checking would be nice: 4d or 5d args vector will work silently but maybe not as indented.

[rsparapa]

More than the dimension returns NA which is similar to the way we handled errors in Quantile, i.e., silently returning nothing of value. What more do we need to do here? This is not even an Rcpp function so doing more seems like over-kill. It just to show other users how to do this without having to spend hours debugging it (like I have sadly). Perhaps it is trivial, but I have wanted this for more than a year so perhaps others do as well. Now, I have explained why Armadillo/Eigen are not inviting... 2020-07-18-accessing-3d-array.txt

[eddelbuettel]

One quick look at RcppArmadillo unit tests leads to cpp/cube.cpp the first test function is easily converted as a one-liner:

R> Rcpp::cppFunction("arma::cube cube_test(const arma::cube& x) { return arma::pow(x, 2); }", \
                     depends="RcppArmadillo")
R> cube_test( array(1:8, dim=c(2,2,2)) )
, , 1

     [,1] [,2]
[1,]    1    9
[2,]    4   16

, , 2

     [,1] [,2]
[1,]   25   49
[2,]   36   64

R> 

So we clearly get arma::cube objects in an out, and hence also have the Armadillo functionality available. Works for me. We always stated relatively clearly (if somewhat apologetically) that you should not try to use Rcpp vectors or matrices for direct linear algebra work.

If I find time, I may try to edit and enhance your draft into something suitable for the Rcpp Gallery. As it stands, it is little too short and lacking motivation. I appreciate that you took the time to write something up, and that you kept the submission straightforward.

[rsparapa]

Sure, go ahead. But, like I said, having an additional dependency is what I am trying to avoid. And I am more of an Eigen fan than Armadillo. Even so, I work on these R packages with a team and we dread dependencies. Furthermore, some favor Eigen and some Armadillo so, currently, we have neither!

[eddelbuettel]

Sure, you can do as you please. If you drop Rcpp, you even have one dependency less. And if you drop R, then you're down to only a C compiler. Better, no?

(That was very tongue-in-cheek. I also believe in lightweight is the right weight but I also believe in using the right tools -- it is a tradeoff everyone has to make. I would take Armadillo a hundred times over home-cooked matrix etc types. Been there, done that, and there is a reason we never found that avenenue all that exciting for Rcpp itself. Your loss I presume.)

[rsparapa]

Well, these are not home-cooked types: these arrays are already present in R and Rcpp (as you obviously know). I'm just trying to make use of them. But, since you are a strong advocate for Rcpp (as you should be and I am right there with you), maybe you meant something else.

[eddelbuettel]

I meant what I wrote. It is not a good idea to do (linear algebra) math with Rcpp because you then end up hand-cooking one-off functions that cost you a day, had no other eyeballs on them, and are likely buggy.

Or you could rely on code used by 100s of thousands of people, that widely tested, robust and (on the margin) likely faster (due to some OpenMP). I do what 750+ other packages on CRAN and enjoy the new zero-depends apart from its header RcppArmadillo package.. Some people like RcppEigen too, I tend(ed) to stick with RcppArmadillo.

On balance, I think, this is where the story should end. There is simply no new Rcpp Gallery piece here. I appreciate all the work you did, but let's all move on with our other projects. We're 20 coments in and not broken any relevant new ground. So maybe time to call it a day.

[rsparapa]

Now, I see what you mean. But, I'm not doing linear algebra. I have these huge MCMC objects that are 3-dimensional that I just want to return to R. Furthermore, these methods generally don't have any linear algebra needs which is why we can get away with not having dependencies on Eigen or Armadillo (like I said, I love Eigen, but it is not necessary here).