Некоторые входные данные из матриц, а другие из векторов, как заменить 2d цикл?

Question

У меня есть функция, которую я бы предпочел не редактировать.Некоторые входы зависят от времени (показаны здесь как векторы), некоторые зависят от времени, а также зависят от другой переменной, Nj.

В настоящее время я выполняю цикл по времени (Ni) и цикл по Nj и вычисляю каждое значение отдельно.Насколько я вижу, семейство функций apply работает только в этом сценарии, когда все входные данные имеют одинаковую размерность.Есть ли другой способ, которым я могу сделать это?

Ni <- 10
Nj <- 10

a <- matrix(1:100/100, Ni, Nj)
b <- matrix(runif(100)*500, Ni, Nj)
c <- runif(Ni)
d <- c + runif(Ni)
e <- runif(1)*100
f <- c(0.3, 0.7)

funky <- function(a, b, c, d, e, f) {

  firstLine <- a / b
  secondLine <- firstLine * c
  thirdLine <- (secondLine + 45) / d
  fourthLine <- thirdLine + e
  result <- c(f[1] * fourthLine, f[2] * fourthLine)

  result

}

resultMatrix1 <- matrix(numeric(), Ni, Nj)
resultMatrix2 <- matrix(numeric(), Ni, Nj)

for (i in 1:Ni) {
  for (j in 1:Nj) {

    result <- funky(a[i, j],
                    b[i, j],
                    c[i],
                    d[i],
                    e, 
                    f
    )

    resultMatrix1[i, j] <- result[1]
    resultMatrix2[i, j] <- result[2]
  }
}

Это некоторый составленный код, который я только что скомбинировал, который показывает, что я имею в виду относительно входных измерений.Проблема в том, что используемая мной функция не очень быстрая, а фактическая таблица результатов, которую я заполняю, составляет около 100 * 150, а запуск занимает около получаса.

Answer 1

Вы можете векторизовать один из циклов. Я буду векторизовать цикл for на i.

set.seed(1234)

Ni <- 10
Nj <- 10

a <- matrix(1:100/100, Ni, Nj)
b <- matrix(runif(100)*500, Ni, Nj)
c <- runif(Ni)
d <- c + runif(Ni)
e <- runif(1)*100
f <- c(0.3, 0.7)

funky <- function(a, b, c, d, e, f) {

  firstLine <- a / b
  secondLine <- firstLine * c
  thirdLine <- (secondLine + 45) / d
  fourthLine <- thirdLine + e
  result <- c(f[1] * fourthLine, f[2] * fourthLine)

  result

}

resultMatrix1 <- matrix(numeric(), Ni, Nj)
resultMatrix2 <- matrix(numeric(), Ni, Nj)

for (i in 1:Ni) {
  for (j in 1:Nj) {

    result <- funky(a[i, j],
                    b[i, j],
                    c[i],
                    d[i],
                    e, 
                    f
    )

    resultMatrix1[i, j] <- result[1]
    resultMatrix2[i, j] <- result[2]
  }
}

resultMatrix1
resultMatrix2


funky2 <- function(a, b, c, d, e, f) {

  firstLine <- a / b
  secondLine <- firstLine * c
  thirdLine <- (secondLine + 45) / d
  fourthLine <- thirdLine + e
  result <- matrix(c(f[1] * fourthLine, f[2] * fourthLine), ncol = 2)

  result
}

rmat <- array(NA, dim = c(2, Ni, Nj))
for(j in 1:Nj) {
  result <- funky2(a[, j], b[, j], c, d, e, f)
  rmat[1, , j] <- result[, 1]
  rmat[2, , j] <- result[, 2]
}


identical(resultMatrix1, rmat[1, , ])
#[1] TRUE
identical(resultMatrix2, rmat[2, , ])
#[1] TRUE


JoePye <- function(a, b, c, d, e, f, Ni, Nj){
  resultMatrix1 <- matrix(numeric(), Ni, Nj)
  resultMatrix2 <- matrix(numeric(), Ni, Nj)

  for (i in 1:Ni) {
    for (j in 1:Nj) {

      result <- funky(a[i, j],
                      b[i, j],
                      c[i],
                      d[i],
                      e, 
                      f
      )

      resultMatrix1[i, j] <- result[1]
      resultMatrix2[i, j] <- result[2]
    }
  }

  list(Mat1 = resultMatrix1,
       Mat2 = resultMatrix2)
}
RuiB <- function(a, b, c, d, e, f, Ni, Nj){
  rmat <- array(NA, dim = c(2, Ni, Nj))
  for(j in 1:Nj) {
    result <- funky2(a[, j], b[, j], c, d, e, f)
    rmat[1, , j] <- result[, 1]
    rmat[2, , j] <- result[, 2]
  }
  rmat
}

library(microbenchmark)

mb10 <- microbenchmark(
  JoePye = JoePye(a, b, c, d, e, f, Ni, Nj),
  RuiB = RuiB(a, b, c, d, e, f, Ni, Nj)
)
mb10
#Unit: microseconds
#   expr     min       lq     mean   median       uq     max neval cld
# JoePye 473.950 479.2840 496.5165 483.6325 495.7015 758.162   100   b
#   RuiB 150.502 157.8205 174.8730 160.1100 165.2400 647.653   100  a 

Это ускорение в 3 раза. Теперь протестируйте с большим вводом.

set.seed(1234)

Ni <- 1e2
Nj <- 2e1

a <- matrix(seq.int(Ni*Nj)/100, Ni, Nj)
b <- matrix(runif(Ni*Nj)*500, Ni, Nj)
c <- runif(Ni)
d <- c + runif(Ni)
e <- runif(1)*100
f <- c(0.3, 0.7)

res1 <- JoePye(a, b, c, d, e, f, Ni, Nj)
res2 <- RuiB(a, b, c, d, e, f, Ni, Nj)

identical(res1[[1]], res2[1, , ])
#[1] TRUE
identical(res1[[2]], res2[2, , ])
#[1] TRUE


mb100 <- microbenchmark(
  JoePye = JoePye(a, b, c, d, e, f, Ni, Nj),
  RuiB = RuiB(a, b, c, d, e, f, Ni, Nj),
  times = 10
)
mb100
#Unit: microseconds
#   expr      min       lq      mean   median       uq       max neval cld
# JoePye 9198.846 9248.114 9421.2359 9352.244 9426.161 10147.642    10   b
#   RuiB  478.564  490.404  533.8198  522.573  594.841   602.938    10  a 

И при большем входе коэффициент ускорения увеличился до 18.