# calibrar package: Automated Calibration for Complex Models --------
#' calibrar package: Automated Calibration for Complex Models
#'
#'\if{html}{\figure{logo_small.png}{options: style='float: right' alt='logo' width='150'}}
#'
#' This package allows the parameter estimation (i.e. calibration) of complex models, including stochastic ones. It implements generic functions that can be used for fitting any type of models, especially those with non-differentiable objective functions, with the same syntax as base::optim.
#' It supports multiple phases estimation (sequential parameter masking), constrained optimization (bounding box restrictions) and automatic parallel computation of numerical gradients.
#' Some common maximum likelihood estimation methods and automated construction of the objective function from simulated model outputs is provided.
#' See <https://roliveros-ramos.github.io/calibrar/> for more details.
#'
#' @title Automated Calibration for Complex Models
#' @description Automated Calibration for Complex Models
#' @name calibrar-package
#' @aliases calibrar-package calibrar
#' @docType package
#' @author Ricardo Oliveros-Ramos Maintainer: Ricardo Oliveros-Ramos
#' <ricardo.oliveros@@gmail.com>
#' @references calibrar: an R package for the calibration of ecological models (Oliveros-Ramos and Shin 2014)
#' @keywords calibration
#' @examples
#' \dontrun{
#' require(calibrar)
#' set.seed(880820)
#' path = NULL # NULL to use the current directory
#' # create the demonstration files
#' demo = calibrar_demo(model="PoissonMixedModel", L=5, T=100)
#' # get calibration information
#' calibrationInfo = calibration_setup(file=demo$path)
#' # get observed data
#' observed = calibration_data(setup=calibrationInfo, path=demo$path)
#' # read forcings for the model
#' forcing = read.csv(file.path(demo$path, "master", "environment.csv"), row.names=1)
#' # Defining 'runModel' function
#' runModel = function(par, forcing) {
#' output = calibrar:::.PoissonMixedModel(par=par, forcing=forcing)
#' # adding gamma parameters for penalties
#' output = c(output, list(gammas=par$gamma))
#' return(output)
#' }
#' # real parameters
#' cat("Real parameters used to simulate data\n")
#' print(demo$par)
#' # objective functions
#' obj = calibration_objFn(model=runModel, setup=calibrationInfo,
#' observed=observed, forcing=forcing)
#' cat("Starting calibration...\n")
#' control = list(weights=calibrationInfo$weights, maxit=3.6e5) # control parameters
#' cat("Running optimization algorithms\n", "\t", date(), "\n")
#' cat("Running optim AHR-ES\n")
#' ahr = calibrate(par=demo$guess, fn=obj, lower=demo$lower, upper=demo$upper, control=control)
#' summary(ahr)
#' }
NULL
# calibrate ---------------------------------------------------------------
#' @title Sequential parameter estimation for the calibration of complex models
#' @description This function performs the optimization of a function, possibly
#' in sequential phases of increasing complexity, and it is designed for the
#' calibration of a model, by minimizing the error function \code{fn} associated to it.
#' @param par A numeric vector or list. The length of the par argument defines the
#' number of parameters to be estimated (i.e. the dimension of the problem).
#' @param fn The function to be minimized.
#' @param gr A function computing the gradient of \code{fn}. If NULL, a numerical approximation
#' of the gradient is used. It can be also a character specifying the method for the computation
#' of the numerical gradient: 'central', 'forward' (the default), 'backward' or 'richardson'.
#' @param \dots Additional parameters to be passed to \code{fn}.
#' @param lower Lower threshold value(s) for parameters. One value or a vector
#' of the same length as par. If one value is provided, it is used for all
#' parameters. \code{NA} means \code{-Inf}. By default \code{-Inf} is used (unconstrained).
#' @param upper Upper threshold value(s) for parameters. One value or a vector
#' of the same length as par. If one value is provided, it is used for all
#' parameters. \code{NA} means \code{Inf}. By default \code{Inf} is used (unconstrained).
#' @param phases An optional vector of the same length as \code{par},
#' indicating the phase at which each parameter becomes active. If omitted,
#' default value is 1 for all parameters, performing a single optimization.
#' @param method The optimization method to be used. The default method
#' is the AHR-ES (Adaptative Hierarchical Recombination Evolutionary Strategy,
#' Oliveros-Ramos & Shin, 2016). See details for the methods available.
#' @param control Parameter for the control of the algorithm itself, see details.
#' @param hessian Logical. Should a numerically differentiated Hessian matrix be returned?
#' Currently not implemented.
#' @param replicates The number of replicates for the evaluation of \code{fn}.
#' The default value is 1. A value greater than 1 is only useful for stochastic
#' functions.
#' @param parallel Logical. Use parallel computation numerical of gradient?
#' @details In the control list, \code{aggFn} is a function to aggregate \code{fn} to
#' a scalar value if the returned value is a vector. Some optimization algorithm can
#' exploite the additional information provided by a vectorial output from \code{fn}.
#' @author Ricardo Oliveros-Ramos
#' @examples
#' calibrate(par=rep(NA, 5), fn=sphereN)
#' \dontrun{
#' calibrate(par=rep(NA, 5), fn=sphereN, replicates=3)
#' calibrate(par=rep(0.5, 5), fn=sphereN, replicates=3, lower=-5, upper=5)
#' calibrate(par=rep(0.5, 5), fn=sphereN, replicates=3, lower=-5, upper=5, phases=c(1,1,1,2,3))
#' calibrate(par=rep(0.5, 5), fn=sphereN, replicates=c(1,1,4), lower=-5, upper=5, phases=c(1,1,1,2,3))
#' }
#' @family optimisers
#' @export
calibrate = function(par, fn, gr, ..., method, lower, upper, phases, control,
hessian, replicates, parallel) {
UseMethod("calibrate", fn)
}
#' @export
#' @rdname calibrate
calibrate.default = function(par, fn, gr = NULL, ..., method = NULL,
lower = NULL, upper = NULL, phases = NULL,
control = list(),
hessian = FALSE, replicates=1, parallel=FALSE) {
if(is.null(method)) {
method = if(all(replicates==1)) "Rvmmin" else "AHR-ES"
}
methods = c("L-BFGS-B", "nlminb", "Rcgmin", "Rvmmin",
"hjn", "spg", "LBFGSB3", "CMA-ES", "genSA", "DE", "soma", "genoud", "PSO",
"hybridPSO", "mads", "hjkb", "nmkb", "bobyqa", "AHR-ES", "Nelder-Mead",
"CG", "BFGS", "SANN")
method = match.arg(method, choices=methods, several.ok = TRUE)
# calibrate should dispatch on 'fn'
if(inherits(fn, "objFn")) {
control$nvar = attr(fn, "nvar")
control$weights = attr(fn, "weights")
names(control$weights) = attr(fn, "variables")
}
# check function and method
if(!all(method %in% multiMethods)) {
if(inherits(fn, "objFn")) {
agg = attr(fn, "aggregate")
if(is.null(agg)) warning("Update your objective function to the last version of the 'calibrar' package.")
if(!isTRUE(agg))
stop(sprintf("Method '%s' does not support multi-objective optimization, use aggregate=TRUE in 'calibration_objFn'.",
method))
}
if(!is.null(control$nvar))
if(control$nvar != 1) stop(sprintf("Method '%s' does not support multi-objective optimization, 'fn' must be scalar.",
method))
control$nvar = 1
}
if(inherits(fn, "objFn")) control$aggFn = attr(fn, "aggFn")
fnx = attr(fn, "fn")
if(!is.null(fnx)) fnx = match.fun(fnx)
# if 'gr' is a character, assume it specify the type of numerical approximation
if(!is.null(gr) & is.character(gr)) {
if(!is.null(control$gr.method))
warning(sprintf("Ignoring provided gr.method='%s', using gr='%s'.",
control$gr.method, gr))
control$gr.method = gr
gr = NULL
}
control = .checkControl_calibrate(control=control, method=method, par=par, fn=fn, ...)
# here, we calculate 'batchsize' according to the method, copyMaster is called just once.
# check for a partial results file to restart from a completed phase
restart = .restartCalibration(control, type="partial")
skeleton = as.relistable(par)
par = unlist(par)
lower = unlist(lower)
upper = unlist(upper)
phases = unlist(phases)
npar = length(par)
# checking conformance of all arguments
if(is.null(names(par))) names(par) = .printSeq(npar, preffix="par")
fn = match.fun(fn)
if(!is.null(gr)) gr = match.fun(gr)
phases = .checkPhases(phases=phases, npar=npar)
bounds = .checkBounds(lower=lower, upper=upper, npar=npar)
guess = .checkOpt(par=par, lower=bounds$lower, upper=bounds$upper)
par = guess
lower = bounds$lower
upper = bounds$upper
nphases = max(phases, na.rm=TRUE)
replicates = .checkReplicates(replicates, nphases)
control$repmax = max(replicates, na.rm=TRUE)
output = if(isTRUE(restart)) .getResults(control=control, type="partial") else list(phase=0)
output$phase = output$phase + 1
msg0 = sprintf("\nStarting calibration from phase %d.\n", output$phase)
msg1 = sprintf("\nRestarting calibration from phase %d.\n", output$phase)
if(isTRUE(restart)) message(msg1)
if(!isTRUE(restart) & isTRUE(control$verbose)) message(msg0)
conv = .checkConvergence(control, nphases, method)
control$parallel = parallel
if(!is.null(control$master)) message(sprintf("Using 'master' directory: %s", control$master))
if(!is.null(control$run)) message(sprintf("Using 'run' directory: %s", control$run))
xtm1 = Sys.time() # timer for all phases
# -------------------------------------------
# start the sequential parameter estimation
# -------------------------------------------
if(nphases > 1)
# message(sprintf("PARAMETER ESTIMATION IN %s PHASES.", toupper(make_number(nphases))))
message(sprintf("Parameter estimation in %s phases.", make_number(nphases)))
for(phase in seq(from=output$phase, to=nphases)) {
if(output$phase > nphases) break
control$maxgen = conv$maxgen[phase]
control$maxit = conv$maxit[phase]
control$convergence = conv$convergence[phase]
control$abstol = conv$abstol[phase]
control$reltol = conv$reltol[phase]
method = conv$method[phase]
active = (phases <= phase) # NAs are corrected in .calibrar
msg2 = sprintf("\n- Phase %d: %d of %d parameters active.",
phase, sum(active, na.rm=TRUE), npar)
msg3 = sprintf("\tStarted at %s.\n", date())
if(nphases > 1) message(msg2)
if(isTRUE(control$verbose)) message(msg3)
if(nphases==1) message(sprintf("Using optimization method '%s'.", method))
if(nphases!=1) message(sprintf("\tUsing optimization method '%s'.", method))
# call optimisers using .calibrar handler
tm1 = Sys.time()
temp = .calibrar(par=par, fn=fn, gr = gr, ..., method = method,
lower = lower, upper = upper, control = control,
hessian = hessian, active=active, skeleton=skeleton,
replicates=replicates[phase])
tm2 = Sys.time()
output$phases[[phase]] = temp # trim?
output$phase = phase
.createOutputFile(output, control, type="partial", phase=phase)
par = temp$par #update parameter guess
control = .updateControl(control=control, opt=temp, method=method) # update CVs?
if(!is.null(control$restart) & !is.null(fnx)) {
wd = getwd()
.setWorkDir(control$run, i=0)
xpar = relist(flesh = par, skeleton = skeleton)
ifile = sprintf(".%s-phase%d.par", control$restart, phase)
saveRDS(xpar, file=ifile)
partial = try(fnx(xpar))
setwd(wd)
if(!inherits(partial, "try-error")) {
ifile = sprintf("%s-phase%d.simulated", control$restart, phase)
saveRDS(partial, file=ifile)
}
}
xmsg = temp$message
if(is.null(xmsg)) xmsg = NA
if(is.na(xmsg)) xmsg = "-"
msg = c(sprintf("Phase %d finished (%s)", phase, format_difftime(tm1, tm2)),
sprintf("Elapsed time: %s", format_difftime(tm1, tm2)),
sprintf("Function value: %g", temp$value),
paste(c("Parameter values:",sprintf("%0.3g", par[which(active)])), collapse=" "))
msg = if(nphases == 1) msg[-1] else paste("\t", msg[-2], sep="")
msg = if(phase == nphases) c(msg, sprintf("\nStatus: %s", xmsg)) else msg
msg = paste(msg, collapse="\n")
message(msg)
} # end of phases
xtm2 = Sys.time() # timer for all phases
isActive = (phases>0) & !is.na(phases)
paropt = output$phases[[nphases]]$par # save parameters of last phase
paropt = relist(paropt, skeleton)
class(paropt) = setdiff(class(paropt), "relistable")
trace = output$phases[[nphases]]$trace # must be a list
if(is.null(trace)) trace = list()
if(!is.null(output$phases[[nphases]]$generations))
trace$generations = output$phases[[nphases]]$generations
trace = c(trace, list(partial=output$phases[[nphases]]$partial), output)
final = list(par = paropt,
value = output$phases[[nphases]]$value,
counts = output$phases[[nphases]]$counts,
convergence = output$phases[[nphases]]$convergence,
message = output$phases[[nphases]]$message,
method = method,
fn = fn,
active = isActive,
elapsed = format_difftime(xtm1, xtm2, value = TRUE),
trace = trace)
# output = c(final, output)
class(final) = c("calibrar.results")
.createOutputFile(final, control, type="results")
return(final)
}
#' @export
calibrate.TMB = function(par, fn, gr = NULL, ..., method = NULL,
lower = NULL, upper = NULL, phases = NULL,
control = list(),
hessian = NULL, replicates=1, parallel=FALSE) {
methods = c("nlminb", "BFGS", "L-BFGS-B", "Rcgmin", "Rvmmin", "spg",
"LBFGSB3", "AHR-ES")
if(!is.null(method)) fn$method = method
method = match.arg(fn$method, choices=methods)
if(!is.null(hessian)) fn$hessian = hessian
out = calibrate(par=par, fn=fn$fn, gr=fn$gr, lower=lower, upper=upper,
phases=phases, method=method, control=control,
hessian=fn$hessian, replicates=1, parallel=parallel)
return(out)
}
# optim2 ------------------------------------------------------------------
#' General-purpose optimization with parallel numerical gradient computation
#'
#' @param active Boolean vector of the same length as par, indicating if the
#' parameter is used in the optimization (TRUE) or hold at a fixed value (FALSE).
#'
#' @return
#' A list with components:
#' \describe{
#' \item{par}{The best set of parameters found.}
#' \item{value}{The value of fn corresponding to par.}
#' \item{counts}{A two-element integer vector giving the number of calls to fn and gr respectively. This excludes those calls needed to compute the Hessian, if requested, and any calls to fn to compute a finite-difference approximation to the gradient.}
#' \item{convergence}{An integer code. 0 indicates successful completion. }
#' \item{message}{A character string giving any additional information returned by the optimizer, or NULL.}
#' \item{hessian}{Only if argument hessian is true. A symmetric matrix giving an estimate of the Hessian at the solution found. Note that this is the Hessian of the unconstrained problem even if the box constraints are active.}
#' }
#' @examples
#' optim2(par=rep(NA, 5), fn=sphereN)
#' @inheritParams calibrate
#' @author Ricardo Oliveros-Ramos
#' @family optimisers
#' @export
optim2 = function(par, fn, gr = NULL, ...,
method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN",
"Brent", "nlm", "nlminb", "Rcgmin", "Rvmmin", "hjn",
"spg", "LBFGSB3", "AHR-ES"),
lower = -Inf, upper = +Inf, active = NULL,
control = list(), hessian = FALSE, parallel=FALSE) {
method = match.arg(method)
out = .optim2(par=par, fn=fn, gr=gr, ..., lower=lower, upper=upper, active=active,
method=method, control=control, hessian=hessian, parallel=parallel)
return(out)
}
# optimh ------------------------------------------------------------------
#' General-purpose optimization using heuristic algorithms
#'
#' @param gr Function to compute the gradient of \code{fn}. Ignored by most methods,
#' added for consistency with other optimization functions.
#' @examples
#' optim2(par=rep(NA, 5), fn=sphereN)
#' @inheritParams optim2
#' @inherit optim2 return author
#' @family optimisers
#' @export
optimh = function(par, fn, gr = NULL, ...,
method = c("AHR-ES", "Nelder-Mead", "SANN", "hjn", "bobyqa",
"CMA-ES", "genSA", "DE", "soma", "genoud", "PSO",
"hybridPSO", "mads", "hjk", "hjkb", "nmk", "nmkb"),
lower = -Inf, upper = +Inf, active = NULL,
control = list(), hessian = FALSE, parallel=FALSE) {
method = match.arg(method)
out = .optim2(par=par, fn=fn, gr=gr, ..., lower=lower, upper=upper, active=active,
method=method, control=control, hessian=hessian, parallel=parallel)
return(out)
}
# AHR-ES ------------------------------------------------------------------
#' @title Adaptative Hierarchical Recombination Evolutionary Strategy (AHR-ES) for
#' derivative-free and black-box optimization
#' @description This function performs the optimization of a function using
#' the Adaptative Hierarchical Recombination Evolutionary Strategy (AHR-ES, Oliveros & Shin, 2015).
#'
#' @author Ricardo Oliveros-Ramos
#' @examples
#' \dontrun{ahres(par=rep(1, 5), fn=sphereN)}
#' @inheritParams calibrate
#' @inheritParams optim2
#' @inherit optim2 return
#' @author Ricardo Oliveros-Ramos
#' @family optimisers
#' @export
ahres = function(par, fn, gr = NULL, ..., lower = -Inf, upper = +Inf, active = NULL,
control = list(), hessian = FALSE, parallel=FALSE) {
# par can be a list
skeleton = as.relistable(par)
par = unlist(par)
lower = unlist(lower)
upper = unlist(upper)
active = unlist(active)
npar = length(par)
# par can be re-scaled
parscale = .checkParscale(control=control, npar=npar)
control$parscale = NULL # reset parscale
par = par/parscale
lower = lower/parscale
upper = upper/parscale
# function can be re-scaled
fnscale = .checkFnscale(control=control)
control$fnscale = NULL # reset fnscale
# Checking par
active = .checkActive(active=active, npar=npar)
bounds = .checkBounds(lower=lower, upper=upper, npar=npar)
guess = .checkOpt(par=par, lower=bounds$lower, upper=bounds$upper)
par = guess
lower = bounds$lower
upper = bounds$upper
isActive = which(active)
activeFlag = isTRUE(all(active))
if(is.null(names(par))) names(par) = .printSeq(npar, preffix="par")
# getting functions
fn = match.fun(fn)
if(!is.null(gr)) warning("Gradient ignored by this method.")
# when 'parscale' and 'fnscale' are supplied,
# optimization is carried out over the scaled function
fn1 = function(par) {
parx = guess
parx[isActive] = par
parx = relist(flesh = parx*parscale, skeleton = skeleton)
output = fn(parx, ...)/fnscale
return(output)
}
# here, make methods explicit (one by one)
output = .ahres(par=par, fn=fn1, gr=NULL, lower=lower, upper=upper, control=control, hessian=hessian)
# reshaping full parameters, un-scaling par
paropt = guess
paropt[isActive] = output$par
paropt = paropt*parscale
output$par = paropt
# value need to be unscaled
output$value = output$value*fnscale
output = c(output, active=list(par=isActive, flag=activeFlag))
class(output) = c("ahres.result", class(output))
return(output)
}
# calibration_setup -------------------------------------------------------
#' Get information to run a calibration using the \code{calibrar} package.
#'
#' A wrapper for \code{read.csv} checking column names and data types
#' for the table with the calibration information.
#'
#' @param file The file with the calibration information, see details.
#' @param control Control arguments for generating the setup. See details.
#' @param \dots Additional arguments to \code{read.csv} function.
#' @return A data.frame with the information for the calibration of a
#' model, to be used with the \code{\link{calibration_objFn}}
#' and \code{\link{calibration_data}}.
#' @author Ricardo Oliveros-Ramos
#' @seealso \code{\link{calibration_objFn}}, \code{\link{calibration_data}}.
#' @export
calibration_setup = function(file, control=list(), ...) {
calibrationInfo = read.csv(file, stringsAsFactors=FALSE, ...)
fullNames = c("variable", "type", "calibrate", "cv", "weight", "use_data", "file", "varid", "col_skip", "nrows")
# mandatory columns
minNames = c("variable", "type", "calibrate", "use_data", "file")
doesNotMatch = !(names(calibrationInfo) %in% fullNames)
dnm = names(calibrationInfo)[doesNotMatch]
isMissing = !(minNames %in% names(calibrationInfo))
im = minNames[isMissing]
sdnm = if(length(dnm)>1) " columns do " else " column does "
sim = if(length(im)>1) " variables are " else " variable is "
msg1 = sprintf("Error in %s (%s %s not match).", file, paste(sapply(dnm, sQuote), collapse=", "), sdnm)
msg2 = sprintf("Error in %s (%s %s missing).", file, paste(sapply(im, sQuote), collapse=", "), sim)
if(any(doesNotMatch)) stop(msg1)
if(any(isMissing)) stop(msg2)
if(is.null(calibrationInfo$weight) & is.null(calibrationInfo$cv))
stop("Either 'cv' or 'weight' columns must be specified.")
# parsing correct data types
calibrationInfo$variable = as.character(calibrationInfo$variable)
calibrationInfo$type = as.character(calibrationInfo$type)
calibrationInfo$calibrate = as.logical(calibrationInfo$calibrate)
if(!is.null(calibrationInfo$weight)) calibrationInfo$weight = as.numeric(calibrationInfo$weight)
if(!is.null(calibrationInfo$cv)) calibrationInfo$cv = as.numeric(calibrationInfo$cv)
calibrationInfo$use_data = as.logical(calibrationInfo$use_data)
calibrationInfo$file = as.character(calibrationInfo$file)
calibrationInfo$varid = as.character(calibrationInfo$varid)
if(is.null(control$col_skip)) control$col_skip = 2
if(is.null(control$nrows)) control$nrows = NA
# optional columns
if(is.null(calibrationInfo$weight))
calibrationInfo$weight = 1/(2*calibrationInfo$cv^2)
if(any(is.na(calibrationInfo$weight))) {
if(is.null(calibrationInfo$cv)) stop("NAs found in 'weight' column.")
isna = is.na(calibrationInfo$weight)
calibrationInfo$weight[isna] = 1/(2*calibrationInfo$cv[isna]^2)
}
if(any(is.na(calibrationInfo$weight))) stop("NAs found in 'weight' column.")
if(any(calibrationInfo$weight <= 0)) stop("All weights must be positive.")
if(is.null(calibrationInfo$varid))
calibrationInfo$varid = NA
if(is.null(calibrationInfo$col_skip))
calibrationInfo$col_skip = control$col_skip
if(is.null(calibrationInfo$nrows))
calibrationInfo$nrows = control$nrows
return(calibrationInfo)
}
# calibration_data --------------------------------------------------------
#' @title Get observed data for the calibration of a model
#'
#' @description Create a list with the observed data with the
#' information provided by its main argument.
#'
#' @param setup A data.frame with the information about the calibration,
#' normally created with the \code{\link{calibration_setup}} function.
#' See details.
#' @param path Path to the directory to look up for the data. Paths in setup are considered relatives to this path.
#' @param file Optional file to save the created object (as an 'rds' file.)
#' @param verbose If TRUE, detailed messages of the process are printed.
#' @param \dots Additional arguments to \code{read.csv} function
#' to read the data files.
#' @return A list with the observed data needed for a calibration, to be used
#' in combination with the \code{\link{calibration_objFn}}.
#' @author Ricardo Oliveros-Ramos
#' @seealso \code{\link{calibration_objFn}}, \code{\link{calibration_setup}}.
#' @export
calibration_data = function(setup, path=".", file=NULL, verbose=TRUE, ...) {
observed = list()
useData = as.logical(setup$use_data)
isActive = as.logical(setup$calibrate)
message("Creating observed data list for calibration...","\n")
readData = useData & isActive
for(iVar in seq_len(nrow(setup))) {
ifile = file.path(path, setup$file[iVar])
varid = setup$varid[iVar]
col_skip = setup$col_skip[iVar]
nrows = setup$nrows[iVar]
if(readData[iVar]) {
observed[[iVar]] = .read_data(file=ifile, col_skip=col_skip, varid=varid, nrows=nrows, ...)
} else {
observed[[iVar]] = NA
}
} # end iVar loop
msg0 = "Loaded observed data for variable: %s.\n"
msg1 = "Loaded observed data for variables: %s.\n"
msg = if(sum(readData, na.rm=TRUE)==1) msg0 else msg1
if(isTRUE(verbose)) message(sprintf(msg, paste(sQuote(setup$variable[readData]), collapse=", ")))
names(observed) = as.character(setup$variable)
if(!is.null(file)) saveRDS(observed, file=file)
return(observed)
}
# calibration_objFn -------------------------------------------------------
#' Create an objective function to be used with optimization routines
#'
#' Create a new function, to be used as the objective function in the
#' calibration, given a function to run the model within R, observed data
#' and information about the comparison with data.
#'
#' @param model Function to run the model and produce a list of outputs.
#' @param setup A data.frame with the information about the calibration,
#' normally created with the \code{\link{calibration_setup}} function.
#' See details.
#' @param observed A list of the observed variables created with the
#' function \code{\link{calibration_data}}
#' @param aggFn A function to aggregate \code{fn} to a scalar value if the
#' returned value is a vector. Some optimization algorithm can explote the
#' additional information provided by a vectorial output from \code{fn}
#' @param aggregate boolean, if TRUE, a scalar value is returned using the
#' \code{aggFn}.
#' @param \dots More arguments passed to the \code{model} function.
#' @return A function, integrating the simulation of the model and the
#' comparison with observed data.
#' @author Ricardo Oliveros-Ramos
#' @seealso \code{\link{calibration_data}}, \code{\link{calibration_setup}}.
#' @export
calibration_objFn = function(model, setup, observed, aggFn=NULL, aggregate=FALSE, ...) {
fn_name = deparse(substitute(model))
fn = match.fun(model)
if(is.null(aggFn)) aggFn = .weighted.sum
aggFn = match.fun(aggFn)
force(observed)
force(setup)
force(aggregate)
weights = setup$weight[setup$calibrate]
msg = sprintf("The '%s' function returned NULL, please check.", fn_name)
# check for names in observed and simulated
fn1 = function(par) {
aggFn = match.fun(aggFn)
simulated = fn(par, ...)
# in case the model produce a NULL, keep moving.
if(is.null(simulated)) {
message(msg)
return(NULL)
}
output = .calculateObjetiveValue(obs=observed, sim=simulated, info=setup)
if(isTRUE(aggregate)) output = aggFn(x=output, w=weights)
return(output)
}
class(fn1) = c(class(fn1), "objFn")
fnx = function(par) {
simulated = fn(par, ...)
return(simulated)
}
attr(fn1, "nvar") = if(aggregate) 1 else sum(setup$calibrate)
attr(fn1, "weights") = weights
attr(fn1, "variables") = setup$variable[setup$calibrate]
attr(fn1, "aggregate") = aggregate
attr(fn1, "fn") = fnx
attr(fn1, "aggFn") = if(aggregate) NULL else aggFn
return(fn1)
}
multiMethods = "AHR-ES" # list of methods supporting multi-objective
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.