[R] Matrix indexing in a loop

[Gabor Grothendieck]

There was an error in the function f -- it only worked correctly if
order was 1.  Here it is with that fixed.  The function f is the
only thing changed from my last post.  It makes use of the
fact that sum(x) == n and sum(x*x) == n*n  only occur when
one element of x is n and the rest are 0.

root3 <- function(x, order = 1) {
	f <- function(x,y) {
		z <- abs(as.numeric(x) - as.numeric(y))
		(sum(z) == order) & (sum(z*z) == order * 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:
> 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]]
> > >
> > > ______________________________________________
> > > R-help at stat.math.ethz.ch mailing list
> > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
> > >
> >
>