[R] How to optimize this codes ?

[William Dunlap]

> -----Original Message-----
> From: William Dunlap 
> Sent: Thursday, December 04, 2008 9:59 AM
> To: 'daren76 at hotmail.com'
> Cc: 'R help'
> Subject: Re: [R] How to optimize this codes ?
> 
>    [R] How to optimize this codes ?
>    Daren Tan daren76 at hotmail.com
>    Thu Dec 4 17:02:49 CET 2008
> 
>    How to optimize the for-loop to be reasonably fast for 
> sample.size=100000000 ? 
>    You may want to change sample.size=1000 to have an idea 
> what I am achieving.  
>  
>    set.seed(143)
>    A <- matrix(sample(0:1, sample.size, TRUE), ncol=10, 
> dimnames=list(NULL, LETTERS[1:10]))
>  
>    B <- list()
>    for(i in 1:10) { 
>      B[[i]] <- apply(combn(LETTERS[1:10], i), 2, function(x) 
> { sum(apply(data.frame(A[,x]), 1, all)) })
>    }
> 
> Instead computing all(A[row,x]) each row of the matrix by looping
> over the rows you could loop over the columns.  That is generally
> quicker if there are many more rows than columns.  S+ and R have
> functions called pmin and pmax that do this for min and max, but
> no pany or pall functions.  In your 0/1 case all is the same as min
> so you can replace
>      apply(data.frame(A[,x]), 1, all)
> with
>      sum(do.call("pmin", unname(A[,x,drop=FALSE])))
> after first converting A to a data.frame before starting to compute B.
> (The do.call/unname combo is a hack to pass all the columns 
> of a data.frame
> as separate arguments to a function.  If pmin took a list 
> that wouldn't
> be necessary.)
> 
> When you do timings, looking at one size of problem often 
> isn't helpful,
> as it doesn't tell you how the time depends on the size of 
> the input.  The
> relationship often is not linear.  E.g., I put your method in 
> a function
> called computeB0 and my modification of it into computeB1 and 
> wrote a makeA
> that makes an A matrix with a given sample.size.  Here are 
> the times, it looks
> like 1000 is below the linear range for this example:
> 
> sample.size       time:computeB0   time:computeB1
>        1000                 5.36             0.98
>       10000                36.82             2.49
>      100000               381.34            18.97
> 
> and here are the functions used
> 
> makeA <-
> function(sample.size){
>    set.seed(143);
>    A<-matrix(sample(0:1,size=sample.size,replace=TRUE), 
> ncol=10, dimnames=list(NULL, LETTERS[1:10]));
>    A
> }
> computeB0 <-
> function(A){B<-list()
>   for(i in 1:10) {
>     B[[i]] <- apply(combn(LETTERS[1:10], i), 2, function(x) { 
> sum(apply(data.frame(A[,x]), 1, all)) })
>   }
>   B
> }
> computeB1 <-
> function(A){
>   A<-as.data.frame(A)
>   B<-list()
>   for(i in 1:10) {
>     B[[i]] <- apply(combn(LETTERS[1:10], i), 2,
>                   FUN=function(x) { sum(do.call("pmin", 
> unname(A[,x,drop=FALSE]))) }
>             )
>   }
>   B
> }
> 
> You could probably same more time by restructuring this as a 
> recursive function,
> still operating on columns, that traversed the binary tree of 
> column inclusion/exclusion.
> I just wanted to point out the advantage of looping over 
> columns instead of looping over
> rows.

A hastily written recursive version is:

computeB3 <-
function(A) {
   B <- rep(list(integer(0)), ncol(A))
   storage.mode(A) <- "logical"
   recurse <- function(thisCol = 1, includedCols = rep(NA, ncol(A)),
allSoFar = rep(TRUE, nrow(A))) {
      if (thisCol < ncol(A)) # skip this column
          Recall(thisCol+1, replace(includedCols, thisCol, FALSE),
allSoFar = allSoFar)
      else if (thisCol > ncol(A))
          return()
      includedCols[thisCol] <- TRUE
      allSoFar <- allSoFar & A[,thisCol]
      nIncludedCols <- sum(includedCols[1:thisCol])
      # cat(thisCol, ":", nIncludedCols, ":", includedCols[1:thisCol],
sum(allSoFar), "\n")
      # Note B<<- in next line uses lexical scoping to change
computeB3:B.  Does not work in S+.
      B[[nIncludedCols]][length(B[[nIncludedCols]])+1] <<- sum(allSoFar)
      Recall(thisCol+1, includedCols, allSoFar = allSoFar)
   }
   recurse()
   B
}

The elements of B[[i]], for each i, are in a different order than they
are in
the other functions, but the histograms are the same.

I added its time to the above table:

sample.size       time:computeB0   time:computeB1    time:computeB3
       1000                 5.36             0.98              0.12    
      10000                36.82             2.49              0.25
     100000               381.34            18.97              1.62
    1000000                    ?           290.21             16.36

Bill Dunlap
TIBCO Software Inc - Spotfire Division
wdunlap tibco.com