Package 'calibrar'

Automated Calibration for Complex Models

Description

Automated Calibration for Complex Models

Details

calibrar package: Automated Calibration for Complex Models

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.

Author(s)

Ricardo Oliveros-Ramos Maintainer: Ricardo Oliveros-Ramos <[email protected]>

References

calibrar: an R package for the calibration of ecological models (Oliveros-Ramos and Shin 2014)

See Also

Useful links:

• https://roliveros-ramos.github.io/calibrar/

• Report bugs at https://github.com/roliveros-ramos/calibrar/issues

Examples

## Not run: 
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)

## End(Not run)

---

Get an specific argument from the command line

Description

Get an specific argument from the command line

Usage

.get_command_argument(
  x,
  argument,
  prefix = "--",
  default = FALSE,
  verbose = FALSE
)

Arguments

x	The command line arguments, from x = commandArgs()
argument	The name of the argument.
prefix	The prefix to any argument of interest, the default is "–"
default	Default value to return is argument is missing, default to FALSE.
verbose	Boolean, if TRUE, shows a warning when the parameter is not found.

Value

The value of the argument, assumed to be followed after '=' or, TRUE if nothing but the argument was found. If the argument is not found, FALSE is returned.

Examples

.get_command_argument(commandArgs(), "interactive")
.get_command_argument(commandArgs(), "RStudio")
.get_command_argument(commandArgs(), "RStudio", prefix="")
.get_command_argument(commandArgs(), "vanilla")
.get_command_argument("--control.file=baz.txt", "control.file")

---

Read a configuration file.

Description

File is expected to have lines of the form 'key SEP value' where key is the name of the parameter, SEP a separator (can be '=' ',', ';') and value the value of the parameter itself. The SEP for each line is determined and parameters values are returned as a list.

Usage

.read_configuration(
  file,
  recursive = TRUE,
  keep.names = TRUE,
  conf.key = NULL,
  ...
)

Arguments

file	File to read the configuration
recursive	Should 'conf.key' keys be read as additional configuration files? Default is TRUE.
keep.names	Should names be kept as they are? By default, are converted to lower case.
conf.key	String indicating the leading key to find an additional configuration file.
...	Additional arguments, not currently in use.

---

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

Usage

ahres(
  par,
  fn,
  gr = NULL,
  ...,
  lower = -Inf,
  upper = +Inf,
  active = NULL,
  control = list(),
  hessian = FALSE,
  parallel = FALSE
)

Arguments

par	A numeric vector or list. The length of par defines the number of parameters to be estimated (i.e., the dimension of the problem).
fn	The objective function to be minimised. It should accept a parameter vector (or list, depending on the wrapper) as first argument and return either a scalar value (single-objective) or a numeric vector (multi-objective).
gr	A function computing the gradient of fn. If NULL, a numerical approximation is used. Alternatively, a character string can specify the numerical gradient scheme: "central", "forward" (default), "backward", or "richardson".
...	Additional arguments passed to fn and gr.
lower	Lower bounds for parameters. One value or a vector of the same length as par. NA is treated as -Inf. Default is unconstrained.
upper	Upper bounds for parameters. One value or a vector of the same length as par. NA is treated as Inf. Default is unconstrained.
active	Boolean vector of the same length as par, indicating if the parameter is used in the optimisation (TRUE) or held at a fixed value (FALSE).
control	A list of control options. Common options include ncores, run, master, verbose, REPORT, restart.file, gradient, and gr.method. Additional solver-specific options may be passed through to the underlying optimiser.
hessian	Logical. Should a numerically differentiated Hessian matrix be returned? Currently not implemented.
parallel	Logical. Enable parallel computation (e.g., for replicated evaluations and/or numerical gradients) using up to control$ncores cores.

Value

A list with components:

par

The best set of parameters found.

value

The value of fn corresponding to par.

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.

convergence

An integer code. 0 indicates successful completion.

message

A character string giving any additional information returned by the optimizer, or NULL.

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.

Author(s)

Ricardo Oliveros-Ramos

See Also

Other optimisers: calibrate(), optim2(), optimh()

Examples

## Not run: ahres(par=rep(1, 5), fn=sphereN)

---

Demos for the calibrar package

Description

Creates demo files able to be processed for a full calibration using the calibrar package

Usage

calibrar_demo(path = NULL, model = NULL, ...)

Arguments

path	Path to create the demo files
model	Model to be used in the demo files, see details.
...	Additional parameters to be used in the construction of the demo files.

Details

Current implemented models are:

PoissonMixedModel

Poisson Autoregressive Mixed model for the dynamics of a population in different sites:

$log(\mu_{i, t+1}) = log(\mu_{i, t}) + \alpha + \beta X_{i, t} + \gamma_t$

where $\mu_{i, t}$ is the size of the population in site $i$ at year $t$, $X_{i, t}$ is the value of an environmental variable in site $i$ at year $t$. The parameters to estimate were $\alpha$, $\beta$, and $\gamma_t$, the random effects for each year, $\gamma_t \sim N(0,\sigma^2)$, and the initial population at each site $\mu_{i, 0}$. We assumed that the observations $N_{i,t}$ follow a Poisson distribution with mean $\mu_{i, t}$.

PredatorPrey

Lotka Volterra Predator-Prey model. The model is defined by a system of ordinary differential equations for the abundance of prey $N$ and predator $P$:

$\frac{dN}{dt} = rN(1-N/K)-\alpha NP$

$\frac{dP}{dt} = -lP + \gamma\alpha NP$

The parameters to estimate are the prey’s growth rate $r$, the predator’s mortality rate $l$, the carrying capacity of the prey $K$ and $\alpha$ and $\gamma$ for the predation interaction. Uses deSolve package for numerical solution of the ODE system.

SIR

Susceptible-Infected-Recovered epidemiological model. The model is defined by a system of ordinary differential equations for the number of susceptible $S$, infected $I$ and recovered $R$ individuals:

$\frac{dS}{dt} = -\beta S I/N$

$\frac{dI}{dt} = \beta S I/N -\gamma I$

$\frac{dR}{dt} = \gamma I$

The parameters to estimate are the average number of contacts per person per time $\beta$ and the instant probability of an infectious individual recovering $\gamma$. Uses deSolve package for numerical solution of the ODE system.

IBMLotkaVolterra

Stochastic Individual Based Model for Lotka-Volterra model. Uses ibm package for the simulation.

Value

A list with the following elements:

path	Path were the files were saved
par	Real value of the parameters used in the demo
setup	Path to the calibration setup file
guess	Values to be provided as initial guess to the calibrate function
lower	Values to be provided as lower bounds to the calibrate function
upper	Values to be provided as upper bounds to the calibrate function
phase	Values to be provided as phases to the calibrate function
constants	Constants used in the demo, any other variable not listed here.
value	NA, set for compatibility with summary methods.
time	NA, set for compatibility with summary methods.
counts	NA, set for compatibility with summary methods.

Author(s)

Ricardo Oliveros–Ramos

References

Oliveros-Ramos and Shin (2014)

Examples

## Not run: 

summary(ahr)
set.seed(880820)
path = NULL # NULL to use the current directory
# create the demonstration files
demo = calibrar_demo(path=path, model="PredatorPrey", T=100) 
# get calibration information
calibration_settings = calibration_setup(file = demo$setup)
# get observed data
observed = calibration_data(setup = calibration_settings, path=demo$path)
# Defining 'run_model' function
run_model = calibrar:::.PredatorPreyModel
# real parameters
cat("Real parameters used to simulate data\n")
print(unlist(demo$par)) # parameters are in a list
# objective functions
obj  = calibration_objFn(model=run_model, setup=calibration_settings, observed=observed, T=demo$T)
obj2 = calibration_objFn(model=run_model, setup=calibration_settings, observed=observed, 
T=demo$T, aggregate=TRUE)
cat("Starting calibration...\n")
cat("Running optimization algorithms\n", "\t")
cat("Running optim AHR-ES\n")
ahr = calibrate(par=demo$guess, fn=obj, lower=demo$lower, upper=demo$upper, phases=demo$phase)
summary(ahr)

## End(Not run)

---

Sequential parameter estimation for the calibration of complex models

Description

calibrate() minimises an objective function fn for the calibration of complex (possibly expensive and stochastic) models. It supports sequential calibration in multiple phases (progressively activating parameters), replicated evaluations for stochastic objectives, restartable runs, and parallel execution patterns for computationally intensive models.

Usage

calibrate(
  par,
  fn,
  gr,
  ...,
  method,
  lower,
  upper,
  phases,
  control,
  hessian,
  replicates,
  parallel
)

## Default S3 method:
calibrate(
  par,
  fn,
  gr = NULL,
  ...,
  method = NULL,
  lower = NULL,
  upper = NULL,
  phases = NULL,
  control = list(),
  hessian = FALSE,
  replicates = 1,
  parallel = FALSE
)

Arguments

par	A numeric vector or list. The length of par defines the number of parameters to be estimated (i.e., the dimension of the problem).
fn	The objective function to be minimised. It should accept a parameter vector (or list, depending on the wrapper) as first argument and return either a scalar value (single-objective) or a numeric vector (multi-objective).
gr	A function computing the gradient of fn. If NULL, a numerical approximation is used. Alternatively, a character string can specify the numerical gradient scheme: "central", "forward" (default), "backward", or "richardson".
...	Additional arguments passed to fn and gr.
method	Optimisation method(s) to be used. Can be a single method name or a vector of method names (e.g., one per phase). If NULL, a default is chosen based on replicates (see Details).
lower	Lower bounds for parameters. One value or a vector of the same length as par. NA is treated as -Inf. Default is unconstrained.
upper	Upper bounds for parameters. One value or a vector of the same length as par. NA is treated as Inf. Default is unconstrained.
phases	Optional integer vector of the same length as par, indicating the phase at which each parameter becomes active. If omitted, all parameters are active in a single phase.
control	A list of control options. Common options include ncores, run, master, verbose, REPORT, restart.file, gradient, and gr.method. Additional solver-specific options may be passed through to the underlying optimiser.
hessian	Logical. Should a numerically differentiated Hessian matrix be returned? Currently not implemented.
replicates	Integer or integer vector controlling the number of replicate evaluations of fn per phase. The default is 1 (deterministic objective).
parallel	Logical. Enable parallel computation (e.g., for replicated evaluations and/or numerical gradients) using up to control$ncores cores.

Details

Sequential phases. When phases is provided, parameters are activated progressively across phases. In phase $k$, parameters with phases <= k are active and estimated, while all others are held fixed. If phases is omitted, all parameters are estimated in a single phase.

Replicated evaluations. The argument replicates controls the number of replicate evaluations of fn per phase. It can be a single integer (applied to all phases) or a vector of length equal to the number of phases. Replication is useful for stochastic models to reduce Monte Carlo noise and/or to target robust parameter sets.

Default method selection. If method is not provided, calibrate() defaults to "Rvmmin" when all(replicates == 1) (deterministic objective), and to "AHR-ES" otherwise (replicated/stochastic objective).

Multi-objective outputs. The objective function fn may return a scalar (single-objective) or a numeric vector of objective components (multi-objective). Multi-objective optimisation is currently supported only by methods in multiMethods (e.g., "AHR-ES"). For single-objective methods, fn must be scalar; alternatively, objective components can be aggregated into a scalar by defining an aggregated objective in calibration_objFn(aggregate = TRUE).

Aggregation. When fn is created via calibration_objFn(), metadata such as the number of components (nvar) and component weights (weights) can be stored as attributes. If aggregate = TRUE, fn is scalar and can be used with any method. If aggregate = FALSE (vector output), only multi-objective methods can use it directly.

Parallel execution and run directories. If parallel = TRUE, calibrate() may distribute replicate evaluations and/or finite-difference gradient computations across multiple cores. The number of cores is controlled by control$ncores. For file-based or externally executed models, control$master and control$run can be used to manage a master/template directory and per-run working directories (created if needed).

Restart. Partial results can be written and used to resume a calibration via control$restart.file. Restart functionality is currently available only for "AHR-ES", "Rvmmin", and "hjn".

Value

An object of class "calibrar.results" with components:

par

Best parameter values found (returned with the same structure as input par).

value

Objective value at par.

counts

Number of calls to fn and gr (where applicable).

convergence

Convergence code returned by the underlying optimiser.

message

Additional information returned by the optimiser, if any.

method

The optimisation method used in the final phase.

fn

The objective function.

active

Logical vector indicating which parameters were active (estimated).

elapsed

Elapsed time for the full calibration run.

trace

Tracing information and per-phase outputs (when available).

Choosing an optimisation method

For smooth deterministic objectives, gradient-based methods such as "Rvmmin", "L-BFGS-B", "LBFGSB3", or "nlminb" are often efficient (with box constraints when needed). For noisy or stochastic objectives (typically when replicates > 1), heuristic/global methods such as "AHR-ES" are generally more appropriate. Not all methods are directly comparable (local vs.\ global, deterministic vs.\ stochastic, scalar vs.\ vector-valued objectives); method choice should reflect the structure of fn and the presence of constraints.

Notes

"SANN" is included for compatibility with stats::optim(), but it is highly sensitive to tuning and often performs poorly on continuous problems under default settings. For stochastic objectives, "AHR-ES" is the intended default within calibrate().

Author(s)

Ricardo Oliveros-Ramos

See Also

calibration_setup, calibration_data, calibration_objFn

Other optimisers: ahres(), optim2(), optimh()

Examples

calibrate(par=rep(NA, 5), fn=sphereN)
## Not run: 
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))

## End(Not run)

---

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.

Usage

calibration_data(setup, path = ".", file = NULL, verbose = TRUE, ...)

Arguments

setup	A data.frame with the information about the calibration, normally created with the calibration_setup function. See details.
path	Path to the directory to look up for the data. Paths in setup are considered relatives to this path.
file	Optional file to save the created object (as an 'rds' file.)
verbose	If TRUE, detailed messages of the process are printed.
...	Additional arguments to read.csv function to read the data files.

Value

A list with the observed data needed for a calibration, to be used in combination with the calibration_objFn.

Author(s)

Ricardo Oliveros-Ramos

See Also

calibration_objFn, calibration_setup.

---

Create an objective function to be used with optimization routines

Description

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.

Usage

calibration_objFn(model, setup, observed, aggFn = NULL, aggregate = FALSE, ...)

Arguments

model	Function to run the model and produce a list of outputs.
setup	A data.frame with the information about the calibration, normally created with the calibration_setup function. See details.
observed	A list of the observed variables created with the function calibration_data
aggFn	A function to aggregate 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 fn
aggregate	boolean, if TRUE, a scalar value is returned using the aggFn.
...	More arguments passed to the model function.

Value

A function, integrating the simulation of the model and the comparison with observed data.

Author(s)

Ricardo Oliveros-Ramos

See Also

calibration_data, calibration_setup.

---

Get information to run a calibration using the calibrar package.

Description

A wrapper for read.csv checking column names and data types for the table with the calibration information.

Usage

calibration_setup(file, control = list(), ...)

Arguments

file	The file with the calibration information, see details.
control	Control arguments for generating the setup. See details.
...	Additional arguments to read.csv function.

Value

A data.frame with the information for the calibration of a model, to be used with the calibration_objFn and calibration_data.

Author(s)

Ricardo Oliveros-Ramos

See Also

calibration_objFn, calibration_data.

---

Calculate a discretization of the 2D Gaussian Kernel

Description

Calculate a discretization of the 2D Gaussian Kernel

Usage

gaussian_kernel(par, lower, upper, n = 10, checkSymmetry = TRUE, ...)

Arguments

par	A list, including the mean and covariance matrix.
lower	A vector, indicating the lower bound for the calculation.
upper	A vector, indicating the upper bound for the calculation.
n	The number of cells for each dimension, can be one or two numbers.
checkSymmetry	TRUE by default, checks if the covariance matrix is symmetric.
...	Additional arguments, currently not used.

Value

A list, with 'x', 'y' and 'z' components.

---

Numerical computation of the gradient, with parallel capabilities

Description

This function calculates the gradient of a function, numerically, including the possibility of doing it in parallel.

Usage

gradient(fn, x, method, control, parallel, ...)

Arguments

fn	The function to calculate the gradient.
x	The value to compute the gradient at.
method	The method used. Currently implemented: central, backward, forward and Richardson. See details.
control	A list of control arguments.
parallel	Boolean, should numerical derivatives be calculated in parallel?
...	Additional arguments to be passed to fn.

Value

The gradient of fn at x.

Examples

gradient(fn=function(x) sum(x^3), x=0)

---

Objective function between observed and simulated data

Description

Compute an objective value (error, negative log-likelihood, and/or penalty) given observed data ('obs') and model outputs ('sim'). This function is a thin dispatcher: it resolves 'FUN' via [match.fun()] and returns 'FUN(obs = obs, sim = sim, ...)'.

Usage

objFn(obs, sim, FUN, ...)

fitness(obs, sim, FUN, ...)

Arguments

obs	Observed data as expected by 'FUN'. Typically a numeric vector, matrix, or array. Missing values ('NA') are permitted; see *Missing values* below.
sim	Simulated data matching 'obs', in the sense expected by 'FUN'. For most pointwise criteria (e.g. 'norm2', 'lnorm2', 'pois'), 'sim' should have the same shape as 'obs'. For composition criteria (e.g. 'multinom'), 'obs' and 'sim' are expected to be matrices with rows representing samples and columns representing classes.
FUN	Objective function to apply. Can be: a function, e.g. 'FUN = norm2'; or a character string naming a function, e.g. 'FUN = "norm2"'. The function must accept arguments named 'obs' and 'sim' (additional arguments may be supplied via '...') and must return a single numeric scalar.
...	Additional arguments forwarded to 'FUN'.

Details

The returned value is intended to be **minimised** by an optimiser: lower values indicate a better match between 'sim' and 'obs' (or a lower penalty).

Value

A numeric scalar: the value of 'FUN(obs = obs, sim = sim, ...)'. By convention this is an objective to be **minimised**.

Conventions and expectations

• **Minimisation**: all provided objectives are formulated so that lower values indicate a better fit (or weaker penalty).

• **Shapes**: 'objFn()' does not reshape or recycle data. It is the caller's responsibility to supply 'obs' and 'sim' in a compatible form for the chosen 'FUN'.

• **Scalar output**: 'FUN' should return a length-one numeric value. Returning vectors is not supported by 'objFn()' itself (although higher level workflows may aggregate vector-valued objectives elsewhere).

Missing values

Most built-in objectives stop if 'all(is.na(obs))'. Otherwise, they ignore missing values using 'na.rm = TRUE' inside 'sum()'/'mean()'. This means:

• partial 'NA's in 'obs' are ignored;

• 'NA's in 'sim' will also be dropped from sums when 'na.rm = TRUE' (possibly masking simulation failures if not checked upstream).

Numerical constraints and stability

Some objectives impose additional constraints:

• **Poisson** ('pois'): uses 'log(sim)'; 'sim' must be strictly positive wherever 'obs > 0' (otherwise '-Inf' may occur). Consider flooring the simulated intensity, e.g. 'sim <- pmax(sim, 1e-12)'.

• **Log-scale** ('lnorm2', 'lnorm3', 'lnorm4', 'lnorm4b'): apply 'log(x + tiny)'; both 'obs + tiny' and 'sim + tiny' must be positive. If your data may contain zeros, 'tiny' should be chosen accordingly.

• **Compositions** ('multinom'): expects matrix inputs and uses row sums to convert counts/weights into proportions. See details below.

Supported built-in objective functions

The following functions are available with their current behaviour. All are formulated as objectives to be minimised.

norm2

Sum of squared errors on the original scale:

$\sum (obs - sim)^2$

Typical use: continuous observations with approximately additive errors.

lnorm2

Sum of squared errors on the log scale:

$\sum (\log(obs + tiny) - \log(sim + tiny))^2$

Arguments: 'tiny' (default '1e-2'), added before the log to avoid 'log(0)'. Typical use: positive-valued data with multiplicative (lognormal) error structure.

lnorm3

Log-scale squared error with an estimated multiplicative scaling factor $q$. Internally:

• compute element-wise ratios 'ratio <- obs/sim';

• set 'NaN' ratios to 'NA';

• estimate 'q <- mean(ratio, na.rm = TRUE)';

• return $\sum (\log(obs+tiny) - \log(sim+tiny) - \log(q))^2$.

Typical use: positive data where an overall multiplicative bias (scale mismatch) is expected and should not be fully penalised.

lnorm4 / lnorm4b

Extensions of 'lnorm3' that add a penalty term to discourage extreme values of the scaling factor $q$ (or extreme per-observation ratios). They rely on the helper rangeq():

• compute 'ratio <- obs/sim', estimate 'q <- mean(ratio, na.rm=TRUE)';

• compute a penalty using parameters 'b' and 'c':

  $pen = n \cdot (\max(|\log_2(q)|, b)^c - b^c)$

  when 'dump = TRUE' (used by 'lnorm4'), or

  $pen = \sum (\max(|\log_2(ratio)|, b)^c - b^c)$

  when 'dump = FALSE' (used by 'lnorm4b'), where 'n' is the number of non-missing ratios.

• add the penalty to the 'lnorm3'-style objective.

Arguments: 'tiny', 'b' (default '1'), 'c' (default '2'). Typical use: log-scale fitting where scale drift must be controlled.

pois

Poisson negative log-likelihood (up to constants):

$-\sum (obs \log(sim) - sim)$

Typical use: counts (or count-like rates) with Poisson observation error. Note: 'sim' must be positive where 'obs > 0'.

multinom

A composition (multinomial-like) objective operating on matrices. Inputs are expected as 'obs' and 'sim' matrices with:

• rows = samples (e.g. time steps, hauls, sites),

• columns = classes (e.g. age/size bins, categories).

Internal steps (high-level):

• Let $A$ be the number of classes ('A <- ncol(sim)').

• Rows of 'sim' that are all zeros (excluding rows that are all 'NA') are replaced by '1' on that row (interpreted as a uniform prior).

• Row sums are used to compute proportions: $Psim = sim/sum(sim)$, $Pobs = obs/sum(obs)$ (row-wise).

• Rows with 'sum(sim) == 0' are set to 'NA' for numerical convenience.

• Rows with 'sum(obs) == 0' are set to 'NA' (interpreted as “no proportion data available”).

• A variance term 'sigma2' and a small stabiliser 'tiny' are used to define an objective that penalises discrepancies between 'Pobs' and 'Psim'.

Arguments: 'size' (default '20') and 'tiny' (default '1e-3'). Interpretation: 'size' plays the role of an effective sample size (larger values typically increase the weight of the composition fit).

normp / re

Pure penalty on simulated values:

$\sum sim^2$

This ignores 'obs' and can be used as a regulariser, or when 'sim' represents a residual vector or deviates already computed upstream. 're' is an alias of 'normp'.

penalty

Scaled quadratic penalty:

$n \cdot mean(sim^2)$

Arguments: 'n' (default '100'). This assumes a fixed sample size and can be used to put the penalty on a comparable scale across datasets.

Writing your own objective function

You can supply any custom function via 'FUN' provided it:

• accepts arguments named 'obs' and 'sim' (plus optional '...');

• returns a length-one numeric scalar to be minimised;

• defines its own parameter checks and missing-value policy.

See Also

match.fun

Examples

## Basic squared-error objective
obs <- c(1, 2, 3, NA, 5)
sim <- c(1.2, 1.9, 2.7, 4.0, 5.1)
objFn(obs, sim, FUN = "norm2")

## Log-scale objective (positive data)
obs <- c(0.1, 1, 10)
sim <- c(0.2, 0.9, 11)
objFn(obs, sim, FUN = lnorm2, tiny = 1e-2)

## Poisson objective (counts) with flooring for numerical safety
obs <- c(0, 3, 10, 2)
sim <- c(0, 2.5, 9.8, 1.9)
sim <- pmax(sim, 1e-12)
objFn(obs, sim, FUN = "pois")

## Composition objective (matrices: rows = samples, cols = classes)
obs <- rbind(c(10, 5, 0),
            c( 0, 0, 0),  # interpreted as “no composition data”
            c( 2, 1, 7))
sim <- rbind(c( 9, 6, 1),
            c( 0, 0, 0),  # replaced internally by sim+1 on that row
            c( 1, 2, 6))
objFn(obs, sim, FUN = "multinom", size = 20, tiny = 1e-3)

## Custom objective function
my_obj <- function(obs, sim, ...) {
  if (all(is.na(obs))) stop("All observed values are NA.")
  sum(abs(obs - sim), na.rm = TRUE)  # L1 error
}
objFn(obs = c(1, 2, NA), sim = c(1.1, 1.7, 3), FUN = my_obj)

---

Unified optimisation interface with structured parameters and parallel numerical gradients

Description

optim2() provides a unified interface to multiple deterministic and stochastic optimisation methods, combining the algorithms available through stats::optim() with a small set of additional solvers accessed via calibrar's internal dispatcher.

Usage

optim2(
  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
)

Arguments

par	A numeric vector or list. The length of par defines the number of parameters to be estimated (i.e., the dimension of the problem).
fn	The objective function to be minimised. It should accept a parameter vector (or list, depending on the wrapper) as first argument and return either a scalar value (single-objective) or a numeric vector (multi-objective).
gr	A function computing the gradient of fn. If NULL, a numerical approximation is used. Alternatively, a character string can specify the numerical gradient scheme: "central", "forward" (default), "backward", or "richardson".
...	Additional arguments passed to fn and gr.
method	Optimisation method(s) to be used. Can be a single method name or a vector of method names (e.g., one per phase). If NULL, a default is chosen based on replicates (see Details).
lower	Lower bounds for parameters. One value or a vector of the same length as par. NA is treated as -Inf. Default is unconstrained.
upper	Upper bounds for parameters. One value or a vector of the same length as par. NA is treated as Inf. Default is unconstrained.
active	Boolean vector of the same length as par, indicating if the parameter is used in the optimisation (TRUE) or held at a fixed value (FALSE).
control	A list of control options. Common options include ncores, run, master, verbose, REPORT, restart.file, gradient, and gr.method. Additional solver-specific options may be passed through to the underlying optimiser.
hessian	Logical. Should a numerically differentiated Hessian matrix be returned? Currently not implemented.
parallel	Logical. Enable parallel computation (e.g., for replicated evaluations and/or numerical gradients) using up to control$ncores cores.

Details

Methods. The current selection includes (i) base stats::optim() methods, (ii) additional gradient-based solvers from external packages, and (iii) heuristic/global methods.

Comparability. Not all methods are directly comparable: some are deterministic local optimisers (e.g., quasi-Newton), others are derivative-free local searches, and others are stochastic/global heuristics. Method choice should reflect the objective function (smooth vs. non-smooth, deterministic vs. stochastic) and the presence of constraints.

Control arguments. optim2() standardises a set of common control arguments (e.g., iteration limits, tolerances, tracing, and scaling) where supported, while still allowing method-specific control parameters to be passed through to the underlying solver.

Parallel numerical gradients. When analytic gradients are not provided, optim2() can compute finite-difference numerical gradients and distribute function evaluations across multiple cores, which can substantially reduce wall-clock time for expensive objective functions.

Value

A list with components:

par

The best set of parameters found.

value

The value of fn corresponding to par.

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.

convergence

An integer code. 0 indicates successful completion.

message

A character string giving any additional information returned by the optimizer, or NULL.

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.

Choosing an optimisation method

For smooth deterministic objectives, quasi-Newton methods (e.g., "BFGS" or "L-BFGS-B" with box constraints) are often efficient when gradients are available or can be reliably approximated. For box-constrained problems, consider "L-BFGS-B" or the extended bounded solvers (e.g., "Rvmmin", "spg"). For noisy or stochastic objectives, heuristic/global methods (e.g., "AHR-ES") may be more appropriate.

Notes

"SANN" is included for compatibility with stats::optim(), but it is highly sensitive to tuning and often performs poorly on continuous problems under default settings. For noisy or rugged objective functions, "AHR-ES" is generally the recommended heuristic alternative within optim2().

Author(s)

Ricardo Oliveros-Ramos

See Also

optim, nlm, nlminb

Other optimisers: ahres(), calibrate(), optimh()

Examples

optim2(par=rep(NA, 5), fn=sphereN)

---

General-purpose optimization using heuristic algorithms

Description

General-purpose optimization using heuristic algorithms

Usage

optimh(
  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
)

Arguments

par	A numeric vector or list. The length of par defines the number of parameters to be estimated (i.e., the dimension of the problem).
fn	The objective function to be minimised. It should accept a parameter vector (or list, depending on the wrapper) as first argument and return either a scalar value (single-objective) or a numeric vector (multi-objective).
gr	Function to compute the gradient of fn. Ignored by most methods, added for consistency with other optimization functions.
...	Additional arguments passed to fn and gr.
method	Optimisation method(s) to be used. Can be a single method name or a vector of method names (e.g., one per phase). If NULL, a default is chosen based on replicates (see Details).
lower	Lower bounds for parameters. One value or a vector of the same length as par. NA is treated as -Inf. Default is unconstrained.
upper	Upper bounds for parameters. One value or a vector of the same length as par. NA is treated as Inf. Default is unconstrained.
active	Boolean vector of the same length as par, indicating if the parameter is used in the optimisation (TRUE) or held at a fixed value (FALSE).
control	A list of control options. Common options include ncores, run, master, verbose, REPORT, restart.file, gradient, and gr.method. Additional solver-specific options may be passed through to the underlying optimiser.
hessian	Logical. Should a numerically differentiated Hessian matrix be returned? Currently not implemented.
parallel	Logical. Enable parallel computation (e.g., for replicated evaluations and/or numerical gradients) using up to control$ncores cores.

Value

A list with components:

par

The best set of parameters found.

value

The value of fn corresponding to par.

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.

convergence

An integer code. 0 indicates successful completion.

message

A character string giving any additional information returned by the optimizer, or NULL.

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.

Author(s)

Ricardo Oliveros-Ramos

See Also

Other optimisers: ahres(), calibrate(), optim2()

Examples

optim2(par=rep(NA, 5), fn=sphereN)

---

Sphere function with random noise

Description

This function calculates the Euclidian distance from a point to the origin after a random displacement of its position.

Usage

sphereN(x, sd = 0.1, aggregate = TRUE)

Arguments

x	The coordinates of the point
sd	The standard deviation of the noise to be added to the position of x, a normal distribution with mean zero is used.
aggregate	If aggregate is TRUE the distance is returned, otherwise the size of the projection of the distance among each axis.

Value

The distance from the point x to the origin after a random displacement.

Author(s)

Ricardo Oliveros–Ramos

Examples

sphereN(rep(0, 10))

---

Predict time-varying parameters using splines.

Description

Predict time-varying parameters using splines.

Usage

spline_par(par, n, knots = NULL, periodic = FALSE, period = NULL)

Arguments

par	Values at knots
n	Number of points. Time (independent variable) is assumed to be between 0 and n with length(par) equidistant points (including 0 and n).
knots	Position of knots. Default, is length(x) equidistant points between 0 and 1. Always are re-scaled to 0 to 1.
periodic	boolean, is the spline periodic?
period	If periodic is TRUE, it specify the time period.

Value

A list with the interpolates values as 'x' and 'time'.

---