tensor() function and sets

[Martin Maechler]

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