Package 'PESTO'

Description

Uses R's linked LAPACK (which on macOS is Apple Accelerate/AMX, on Linux is typically OpenBLAS or MKL) for hardware-optimised SVD computation.

Usage

accelerate_svd(A, thin = TRUE)

Arguments

Value

A list with components U, d, V.

Examples

set.seed(1L)
A <- matrix(rnorm(8 * 5), nrow = 8, ncol = 5)
res <- accelerate_svd(A, thin = TRUE)
length(res$d)
all.equal(sort(res$d, decreasing = TRUE), svd(A)$d)

Description

Dynamically determines the optimal ensemble size based on convergence diagnostics and information-theoretic criteria.

Usage

adaptive_ensemble_size(
  phi_values,
  current_size,
  min_size = 20L,
  max_size = 500L,
  cv_target = 0.3
)

Arguments

Details

Uses the effective sample size (ESS) and coefficient of variation of phi to determine whether the ensemble is too large (wasting compute) or too small (poor UQ coverage).

Value

A list with recommended_size, reasoning, and diagnostics.

Examples

set.seed(1L)
phi_values <- rnorm(50L, mean = 100, sd = 20)^2
res <- adaptive_ensemble_size(
  phi_values   = phi_values,
  current_size = 50L
)
res$recommended_size
res$cv_phi
res$ess_ratio

Description

Automatically selects the SVD algorithm from the matrix dimensions:

Usage

adaptive_svd(A, k = 0L, method = "auto")

Arguments

Details

• Randomised SVD (Halko et al., 2011): when target rank k << min(m, n)

• Apple Accelerate / LAPACK: a full dense decomposition otherwise

• Eigen BDCSVD: cross-platform divide-and-conquer fallback

Value

A list with components U (m x k), d (k), V (n x k), plus method_used and time_ms.

Examples

set.seed(1L)
A <- matrix(rnorm(20 * 12), nrow = 20, ncol = 12)
res <- adaptive_svd(A, k = 5L, method = "auto")
length(res$d)
is.character(res$method_used)

Description

Builds an in-process forward-model closure for pesto_ies_callback() that drives APSIM (Next Gen .apsimx or Classic .apsim) through the apsimx package without going via the .pst-file path. Used as Year-1 §D4 of the UQ ag-stack roadmap (uq_ag_stack_roadmap_v0.md).

Usage

apsim_callback(
  template,
  param_map,
  output_extractor,
  workdir = tempfile("apsim_cb_"),
  param_writer = NULL,
  simulation_runner = NULL,
  verbose = FALSE
)

Arguments

Details

The returned closure has signature ⁠function(theta) -> obs⁠, where theta is an ⁠nreal x npar⁠ matrix with column names matching names(param_map), and obs is an ⁠nreal x nobs⁠ matrix. Each row is produced by:

1. Copying template into a fresh per-realisation file under workdir.

2. For each ⁠(par_name, node_path)⁠ in param_map, calling the appropriate ⁠apsimx::edit_*⁠ function to write theta[i, par_name] to that node.

3. Calling apsimx::apsimx() (or apsim() for Classic) and passing the returned simulation object to output_extractor(), which must return a length-nobs numeric vector.

Per-realisation failures (APSIM crash, missing output, extractor error) populate an NA row; pesto_ies_callback() then carries that realisation forward unchanged or aborts, depending on its on_failure setting.

Value

A closure of signature ⁠function(theta) -> obs⁠ suitable for the ⁠forward_model =⁠ argument of pesto_ies_callback(). The closure carries an "apsim_version" attribute recording the in-use APSIM binary version (read from Models --version on the configured engine, NA_character_ if it cannot be determined), so a calibrated run can be grounded to the exact simulator that produced it. Thread it into the manifest with as_manifest(fit, apsim_version = attr(fm, "apsim_version")), so a downstream consumer can refuse to compare two manifests built against incompatible APSIM major versions.

APSIM version compatibility

Tested against the apsimx package API as of CRAN 2.7.x. The exact editor function differs between APSIM Next Gen (⁠edit_apsimx*⁠) and APSIM Classic (edit_apsim); selection is by file extension of template. If your installed apsimx version exposes a different editor signature, supply param_writer to override the default per-parameter writer.

Concurrency

The returned closure evaluates whatever block of realisations it is handed and writes each to a unique per-realisation file under workdir, so it is safe to drive in parallel. To run an ensemble concurrently, wrap the closure in a pesto_forward_model() with parallel = "multicore" (or a custom map_fn); the PESTO evaluation engine then dispatches realisations across forked workers, each invoking APSIM on its own input file. Called directly (the bulk path), realisations run serially. apsimx's own thread-safety under heavy ensemble load has not been independently verified, so start with a modest n_cores.

See Also

pesto_ies_callback() for the IES driver; pesto_ies() for the classic .pst-file path.

Examples

## Not run: 
# Requires apsimx and a working APSIM installation
fm <- apsim_callback(
  template  = "wheat_wagga.apsimx",
  param_map = list(
    RUE       = "Wheat.Leaf.Photosynthesis.RUE.FixedValue",
    CN2       = "Soil.SoilWater.CN2Bare"
  ),
  output_extractor = function(sim) {
    # sim is a data.frame; extract end-of-season yield trajectory
    as.numeric(sim$Wheat.Grain.Total.Wt)
  }
)
prior <- matrix(c(runif(40, 1.0, 2.0), runif(40, 60, 90)),
                ncol = 2, dimnames = list(NULL, c("RUE", "CN2")))
fit <- pesto_ies_callback(
  forward_model  = fm,
  prior_ensemble = prior,
  obs            = c(y1 = 4500, y2 = 5200),
  obs_sd         = 200,
  noptmax        = 4
)

## End(Not run)

Description

Generic used by the IES driver so a caller can pass either a bare ⁠function(theta) -> obs⁠ or a fully-specified pesto_forward_model(). A bare function is wrapped with the supplied policy arguments; an existing pesto_forward_model is returned unchanged (its own policy is authoritative and the ... are ignored).

Usage

as_forward_model(x, ...)

Arguments

Value

A pesto_forward_model S7 object.

See Also

pesto_forward_model(), pesto_evaluate().

Examples

fm <- as_forward_model(function(theta) theta, on_failure = "stop")
S7::S7_inherits(fm, pesto_forward_model)

Description

Wraps the plain-list result returned by pesto_ies_callback() (and, eventually, pesto_ies()) in the S7 manifest contract object so downstream packages can consume it via S7 dispatch without reaching into PESTO-specific list internals.

Usage

as_manifest(x, ...)

Arguments

Value

A pesto_ensemble_manifest S7 object.

Description

Issues a warning when the ratio of training points to parameters falls below an empirical threshold, where the Gaussian-process surrogate inside pesto_surrogate_ies() and surrogate_ensemble_update() is unlikely to repay its training cost. The check is a soft guardrail — it does not modify the run, only flags an unfavourable regime so the caller can decide whether to fall back to pure IES.

Usage

check_surrogate_regime(n_params, n_train, threshold = 5)

Arguments

Details

The default threshold of 5 corresponds to the soft floor n_train >= 5 * n_params documented in vignette("surrogate-ies", package = "PESTO"). Below that floor the GP posterior variance typically stays above the uncertainty-driven switching threshold and surrogate savings collapse to near zero.

This is exposed as a stand-alone helper so users can call it explicitly before scheduling an expensive ensemble. It is not invoked automatically by pesto_surrogate_ies() in the current release; that wiring is tracked as a v0.2 enhancement candidate.

Value

Invisibly returns TRUE when the regime is favourable (n_train >= threshold * n_params) and FALSE otherwise. Called for the warning side-effect.

References

Rasmussen, C. E. & Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press.

See Also

pesto_surrogate_ies(), surrogate_ensemble_update(), vignette("surrogate-ies", package = "PESTO")

Examples

# Favourable regime: 100 training points for 10 parameters.
check_surrogate_regime(n_params = 10L, n_train = 100L)

# Unfavourable regime: 30 training points for 30 parameters
# (the curse-of-dimensionality case from Scenario C of the
# comparison-and-simulation vignette). Emits a warning.
suppressWarnings(
  check_surrogate_regime(n_params = 30L, n_train = 30L)
)

# Custom threshold for users with a smoother forward model.
check_surrogate_regime(n_params = 20L, n_train = 60L, threshold = 3)

Description

Calculates the weighted sum of squared residuals for each realisation in the ensemble.

Usage

compute_phi(residuals, weights)

Arguments

Value

Numeric vector (nreal). Phi value per realisation.

Examples

set.seed(1L)
residuals <- matrix(rnorm(5 * 4), 5, 4)
weights   <- rep(1, 5)
phi <- compute_phi(residuals, weights)
length(phi)
phi

Description

Builds an npar x nobs localisation taper directly from the ensemble, with no parameter / observation coordinates required. This is the iterative-ensemble-smoother-native localisation of Luo & Bhakta (2020): spurious sample correlations between a parameter and an observation (an artefact of finite ensemble size) are damped, while genuine correlations that stand above an estimated noise floor are retained.

Usage

correlation_localisation(
  par_diff,
  obs_diff,
  threshold = -1,
  taper = "hard",
  n_shuffle = 1L,
  quantile = 0.95
)

Arguments

Details

The sample correlation $\rho_{ij}$ between parameter-anomaly row $i$ and observation-anomaly row $j$ is compared against a noise floor $\theta$. When threshold < 0 the floor is estimated by destroying the parameter-observation link — the realisation order of the observation anomalies is randomly permuted, independently per replicate, and the floor is taken as a high quantile (default 0.95) of the resulting spurious $|\rho|$ values. The permutation uses R's RNG, so the estimate is reproducible under set.seed().

Two tapers are offered. "hard" keeps correlations above the floor unchanged (weight 1) and zeroes the rest. "soft" applies a smooth, monotone ramp $w_{ij} = \mathrm{clip}((|\rho_{ij}| - \theta) / (1 - \theta), 0, 1)$, which downweights near-floor correlations rather than thresholding them sharply.

Value

A list with rho (the npar x nobs taper), threshold (the floor used), n_active (count of entries with non-zero weight), and frac_active (that count over npar * nobs).

References

Luo, X. & Bhakta, T. (2020). Automatic and adaptive localization for ensemble-based history matching. Journal of Petroleum Science and Engineering, 184, 106559.

Examples

set.seed(1L)
npar <- 8L; nobs <- 5L; nreal <- 40L
pd <- matrix(rnorm(npar * nreal), npar, nreal)
# Make parameter 1 genuinely correlated with observation 1.
od <- matrix(rnorm(nobs * nreal), nobs, nreal)
od[1L, ] <- od[1L, ] + 2 * pd[1L, ]
loc <- correlation_localisation(pd, od)
loc$threshold
loc$rho[1L, 1L]

Description

Builds a pesto_pst object from R data structures, without requiring an existing .pst file.

Usage

create_pest_scenario(
  parameters,
  observations,
  model_command,
  template_files = NULL,
  instruction_files = NULL,
  pestpp_options = list()
)

Arguments

Value

A pesto_pst object.

Examples

pars <- data.table::data.table(
  parnme = c("k1", "k2"),
  partrans = "log",
  parchglim = "factor",
  parval1 = c(1.0, 0.5),
  parlbnd = c(0.01, 0.001),
  parubnd = c(100, 50),
  pargp = "hydraulic"
)
obs <- data.table::data.table(
  obsnme = c("h1", "h2"),
  obsval = c(1.0, 2.0),
  weight = c(1.0, 1.0),
  obgnme = "head"
)
pst <- create_pest_scenario(
  parameters    = pars,
  observations  = obs,
  model_command = "python model.py"
)
inherits(pst, "pesto_pst")
pst$control_data$npar

Description

A ready-to-use ode_forward_model() specialisation for the canonical single-state crop biomass-accumulation curve. Above-ground biomass $B$ follows a logistic growth law

$\frac{dB}{dt} = r\,B\left(1 - \frac{B}{B_{\max}}\right),$

the standard sigmoid description of a crop's dry-matter accumulation over a season: an early near-exponential phase at relative growth rate $r$, decelerating to a canopy- and resource-limited ceiling $B_{\max}$ (Goudriaan & Monteith 1990). The calibration parameters are the relative growth rate r, the asymptotic biomass b_max, and the initial biomass b0.

Usage

crop_growth_forward_model(
  times,
  solver = c("rk4", "desolve"),
  n_steps = 10L,
  ...
)

Arguments

Details

The forward map integrates the logistic ODE across times and returns the modelled biomass at every time after the first – exactly the shape a destructive-harvest or remote-sensing biomass series takes, so the returned object calibrates directly against an observed biomass-over-time vector through pesto_ies_callback(). The single-state logistic form is deliberately the simplest defensible crop-growth template; richer multi-organ partitioning models compose through the same ode_forward_model() entry point by supplying a vector-valued derivs.

Value

A pesto_forward_model() S7 object with param_names c("r", "b_max", "b0") and n_obs = length(times) - 1L.

References

Goudriaan, J. & Monteith, J. L. (1990). A mathematical function for crop growth based on light interception and leaf area expansion. Annals of Botany, 66(6), 695–701.

See Also

ode_forward_model() for the generic builder; seir_forward_model() for the compartmental epidemic template; pesto_ies_callback() for calibration.

Examples

# Simulate a biomass series at a known parameter, then recover it.
times <- seq(0, 120, by = 15)
fm <- crop_growth_forward_model(times = times)
truth <- matrix(c(0.06, 1400, 20), nrow = 1L,
                dimnames = list(NULL, c("r", "b_max", "b0")))
biomass <- as.numeric(pesto_evaluate(fm, truth))
round(biomass)

Description

Implements the core IES ensemble update equation using the Gauss-Levenberg-Marquardt (GLM) formulation from Chen & Oliver (2013). This is the computational hotspot of the iterative ensemble smoother.

Usage

ensemble_solution(
  par_diff,
  obs_diff,
  obs_resid,
  par_resid,
  weights,
  parcov_inv,
  Am,
  cur_lam,
  eigthresh = 1e-06,
  use_approx = TRUE,
  use_prior_scaling = FALSE,
  iter = 1L,
  reg_factor = -1
)

Arguments

Details

The update equation solves:

$\Delta\theta = -\Delta\theta' V s (s^2 + \lambda I)^{-1} U^T r$

where the SVD is performed on the scaled observation difference matrix.

Value

Matrix (nreal x npar). Parameter upgrade vectors (one row per realisation). The returned matrix is the negative-direction step from the Chen-Oliver 2013 GLM update formula. To advance the ensemble, apply by subtraction: par_new = par_old - upgrade. The R-side driver pesto_ies_callback() handles this convention internally.

Examples

set.seed(1L)
npar  <- 4L
nreal <- 20L
nobs  <- 30L
par_diff  <- matrix(rnorm(npar * nreal), npar, nreal)
obs_diff  <- matrix(rnorm(nobs * nreal), nobs, nreal)
obs_resid <- matrix(rnorm(nobs * nreal, sd = 0.5), nobs, nreal)
par_resid <- matrix(rnorm(npar * nreal, sd = 0.1), npar, nreal)
weights    <- rep(1, nobs)
parcov_inv <- rep(1, npar)
Am         <- matrix(0, 0, 0)
upgrade <- ensemble_solution(
  par_diff   = par_diff,
  obs_diff   = obs_diff,
  obs_resid  = obs_resid,
  par_resid  = par_resid,
  weights    = weights,
  parcov_inv = parcov_inv,
  Am         = Am,
  cur_lam    = 1.0
)
dim(upgrade)

Description

A variant of ensemble_solution() that selects the SVD backend automatically – a randomised SVD for low-rank problems, otherwise a dense LAPACK / Accelerate decomposition – and returns the upgrade together with timing diagnostics. It is the convenient entry point when the best backend for a given problem size is not known in advance. All computation is on the CPU.

Usage

ensemble_solution_adaptive(
  par_diff,
  obs_diff,
  obs_resid,
  par_resid,
  weights,
  parcov_inv,
  Am,
  cur_lam,
  eigthresh = 1e-06,
  use_approx = TRUE,
  use_prior_scaling = FALSE,
  iter = 1L,
  reg_factor = -1,
  svd_method = "auto",
  target_rank = 0L
)

Arguments

Value

A list with upgrade matrix and performance diagnostics.

Examples

set.seed(1L)
npar  <- 4L
nreal <- 20L
nobs  <- 30L
par_diff  <- matrix(rnorm(npar * nreal), npar, nreal)
obs_diff  <- matrix(rnorm(nobs * nreal), nobs, nreal)
obs_resid <- matrix(rnorm(nobs * nreal, sd = 0.5), nobs, nreal)
par_resid <- matrix(rnorm(npar * nreal, sd = 0.1), npar, nreal)
weights    <- rep(1, nobs)
parcov_inv <- rep(1, npar)
Am         <- matrix(0, 0, 0)
res <- ensemble_solution_adaptive(
  par_diff   = par_diff,
  obs_diff   = obs_diff,
  obs_resid  = obs_resid,
  par_resid  = par_resid,
  weights    = weights,
  parcov_inv = parcov_inv,
  Am         = Am,
  cur_lam    = 1.0,
  svd_method = "auto"
)
dim(res$upgrade)

Description

ensemble_solution_gpu() was renamed to ensemble_solution_adaptive(). The original name implied GPU computation; the function performs CPU adaptive-SVD backend selection (randomised SVD versus a dense LAPACK / Accelerate decomposition). The alias is retained for backward compatibility and will be removed in a future release.

Usage

ensemble_solution_gpu(...)

Arguments

Value

The list returned by ensemble_solution_adaptive().

See Also

ensemble_solution_adaptive()

Examples

set.seed(1L)
np <- 4L; no <- 30L; nr <- 20L
pd <- matrix(rnorm(np * nr), np, nr)
od <- matrix(rnorm(no * nr), no, nr)
or_ <- matrix(rnorm(no * nr, sd = 0.5), no, nr)
Am <- matrix(0, 0, 0)
suppressWarnings(
  res <- ensemble_solution_gpu(pd, od, or_, pd, rep(1, no), rep(1, np),
                               Am, cur_lam = 1.0)
)
dim(res$upgrade)

Description

Computes the IES Gauss-Levenberg-Marquardt update with state-space covariance localisation applied as a Schur (elementwise) product on the explicit Kalman gain. The standard SVD kernel ensemble_solution() works in the reduced observation-anomaly subspace and never forms the npar x nobs gain, so it cannot host localisation; this kernel reconstructs the gain

$K = \Delta\Theta\, V\, \mathrm{diag}(s)\, \mathrm{diag}((s^2 + (\lambda+1))^{-1})\, U^{T},$

(with $U s V^{T}$ the thin SVD of the weight-scaled observation anomalies) tapers it as $K \circ \rho$, and applies it to the weighted residuals.

Usage

ensemble_solution_localised(
  par_diff,
  obs_diff,
  obs_resid,
  weights,
  rho,
  cur_lam,
  eigthresh = 1e-06
)

Arguments

Details

When $\rho \equiv 1$ the result is identical (to truncation tolerance) to ensemble_solution() with use_approx = TRUE; the prior-scaling null-space correction (upgrade_2) is not part of the localised path. The returned matrix follows the same sign convention as ensemble_solution() — it is the negative-direction step, applied to the ensemble by subtraction (par_new = par_old - upgrade).

Value

Matrix (nreal x npar). Negative-direction parameter upgrade, applied by subtraction.

References

Chen, Y. & Oliver, D.S. (2013). Levenberg-Marquardt forms of the iterative ensemble smoother for efficient history matching and uncertainty quantification. Computational Geosciences, 17(4), 689–703.

Examples

set.seed(1L)
npar <- 4L; nreal <- 20L; nobs <- 6L
par_diff  <- matrix(rnorm(npar * nreal), npar, nreal)
obs_diff  <- matrix(rnorm(nobs * nreal), nobs, nreal)
obs_resid <- matrix(rnorm(nobs * nreal, sd = 0.5), nobs, nreal)
weights   <- rep(1, nobs)
rho       <- matrix(1, npar, nobs)          # no localisation
upg <- ensemble_solution_localised(
  par_diff, obs_diff, obs_resid, weights, rho, cur_lam = 1.0
)
dim(upg)

Description

Implements the Multiple Data Assimilation (MDA) update from Evensen (2018), Section 14.3.2. Uses low-rank representation of the error covariance.

Usage

ensemble_solution_mda(
  par_diff,
  obs_diff,
  obs_resid,
  obs_err,
  cur_lam = 1,
  eigthresh = 1e-06
)

Arguments

Value

Matrix (nreal x npar). Parameter upgrade vectors.

Examples

set.seed(1L)
npar  <- 4L
nreal <- 20L
nobs  <- 30L
par_diff  <- matrix(rnorm(npar * nreal), npar, nreal)
obs_diff  <- matrix(rnorm(nobs * nreal), nobs, nreal)
obs_resid <- matrix(rnorm(nobs * nreal, sd = 0.5), nobs, nreal)
obs_err   <- matrix(rnorm(nobs * nreal, sd = 0.5), nobs, nreal)
upgrade <- ensemble_solution_mda(
  par_diff  = par_diff,
  obs_diff  = obs_diff,
  obs_resid = obs_resid,
  obs_err   = obs_err,
  cur_lam   = 1.0
)
dim(upgrade)

Description

Diagnoses ensemble collapse (under-dispersion) by the participation ratio of the parameter-anomaly covariance eigenspectrum. Given the parameter anomalies $\Delta\Theta$ (deviations from the ensemble mean, npar x nreal), the anomaly covariance is $C = \Delta\Theta \Delta\Theta^{T} / (N - 1)$ with eigenvalues $\lambda_i = s_i^2 / (N - 1)$, where $s_i$ are the singular values of $\Delta\Theta$.

Usage

ensemble_spread_ess(par_diff)

Arguments

Details

The spectral spread-ESS is the participation ratio

$\mathrm{ESS} = \frac{(\sum_i \lambda_i)^2}{\sum_i \lambda_i^2} = \frac{(\sum_i s_i^2)^2}{\sum_i s_i^4},$

the effective number of directions carrying variance. It is bounded in $[1, r_{\max}]$ with $r_{\max} = \min(\mathrm{npar}, N - 1)$: equal to $r_{\max}$ when variance is spread isotropically across all modes, and approaching 1 when the ensemble collapses onto a single direction. Because the ratio is invariant to a global rescaling of the anomalies, it isolates the shape of the collapse (directional degeneracy) from its magnitude; magnitude is tracked separately by the R-side spread-retention ratio.

Value

A list with components ess (the spectral spread-ESS), r_max (the maximum attainable value $\min(\mathrm{npar}, N - 1)$), and ess_ratio (ess / r_max, in $(0, 1]$).

References

Bretherton, C.S., Widmann, M., Dymnikov, V.P., Wallace, J.M. & Blade, I. (1999). The effective number of spatial degrees of freedom of a time-varying field. Journal of Climate, 12(7), 1990–2009.

Examples

set.seed(1L)
# Healthy isotropic spread -> ESS near r_max
good <- matrix(rnorm(6L * 40L), 6L, 40L)
ensemble_spread_ess(good)$ess_ratio
# Collapsed onto one direction -> ESS near 1
v <- rnorm(6L)
bad <- outer(v, rnorm(40L)) + matrix(rnorm(6L * 40L, sd = 1e-3), 6L, 40L)
ensemble_spread_ess(bad)$ess_ratio

Description

Evaluates the Gaspari & Cohn (1999) fifth-order piecewise-rational compactly-supported correlation function on a matrix of distances. This is the classical distance-based localisation taper: a smooth bump that equals 1 at zero distance, decays to 0 at twice the localisation radius, and is identically 0 beyond. Used to taper the Kalman gain when the parameters and observations carry a spatial (or otherwise metric) coordinate.

Usage

gaspari_cohn(distances, radius)

Arguments

Details

With $z = d / c$ (distance over localisation radius $c$):

$G(z) = \begin{cases} -\tfrac{1}{4}z^5 + \tfrac{1}{2}z^4 + \tfrac{5}{8}z^3 - \tfrac{5}{3}z^2 + 1 & 0 \le z \le 1 \\ \tfrac{1}{12}z^5 - \tfrac{1}{2}z^4 + \tfrac{5}{8}z^3 + \tfrac{5}{3}z^2 - 5z + 4 - \tfrac{2}{3}z^{-1} & 1 < z \le 2 \\ 0 & z > 2. \end{cases}$

Value

Matrix (npar x nobs) of taper weights in $[0, 1]$.

References

Gaspari, G. & Cohn, S.E. (1999). Construction of correlation functions in two and three dimensions. Quarterly Journal of the Royal Meteorological Society, 125(554), 723–757.

Examples

d <- matrix(c(0, 0.5, 1, 1.5, 2, 3), nrow = 2L)
gaspari_cohn(d, radius = 1.0)

Description

Debiases a cheap-fidelity output ensemble against a sparse expensive sample, per observation dimension. For each observation j it fits the first-order autoregressive control variate (the linear term of the Kennedy-O'Hagan multi-fidelity model) high_j ~ a_j + b_j * low_j on the paired subset, with b_j = cov(high_j, low_j) / var(low_j) the variance-minimising coefficient, then predicts the corrected high-fidelity output for every realisation as a_j + b_j * low_all_j.

Usage

mf_control_variate(low_all, high_sub, low_sub)

Arguments

Details

The estimator degrades gracefully: where the cheap output has zero variance on the subset it falls back to the expensive subset mean (b_j = 0), and where the two fidelities are weakly correlated the correction shrinks toward that mean rather than amplifying noise.

This is the plug-in primitive for surrogate cascades: a cascade runs the cheap level over the full ensemble, the expensive level over a chosen subset, and calls this to lift the cheap ensemble toward the expensive one at a fraction of the cost.

Value

A ⁠nreal x nobs⁠ matrix of bias-corrected outputs, with attributes "intercept" (a_j), "slope" (b_j), and "subset_cor" (per-dimension subset correlation; NA for a degenerate dimension).

References

Kennedy, M. C. & O'Hagan, A. (2000). Predicting the output from a complex computer code when fast approximations are available. Biometrika, 87(1), 1–13.

See Also

pesto_multifidelity_model().

Examples

set.seed(1L)
low_all  <- matrix(rnorm(40L), ncol = 2L)
sub      <- 1:5
low_sub  <- low_all[sub, , drop = FALSE]
high_sub <- 0.3 + 1.2 * low_sub + matrix(rnorm(10L, sd = 0.01), ncol = 2L)
corrected <- mf_control_variate(low_all, high_sub, low_sub)
attr(corrected, "slope")

Description

Each calibration parameter is a named entry of the theta matrix supplied by the IES driver. param_names declares which columns the model consumes and in what order. For each realisation the template integrates

$\frac{d\mathbf{y}}{dt} = f(t, \mathbf{y}, \theta)$

from times[1] with initial state y0 – itself possibly a function of theta – and applies observe() to the state trajectory to produce a length-n_obs numeric vector.

Usage

ode_forward_model(
  derivs,
  y0,
  times,
  param_names = character(0),
  observe = NULL,
  solver = c("rk4", "desolve"),
  n_steps = 10L,
  desolve_method = "lsoda",
  n_obs = NA_integer_,
  ...
)

Arguments

Details

Wraps a user-supplied ODE right-hand side in a typed pesto_forward_model() whose forward map integrates the system over a fixed time grid and reads the observation vector off the resulting state trajectory. This is the generic differential_equations template that crop_growth_forward_model() and seir_forward_model() specialise; use it directly for any compartmental or mechanistic ODE.

The derivative function derivs has signature ⁠function(t, y, theta) -> dydt⁠, returning a numeric vector the same length as the state y. The initial state y0 may be a plain numeric vector (shared across realisations) or a function ⁠function(theta) -> y0⁠ when the starting condition is itself calibrated (for example an unknown initial inoculum). The observation map observe has signature ⁠function(traj, theta) -> obs⁠, where traj is an ⁠length(times) x length(y0)⁠ matrix of states at each integration time (column order matching y0); the default reads the first state variable at every time after the first.

Two integrators are available. The default solver = "rk4" is a self-contained classical fourth-order Runge-Kutta with n_steps fixed sub-steps between successive observation times – no external dependency, deterministic, and adequate for the smooth non-stiff systems these templates target. solver = "desolve" delegates to deSolve::ode() (an optional Suggests dependency) with method desolve_method, which brings adaptive step control and stiff solvers ("lsoda", the default, switches between non-stiff and stiff automatically). A realisation whose integration fails or returns a non-finite trajectory is reported as an NA row, which pesto_evaluate() and pesto_ies_callback() handle under their failure policy.

Value

A pesto_forward_model() S7 object whose forward map integrates the ODE system. Pass it to pesto_ies_callback() as forward_model, evaluate it directly with pesto_evaluate(), or bundle several across fidelity levels with pesto_multifidelity_model().

References

Soetaert, K., Petzoldt, T. & Setzer, R. W. (2010). Solving differential equations in R: package deSolve. Journal of Statistical Software, 33(9), 1–25. doi:10.18637/jss.v033.i09

See Also

crop_growth_forward_model() and seir_forward_model() for ready-made specialisations; pesto_forward_model() for the contract; pesto_ies_callback() for the IES driver.

Examples

# Exponential decay dy/dt = -k y, observed at five times. Calibrate k.
fm <- ode_forward_model(
  derivs      = function(t, y, theta) -theta[["k"]] * y,
  y0          = c(y = 1),
  times       = seq(0, 4, by = 1),
  param_names = "k"
)
theta <- matrix(c(0.5, 1.0), ncol = 1L, dimnames = list(NULL, "k"))
pesto_evaluate(fm, theta)

Description

A versioned, hashed, provenance-tracked container for the result of a PESTO ensemble run. It is a portable data contract: a documented, versioned format that any downstream tool can read without depending on PESTO internals, carrying enough provenance to make a run independently reproducible and independently verifiable (via verify_manifest()).

Arguments

Details

Construct directly via the class constructor, or — preferred — via as_manifest() applied to a pesto_ies_callback_result.

Persistence

See write_manifest() / read_manifest() for YAML serialisation and verify_manifest() for integrity checking. The YAML file is the authoritative human-readable entry; params, outputs, and the assimilation context (weights + obs_target) are emitted as sidecar files whose relative paths are recorded inside the YAML. Three sidecar modes:

• format = "rds" (default, verifiable) — RDS only. IEEE 754 doubles round-trip bit-exactly; SHA-256 integrity check holds.

• format = "both" (verifiable) — RDS sidecars plus parallel ⁠*_inspection.csv⁠ files. Hash bound to RDS; CSVs are decorative.

• format = "csv_unverified" (not verifiable, renamed from "csv" in PESTO 0.3.2) — CSV-only sidecars. Hash is recorded but cannot be recomputed from disk (CSV formatter precision loss ~1 ULP). verify_manifest() returns ok = NA + message. Use only for one-way export to non-R analysts; round-trip integrity is by construction broken.

The YAML carries an explicit ⁠integrity:⁠ field ("verifiable" / "not_verifiable") derived from format so downstream non-R tools can branch on the integrity contract without needing to know the PESTO-specific format vocabulary.

Examples

set.seed(1)
npar <- 2L; nobs <- 3L; nreal <- 30L
G <- matrix(stats::rnorm(nobs * npar), nobs, npar)
y <- as.numeric(G %*% c(0.5, -1.0)) + stats::rnorm(nobs, sd = 0.05)
names(y) <- paste0("o", seq_len(nobs))
prior <- matrix(stats::rnorm(nreal * npar), nreal, npar,
                dimnames = list(NULL, c("a", "b")))
fit <- pesto_ies_callback(function(t) t %*% t(G), prior, y, 0.05,
                          noptmax = 3, verbose = FALSE)
m <- as_manifest(fit, seed = 1L)
verify_manifest(m)$ok

Description

Runs the forward model on a parameter matrix under its own failure policy and concurrency strategy, returning a shape-guaranteed ⁠nreal x nobs⁠ observation matrix. Failed realisations populate NA rows (when on_failure = "na"). The returned matrix carries two attributes: "n_failures" (integer count of NA rows) and "fail_idx" (integer realisation indices that failed).

Usage

pesto_evaluate(model, theta, ...)

Arguments

Value

An ⁠nreal x nobs⁠ numeric matrix with attributes "n_failures" and "fail_idx".

See Also

pesto_forward_model(), pesto_multifidelity_model().

Examples

fm <- pesto_forward_model(fn = function(theta) theta[, 1, drop = FALSE],
                          n_obs = 1L)
pesto_evaluate(fm, matrix(c(1, 2, 3), ncol = 1L))

Description

A pesto_forward_model wraps a user callable of signature ⁠function(theta) -> obs⁠ – where theta is an ⁠nreal x npar⁠ numeric matrix and obs is an ⁠nreal x nobs⁠ numeric matrix – in a typed object that owns the evaluation contract: output dimensionality, expected parameter names, the failure policy, the concurrency strategy, and a fidelity tag.

Usage

pesto_forward_model(
  fn,
  n_obs = NA_integer_,
  param_names = character(0),
  on_failure = "na",
  parallel = "serial",
  n_cores = NA_integer_,
  map_fn = NULL,
  max_fail_frac = 1,
  fidelity = 0L,
  label = NA_character_
)

Arguments

Details

This is the single contract both PESTO adapter modes honour. The native callback driver pesto_ies_callback() accepts either a bare function (auto-wrapped via as_forward_model() with that driver's failure policy) or a pesto_forward_model built here; the apsimx adapter apsim_callback() can emit one directly. Evaluation is via pesto_evaluate(), which guarantees the returned shape and accounts for failed realisations as NA rows.

Value

A pesto_forward_model S7 object.

Concurrency

With parallel = "serial" (the default) and no map_fn, evaluation attempts a single bulk call fn(theta) and falls back to a serial per-realisation loop only if the bulk call errors – preserving the fast path for vectorised forward models. With parallel = "multicore" realisations are dispatched per row through parallel::mclapply() (fork-based; silently serial on Windows). A custom map_fn (an lapply-shaped ⁠function(X, FUN, ...)⁠) overrides both and lets callers plug in future.apply::future_lapply, mirai, or a cluster backend. For reproducible parallel runs set RNGkind("L'Ecuyer-CMRG") and set.seed() before evaluating; parallel::mclapply() then draws independent streams per realisation.

Fidelity

fidelity is an integer level tag (0L = base / cheapest by convention; higher = more expensive / higher resolution). A single-fidelity model leaves it at 0L. The tag is what pesto_multifidelity_model() uses to order levels and what the manifest emitter threads into provenance.

See Also

pesto_evaluate() to evaluate one; as_forward_model() to coerce a bare function; pesto_multifidelity_model() to compose several across fidelity levels; pesto_ies_callback() for the IES driver that consumes it.

Examples

# A vectorised linear forward model wrapped as a contract object.
G  <- matrix(c(1, 0, 0, 1, 1, -1), nrow = 3L, byrow = TRUE)
fm <- pesto_forward_model(
  fn          = function(theta) theta %*% t(G),
  n_obs       = 3L,
  param_names = c("a", "b")
)
theta <- matrix(c(1, 2, 3, 4), nrow = 2L, byrow = TRUE,
                dimnames = list(NULL, c("a", "b")))
pesto_evaluate(fm, theta)

Description

Executes the pestpp-glm algorithm for deterministic parameter estimation.

Usage

pesto_glm(
  pst_file,
  exe = NULL,
  noptmax = NULL,
  extra_args = list(),
  working_dir = NULL,
  verbose = TRUE
)

Arguments

Value

A list of class pesto_glm_result.

Examples

if (nzchar(Sys.which("pestpp-glm"))) {
  pars <- data.table::data.table(
    parnme = c("k1", "k2"), partrans = "log", parchglim = "factor",
    parval1 = c(1.0, 0.5), parlbnd = c(0.01, 0.001),
    parubnd = c(100, 50), pargp = "hydraulic"
  )
  obs <- data.table::data.table(
    obsnme = c("h1", "h2"), obsval = c(1.0, 2.0),
    weight = c(1.0, 1.0), obgnme = "head"
  )
  pst <- create_pest_scenario(pars, obs, model_command = "echo run")
  tf <- tempfile(fileext = ".pst")
  on.exit(unlink(tf), add = TRUE)
  write_pst(pst, tf)
  res <- pesto_glm(tf, noptmax = 1, verbose = FALSE)
  res$exit_code
}

Description

Executes the pestpp-ies algorithm on a PEST control file. This is the primary method for ensemble-based parameter estimation and uncertainty quantification.

Usage

pesto_ies(
  pst_file,
  exe = NULL,
  num_reals = 50,
  noptmax = NULL,
  lambda_scale_fac = c(0.1, 0.5, 1),
  ies_par_en = NULL,
  extra_args = list(),
  working_dir = NULL,
  verbose = TRUE
)

Arguments

Value

A list of class pesto_ies_result containing:

phi

data.table of objective function values per iteration

par_ensemble

Final parameter ensemble (data.table)

obs_ensemble

Final observation ensemble (data.table)

exit_code

Integer exit code from pestpp-ies

runtime_seconds

Total wall-clock runtime

Examples

if (nzchar(Sys.which("pestpp-ies"))) {
  pars <- data.table::data.table(
    parnme = c("k1", "k2"), partrans = "log", parchglim = "factor",
    parval1 = c(1.0, 0.5), parlbnd = c(0.01, 0.001),
    parubnd = c(100, 50), pargp = "hydraulic"
  )
  obs <- data.table::data.table(
    obsnme = c("h1", "h2"), obsval = c(1.0, 2.0),
    weight = c(1.0, 1.0), obgnme = "head"
  )
  pst <- create_pest_scenario(pars, obs, model_command = "echo run")
  tf <- tempfile(fileext = ".pst")
  on.exit(unlink(tf), add = TRUE)
  write_pst(pst, tf)
  res <- pesto_ies(tf, num_reals = 3, noptmax = 1, verbose = FALSE)
  res$exit_code
}

Description

Drives an Iterative Ensemble Smoother entirely in R, using a user-supplied forward model callable instead of the PEST++ .pst-file invocation cycle. Each iteration:

Usage

pesto_ies_callback(
  forward_model,
  prior_ensemble,
  obs,
  obs_sd,
  noptmax = 4L,
  lambda = 1,
  fidelity_schedule = NULL,
  parcov = NULL,
  eigthresh = 1e-06,
  use_approx = TRUE,
  inflation = NULL,
  localisation = NULL,
  on_failure = c("na", "stop"),
  verbose = TRUE,
  phi_tol = NULL
)

Arguments

Details

1. Evaluates forward_model(par_ensemble) to obtain simulated observations.

2. Forms parameter / observation anomalies and residuals.

3. Calls the C++ kernel ensemble_solution() (GLM form, Chen & Oliver 2013) to compute an ⁠nreal x npar⁠ upgrade.

4. Adds the upgrade to the current ensemble.

The classic .pst-file path remains available via pesto_ies() for full PEST++ compatibility. Use this callback driver when the forward model is itself an R function (e.g. an apsimx wrapper from apsim_callback(), a Python bridge, or a synthetic test problem) and the per-realisation file-I/O overhead of the .pst path is the bottleneck.

Phase-1 behaviour: single lambda per iteration (or a user-supplied schedule). A line-search over ies_lambda_mults matching pestpp-ies is a planned Phase-2 enhancement; for the common case of a well-behaved forward model with lambda = 1, the GLM update reduces phi reliably (see vignette apsim-callback).

Value

A list of class c("pesto_ies_callback_result", "pesto_ies_result") with components:

phi

data.table of per-realisation phi by iteration.

par_ensemble

Final parameter ensemble (data.table).

obs_ensemble

Final simulated-observation ensemble (data.table).

iterations

List of per-iteration metadata: lambda, mean_phi, n_failures, and the dispersion diagnostics spread_ess / spread_ess_ratio (ensemble_spread_ess()), plus inflation_method / inflation_factor / retention and localisation / loc_threshold / loc_frac_active when those countermeasures are active.

runtime_seconds

Total wall-clock runtime.

n_forward_evals

Total number of realisation-level forward evaluations across all iterations (including the final refresh).

failure_rate

Fraction of forward evaluations that returned NA.

converged

Logical: TRUE if the run stopped early on the phi_tol convergence criterion; FALSE when phi_tol is NULL or the full noptmax was reached.

n_iterations

Number of IES iterations actually run (fewer than noptmax if the convergence checker stopped the run early).

fidelity

For a pesto_multifidelity_model() run, a provenance list list(type, schedule, final_level, n_levels, costs) recording the realised per-iteration fidelity schedule; NULL for a single-fidelity run. Consumed by as_manifest() to populate the manifest contract.

Multi-fidelity

When forward_model is a pesto_multifidelity_model(), fidelity_schedule selects which fidelity level each iteration evaluates – a recycled / padded integer vector, one entry per iteration (default: the highest fidelity every iteration, i.e. exactly the single-fidelity behaviour). This supports fidelity ramping (cheap early iterations, expensive late ones); the final ensemble refresh always uses the highest fidelity so the returned posterior is at full resolution. The control-variate combiner mf_control_variate() is the plug-in point for bias-corrected surrogate cascades.

References

Chen, Y. & Oliver, D.S. (2013). Levenberg-Marquardt forms of the iterative ensemble smoother for efficient history matching and uncertainty quantification. Computational Geosciences, 17(4), 689–703.

White, J.T. (2018). A model-independent iterative ensemble smoother for efficient history-matching and uncertainty quantification in very high dimensions. Environmental Modelling & Software, 109, 191–201.

See Also

pesto_ies() for the .pst-file path; apsim_callback() for the apsimx adapter.

Examples

# Linear-Gaussian recovery toy
set.seed(1)
npar <- 3; nobs <- 6; nreal <- 80
G <- matrix(rnorm(nobs * npar), nobs, npar)
theta_true <- c(1.0, -0.5, 2.0)
y <- as.numeric(G %*% theta_true) + rnorm(nobs, sd = 0.05)
f <- function(theta) theta %*% t(G)
prior <- matrix(rnorm(nreal * npar), nreal, npar,
                dimnames = list(NULL, paste0("p", 1:npar)))
fit <- pesto_ies_callback(
  forward_model = f, prior_ensemble = prior,
  obs = setNames(y, paste0("o", 1:nobs)), obs_sd = 0.05,
  noptmax = 5, verbose = FALSE
)
colMeans(as.matrix(fit$par_ensemble[, -1]))  # should approach theta_true

Description

Assimilates observations in time-ordered windows against a static parameter ensemble. Each window runs an IES Gauss-Levenberg-Marquardt update using only that window's observation block; the updated ensemble is carried forward as the prior for the next window. This is the filtering analogue of pesto_ies_callback(): instead of one batch assimilation of all observations, the posterior is refined window by window and is available after each.

Usage

pesto_ies_filter(
  forward_model,
  prior_ensemble,
  obs,
  obs_sd,
  windows,
  window_noptmax = 1L,
  lambda = 1,
  fidelity_schedule = NULL,
  parcov = NULL,
  eigthresh = 1e-06,
  use_approx = TRUE,
  inflation = NULL,
  localisation = NULL,
  on_failure = c("na", "stop"),
  verbose = TRUE
)

Arguments

Value

A list of class c("pesto_ies_filter_result", "pesto_ies_result") with components:

phi

data.table of per-realisation phi by window (on that window's observation block).

par_ensemble

Final parameter ensemble (data.table).

obs_ensemble

Final simulated-observation ensemble, full nobs columns (data.table).

windows

List of per-window metadata: assimilated indices, lambda, mean phi, per-parameter ensemble mean and standard deviation (the sd trace shows the posterior tightening), and failure count.

runtime_seconds, n_forward_evals, failure_rate

Run totals.

fidelity

Multi-fidelity provenance (or NULL), as in pesto_ies_callback().

Filter vs smoother

pesto_ies_callback() forms parameter residuals against the original prior mean and assimilates every observation together. Here each window forms residuals against the current (carried-forward) ensemble mean and assimilates only its own block, so information accrues sequentially. With the default use_approx = TRUE the carried-forward ensemble itself encodes the background covariance, so the sequential behaviour comes from propagating the ensemble, not from an explicit prior term.

Multi-fidelity

When forward_model is a pesto_multifidelity_model(), fidelity_schedule selects the fidelity level evaluated at each window (recycled / padded to the number of windows; default: highest fidelity throughout). The final ensemble refresh always uses the highest fidelity.

References

Chen, Y. & Oliver, D.S. (2013). Levenberg-Marquardt forms of the iterative ensemble smoother for efficient history matching and uncertainty quantification. Computational Geosciences, 17(4), 689–703.

See Also

pesto_ies_callback() for the batch smoother; pesto_multifidelity_model() for fidelity stacks; as_manifest() to wrap the result in the ensemble-manifest contract.

Examples

# Linear-Gaussian recovery, assimilated in three observation windows.
set.seed(1)
npar <- 3; nobs <- 9; nreal <- 80
G <- matrix(rnorm(nobs * npar), nobs, npar)
theta_true <- c(1.0, -0.5, 2.0)
y <- as.numeric(G %*% theta_true) + rnorm(nobs, sd = 0.05)
f <- function(theta) theta %*% t(G)
prior <- matrix(rnorm(nreal * npar), nreal, npar,
                dimnames = list(NULL, paste0("p", 1:npar)))
fit <- pesto_ies_filter(
  forward_model  = f, prior_ensemble = prior,
  obs = setNames(y, paste0("o", 1:nobs)), obs_sd = 0.05,
  windows = list(1:3, 4:6, 7:9), verbose = FALSE
)
# Posterior sd should shrink window over window:
vapply(fit$windows, function(w) mean(w$par_sd), numeric(1))

Description

print() gives a one-screen summary of a fitted iterative ensemble-smoother run (driver, problem dimensions, the phi-convergence trace, the spread-ESS dispersion diagnostic and the failure rate). plot() draws the objective-function (phi) convergence trace by delegating to plot_phi(). For posterior parameter distributions (prior vs posterior) use plot_ensemble().

Usage

## S3 method for class 'pesto_ies_result'
print(x, ...)

## S3 method for class 'pesto_ies_result'
plot(x, ...)

Arguments

Value

print() returns x invisibly; plot() returns a ggplot2 object.

Examples

set.seed(1)
G <- matrix(rnorm(18L), 6L, 3L)
prior <- matrix(rnorm(150L), 50L, 3L, dimnames = list(NULL, paste0("p", 1:3)))
obs <- stats::setNames(as.numeric(G %*% c(1, -0.5, 2)) + rnorm(6L, sd = 0.05),
                       paste0("o", 1:6))
fit <- pesto_ies_callback(function(t) t %*% t(G), prior, obs,
                          obs_sd = 0.05, noptmax = 4L, verbose = FALSE)
print(fit)
p <- plot(fit)        # phi convergence (a ggplot2 object)

Description

Builds a control object that tells pesto_ies_callback() and pesto_ies_filter() how to counteract ensemble under-dispersion — the progressive collapse of posterior spread that an iterative ensemble smoother suffers with a finite ensemble. Inflation re-expands the post-update parameter spread each iteration; the default method = "none" leaves the update byte-identical to the un-inflated smoother.

Usage

pesto_inflation(
  method = c("none", "rtps", "adaptive", "multiplicative"),
  alpha = 0.5,
  factor = 1,
  retention_floor = 0.5,
  max_factor = 5
)

Arguments

Details

Four methods are offered. "rtps" is relaxation-to-prior-spread (Whitaker & Hamill 2012): each parameter's posterior anomalies are rescaled by $\alpha(\sigma^{b} - \sigma^{a})/\sigma^{a} + 1$, where $\sigma^{b}$ and $\sigma^{a}$ are the background (pre-update) and analysis (post-update) standard deviations. Being per-parameter, it re-inflates the directions that collapsed hardest, so it is the spectrally-aware workhorse. "adaptive" is a global, magnitude-targeting scheme: it measures the mean spread-retention ratio $q = \mathrm{mean}_j(\sigma^{a}_j/\sigma^{b}_j)$ and, when q falls below retention_floor, applies a single multiplicative factor $\min(\texttt{max\_factor}, \texttt{retention\_floor}/q)$ to restore the lost variance magnitude. "multiplicative" applies a fixed factor every iteration. "none" disables inflation.

The companion diagnostic is the spectral spread-ESS (ensemble_spread_ess()), recorded each iteration regardless of method: it reports the effective number of variance-carrying directions and is what detects directional collapse. Because that participation ratio is invariant to a global rescaling, a global ("multiplicative" / "adaptive") inflation restores variance magnitude but not the spectral shape; "rtps" is the method that reshapes the spectrum. The two compose well.

Value

An object of class "pesto_inflation".

References

Whitaker, J.S. & Hamill, T.M. (2012). Evaluating methods to account for system errors in ensemble data assimilation. Monthly Weather Review, 140(9), 3078–3089.

See Also

pesto_localisation(), ensemble_spread_ess(), pesto_ies_callback().

Examples

inf <- pesto_inflation("rtps", alpha = 0.5)
inf

Description

Builds a control object that tells pesto_ies_callback() and pesto_ies_filter() how to taper the ensemble Kalman gain, suppressing the spurious long-range parameter-observation correlations that a finite ensemble manufactures. Localisation is applied as a Schur (elementwise) product on the explicit gain inside ensemble_solution_localised(); the default method = "none" leaves the standard SVD update untouched.

Usage

pesto_localisation(
  method = c("none", "correlation", "distance"),
  taper = c("hard", "soft"),
  threshold = -1,
  n_shuffle = 1L,
  quantile = 0.95,
  distances = NULL,
  par_coords = NULL,
  obs_coords = NULL,
  radius = NULL
)

Arguments

Details

"correlation" is the iterative-ensemble-smoother-native automatic localisation of Luo & Bhakta (2020): it needs no parameter or observation coordinates, estimating a noise floor from the ensemble itself and damping sample correlations that fall below it (see correlation_localisation()). This is the recommended default for parameter-estimation problems whose parameters carry no spatial metric. "distance" is classical distance-based localisation: a Gaspari-Cohn taper (gaspari_cohn()) of a parameter-to-observation distance matrix, for problems where such a metric exists — supply either distances directly or par_coords + obs_coords (Euclidean distances are then computed), together with radius.

Value

An object of class "pesto_localisation".

References

Luo, X. & Bhakta, T. (2020). Automatic and adaptive localization for ensemble-based history matching. Journal of Petroleum Science and Engineering, 184, 106559.

See Also

pesto_inflation(), correlation_localisation(), gaspari_cohn().

Examples

loc <- pesto_localisation("correlation", taper = "soft")
loc

Description

Bundles an ordered stack of pesto_forward_model() levels – cheapest (level = 0) to most expensive (level = n - 1) – with their relative per-evaluation costs. This is the first-class form of the bridge's ⁠(cheap, expensive)⁠ fidelity vector: the IES driver pesto_ies_callback() selects a level per iteration via fidelity_schedule, and mf_control_variate() debiases a cheap level against a sparse expensive sample for surrogate cascades.

Usage

pesto_multifidelity_model(levels, costs = NULL, label = NA_character_)

Arguments

Details

Each element of levels may be a bare ⁠function(theta) -> obs⁠ or a fully-specified pesto_forward_model(); bare functions are coerced via as_forward_model(). Levels must be ordered by ascending fidelity (cheapest first) – the convention the level index and the costs vector both follow.

Value

A pesto_multifidelity_model S7 object.

See Also

pesto_forward_model(), pesto_evaluate(), mf_control_variate(), pesto_ies_callback().

Examples

cheap     <- function(theta) theta %*% c(1, 1)        # fast, biased
expensive <- function(theta) theta %*% c(1, 1) + 0.5  # slow, truth
mf <- pesto_multifidelity_model(
  levels = list(
    pesto_forward_model(fn = cheap,     n_obs = 1L, fidelity = 0L),
    pesto_forward_model(fn = expensive, n_obs = 1L, fidelity = 1L)
  ),
  costs = c(1, 25)
)
theta <- matrix(c(1, 0, 0, 1), nrow = 2L, byrow = TRUE)
pesto_evaluate(mf, theta, level = 0L)  # cheap
pesto_evaluate(mf, theta, level = 1L)  # expensive

Description

Constructs the optional obs_schema slot of a pesto_ensemble_manifest: a machine-checkable statement of what each output and parameter column means (its physical quantity and unit), so a downstream consumer can verify two manifests are commensurable by name rather than trusting a positional convention. This is the single Independent-Oracle anchor the manifest format previously lacked.

Usage

pesto_obs_schema(outputs = NULL, params = NULL)

Arguments

Details

Each row optionally carries provenance: verified_on (a Date, or NA for an unverified fact), oracle_kind, and evidence_path. A column is grounded only when verified_on is a non-NA date.

Value

A validated list(outputs = , params = ) suitable for the obs_schema argument of pesto_ensemble_manifest. Either side may be omitted (NULL).

See Also

The manifest format specification in inst/manifest_contract.md.

Examples

# Describe two output columns and one parameter, grounding the yield
# column against a dated oracle and leaving the rest unverified.
schema <- pesto_obs_schema(
  outputs = data.frame(
    name        = c("yield", "biomass"),
    quantity    = c("grain yield", "above-ground biomass"),
    unit        = c("kg/ha", "kg/ha"),
    verified_on = c(as.Date("2026-06-01"), as.Date(NA)),
    stringsAsFactors = FALSE
  ),
  params = data.frame(
    name      = "rue",
    apsim_node = "Wheat.Leaf.Photosynthesis.RUE.FixedValue",
    unit      = "g/MJ",
    stringsAsFactors = FALSE
  )
)
str(schema, max.level = 2L)

Description

A textbook, pure-R implementation of one Iterative Ensemble Smoother (IES) parameter upgrade as published in Chen & Oliver (2013) eq. 12. This function is independent of ensemble_solution() (the C++ kernel) and exists for two purposes:

Usage

pesto_reference_ies(
  par_ensemble,
  obs_ensemble,
  obs_target,
  weights,
  lambda = 1
)

Arguments

Details

1. Independent numerical validation. PESTO's C++ kernel can be cross-checked against this reference at machine precision (the package's regression test tests/testthat/test-ensemble-solution-sign.R and the paired-seed protocol in inst/scripts/i2_paired_seed_check.R both rely on this).

2. Self-contained pestpp-ies comparison. Vignettes can compare PESTO native IES to this textbook reference without requiring the upstream pestpp-ies binary, making the comparison story reproducible on CRAN check farms and other binary-free environments.

The implementation is deliberately pedagogical — readability over speed. For production work use ensemble_solution() (the C++ kernel; $\sim 35\times$ faster on production-scale ensembles).

Value

A numeric matrix ⁠(n_par x n_real)⁠: the parameter upgrade (add to par_ensemble to get the next iterate).

Sign convention

This reference uses the textbook sign obs_resid = obs - sim with a positive leading sign in the upgrade. PESTO's C++ kernel uses the equivalent obs_resid = sim - obs with a leading negative sign (see ?ensemble_solution). Both yield identical upgrades; the difference is purely cosmetic. The two are cross-validated to machine precision in inst/scripts/i2_paired_seed_check.R.

References

Chen, Y. & Oliver, D. S. (2013). Levenberg-Marquardt forms of the iterative ensemble smoother for efficient history matching and uncertainty quantification. Computational Geosciences, 17(4), 689–703. doi:10.1007/s10596-013-9351-5

Evensen, G. (2018). Analysis of iterative ensemble smoothers for solving inverse problems. Computational Geosciences, 22(3), 885–908. doi:10.1007/s10596-018-9731-y

See Also

ensemble_solution() for the production C++ kernel.

Examples

# Simple linear inverse problem: identify the true theta given noisy obs.
set.seed(20260425L)
n_par <- 4L; n_obs <- 12L; n_real <- 30L
G <- matrix(rnorm(n_obs * n_par) / sqrt(n_par), n_obs, n_par)
theta_true <- rnorm(n_par)
y_obs <- as.numeric(G %*% theta_true) + rnorm(n_obs, sd = 0.05)
weights <- rep(1 / 0.05, n_obs)

par_ens <- matrix(rnorm(n_par * n_real, sd = 1), n_par, n_real)
obs_ens <- G %*% par_ens

# One textbook IES upgrade.
upg <- pesto_reference_ies(
  par_ensemble = par_ens,
  obs_ensemble = obs_ens,
  obs_target   = y_obs,
  weights      = weights,
  lambda       = 1.0
)
dim(upg)

# Apply the upgrade and check phi reduces.
par_next <- par_ens + upg
obs_next <- G %*% par_next
Y <- matrix(rep(y_obs, n_real), n_obs, n_real)
phi0 <- mean(compute_phi(Y - obs_ens,  weights))
phi1 <- mean(compute_phi(Y - obs_next, weights))
phi1 < phi0

Description

Executes pestpp-sen for Morris or Sobol sensitivity analysis.

Usage

pesto_sensitivity(
  pst_file,
  method = c("morris", "sobol"),
  exe = NULL,
  extra_args = list(),
  working_dir = NULL,
  verbose = TRUE
)

Arguments

Value

A list of class pesto_sen_result.

Examples

if (nzchar(Sys.which("pestpp-sen"))) {
  pars <- data.table::data.table(
    parnme = c("k1", "k2"), partrans = "log", parchglim = "factor",
    parval1 = c(1.0, 0.5), parlbnd = c(0.01, 0.001),
    parubnd = c(100, 50), pargp = "hydraulic"
  )
  obs <- data.table::data.table(
    obsnme = c("h1", "h2"), obsval = c(1.0, 2.0),
    weight = c(1.0, 1.0), obgnme = "head"
  )
  pst <- create_pest_scenario(pars, obs, model_command = "echo run")
  tf <- tempfile(fileext = ".pst")
  on.exit(unlink(tf), add = TRUE)
  write_pst(pst, tf)
  res <- pesto_sensitivity(tf, method = "morris", verbose = FALSE)
  res$method
}

Description

Performs a single IES iteration using a Gaussian Process surrogate to reduce the number of expensive full-model evaluations.

Usage

pesto_surrogate_ies(
  par_ensemble,
  obs_ensemble,
  obs_target,
  weights,
  parcov_inv,
  lambda = 1,
  uncertainty_threshold = 0.1,
  eigthresh = 1e-06
)

Arguments

Details

How it works:

1. A GP surrogate is trained on existing parameter-observation pairs

2. The surrogate predicts model outputs for all ensemble members

3. Only members with high prediction uncertainty trigger full model runs

4. Control-variate bias correction blends surrogate and model results

5. Standard IES update is computed on the blended ensemble

This typically saves 50-90% of model evaluations per iteration.

Value

A list containing:

upgrade

Matrix of parameter upgrades

n_model_runs

Number of full model evaluations needed

n_surrogate_runs

Number of surrogate-only evaluations

savings_pct

Percentage of model runs saved

gp_diagnostics

GP training diagnostics

References

Rasmussen, C.E. & Williams, C.K.I. (2006). Gaussian Processes for Machine Learning. MIT Press.

Liu, F. & Guillas, S. (2017). Dimension reduction for Gaussian process emulation. Statistics and Computing, 27(3), 785-802.

Examples

set.seed(7L)
n_real <- 15L; n_par <- 5L; n_obs <- 8L
par_ens <- matrix(rnorm(n_real * n_par), n_real, n_par,
                  dimnames = list(NULL, paste0("k", 1:n_par)))
obs_ens <- matrix(rnorm(n_real * n_obs), n_real, n_obs,
                  dimnames = list(NULL, paste0("h", 1:n_obs)))
obs_target <- rnorm(n_obs)
weights    <- rep(1.0, n_obs)
parcov_inv <- rep(1.0, n_par)
res <- pesto_surrogate_ies(
  par_ensemble = par_ens,
  obs_ensemble = obs_ens,
  obs_target   = obs_target,
  weights      = weights,
  parcov_inv   = parcov_inv,
  lambda       = 1.0
)
res$savings_pct

Description

Executes pestpp-swp for embarrassingly parallel model runs across a parameter ensemble.

Usage

pesto_sweep(
  pst_file,
  par_ensemble,
  exe = NULL,
  working_dir = NULL,
  verbose = TRUE
)

Arguments

Value

A list containing observation outputs for each realisation.

Examples

if (nzchar(Sys.which("pestpp-swp"))) {
  pars <- data.table::data.table(
    parnme = c("k1", "k2"), partrans = "log", parchglim = "factor",
    parval1 = c(1.0, 0.5), parlbnd = c(0.01, 0.001),
    parubnd = c(100, 50), pargp = "hydraulic"
  )
  obs <- data.table::data.table(
    obsnme = c("h1", "h2"), obsval = c(1.0, 2.0),
    weight = c(1.0, 1.0), obgnme = "head"
  )
  pst <- create_pest_scenario(pars, obs, model_command = "echo run")
  tf <- tempfile(fileext = ".pst")
  on.exit(unlink(tf), add = TRUE)
  write_pst(pst, tf)
  par_ens <- data.table::data.table(
    real_name = c("r1", "r2", "r3"),
    k1 = c(0.8, 1.0, 1.2),
    k2 = c(0.4, 0.5, 0.6)
  )
  res <- pesto_sweep(tf, par_ensemble = par_ens, verbose = FALSE)
  res$exit_code
}

Description

Returns the PESTO package version, plus the version of the PEST++ install PESTO resolves (see pestpp_available() for how it is found). PESTO bundles no binaries, so pestpp_version reports whatever install is configured, or "not found" when there is none.

Usage

pesto_version()

Value

A list with version strings: pesto_version, pestpp_version, platform, and r_version.

Examples

v <- pesto_version()
v$pesto_version
v$platform

Description

A non-erroring capability probe for a PEST++ family executable. It is the documented way for examples, vignettes, and conditional tests to skip gracefully when no binary is installed: every PESTO algorithm runs natively in R, so the external binaries are needed only for the optional cross-checking and .pst-file paths.

Usage

pestpp_available(which = "pestpp-ies")

Arguments

Details

The probe resolves the executable exactly as pesto_ies() and friends do – the per-tool environment variable (e.g. PESTPP_IES_EXE_PATH), then PESTPP_BIN_DIR, then the system PATH – but never throws. PESTO bundles no binaries: PEST++ ships no installer and is not on the PATH by default, so pointing at an existing install by environment variable is the normal way to reach it from R.

Value

A length-one logical: TRUE if the named executable is resolvable, FALSE otherwise. Never errors.

See Also

pesto_ies(), pesto_glm()

Examples

# FALSE on a machine without PEST++ installed -- and that is fine:
pestpp_available()
pestpp_available("pestpp-glm")

Description

Visualises the prior and/or posterior parameter ensemble distributions as violin plots or density ridges.

Usage

plot_ensemble(
  ensemble,
  parameters = NULL,
  prior_ensemble = NULL,
  max_params = 20L,
  title = "Parameter Ensemble Distributions"
)

Arguments

Value

A ggplot2 object.

Examples

posterior <- data.table::data.table(
  real_name = sprintf("real_%02d", 1:50),
  k1 = rnorm(50, 1.0, 0.15),
  k2 = rnorm(50, 0.5, 0.05),
  k3 = rnorm(50, 2.0, 0.40),
  k4 = rnorm(50, 0.1, 0.02)
)
prior <- data.table::data.table(
  real_name = sprintf("real_%02d", 1:50),
  k1 = rnorm(50, 1.0, 0.50),
  k2 = rnorm(50, 0.5, 0.20),
  k3 = rnorm(50, 2.0, 1.00),
  k4 = rnorm(50, 0.1, 0.08)
)
p <- plot_ensemble(posterior, prior_ensemble = prior)
inherits(p, "ggplot")

Description

Creates a ranked lollipop plot of parameter identifiability based on the singular value decomposition of the Jacobian matrix. Accepts either a numeric Jacobian matrix in memory or a path to a .jco (PEST binary) file. High-dimensional problems are kept legible by showing only the most identifiable parameters (see top_n).

Usage

plot_identifiability(
  jacobian = NULL,
  jco_file = NULL,
  pst = NULL,
  n_sv = NULL,
  top_n = NULL,
  title = "Parameter Identifiability"
)

Arguments

Value

A ggplot2 object.

Examples

J <- matrix(rnorm(30 * 8), nrow = 30, ncol = 8)
J[, 7] <- 0.5 * J[, 1] + 0.5 * J[, 2]
J[, 8] <- 1e-6 * rnorm(30)
colnames(J) <- paste0("k", 1:8)
p <- plot_identifiability(jacobian = J)
inherits(p, "ggplot")

Description

Creates a publication-quality plot of objective function values across iterations, showing mean, min, max, and individual realisation traces.

Usage

plot_phi(
  result,
  log_scale = TRUE,
  show_reals = FALSE,
  title = "Objective Function Convergence"
)

Arguments

Value

A ggplot2 object.

Examples

phi_dt <- data.table::data.table(
  iteration = 0:4,
  total_runs = c(50L, 100L, 150L, 200L, 250L),
  mean = c(1200, 450, 180, 95, 72),
  min  = c(900, 320, 130, 70, 55),
  max  = c(1700, 680, 260, 140, 105),
  median = c(1180, 440, 175, 92, 70),
  std    = c(180, 80, 35, 18, 12)
)
p <- plot_phi(phi_dt, log_scale = TRUE,
              title = "Synthetic Phi Convergence")
inherits(p, "ggplot")

Description

Visualises the surrogate-accelerated IES performance including model savings, uncertainty distribution, and GP quality metrics.

Usage

plot_surrogate_diagnostics(results, title = "Surrogate IES Diagnostics")

Arguments

Value

A ggplot2 object.

Examples

iter1 <- list(n_model_runs = 12L, n_surrogate_runs = 38L,
              savings_pct = 76.0, mean_uncertainty = 0.18)
iter2 <- list(n_model_runs = 8L,  n_surrogate_runs = 42L,
              savings_pct = 84.0, mean_uncertainty = 0.11)
iter3 <- list(n_model_runs = 5L,  n_surrogate_runs = 45L,
              savings_pct = 90.0, mean_uncertainty = 0.07)
p <- plot_surrogate_diagnostics(list(iter1, iter2, iter3))
inherits(p, "ggplot")

Description

Generates predictions and prediction uncertainties for new parameter sets using a trained GP surrogate. The uncertainty estimates are crucial for the adaptive switching criterion.

Usage

predict_gp_surrogate(gp, X_new)

Arguments

Value

A list with:

mean

Matrix (m x nobs). Predicted observations.

variance

Matrix (m x nobs). Prediction variance per output.

uncertainty

Numeric vector (m). Mean prediction uncertainty per realisation.

Examples

set.seed(1L)
X_train <- matrix(rnorm(20 * 4), 20, 4)
Y_train <- matrix(rnorm(20 * 6), 20, 6)
gp      <- train_gp_surrogate(X_train, Y_train)
X_new   <- matrix(rnorm(5 * 4), 5, 4)
pred    <- predict_gp_surrogate(gp, X_new)
dim(pred$mean)
length(pred$uncertainty)

Description

Companion predictor for train_gp_surrogate_tuned(). It reapplies the per-axis pre-scaling the tuner stored (so an anisotropic surrogate predicts on the same geometry it was fitted on) and adds the centred response mean back, returning predictions on the original response scale.

Usage

predict_gp_surrogate_tuned(gp, X_new)

Arguments

Value

The list predict_gp_surrogate() returns (mean, variance, uncertainty), with mean on the original response scale.

See Also

train_gp_surrogate_tuned().

Examples

set.seed(1L)
X <- matrix(runif(40L * 2L), 40L, 2L)
y <- sin(3 * X[, 1]) + 0.5 * X[, 2]^2
gp <- train_gp_surrogate_tuned(X, matrix(y, ncol = 1L))
pred <- predict_gp_surrogate_tuned(gp, X)

Description

Predict with RFF Sparse GP Surrogate

Usage

predict_rff_surrogate(rff, X_new)

Arguments

Value

A list with mean predictions and approximate uncertainties.

Examples

set.seed(1L)
X_train <- matrix(rnorm(30 * 4), 30, 4)
Y_train <- matrix(rnorm(30 * 6), 30, 6)
rff     <- train_rff_surrogate(X_train, Y_train, n_features = 100L)
X_new   <- matrix(rnorm(5 * 4), 5, 4)
pred    <- predict_rff_surrogate(rff, X_new)
dim(pred$mean)
length(pred$uncertainty)

Description

Print method for pesto_pst objects

Usage

## S3 method for class 'pesto_pst'
print(x, ...)

Arguments

Value

Invisibly returns x. Called for the side effect of printing.

Examples

pars <- data.table::data.table(
  parnme = c("k1", "k2"), partrans = "log", parchglim = "factor",
  parval1 = c(1.0, 0.5), parlbnd = c(0.01, 0.001),
  parubnd = c(100, 50), pargp = "hydraulic"
)
obs <- data.table::data.table(
  obsnme = c("h1", "h2"), obsval = c(1.0, 2.0),
  weight = c(1.0, 1.0), obgnme = "head"
)
pst <- create_pest_scenario(pars, obs, model_command = "echo run")
print(pst)

Description

Reads PEST++ ensemble files in CSV or binary (.jcb/.jco) format.

Usage

read_ensemble(file, format = c("csv", "binary"))

Arguments

Value

A data.table with realisations as rows and parameters/observations as columns. The first column real_name contains realisation names.

See Also

write_ensemble()

Examples

ens <- data.table::data.table(
  real_name = sprintf("real_%02d", 1:10),
  k1 = rnorm(10, mean = 1.0, sd = 0.2),
  k2 = rnorm(10, mean = 0.5, sd = 0.1),
  k3 = rnorm(10, mean = 2.0, sd = 0.3)
)
tf <- tempfile(fileext = ".csv")
on.exit(unlink(tf), add = TRUE)
write_ensemble(ens, tf)
ens_back <- read_ensemble(tf, format = "csv")
identical(names(ens_back), names(ens))
nrow(ens_back) == nrow(ens)

Description

Inverse of write_manifest(). Reads the YAML, loads the three sidecar data files (paths resolved relative to the YAML file), and reconstructs the pesto_ensemble_manifest S7 object. The file extensions in the YAML's ⁠artefacts:⁠ block determine the read path (.rds via readRDS, .csv via utils::read.csv).