[R] multiple dimensional diag()

[Robin Hankin]

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)