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 "uglobal" (default), "ulocal", "clocal", "cglobal" or "sce" to choose the approach for the nonlinear optimisation, see Details and soilhypfitIntro.

method

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.

hessian

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).

nloptr

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.

sce

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.

wrc_model

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.

hcc_model

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.

initial_param

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".

approximation_alpha_k0

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.

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 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.

gam_k

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.

gam_n_newdata

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.

precBits

an integer scalar defining the default precision (in bits) to be used in high-precision computations by mpfr, see Details.

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 (FALSE, default) or kept (TRUE).

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_boundf which generates this list, see Details.

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 param_transf, which generates this vector, see Details.

fwd_tf

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.

deriv_fwd_tfd

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.

bwd_tf

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.

pcmp

a list to control parallel computations or a function such as control_pcmp that generates this list, see control_pcmp for allowed arguments.

alpha

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.

n

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.

tau

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.

thetar

a numeric vector of length 2 defining the allowed lower and upper bounds for the linear parameter \theta_r (param_boundf), see Details.

thetas

a numeric vector of length 2 defining the allowed lower and upper bounds for the linear parameter \theta_s (param_boundf), see Details.

k0

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.

reltol

a numeric scalar defining (one possible) convergence criterion for the optimiser SCEoptim, see argument control of SCEoptim for details.

maxtime

a numeric scalar defining the maximum duration of optimisation in seconds by the optimiser SCEoptim, see see argument control of SCEoptim for details.

ncores

an integer defining the number of cores for parallel computations. Defaults to the number of available cores minus one. ncores = 1 suppresses parallel computations.

fork

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.

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:

  1. 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.

  2. 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.

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":

\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:

  1. 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),

  2. 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 \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.

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 (:

  1. 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.

  2. 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.

  3. Unconstrained, local optimisation (settings = "ulocal"): 

    nloptr = control_nloptr(
  algorithm = "NLOPT_LN_BOBYQA",
  xtol_rel = -1,
  ftol_rel = 1.e-8,
  maxeval = 250,
  maxtime = -1)

  4. 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.

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)
[Package soilhypfit version 0.1-7 Index]