[Rd] [R] unvectorized option for outer()

[Tony Plate]

When I read the preface to The Blue Book (The New S Language, Becker, 
Chambers & Wilks) I see comments along the lines of "high-level 
language", "primary goal of the S environment is to enable and encourage 
good data analysis", etc.  While vectorization is a great feature of S 
(and R), I don't see it, or programming efficiency, mentioned there at all.

Nonetheless, Peter's suggestion of a general Vectorize() function is 
intriguing, and could be useful with other functions that trip users up 
in the same way.  (Also, I do apprecicate Peter pointing out that not 
all functions vectorize naturally, as in his example density 
calculations over grids of parameters).

So, here's a first pass at a general Vectorize() function:

Vectorize <- function(FUN, vectorize.args) {
     if (!all(is.element(vectorize.args, names(formals(FUN)))))
         stop("some args to vectorize are not args of FUN")
     FUNV <- eval(substitute(function(x, ...) mapply(FUN, x, 
MoreArgs=list(...)), list(FUN=FUN)))
     formals(FUNV) <- formals(FUNV)[c(rep(1, length(vectorize.args)), 2)]
     names(formals(FUNV))[seq(along=vectorize.args)] <- vectorize.args
     body(FUNV) <- body(FUNV)[c(1, 2, rep(3, length(vectorize.args)), 4)]
     body(FUNV)[seq(3,len=length(vectorize.args))] <- 
lapply(vectorize.args, as.name)
     FUNV
}

ssd <- function(A,alpha,Y,t) sum((Y - A*exp(-alpha*t))2)
# SSD is a vectorized version of ssd
SSD <- function(Avec, alphavec, ...) mapply(ssd, Avec, alphavec, 
MoreArgs=list(...))
# Vectorize(ssd, c("A", "alpha")) should produce
# function(A, alpha, ...) mapply(ssd, A, alpha, MoreArgs=list(...))
Y <- 1:5; t <- 3
outer(1:3, 1:2, SSD, Y, t)
outer(1:3, 1:2, Vectorize(ssd, c("A", "alpha")), Y, t)

 > # transcript of running the above commands
 > Vectorize(ssd, c("A", "alpha"))
function (A, alpha, ...)
mapply(function (A, alpha, Y, t)
sum((Y - A * exp(-alpha * t))^2), A, alpha, MoreArgs = list(...))
<environment: 0x1361f40>
 > Y <- 1:5; t <- 3
 > outer(1:3, 1:2, SSD, Y, t)
          [,1]     [,2]
[1,] 53.51878 54.92567
[2,] 52.06235 54.85140
[3,] 50.63071 54.77719
 > outer(1:3, 1:2, Vectorize(ssd, c("A", "alpha")), Y, t)
          [,1]     [,2]
[1,] 53.51878 54.92567
[2,] 52.06235 54.85140
[3,] 50.63071 54.77719
 >

[There are a couple of minor design issues around syntax -- what is the 
best way of specifying the arguments to vectorize? (e.g., what about an 
interface that allowed Vectorize(ssd ~ A * alpha)?), and should the 
function name rather than the definition appear in the result of 
Vectorize()?  But those are issues of secondary importance.]

I have to confess I don't really understand how environments work with 
functions, so I don't know if this Vectorize() function will work in 
general.  What is the appropriate environment for returned value of 
Vectorize()?  Is this approach to creating a Vectorize() function on the 
right tack at all?  Any other improvements or fixes?

-- Tony Plate


Peter Dalgaard wrote:
> Thomas Lumley <tlumley at u.washington.edu> writes:
> 
> 
>>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.
> 
> 
> In fairness, it should probably be said that not all problems
> vectorize naturally. One example is
> 
>   ssd <- function(A,alpha) sum((Y - A*exp(-alpha*t))^2)
> 
> However, it should be worth noting that with the mapply() function at
> hand, it is pretty easy to turn a non-vectorized function into a
> vectorized one. 
> 
>   SSD <- function(A,alpha) mapply(ssd, A, alpha)
> 
> (Anybody want to try their hand on writing a general Vectorize()
> function? I.e. one that allowed
> 
>    outer(Avec, alphavec, Vectorize(ssd))
> 
> to work.)