[R] multiple dimensional diag()
Robin Hankin
rksh at soc.soton.ac.uk
Fri Oct 1 13:24:56 CEST 2004
Hi
I have two arbitrarily dimensioned arrays, "a" and "b", with
length(dim(a))==length(dim(b)). I want to form a sort of
"corner-to-corner" version of abind(), or a multidimensional version
of blockdiag().
In the case of matrices, the function is easy to write and if
a=matrix(1,3,4) and b=matrix(2,2,2), then adiag(a,b) would return:
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 1 1 1 0 0
[2,] 1 1 1 1 0 0
[3,] 1 1 1 1 0 0
[4,] 0 0 0 0 2 2
[5,] 0 0 0 0 2 2
I am trying to generalize this to two higher dimensional arrays.
If x <- adiag(a,b) then I want all(dim(x)==dim(a)+dim(b)); and if
dim(a)=c(a_1, a_2,...a_d) then x[1:a_1,1:a_2,...,1:a_d]=a, and
x[(a_1+1):(a_1+b_1),...,(a_d+1):(a_d+b_d)]=b. Other elements of x are
zero.
The fact that I'm having difficulty expressing this succinctly makes
me think I'm missing something basic.
If a and b have identical dimensions [ie all(dim(a)==dim(b)) ], the
following ghastly kludge (which is one of many) works:
adiag <- function(a,b) {
if(any(dim(a) != dim(b))){stop("a and b must have identical dimensions")}
jj <- array(0,rep(2,length(dim(a))))
jj[1] <- 1
jj[length(jj)] <- 1
jj <- kronecker(jj,b)
f <- function(i){1:i}
do.call("[<-",c(list(jj),sapply(dim(a),f,simplify=FALSE),list(a)))
}
Then "adiag(array(1:8,rep(2,3)),array(-1,rep(2,3)))" is OK. What is
the best way to bind arbitrarily dimensioned arrays together
corner-to-corner?
--
Robin Hankin
Uncertainty Analyst
Southampton Oceanography Centre
SO14 3ZH
tel +44(0)23-8059-7743
initialDOTsurname at soc.soton.ac.uk (edit in obvious way; spam precaution)
More information about the R-help
mailing list