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

______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Reply via email to