vectorizing loops

General
1

i am beginner with R and i am trying to avoid loops. This is my code

NII <-280
NJJ<-237
#CUTK  is a nonegative matrix 282x239
ADVDOM <- array(0, dim=c(282,239))
for(i in seq(1,NII + 2){
  for(j in seq(1,NJJ + 2)){
    if(CUTK[i,j] > 0){
      iplus <- trunc(15* CUTK[i,j] / 4)
      jplus <- trunc(15* CUTK[i,j] / 4)
     
      for(i1 in seq(i - iplus, i + iplus + 1)){
        for(j1 in seq(j - jplus, j+ jplus)){
          if(i1 <= NII + 1 && i1 > 1 && j1 <= NJJ + 1 && j1 > 1){
            ADVDOM[i1,j1] <- 1
          }
        }
      }
    }
  }
}

These are my attemps, but they don't work. I can't spot the mistakes.

#1
i <- 1 : (NII + 2)
j <- 1 : (NJJ + 2)
outer(i,j, Vectorize(function(i,j){
  if(CUTK[i,j] > 0){
    iplus <- trunc(15* CUTK[i,j] / 4)
    jplus <- trunc(15* CUTK[i,j] / 4)
    
    i1 <- (i - iplus) : (i + iplus + 1)
    j1 <- (j - jplus) : (j + jplus)
    
    outer(i1,j1,Vectorize(function(i1,j1){
      if(i1 <= NII + 1 && i1 > 1 && j1 <= NJJ + 1 && j1 > 1){
        ADVDOM[i1,j1] <- 1
      }
    }))
  }
}))
#2
sapply(1: NII + 2 , function(i){
  sapply(1: NJJ + 2, function(j){
    if(CUTK[i,j] > 0){
      iplus <- trunc(15* CUTK[i,j] / 4)
      jplus <- trunc(15* CUTK[i,j] / 4)
      
      sapply(i - iplus: i + iplus + 1 , function(i1){
        sapply(j - jplus : j + jplus, function(j1){
          if(i1 <= NII + 1 && i1 > 1 && j1 <= NJJ + 1 && j1 > 1){
            ADVDOM[i1,j1] <- 1
          }
        })
      })
    }
  })
})
#3

ij <- expand.grid(i=1:(NII + 2), j=1:(NJJ + 2))
invisible(apply(ij,1, function(ij){
  if(CUTK[ij[1],ij[2]] > 0){
    iplus <- trunc(PrognosticSubDomainFactor * CUTK[i,j] / DXK)
    jplus <- trunc(PrognosticSubDomainFactor * CUTK[i,j] / DYK)
    
    ijplus <- expand.grid(ip = (ij[1] - iplus):(ij[1] + iplus + 1), jp = (ij[2] - jplus):(ij[2] + jplus) )
    apply(ijplus, 1 , function(ijplus){
      if(ijplus[1] <= NII + 1 && ijplus[1] > 1 && ijplus[2] <= NJJ + 1 && ijplus[2] > 1){
        ADVDOM[ijplus[1],ijplus[2]] <- 1
      }
    })
  }
}))
2

Is this your own algorithm ? What is it about ?
I can't even run your for loop version as I don't have the objects you reference....

3

The only variable I can't write is CUTK, but you can use a random nonnegative matrix 282x239.
The for loop version works fine, but it takes so much time...

4

Yeah, it seems your algorithm with 4 nested loops inherently scale poorly. Vectorising has technical benefits around reducing memory allocation costs but that's more a linear cost saving and probably not noticeable when you are working on O4.

Can you not give any context to the purpose of your code. Perhaps a clever mathematician or someone that knows good heuristics about your application can find a better algorithm.

5

It's part of bigger project about the calculation of pollution in a real city.
CUTK is the height of buildings ( so with the first two loops, you are considering all the buildings). In the other two loops you are marking the building in a certain area (iplus*jplus), where this area is centred in the i,j-building.

6

I think you are generate large innerloops that you mostly throw away.
I would like to optimise this process for you, but im distracted by the fact that, im assuming you did it for example purposes made the core modification of ADVDOM to be ADVDOM[i1,j1] <- 1

Are you perhaps intending to be more sophisticated at that point , if so, do you have a method of debiasing the fact that ADVDOM[i1,j1] would be modified more or less by larger and smaller iplus and jplus coefficients ?

7

non.zeros <- which(CUTK > 0, arr.ind = TRUE)
iplus <- trunc(PrognosticSubDomainFactor * CUTK[CUTK > 0] / DXK)
jplus <- trunc(PrognosticSubDomainFactor * CUTK[CUTK > 0] / DYK)
i1 <- (non.zeros[,"row"] - iplus)
i2 <- (non.zeros[,"row"] + iplus + 1)
j1 <- (non.zeros[,"col"] - jplus)
j2 <- (non.zeros[,"col"] + jplus)
for(z in seq_along(i1)){
il <- pmin(pmax(1,i1[z]:i2[z]), NII + 1)
jl <- pmin(pmax(1,j1[z]:j2[z]), NJJ + 1)
ADVDOM[il ,jl] <- 1
}

8

seems like a good direction, but this does give different results than the first approach.
Though I can't comment on which of the two is accurate...



calc1 <- function(x,y,PrognosticSubDomainFactor,DXK,DYK){
  set.seed(42)
  #CUTK  is a nonegative matrix e.g. 282x239
  CUTK <- array(runif(n=(y+2)*(x+2),min = 0,max=10), dim=c(y+2,x+2))
  NII <-y
  NJJ<-x
  ADVDOM <- array(0, dim=c(y+2,x+2))
  for(i in seq(1,NII + 2)){
    for(j in seq(1,NJJ + 2)){
      if(CUTK[i,j] > 0){
        iplus <- trunc(PrognosticSubDomainFactor* CUTK[i,j] / DXK)
        jplus <- trunc(PrognosticSubDomainFactor* CUTK[i,j] / DYK)
        
        for(i1 in seq(i - iplus, i + iplus + 1)){
          for(j1 in seq(j - jplus, j+ jplus)){
            if(i1 <= NII + 1 && i1 > 1 && j1 <= NJJ + 1 && j1 > 1){
              if( ADVDOM[i1,j1] ==0){
              }
              ADVDOM[i1,j1] <- 1
            }
          }
        }
      }
    }
  }
  return(ADVDOM)
}


calc2 <- function(x,y,PrognosticSubDomainFactor,DXK,DYK){
  set.seed(42)
  #CUTK  is a nonegative matrix e.g. 282x239
  CUTK <- array(runif(n=(y+2)*(x+2),min = 0,max=10), dim=c(y+2,x+2))
  NII <-y
  NJJ<-x
  ADVDOM <- array(0, dim=c(y+2,x+2))
  
  non.zeros <- which(CUTK > 0, arr.ind = TRUE)
  iplus <- trunc(PrognosticSubDomainFactor * CUTK[CUTK > 0] / DXK)
  jplus <- trunc(PrognosticSubDomainFactor * CUTK[CUTK > 0] / DYK)
  i1 <- (non.zeros[,"row"] - iplus)
  i2 <- (non.zeros[,"row"] + iplus + 1)
  j1 <- (non.zeros[,"col"] - jplus)
  j2 <- (non.zeros[,"col"] + jplus)
  for(z in seq_along(i1)){
    il <- pmin(pmax(1,i1[z]:i2[z]), NII + 1)
    jl <- pmin(pmax(1,j1[z]:j2[z]), NJJ + 1)
    ADVDOM[il ,jl] <- 1
  }
  return(ADVDOM)
}


#dummy process
calc3<- function(x,y)
{
  ADVDOM <- array(1, dim=c(y+2,x+2))
  ADVDOM[,1] <- 0
  ADVDOM[1,] <- 0
  ADVDOM[,x+2] <- 0
  ADVDOM[x+2,] <- 0
  return(ADVDOM)
}

library(microbenchmark)


microbenchmark(
  res_1 <- calc1(5,5,15,4,4),
  res_2 <- calc2(5,5,15,4,4),
  res_3 <- calc3(5,5)
  , times = 1L)

identical(res_1,res_2)
identical(res_1,res_3)

microbenchmark(
  res_1 <- calc1(300,300,15,4,4),
  res_2 <- calc2(300,300,15,4,4),
  res_3 <- calc3(300,300)
  , times = 1L)

identical(res_1,res_2)
identical(res_1,res_3)

Open question for me if there are particular numbers of PrognosticSubDomainFactor , DXK, DYK, that would lead a difference from the uniform behaviour so far observed and so would differentiate res_1 from res_3

9

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.