[R] multiple dimensional diag()

[Peter Wolf]

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