[R] multiple dimensional diag()
Peter Wolf
s-plus at wiwi.uni-bielefeld.de
Tue Oct 5 19:14:06 CEST 2004
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
>>
>>>>
>>>>
>
More information about the R-help
mailing list