Merge branch 'feature/install-execute' into feature/install-and-execute

rogiersbart · Jan 30, 2020 · 79bec06 · 79bec06
2 parents fde12ae + d90230d
commit 79bec06
Show file tree

Hide file tree

Showing 4 changed files with 246 additions and 177 deletions.
diff --git a/.gitignore b/.gitignore
@@ -2,4 +2,5 @@
 .Rhistory
 .RData
 inst/doc
-code
+code
+inst/code
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -19,8 +19,6 @@ Imports:
     animation,
     tools,
     lubridate,
-    xml2,
-    rvest,
     sf,
     tibble,
     dplyr,

diff --git a/R/execute.R b/R/execute.R
@@ -1,6 +1,6 @@
 #' Run a MODFLOW model
 #' 
-#' \code{execute} runs a MODFLOW model.
+#' \code{rmf_execute} runs a MODFLOW model.
 #' 
 #' @param path path to name file; typically '*.nam'
 #' @param code name of the MODFLOW variant to use
@@ -12,6 +12,7 @@
 #' @rdname rmf_execute
 #' @method rmf_execute character
 #' @export
+#' @seealso \code{\link{rmf_install}}
 rmf_execute.character <- function(
   path,
   code = 2005,
@@ -27,12 +28,21 @@ rmf_execute.character <- function(
     dir <- dirname(path)
     file <- basename(path)
 
-  # set initial parameters if provided  
+    # set initial parameters if provided  
     if(!is.null(par)) {
-      nam <- rmf_read_nam(paste0(dir,'/',file))
-      pvl <- rmf_read_pvl(paste0(dir,'/',nam$fname[which(nam$ftype=='PVAL')]))
-      pvl$parval <- par
-      rmf_write_pvl(pvl, file=paste0(dir,'/',nam$fname[which(nam$ftype=='PVAL')]))
+      nam <- rmf_read_nam(path)
+      if('PVAL' %in% nam$ftype) {
+        pvl <- rmf_read_pvl(file.path(dir, nam$fname[which(nam$ftype == 'PVAL')]))
+
+        # copy and write original pvl
+        rmf_write_pvl(pvl, file = file.path(dir, paste('original', nam$fname[which(nam$ftype=='PVAL')], sep = '_')))
+
+        pvl$parval <- par
+        rmf_write_pvl(pvl, file = file.path(dir, nam$fname[which(nam$ftype == 'PVAL')]))
+
+      } else {
+        warning('Model does not have PVL input file. Ignoring par argument', call. = FALSE)
+      }
     }
 
   # run modflow
@@ -95,60 +105,71 @@ rmf_execute <- function(...) {
   UseMethod('rmf_execute')
 }
 
-#' Run a MODFLOW model optimization, based on the parameter value file
+#' Run a MODFLOW model optimization, based on the parameter value file and the head predictions output file
 #' 
 #' \code{rmf_optimize} runs a MODFLOW optimization.
 #' 
 #' @param path path to name file; typically '*.nam'
 #' @param code name of the MODFLOW variant to use
-#' @param par initial parameter values (for all or only included parameters); current parameter value file values are used if par is not provided
+#' @param par initial parameter values (for all or only included parameters); current parameter value file values are used if par is not provided (default)
 #' @param include logical vector indicating which parameters in the parameter value file to include in the optimization
-#' @param trans vector of transformations; currently only 'log' is supported
-#' @param method optimization method: 'Nelder-Mead','BFGS','CG','SANN','Brent','L-BGFS-B'
 #' @param lower lower parameter bounds
 #' @param upper upper parameter bounds
 #' @param control list of control arguments
+#' @param logtrans logical vector indicating which parameters should be logtransformed
+#' @param verbose logical; should the terminal output be printed to the R console? Defaults to TRUE
+#' @param cost character denoting which statistic should be used to compute the cost. Possible values are those returned by \code{\link{rmf_performance}}. Defaults to 'ssq'
+#' @param out_file optional file path to write cost value, convergence and parameter values to for each iteration
+#' @param cvg_message character denoting the message in the terminal output used to check for convergence
 #' @param ... further arguments provided to \code{optim}
+#'
+#' @details Only works with models using MODFLOW parameters and having a head predictions output file as defined in the HOB object
+#'          The only method currently available is "L-BFGS-B". See \code{\link{optim}}
 #' @return \code{optim} results with the full list of parameters
 #' @export
+#' @seealso \code{\link{optim}}, \code{\link{rmf_execute}}, \code{\link{rmf_create_pvl}} & \code{\link{rmf_create_hob}} 
 rmf_optimize <- function(
   path,
   code = 2005,
   par = NULL,
   include = NULL,
-  trans = NULL,
-  method = 'Nelder-Mead',
-  lower = -Inf,
+  logtrans = NULL,
+  lower = 0, # TODO think about this default value ...
   upper = Inf,
+  verbose = TRUE,
+  cost = 'ssq',
+  out_file = NULL,
   control = list(),
+  cvg_message = 'Normal termination',
   ...
 ) {
   code <- rmf_find(code)
   dir <- dirname(path)
-  file <- basename(path)
-  nam <- rmf_read_nam(paste0(dir,'/',file))
-  pvl <- rmf_read_pvl(paste0(dir,'/',nam$fname[which(nam$ftype=='PVAL')]))
-  hob <- rmf_read_hob(paste0(dir,'/',nam$fname[which(nam$ftype=='HOB')]))
+  nam <- rmf_read_nam(path)
+  if(!('PVAL' %in% nam$ftype) || !('HOB' %in% nam$ftype)) stop('rmf_optimize only works with models having PVL and HOB input', call. = FALSE)
+  pvl <- rmf_read_pvl(file.path(dir, nam$fname[which(nam$ftype=='PVAL')]))
+  hob <- rmf_read_hob(file.path(dir, nam$fname[which(nam$ftype=='HOB')]))
+  if(hob$iuhobsv == 0 || toupper(nam$ftype[which(nam$nunit == hob$iuhobsv)]) != 'DATA') {
+    stop('rmf_optimize only works with models having a head predictions output file. This can be created by setting the iuhobsv argument of the HOB file to a non-zero value', call. = FALSE)
+  }
+
+  # copy and write original pvl
+  rmf_write_pvl(pvl, file = file.path(dir, paste('original', nam$fname[which(nam$ftype=='PVAL')], sep = '_')))
+
   if(is.null(par)) par <- pvl$parval
   if(is.null(include)) include <- rep(TRUE,length(par))
-  if(length(par)!=length(pvl$parval)) 
-  {
-    par2 <- par
-    par <- pvl$parval
-    par[which(include)] <- par2
-  }
+
   if(is.infinite(lower)) lower <- rep(lower,pvl$np)
   if(is.infinite(upper)) upper <- rep(upper,pvl$np)
-  if(!is.null(trans))
-  {
-    par[which(trans=='log')] <- log(par[which(trans=='log')])
-    lower[which(trans=='log')] <- log(lower[which(trans=='log')])
-    upper[which(trans=='log')] <- log(upper[which(trans=='log')])
-    if('parscale' %in% names(control))
-    {
-      control$parscale[which(trans[include]=='log')] <- log(control$parscale[which(trans[include]=='log')])
+  if(!is.null(logtrans)) {
+    par[which(logtrans)] <- log(par[which(logtrans)])
+    lower[which(logtrans)] <- log(lower[which(logtrans)])
+    upper[which(logtrans)] <- log(upper[which(logtrans)])
+    if('parscale' %in% names(control)) {
+      control$parscale[which(trans[include])] <- log(control$parscale[which(trans[include])])
     }
   }
+<<<<<<< HEAD
   optim_modflow <- function(par_include)
   {
     pvl$parval <- par
@@ -175,55 +196,113 @@ rmf_optimize <- function(
   opt$par <- par
   opt$par[which(include)] <- par2
   if(!is.null(trans)) opt$par[which(trans=='log')] <- exp(opt$par[which(trans=='log')])
+=======
+
+  # if par only has values for include = TRUE
+  if(length(par) != length(pvl$parval)) {
+    par2 <- par
+    par <- pvl$parval
+    par[which(include)] <- par2
+  }
+
+  optim_modflow <- function(par_include) {
+    # adjust values
+    pvl$parval <- par
+    pvl$parval[which(include)] <- par_include
+    if(!is.null(logtrans)) pvl$parval[which(logtrans)] <- exp(pvl$parval[which(logtrans)])
+
+    # write values
+    rmf_write_pvl(pvl, file = file.path(dir, nam$fname[which(nam$ftype=='PVAL')]))
+
+    # run modflow
+    cvg <- rmf_execute(file = file, executable = executable, version = version, verbose = verbose, cvg_message = cvg_message, convergence = TRUE)
+
+    # get cost
+    if(cvg) {
+      hpr <- rmf_read_hpr(file.path(dir, nam$fname[which(nam$nunit==hob$iuhobsv)]))
+      cost_value <- rmf_performance(hpr)[cost]
+    } else { # return large cost when not converging
+      cost_value <- ifelse(!is.null(control$fnscale) && control$fnscale < 0, -1e100, 1e100) # maximization/minimization
+    }
+
+    # print 
+    if(verbose) cat(paste(cost, '=', format(cost_value,scientific=TRUE,digits=4), 'converged =', as.character(cvg), 'parval =', paste(format(pvl$parval[include],scientific=TRUE,digits=4), collapse=' '),'\n'))
+    if(!is.null(out_file)) cat(paste(cost, '=', format(cost_value,scientific=TRUE,digits=4), 'converged =', as.character(cvg), 'parval =', paste(format(pvl$parval[include],scientific=TRUE,digits=4), collapse=' '),'\n'), file = out_file, append = TRUE)
+
+    return(cost_value)
+  }
+
+  # optimize; L-BGFS-B only
+  opt <- optim(par[which(include)], optim_modflow, method = 'L-BFGS-B', lower = lower[which(include)], upper = upper[which(include)], control = control, ...)
+>>>>>>> install-execute
   opt$included <- include
   return(opt)
 }
 
-#' Run a MODFLOW model sensitivity analysis, based on the parameter value file
+#' Run a MODFLOW model sensitivity analysis, based on the parameter value file and the head predictions output file
 #' 
 #' \code{rmf_analyze} performs a MODFLOW model sensitivity analysis.
 #' 
 #' @param path path to name file; typically '*.nam'
 #' @param code name of the MODFLOW variant to use
-#' @param par central parameter values (for all or only included parameters); parameter value file values are used if par is not provided
+#' @param par central parameter values (for all or only included parameters); parameter value file values are used if par is not provided (default)
 #' @param include logical vector indicating which parameters in the parameter value file to include in the sensitivity analysis
-#' @return sensitivity analysis results
+#' @param ... optional arguments passed to \code{rmf_execute}
+#' @details Only works with models using MODFLOW parameters and having a head predictions output file as defined in the HOB object
+#' @return object of class \code{sen} which is a list with dimensionless scaled sensitivities (dss) and composite scaled sensitivies (css)
 #' @export
-rmf_analyze <- function(
-  path,
-  code = 'mf2005',
-  par = NULL,
-  include = NULL
-) {
+#' @seealso \code{\link{rmf_execute}}, \code{\link{rmf_create_pvl}} & \code{\link{rmf_create_hob}} 
+rmf_analyze <- function(path,
+                        code = "2005",
+                        par = NULL,
+                        include = NULL,
+                        ...) {
   code <- rmf_find(code)
   dir <- dirname(path)
-  file <- basename(path)
-  nam <- rmf_read_nam(paste0(dir,'/',file))
-  pvl <- rmf_read_pvl(paste0(dir,'/',nam$fname[which(nam$ftype=='PVAL')]))
-  hob <- rmf_read_hob(paste0(dir,'/',nam$fname[which(nam$ftype=='HOB')]))
+  nam <- rmf_read_nam(path)
+  if(!('PVAL' %in% nam$ftype) || !('HOB' %in% nam$ftype)) stop('rmf_analyze only works with models having PVL and HOB input', call. = FALSE)
+  pvl <- rmf_read_pvl(file.path(dir, nam$fname[which(nam$ftype=='PVAL')]))
+  hob <- rmf_read_hob(file.path(dir, nam$fname[which(nam$ftype=='HOB')]))
+  if(hob$iuhobsv == 0 || toupper(nam$ftype[which(nam$nunit == hob$iuhobsv)]) != 'DATA') {
+    stop('rmf_analyze only works with models having a head predictions output file. This can be created by setting the iuhobsv argument of the HOB file to a non-zero value', call. = FALSE)
+  }
+  pvl_org <- pvl
   if(is.null(par)) par <- pvl$parval
   if(is.null(include)) include <- rep(TRUE,length(par))
-  if(length(par)!=length(pvl$parval)) 
-  {
+
+  # if par only has values for include = TRUE
+  if(length(par) != length(pvl$parval)) {
     par2 <- par
     par <- pvl$parval
     par[which(include)] <- par2
   } 
-  rmf_execute(file = paste0(dir,'/',file), code = code,par)
-  hpr_orig <- rmf_read_hpr(paste0(dir,'/',nam$fname[which(nam$nunit==hob$iuhobsv)]))
+
+  # copy and write original pvl
+  rmf_write_pvl(pvl_org, file = file.path(dir, paste('original', nam$fname[which(nam$ftype=='PVAL')], sep = '_')))
+
+  # initial run
+  rmf_execute(path = path, code = code, par = par, version = version, convergence = FALSE, ...)
+  hpr_orig <- rmf_read_hpr(file.path(dir, nam$fname[which(nam$nunit == hob$iuhobsv)]))
+
   sens <- list()
-  sens$dss <- matrix(NA,nrow=length(hob$obsnam),ncol=length(pvl$parval))
+  sens$dss <- matrix(NA, nrow = length(hob$obsnam), ncol = length(pvl$parval))
   sens$css <- pvl$parval*NA
-  for(i in which(include))
-  {
+
+  # dss & css
+  for(i in which(include)) {
     pvl$parval <- par
     pvl$parval[i] <- pvl$parval[i]*1.01
-    rmf_write_pvl(pvl, file=paste0(dir,'/',nam$fname[which(nam$ftype=='PVAL')]))
-    rmf_execute(file = paste0(dir,'/',file), code = code)
-    hpr <- rmf_read_hpr(paste0(dir,'/',nam$fname[which(nam$nunit==hob$iuhobsv)]))
-    sens$dss[,i] <- (hpr$simulated_equivalent-hpr_orig$simulated_equivalent)/(0.01)
+    rmf_write_pvl(pvl, file = file.path(dir, nam$fname[which(nam$ftype == 'PVAL')]))
+    rmf_execute(file = file, executable = executable, par = NULL, version = version, convergence = FALSE, ...)
+    hpr <- rmf_read_hpr(file.path(dir, nam$fname[which(nam$nunit == hob$iuhobsv)]))
+    sens$dss[,i] <- (hpr$simulated - hpr_orig$simulated)/(0.01)
     sens$css[i] <- sqrt(sum(sens$dss[,i]^2)/hob$nh)
   }
+
+  #  write original pvl and run
+  rmf_write_pvl(pvl_org, file = file.path(dir, nam$fname[which(nam$ftype=='PVAL')]))
+  rmf_execute(file = file, executable = executable, par = NULL, version = version, convergence = FALSE, ...)
+
   sens$parnam <- pvl$parnam
   class(sens) <- 'sen'
   return(sens)
@@ -253,7 +332,9 @@ rmf_find <- function(
   precision = "single"
 ) {
   # TODO automatically select 32 bit executables if available, otherwise throw
-  # error on 32 bit systems
+  #      error on 32 bit systems
+  # TODO check if executable provided in name exists in current wd, not name
+  #      of executable that we are using in RMODFLOW!
   if (grepl("2005", name)) {
     if (grepl("dbl", name)) precision <- "double"
     name <- "MODFLOW-2005"
-Original file line number
+Diff line change
@@ Expand Up / @@ -19,8 +19,6 @@ Imports: @@
         animation,
         tools,
         lubridate,
-        xml2,
-        rvest,
         sf,
         tibble,
         dplyr,
@@ Expand Down @@