Re: [S] answers to trying to avoid looping

[Bill Venables <Bill.Venables@cmis.csiro.au>]

 At 17:33 15/09/00 -0400, Joseph S. Verducci wrote:
>
>
>There seems to be a minor problem with the first solution
> and the second solution is limited by the stacking of
> the recursion.  For standard settings, the limit is 15 columns:
>
>> unix.time(Subsets(15,2))
>[1] 0.48 0.01 0.00 0.00 0.00
>
>> unix.time(Subsets(16,2))
>Error: Expressions nested beyond limit (256)
>Timing stopped at: 2.41999999999996 0.00999999999999979 3 0 0 

No, I would not regard those functions as being at all appropriate here.  They were written for quite different purposes.

>
>Combining modifications of both solutions seems to work 
>reasonable well:
>
>Modifying Gardar's function --
>
>Pairs2 <- function(n){
>#
>#  returns an {n(n+1)/2} x 2 matrix of paired integers
>#
>a <- matrix(0,n,n)
>nr <- row(a)
>nc <- col(a)
>ut <- !lower.tri(a)
>cbind(nc[ut],nr[ut])}
>
>Then using Steve's idea
>
>PamG <- function(x){
>#
>#  returns an {nc(nc+1)/2} x nr} matrix whose columns
>#     give the required products 
>#
>apply(Pairs2(ncol(x)),1,function(y) {x[,y[1]] * x[,y[2]]})}

Despite appearances this is really a loop solution.  I suspect your best solution is going to be an even more explicit loop, though.  Here's my suggestion:

> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > Quad <- function(X) {
> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > nc <- ncol(X)
> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > X2 <- matrix(0, nrow(X), (nc * (nc + 1))/2)
> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > m <- 1
> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > for(j in 1:nc) {
> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > X2[, m:(m + nc - j)] <- X[, j] * X[, j:nc]
> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > m <- m + nc - j + 1
> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > }
> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > X2
> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > }

Call as


X2 <- Quad(X)

of course.

This works if X is a data frame or a matrix, but the value, X2, is always a matrix.  If speed were a consideration I would suggest you work with a simple, unadorned matrices X and X2, though, (i.e. no names or dimnames attributes).  

The most important line here is the one setting up the initial X2, namely 

> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > X2 <- matrix(0, nrow(X), (nc * (nc + 1))/2)

It sets aside enough storage to hold the final answer, so you need only make one trip to the allocator.  The loop itself is then a very tiny overhead.

Bill Venables.
Bill Venables, Statistician           Tel. +61 7 3826 7251 
CSIRO Marine Laboratories,                          Fax. +61 7 3826  7304
Cleveland, Qld, 4163                  Email:  Bill.Venables@cmis.csiro.au
AUSTRALIA                        http://www.cmis.csiro.au/bill.venables/ ----------------------------------------------------------------------- This message was distributed by s-news@wubios.wustl.edu.  To unsubscribe send e-mail to s-news-request@wubios.wustl.edu with the BODY of the message:  unsubscribe s-news