Matrix indexing in a loop

How do you specify matrix location a[i,j] (or a[i-1,j], etc.) in a "for"
loop?

I am looking for a flexible method of indexing neighbors over a series
of lags (1,2,3...) and I may wish to extend this method to 3D arrays.


Example:

Data matrix

fun

    [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12


For each element a[i,j] in "fun", sum the 1st order (Rook's) neighbors:

a[i-1,j]

a[i+1,j]

a[i,j-1]

a[i,j+1]

Then divide by the number of elements included as neighbors-- this
number depends on the location of a[i,j] in the matrix.


Insert the product of the neighbor calculation for each a[i,j] into the
corresponding position b[i,j] in an empty matrix with the same
dimensions as "fun".


For example, element [2,2] in "fun" should yield element [2,2] in a new
matrix equal to 24/4=6.  Of course, element [1,1] in the new matrix
should be the product of only two numbers.


Thanks

J. Mills

       [[alternative HTML version deleted]]

For 2d here is a solution based on zoo.  It turns the matrix into
a time series and lags it forwards and backwards and does the
same for its transpose in order to avoid index machinations.
The function is called rook2 and it first defines three local
functions, one that converts NAs to zero, one that does
a lag using na.pad = TRUE and one to invoke Lag and
add up the lags:

library(zoo)
rook2 <- function(x, i = 1) {
  na2zero <- function(x) ifelse(is.na(x), 0, x)
  Lag <- function(x, i) na2zero(lag(zoo(x), i, na.pad = TRUE))
  Rook <- function(x, i) Lag(x, i) + Lag(x, -i) + t(Lag(t(x), i) +
Lag(t(x), -i))
  Rook(x, i) / Rook(1+0*x, i)
}

# test
m <- matrix(1:24, 6)
rook2(m)


On 2/17/06, Mills, Jason <Jason.Mills at afhe.ualberta.ca> wrote:

How do you specify matrix location a[i,j] (or a[i-1,j], etc.) in a "for"
loop?

I am looking for a flexible method of indexing neighbors over a series
of lags (1,2,3...) and I may wish to extend this method to 3D arrays.


Example:

Data matrix

fun

    [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12


For each element a[i,j] in "fun", sum the 1st order (Rook's) neighbors:

a[i-1,j]

a[i+1,j]

a[i,j-1]

a[i,j+1]

Then divide by the number of elements included as neighbors-- this
number depends on the location of a[i,j] in the matrix.


Insert the product of the neighbor calculation for each a[i,j] into the
corresponding position b[i,j] in an empty matrix with the same
dimensions as "fun".


For example, element [2,2] in "fun" should yield element [2,2] in a new
matrix equal to 24/4=6.  Of course, element [1,1] in the new matrix
should be the product of only two numbers.


Thanks

J. Mills

       [[alternative HTML version deleted]]

Thought about this some more and here is a solution that
works with 2d and 3d (and higher dimensions).

inner is a generalized inner product similar to a function
I have posted previously and f(x,y) is an inner product such
that f(x,y) is TRUE if abs(x - y) == order (after converting
both x and y to numeric) and FALSE otherwise.  Root then
uses as.data.frame.table to turn the array into data frames
with the last column having
the data and the other columns representing the indexes.
We then perform the inner product and use the resulting
matrix to multiply c(x) which is the input strung out
into a vector reshaping it back into the same shape as x.
Finally we divide each cell by the number of surrounding
cells.

root3 <- function(x, order = 1) {
       f <- function(x,y) sum(abs(as.numeric(x) - as.numeric(y))) == order
       inner <- function(a,b,f)
                       apply(b,1,function(x)apply(a,1,function(y)f(x,y)))

       Root <- function(x) {
               n <- length(dim(x))
               dd <- sapply(as.data.frame.table(x)[,-n-1], as.numeric)
               structure(inner(dd, dd, f) %*% c(x), .Dim = dim(x))
       }
       Root(x) / Root(1+0*x)
}

# tests
m <- matrix(1:24, 6) # 2d
root3(m)
mm <- array(1:24, 2:4)  # 3d
root3(mm)


On 2/17/06, Gabor Grothendieck <ggrothendieck at gmail.com> wrote:

For 2d here is a solution based on zoo.  It turns the matrix into
a time series and lags it forwards and backwards and does the
same for its transpose in order to avoid index machinations.
The function is called rook2 and it first defines three local
functions, one that converts NAs to zero, one that does
a lag using na.pad = TRUE and one to invoke Lag and
add up the lags:

library(zoo)
rook2 <- function(x, i = 1) {
  na2zero <- function(x) ifelse(is.na(x), 0, x)
  Lag <- function(x, i) na2zero(lag(zoo(x), i, na.pad = TRUE))
  Rook <- function(x, i) Lag(x, i) + Lag(x, -i) + t(Lag(t(x), i) +
Lag(t(x), -i))
  Rook(x, i) / Rook(1+0*x, i)
}

# test
m <- matrix(1:24, 6)
rook2(m)


On 2/17/06, Mills, Jason <Jason.Mills at afhe.ualberta.ca> wrote:

How do you specify matrix location a[i,j] (or a[i-1,j], etc.) in a "for"
loop?

I am looking for a flexible method of indexing neighbors over a series
of lags (1,2,3...) and I may wish to extend this method to 3D arrays.


Example:

Data matrix

fun

    [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12


For each element a[i,j] in "fun", sum the 1st order (Rook's) neighbors:

a[i-1,j]

a[i+1,j]

a[i,j-1]

a[i,j+1]

Then divide by the number of elements included as neighbors-- this
number depends on the location of a[i,j] in the matrix.


Insert the product of the neighbor calculation for each a[i,j] into the
corresponding position b[i,j] in an empty matrix with the same
dimensions as "fun".


For example, element [2,2] in "fun" should yield element [2,2] in a new
matrix equal to 24/4=6.  Of course, element [1,1] in the new matrix
should be the product of only two numbers.


Thanks

J. Mills

       [[alternative HTML version deleted]]