[R] Using apply with more than one matrix

William Dunlap wdunlap at tibco.com
Thu May 1 18:45:38 CEST 2014


> Thank you, A.K.  I learned from both of your solutions.  I find the one that uses alply easier to follow and understand intuitively.

Another approach is to only loop over the 3rd dimensional slices of
a1.  Your original code, converted to a function so it is easier to
test and think about is

f0 <- function(a1, m1, m2) {
    stopifnot(length(dim(a1))==3, length(dim(m1))==2, length(dim(m2))==2,
              all(dim(m1)==dim(m2)), all(dim(m1)==dim(a1)[1:2]),
              dim(a1)[3] > 1)
    condition1 <- !is.na(m1)& m1 > m2
    ans <- array(NA_real_, dim(m1)) # initialize
    for(i in seq_len(dim(m1)[1])) {
        for(j in seq_len(dim(m1)[2])) {
            ind.not.na <- which(!is.na(a1[i,j,]))
            if(condition1[i,j] && length(ind.not.na) > 0) {
                ans[i,j] <- a1[i,j,ind.not.na[1]] + m2[i,j]
            }
        }
    }
    ans
}

I believe the following is an equivalent function, but I haven't
checked it much.
Note that the only loop is over the third dimension of a1, which I assume is not
too big.  That updating loop could probably be replaced by a call to Reduce.

f1 <- function(a1, m1, m2) {
    stopifnot(length(dim(a1))==3, length(dim(m1))==2, length(dim(m2))==2,
              all(dim(m1)==dim(m2)), all(dim(m1)==dim(a1)[1:2]),
              dim(a1)[3] > 1)

    firstGoodInA1 <- array(NA_real_, dim(a1)[1:2])
    for(k in seq_len(dim(a1)[3])) {
        isStillBad <- is.na(firstGoodInA1)
        firstGoodInA1[isStillBad] <- a1[, ,k][isStillBad]
    }
    condition1 <- !is.na(firstGoodInA1) & !is.na(m1) & (m1 > m2)
    ans <- array(NA_real_, dim(m1))
    ans[condition1] <- firstGoodInA1[condition1] + m2[condition1]
    ans
}


Bill Dunlap
TIBCO Software
wdunlap tibco.com


On Thu, May 1, 2014 at 8:30 AM, Waichler, Scott R
<Scott.Waichler at pnnl.gov> wrote:
> Thank you, A.K.  I learned from both of your solutions.  I find the one that uses alply easier to follow and understand intuitively.  I guess I'll want to learn more about what plyr can do.  I've been using R for years but hadn't pushed vectorization far enough in my code.  Now my computing will be more efficient.
>
> Thanks also to David and the others who responded to my question.  It all helped.  --Scott Waichler
>
>> -----Original Message-----
>> From: arun [mailto:smartpink111 at yahoo.com]
>> Sent: Thursday, May 01, 2014 6:24 AM
>> To: R. Help
>> Cc: Waichler, Scott R
>> Subject: Re: [R] Using apply with more than one matrix
>>
>> Hi,
>> Sorry, a typo in the previous function:
>> ---------------------------------
>> if (condition1[i] && !is.na(indx3)) {
>>             arr[x1][indx3] + m2[i]  ###### should be mat2[i]
>>         } else NA
>>
>> -----------
>>
>> Also, you can try:
>> library(plyr)
>> evaluateNew <- function(arr, mat1, mat2){ if (!all(dim(mat1) ==
>> dim(mat2))) {
>>         stop("both matrices should have equal dimensions")
>>     }
>>  indx1 <- unlist(alply(arr, c(1,2), function(x) {x[!is.na(x)][1]}))
>> condition1 <- !is.na(mat1) & mat1 > mat2 ifelse(condition1, indx1+mat2,
>> NA) } evaluateNew(a1, m1, m2)
>>
>> #    [,1] [,2]
>> #[1,]   NA 1.66
>> #[2,] 3.14   NA
>> A.K.
>>
>>
>>
>>
>> On Thursday, May 1, 2014 5:49 AM, arun <smartpink111 at yahoo.com> wrote:
>>
>>
>> Hi,
>>
>> You may try:
>> evaluate <- function(arr, mat1, mat2) {
>>     if (!all(dim(mat1) == dim(mat2))) {
>>         stop("both matrices should have equal dimensions")
>>     }
>>     indx1 <- as.matrix(do.call(expand.grid, lapply(dim(arr), sequence)))
>>     indx2 <- paste0(indx1[, 1], indx1[, 2])
>>     condition1 <- !is.na(mat1) & mat1 > mat2
>>     ans <- sapply(seq_along(unique(indx2)), function(i) {
>>         x1 <- indx1[indx2 %in% unique(indx2)[i], ]
>>         indx3 <- which(!is.na(arr[x1]))[1]
>>         if (condition1[i] && !is.na(indx3)) {
>>             arr[x1][indx3] + m2[i]
>>         } else NA
>>     })
>>     dim(ans) <- dim(mat1)
>>     ans
>> }
>> evaluate(a1,m1,m2)
>> #     [,1] [,2]
>> #[1,]   NA 1.66
>> #[2,] 3.14   NA
>>
>> A.K.
>>
>>
>>
>> On Thursday, May 1, 2014 2:53 AM, "Waichler, Scott R"
>> <Scott.Waichler at pnnl.gov> wrote:
>> > I would ask you to look at this loop-free approach and ask if this is
>> > not equally valid?
>> >
>> > ans <- matrix(NA, ncol=2, nrow=2)
>> > ind.not.na <- which(!is.na(a1))
>> > ans[] <- condition1*a1[,,ind.not.na[1]]+ m2  # two matrices of equal
>> >dimensions, one logical.
>> >  ans
>> >      [,1] [,2]
>> > [1,]   NA 1.66
>> > [2,] 2.74   NA
>>
>> Thanks, I am learning something.  I didn't know you could multiply a
>> logical object by a numerical one.  But notice the answer is not the same
>> as mine, because I am doing an operation on the vector of values a1[i,j,]
>> first.
>>
>> I tried a modification on sapply below, but it doesn't work  because I
>> haven't referenced the 3d array a1 properly.  So I guess I must try to get
>> a 2d result from a1 first, then use that in matrix arithmetic.
>>
>> Sapply or mapply may work, I haven't used these much and will try to learn
>> better how to use them.  Your use of sapply looks good; but I'm trying to
>> understand if and how I can bring in the operation on a1.  This doesn't
>> work:
>>
>> evaluate <- function(idx) {
>>   ind.not.na <- which(!is.na(a1[idx,])) ]))  # doesn't work; improper
>> indexing for a1
>>   if(length(ind.not.na) > 0) {
>>     return(condition1*(a1[idx,ind.not.na[1]] + m2[idx]))  # doesn't work;
>> improper indexing for a1
>>   }
>> }
>> vec <- sapply(seq(length(m2)), evaluate)
>>
>> Scott Waichler
>>
>>
>> > -----Original Message-----
>> > From: David Winsemius [mailto:dwinsemius at comcast.net]
>> > Sent: Wednesday, April 30, 2014 8:46 PM
>> > To: Waichler, Scott R
>> > Cc: Bert Gunter; r-help at r-project.org
>> > Subject: Re: [R] Using apply with more than one matrix
>> >
>> >
>> > On Apr 30, 2014, at 6:03 PM, Waichler, Scott R wrote:
>> >
>> > > Here is a working example with no random parts.  Thanks for your
>> > patience and if I'm still off the mark with my presentation I'll drop
>> > the matter.
>> > >
>> > > v <- c(NA, 1.5, NA, NA,
>> > >       NA, 1.1, 0.5, NA,
>> > >       NA, 1.3, 0.4, 0.9)
>> > > a1 <- array(v, dim=c(2,2,3))
>> > > m1 <- matrix(c(NA, 1.5, 2.1, NA), ncol=2, byrow=T)
>> > > m2 <- matrix(c(1.56, 1.64, 1.16, 2.92), ncol=2)
>> > > condition1 <- !is.na(m1)& m1 > m2
>> > >
>> > > ans <- matrix(NA, ncol=2, nrow=2) # initialize for(i in 1:2) {
>> > >for(j  in 1:2) {
>> > >    ind.not.na <- which(!is.na(a1[i,j,]))
>> > >    if(condition1[i,j] && length(ind.not.na) > 0) ans[i,j] <-
>> > >a1[i,j,ind.not.na[1]] + m2[i,j]  } } ans
>> > >     [,1] [,2]
>> > > [1,]   NA 1.66
>> > > [2,] 3.14   NA
>> >
>> > I would ask you to look at this loop-free approach and ask if this is
>> > not equally valid?
>> >
>> > ans <- matrix(NA, ncol=2, nrow=2)
>> > ind.not.na <- which(!is.na(a1))
>> > ans[] <- condition1*a1[,,ind.not.na[1]]+ m2  # two matrices of equal
>> >dimensions, one logical.
>> >  ans
>> >      [,1] [,2]
>> > [1,]   NA 1.66
>> > [2,] 2.74   NA
>> > >
>> > > Let me try asking again in words.  If I have multiple matrices or
>> > > slices
>> > of 3d arrays that are the same dimension, is there a way to pass them
>> > all to apply, and have apply take care of looping through i,j?
>> >
>> > I don't think `apply` is the correct function for this. Either
>> > `mapply` or basic matrix operation seem more likely to deliver correct
>> results:
>> >
>> >
>> > > I understand that apply has just one input object x.  I want to work
>> > > on
>> > more than one array object at once using a custom function that has
>> > this
>> > characteristic:  in order to compute the answer at i,j I need a result
>> > from higher order array at the same i,j.
>> >
>> > If you want to iterate over matrix indices you can either use the
>> > vector version e.g. m2[3] or the matrix version, m2[2,1[.
>> >
>> > vec <- sapply(seq(length(m2) , function(idx) m2[idx]*condition1[idx] )
>> >
>> >
>> >
>> > > This is what I tried to demonstrate in my example above.
>> > >
>> > > Thanks,
>> > > Scott
>> >
>> > David Winsemius
>> > Alameda, CA, USA
>>
>> ______________________________________________
>> R-help at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-
>> guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.



More information about the R-help mailing list