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)
Here is a function that does what you want (perhaps):
<<*>>a.block.diag <- function(a,b) {
# a.block.daig combines arrays a and b and builds a new array which has
# a and b as blocks on its diagonal. pw 10/2004
if(length(dim.a<-dim(a))!= length(dim.b<-dim(b))){
stop("a and b must have identical number of dimensions")}
s<-array(0,dim.a+dim.b)
s<-do.call("[<-",c(list(s),lapply(dim.a,seq),list(a)))
ind<-lapply(seq(dim.b),function(i)seq(dim.b[[i]])+dim.a[[i]])
do.call("[<-",c(list(s),ind,list(b)))
}
@
Try:
<<*>>a=matrix(1,3,4); b=matrix(2,2,2)
a.block.diag(a,b)
@
output-start
Fri Oct 1 17:20:26 2004
[,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
output-end
and an another example:
<<*>>xx<-array(1:8,c(2,rep(2,2)));
yy<-array(-1*(1:9),c(2,rep(3,2)))
a.block.diag(xx,yy)
@
output-start
Fri Oct 1 17:21:38 2004
, , 1
[,1] [,2] [,3] [,4] [,5]
[1,] 1 3 0 0 0
[2,] 2 4 0 0 0
[3,] 0 0 0 0 0
[4,] 0 0 0 0 0
, , 2
[,1] [,2] [,3] [,4] [,5]
[1,] 5 7 0 0 0
[2,] 6 8 0 0 0
[3,] 0 0 0 0 0
[4,] 0 0 0 0 0
, , 3
[,1] [,2] [,3] [,4] [,5]
[1,] 0 0 0 0 0
[2,] 0 0 0 0 0
[3,] 0 0 -1 -3 -5
[4,] 0 0 -2 -4 -6
, , 4
[,1] [,2] [,3] [,4] [,5]
[1,] 0 0 0 0 0
[2,] 0 0 0 0 0
[3,] 0 0 -7 -9 -2
[4,] 0 0 -8 -1 -3
, , 5
[,1] [,2] [,3] [,4] [,5]
[1,] 0 0 0 0 0
[2,] 0 0 0 0 0
[3,] 0 0 -4 -6 -8
[4,] 0 0 -5 -7 -9
output-end
a.block.diag is not always the best solution. In case of
x<-array(1:8,rep(2,3)); y<-array(-1,rep(2,3))
the function adiag will be a little bit faster.
Peter Wolf
Robin Hankin wrote:
> 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 wrote:> hi again Peter > > thanks for this. Now we've got the legal stuff sorted, on to business! > > One of the reasons I wanted the function was to to multidimensional > moving- > window averaging, and for this I needed to be able to call > a.block.diag(a,b) when a and b might > have one or more dimensions of zero extent. > > I also figure that a scalar argument should be interpreted as having > dimensions > rep(1,d) where d=length(dim(<other matrix>)), and that a.block.diag(n,m) > should return plain old diag(n,m) if both n and m are scalars. > > I've modified your function do to this, and added a padding value > argument (function > and a couple of calls pasted below) > > what do you think?perfect! ... and if you want to combine more than two arrays you can add a function like mbd mbd<-function(...,pad=0){ result<-(args<-list(...))[[1]] for(li in args[-1]) result<-a.block.diag(result,li,pad=pad) return(result) } a <- matrix(1:4,2,2) mbd(a,a,a,a,pad=10) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [1,] 1 3 10 10 10 10 10 10 [2,] 2 4 10 10 10 10 10 10 [3,] 10 10 1 3 10 10 10 10 [4,] 10 10 2 4 10 10 10 10 [5,] 10 10 10 10 1 3 10 10 [6,] 10 10 10 10 2 4 10 10 [7,] 10 10 10 10 10 10 1 3 [8,] 10 10 10 10 10 10 2 4 ... or integrate the loop into a.block.diag Peter Wolf> > > > best wishes > > Robin > > > > a.block.diag" <- function(a,b,pad=0) { > ## a.block.daig combines arrays a and b and builds a new array which > has > ## a and b as blocks on its diagonal. pw 10/2004 > > if( (length(a)==1) & (length(b)==1) ){return(diag(c(a,b)))} > if(length(a)==1){dim(a) <- rep(1,length(dim(b)))} > if(length(b)==1){dim(b) <- rep(1,length(dim(a)))} > > if(length(dim.a <- dim(a)) != length(dim.b <- dim(b))){ > stop("a and b must have identical number of dimensions") > } > > seq.new <- function(i){if(i==0){return(NULL)} else {return(1:i)}} > s <- array(pad, dim.a+dim.b) > s <- do.call("[<-",c(list(s),lapply(dim.a,seq.new),list(a))) > ind <- lapply(seq(dim.b),function(i)seq.new(dim.b[[i]])+dim.a[[i]]) > do.call("[<-",c(list(s),ind,list(b))) > } > > > a <- matrix(1:4,2,2) > >> a > > [,1] [,2] > [1,] 1 3 > [2,] 2 4 > >> b <- array(1e44,c(0,3)) >> b > > [,1] [,2] [,3] > >> a.block.diag(a,b) > > [,1] [,2] [,3] [,4] [,5] > [1,] 1 3 0 0 0 > [2,] 2 4 0 0 0 > >> >> > >>> >>> >>> Before I go any further, I need to check that it's OK with you for >>> me to put this >>> function (or modifications of it) into my package, which is GPL. >>> Full authorship credit given. >>> Is this OK? >> >> >> Hallo Robin, >> >> you are welcome to put a.block.diag() in your package. >> The function is free software; you can redistribute it and/or modify >> it under the terms of the GNU... >> >> Peter Wolf >> >>>> >>>> >