Re: [R] Matrix indexing in a loop

From: Gabor Grothendieck <ggrothendieck_at_gmail.com>
Date: Sun 19 Feb 2006 - 01:08:36 EST

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



R-help@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 Received on Sun Feb 19 01:12:17 2006

This archive was generated by hypermail 2.1.8 : Mon 20 Feb 2006 - 14:08:44 EST