From andy_liaw at merck.com  Mon Oct 31 14:12:35 2005
From: andy_liaw at merck.com (Liaw, Andy)
Date: Mon, 31 Oct 2005 08:12:35 -0500
Subject: [Rd] [R] unvectorized option for outer()
Message-ID: <39B6DDB9048D0F4DAD42CB26AAFF0AFAFED53F@usctmx1106.merck.com>

> From: Thomas Lumley
> 
> On Sun, 30 Oct 2005, Jonathan Rougier wrote:
> 
> > I'm not sure about this.  Perhaps I am a dinosaur, but my feeling is
> > that if people are writing functions in R that might be subject to
> > simple operations like outer products, then they ought to be writing
> > vectorised functions!
> 
> I would agree.  How about an oapply() function that does 
> multiway (rather 
> than just two-way) outer products.  Basing the name on "apply" would 
> emphasize the similarity to other flexible, not particularly 
> optimized 
> second-order functions.
> 
>  	-thomas

I'll toss in my $0.02:  The following is my attempt at creating a "general
outer" that works with more than two arguments.

gouter <- function (x, FUN, ...) {
    xgrid <- as.list(do.call("expand.grid", x))
    names(xgrid) <- NULL
    array(do.call(deparse(substitute(FUN)), c(xgrid, list(...))), 
        dim = sapply(x, length), dimnames = x)
}


Here's an example:

> example(gouter)

gouter> gouter(list(letters[1:2], 1:3, 2:4), paste, sep = "")
, , 2

  1     2     3    
a "a12" "a22" "a32"
b "b12" "b22" "b32"

, , 3

  1     2     3    
a "a13" "a23" "a33"
b "b13" "b23" "b33"

, , 4

  1     2     3    
a "a14" "a24" "a34"
b "b14" "b24" "b34"

Andy

> 
> 
> >		Maybe it's not possible to hold this line, and
> > maybe "outer" is not the right place to draw it, but I 
> think we ought to
> > respect the "x is a vector" mindset as much as possible in the base
> > package.  As Tony notes, the documentation does try to be 
> clear about
> > what outer actually does, and how it can be used.
> >
> > So I am not a fan of the VECTORIZED argument, and 
> definitely not a fan
> > of the VECTORIZED=FALSE default.
> >
> > Jonathan.
> >
> > Gabor Grothendieck wrote:
> >> If the default were changed to VECTORIZED=FALSE then it would
> >> still be functionally compatible with what we have now so 
> all existing
> >> software would continue to run correctly yet would not cause
> >> problems for the unwary.  Existing software would not have 
> to be changed
> >> to add VECTORIZED=TRUE except for those, presumbly few, cases
> >> where outer performance is critical.   One optimization might be to
> >> have the default be TRUE if the function is * or perhaps 
> if its specified
> >> as a single character and FALSE otherwise.
> >>
> >> Having used APL I always regarded the original design of 
> outer in R as
> >> permature performance optimization and this would be a 
> chance to get
> >> it right.
> >>
> >> On 10/28/05, Tony Plate <tplate at acm.org> wrote:
> >>
> >>> [following on from a thread on R-help, but my post here seems more
> >>> appropriate to R-devel]
> >>>
> >>> Would a patch to make outer() work with non-vectorized 
> functions be
> >>> considered?  It seems to come up moderately often on the 
> list, which
> >>> probably indicates that many many people get bitten by the same
> >>> incorrect expectation, despite the documentation and the 
> FAQ entry.  It
> >>> looks pretty simple to modify outer() appropriately: one 
> extra function
> >>> argument and an if-then-else clause to call mapply(FUN, 
> ...) instead of
> >>> calling FUN directly.
> >>>
> >>> Here's a function demonstrating this:
> >>>
> >>> outer2 <- function (X, Y, FUN = "*", ..., VECTORIZED=TRUE)
> >>> {
> >>>    no.nx <- is.null(nx <- dimnames(X <- as.array(X)))
> >>>    dX <- dim(X)
> >>>    no.ny <- is.null(ny <- dimnames(Y <- as.array(Y)))
> >>>    dY <- dim(Y)
> >>>    if (is.character(FUN) && FUN == "*") {
> >>>        robj <- as.vector(X) %*% t(as.vector(Y))
> >>>        dim(robj) <- c(dX, dY)
> >>>    }
> >>>    else {
> >>>        FUN <- match.fun(FUN)
> >>>        Y <- rep(Y, rep.int(length(X), length(Y)))
> >>>        if (length(X) > 0)
> >>>            X <- rep(X, times = ceiling(length(Y)/length(X)))
> >>>        if (VECTORIZED)
> >>>            robj <- FUN(X, Y, ...)
> >>>        else
> >>>            robj <- mapply(FUN, X, Y, MoreArgs=list(...))
> >>>        dim(robj) <- c(dX, dY)
> >>>    }
> >>>    if (no.nx)
> >>>        nx <- vector("list", length(dX))
> >>>    else if (no.ny)
> >>>        ny <- vector("list", length(dY))
> >>>    if (!(no.nx && no.ny))
> >>>        dimnames(robj) <- c(nx, ny)
> >>>    robj
> >>> }
> >>> # Some examples
> >>> f <- function(x, y, p=1) {cat("in f\n"); (x*y)^p}
> >>> outer2(1:2, 3:5, f, 2)
> >>> outer2(numeric(0), 3:5, f, 2)
> >>> outer2(1:2, numeric(0), f, 2)
> >>> outer2(1:2, 3:5, f, 2, VECTORIZED=F)
> >>> outer2(numeric(0), 3:5, f, 2, VECTORIZED=F)
> >>> outer2(1:2, numeric(0), f, 2, VECTORIZED=F)
> >>>
> >>> # Output on examples
> >>>
> >>>> f <- function(x, y, p=1) {cat("in f\n"); (x*y)^p}
> >>>> outer2(1:2, 3:5, f, 2)
> >>>
> >>> in f
> >>>     [,1] [,2] [,3]
> >>> [1,]    9   16   25
> >>> [2,]   36   64  100
> >>>
> >>>> outer2(numeric(0), 3:5, f, 2)
> >>>
> >>> in f
> >>>     [,1] [,2] [,3]
> >>>
> >>>> outer2(1:2, numeric(0), f, 2)
> >>>
> >>> in f
> >>>
> >>> [1,]
> >>> [2,]
> >>>
> >>>> outer2(1:2, 3:5, f, 2, VECTORIZED=F)
> >>>
> >>> in f
> >>> in f
> >>> in f
> >>> in f
> >>> in f
> >>> in f
> >>>     [,1] [,2] [,3]
> >>> [1,]    9   16   25
> >>> [2,]   36   64  100
> >>>
> >>>> outer2(numeric(0), 3:5, f, 2, VECTORIZED=F)
> >>>
> >>>     [,1] [,2] [,3]
> >>>
> >>>> outer2(1:2, numeric(0), f, 2, VECTORIZED=F)
> >>>
> >>> [1,]
> >>> [2,]
> >>>
> >>> If a patch to add this feature would be considered, I'd 
> be happy to
> >>> submit one (including documentation).  If so, and if there are any
> >>> potential traps I should bear in mind, please let me know!
> >>>
> >>> -- Tony Plate
> >>>
> >>> Rau, Roland wrote:
> >>>
> >>>> Dear all,
> >>>>
> >>>> a big thanks to Thomas Lumley, James Holtman and Tony 
> Plate for their
> >>>> answers. They all pointed in the same direction => I 
> need a vectorized
> >>>> function to be applied. Hence, I will try to work with a 
> 'wrapper'
> >>>> function as described in the FAQ.
> >>>>
> >>>> Thanks again,
> >>>> Roland
> >>>>
> >>>>
> >>>>
> >>>>
> >>>>> -----Original Message-----
> >>>>> From: Thomas Lumley [mailto:tlumley at u.washington.edu]
> >>>>> Sent: Thursday, October 27, 2005 11:39 PM
> >>>>> To: Rau, Roland
> >>>>> Cc: r-help at stat.math.ethz.ch
> >>>>> Subject: Re: [R] outer-question
> >>>>>
> >>>>>
> >>>>> You want FAQ 7.17 Why does outer() behave strangely 
> with my function?
> >>>>>
> >>>>>     -thomas
> >>>>>
> >>>>> On Thu, 27 Oct 2005, Rau, Roland wrote:
> >>>>>
> >>>>>
> >>>>>
> >>>>>> Dear all,
> >>>>>>
> >>>>>> This is a rather lengthy message, but I don't know what I
> >>>>>
> >>>>> made wrong in
> >>>>>
> >>>>>
> >>>>>> my real example since the simple code works.
> >>>>>> I have two variables a, b and a function f for which I 
> would like to
> >>>>>> calculate all possible combinations of the values of a and b.
> >>>>>> If f is multiplication, I would simply do:
> >>>>>>
> >>>>>> a <- 1:5
> >>>>>> b <- 1:5
> >>>>>> outer(a,b)
> >>>>>>
> >>>>>> ## A bit more complicated is this:
> >>>>>> f <- function(a,b,d) {
> >>>>>>    return(a*b+(sum(d)))
> >>>>>> }
> >>>>>> additional <- runif(100)
> >>>>>> outer(X=a, Y=b, FUN=f, d=additional)
> >>>>>>
> >>>>>> ## So far so good. But now my real example. I would 
> like to plot the
> >>>>>> ## log-likelihood surface for two parameters alpha and beta of
> >>>>>> ## a Gompertz distribution with given data
> >>>>>>
> >>>>>> ### I have a function to generate random-numbers from a
> >>>>>> Gompertz-Distribution
> >>>>>> ### (using the 'inversion method')
> >>>>>>
> >>>>>> random.gomp <- function(n, alpha, beta) {
> >>>>>>          return( (log(1-(beta/alpha*log(1-runif(n)))))/beta)
> >>>>>> }
> >>>>>>
> >>>>>> ## Now I generate some 'lifetimes'
> >>>>>> no.people <- 1000
> >>>>>> al <- 0.1
> >>>>>> bet <- 0.1
> >>>>>> lifetimes <- random.gomp(n=no.people, alpha=al, beta=bet)
> >>>>>>
> >>>>>> ### Since I neither have censoring nor truncation in this
> >>>>>
> >>>>> simple case,
> >>>>>
> >>>>>
> >>>>>> ### the log-likelihood should be simply the sum of the 
> log of the
> >>>>>> ### the densities (following the parametrization of
> >>>>>
> >>>>> Klein/Moeschberger
> >>>>>
> >>>>>
> >>>>>> ### Survival Analysis, p. 38)
> >>>>>>
> >>>>>> loggomp <- function(alphas, betas, timep) {
> >>>>>> return(sum(log(alphas) + betas*timep + (alphas/betas *
> >>>>>> (1-exp(betas*timep)))))
> >>>>>> }
> >>>>>>
> >>>>>> ### Now I thought I could obtain a matrix of the
> >>>>>
> >>>>> log-likelihood surface
> >>>>>
> >>>>>
> >>>>>> ### by specifying possible values for alpha and beta 
> with the given
> >>>>>> data.
> >>>>>> ### I was able to produce this matrix with two for-loops.
> >>>>>
> >>>>> But I thought
> >>>>>
> >>>>>
> >>>>>> ### I could use also 'outer' in this case.
> >>>>>> ### This is what I tried:
> >>>>>>
> >>>>>> possible.alphas <- seq(from=0.05, to=0.15, length=30)
> >>>>>> possible.betas <- seq(from=0.05, to=0.15, length=30)
> >>>>>>
> >>>>>> outer(X=possible.alphas, Y=possible.betas, FUN=loggomp,
> >>>>>
> >>>>> timep=lifetimes)
> >>>>>
> >>>>>
> >>>>>> ### But the result is:
> >>>>>>
> >>>>>>
> >>>>>>> outer(X=possible.alphas, Y=possible.betas, FUN=loggomp,
> >>>>>>
> >>>>>> timep=lifetimes)
> >>>>>> Error in outer(X = possible.alphas, Y = possible.betas, FUN
> >>>>>
> >>>>> = loggomp,
> >>>>>
> >>>>>
> >>>>>> :
> >>>>>>      dim<- : dims [product 900] do not match the length
> >>>>>
> >>>>> of object [1]
> >>>>>
> >>>>>
> >>>>>> In addition: Warning messages:
> >>>>>> ...
> >>>>>>
> >>>>>> ### Can somebody give me some hint where the problem is?
> >>>>>> ### I checked my definition of 'loggomp' but I thought this
> >>>>>
> >>>>> looks fine:
> >>>>>
> >>>>>
> >>>>>> loggomp(alphas=possible.alphas[1], betas=possible.betas[1],
> >>>>>> timep=lifetimes)
> >>>>>> loggomp(alphas=possible.alphas[4], betas=possible.betas[10],
> >>>>>> timep=lifetimes)
> >>>>>> loggomp(alphas=possible.alphas[3], betas=possible.betas[11],
> >>>>>> timep=lifetimes)
> >>>>>>
> >>>>>>
> >>>>>> ### I'd appreciate any kind of advice.
> >>>>>> ### Thanks a lot in advance.
> >>>>>> ### Roland
> >>>>>>
> >>>>>>
> >>>>>> +++++
> >>>>>> This mail has been sent through the MPI for Demographic
> >>>>>
> >>>>> Rese...{{dropped}}
> >>>>>
> >>>>>
> >>>>>> ______________________________________________
> >>>>>> R-help at stat.math.ethz.ch mailing list
> >>>>>> https://stat.ethz.ch/mailman/listinfo/r-help
> >>>>>> PLEASE do read the posting guide!
> >>>>>
> >>>>> http://www.R-project.org/posting-guide.html
> >>>>>
> >>>>> Thomas Lumley                 Assoc. Professor, Biostatistics
> >>>>> tlumley at u.washington.edu      University of Washington, Seattle
> >>>>>
> >>>>
> >>>>
> >>>> +++++
> >>>> This mail has been sent through the MPI for Demographic 
> Research.  Should you receive a mail that is apparently from 
> a MPI user without this text displayed, then the address has 
> most likely been faked. If you are uncertain about the 
> validity of this message, please check the mail header or ask 
> your system administrator for assistance.
> >>>>
> >>>>
> >>>
> >>> ______________________________________________
> >>> R-help at stat.math.ethz.ch mailing list
> >>> https://stat.ethz.ch/mailman/listinfo/r-help
> >>> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html
> >>>
> >>
> >>
> >> 
> ______________________________________________
> >> R-devel at r-project.org mailing list
> >> https://stat.ethz.ch/mailman/listinfo/r-devel
> >
> > -- 
> > Jonathan Rougier                       Science Laboratories
> > Department of Mathematical Sciences    South Road
> > University of Durham                   Durham DH1 3LE
> > tel: +44 (0)191 334 3111, fax: +44 (0)191 334 3051
> > http://www.maths.dur.ac.uk/stats/people/jcr/jcr.html
> >
> > ______________________________________________
> > R-devel at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-devel
> >
> 
> Thomas Lumley			Assoc. Professor, Biostatistics
> tlumley at u.washington.edu	University of Washington, Seattle
> 
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
> 
>



[Rd] [R] unvectorized option for outer()
Liaw, Andy