update parallel example

DESCRIPTION

@@ -42,6 +42,7 @@ Imports:
     stats,
     utils,
     lubridate,
+    magrittr,
     Rcpp,
     foreach
 Suggests:

NAMESPACE

@@ -11,6 +11,7 @@ export(vic)
 export(vic_param)
 export(vic_params)
 export(vic_version)
+import(magrittr)
 importFrom(Rcpp,sourceCpp)
 importFrom(lubridate,day)
 importFrom(lubridate,month)

R/tools.R

@@ -12,3 +12,17 @@ listk <- function(...) {
   x <- setNames(list(...), vars)
   return(x)
 }
+
+
+`%||%` <- function(x, y) {
+  if (is.null(x)) {
+    y
+  } else {
+    x
+  }
+}
+
+set_dimnames <- function(x, value) {
+  dimnames(x) <- value
+  x
+}

R/vic.R

@@ -248,149 +248,152 @@ vic <- function(forcing, soil, veg,
                 output_info = NULL,
                 veglib = NULL, snowband = NULL, lake = NULL,
                 forcing_veg = NULL,
-                veg_force_types = c('albedo', 'LAI', 'fcanopy'),
+                veg_force_types = c("albedo", "LAI", "fcanopy"),
                 parall = FALSE) {
-
   tp1 <- proc.time()
 
-  if(!is.list(output_info) || length(output_info) == 0) {
-    output_info <- list(
-      fluxes = list(timescale = 'day', aggpar = 1,
-                    outvars = c('OUT_PREC', 'OUT_EVAP', 'OUT_RUNOFF',
-                                'OUT_BASEFLOW', 'OUT_WDEW', 'OUT_SOIL_MOIST')
-                    ),
-      snow = list(timescale = 'day', aggpar = 1,
-                    outvars = c('OUT_SWE', 'OUT_SNOW_DEPTH', 'OUT_SNOWF',
-                                'OUT_SNOW_MELT', 'OUT_SNOW_SURF_TEMP',
-                                'OUT_SNOW_PACK_TEMP')
-                    )
-    )
+  if (!is.list(output_info) || length(output_info) == 0) {
+    output_info <- default_output_info
   }
-  if(!is.list(output_info[[1]]))
-    output_info <- list(output = output_info)
 
+  if (!is.list(output_info[[1]])) output_info %<>% list(output = .)
+
   output_info <- deal_output_info(output_info)
   n_outputs <- length(output_info)
   out_names <- names(output_info)
 
-  if(is.null(veglib)) {
+  if (is.null(veglib)) {
     veglib <- veglib_LDAS
   } else {
-    if(is.data.frame(veglib)) {
-      veglib <- veglib[, !sapply(veglib, is.character)]
-      veglib <- as.matrix(veglib)
+    if (is.data.frame(veglib)) {
+      veglib <- veglib[, !sapply(veglib, is.character)] %>% as.matrix()
     }
   }
 
-  if(is.vector(soil))
+  if (is.vector(soil)) {
     soil <- t(soil)
+  }
   soil <- as.matrix(soil)
+
   ncell <- nrow(soil)
   cellid <- soil[, 2]
 
   # Check length of forcing data.
   minfl <- get_forclen()
-  for(ft in names(forcing)) {
+  for (ft in names(forcing)) {
     fl <- nrow(forcing[[ft]])
-    if(is.null(nrow(forcing[[ft]]))) fl <- length(forcing[[ft]])
-    if(fl < minfl) {
-      stop(sprintf('Length of forcing data "%s" (%d) is too short for model require (%d).',
-                   ft, fl, minfl))
+    if (is.null(nrow(forcing[[ft]]))) fl <- length(forcing[[ft]])
+    if (fl < minfl) {
+      stop(sprintf(
+        'Length of forcing data "%s" (%d) is too short for model require (%d).',
+        ft, fl, minfl
+      ))
     }
     fc <- ncol(forcing[[ft]])
-    if(is.null(fc) && ncell > 1 || fc < ncell)
+    if (is.null(fc) && ncell > 1 || fc < ncell) {
       stop(sprintf('Columns of forcing data "%s" (%d) are too few for model require. It should no less than number of gridcells (%d).', ft, fc, ncell))
+    }
   }
 
-  if(!is.list(veg)) veg <- list(veg)
+  if (!is.list(veg)) veg <- list(veg)
 
-  if(is.vector(lake)) {
+  if (is.vector(lake)) {
     lake <- t(lake)
-  } else if(!is.null(lake)) {
-    lake <- as.matrix(lake)
+  } else if (!is.null(lake)) {
+    lake %<>% as.matrix()
   }
 
-  if(is.vector(snowband)) {
-    snowband <- t(snowband)
-  } else if(!is.null(snowband)) {
-    snowband <- as.matrix(snowband)
+  if (is.vector(snowband)) {
+    snowband %<>% t()
+  } else if (!is.null(snowband)) {
+    snowband %<>% as.matrix()
   }
 
   forc_types <- c("PREC", "TEMP", "SW", "LW", "PRESS", "VP", "WIND")
   # The order of forc_types must not be change.
 
   forc_lack <- forc_types[!(forc_types %in% names(forcing))]
-  celev <- getOption('VIC_global_params')[['nlayers']]*4+10
-  if("LW" %in% forc_lack) J <- get_J()
+  celev <- getOption("VIC_global_params")[["nlayers"]] * 4 + 10
+  if ("LW" %in% forc_lack) J <- get_J()
 
-  if(ncell == 1) forcing <- data.frame(forcing[forc_types])
-  
-  globalopt <- getOption('VIC_global_params')
+  if (ncell == 1) forcing <- data.frame(forcing[forc_types])
+
+  globalopt <- getOption("VIC_global_params")
   `%dof%` <- ifelse(parall, foreach::`%dopar%`, foreach::`%do%`)
 
   tp2 <- proc.time()
 
-  tmp_out <- foreach::foreach(i=1:ncell) %dof% {
+  # 需要传递全局变量
+  tmp_out <- foreach::foreach(i = 1:ncell) %dof% {
+    # tmp_out <- list()
+    # for (i in 1:ncell) {
     # running(i, 100)
     forc_veg <- matrix() # default
-    if(!is.null(forcing_veg)) {
+    if (!is.null(forcing_veg)) {
       forc_veg <- forcing_veg[[i]]
       attr(forc_veg, "types") <- veg_force_types
     }
 
-    iband <- if (is.null(snowband)) -1 else snowband[i, ]
-    ilake <- if(is.null(lake)) -1 else lake[i, ]
+    iband <- snowband[i, ] %||% -1
+    ilake <- lake[i, ] %||% -1
 
-    if(is.null(nrow(veg[[i]])) || ncol(veg[[i]]) <= 1){
+    if (is.null(nrow(veg[[i]])) || ncol(veg[[i]]) <= 1) {
       iveg <- matrix(nrow = 0, ncol = 8)
     } else {
       iveg <- veg[[i]]
     }
 
-    forc <- sapply(forc_types, function(ft){
-      if(ft %in% forc_lack) {
+    forc <- sapply(forc_types, function(ft) {
+      if (ft %in% forc_lack) {
         rep(0, minfl)
       } else {
         forcing[[ft]][1:minfl, i]
       }
     })
 
     # Estimate forcing data that not supplied.
-    if("PRESS" %in% forc_lack) {
+    if ("PRESS" %in% forc_lack) {
       elev <- soil[i, celev]
-      forc[,5] <- 101.3*exp(-elev*9.81/(287*(273.15+forc[,2]-0.5*elev*0.0065)))
+      forc[, 5] <- cal_Pa(forc[, 2], elev)
     }
-    if("LW" %in% forc_lack) {
+    if ("LW" %in% forc_lack) {
       lat <- soil[i, 3]
-      forc[, 4] <- cal_lw(forc[,2], forc[,6], forc[,3], lat, J)
+      # cal_lw(temp, vp, rdsn. lat, J)
+      forc[, 4] <- cal_lw(forc[, 2], forc[, 6], forc[, 3], lat, J)
       # temp, vp, rsds, lat, J
     }
-
-    p <- listk(globalopt,
-             forc, soil[i, ], iband, iveg,
-             ilake, forc_veg, veglib, output_info)
+
+    # p <- listk(
+    #   globalopt,
+    #   forc, soil[i, ], iband, iveg,
+    #   ilake, forc_veg, veglib, output_info
+    # )
     # print(str(p))
-    iout <- vic_run_cell(globalopt,
-                         forc, soil[i, ], iband, iveg,
-                         ilake, forc_veg, veglib, output_info)
+
+    iout <- vic_run_cell(
+      globalopt,
+      forc, soil[i, ], iband, iveg,
+      ilake, forc_veg, veglib, output_info
+    )
+    # tmp_out[[i]] = iout
     iout
   }
 
   tp3 <- proc.time()
-
   # Run time in C
-  #cpp_time <- rbind(
+  # cpp_time <- rbind(
   #  c(rowSums(sapply(tmp_out, function(x)x[[n_outputs + 1]])), 0),
   #  c(rowSums(sapply(tmp_out, function(x)x[[n_outputs + 2]])), 0),
   #  c(rowSums(sapply(tmp_out, function(x)x[[n_outputs + 3]])), 0))
   out <- list()
-  for(i in 1:n_outputs) {
-    out_header <- attr(tmp_out[[1]][[i]], 'header')
-    out_time <- attr(tmp_out[[1]][[i]], 'time')
+
+  for (i in 1:n_outputs) {
+    out_header <- attr(tmp_out[[1]][[i]], "header")
+    out_time <- attr(tmp_out[[1]][[i]], "time")
+
     iout <- lapply(1:length(out_header), function(h) {
-      itype <- sapply(1:ncell, function(g) tmp_out[[g]][[i]][, h])
-      rownames(itype) <- out_time
-      colnames(itype) <- cellid
+      itype <- sapply(1:ncell, function(g) tmp_out[[g]][[i]][, h]) %>% 
+        set_dimnames(list(out_time, cellid))
       itype
     })
     names(iout) <- out_header
@@ -399,14 +402,35 @@ vic <- function(forcing, soil, veg,
 
   tp4 <- proc.time()
 
-  time_table <- rbind(tp2-tp1, tp3-tp2, tp4-tp3, tp4-tp1)[,1:3]
-  #time_table <- rbind(time_table, cpp_time)
-  #rownames(time_table) <- c('init_time', 'run_time', 'final_time', 'total_time',
+  time_table <- rbind(tp2 - tp1, tp3 - tp2, tp4 - tp3, tp4 - tp1)[, 1:3] %>% 
+    set_rownames(c("init_time", "run_time", "final_time", "total_time")) %>% 
+    round(4)
+
+  # time_table <- rbind(time_table, cpp_time)
+  # rownames(time_table) <- c('init_time', 'run_time', 'final_time', 'total_time',
   #                       'init_time(cpp)', 'run_time(cpp)', 'final_time(cpp)')
-  rownames(time_table) <- c('init_time', 'run_time', 'final_time', 'total_time')
-  out[['timetable']] <- round(time_table, 4)
-  out[['global_options']] <- globalopt
-  out[['output_info']] <- output_info
-  class(out) <- 'vic_output'
+  out[["timetable"]] <- time_table
+  out[["global_options"]] <- globalopt
+  out[["output_info"]] <- output_info
+  class(out) <- "vic_output"
   out
 }
+
+#' @import magrittr
+default_output_info <- list(
+  fluxes = list(
+    timescale = "day", aggpar = 1,
+    outvars = c(
+      "OUT_PREC", "OUT_EVAP", "OUT_RUNOFF",
+      "OUT_BASEFLOW", "OUT_WDEW", "OUT_SOIL_MOIST"
+    )
+  ),
+  snow = list(
+    timescale = "day", aggpar = 1,
+    outvars = c(
+      "OUT_SWE", "OUT_SNOW_DEPTH", "OUT_SNOWF",
+      "OUT_SNOW_MELT", "OUT_SNOW_SURF_TEMP",
+      "OUT_SNOW_PACK_TEMP"
+    )
+  )
+)

R/vic_forcing.R

@@ -20,6 +20,11 @@ get_J <- function() {
   J
 }
 
+# air pressure in kPa
+cal_Pa <- function(Tair, elev = 0) {
+  101.3 * exp(-elev * 9.81 / (287 * (273.15 + Tair - 0.5 * elev * 0.0065)))
+}
+
 cal_lw <- function(temp, vp, rsds, lat, J) {
   # Calc s
   lat <- lat * pi/360