Как управлять трехмерными массивами в агентных моделях с помощью R?

Question

Я строю агентную модель с помощью R, но у меня возникают проблемы с памятью, когда я пытаюсь работать с большими трехмерными массивами.В трехмерных массивах первое измерение соответствует времени (от 1 до 3650 дней), второе измерение определяет свойства отдельных лиц или ячеек ландшафта, а третье измерение представляет каждую отдельную ячейку или ячейку ландшафта.На каждом временном шаге (1 день) каждый трехмерный массив заполняется с использованием нескольких функций.В идеале я хотел бы использовать свою ПРО для больших ландшафтов (например, 90000 ячеек), содержащих большое количество людей (например, 720000).На самом деле это невозможно из-за проблем с памятью.

В настоящее время трехмерные массивы определяются при инициализации, поэтому данные сохраняются в массиве на каждом временном шаге.Однако, чтобы заполнить один трехмерный массив в t из модели, мне нужно хранить только данные в t - 1 и t - tf - 1, где tf - фиксированный параметр длительности (например, tf = 320 дней).Вот пример одной модельной функции, которая используется для заполнения одного трехмерного массива (есть 3 параметра продолжительности):

    library(ff)
    library(magrittr)
    library(dplyr)
    library(tidyr)
    library(gtools)

    ## Define parameters
    tf1 <- 288
    tf2 <- 292
    tf3 <- 150

    ## Define 3D array
    col_array <- c(letters[seq( from = 1, to = 9 )])
    s_array <- ff(-999, dim=c(3650, 9, 2500), dimnames=list(NULL, col_array, as.character(seq(1, 2500, 1))), 
                              filename="s_array.ffd", vmode="double", overwrite = t) ## 3th dimension = patch ID

    ## Define initial array
    initial_s_array <- matrix(sample.int(100, size = 2500*9, replace = TRUE), nrow = 2500, ncol = 9, dimnames=list(NULL, col_array))

   ## Loop over time 
    line <- 1
    for(t in 1:3650){
      print(t)
      s_array[t,c("a", "b", "c", "d", "e", "f", "g", "h", "i"),] <- func(v_t_1 = ifelse((t - 1) >= 1, list(s_array[(t - 1),,]), NA)[[1]], 
                                                                         v_t_tf1_1 = ifelse((t - tf1 - 1) >= 1, list(s_array[(t - tf1 - 1),,]), NA)[[1]], 
                                                                         v_t_tf2_1 = ifelse((t - tf2 - 1) >= 1, list(s_array[(t - tf2 - 1),,]), NA)[[1]], 
                                                                         v_t_tf3_1 = ifelse((t - tf3 - 1) >= 1, list(s_array[(t - tf3 - 1),,]), NA)[[1]], 
                                                                         v_t_0 = initial_s_array, columns_names = col_array)
      line <- line + 1
    }


    func <- function(v_t_1, v_t_tf1_1, v_t_tf2_1, v_t_tf3_1, v_t_0, columns_names){

      ## Data at t-1
      dt_t_1 <- (ifelse(!(all(is.na(v_t_1))), 
                        list(v_t_1 %>% 
                               as.data.frame.table(stringsAsFactors = FALSE) %>% 
                               dplyr::mutate_all(as.character)), NA))[[1]]

      ## Data at t-tf1-1
      dt_t_tf1_1 <- (ifelse(!(all(is.na(v_t_tf1_1))), 
                             list(v_t_tf1_1 %>% 
                                    as.data.frame.table(stringsAsFactors = FALSE) %>% 
                                    dplyr::mutate_all(as.character)), NA))[[1]]

      ## Data at t-tf2-1
      dt_t_tf2_1 <- (ifelse(!(all(is.na(v_t_tf2_1))), 
                             list(v_t_tf2_1 %>% 
                                    as.data.frame.table(stringsAsFactors = FALSE) %>% 
                                    dplyr::mutate_all(as.character)), NA))[[1]]

      ## Data at t-tf3-1
      dt_t_tf3_1 <- (ifelse(!(all(is.na(v_t_tf3_1))), 
                             list(v_t_tf3_1 %>% 
                                    as.data.frame.table(stringsAsFactors = FALSE) %>% 
                                    dplyr::mutate_all(as.character)), NA))[[1]]

      ## Format data at t-1
      dt_t_1_reshape <- (ifelse(!(all(is.na(dt_t_1))), 
                                list(dt_t_1 %>%
                                       dplyr::rename(ID = Var2) %>%
                                       tidyr::spread(Var1, Freq) %>%
                                       dplyr::select(ID, columns_names) %>%
                                       dplyr::arrange(match(ID, mixedsort(colnames(v_t_1))))), NA))[[1]]

      ## Format data at t-tf1-1
      dt_t_tf1_1_reshape <- (ifelse(!(all(is.na(dt_t_tf1_1))), 
                                     list(dt_t_tf1_1 %>%
                                            dplyr::rename(ID = Var2) %>%
                                            tidyr::spread(Var1, Freq) %>%
                                            dplyr::select(ID, columns_names) %>%
                                            dplyr::arrange(match(ID, mixedsort(colnames(v_t_tf1_1))))), NA))[[1]]

      ## Format data at t-tf2-1
      dt_t_tf2_1_reshape <- (ifelse(!(all(is.na(dt_t_tf2_1))), 
                                     list(dt_t_tf2_1 %>%
                                            dplyr::rename(ID = Var2) %>%
                                            tidyr::spread(Var1, Freq) %>%
                                            dplyr::select(ID, columns_names) %>%
                                            dplyr::arrange(match(ID, mixedsort(colnames(v_t_tf2_1))))), NA))[[1]]

      ## Format data at t-tf3-1
      dt_t_tf3_1_reshape <- (ifelse(!(all(is.na(dt_t_tf3_1))), 
                                     list(dt_t_tf3_1 %>%
                                            dplyr::rename(ID = Var2) %>%
                                            tidyr::spread(Var1, Freq) %>%
                                            dplyr::select(ID, columns_names) %>%
                                            dplyr::arrange(match(ID, mixedsort(colnames(v_t_tf3_1))))), NA))[[1]]

      ## Retrieve data
      a_t_1 <- (ifelse((!(all(is.na(dt_t_1_reshape)))), list(as.numeric(dt_t_1_reshape[,c("a")])), list(v_t_0[,c("a")])))[[1]] 
      d_t_1 <- (ifelse((!(all(is.na(dt_t_1_reshape)))), list(as.numeric(dt_t_1_reshape[,c("d")])), list(v_t_0[,c("d")])))[[1]]   
      g_t_1 <- (ifelse((!(all(is.na(dt_t_1_reshape)))), list(as.numeric(dt_t_1_reshape[,c("f")])), list(v_t_0[,c("f")])))[[1]] 

      a_t_tf1_1 <- (ifelse(!(all(is.na(dt_t_tf1_1_reshape))), list(as.numeric(dt_t_tf1_1_reshape[,c("a")])), 0))[[1]] 
      d_t_tf2_1 <- (ifelse(!(all(is.na(dt_t_tf2_1_reshape))), list(as.numeric(dt_t_tf2_1_reshape[,c("d")])), 0))[[1]] 
      g_t_tf3_1 <- (ifelse(!(all(is.na(dt_t_tf3_1_reshape))), list(as.numeric(dt_t_tf3_1_reshape[,c("f")])), 0))[[1]] 

      b_t_1 <- (ifelse((!(all(is.na(dt_t_1_reshape)))), list(as.numeric(dt_t_1_reshape[,c("b")])), list(v_t_0[,c("b")])))[[1]] 
      e_t_1 <- (ifelse((!(all(is.na(dt_t_1_reshape)))), list(as.numeric(dt_t_1_reshape[,c("e")])), list(v_t_0[,c("e")])))[[1]] 
      h_t_1 <- (ifelse((!(all(is.na(dt_t_1_reshape)))), list(as.numeric(dt_t_1_reshape[,c("h")])), list(v_t_0[,c("h")])))[[1]] 

      b_t_tf1_1 <- (ifelse(!(all(is.na(dt_t_tf1_1_reshape))), list(as.numeric(dt_t_tf1_1_reshape[,c("b")])), 0))[[1]] 
      e_t_tf2_1 <- (ifelse(!(all(is.na(dt_t_tf2_1_reshape))), list(as.numeric(dt_t_tf2_1_reshape[,c("e")])), 0))[[1]] 
      h_t_tf3_1 <- (ifelse(!(all(is.na(dt_t_tf3_1_reshape))), list(as.numeric(dt_t_tf3_1_reshape[,c("h")])), 0))[[1]] 

      ## Define discrete equations
      a_t <- round(0.4*a_t_1 + 0.5*a_t_tf1_1)
      b_t <- round(0.5*b_t_1 + 0.6*b_t_tf1_1)
      c_t <- a_t + b_t
      d_t <- round(0.7*d_t_1 + 0.7*d_t_tf2_1)
      e_t <- round(0.9*e_t_1 + 0.4*e_t_tf2_1)
      f_t <- d_t + e_t
      g_t <- round(0.3*g_t_1 + 0.2*g_t_tf3_1)
      h_t <- round(0.5*h_t_1 + 0.1*h_t_tf3_1)
      i_t <- g_t + h_t

      ## Update the values
      dt_array <- as.matrix(cbind(a_t, b_t, c_t, d_t, e_t, f_t, g_t, h_t, i_t))
      ## print(dt_array)

      ## Build the output matrix         
      dt_array <- t(dt_array)

      return(dt_array)

    }

Функция «func» принимает в качестве данных аргумента при t - 1 и t - tf -1, предоставляя один трехмерный массив «s_array».Функция возвращает фрейм данных, который используется для заполнения трехмерного массива.

Я думаю, что я мог бы уменьшить первое измерение моих массивов, сохраняя только данные в t - 1 и t - tf - 1 (вместохранить данные на каждом временном шаге от 1 до 3650 дней).Однако я не знаю, как управлять этими новыми трехмерными массивами в ПРО на каждом временном шаге (т. Е. Как инициализировать трехмерные массивы и хранить данные только в t - 1 и t - tf - 1)?

Редактировать: я протестировал пример с 90000 наблюдениями для 3-го измерения.Количество строк (т.е. 3650) в каждом массиве слишком велико.

> s_array <- ff(-999, dim=c(3650, 9, 90000), dimnames=list(NULL, col_array, as.character(seq(1, 90000, 1))), 
+               filename="s_array.ffd", vmode="double", overwrite = TRUE) 
Error in if (length < 0 || length > .Machine$integer.max) stop("length must be between 1 and .Machine$integer.max") : 
  missing value where TRUE/FALSE needed
In addition: Warning message:
In ff(-999, dim = c(3650, 9, 90000), dimnames = list(NULL, col_array,  :
  NAs introduced by coercion to integer range

Есть ли способ уменьшить количество строк и применить функцию, которая используется для заполнения массивов?

Answer 1

Причина, по которой я сказал, что R, вероятно, не идеален, заключается в его семантике «копировать при изменении», поэтому каждый раз, когда вы что-то изменяете в любом массиве / матрице / фрейме данных, должна быть сделана копия.Я думаю, что R может делать некоторые умные вещи с помощью управления памятью, но все же.

Использование ff, чтобы избежать одновременного хранения всех временных интервалов в ОЗУ, вероятно, действительно выгодно, но ваш код, вероятно, меняет памятьформатирует слишком часто, манипулируя структурами данных взад-впередЯ думаю Я упаковал логику в пару R6 классов (вместе с некоторыми улучшениями), возможно, это поможет вам начать улучшать использование памяти вашего кода:

suppressPackageStartupMessages({
  library(R6)
  library(ff)
})

SArray <- R6::R6Class(
  "SArray",
  public = list(
    time_slices = NULL,

    initialize = function(sdim, sdimnames) {
      self$time_slices <- lapply(1L:sdim[1L], function(ignored) {
        ff(NA_real_, vmode="double", dim=sdim[-1L], dimnames=sdimnames[-1L])#, FF_RETURN=FALSE)
      })

      names(self$time_slices) <- sdimnames[[1L]]
    }
  )
)

`[.SArray` <- function(s_array, i, j, ...) {
  s_array$time_slices[[i]][j,]
}

`[<-.SArray` <- function(s_array, i, j, ..., value) {
  s_array$time_slices[[i]][j,] <- value
  s_array
}

dim.SArray <- function(x) {
  c(length(x$time_slices), dim(x$time_slices[[1L]]))
}

ABM <- R6::R6Class(
  "ABM",
  public = list(
    s_array = NULL,
    tf1 = NULL,
    tf2 = NULL,
    tf3 = NULL,

    initialize = function(sdim, sdimnames, tfs) {
      self$tf1 <- tfs[1L]
      self$tf2 <- tfs[2L]
      self$tf3 <- tfs[3L]

      self$s_array <- SArray$new(sdim, sdimnames)
    },

    init_abm = function(seed = NULL) {
      set.seed(seed)
      sdim <- dim(self$s_array)
      s_init <- matrix(sample.int(100L, size = 6L * sdim[3L], replace=TRUE),
                       nrow=6L, ncol=sdim[3L],
                       dimnames=list(c("a", "b", "d", "e", "g", "h"),
                                     as.character(1:sdim[3L])))

      self$a(1L, s_init["a", ])
      self$b(1L, s_init["b", ])
      self$c(1L)
      self$d(1L, s_init["d", ])
      self$e(1L, s_init["e", ])
      self$f(1L)
      self$g(1L, s_init["g", ])
      self$h(1L, s_init["h", ])
      self$i(1L)

      private$t <- 1L

      invisible()
    },

    can_advance = function() {
      private$t < dim(self$s_array)[1L]
    },

    advance = function(verbose = FALSE) {
      t <- private$t + 1L
      if (verbose) print(t)

      self$a(t)
      self$b(t)
      self$c(t)
      self$d(t)
      self$e(t)
      self$f(t)
      self$g(t)
      self$h(t)
      self$i(t)

      private$t <- t

      invisible()
    },

    # get time slice at t - tf - 1 for given letter
    s_tf = function(t, tf, letter) {
      t_tf_1 <- t - tf - 1L

      if (t_tf_1 > 0L)
        self$s_array[t_tf_1, letter, ]
      else
        0
    },

    # discrete equations

    a = function(t, t_0) {
      if (t < 2L) {
        t <- 1L
        t_1 <- t_0
        t_tf_1 <- 0
      }
      else {
        t_1 <- self$s_array[t - 1L, "a", ]
        t_tf_1 <- self$s_tf(t, self$tf1, "a")
      }

      self$s_array[t, "a", ] <- round(0.4 * t_1 + 0.5 * t_tf_1)
      invisible()
    },

    b = function(t, t_0) {
      if (t < 2L) {
        t <- 1L
        t_1 <- t_0
        t_tf_1 <- 0
      }
      else {
        t_1 <- self$s_array[t - 1L, "b", ]
        t_tf_1 <- self$s_tf(t, self$tf1, "b")
      }

      self$s_array[t, "b", ] <- round(0.5 * t_1 + 0.6 * t_tf_1)
      invisible()
    },

    c = function(t) {
      if (t < 1L) stop("t must be positive")
      a_t <- self$s_array[t, "a", ]
      b_t <- self$s_array[t, "b", ]
      self$s_array[t, "c", ] <- a_t + b_t
      invisible()
    },

    d = function(t, t_0) {
      if (t < 2L) {
        t <- 1L
        t_1 <- t_0
        t_tf_1 <- 0
      }
      else {
        t_1 <- self$s_array[t - 1L, "d", ]
        t_tf_1 <- self$s_tf(t, self$tf2, "d")
      }

      self$s_array[t, "d", ] <- round(0.7 * t_1 + 0.7 * t_tf_1)
      invisible()
    },

    e = function(t, t_0) {
      if (t < 2L) {
        t <- 1L
        t_1 <- t_0
        t_tf_1 <- 0
      }
      else {
        t_1 <- self$s_array[t - 1L, "e", ]
        t_tf_1 <- self$s_tf(t, self$tf2, "e")
      }

      self$s_array[t, "e", ] <- round(0.9 * t_1 + 0.4 * t_tf_1)
      invisible()
    },

    f = function(t) {
      if (t < 1L) stop("t must be positive")
      d_t <- self$s_array[t, "d", ]
      e_t <- self$s_array[t, "e", ]
      self$s_array[t, "f", ] <- d_t + e_t
      invisible()
    },

    g = function(t, t_0) {
      if (t < 2L) {
        t <- 1L
        t_1 <- t_0
        t_tf_1 <- 0
      }
      else {
        t_1 <- self$s_array[t - 1L, "g", ]
        t_tf_1 <- self$s_tf(t, self$tf3, "g")
      }

      self$s_array[t, "g", ] <- round(0.3 * t_1 + 0.2 * t_tf_1)
      invisible()
    },

    h = function(t, t_0) {
      if (t < 2L) {
        t <- 1L
        t_1 <- t_0
        t_tf_1 <- 0
      }
      else {
        t_1 <- self$s_array[t - 1L, "h", ]
        t_tf_1 <- self$s_tf(t, self$tf3, "h")
      }

      self$s_array[t, "h", ] <- round(0.5 * t_1 + 0.1 * t_tf_1)
      invisible()
    },

    i = function(t) {
      if (t < 1L) stop("t must be positive")
      g_t <- self$s_array[t, "g", ]
      h_t <- self$s_array[t, "h", ]
      self$s_array[t, "i", ] <- g_t + h_t
      invisible()
    }
  ),
  private = list(
    t = NULL
  )
)

max_t <- 10
abm <- ABM$new(c(max_t, 9, 2500),
               list(NULL, letters[1:9], as.character(1:2500)),
               c(288L, 292L, 150L))

abm$init_abm()

while (abm$can_advance()) {
  abm$advance(TRUE)
}

anyNA(abm$s_array[])
# FALSE

Некоторые функции в рамках дискретных уравнений инкапсулируют логику для инициализации, когда t < 2L.Класс SArray разделяет трехмерный массив на список двумерных массивов, чтобы обойти ограничение .Machine$integer.max.