BO(f, modelSettings, optimizationSettings, warmStart)

Performs a full Bayesian Optimization run. ... (resten av docstringen) ...

Returns

• A Tuple containing:
  • gp: The final, trained GP model.
  • X: The full n x d matrix of all evaluated points.
  • y: The full n-element vector of all observations.
  • objectMaximizer: The location x of the global maximum of the final posterior mean.
  • objectMaximizerY: The predicted posterior mean value at objectMaximizer. # <– NEW
  • postMaxObserved: The location x of the observed point with the highest posterior mean.
  • postMaxObservedY: The predicted posterior mean value at postMaxObserved.

...

ExpectedMaxGaussian(μ::Vector{Float64}, σ::Vector{Float64})

Analytically computes the expectation E[max(μ + σZ)] where Z ~ N(0,1).

This is the core for knowledgeGradientDiscrete. It calculates the expected maximum of a set of correlated Gaussian random variables that share a single source of randomness, Z. The function is robust to cases where slopes (σ) are equal or nearly equal.

Arguments

• μ::Vector{Float64}: A vector of the current means of the random variables.
• σ::Vector{Float64}: A vector of the sensitivities ("slopes") of each random variable with respect to the common random factor Z.

Returns

• Float64: The analytically computed expected maximum value.

estimate_integral_wsabi(gp, bounds; n_samples=100_000, y_mean=0.0, y_std=1.0)

Estimates the integral of the original function f(x) using the final GP model trained on warped data g(x) = log(f(x)).

It uses Monte Carlo integration on the posterior expectation of f(x). E[f(x)] = exp(μg(x) + σ²g(x)/2), where μg and σ²g are the posterior mean and variance of the GP fitted to the log-transformed data.

Arguments

• gp: The final trained GP model.
• bounds: A tuple (lo, hi) defining the integration domain.
• n_samples: Number of Monte Carlo samples for the approximation.
• y_mean, y_std: The mean and std dev used to standardize the warped y-values, needed to un-scale the GP's predictions.

Returns

• Float64: The estimated value of the integral.

  expected_improvement(gp, xnew, fMax; ξ = 0.01)

Computes the expected improvement given a Gaussian process (gp) object, and the earliest best evaluation (fMax) at the new evaluation point xnew). ξ is a tuning parameter that controls the exploration-exploitation trade-off. A large value of ξ encourages exploration and vice versa.

Returns the expected improvement at xnew.

Examples

julia> Set up GP model.
julia> X_train = [1.0, 2.5, 4.0]; julia> y_train = [sin(x) for x in X_train];
julia> gp_model = GP(X_train', y_train, MeanZero(), SE(0.0, 0.0));
julia> optimize!(gp_model);

julia> # 2. Define the best observed value and a candidate point
julia> fMax = maximum(y_train);
julia> x_candidate = 3.0;

julia> # 3. Compute Expected Improvement
julia> ei = expected_improvement(gp_model, x_candidate, fMax; ξ=0.01)
≈ 4.89.

expected_improvement_boundary(gp, xnew, ybest; ξ = 0.10, bounds=nothing)

Computes the Expected Improvement (EI) with an optional penalty for points near the boundary of the search space. VERY EXPERIMENTAL.

The penalty is a multiplicative weight w = prod(1 - x_norm^2) that smoothly pushes the search away from the edges of the scaled [-1, 1] domain. This can be useful to prevent the optimizer from selecting points at the very edge of the feasible region. The idea comes from practical experience that EI otherwise spent too much time exploring the boundaries.

Arguments

• gp: The trained Gaussian Process model.
• xnew: The candidate point at which to evaluate the EI.
• ybest: The current best observed value.
• ξ: The exploration-exploitation trade-off parameter.
• bounds: A tuple (lower, upper) defining the original search space, required for the penalty.

Returns

• Float64: The (potentially penalized) Expected Improvement value.

    first_func(x)

A first dummy function that returns the log of the square of the input value.

julia> T = 500; ϕ = [sin.(2π*(1:T)/T) -0.1*ones(T)]; θ = zeros(T); σ = 1; μ = 2;
julia> y = simTvARMA(ϕ, θ, μ, σ);
julia>  tvPeriodogram(y, 25, 15)

inv_rescale(X_scaled, lo, hi)

Performs the inverse of rescale, scaling an n x d matrix X_scaled from the hypercube [-1, 1]^d back to the original domain.

Arguments

• X_scaled: An n x d matrix where each row is a point in [-1, 1]^d.
• lo: A d-element vector of lower bounds for the original domain.
• hi: A d-element vector of upper bounds for the original domain.

Returns

• The un-scaled n x d matrix in the original domain.

Examples

julia> lo = [0.0, 10.0];
julia> hi = [1.0, 20.0];
julia> X_scaled = [0.0 0.0; -1.0 -1.0];
julia> inv_rescale(X_scaled, lo, hi)
2×2 Matrix{Float64}:
 0.5  15.0
 0.0  10.0

knowledgeGradientDiscrete(gp, xnew, domain_points)

Computes the Knowledge Gradient (KG) acquisition function where the future maximum is constrained to occur on a fixed, discrete set of points.

This function is the analytical engine for the knowledgeGradientHybrid heuristic.

Arguments

• gp: The trained Gaussian Process model.
• xnew: The candidate point at which to evaluate the KG.
• domain_points: A d x M matrix of M discrete points in the scaled [-1, 1] space where the future maximum is sought.

Returns

• Float64: The discrete Knowledge Gradient value.

knowledgeGradientHybrid(gp, xnew, lower, upper; n_z=5)

Calculates the Hybrid Knowledge Gradient (KGh) acquisition function for a candidate point xnew.

The Hybrid KG combines the strengths of the Monte-Carlo and Discrete KG methods. It uses a small, deterministic set of n_z fantasy scenarios to identify a high-potential set of future maximizers (X_MC), and then uses the fast, analytical Discrete KG algorithm on that small set.

Arguments

• gp::GPE: The trained Gaussian Process model.
• xnew: The candidate point at which to evaluate the KG.
• lower, upper: The bounds of the search domain for finding the fantasy maximizers.
• n_z::Int: The number of deterministic Z-samples to use (default is 5, as in the paper).

Returns

• Float64: The positive Hybrid Knowledge Gradient value.

knowledgeGradientMonteCarlo(gp, xnew, lower, upper; n_samples=20)

Calculates the Knowledge Gradient (KG) acquisition function for a candidate point xnew using Monte Carlo. This version gives a noisy surface to the optimization landscape due to MC variation. Try the "quadrature methods below".

The Knowledge Gradient quantifies the expected increase in the maximum estimated value after sampling at xnew.

Arguments

• gp::GPE: The trained Gaussian Process model.
• xnew: The candidate point at which to evaluate the KG.
• n_samples::Int: The number of Monte Carlo samples to use.

Returns

• Float64: The positive Knowledge Gradient value. This function returns the natural KG score, which should be maximized by the optimization loop.

multi_start_maximize(f, lower, upper; n_starts=20)

Finds the global maximum of a function f within a given box [lower, upper] by using multiple starting points. This is a wrapper around multi_start_minimize that simply minimizes the negative of the function.

Arguments

• f: The function to maximize.
• lower: A vector of lower bounds.
• upper: A vector of upper bounds.
• n_starts::Int: The number of distinct starting points to use.

Returns

• Tuple{Float64, Vector{Float64}}: A tuple (best_max, best_argmax) containing the maximum value found and the location where it was found.

multi_start_minimize(f, lower, upper; n_starts=20)

Finds the global minimum of a function f within a given box [lower, upper] by using multiple starting points.

This function uses the L-BFGS optimizer from Optim.jl together with autodiff, starting from n_starts evenly spaced points within the domain to increase the probability of finding a global, rather than local, minimum.

Arguments

• f: The function to minimize.
• lower: A vector of lower bounds.
• upper: A vector of upper bounds.
• n_starts::Int: The number of distinct starting points to use.

Returns

• Tuple{Float64, Vector{Float64}}: A tuple (best_min, best_argmin) containing the minimum value found and the location where it was found.

posteriorMax(gp; n_starts=20)

Finds the global maximum of the GP's posterior mean over the entire continuous, scaled domain [-1, 1]^d using multi-start optimization. This maximum can occur at a previously unobserved point.

Arguments

• gp: The trained Gaussian Process model.
• n_starts::Int: The number of restarts for the optimization.

Returns

• NamedTuple{(:fX_max, :X_max)}: A named tuple containing the maximum value of the posterior mean (fX_max) and its location (X_max) in the scaled space.

posteriorMaxObs(gp, X_scaled)

Finds the maximum of the GP's posterior mean evaluated only at the points that have already been observed.

Arguments

• gp: The trained Gaussian Process model.
• X_scaled: A d x n matrix of the n locations already observed, in the scaled space.

Returns

• Float64: The maximum value of the posterior mean among the observed points.

posterior_cov(gp::GPE, X1::AbstractMatrix, X2::AbstractMatrix)

Computes the posterior covariance matrix k_n(X1, X2) between two sets of points X1 and X2, given the GP's training data.

The calculation uses the formula k(X1, X2) - k(X1, X) * K_inv * k(X, X2), using the pre-computed Cholesky factorization of the GP's kernel matrix for high efficiency.

Arguments

• gp::GPE: The trained Gaussian Process model.
• X1::AbstractMatrix: A d x n1 matrix of n1 points.
• X2::AbstractMatrix: A d x n2 matrix of n2 points.

Returns

• Matrix{Float64}: The n1 x n2 posterior covariance matrix.

posterior_variance(gp, xnew)

Computes the posterior variance of the Gaussian Process gp at a new point xnew.

This acquisition function directs sampling to regions of highest uncertainty in the model. It is the primary acquisition function for standard Bayesian Quadrature methods like WSABI, but can also be used for pure exploration in Bayesian Optimization.

Arguments

• gp: The trained Gaussian Process model.
• xnew: The candidate point at which to evaluate the posterior variance.

Returns

• Float64: The posterior variance at xnew.

Examples

```julia-repl julia> Xtrain = [1.0, 5.0]; julia> ytrain = [sin(x) for x in Xtrain]; julia> gp = GP(Xtrain', y_train, MeanZero(), SE(0.0, 0.0)); julia> optimize!(gp);

julia> # The point of highest uncertainty is halfway between the samples julia> xcandidate = 3.0; julia> posteriorvariance(gp, x_candidate) 1.0000010000000002

propose_next(gp::GPE, f_max; n_restarts::Int, acq_config::AcquisitionConfig)

Optimizes the acquisition function to find the best next point to sample. The optimization is performed in the scaled [-1, 1]^d space.

Arguments

• gp::GPE: The current trained Gaussian Process model.
• f_max: The current best observed value (or posterior mean max), used by some acquisition functions like EI.
• n_restarts::Int: The number of restarts for the multi-start optimization of the acquisition function.
• acq_config::AcquisitionConfig: A struct holding the configuration for the chosen acquisition function (e.g., EIConfig, UCBConfig, KGQConfig).

Returns

• Vector{Float64}: The coordinates of the next best point to sample, in the scaled [-1, 1]^d space.

rescale(X, lo, hi)

Scales an n x d matrix X from an original domain to the hypercube [-1, 1]^d.

Arguments

• X: An n x d matrix where each row is a point in the original domain.
• lo: A d-element vector of lower bounds for the original domain.
• hi: A d-element vector of upper bounds for the original domain.

Returns

• The scaled n x d matrix where all values are in [-1, 1].

Examples

```julia-repl julia> lo = [0.0, 10.0]; julia> hi = [1.0, 20.0]; julia> X = [0.5 15.0; 0.0 10.0]; julia> rescale(X, lo, hi) 2×2 Matrix{Float64}: 0.0 0.0 -1.0 -1.0

  upper_confidence_bound(gp, xnew; κ = 2.0)

Computes the upper confidence bound criteria at the new point xnew given a Gaussian process (gp) object, and the exploitation/exploitation parameter κ. High values of κ encourage exploration.

Returns the upper confidence bound criteria at xnew.

Examples

julia> Set up GP model.
julia> X_train = [1.0, 2.5, 4.0]; julia> y_train = [sin(x) for x in X_train];
julia> gp_model = GP(X_train', y_train, MeanZero(), SE(0.0, 0.0));
julia> optimize!(gp_model);
julia> x_candidate = 3.0;

julia> ucb = upper_confidence_bound(gp_model, x_candidate, κ = 2.0)
-0.22727575567547253   # Updete this example, it is from when i used minimum instead of maximum.