This page documents options to control fit_wrc_hcc. It
describes the arguments of the functions control_fit_wrc_hcc,
param_boundf, param_transf, fwd_transf,
dfwd_transf, bwd_transf, control_nloptr,
control_sce and control_pcmp, which all serve to steer
fit_wrc_hcc.
control_fit_wrc_hcc(
settings = c("uglobal", "ulocal", "clocal", "cglobal", "sce"),
method = c("ml", "mpd", "wls"), hessian,
nloptr = control_nloptr(), sce = control_sce(),
wrc_model = "vg", hcc_model = "vgm",
initial_param = c(alpha = 2., n = 1.5, tau = 0.5),
approximation_alpha_k0 =
c(c0 = 1, c1 = 5.55, c2 = 1.204, c3 = 2.11, c4 = 1.71),
variable_weight = c(wrc = 1, hcc = 1),
gam_k = 6, gam_n_newdata = 101, precBits = 256,
min_nobs_wc = 5, min_nobs_hc = 5,
keep_empty_fits = FALSE,
param_bound = param_boundf(), param_tf = param_transf(),
fwd_tf = fwd_transf(), deriv_fwd_tfd = dfwd_transf(), bwd_tf = bwd_transf(),
pcmp = control_pcmp())param_boundf(alpha = c(0 + 10 * sqrt(.Machine$double.eps), 500),
n = c(1 + 10 * sqrt(.Machine$double.eps), 20), tau = c(-2, 20),
thetar = c(0, 1), thetas = c(0, 1), k0 = c(0, Inf))
param_transf(alpha = c("log", "identity"), n = c("log1", "identity"),
tau = c("logitlu", "identity"), k0 = c("log", "identity"))
fwd_transf(...)
dfwd_transf(...)
bwd_transf(...)
control_nloptr(...)
control_sce(reltol = 1.e-8, maxtime = 20, ...)
control_pcmp(ncores = detectCores() - 1L,
fork = !identical(.Platform[["OS.type"]], "windows"))
control_fit_wrc_hcc creates a list with components
settings, hessian, method, nloptr,
sce, wrc_model, hcc_model, initial_param,
approximation_alpha_k0, variable_weight, gam_k,
gam_n_newdata, precBits, min_nobs_wc,
min_nobs_hc, keep_empty_fits, param_bound,
param_tf, fwd_tf, deriv_fwd_tfd, bwd_tf,
pcmp corresponding to its arguments and some further components
(delta_sat_0, grad_eps, use_derivative) that cannot
be changed by the user.
control_nloptr and control_sce
create lists with control parameters passed to
nloptr and SCEoptim,
respectively, see Details.
param_boundf generates a list with allowed lower and upper bounds of
the model parameters.
param_transf generates a list with keywords that
define what transformations are used for estimating the model
parameters, and fwd_transf, bwd_transf and
dfwd_transf return lists of functions with forward and backward
transformations and the first derivatives of the forward
transformations, see Details.
control_pcmp generates a list with control parameters for parallel
computations.
a keyword with possible values "uglobal"
(default), "ulocal", "clocal", "cglobal" or
"sce" to choose the approach for the nonlinear optimisation, see
Details and soilhypfitIntro.
a keyword with possible values "ml" (maximum
likelihood, default), "mpd" (maximum posterior density) or
"wls" (weighted least squares) to choose the method for estimating
the nonlinear parameters \(\boldsymbol{\nu}^\mathrm{} \), see
soilhypfitIntro.
a logical scalar controlling whether the Hessian matrix of
the objective function should be computed with respect to the possibly
transformed nonlinear parameters \(\boldsymbol{\nu}^\mathrm{} \) at the solution
(default: TRUE if settings %in% c("uglobal", "ulocal",
"sce") && method %in% c("mpd", "ml") and FALSE otherwise).
a list of arguments passed to the optimiser
nloptr by its argument opts, or a function
such as control_nloptr that generates such a list.
Note that
control_fit_wrc_hcc chooses sensible default values for the
various components of the list in dependence of the value chosen for
settings, but these defaults can be overridden by the arguments of
control_nloptr, see Details.
a list of arguments passed to the optimiser
SCEoptim by its argument control, or a
function such as control_sce that generates such a list.
Note
that control_fit_wrc_hcc chooses sensible default values for the
various components of the list, but these defaults can be overridden by
the arguments of control_sce, see Details.
a keyword choosing the parametrical model for the
water retention curve. Currently, only the Van Genuchten model
("vg") is implemented, see wc_model
and soilhypfitIntro.
a keyword choosing the parametrical model for the
hydraulic conductivity function. Currently, only the Van
Genuchten-Mualem model ("vgm") is implemented, see
hc_model and soilhypfitIntro.
a named numeric vector with default initial values
for the nonlinear parameters \(\boldsymbol{\nu}^\mathrm{} \)
of the models for the water retention curve and/or hydraulic conductivity
function. Currently, initial values must be defined for the parameters
\(\alpha\) (default 1.5 [\(\mathrm{m}^{-1}\)]), \(n\)
(default 2 [-]) and \(\tau\) (default 0.5 [-]), hence the elements of
initial_param must be named "alpha", "n" and
"tau".
a named numeric vector with constants to
approximate the parameter \(\alpha\) and the saturated hydraulic
conductivity \(K_0\) when constraining the estimated nonlinear
parameters \(\boldsymbol{\nu}^\mathrm{} \) of the Van
Genuchten-Mualem model by the characteristic evaporative length
(Lehmann et al., 2020), see evaporative-length and
soilhypfitIntro. For consistency with other quantities,
the following units should be used for the constants:
c1: \(\mathrm{m}^{-1}\),
c3: \(\mathrm{m}\,\mathrm{d}^{-1}\).
The remaining constants are dimensionless.
a named numeric vector of length 2 that defines
the weights of water content and hydraulic conductivity measurements in
the objective function for method = "wls". If equal to 1
(default) the weights of the variables are equal to the inverse variances
of the water content and (possibly log-transformed) hydraulic
conductivity measurements. If different from 1, then the inverse
variances are multiplied by variable_weight to get the variable
weights \(w_\theta\) and \(w_K\), see fit_wrc_hcc and
soilhypfitIntro. Note that for estimation methods mpd and
ml the variable weights are equal to 1 but the case weights
\(w^\prime_{\theta,i}\) and \(w^\prime_{K,j}\) may differ from 1.
the dimension of the basis of the additive model for
the water retention data when computing initial values of the
parameters \(\alpha\) and \(n\), see s and
soilhypfitIntro.
the number of evaluation points of the additive
model for the water retention data when computing initial values of the
parameters \(\alpha\) and \(n\), see soilhypfitIntro.
an integer scalar defining the default precision (in
bits) to be used in high-precision computations by
mpfr, see Details.
an integer scalar defining the minimum number of water content measurements per sample required for fitting a model to the water content data.
an integer scalar defining the minimum number of hydraulic conductivity measurements per sample required for fitting a model to hydraulic conductivity data.
a logical scalar controlling whether missing
fitting results should be dropped for samples without any data or failed
fits (FALSE, default) or kept (TRUE).
a named list of numeric vectors of length 2 that
define the allowed lower and upper bounds (box constraints) for the
parameters of the models or a function such as param_boundf which
generates this list, see Details.
a named vector of keywords that define the
transformations to be applied to the model parameters before estimation
or a function such as param_transf, which generates this vector,
see Details.
a named list of invertible functions to be used to
transform model parameters or a function such as fwd_transf, which
generates this list, see Details.
a named list of functions corresponding to the first
derivatives of the functions defined in fwd_tf or a function such
as dfwd_transf, which generates this list, see Details.
a named list of inverse functions corresponding the
functions defined in fwd_tf or a function such as bwd_transf,
which generates this list, see Details.
a list to control parallel computations or a function such as
control_pcmp that generates this list, see control_pcmp for
allowed arguments.
either a numeric vector of length 2 defining the allowed
lower and upper bounds for the nonlinear parameter \(\alpha\)
(param_boundf), or a keyword defining the transformation
to be applied to the parameter \(\alpha\) before estimation
(param_transf), see Details.
either a numeric vector of length 2 defining the allowed lower
and upper bounds for the nonlinear parameter \(n\)
(param_boundf), or a keyword defining the transformation
to be applied to the parameter \(n\) before estimation
(param_transf), see Details.
either a numeric vector of length 2 defining the allowed
lower and upper bounds for the nonlinear parameter \(\tau\)
(param_boundf), or a keyword defining the transformation
to be applied to the parameter \(\tau\) before estimation
(param_transf), see Details.
a numeric vector of length 2 defining the allowed
lower and upper bounds for the linear parameter \(\theta_r\)
(param_boundf), see Details.
a numeric vector of length 2 defining the allowed
lower and upper bounds for the linear parameter \(\theta_s\)
(param_boundf), see Details.
either a numeric vector of length 2 defining the allowed
lower and upper bounds for the linear parameter \(K_0\)
(param_boundf), or a keyword defining the transformation
to be applied to the parameter \(K_0\) before estimation
(param_transf), see Details.
a numeric scalar defining (one possible) convergence
criterion for the optimiser SCEoptim, see argument control
of SCEoptim for details.
a numeric scalar defining the maximum duration of
optimisation in seconds by the optimiser SCEoptim, see see
argument control of SCEoptim for details.
an integer defining the number of cores for
parallel computations. Defaults to the number of available cores
minus one. ncores = 1 suppresses parallel computations.
a logical scalar controlling whether forking should be used
for parallel computations (default: TRUE on Unix and MacOS and
FALSE on Windows operating systems). Note that stetting
fork = TRUE on Windows suppresses parallel computations.
further arguments, such as details on parameter
transformations (fwd_transf, dfwd_transf,
bwd_transf) or control options passed to the
optimisers nloptr and
SCEoptim, see Details.
Andreas Papritz papritz@retired.ethz.ch.
Parameters of models for the water retention curve and the hydraulic
conductivity function may vary only within certain bounds (see
param_boundf for allowed ranges). fit_wrc_hcc
uses two mechanisms to constrain parameter estimates to permissible
ranges:
Parameter transformations
If a local algorithm is used for nonlinear optimisation
(settings = "ulocal" or settings = "clocal") and a
transformation not equal to "identity" is specified in
param_tf for any of the nonlinear parameters
\(\boldsymbol{\nu}^\mathrm{} \), then the elements of \(\boldsymbol{\nu}^\mathrm{} \) are
transformed by the functions given in param_tf. The values
of the transformed parameters vary then over the whole real line,
and an unconstrained algorithm can be used for nonlinear
optimisation.
Note that the linear parameters \(\theta_r\) (residual)
and \(\theta_s\) (saturated water content) are never transformed
and for the saturated hydraulic conductivity, \(K_0\), only
"log" (default) or "identity" can be chosen.
Quadratic programming (see solve.QP) is
employed to enforce the box constraints specified in the argument
param_bound for \(\theta_r\) and \(\theta_s\).
Quadratic programming is also used to enforce the positivity
constraint on \(K_0\) if \(K_0\) is not log-transformed
("identity"). Otherwise, the logarithm of \(K_0\) is
estimated unconstrainedly, see soilhypfitIntro for
further details.
Box constraints
If a global algorithm is used for the optimisation (settings
equal to "uglobal" "cglobal" or "sce") or if
"identity" transformations are specified for all elements of
\(\boldsymbol{\nu}^\mathrm{} \), then an optimisation algorithm is deployed that
respects the box constraints given in param_bound. If
parameters are estimated for several soil samples in a single call
of fit_wrc_hcc and if sample-specific box constraints
should be used then the lower and upper bounds of the
box-constraints must be specified by the arguments
lower_param and upper_param of the function
fit_wrc_hcc, see explanations there.
Further note that the transformations specified by param_tf
for the nonlinear parameters \(\boldsymbol{\nu}^\mathrm{} \) are ignored when a
global optimisation algorithm is used.
The arguments param_tf, fwd_tf, deriv_fwd_tfd,
bwd_tf define how the model parameters are transformed for
estimation by local optimisation algortihms (see above and
soilhypfitIntro). The following transformations are
currently available:
"log":\(\log(x)\),
"log1":\(\log(x-1)\),
"logitlu":\(\log((x - l) / (u - x))\) with \(l\)
and \(u\) the allowed lower and upper bounds for a parameter, see
param_boundf,
"identity":no transformation.
These are the possible values that the various arguments of the
function param_transf accept (as quoted character strings), and
these are the names of the list components returned by
fwd_transf, dfwd_transf and bwd_transf.
Additional transformations can be implemented by:
Extending the function definitions by arguments like
fwd_tf = fwd_transf(my_fun = function(x) your transformation),
deriv_fwd_tfd = dfwd_transf(my_fun = function(x) your derivative),
bwd_tf = bwd_transf(my_fun = function(x) your back-transformation),
Assigning to a given argument of param_transf the name
of the new function, e.g.
alpha = "my_fun".
Note that the values given to the arguments of param_transf must
match the names of the functions returned by fwd_transf,
dfwd_transf and bwd_transf.
Estimation of \(\log(K_0)\) is somewhat delicate for large
values of the shape parameter \(n\) and/or small values of
\(\alpha\). The water saturation and the relative conductivity are
then close to zero for capillary pressure head exceeding
\(1/\alpha\). To avoid numerical problems caused by limited accuracy
of double precision numbers, fit_wrc_hcc uses the
function mpfr of the package Rmpfr for
high-accuracy computations. The argument precBits of
control_fit_wrc_hcc controls the accuracy. Increase its
value if computation fails with a respective diagnostic message.
The argument settings defines sets of default options to control
the optimisers. The following settings are currently defined:
"uglobal":unconstrained optimisation by any of the
global algorithms (named "NLOPT_G...") of the NLopt library.
"cglobal":constrained optimisation by the global
algorithm "NLOPT_GN_ISRES" of NLopt that allows for inequality
constraints.
"ulocal":unconstrained optimisation by any of the
local algorithms
(named "NLOPT_L...") of
NLopt.
"clocal":constrained optimisation by any of the local
algorithms
("NLOPT_LN_COBYLA", "NLOPT_LN_AUGLAG",
"NLOPT_LD_AUGLAG", "NLOPT_LD_SLSQP",
"NLOPT_LD_MMA"), "NLOPT_LD_CCSAQ") of NLopt that allow
for inequality constraints.
"sce":unconstrained optimisation by the global
algorithm implemented in SCEoptim.
The functions control_nloptr and control_sce allow finer
control of the optimisers.
control_nloptr and control_sce take any
argument available to pass controlling options to the optimisers
nloptr (by its argument opts) and
SCEoptim (by its argument control),
respectively.
The function nloptr.print.options prints all
options available to control nloptr by its
argument opts. Detailed information on the options can be
found in the
NLopt
documentation.
The function control_fit_wrc_hcc sets meaningful defaults for
opts in dependence of the chosen optimisation approach as
specified by the argument settings, and it checks the
consistency of the arguments of control_nloptr if they are
explicitly passed to fit_wrc_hcc.
The following defaults are set by control_fit_wrc_hcc for the
argument opts of nloptr (:
Unconstrained, global optimisation (settings =
"uglobal"):
nloptr = control_nloptr(
algorithm = "NLOPT_GN_MLSL_LDS",
local_opts = list(
algorithm = "NLOPT_LN_BOBYQA",
xtol_rel = -1.,
ftol_rel = 1.e-6
),
xtol_rel = -1,
ftol_rel = -1,
maxeval = 125,
maxtime = -1)
In addition, any parameter transformations specified by
param_tf are overridden and the untransformed parameters
("identity") are estimated when settings = "uglobal" is
chosen.
Constrained, global optimisation (settings =
"cglobal"):
nloptr = control_nloptr(
algorithm = "NLOPT_GN_ISRES",
xtol_rel = -1,
ftol_rel = -1,
maxeval = 1000,
maxtime = -1)
In addition, any parameter transformations specified by
param_tf are overridden and the untransformed parameters
("identity") are estimated when settings = "cglobal"
is chosen.
Unconstrained, local optimisation (settings =
"ulocal"):
nloptr = control_nloptr(
algorithm = "NLOPT_LN_BOBYQA",
xtol_rel = -1,
ftol_rel = 1.e-8,
maxeval = 250,
maxtime = -1)
Constrained, local optimisation (settings = "clocal"):
nloptr = control_nloptr(
algorithm = "NLOPT_LD_CCSAQ",
xtol_rel = -1,
ftol_rel = 1.e-8,
maxeval = 1000,
maxtime = -1)
If the algorithm "NLOPT_LD_AUGLAG" is used for constrained,
local optimisation then
nloptr = control_nloptr(
algorithm = "NLOPT_LD_AUGLAG",
local_opts = list(
algorithm = "NLOPT_LD_LBFGS",
xtol_rel = -1.,
ftol_rel = 1.e-6
),
xtol_rel = -1,
ftol_rel = 1.e-8,
maxeval = 1000,
maxtime = -1)
For other, unspecified elements of opts default values as
listed by nloptr.print.options are used.
The function control_sce sets meaningful defaults for the
argument control of SCEoptim.
Currently, the following defaults are defined:
sce = control_sce(
reltol = 1e-08,
maxtime = 20)
In addition, any parameter transformations specified by
param_tf are overridden and the untransformed parameters
("identity") are estimated when settings = "sce" is
chosen.
Johnson, S.G. The NLopt nonlinear-optimisation package. https://github.com/stevengj/nlopt.
Lehmann, P., Assouline, S., Or, D. (2008) Characteristic lengths affecting evaporative drying of porous media. Physical Review E, 77, 056309, tools:::Rd_expr_doi("10.1103/PhysRevE.77.056309").
Lehmann, P., Bickel, S., Wei, Z., Or, D. (2020) Physical Constraints for Improved Soil Hydraulic Parameter Estimation by Pedotransfer Functions. Water Resources Research 56, e2019WR025963, tools:::Rd_expr_doi("10.1029/2019WR025963").
soilhypfitIntro for a description of the models and a brief
summary of the parameter estimation approach;
fit_wrc_hcc for (constrained) estimation of parameters of
models for soil water retention and hydraulic conductivity data;
soilhypfitmethods for common S3 methods for class
fit_wrc_hcc;
vcov for computing (co-)variances of the estimated
nonlinear parameters;
prfloglik_sample for profile loglikelihood
computations;
wc_model and hc_model for currently
implemented models for soil water retention curves and hydraulic
conductivity functions;
evaporative-length for physically constraining parameter
estimates of soil hydraulic material functions.
# \donttest{
# use of \donttest{} because execution time exceeds 5 seconds
data(sim_wrc_hcc)
# estimate parameters for a single soil sample by maximizing loglikelihood ...
# ... with unconstrained, global optimisation algorithm NLOPT_GN_MLSL
coef(
fit1 <- fit_wrc_hcc(
wrc_formula = wc ~ head, hcc_formula = hc ~ head,
data = sim_wrc_hcc, subset = id == 2
), gof = TRUE)
# ... as fit1 but fitting parameter tau as well
coef(
fit2 <- update(fit1,
fit_param = default_fit_param(tau = TRUE)
), gof = TRUE)
plot(fit1, y = fit2)
# ... with unconstrained, local optimisation algorithm NLOPT_LN_BOBYQA,
# initial values for alpha and n are computed from data and
# transformed nonlinear parameters are estimated without box-constraints
coef(
fit3 <- update(
fit2,
control = control_fit_wrc_hcc(settings = "ulocal"),
verbose = 2), gof = TRUE)
# estimate parameters by unconstrained, weighted least squares minimisation with
# algorithm NLOPT_LD_LBFGS, giving larger weight to conductivity data,
# using specified initial values for alpha and n and
# fitting untransformed nonlinear parameters with default box constraints
# defined by param_boundf()
# diagnostic output directly from nloptr
coef(
fit4 <- update(
fit2,
param = c(alpha = 1.7, n = 2),
control = control_fit_wrc_hcc(
settings = "ulocal", method = "wls",
variable_weight = c(wrc = 1, hcc = 2),
nloptr = control_nloptr(algorithm = "NLOPT_LD_LBFGS", print_level = 3),
param_tf = param_transf(alpha = "identity", n = "identity", tau = "identity")
), verbose = 0), gof = TRUE)
# ... as fit4 but giving even larger weight to conductivity data
coef(
fit5 <- update(
fit4,
control = control_fit_wrc_hcc(
settings = "ulocal", method = "wls",
variable_weight = c(wrc = 1, hcc = 5),
nloptr = control_nloptr(algorithm = "NLOPT_LD_LBFGS", print_level = 3),
param_tf = param_transf(alpha = "identity", n = "identity", tau = "identity")
), verbose = 0), gof = TRUE)
plot(fit4, y = fit5)
# }Run the code above in your browser using DataLab