( 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
_______________________________________________
Rcpp-devel mailing list
[email protected]
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
