>>>>> On Tue, 20 Jul 1999 15:26:35 +0100 (BST), Jonathan Rougier
<J.C.Rougier@durham.ac.uk> said:
JonR> Hi Everyone,
JonR> To complete the outer() and kronecker() functions in the base, may
I
JonR> suggest the following tensor() function, which allows the
multiplication
JonR> of arrays through sets of conformable dimensions. I am happy to
write a
JonR> help page if required.
JonR> The code also needs a setdiff() function which prompts me to ask:
JonR> what about simple set functions? I expect many of us have written
our own
JonR> (Brian has a setdiff() in drop1.lm(), for example), which seems
like a
JonR> good reason for putting versions in the base. I would be happy to
provide
JonR> mine for general scrutiny.
Actually your version of "setdiff" is in
example("%in%")
"%w/o%" <- function(x,y) x[!x %in% y] #-- x without y
Maybe we should make this part of core R, since it allows the nicer syntax
(1:10) %w/o% c(3,7,12).
One problem with all these might be that you forget to think about efficiency
{ %w/o% -> %in% -> match(.) }
but that might be ok.
union
and
intersect
are two other set functions, the first of which is also contained in
example(match) { == example("%in%") }:
intersect <- function(x, y) y[match(x, y, nomatch = 0)]
What do others think? Should (some of) these be put into "core R" ?
JonR> Jonathan Rougier Science Laboratories
JonR> Department of Mathematical Sciences South Road
JonR> University of Durham Durham DH1 3LE
JonR> "setdiff" <-
JonR> function (x, y) x[!(x %in% y)]
JonR> "tensor" <-
JonR> function (A, B, da, db)
JonR> {
JonR> # tensor product of A and B through da and db
JonR> no.na <- is.null(na <- dimnames(A <- as.array(A)))
JonR> dima <- dim(A)
JonR> no.nb <- is.null(nb <- dimnames(B <- as.array(B)))
JonR> dimb <- dim(B)
JonR> if (any(dima[da] != dimb[db]))
JonR> stop("Mismatched dimensions")
JonR> kpa <- setdiff(seq(along = dima), da)
JonR> kpb <- setdiff(seq(along = dimb), db)
JonR> # fix up the dimnames (see outer)
JonR> if (no.na && no.nb)
JonR> nms <- NULL
JonR> else {
JonR> if (no.na)
JonR> na <- vector("list", length(dima))
JonR> else if (no.nb)
JonR> nb <- vector("list", length(dimb))
JonR> nms <- c(na[kpa], nb[kpb])
JonR> }
JonR> A <- matrix(aperm(A, c(kpa, da)), ncol = prod(dima[da]))
JonR> B <- matrix(aperm(B, c(db, kpb)), nrow = prod(dimb[db]))
JonR> array(A %*% B, c(dima[kpa], dimb[kpb]), dimnames = nms)
JonR> }
JonR> # examples
JonR> A <- array(1:6, 2:3, list(LETTERS[1:2], LETTERS[3:5]))
JonR> B <- 1:3
JonR> tensor(A, B, 2, 1) # same as drop(A %*% B)
JonR> A <- outer(A, B) # A now 2 by 3 by 3
JonR> tensor(A, B, 2, 1); tensor(A, B, 3, 1) # both 2 by 3 (nb dimnames)
JonR> B <- outer(1:2, B) # B now 2 by 3
JonR> tensor(A, B, c(1, 3), 1:2) # must be 3 vector
JonR> B <- outer(1:3, B) # B now 3 by 2 by 3
JonR> tensor(A, B, 1:3, c(2, 1, 3)) # must be a scalar
JonR> sum(A * aperm(B, c(2, 1, 3))) # same scalar value
Looks interesting, however is it really useful
{for doing statistics, graphics, numerical computing, I mean} ?
Yes, if you'd provide a *.Rd help file (use ``prompt(tensor)''),
we'd be even more inclined to include the tensor function....
Also, I'd like to see a "real" use for it.
Is
tensor(A, B, 1:3, c(2, 1, 3))
really easier than
sum(A * aperm(B, c(2, 1, 3)))
[probably yes, you have to think less, but still:
couldn't "da" and "db" (3rd and 4th arg) get
reasonable default values ?]
----
Martin
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To:
r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._