Как векторизовать цикл for

Question

Я начинающий с R и стараюсь избегать циклов. Это мой код:

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

Это мои попытки, но они не работают. Я не могу определить ошибки.

#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
      }
    })
  }
}))

Я хотел бы понять эту технику векторизации в R. У меня есть несколько второстепенных вопросов:

1. как можно избежать печать после решения № 1 и № 2? Я пробовал invisible() в # 3, но я не знаю, есть ли более простое решение.

2. если бы я использовал mapply?

Answer 1

Следующий код намного проще, имеет только один for l oop и дает тот же результат. Он начинается с создания индексной матрицы, сообщающей положительное значение CUTK, а затем использует эту матрицу в остальной части кода.

ADVDOM2 <- array(0, dim = c(282, 239))

non.zeros <- which(CUTK > 0, arr.ind = TRUE)
iplus <- trunc(15*CUTK[non.zeros]/4)
jplus <- trunc(15*CUTK[non.zeros]/4)
ip1 <- pmax(2, non.zeros[, 'row'] - iplus)
ip2 <- pmin(NII + 1, non.zeros[, 'row'] + iplus + 1)
jp1 <- pmax(2, non.zeros[, 'col'] - jplus)
jp2 <- pmin(NJJ + 1, non.zeros[, 'col'] + jplus)
for(k in seq_along(iplus)){
  i <- ip1[k]:ip2[k]
  j <- jp1[k]:jp2[k]
  ADVDOM2[i, j] <- 1
}

identical(ADVDOM, ADVDOM2)
#[1] TRUE

Answer 2

В отсутствие CUTK и в интересах указать вам только правильное направление, я просто хочу быстро сказать: вам не нужны циклы для доступа к каждому элементу матрицы. Например, посмотрите, как R работает для автоматического доступа к каждому элементу:

a <- matrix(c(1:9),nrow=3, ncol=3) 
#     [,1] [,2] [,3]
#[1,]    1    4    7
#[2,]    2    5    8
#[3,]    3    6    9

ifelse(a<5,0,1)      # assign to a named var if no output wanted
#    [,1] [,2] [,3]
#[1,]    0    0    1
#[2,]    0    1    1
#[3,]    0    1    1

Поскольку вы обращаетесь ко второй матрице и можете подумать, что это проблема, это не так. R заботится об индексации (убедитесь, что ранги одинаковые или кратные):

b <- t(a)           # t() is the transpose function
ifelse(a>b,1,0)     # assign to a named var if no output wanted
#     [,1] [,2] [,3]
#[1,]    0    1    1
#[2,]    0    0    1
#[3,]    0    0    0