From murdoch at stats.uwo.ca  Mon Oct 31 21:52:32 2005
From: murdoch at stats.uwo.ca (Duncan Murdoch)
Date: Mon, 31 Oct 2005 15:52:32 -0500
Subject: [Rd] [R] unvectorized option for outer()
In-Reply-To: <43666D4D.7070602@acm.org>
References: <8B08A3A1EA7AAC41BE24C750338754E69FDF14@HERMES.demogr.mpg.de>	<43629490.9060204@acm.org>	<971536df0510281643u50929fa6rf097c40948700804@mail.gmail.com>	<4364A4FC.2030400@durham.ac.uk>	<Pine.LNX.4.63a.0510301300260.28148@homer23.u.washington.edu>	<x2y84amyv4.fsf@turmalin.kubism.ku.dk>
	<43666D4D.7070602@acm.org>
Message-ID: <43668410.8060403@stats.uwo.ca>

Here's my entry for Vectorize:

Vectorize <- function(FUN, vectorize.args = names(fargs)) {

      fargs <- formals(FUN)
      dovec <- match(vectorize.args, names(fargs))

      FUNV <- function() { # will set the formals below
         args <- as.list(match.call())[-1]
         args <- lapply(args, eval, parent.frame())
      	mapply.args <- c(FUN = FUN,
      			    args[dovec],
      			    MoreArgs = list(args[-dovec]))
         do.call("mapply", mapply.args)
      }
      formals(FUNV) <- fargs
      FUNV
  }


I *think* this works fairly generally, but I'm not 100% sure of that, 
especially about how it handles default parameter values, e.g. in a 
weird case like

 > ssd <- function(A,alpha,Y,t=Time) {
+     Time <- 3
+     sum((Y - A*exp(-alpha*t))^2)
+ }
 >
 > SSD <- Vectorize(ssd, c("A", "alpha"))

It appears to give the right answer:

 > outer(1:3, 1:2, SSD, Y)
          [,1]     [,2]
[1,] 53.51878 54.92567
[2,] 52.06235 54.85140
[3,] 50.63071 54.77719

but that's a bit of a surprise, and I wouldn't be too surprised if there 
are other weird examples where it screws up.

Duncan Murdoch





On 10/31/2005 2:15 PM, Tony Plate wrote:
> 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.)
> 
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel



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