> 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.waich...@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:smartpink...@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 <smartpink...@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.waich...@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:dwinsem...@comcast.net] >> > Sent: Wednesday, April 30, 2014 8:46 PM >> > To: Waichler, Scott R >> > Cc: Bert Gunter; r-help@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@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@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@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.