unvectorized option for outer()

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

example(gouter)

	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}}

+++++
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.