Matrix indexing in a loop

classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

Matrix indexing in a loop

Mills, Jason

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

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: Matrix indexing in a loop

Pontarelli, Brett
Do you have to use a loop?  The following function should do what you want for the 1st order:

rook = function(Y) {
        rsub = function(Z) {
                X = matrix(0,nrow(Z),ncol(Z));
                X[1:(N-1),1:M] = X[1:(N-1),1:M] + Z[2:N,1:M];
                X[2:N,1:M] = X[2:N,1:M] + Z[1:(N-1),1:M];
                X[1:N,1:(M-1)] = X[1:N,1:(M-1)] + Z[1:N,2:M];
                X[1:N,2:M] = X[1:N,2:M] + Z[1:N,1:(M-1)];
                return(X);
        }
        return(rsub(Y)/rsub(matrix(1,nrow(Y),ncol(Y))));
}

I'm not sure I understand how the higher orders work.  For example, an interior element for the 1st order is always divided by 4.  Is an interior element for a 3rd order divided by 4 or 8 or something else?  Also, how are you implementing your 3D matrices?

--Brett


-----Original Message-----
From: [hidden email] [mailto:[hidden email]] On Behalf Of Mills, Jason
Sent: Friday, February 17, 2006 1:36 PM
To: [hidden email]
Subject: [R] 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]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: Matrix indexing in a loop

Pontarelli, Brett
In reply to this post by Mills, Jason
You're right I had N and M defined outside of the function and rook and rsub were picking up on that.  The following is a bit better and more cleaned up version with i-th order option:

rook = function(Y,i) {
        N = nrow(Y);
        M = ncol(Y);
        rsub = function(Z,i) {
                X = matrix(0,N,M);
                X[1:(N-i),] = X[1:(N-i),] + Z[(1+i):N,];
                X[(1+i):N,] = X[(1+i):N,] + Z[1:(N-i),];
                X[,1:(M-i)] = X[,1:(M-i)] + Z[,(1+i):M];
                X[,(1+i):M] = X[,(1+i):M] + Z[,1:(M-i)];
                return(X);
        }
        return(rsub(Y,i)/rsub(matrix(1,N,M),i));
}

Notice that the variable "i" can be passed any value even one that causes an error.

--Brett


-----Original Message-----
From: Mills, Jason [mailto:[hidden email]]
Sent: Friday, February 17, 2006 4:11 PM
To: Pontarelli, Brett
Subject: RE: [R] Matrix indexing in a loop

Hi Brett, thanks for the tip.

I tried this function on my sample matrix and got the error message "Error in rsub(fun) : object "N" not found".  Your code looks like it should work, so I must be doing some wrong.  I will continue to experiment.

As for the neighbor pattern, the convention follows the rules of chess.
I would consider a 2nd order rook case to include only 4 elements:

a[i-2,j]
a[i+2,j]
a[i,j-2]
a[i,j+2]

Even 3rd order Rook would still only include 4 elements.  As another example, a Queen neighborhood includes 8 elements, independent of the lag order.


I haven't started using 3D arrays yet.  I am attempting spatio-temporal analysis and I have thought about representing my landscape (two
dimensions) over time (the third dimension).  For now, I'm trying to get a handle on working in two dimensions.

Thanks.  

Jason






-----Original Message-----
From: Pontarelli, Brett [mailto:[hidden email]]
Sent: Friday, February 17, 2006 4:04 PM
To: Mills, Jason; [hidden email]
Subject: RE: [R] Matrix indexing in a loop

Do you have to use a loop?  The following function should do what you want for the 1st order:

rook = function(Y) {
        rsub = function(Z) {
                X = matrix(0,nrow(Z),ncol(Z));
                X[1:(N-1),1:M] = X[1:(N-1),1:M] + Z[2:N,1:M];
                X[2:N,1:M] = X[2:N,1:M] + Z[1:(N-1),1:M];
                X[1:N,1:(M-1)] = X[1:N,1:(M-1)] + Z[1:N,2:M];
                X[1:N,2:M] = X[1:N,2:M] + Z[1:N,1:(M-1)];
                return(X);
        }
        return(rsub(Y)/rsub(matrix(1,nrow(Y),ncol(Y))));
}

I'm not sure I understand how the higher orders work.  For example, an interior element for the 1st order is always divided by 4.  Is an interior element for a 3rd order divided by 4 or 8 or something else?
Also, how are you implementing your 3D matrices?

--Brett


-----Original Message-----
From: [hidden email]
[mailto:[hidden email]] On Behalf Of Mills, Jason
Sent: Friday, February 17, 2006 1:36 PM
To: [hidden email]
Subject: [R] 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]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: Matrix indexing in a loop

Gabor Grothendieck
In reply to this post by Mills, Jason
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 <[hidden email]> 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]]
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: Matrix indexing in a loop

Gabor Grothendieck
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 <[hidden email]> 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 <[hidden email]> 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]]
> >
> > ______________________________________________
> > [hidden email] mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
> >
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: 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 <[hidden email]> 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 <[hidden email]> 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 <[hidden email]> 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]]
> > >
> > > ______________________________________________
> > > [hidden email] mailing list
> > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
> > >
> >
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html