[Rcpp-devel] [R] Call by reference or suggest workaround

Romain Francois romain at r-enthusiasts.com
Sat Jun 19 17:32:34 CEST 2010


( Reposting this here to Rcpp-devel since this might be of interest to 
the readers of this list. )

Le 19/06/10 16:32, Chidambaram Annamalai a écrit :
 >
 > I have written code to compute multi-indices in R [1] and due to the
 > recursive nature of the computation I need to pass around the *same*
 > matrix object (where each row corresponds to one multi-index). As pass
 > by reference wasn't the default behavior I declared a global matrix
 > (mat) and used the<<- operator to write to the global matrix. So the
 > usage would be to call genMultiIndices(3,2) for side effects to
 > generate all multi-indices of length 3 and sum 2. And then access the
 > global matrix.
 >
 > However, after coding this I can't seem to export the global matrix
 > object (in the NAMESPACE file) and still retain mutability since its
 > binding is "locked" (R throws an error). Can I somehow unlock this?
 >
 > Ideally I would want to pass around the same matrix to the recursive
 > function. Is that possible? If not, could someone please suggest a
 > workaround to use the code in an R package?
 >
 > [1]: http://dpaste.com/209186/

Hi,

You can use lexical scoping and you might like ?Recall as well.

genMultiIndices <- function(N, V) {
     mat <- matrix(nrow=choose(N + V - 1, V), ncol=N)
     fillMultiIndices <- function(i, j, n, v) {
         if (n == 1) {
             mat[i, j] <<- v
         }
         else if (v == 0) {
             mat[i, j:(j + n - 1)] <<- 0L
         }
         else {
             rowOffset <- 0
             # the first element of each multi-index can be any of 0, 1, 
..., v
             for (k in v:0) {
                 times <- choose((n - 1) + (v - k) - 1, (v - k))
                 mat[(i + rowOffset):(i + rowOffset + times - 1), j] <<- k
                 Recall(i + rowOffset, j + 1, n - 1, v - k)
                 rowOffset <- rowOffset + times
             }
         }
     }
     fillMultiIndices(1, 1, N, V)
     mat
}



Also, you can consider writing your code in a language that supports 
references, e.g. C++. Here is a start with inline/Rcpp :

require( inline )
require( Rcpp )

genMultiIndices_internal <-  local({
inc <- '
void fillMultiIndices( Rcpp::IntegerMatrix& mat, int i, int j, int n, 
int v ){
     if( n == 1 ){
         mat( (i-1), (j-1) ) = v ;
     } else if( v == 0 ){
         for( int k=j; k < j+n; k++){
             mat( (i-1), (k-1) ) = 0 ;
         }
     } else {

         // using the R function
         // I leave it to you to use a C implementation
         Function choose( "choose" ) ;

         int rowOffset = 0 ;
         int times ;
         for( int k=v; k>=0; k--){
             times = as<int>( choose( (n-1) + (v-k) - 1, (v-k) ) );
             int start = i + rowOffset ;
             int end   = i + rowOffset + times ;
             for( int z = start; z < end; z++ ){
                 mat( z-1 , j-1 ) = k ;
             }
             fillMultiIndices( mat, i + rowOffset, j+1, n-1, v-k ) ;
             rowOffset += times ;
         }
     }
}
'
code <- '
     int N  = as<int>( N_) ;
     int V  = as<int>( V_) ;
     int NR = as<int>( NR_) ;
     Rcpp::IntegerMatrix mat( NR, N ) ;
     fillMultiIndices( mat, 1, 1, N, V ) ;
     return mat ;
'
     .genMultiIndices <- cxxfunction(
         signature( N_ = "integer", V_ = "integer", NR_ = "integer" ),
         code, include = inc, plugin = "Rcpp" )
     function( N, V){
         .genMultiIndices( N, V, choose(N + V - 1, V) )
     }
} )


I've been lazy here and I am using the choose function from R, so there 
is room for some improvement.

 > ( x <- genMultiIndices( 3L , 2L )          )
      [,1] [,2] [,3]
[1,]    2    0    0
[2,]    1    1    0
[3,]    1    0    1
[4,]    0    2    0
[5,]    0    1    1
[6,]    0    0    2
 > ( y <- genMultiIndices_internal( 3L, 2L )  )
      [,1] [,2] [,3]
[1,]    2    0    0
[2,]    1    1    0
[3,]    1    0    1
[4,]    0    2    0
[5,]    0    1    1
[6,]    0    0    2
 > identical( x, y )
[1] TRUE

Romain
-- 
Romain Francois
Professional R Enthusiast
+33(0) 6 28 91 30 30
http://romainfrancois.blog.free.fr
|- http://bit.ly/98Uf7u : Rcpp 0.8.1
|- http://bit.ly/c6YnCi : graph gallery collage
`- http://bit.ly/bZ7ltC : inline 0.3.5



More information about the Rcpp-devel mailing list