control_fit_wrc_hcc {soilhypfit} | R Documentation |
Controlling fit_wrc_hcc
Description
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
.
Usage
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"))
Arguments
settings |
a keyword with possible values |
method |
a keyword with possible values |
hessian |
a logical scalar controlling whether the Hessian matrix of
the objective function should be computed with respect to the possibly
transformed nonlinear parameters |
nloptr |
a list of arguments passed to the optimiser
|
sce |
a list of arguments passed to the optimiser
|
wrc_model |
a keyword choosing the parametrical model for the
water retention curve. Currently, only the Van Genuchten model
( |
hcc_model |
a keyword choosing the parametrical model for the
hydraulic conductivity function. Currently, only the Van
Genuchten-Mualem model ( |
initial_param |
a named numeric vector with default initial values
for the nonlinear parameters |
approximation_alpha_k0 |
a named numeric vector with constants to
approximate the parameter
The remaining constants are dimensionless. |
variable_weight |
a named numeric vector of length 2 that defines
the weights of water content and hydraulic conductivity measurements in
the objective function for |
gam_k |
the dimension of the basis of the additive model for
the water retention data when computing initial values of the
parameters |
gam_n_newdata |
the number of evaluation points of the additive
model for the water retention data when computing initial values of the
parameters |
precBits |
an integer scalar defining the default precision (in
bits) to be used in high-precision computations by
|
min_nobs_wc |
an integer scalar defining the minimum number of water content measurements per sample required for fitting a model to the water content data. |
min_nobs_hc |
an integer scalar defining the minimum number of hydraulic conductivity measurements per sample required for fitting a model to hydraulic conductivity data. |
keep_empty_fits |
a logical scalar controlling whether missing
fitting results should be dropped for samples without any data or failed
fits ( |
param_bound |
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_tf |
a named vector of keywords that define the
transformations to be applied to the model parameters before estimation
or a function such as |
fwd_tf |
a named list of invertible functions to be used to
transform model parameters or a function such as |
deriv_fwd_tfd |
a named list of functions corresponding to the first
derivatives of the functions defined in |
bwd_tf |
a named list of inverse functions corresponding the
functions defined in |
pcmp |
a list to control parallel computations or a function such as
|
alpha |
either a numeric vector of length 2 defining the allowed
lower and upper bounds for the nonlinear parameter |
n |
either a numeric vector of length 2 defining the allowed lower
and upper bounds for the nonlinear parameter |
tau |
either a numeric vector of length 2 defining the allowed
lower and upper bounds for the nonlinear parameter |
thetar |
a numeric vector of length 2 defining the allowed
lower and upper bounds for the linear parameter |
thetas |
a numeric vector of length 2 defining the allowed
lower and upper bounds for the linear parameter |
k0 |
either a numeric vector of length 2 defining the allowed
lower and upper bounds for the linear parameter |
reltol |
a numeric scalar defining (one possible) convergence
criterion for the optimiser |
maxtime |
a numeric scalar defining the maximum duration of
optimisation in seconds by the optimiser |
ncores |
an integer defining the number of cores for
parallel computations. Defaults to the number of available cores
minus one. |
fork |
a logical scalar controlling whether forking should be used
for parallel computations (default: |
... |
further arguments, such as details on parameter
transformations ( |
Details
Enforcing bounds on the estimated parameters
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"
orsettings = "clocal"
) and a transformation not equal to"identity"
is specified inparam_tf
for any of the nonlinear parameters, then the elements of
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
(residual) and
(saturated water content) are never transformed and for the saturated hydraulic conductivity,
, only
"log"
(default) or"identity"
can be chosen. Quadratic programming (seesolve.QP
) is employed to enforce the box constraints specified in the argumentparam_bound
forand
. Quadratic programming is also used to enforce the positivity constraint on
if
is not log-transformed (
"identity"
). Otherwise, the logarithm ofis 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, 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 offit_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 argumentslower_param
andupper_param
of the functionfit_wrc_hcc
, see explanations there.Further note that the transformations specified by
param_tf
for the nonlinear parametersare ignored when a global optimisation algorithm is used.
Parameter transformations
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"
:,
"log1"
:,
"logitlu"
:with
and
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
.
High-precision numerical computations
Estimation of is somewhat delicate for large
values of the shape parameter
and/or small values of
. The water saturation and the relative conductivity are
then close to zero for capillary pressure head exceeding
. 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.
Options to choose the approach for nonlinear optimisation
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.
Controlling nloptr
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 whensettings = "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 whensettings = "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 thennloptr = 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.
Controlling SCEoptim
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.
Value
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.
Author(s)
Andreas Papritz papritz@retired.ethz.ch.
References
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, 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, doi:10.1029/2019WR025963.
See Also
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.
Examples
# 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)