From a930770b7ff25a45693b8c03a58797cc7a6a7c65 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Wed, 29 Jul 2026 12:33:36 -0400 Subject: [PATCH 01/10] Initial --- Project.toml | 2 + ext/InfiniteDisjunctiveProgramming.jl | 104 +++++++++++++++++++++++--- 2 files changed, 97 insertions(+), 9 deletions(-) diff --git a/Project.toml b/Project.toml index bdec9be..5de25a7 100644 --- a/Project.toml +++ b/Project.toml @@ -4,7 +4,9 @@ authors = ["hdavid16 "] version = "0.6.1" [deps] +AbstractGPs = "99985d1d-32ba-4be9-9821-2ec096f28918" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +KernelFunctions = "ec8451be-7e33-11e9-00cf-bbf324bd1392" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" [weakdeps] diff --git a/ext/InfiniteDisjunctiveProgramming.jl b/ext/InfiniteDisjunctiveProgramming.jl index 5f77dce..a251067 100644 --- a/ext/InfiniteDisjunctiveProgramming.jl +++ b/ext/InfiniteDisjunctiveProgramming.jl @@ -3,6 +3,7 @@ module InfiniteDisjunctiveProgramming import JuMP.MOI as _MOI import InfiniteOpt, JuMP import DisjunctiveProgramming as DP +import AbstractGPs, KernelFunctions ################################################################################ # MODEL @@ -263,9 +264,95 @@ function _interpolate_at( ) end -# Transcribe mini_expr, solve per support on the transcribed JuMP -# model, and aggregate to a scalar if uniform, else to a parameter -# function on main. +# ------ GP active-learning M(d) (hard-coded for the mbm_gp experiments) ------ +# Coordinates of every support, normalized to [0,1]^d so one isotropic +# lengthscale works across dimensions. +function _gp_support_coords(grids) + idxs = CartesianIndices(length.(grids)) + los = [minimum(g) for g in grids] + rng = [max(maximum(g) - minimum(g), eps()) for g in grids] + X = [[(grids[d][I[d]] - los[d]) / rng[d] for d in 1:length(grids)] + for I in idxs] + return X, collect(idxs) +end + +# Posterior mean and sd at all coords, given solved (index => M) samples. +# Outputs are standardized so the k*sd term scales with the M spread. +function _gp_mean_sd(X, solved) + lis = collect(keys(solved)) + yt = [solved[li] for li in lis] + ybar = sum(yt) / length(yt) + ystd = max(sqrt(sum(abs2, yt .- ybar) / max(length(yt) - 1, 1)), 1e-8) + kern = KernelFunctions.with_lengthscale( + KernelFunctions.SqExponentialKernel(), 0.1) + post = AbstractGPs.posterior( + AbstractGPs.GP(kern)(X[lis], 1e-8), (yt .- ybar) ./ ystd) + mz = AbstractGPs.mean(post, X) + vz = max.(AbstractGPs.var(post, X), 0.0) + return mz .* ystd .+ ybar, sqrt.(vz) .* ystd +end + +# Count of per-support M subproblems solved (for the grid-vs-GP comparison). +const _M_SOLVE_COUNT = Ref(0) + +# Original workflow: solve M at every support (grid). Kept for the mbm_gp +# vs grid comparison, gated by ENV["DP_MBM_GRID"]. +function _grid_M_vals(objectives, inner_sub, method) + M_vals = Array{Float64}(undef, size(objectives)) + for I in eachindex(objectives) + _M_SOLVE_COUNT[] += 1 + m = DP.raw_M(inner_sub, objectives[I], method) + m === nothing && return nothing + M_vals[I] = m + end + return M_vals +end + +# Solve M at actively-selected supports (max-UCB acquisition) and fill the +# rest with the UCB (mean + k*sd), a valid over-estimate. Returns a scalar +# when M is uniform (e.g. dependent multi-dim parameters where M does not +# vary), matching the grid workflow's early return. +function _gp_active_M_vals(objectives, inner_sub, method, sub, mini_expr) + idxs = collect(CartesianIndices(objectives)) + n = length(idxs) + k = 2.5 + budget = min(n, max(6, cld(n, 4))) + solved = Dict{Int, Float64}() + solve_at!(li) = begin + _M_SOLVE_COUNT[] += 1 + m = DP.raw_M(inner_sub, objectives[idxs[li]], method) + m === nothing && return false + solved[li] = m + true + end + for s in unique([1, cld(n + 1, 2), n]) + solve_at!(s) || return nothing + end + seed = collect(values(solved)) + all(==(first(seed)), seed) && return first(seed) # uniform M -> scalar + mini_prefs = InfiniteOpt.parameter_refs(mini_expr) + reverse_map = Dict(ws[1] => v for (v, ws) in sub.fwd_map) + prefs = Tuple(reverse_map[p] for p in mini_prefs) + grids = Tuple(InfiniteOpt.supports(p) for p in prefs) + X, _ = _gp_support_coords(grids) + while length(solved) < budget + ms, ss = _gp_mean_sd(X, solved) + acq = ms .+ k .* ss + for li in keys(solved) + acq[li] = -Inf + end + solve_at!(argmax(acq)) || return nothing + end + ms, ss = _gp_mean_sd(X, solved) + M_vals = Array{Float64}(undef, size(objectives)) + for (li, I) in enumerate(idxs) + M_vals[I] = get(solved, li, ms[li] + k * ss[li]) + end + return M_vals +end + +# Transcribe mini_expr, then approximate M(d) with an actively-sampled GP +# instead of solving at every support; aggregate to a scalar if uniform. function DP.raw_M( sub::DP.GDPSubmodel{<:InfiniteOpt.InfiniteModel}, mini_expr::JuMP.AbstractJuMPScalar, @@ -276,12 +363,11 @@ function DP.raw_M( inner_sub = DP.GDPSubmodel(transcribed,JuMP.VariableRef[], Dict{JuMP.VariableRef, Vector{JuMP.VariableRef}}() ) - M_vals = Array{typeof(method.default_M)}(undef, size(objectives)) - for I in eachindex(objectives) - m = DP.raw_M(inner_sub, objectives[I], method) - m === nothing && return nothing - M_vals[I] = m - end + M_vals = haskey(ENV, "DP_MBM_GRID") ? + _grid_M_vals(objectives, inner_sub, method) : + _gp_active_M_vals(objectives, inner_sub, method, sub, mini_expr) + M_vals === nothing && return nothing + M_vals isa Number && return M_vals all(==(first(M_vals)), M_vals) && return first(M_vals) mini_prefs = InfiniteOpt.parameter_refs(mini_expr) reverse_map = Dict(ws[1] => v for (v, ws) in sub.fwd_map) From ed08522cbfc1466d71845ffcda6eb85b07174166 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Sat, 1 Aug 2026 21:54:27 -0400 Subject: [PATCH 02/10] Simplify and code coverage --- Project.toml | 9 +- README.md | 3 + ext/InfiniteDisjunctiveProgramming.jl | 124 ++++--------- ext/InfiniteGPDisjunctiveProgramming.jl | 135 +++++++++++++++ src/datatypes.jl | 17 +- src/extension_api.jl | 50 ++++++ src/mbm.jl | 3 +- .../InfiniteDisjunctiveProgramming.jl | 33 +++- .../InfiniteGPDisjunctiveProgramming.jl | 163 ++++++++++++++++++ test/runtests.jl | 1 + 10 files changed, 438 insertions(+), 100 deletions(-) create mode 100644 ext/InfiniteGPDisjunctiveProgramming.jl create mode 100644 test/extensions/InfiniteGPDisjunctiveProgramming.jl diff --git a/Project.toml b/Project.toml index 5de25a7..ba1a492 100644 --- a/Project.toml +++ b/Project.toml @@ -4,20 +4,23 @@ authors = ["hdavid16 "] version = "0.6.1" [deps] -AbstractGPs = "99985d1d-32ba-4be9-9821-2ec096f28918" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" -KernelFunctions = "ec8451be-7e33-11e9-00cf-bbf324bd1392" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" [weakdeps] +AbstractGPs = "99985d1d-32ba-4be9-9821-2ec096f28918" InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57" +KernelFunctions = "ec8451be-7e33-11e9-00cf-bbf324bd1392" [extensions] InfiniteDisjunctiveProgramming = "InfiniteOpt" +InfiniteGPDisjunctiveProgramming = ["InfiniteOpt", "AbstractGPs", "KernelFunctions"] [compat] +AbstractGPs = "0.5" Aqua = "0.8" JuMP = "1.18" +KernelFunctions = "0.10" Reexport = "1" julia = "1.10" Juniper = "0.9.3" @@ -32,4 +35,4 @@ Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" Juniper = "2ddba703-00a4-53a7-87a5-e8b9971dde84" [targets] -test = ["Aqua", "HiGHS", "Test", "Juniper", "Ipopt", "InfiniteOpt"] +test = ["Aqua", "HiGHS", "Test", "Juniper", "Ipopt", "InfiniteOpt", "AbstractGPs", "KernelFunctions"] diff --git a/README.md b/README.md index 543d1f7..c3cd037 100644 --- a/README.md +++ b/README.md @@ -176,6 +176,7 @@ The following reformulation methods are currently supported: - `optimizer`: Optimizer to use when solving subproblems to determine M values. This is a required value. - `default_M`: Default big-M value to use if no big-M is specified for a logical variable (1e9). + - `M_sampler`: Strategy for computing M values across the supports of an infinite model. Default: `:auto`, which uses a Gaussian-process sampler (`GPSampler`) when AbstractGPs is loaded and otherwise solves an M subproblem at every support (`:exact`). Ignored for finite models. 5. [P-Split](https://arxiv.org/abs/2202.05198): This method reformulates each disjunct constraint into P constraints, each with a partitioned group defined by the user. This method requires that terms in the constraint be convex additively seperable with respect to each variable. The `PSplit` struct is created with the following required arguments: @@ -223,6 +224,8 @@ optimize!(model, gdp_method = Hull()) value(W) ``` +When the `MBM` reformulation is used on an infinite model, an M subproblem is solved at every support by default. Loading [AbstractGPs.jl](https://github.com/JuliaGaussianProcesses/AbstractGPs.jl) (`using AbstractGPs`) activates an additional extension that instead solves M at a subset of the supports and fills the rest with a conservative Gaussian-process estimate, which can substantially reduce the number of subproblem solves. The filled values are heuristic upper estimates rather than certificates; see the `GPSampler` docstring for the tuning knobs (`kappa`, `budget`) and use `MBM(optimizer, M_sampler = :exact)` to force exact solves. + ## Release Notes Prior to `v0.4.0`, the package did not leverage the JuMP extension capabilities and was not as robust. For these earlier releases, refer to [Perez, Joshi, and Grossmann, 2023](https://arxiv.org/abs/2304.10492v1) and the following [JuliaCon 2022 Talk](https://www.youtube.com/watch?v=AMIrgTTfUkI). diff --git a/ext/InfiniteDisjunctiveProgramming.jl b/ext/InfiniteDisjunctiveProgramming.jl index a251067..2055422 100644 --- a/ext/InfiniteDisjunctiveProgramming.jl +++ b/ext/InfiniteDisjunctiveProgramming.jl @@ -3,7 +3,6 @@ module InfiniteDisjunctiveProgramming import JuMP.MOI as _MOI import InfiniteOpt, JuMP import DisjunctiveProgramming as DP -import AbstractGPs, KernelFunctions ################################################################################ # MODEL @@ -264,116 +263,61 @@ function _interpolate_at( ) end -# ------ GP active-learning M(d) (hard-coded for the mbm_gp experiments) ------ -# Coordinates of every support, normalized to [0,1]^d so one isotropic -# lengthscale works across dimensions. -function _gp_support_coords(grids) - idxs = CartesianIndices(length.(grids)) - los = [minimum(g) for g in grids] - rng = [max(maximum(g) - minimum(g), eps()) for g in grids] - X = [[(grids[d][I[d]] - los[d]) / rng[d] for d in 1:length(grids)] - for I in idxs] - return X, collect(idxs) +# Resolve `:auto` to the GP sampler when the GP extension is loaded +function _resolve_M_sampler(sampler) + sampler === :auto || return sampler + gp = Base.get_extension(DP, :InfiniteGPDisjunctiveProgramming) + return isnothing(gp) ? :exact : DP.GPSampler() end -# Posterior mean and sd at all coords, given solved (index => M) samples. -# Outputs are standardized so the k*sd term scales with the M spread. -function _gp_mean_sd(X, solved) - lis = collect(keys(solved)) - yt = [solved[li] for li in lis] - ybar = sum(yt) / length(yt) - ystd = max(sqrt(sum(abs2, yt .- ybar) / max(length(yt) - 1, 1)), 1e-8) - kern = KernelFunctions.with_lengthscale( - KernelFunctions.SqExponentialKernel(), 0.1) - post = AbstractGPs.posterior( - AbstractGPs.GP(kern)(X[lis], 1e-8), (yt .- ybar) ./ ystd) - mz = AbstractGPs.mean(post, X) - vz = max.(AbstractGPs.var(post, X), 0.0) - return mz .* ystd .+ ybar, sqrt.(vz) .* ystd -end - -# Count of per-support M subproblems solved (for the grid-vs-GP comparison). -const _M_SOLVE_COUNT = Ref(0) - -# Original workflow: solve M at every support (grid). Kept for the mbm_gp -# vs grid comparison, gated by ENV["DP_MBM_GRID"]. -function _grid_M_vals(objectives, inner_sub, method) +# Solve the M subproblem exactly at every support +function DP.sample_M_values( + sampler::Symbol, + objectives::AbstractArray, + sub::DP.GDPSubmodel, + method::DP._MBM, + grids::Tuple + ) + sampler === :exact || error( + "Unrecognized `M_sampler` `$(repr(sampler))` for MBM on an " * + "infinite model. Use `:auto`, `:exact`, or a `GPSampler`.") M_vals = Array{Float64}(undef, size(objectives)) for I in eachindex(objectives) - _M_SOLVE_COUNT[] += 1 - m = DP.raw_M(inner_sub, objectives[I], method) + m = DP.raw_M(sub, objectives[I], method) m === nothing && return nothing M_vals[I] = m end return M_vals end -# Solve M at actively-selected supports (max-UCB acquisition) and fill the -# rest with the UCB (mean + k*sd), a valid over-estimate. Returns a scalar -# when M is uniform (e.g. dependent multi-dim parameters where M does not -# vary), matching the grid workflow's early return. -function _gp_active_M_vals(objectives, inner_sub, method, sub, mini_expr) - idxs = collect(CartesianIndices(objectives)) - n = length(idxs) - k = 2.5 - budget = min(n, max(6, cld(n, 4))) - solved = Dict{Int, Float64}() - solve_at!(li) = begin - _M_SOLVE_COUNT[] += 1 - m = DP.raw_M(inner_sub, objectives[idxs[li]], method) - m === nothing && return false - solved[li] = m - true - end - for s in unique([1, cld(n + 1, 2), n]) - solve_at!(s) || return nothing - end - seed = collect(values(solved)) - all(==(first(seed)), seed) && return first(seed) # uniform M -> scalar - mini_prefs = InfiniteOpt.parameter_refs(mini_expr) - reverse_map = Dict(ws[1] => v for (v, ws) in sub.fwd_map) - prefs = Tuple(reverse_map[p] for p in mini_prefs) - grids = Tuple(InfiniteOpt.supports(p) for p in prefs) - X, _ = _gp_support_coords(grids) - while length(solved) < budget - ms, ss = _gp_mean_sd(X, solved) - acq = ms .+ k .* ss - for li in keys(solved) - acq[li] = -Inf - end - solve_at!(argmax(acq)) || return nothing - end - ms, ss = _gp_mean_sd(X, solved) - M_vals = Array{Float64}(undef, size(objectives)) - for (li, I) in enumerate(idxs) - M_vals[I] = get(solved, li, ms[li] + k * ss[li]) - end - return M_vals -end - -# Transcribe mini_expr, then approximate M(d) with an actively-sampled GP -# instead of solving at every support; aggregate to a scalar if uniform. +# Transcribe mini_expr, compute the per-support M values with the +# resolved M sampler, and aggregate to a scalar if uniform, else to a +# parameter function on main. function DP.raw_M( sub::DP.GDPSubmodel{<:InfiniteOpt.InfiniteModel}, mini_expr::JuMP.AbstractJuMPScalar, method::DP._MBM ) objectives = InfiniteOpt.transformation_expression(mini_expr) + # transcription orders the dimensions by parameter group, which is + # not the ascending order `parameter_refs` gives the grids below + group_idxs = InfiniteOpt.parameter_group_int_indices(mini_expr) + if length(group_idxs) > 1 && ndims(objectives) == length(group_idxs) + objectives = permutedims(objectives, sortperm(group_idxs)) + end transcribed = InfiniteOpt.transformation_model(sub.model) - inner_sub = DP.GDPSubmodel(transcribed,JuMP.VariableRef[], - Dict{JuMP.VariableRef, Vector{JuMP.VariableRef}}() - ) - M_vals = haskey(ENV, "DP_MBM_GRID") ? - _grid_M_vals(objectives, inner_sub, method) : - _gp_active_M_vals(objectives, inner_sub, method, sub, mini_expr) - M_vals === nothing && return nothing - M_vals isa Number && return M_vals - all(==(first(M_vals)), M_vals) && return first(M_vals) + inner_sub = DP.GDPSubmodel(transcribed, JuMP.VariableRef[], + Dict{JuMP.VariableRef, Vector{JuMP.VariableRef}}()) mini_prefs = InfiniteOpt.parameter_refs(mini_expr) reverse_map = Dict(ws[1] => v for (v, ws) in sub.fwd_map) prefs = Tuple(reverse_map[p] for p in mini_prefs) - main = JuMP.owner_model(first(prefs)) grids = Tuple(InfiniteOpt.supports(p) for p in prefs) + sampler = _resolve_M_sampler(method.M_sampler) + M_vals = DP.sample_M_values(sampler, objectives, inner_sub, method, grids) + M_vals === nothing && return nothing + M_vals isa Number && return M_vals + all(==(first(M_vals)), M_vals) && return first(M_vals) + main = JuMP.owner_model(first(prefs)) param_func = InfiniteOpt.build_parameter_function( error, _interpolate(grids, M_vals), prefs) return InfiniteOpt.add_parameter_function(main, param_func) diff --git a/ext/InfiniteGPDisjunctiveProgramming.jl b/ext/InfiniteGPDisjunctiveProgramming.jl new file mode 100644 index 0000000..0e9e136 --- /dev/null +++ b/ext/InfiniteGPDisjunctiveProgramming.jl @@ -0,0 +1,135 @@ +module InfiniteGPDisjunctiveProgramming + +import InfiniteOpt, JuMP +import AbstractGPs, KernelFunctions +import DisjunctiveProgramming as DP + +################################################################################ +# GP SAMPLER +################################################################################ +# See the `GPSampler` docstring in `src/extension_api.jl` +struct GPSampler{K} + kappa::Float64 + budget::Float64 + min_solves::Int + kernel::K +end + +function DP.GPSampler(; + kappa::Real = 2.5, + budget::Real = 0.25, + min_solves::Int = 6, + kernel = nothing + ) + kappa >= 0 || error("`kappa` must be nonnegative.") + 0 < budget <= 1 || error("`budget` must be in `(0, 1]`.") + min_solves >= 1 || error("`min_solves` must be at least 1.") + return GPSampler(Float64(kappa), Float64(budget), min_solves, kernel) +end + +################################################################################ +# GP FITTING +################################################################################ +# Lengthscale candidates on the [0, 1]-normalized support coordinates +const _LENGTHSCALES = (0.05, 0.1, 0.2, 0.4, 0.8) + +# Observation jitter for the GP fit +const _JITTER = 1e-8 + +# Coordinates of every support in linear index order, normalized to +# [0, 1]^d so one isotropic lengthscale works across dimensions +function _support_coords(grids) + idxs = CartesianIndices(length.(grids)) + los = [minimum(g) for g in grids] + rng = [max(maximum(g) - minimum(g), eps()) for g in grids] + return [[(grids[d][I[d]] - los[d]) / rng[d] for d in 1:length(grids)] + for I in vec(idxs)] +end + +# Fit the GP posterior on the solved coordinates; `y` is standardized +# by the caller. With no user kernel, select the lengthscale of a +# squared exponential kernel by maximizing the marginal likelihood. +function _fit_posterior(sampler::GPSampler, X, y) + isnothing(sampler.kernel) || return AbstractGPs.posterior( + AbstractGPs.GP(sampler.kernel)(X, _JITTER), y) + best_post, best_lp = nothing, -Inf + for ls in _LENGTHSCALES + kern = KernelFunctions.with_lengthscale( + KernelFunctions.SqExponentialKernel(), ls) + fx = AbstractGPs.GP(kern)(X, _JITTER) + lp = AbstractGPs.logpdf(fx, y) + if lp > best_lp + best_post, best_lp = AbstractGPs.posterior(fx, y), lp + end + end + return best_post +end + +# Posterior mean and sd at all coords given the solved (index => M) +# samples, destandardized back to M units +function _mean_sd(sampler::GPSampler, X, solved) + lis = collect(keys(solved)) + y = [solved[li] for li in lis] + ybar = sum(y) / length(y) + ystd = max(sqrt(sum(abs2, y .- ybar) / max(length(y) - 1, 1)), 1e-8) + post = _fit_posterior(sampler, X[lis], (y .- ybar) ./ ystd) + mz = AbstractGPs.mean(post, X) + vz = max.(AbstractGPs.var(post, X), 0.0) + return mz .* ystd .+ ybar, sqrt.(vz) .* ystd +end + +################################################################################ +# M VALUE SAMPLING +################################################################################ +# Solve M at actively-selected supports (max-UCB acquisition) and fill +# the rest with the upper confidence bound mean + kappa * sd, a +# heuristic over-estimate. Returns a scalar when the seed M values are +# uniform (e.g. M does not vary over the supports). +function DP.sample_M_values( + sampler::GPSampler, + objectives::AbstractArray, + sub::DP.GDPSubmodel, + method::DP._MBM, + grids::Tuple + ) + idxs = collect(CartesianIndices(objectives)) + n = length(idxs) + solved = Dict{Int, Float64}() + solve_at(li) = begin + m = DP.raw_M(sub, objectives[idxs[li]], method) + m === nothing && return false + solved[li] = m + return true + end + for s in unique([1, cld(n + 1, 2), n]) + solve_at(s) || return nothing + end + seed = collect(values(solved)) + all(==(first(seed)), seed) && return first(seed) + budget = clamp( + ceil(Int, sampler.budget * n), min(sampler.min_solves, n), n) + X = _support_coords(grids) + while length(solved) < budget + ms, ss = _mean_sd(sampler, X, solved) + acq = ms .+ sampler.kappa .* ss + for li in keys(solved) + acq[li] = -Inf + end + solve_at(argmax(acq)) || return nothing + end + M_vals = Array{Float64}(undef, size(objectives)) + if length(solved) == n # nothing left to estimate + for (li, I) in enumerate(idxs) + M_vals[I] = solved[li] + end + return M_vals + end + ms, ss = _mean_sd(sampler, X, solved) + for (li, I) in enumerate(idxs) + # exact M values are nonnegative, so the fill is too + M_vals[I] = get(solved, li, max(ms[li] + sampler.kappa * ss[li], 0.0)) + end + return M_vals +end + +end diff --git a/src/datatypes.jl b/src/datatypes.jl index bdcd4c2..2bb9038 100644 --- a/src/datatypes.jl +++ b/src/datatypes.jl @@ -368,21 +368,28 @@ struct BigM{T} <: AbstractReformulationMethod end """ - MBM{O, T, L <: LogicalVariableRef} <: AbstractReformulationMethod + MBM{O, T} <: AbstractReformulationMethod A type for using the multiple big-M reformulation approach for disjunctive constraints. **Fields** - `optimizer::O`: Optimizer to use when solving mini-models (required). - `default_M::T`: Default big-M value to use if no big-M is specified for a logical variable (1e9). +- `M_sampler::Any`: Strategy for computing M values across the supports + of an infinite model (`:auto`). `:auto` uses [`GPSampler`](@ref) when + the AbstractGPs extension is loaded and `:exact` otherwise; `:exact` + solves an M subproblem at every support. Ignored for finite models. """ mutable struct MBM{O, T} <: AbstractReformulationMethod optimizer::O default_M::T - + M_sampler::Any + # Constructor with optimizer (required) and optional default_M - function MBM(optimizer::O, default_M::T = 1e9) where {O, T} - new{O, T}(optimizer, default_M) + function MBM( + optimizer::O, default_M::T = 1e9; M_sampler = :auto + ) where {O, T} + new{O, T}(optimizer, default_M, M_sampler) end end @@ -390,6 +397,7 @@ mutable struct _MBM{O, T, M <: JuMP.AbstractModel} <: AbstractReformulationMetho optimizer::O M::Dict{LogicalVariableRef{M}, Any} default_M::T + M_sampler::Any subproblem_indicators::Vector{LogicalVariableRef{M}} # Cached submodels: indicator => GDPSubmodel. # Typed Any so extensions can store different types. @@ -400,6 +408,7 @@ mutable struct _MBM{O, T, M <: JuMP.AbstractModel} <: AbstractReformulationMetho method.optimizer, Dict{LogicalVariableRef{M}, Any}(), method.default_M, + method.M_sampler, Vector{LogicalVariableRef{M}}(), Dict{LogicalVariableRef{M}, Any}() ) diff --git a/src/extension_api.jl b/src/extension_api.jl index ad3566f..6bfa129 100644 --- a/src/extension_api.jl +++ b/src/extension_api.jl @@ -38,3 +38,53 @@ Y(t, x) ``` """ function InfiniteLogical end + +""" + GPSampler(; kappa = 2.5, budget = 0.25, min_solves = 6, + kernel = nothing) + +Creates a Gaussian-process M sampler for [`MBM`](@ref) on infinite +models. Instead of solving an M subproblem at every support of the +infinite parameters, the sampler solves a subset of the supports +selected by an upper-confidence-bound acquisition and fills the +remaining supports with the posterior upper confidence bound +`mean + kappa * sd`. The filled values are heuristic upper estimates +of the exact M values, not certificates; increase `kappa` for more +conservative estimates or use `MBM(...; M_sampler = :exact)` to solve +every support. This requires that InfiniteOpt and AbstractGPs be +imported first, in which case it is also the default M sampler (see +the `M_sampler` field of [`MBM`](@ref)). + +**Keyword Arguments** +- `kappa::Real`: Upper-confidence-bound multiplier used to select the + next support to solve and to fill unsolved supports (2.5). +- `budget::Real`: Fraction of the supports to solve exactly, in + `(0, 1]` (0.25). +- `min_solves::Int`: Minimum number of exactly solved supports (6). +- `kernel`: Covariance kernel for the GP fit. Defaults to a squared + exponential kernel whose lengthscale is selected by maximizing the + marginal likelihood; pass a `KernelFunctions` kernel to override. + +**Example** +```julia +julia> using DisjunctiveProgramming, InfiniteOpt, AbstractGPs, HiGHS + +julia> method = MBM(HiGHS.Optimizer, M_sampler = GPSampler(kappa = 4.0)) +``` +""" +function GPSampler end + +""" + sample_M_values(sampler, objectives, sub, method, grids) + +Compute the MBM M values at the transcription supports of an infinite +model. `objectives` is the array of per-support objective expressions, +`sub` is the transcribed submodel wrapped as a `GDPSubmodel`, `method` +is the `_MBM` data, and `grids` are the support vectors of the +infinite parameters. Returns an array of M values shaped like +`objectives`, a scalar when M is uniform across the supports, or +`nothing` if an M subproblem is infeasible. Extensions implement +methods for their sampler types (e.g. [`GPSampler`](@ref)); the +`:exact` sampler solves every support. +""" +function sample_M_values end diff --git a/src/mbm.jl b/src/mbm.jl index 23c5096..4441d95 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -87,7 +87,8 @@ function reformulate_disjunct_constraint( }, method::_MBM ) - ref_cons = reformulate_disjunction(model, con, MBM(method.optimizer)) + ref_cons = reformulate_disjunction(model, con, MBM( + method.optimizer, method.default_M, M_sampler = method.M_sampler)) new_ref_cons = Vector{JuMP.AbstractConstraint}() for ref_con in ref_cons append!(new_ref_cons, diff --git a/test/extensions/InfiniteDisjunctiveProgramming.jl b/test/extensions/InfiniteDisjunctiveProgramming.jl index 47bb884..a2daf8e 100644 --- a/test/extensions/InfiniteDisjunctiveProgramming.jl +++ b/test/extensions/InfiniteDisjunctiveProgramming.jl @@ -376,7 +376,7 @@ function test_raw_M_infinite_scalar() @constraint(model, con, x >= 5, Disjunct(Y[1])) @constraint(model, con2, x <= 3, Disjunct(Y[2])) @disjunction(model, Y) - mbm = DP._MBM(MBM(HiGHS.Optimizer), model) + mbm = DP._MBM(MBM(HiGHS.Optimizer, M_sampler = :exact), model) sub = DP.copy_model_with_constraints( model, DP.DisjunctConstraintRef[con2], mbm) obj = DP.prepare_max_M_objective( @@ -399,7 +399,7 @@ function test_raw_M_infinite_param_function() @constraint(model, con, x <= f, Disjunct(Y[1])) @constraint(model, con2, x >= 0.5, Disjunct(Y[2])) @disjunction(model, Y) - mbm = DP._MBM(MBM(HiGHS.Optimizer), model) + mbm = DP._MBM(MBM(HiGHS.Optimizer, M_sampler = :exact), model) sub = DP.copy_model_with_constraints( model, DP.DisjunctConstraintRef[con2], mbm) obj = DP.prepare_max_M_objective( @@ -413,6 +413,34 @@ function test_raw_M_infinite_param_function() end end +# raw_M over two infinite parameters with different support counts. +# Transcription orders the objective dimensions by parameter group, +# which need not be the ascending order of the grids, so the M values +# must be permuted to line up. Setup: x(t, s) in [0, 10], +# disj1: x <= t + s, disj2: x >= 0.5. Slack r(x) = x - t - s +# maximized over x in [0.5, 10]: 10 - t - s. +function test_raw_M_infinite_two_params() + model = InfiniteGDPModel() + @infinite_parameter(model, t ∈ [0, 1], supports = [0.0, 0.5, 1.0]) + @infinite_parameter(model, s ∈ [0, 1], supports = [0.0, 1.0]) + @variable(model, 0 <= x <= 10, Infinite(t, s)) + @variable(model, Y[1:2], InfiniteLogical(t, s)) + @constraint(model, con, x <= t + s, Disjunct(Y[1])) + @constraint(model, con2, x >= 0.5, Disjunct(Y[2])) + @disjunction(model, Y) + mbm = DP._MBM(MBM(HiGHS.Optimizer, M_sampler = :exact), model) + sub = DP.copy_model_with_constraints( + model, DP.DisjunctConstraintRef[con2], mbm) + obj = DP.prepare_max_M_objective( + model, JuMP.constraint_object(con), sub) + M = DP.raw_M(sub, obj, mbm) + @test M isa InfiniteOpt.GeneralVariableRef + raw_fn = InfiniteOpt.raw_function(M) + for t_val in [0.0, 0.5, 1.0], s_val in [0.0, 1.0] + @test raw_fn(t_val, s_val) >= 10.0 - t_val - s_val - 1e-6 + end +end + # Piecewise-constant max-of-corners: returns the maximum value over # the 2^n corners of the cell containing the query. function test_interpolate() @@ -834,6 +862,7 @@ end test_interpolate() test_raw_M_infinite_scalar() test_raw_M_infinite_param_function() + test_raw_M_infinite_two_params() test_mbm_finite_and_integer_var() test_mbm_infinite_simple() test_mbm_infinite_param_dependent() diff --git a/test/extensions/InfiniteGPDisjunctiveProgramming.jl b/test/extensions/InfiniteGPDisjunctiveProgramming.jl new file mode 100644 index 0000000..8b702df --- /dev/null +++ b/test/extensions/InfiniteGPDisjunctiveProgramming.jl @@ -0,0 +1,163 @@ +using InfiniteOpt, HiGHS, AbstractGPs, KernelFunctions +import DisjunctiveProgramming as DP + +# Helpers to access internal functions of the two extensions +const IGDP = Base.get_extension(DP, :InfiniteGPDisjunctiveProgramming) +const IDP = Base.get_extension(DP, :InfiniteDisjunctiveProgramming) + +function test_gp_sampler_creation() + sampler = GPSampler() + @test sampler isa IGDP.GPSampler + @test sampler.kappa == 2.5 + @test sampler.budget == 0.25 + @test sampler.min_solves == 6 + @test isnothing(sampler.kernel) + kern = with_lengthscale(SqExponentialKernel(), 0.3) + sampler = GPSampler( + kappa = 4.0, budget = 0.1, min_solves = 3, kernel = kern) + @test sampler.kappa == 4.0 + @test sampler.budget == 0.1 + @test sampler.min_solves == 3 + @test sampler.kernel === kern + @test_throws ErrorException GPSampler(kappa = -1) + @test_throws ErrorException GPSampler(budget = 0) + @test_throws ErrorException GPSampler(budget = 1.5) + @test_throws ErrorException GPSampler(min_solves = 0) +end + +function test_gp_sampler_resolution() + # the GP extension is loaded, so :auto resolves to a GPSampler + @test IDP._resolve_M_sampler(:auto) isa IGDP.GPSampler + @test IDP._resolve_M_sampler(:exact) === :exact + sampler = GPSampler(kappa = 3.0) + @test IDP._resolve_M_sampler(sampler) === sampler + @test MBM(HiGHS.Optimizer).M_sampler === :auto + @test MBM(HiGHS.Optimizer, M_sampler = :exact).M_sampler === :exact +end + +# Mirror of test_raw_M_infinite_scalar: uniform seed M values collapse +# to the exactly-solved scalar under the GP sampler +function test_gp_raw_M_scalar() + model = InfiniteGDPModel() + @infinite_parameter(model, t ∈ [0, 1], num_supports = 5) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, Y[1:2], InfiniteLogical(t)) + @constraint(model, con, x >= 5, Disjunct(Y[1])) + @constraint(model, con2, x <= 3, Disjunct(Y[2])) + @disjunction(model, Y) + mbm = DP._MBM(MBM(HiGHS.Optimizer), model) + sub = DP.copy_model_with_constraints( + model, DP.DisjunctConstraintRef[con2], mbm) + obj = DP.prepare_max_M_objective( + model, JuMP.constraint_object(con), sub) + @test DP.raw_M(sub, obj, mbm) == 5.0 +end + +# With few supports the budget floor covers every support, so the GP +# sampler solves all of them exactly and must reproduce the exact +# grid parameter function +function test_gp_raw_M_matches_exact() + function pfunc_values(M_sampler, supports) + model = InfiniteGDPModel() + @infinite_parameter(model, t ∈ [0, 1], supports = supports) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, Y[1:2], InfiniteLogical(t)) + @parameter_function(model, f == t -> 2*t) + @constraint(model, con, x <= f, Disjunct(Y[1])) + @constraint(model, con2, x >= 0.5, Disjunct(Y[2])) + @disjunction(model, Y) + mbm = DP._MBM( + MBM(HiGHS.Optimizer, M_sampler = M_sampler), model) + sub = DP.copy_model_with_constraints( + model, DP.DisjunctConstraintRef[con2], mbm) + obj = DP.prepare_max_M_objective( + model, JuMP.constraint_object(con), sub) + M = DP.raw_M(sub, obj, mbm) + @test M isa InfiniteOpt.GeneralVariableRef + return [InfiniteOpt.raw_function(M)(t_val) for t_val in supports] + end + supports = [0.0, 0.25, 0.5, 0.75, 1.0] + exact_vals = pfunc_values(:exact, supports) + @test pfunc_values(GPSampler(), supports) == exact_vals + # a user kernel skips the lengthscale fit but solves the same supports + kern = with_lengthscale(SqExponentialKernel(), 0.2) + @test pfunc_values(GPSampler(kernel = kern), supports) == exact_vals +end + +# an empty disjunct region makes the M subproblems infeasible; both +# samplers propagate that up to the reformulation error +function test_gp_infeasible_disjunct() + function build() + model = InfiniteGDPModel(HiGHS.Optimizer) + set_silent(model) + @infinite_parameter(model, t ∈ [0, 1], num_supports = 5) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, Y[1:2], InfiniteLogical(t)) + @parameter_function(model, f == t -> 2*t) + @constraint(model, x <= f, Disjunct(Y[1])) + @constraint(model, x >= 8, Disjunct(Y[2])) + @constraint(model, x <= 3, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, 𝔼(x, t)) + return model + end + for sampler in (:exact, GPSampler()) + model = build() + @test_throws ErrorException optimize!(model, + gdp_method = MBM(HiGHS.Optimizer, M_sampler = sampler)) + end +end + +# optimum (10) needs M(t) >= 10 - 2t pointwise; the GP fill is heuristic +function test_gp_mbm_solve_equivalence() + function solve_with(M_sampler) + model = InfiniteGDPModel(HiGHS.Optimizer) + set_silent(model) + @infinite_parameter(model, t ∈ [0, 1], num_supports = 20) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, Y[1:2], InfiniteLogical(t)) + @parameter_function(model, f == t -> 2*t) + @constraint(model, x <= f, Disjunct(Y[1])) + @constraint(model, x >= 0.5, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, 𝔼(x, t)) + optimize!(model, + gdp_method = MBM(HiGHS.Optimizer, M_sampler = M_sampler)) + @test termination_status(model) == MOI.OPTIMAL + return objective_value(model) + end + obj_exact = solve_with(:exact) + obj_auto = solve_with(:auto) + obj_gp = solve_with(GPSampler(kappa = 4.0, budget = 0.2)) + @test obj_exact ≈ 10.0 atol = 1e-4 + # over-M can't raise the optimum, under-M can only shave it a bit + @test obj_auto <= obj_exact + 1e-6 + @test obj_auto ≈ obj_exact atol = 1e-2 + @test obj_gp <= obj_exact + 1e-6 + @test obj_gp ≈ obj_exact atol = 1e-2 +end + +function test_gp_unknown_sampler_error() + model = InfiniteGDPModel(HiGHS.Optimizer) + set_silent(model) + @infinite_parameter(model, t ∈ [0, 1], num_supports = 5) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, Y[1:2], InfiniteLogical(t)) + @parameter_function(model, f == t -> 2*t) + @constraint(model, x <= f, Disjunct(Y[1])) + @constraint(model, x >= 0.5, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, 𝔼(x, t)) + @test_throws ErrorException optimize!(model, + gdp_method = MBM(HiGHS.Optimizer, M_sampler = :grid)) +end + +@testset "InfiniteGPDisjunctiveProgramming" begin + test_gp_sampler_creation() + test_gp_sampler_resolution() + test_gp_raw_M_scalar() + test_gp_raw_M_matches_exact() + test_gp_mbm_solve_equivalence() + test_gp_unknown_sampler_error() + test_gp_infeasible_disjunct() +end diff --git a/test/runtests.jl b/test/runtests.jl index 06e8813..569a12b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -24,4 +24,5 @@ include("constraints/disjunction.jl") include("print.jl") include("solve.jl") include("extensions/InfiniteDisjunctiveProgramming.jl") +include("extensions/InfiniteGPDisjunctiveProgramming.jl") From fa4a7ca7530537442e8d56802986b8f5102701a5 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Sat, 1 Aug 2026 23:39:33 -0400 Subject: [PATCH 03/10] Don't compute grids until we need them --- ext/InfiniteDisjunctiveProgramming.jl | 20 ++++++++++++----- ext/InfiniteGPDisjunctiveProgramming.jl | 4 ++-- src/extension_api.jl | 9 +++++--- .../InfiniteDisjunctiveProgramming.jl | 22 +++++++++++++++++++ 4 files changed, 44 insertions(+), 11 deletions(-) diff --git a/ext/InfiniteDisjunctiveProgramming.jl b/ext/InfiniteDisjunctiveProgramming.jl index 2055422..37eb8ff 100644 --- a/ext/InfiniteDisjunctiveProgramming.jl +++ b/ext/InfiniteDisjunctiveProgramming.jl @@ -263,6 +263,16 @@ function _interpolate_at( ) end +# The infinite parameters of `mini_expr` and their supports, in the +# ascending order of `parameter_refs`. Only defined when M varies over +# the supports, so it is deferred until an M sampler needs it. +function _support_grids(sub::DP.GDPSubmodel, mini_expr) + reverse_map = Dict(ws[1] => v for (v, ws) in sub.fwd_map) + prefs = Tuple(reverse_map[p] + for p in InfiniteOpt.parameter_refs(mini_expr)) + return prefs, Tuple(InfiniteOpt.supports(p) for p in prefs) +end + # Resolve `:auto` to the GP sampler when the GP extension is loaded function _resolve_M_sampler(sampler) sampler === :auto || return sampler @@ -276,7 +286,7 @@ function DP.sample_M_values( objectives::AbstractArray, sub::DP.GDPSubmodel, method::DP._MBM, - grids::Tuple + support_grids ) sampler === :exact || error( "Unrecognized `M_sampler` `$(repr(sampler))` for MBM on an " * @@ -308,15 +318,13 @@ function DP.raw_M( transcribed = InfiniteOpt.transformation_model(sub.model) inner_sub = DP.GDPSubmodel(transcribed, JuMP.VariableRef[], Dict{JuMP.VariableRef, Vector{JuMP.VariableRef}}()) - mini_prefs = InfiniteOpt.parameter_refs(mini_expr) - reverse_map = Dict(ws[1] => v for (v, ws) in sub.fwd_map) - prefs = Tuple(reverse_map[p] for p in mini_prefs) - grids = Tuple(InfiniteOpt.supports(p) for p in prefs) sampler = _resolve_M_sampler(method.M_sampler) - M_vals = DP.sample_M_values(sampler, objectives, inner_sub, method, grids) + M_vals = DP.sample_M_values(sampler, objectives, inner_sub, method, + () -> _support_grids(sub, mini_expr)[2]) M_vals === nothing && return nothing M_vals isa Number && return M_vals all(==(first(M_vals)), M_vals) && return first(M_vals) + prefs, grids = _support_grids(sub, mini_expr) main = JuMP.owner_model(first(prefs)) param_func = InfiniteOpt.build_parameter_function( error, _interpolate(grids, M_vals), prefs) diff --git a/ext/InfiniteGPDisjunctiveProgramming.jl b/ext/InfiniteGPDisjunctiveProgramming.jl index 0e9e136..4e93519 100644 --- a/ext/InfiniteGPDisjunctiveProgramming.jl +++ b/ext/InfiniteGPDisjunctiveProgramming.jl @@ -90,7 +90,7 @@ function DP.sample_M_values( objectives::AbstractArray, sub::DP.GDPSubmodel, method::DP._MBM, - grids::Tuple + support_grids ) idxs = collect(CartesianIndices(objectives)) n = length(idxs) @@ -108,7 +108,7 @@ function DP.sample_M_values( all(==(first(seed)), seed) && return first(seed) budget = clamp( ceil(Int, sampler.budget * n), min(sampler.min_solves, n), n) - X = _support_coords(grids) + X = _support_coords(support_grids()) while length(solved) < budget ms, ss = _mean_sd(sampler, X, solved) acq = ms .+ sampler.kappa .* ss diff --git a/src/extension_api.jl b/src/extension_api.jl index 6bfa129..38273ba 100644 --- a/src/extension_api.jl +++ b/src/extension_api.jl @@ -75,13 +75,16 @@ julia> method = MBM(HiGHS.Optimizer, M_sampler = GPSampler(kappa = 4.0)) function GPSampler end """ - sample_M_values(sampler, objectives, sub, method, grids) + sample_M_values(sampler, objectives, sub, method, support_grids) Compute the MBM M values at the transcription supports of an infinite model. `objectives` is the array of per-support objective expressions, `sub` is the transcribed submodel wrapped as a `GDPSubmodel`, `method` -is the `_MBM` data, and `grids` are the support vectors of the -infinite parameters. Returns an array of M values shaped like +is the `_MBM` data, and `support_grids` is a function returning the +support vectors of the infinite parameters. It is a function because +the supports are only well defined once M is known to vary over them, +so samplers that return early (or never need coordinates) must not +call it. Returns an array of M values shaped like `objectives`, a scalar when M is uniform across the supports, or `nothing` if an M subproblem is infeasible. Extensions implement methods for their sampler types (e.g. [`GPSampler`](@ref)); the diff --git a/test/extensions/InfiniteDisjunctiveProgramming.jl b/test/extensions/InfiniteDisjunctiveProgramming.jl index a2daf8e..3526d62 100644 --- a/test/extensions/InfiniteDisjunctiveProgramming.jl +++ b/test/extensions/InfiniteDisjunctiveProgramming.jl @@ -441,6 +441,27 @@ function test_raw_M_infinite_two_params() end end +# Dependent parameters have no per-parameter support grid, so raw_M +# must not need one when M does not vary. Setup as in +# test_raw_M_infinite_scalar, over a dependent parameter array. +function test_raw_M_infinite_dependent_params() + model = InfiniteGDPModel() + @infinite_parameter(model, ξ[1:2] ∈ [0, 1], num_supports = 4) + @variable(model, 0 <= x <= 10, Infinite(ξ)) + @variable(model, Y[1:2], InfiniteLogical(ξ)) + @constraint(model, con, x >= 5, Disjunct(Y[1])) + @constraint(model, con2, x <= 3, Disjunct(Y[2])) + @disjunction(model, Y) + for sampler in (:exact, :auto) + mbm = DP._MBM(MBM(HiGHS.Optimizer, M_sampler = sampler), model) + sub = DP.copy_model_with_constraints( + model, DP.DisjunctConstraintRef[con2], mbm) + obj = DP.prepare_max_M_objective( + model, JuMP.constraint_object(con), sub) + @test DP.raw_M(sub, obj, mbm) == 5.0 + end +end + # Piecewise-constant max-of-corners: returns the maximum value over # the 2^n corners of the cell containing the query. function test_interpolate() @@ -863,6 +884,7 @@ end test_raw_M_infinite_scalar() test_raw_M_infinite_param_function() test_raw_M_infinite_two_params() + test_raw_M_infinite_dependent_params() test_mbm_finite_and_integer_var() test_mbm_infinite_simple() test_mbm_infinite_param_dependent() From 5a010118e6327aebf71a209d5c0ebea20f05ae09 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Sun, 2 Aug 2026 13:26:48 -0400 Subject: [PATCH 04/10] detect_uniform_M arg --- ext/InfiniteDisjunctiveProgramming.jl | 7 +- ext/InfiniteGPDisjunctiveProgramming.jl | 42 +++++----- src/extension_api.jl | 10 ++- .../InfiniteGPDisjunctiveProgramming.jl | 82 ++++++++++++++++++- 4 files changed, 114 insertions(+), 27 deletions(-) diff --git a/ext/InfiniteDisjunctiveProgramming.jl b/ext/InfiniteDisjunctiveProgramming.jl index 37eb8ff..38d0a82 100644 --- a/ext/InfiniteDisjunctiveProgramming.jl +++ b/ext/InfiniteDisjunctiveProgramming.jl @@ -268,8 +268,11 @@ end # the supports, so it is deferred until an M sampler needs it. function _support_grids(sub::DP.GDPSubmodel, mini_expr) reverse_map = Dict(ws[1] => v for (v, ws) in sub.fwd_map) - prefs = Tuple(reverse_map[p] - for p in InfiniteOpt.parameter_refs(mini_expr)) + prefs = Tuple(get(reverse_map, p) do + error("MBM cannot build a support grid over `$p`, which " * + "is a group of dependent infinite parameters, so M " * + "must not vary over its supports.") + end for p in InfiniteOpt.parameter_refs(mini_expr)) return prefs, Tuple(InfiniteOpt.supports(p) for p in prefs) end diff --git a/ext/InfiniteGPDisjunctiveProgramming.jl b/ext/InfiniteGPDisjunctiveProgramming.jl index 4e93519..887f837 100644 --- a/ext/InfiniteGPDisjunctiveProgramming.jl +++ b/ext/InfiniteGPDisjunctiveProgramming.jl @@ -13,31 +13,30 @@ struct GPSampler{K} budget::Float64 min_solves::Int kernel::K + detect_uniform_M::Bool end function DP.GPSampler(; kappa::Real = 2.5, budget::Real = 0.25, min_solves::Int = 6, - kernel = nothing + kernel = nothing, + detect_uniform_M::Bool = true ) kappa >= 0 || error("`kappa` must be nonnegative.") 0 < budget <= 1 || error("`budget` must be in `(0, 1]`.") min_solves >= 1 || error("`min_solves` must be at least 1.") - return GPSampler(Float64(kappa), Float64(budget), min_solves, kernel) + return GPSampler(Float64(kappa), Float64(budget), min_solves, + kernel, detect_uniform_M) end ################################################################################ # GP FITTING ################################################################################ -# Lengthscale candidates on the [0, 1]-normalized support coordinates const _LENGTHSCALES = (0.05, 0.1, 0.2, 0.4, 0.8) - -# Observation jitter for the GP fit const _JITTER = 1e-8 -# Coordinates of every support in linear index order, normalized to -# [0, 1]^d so one isotropic lengthscale works across dimensions +# Normalized to [0, 1]^d so one lengthscale works across dimensions function _support_coords(grids) idxs = CartesianIndices(length.(grids)) los = [minimum(g) for g in grids] @@ -46,9 +45,7 @@ function _support_coords(grids) for I in vec(idxs)] end -# Fit the GP posterior on the solved coordinates; `y` is standardized -# by the caller. With no user kernel, select the lengthscale of a -# squared exponential kernel by maximizing the marginal likelihood. +# Lengthscale by marginal likelihood unless the user gave a kernel function _fit_posterior(sampler::GPSampler, X, y) isnothing(sampler.kernel) || return AbstractGPs.posterior( AbstractGPs.GP(sampler.kernel)(X, _JITTER), y) @@ -65,13 +62,13 @@ function _fit_posterior(sampler::GPSampler, X, y) return best_post end -# Posterior mean and sd at all coords given the solved (index => M) -# samples, destandardized back to M units function _mean_sd(sampler::GPSampler, X, solved) lis = collect(keys(solved)) y = [solved[li] for li in lis] ybar = sum(y) / length(y) - ystd = max(sqrt(sum(abs2, y .- ybar) / max(length(y) - 1, 1)), 1e-8) + # floored so near-equal solved values still cushion the filled ones + ystd = max(sqrt(sum(abs2, y .- ybar) / max(length(y) - 1, 1)), + 1e-2 * abs(ybar), 1e-8) post = _fit_posterior(sampler, X[lis], (y .- ybar) ./ ystd) mz = AbstractGPs.mean(post, X) vz = max.(AbstractGPs.var(post, X), 0.0) @@ -81,10 +78,7 @@ end ################################################################################ # M VALUE SAMPLING ################################################################################ -# Solve M at actively-selected supports (max-UCB acquisition) and fill -# the rest with the upper confidence bound mean + kappa * sd, a -# heuristic over-estimate. Returns a scalar when the seed M values are -# uniform (e.g. M does not vary over the supports). +# Solve M at max-UCB selected supports, fill the rest with the bound function DP.sample_M_values( sampler::GPSampler, objectives::AbstractArray, @@ -101,11 +95,16 @@ function DP.sample_M_values( solved[li] = m return true end - for s in unique([1, cld(n + 1, 2), n]) + # golden fractions; even spacing aliases with a periodic M + for s in unique([1, n, 1 + floor(Int, 0.618 * (n - 1)), + 1 + floor(Int, 0.382 * (n - 1))]) solve_at(s) || return nothing end - seed = collect(values(solved)) - all(==(first(seed)), seed) && return first(seed) + if sampler.detect_uniform_M + # a uniform M needs no fit, and so no support grid either + probes = collect(values(solved)) + all(==(first(probes)), probes) && return first(probes) + end budget = clamp( ceil(Int, sampler.budget * n), min(sampler.min_solves, n), n) X = _support_coords(support_grids()) @@ -125,8 +124,7 @@ function DP.sample_M_values( return M_vals end ms, ss = _mean_sd(sampler, X, solved) - for (li, I) in enumerate(idxs) - # exact M values are nonnegative, so the fill is too + for (li, I) in enumerate(idxs) # exact M values are nonnegative M_vals[I] = get(solved, li, max(ms[li] + sampler.kappa * ss[li], 0.0)) end return M_vals diff --git a/src/extension_api.jl b/src/extension_api.jl index 38273ba..7c220b2 100644 --- a/src/extension_api.jl +++ b/src/extension_api.jl @@ -41,7 +41,7 @@ function InfiniteLogical end """ GPSampler(; kappa = 2.5, budget = 0.25, min_solves = 6, - kernel = nothing) + kernel = nothing, detect_uniform_M = true) Creates a Gaussian-process M sampler for [`MBM`](@ref) on infinite models. Instead of solving an M subproblem at every support of the @@ -64,6 +64,14 @@ the `M_sampler` field of [`MBM`](@ref)). - `kernel`: Covariance kernel for the GP fit. Defaults to a squared exponential kernel whose lengthscale is selected by maximizing the marginal likelihood; pass a `KernelFunctions` kernel to override. +- `detect_uniform_M::Bool`: If `true` (the default), M values that + agree at the first few supports are taken to be uniform and used + for every support. This is cheap and is what makes infinite + parameters without a support grid (e.g. dependent ones) workable, + but it assumes M does not vary elsewhere. Set it to `false` to + always fit the GP, which leaves the usual `kappa * sd` cushion on + the unsolved supports at the cost of the extra solves, and which + requires that every infinite parameter have a support grid. **Example** ```julia diff --git a/test/extensions/InfiniteGPDisjunctiveProgramming.jl b/test/extensions/InfiniteGPDisjunctiveProgramming.jl index 8b702df..796507b 100644 --- a/test/extensions/InfiniteGPDisjunctiveProgramming.jl +++ b/test/extensions/InfiniteGPDisjunctiveProgramming.jl @@ -12,13 +12,15 @@ function test_gp_sampler_creation() @test sampler.budget == 0.25 @test sampler.min_solves == 6 @test isnothing(sampler.kernel) + @test sampler.detect_uniform_M kern = with_lengthscale(SqExponentialKernel(), 0.3) - sampler = GPSampler( - kappa = 4.0, budget = 0.1, min_solves = 3, kernel = kern) + sampler = GPSampler(kappa = 4.0, budget = 0.1, min_solves = 3, + kernel = kern, detect_uniform_M = false) @test sampler.kappa == 4.0 @test sampler.budget == 0.1 @test sampler.min_solves == 3 @test sampler.kernel === kern + @test !sampler.detect_uniform_M @test_throws ErrorException GPSampler(kappa = -1) @test_throws ErrorException GPSampler(budget = 0) @test_throws ErrorException GPSampler(budget = 1.5) @@ -137,6 +139,79 @@ function test_gp_mbm_solve_equivalence() @test obj_gp ≈ obj_exact atol = 1e-2 end +# A periodic M must not read as uniform. With f(t) = 2|cos(2*pi*t)| +# on these supports, M = 10 - f is 8 at supports 1, 3, 5 and 10 at +# supports 2, 4, so evenly spaced probes alias and collapse M to 8, +# which caps x at 8 and cuts the optimum from 10 down to 9. +function test_gp_periodic_M_not_uniform() + supports = [0.0, 0.25, 0.5, 0.75, 1.0] + function solve_with(M_sampler) + model = InfiniteGDPModel(HiGHS.Optimizer) + set_silent(model) + @infinite_parameter(model, t ∈ [0, 1], supports = supports) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, Y[1:2], InfiniteLogical(t)) + @parameter_function(model, f == t -> 2 * abs(cos(2 * pi * t))) + @constraint(model, x <= f, Disjunct(Y[1])) + @constraint(model, x >= 0.5, Disjunct(Y[2])) + @disjunction(model, Y) + @objective(model, Max, 𝔼(x, t)) + optimize!(model, gdp_method = MBM( + HiGHS.Optimizer, M_sampler = M_sampler)) + return objective_value(model) + end + @test solve_with(:exact) ≈ 10.0 atol = 1e-6 + @test solve_with(GPSampler()) ≈ 10.0 atol = 1e-6 +end + +# With detection off the uniform M is not collapsed to a scalar: the +# GP is fit and the unsolved supports keep their kappa * sd cushion, +# which must sit above the M that detection would have returned. +function test_gp_detect_uniform_M_off() + function raw_M_with(detect) + model = InfiniteGDPModel() + @infinite_parameter(model, t ∈ [0, 1], num_supports = 20) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, Y[1:2], InfiniteLogical(t)) + @constraint(model, con, x >= 5, Disjunct(Y[1])) + @constraint(model, con2, x <= 3, Disjunct(Y[2])) + @disjunction(model, Y) + mbm = DP._MBM(MBM(HiGHS.Optimizer, + M_sampler = GPSampler(detect_uniform_M = detect)), model) + sub = DP.copy_model_with_constraints( + model, DP.DisjunctConstraintRef[con2], mbm) + obj = DP.prepare_max_M_objective( + model, JuMP.constraint_object(con), sub) + return DP.raw_M(sub, obj, mbm) + end + @test raw_M_with(true) == 5.0 + M = raw_M_with(false) + @test M isa InfiniteOpt.GeneralVariableRef + raw_fn = InfiniteOpt.raw_function(M) + vals = [raw_fn(t) for t in range(0, 1, length = 20)] + @test all(vals .>= 5.0 - 1e-6) + @test maximum(vals) > 5.0 +end + +# Dependent parameters have no support grid, so turning detection off +# leaves the GP with nothing to fit over +function test_gp_detect_uniform_M_off_dependent() + model = InfiniteGDPModel() + @infinite_parameter(model, ξ[1:2] ∈ [0, 1], num_supports = 4) + @variable(model, 0 <= x <= 10, Infinite(ξ)) + @variable(model, Y[1:2], InfiniteLogical(ξ)) + @constraint(model, con, x >= 5, Disjunct(Y[1])) + @constraint(model, con2, x <= 3, Disjunct(Y[2])) + @disjunction(model, Y) + mbm = DP._MBM(MBM(HiGHS.Optimizer, + M_sampler = GPSampler(detect_uniform_M = false)), model) + sub = DP.copy_model_with_constraints( + model, DP.DisjunctConstraintRef[con2], mbm) + obj = DP.prepare_max_M_objective( + model, JuMP.constraint_object(con), sub) + @test_throws ErrorException DP.raw_M(sub, obj, mbm) +end + function test_gp_unknown_sampler_error() model = InfiniteGDPModel(HiGHS.Optimizer) set_silent(model) @@ -158,6 +233,9 @@ end test_gp_raw_M_scalar() test_gp_raw_M_matches_exact() test_gp_mbm_solve_equivalence() + test_gp_periodic_M_not_uniform() + test_gp_detect_uniform_M_off() + test_gp_detect_uniform_M_off_dependent() test_gp_unknown_sampler_error() test_gp_infeasible_disjunct() end From 55e8b3fadd6b92325ef463a60046e1a9622858d7 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Thu, 13 Aug 2026 17:31:55 -0400 Subject: [PATCH 05/10] Bound info for parameter functions in disjunct constraints --- ext/InfiniteDisjunctiveProgramming.jl | 55 +++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) diff --git a/ext/InfiniteDisjunctiveProgramming.jl b/ext/InfiniteDisjunctiveProgramming.jl index 38d0a82..1a79fa0 100644 --- a/ext/InfiniteDisjunctiveProgramming.jl +++ b/ext/InfiniteDisjunctiveProgramming.jl @@ -38,6 +38,61 @@ function DP.requires_disaggregation(vref::InfiniteOpt.GeneralVariableRef) return !_is_parameter(vref) end +# Bound info for parameter refs in disjunct constraints: parameter +# functions report their support extrema, finite parameters a point, +# other parameters (-Inf, Inf). Returns nothing for variables. +function _parameter_bound_info(vref::InfiniteOpt.GeneralVariableRef) + _is_parameter(vref) || return nothing + dvref = InfiniteOpt.dispatch_variable_ref(vref) + if dvref isa InfiniteOpt.ParameterFunctionRef + prefs = InfiniteOpt.parameter_list(dvref) + length(prefs) == 1 || return (-Inf, Inf) + supps = InfiniteOpt.supports(first(prefs)) + isempty(supps) && return (-Inf, Inf) + f = InfiniteOpt.raw_function(dvref) + vals = [f(s) for s in supps] + return (minimum(vals), maximum(vals)) + elseif dvref isa InfiniteOpt.FiniteParameterRef + v = InfiniteOpt.parameter_value(dvref) + return (v, v) + end + return (-Inf, Inf) +end + +function DP.set_variable_bound_info( + vref::InfiniteOpt.GeneralVariableRef, ::DP.BigM) + info = _parameter_bound_info(vref) + info === nothing || return info + lb = JuMP.has_lower_bound(vref) ? JuMP.lower_bound(vref) : -Inf + ub = JuMP.has_upper_bound(vref) ? JuMP.upper_bound(vref) : Inf + return lb, ub +end + +# Hull and PSplit require finite bounds that include 0 +function _clamped_bound_info(vref, info, method_name) + info === nothing || return (min(0, info[1]), max(0, info[2])) + if !JuMP.has_lower_bound(vref) || !JuMP.has_upper_bound(vref) + error("Variable $vref must have both lower and upper " * + "bounds defined when using the $method_name " * + "reformulation.") + end + return (min(0, JuMP.lower_bound(vref)), + max(0, JuMP.upper_bound(vref))) +end + +function DP.set_variable_bound_info( + vref::InfiniteOpt.GeneralVariableRef, ::DP.Hull) + return _clamped_bound_info(vref, + _parameter_bound_info(vref), "Hull") +end + +function DP.set_variable_bound_info( + vref::InfiniteOpt.GeneralVariableRef, + ::Union{DP.PSplit, DP._PSplit}) + return _clamped_bound_info(vref, + _parameter_bound_info(vref), "PSplit") +end + function DP.VariableProperties(vref::InfiniteOpt.GeneralVariableRef) info = DP.get_variable_info(vref) name = JuMP.name(vref) From 6d3909af4bdea02d018ef4902fe4ccd197284960 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Mon, 17 Aug 2026 13:18:28 -0400 Subject: [PATCH 06/10] New API --- Project.toml | 6 +- README.md | 4 +- ext/InfiniteDisjunctiveProgramming.jl | 21 +--- ext/InfiniteGPDisjunctiveProgramming.jl | 60 +++------ src/datatypes.jl | 53 ++++++-- src/extension_api.jl | 77 ++++-------- src/mbm.jl | 5 +- .../InfiniteDisjunctiveProgramming.jl | 20 ++- .../InfiniteGPDisjunctiveProgramming.jl | 117 ++++++++---------- 9 files changed, 159 insertions(+), 204 deletions(-) diff --git a/Project.toml b/Project.toml index ba1a492..4fa94e8 100644 --- a/Project.toml +++ b/Project.toml @@ -10,17 +10,15 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" [weakdeps] AbstractGPs = "99985d1d-32ba-4be9-9821-2ec096f28918" InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57" -KernelFunctions = "ec8451be-7e33-11e9-00cf-bbf324bd1392" [extensions] InfiniteDisjunctiveProgramming = "InfiniteOpt" -InfiniteGPDisjunctiveProgramming = ["InfiniteOpt", "AbstractGPs", "KernelFunctions"] +InfiniteGPDisjunctiveProgramming = ["InfiniteOpt", "AbstractGPs"] [compat] AbstractGPs = "0.5" Aqua = "0.8" JuMP = "1.18" -KernelFunctions = "0.10" Reexport = "1" julia = "1.10" Juniper = "0.9.3" @@ -35,4 +33,4 @@ Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" Juniper = "2ddba703-00a4-53a7-87a5-e8b9971dde84" [targets] -test = ["Aqua", "HiGHS", "Test", "Juniper", "Ipopt", "InfiniteOpt", "AbstractGPs", "KernelFunctions"] +test = ["Aqua", "HiGHS", "Test", "Juniper", "Ipopt", "InfiniteOpt", "AbstractGPs"] diff --git a/README.md b/README.md index c3cd037..bbf80b4 100644 --- a/README.md +++ b/README.md @@ -176,7 +176,7 @@ The following reformulation methods are currently supported: - `optimizer`: Optimizer to use when solving subproblems to determine M values. This is a required value. - `default_M`: Default big-M value to use if no big-M is specified for a logical variable (1e9). - - `M_sampler`: Strategy for computing M values across the supports of an infinite model. Default: `:auto`, which uses a Gaussian-process sampler (`GPSampler`) when AbstractGPs is loaded and otherwise solves an M subproblem at every support (`:exact`). Ignored for finite models. + - `gp`: Gaussian process (or kernel) used to estimate the M values across the supports of an infinite model. Default: `nothing`, which solves an M subproblem at every support. Pass an `AbstractGPs.AbstractGP` or a `KernelFunctions.Kernel` to instead solve a subset of the supports and fill the rest with a conservative GP estimate. Ignored for finite models, as are the tuning keywords `kappa`, `budget`, `min_solves`, and `detect_uniform_M` (see the `MBM` docstring). 5. [P-Split](https://arxiv.org/abs/2202.05198): This method reformulates each disjunct constraint into P constraints, each with a partitioned group defined by the user. This method requires that terms in the constraint be convex additively seperable with respect to each variable. The `PSplit` struct is created with the following required arguments: @@ -224,7 +224,7 @@ optimize!(model, gdp_method = Hull()) value(W) ``` -When the `MBM` reformulation is used on an infinite model, an M subproblem is solved at every support by default. Loading [AbstractGPs.jl](https://github.com/JuliaGaussianProcesses/AbstractGPs.jl) (`using AbstractGPs`) activates an additional extension that instead solves M at a subset of the supports and fills the rest with a conservative Gaussian-process estimate, which can substantially reduce the number of subproblem solves. The filled values are heuristic upper estimates rather than certificates; see the `GPSampler` docstring for the tuning knobs (`kappa`, `budget`) and use `MBM(optimizer, M_sampler = :exact)` to force exact solves. +When the `MBM` reformulation is used on an infinite model, an M subproblem is solved at every support by default. Loading [AbstractGPs.jl](https://github.com/JuliaGaussianProcesses/AbstractGPs.jl) (`using AbstractGPs`) enables an additional extension that instead solves M at a subset of the supports and fills the rest with a conservative Gaussian-process estimate, which can substantially reduce the number of subproblem solves. To opt in, pass a GP or kernel via the `gp` keyword, e.g. `MBM(optimizer, gp = SqExponentialKernel())` (lengthscale selected by marginal likelihood) or `MBM(optimizer, gp = GP(with_lengthscale(Matern52Kernel(), 0.2)))` (used as given). The filled values are heuristic upper estimates rather than certificates; see the `MBM` docstring for the tuning knobs (`kappa`, `budget`, `min_solves`, `detect_uniform_M`). ## Release Notes diff --git a/ext/InfiniteDisjunctiveProgramming.jl b/ext/InfiniteDisjunctiveProgramming.jl index 1a79fa0..be572ec 100644 --- a/ext/InfiniteDisjunctiveProgramming.jl +++ b/ext/InfiniteDisjunctiveProgramming.jl @@ -320,7 +320,7 @@ end # The infinite parameters of `mini_expr` and their supports, in the # ascending order of `parameter_refs`. Only defined when M varies over -# the supports, so it is deferred until an M sampler needs it. +# the supports, so it is deferred until sample_M_values needs it. function _support_grids(sub::DP.GDPSubmodel, mini_expr) reverse_map = Dict(ws[1] => v for (v, ws) in sub.fwd_map) prefs = Tuple(get(reverse_map, p) do @@ -331,24 +331,14 @@ function _support_grids(sub::DP.GDPSubmodel, mini_expr) return prefs, Tuple(InfiniteOpt.supports(p) for p in prefs) end -# Resolve `:auto` to the GP sampler when the GP extension is loaded -function _resolve_M_sampler(sampler) - sampler === :auto || return sampler - gp = Base.get_extension(DP, :InfiniteGPDisjunctiveProgramming) - return isnothing(gp) ? :exact : DP.GPSampler() -end - # Solve the M subproblem exactly at every support function DP.sample_M_values( - sampler::Symbol, + gp::Nothing, objectives::AbstractArray, sub::DP.GDPSubmodel, method::DP._MBM, support_grids ) - sampler === :exact || error( - "Unrecognized `M_sampler` `$(repr(sampler))` for MBM on an " * - "infinite model. Use `:auto`, `:exact`, or a `GPSampler`.") M_vals = Array{Float64}(undef, size(objectives)) for I in eachindex(objectives) m = DP.raw_M(sub, objectives[I], method) @@ -359,7 +349,7 @@ function DP.sample_M_values( end # Transcribe mini_expr, compute the per-support M values with the -# resolved M sampler, and aggregate to a scalar if uniform, else to a +# method's gp, and aggregate to a scalar if uniform, else to a # parameter function on main. function DP.raw_M( sub::DP.GDPSubmodel{<:InfiniteOpt.InfiniteModel}, @@ -376,9 +366,8 @@ function DP.raw_M( transcribed = InfiniteOpt.transformation_model(sub.model) inner_sub = DP.GDPSubmodel(transcribed, JuMP.VariableRef[], Dict{JuMP.VariableRef, Vector{JuMP.VariableRef}}()) - sampler = _resolve_M_sampler(method.M_sampler) - M_vals = DP.sample_M_values(sampler, objectives, inner_sub, method, - () -> _support_grids(sub, mini_expr)[2]) + M_vals = DP.sample_M_values(method.gp, objectives, inner_sub, + method, () -> _support_grids(sub, mini_expr)[2]) M_vals === nothing && return nothing M_vals isa Number && return M_vals all(==(first(M_vals)), M_vals) && return first(M_vals) diff --git a/ext/InfiniteGPDisjunctiveProgramming.jl b/ext/InfiniteGPDisjunctiveProgramming.jl index 887f837..1c46687 100644 --- a/ext/InfiniteGPDisjunctiveProgramming.jl +++ b/ext/InfiniteGPDisjunctiveProgramming.jl @@ -1,35 +1,10 @@ module InfiniteGPDisjunctiveProgramming import InfiniteOpt, JuMP -import AbstractGPs, KernelFunctions +import AbstractGPs +import AbstractGPs.KernelFunctions import DisjunctiveProgramming as DP -################################################################################ -# GP SAMPLER -################################################################################ -# See the `GPSampler` docstring in `src/extension_api.jl` -struct GPSampler{K} - kappa::Float64 - budget::Float64 - min_solves::Int - kernel::K - detect_uniform_M::Bool -end - -function DP.GPSampler(; - kappa::Real = 2.5, - budget::Real = 0.25, - min_solves::Int = 6, - kernel = nothing, - detect_uniform_M::Bool = true - ) - kappa >= 0 || error("`kappa` must be nonnegative.") - 0 < budget <= 1 || error("`budget` must be in `(0, 1]`.") - min_solves >= 1 || error("`min_solves` must be at least 1.") - return GPSampler(Float64(kappa), Float64(budget), min_solves, - kernel, detect_uniform_M) -end - ################################################################################ # GP FITTING ################################################################################ @@ -45,14 +20,15 @@ function _support_coords(grids) for I in vec(idxs)] end -# Lengthscale by marginal likelihood unless the user gave a kernel -function _fit_posterior(sampler::GPSampler, X, y) - isnothing(sampler.kernel) || return AbstractGPs.posterior( - AbstractGPs.GP(sampler.kernel)(X, _JITTER), y) +# A user GP is used as the prior directly; a kernel gets its +# lengthscale selected by marginal likelihood +function _fit_posterior(gp::AbstractGPs.AbstractGP, X, y) + return AbstractGPs.posterior(gp(X, _JITTER), y) +end +function _fit_posterior(kernel::KernelFunctions.Kernel, X, y) best_post, best_lp = nothing, -Inf for ls in _LENGTHSCALES - kern = KernelFunctions.with_lengthscale( - KernelFunctions.SqExponentialKernel(), ls) + kern = KernelFunctions.with_lengthscale(kernel, ls) fx = AbstractGPs.GP(kern)(X, _JITTER) lp = AbstractGPs.logpdf(fx, y) if lp > best_lp @@ -62,14 +38,14 @@ function _fit_posterior(sampler::GPSampler, X, y) return best_post end -function _mean_sd(sampler::GPSampler, X, solved) +function _mean_sd(gp, X, solved) lis = collect(keys(solved)) y = [solved[li] for li in lis] ybar = sum(y) / length(y) # floored so near-equal solved values still cushion the filled ones ystd = max(sqrt(sum(abs2, y .- ybar) / max(length(y) - 1, 1)), 1e-2 * abs(ybar), 1e-8) - post = _fit_posterior(sampler, X[lis], (y .- ybar) ./ ystd) + post = _fit_posterior(gp, X[lis], (y .- ybar) ./ ystd) mz = AbstractGPs.mean(post, X) vz = max.(AbstractGPs.var(post, X), 0.0) return mz .* ystd .+ ybar, sqrt.(vz) .* ystd @@ -80,7 +56,7 @@ end ################################################################################ # Solve M at max-UCB selected supports, fill the rest with the bound function DP.sample_M_values( - sampler::GPSampler, + gp::Union{AbstractGPs.AbstractGP, KernelFunctions.Kernel}, objectives::AbstractArray, sub::DP.GDPSubmodel, method::DP._MBM, @@ -100,17 +76,17 @@ function DP.sample_M_values( 1 + floor(Int, 0.382 * (n - 1))]) solve_at(s) || return nothing end - if sampler.detect_uniform_M + if method.detect_uniform_M # a uniform M needs no fit, and so no support grid either probes = collect(values(solved)) all(==(first(probes)), probes) && return first(probes) end budget = clamp( - ceil(Int, sampler.budget * n), min(sampler.min_solves, n), n) + ceil(Int, method.budget * n), min(method.min_solves, n), n) X = _support_coords(support_grids()) while length(solved) < budget - ms, ss = _mean_sd(sampler, X, solved) - acq = ms .+ sampler.kappa .* ss + ms, ss = _mean_sd(gp, X, solved) + acq = ms .+ method.kappa .* ss for li in keys(solved) acq[li] = -Inf end @@ -123,9 +99,9 @@ function DP.sample_M_values( end return M_vals end - ms, ss = _mean_sd(sampler, X, solved) + ms, ss = _mean_sd(gp, X, solved) for (li, I) in enumerate(idxs) # exact M values are nonnegative - M_vals[I] = get(solved, li, max(ms[li] + sampler.kappa * ss[li], 0.0)) + M_vals[I] = get(solved, li, max(ms[li] + method.kappa * ss[li], 0.0)) end return M_vals end diff --git a/src/datatypes.jl b/src/datatypes.jl index 2bb9038..3196853 100644 --- a/src/datatypes.jl +++ b/src/datatypes.jl @@ -375,21 +375,48 @@ A type for using the multiple big-M reformulation approach for disjunctive const **Fields** - `optimizer::O`: Optimizer to use when solving mini-models (required). - `default_M::T`: Default big-M value to use if no big-M is specified for a logical variable (1e9). -- `M_sampler::Any`: Strategy for computing M values across the supports - of an infinite model (`:auto`). `:auto` uses [`GPSampler`](@ref) when - the AbstractGPs extension is loaded and `:exact` otherwise; `:exact` - solves an M subproblem at every support. Ignored for finite models. +- `gp::Any`: Gaussian process (or kernel) used to estimate the M values + across the supports of an infinite model (`nothing`). `nothing` solves + an M subproblem at every support; an `AbstractGPs.AbstractGP` or a + `KernelFunctions.Kernel` solves a subset of the supports and fills the + rest with a conservative GP estimate (see [`sample_M_values`](@ref)). + Ignored for finite models, as are the remaining fields. +- `kappa::Float64`: Upper-confidence-bound multiplier used to select the + next support to solve and to fill unsolved supports (2.5). +- `budget::Float64`: Fraction of the supports to solve exactly, in + `(0, 1]` (0.25). +- `min_solves::Int`: Minimum number of exactly solved supports (6). +- `detect_uniform_M::Bool`: If `true` (the default), M values that agree + at the first few supports are taken to be uniform and used for every + support. Set it to `false` to always fit the GP, which leaves the + usual `kappa * sd` cushion on the unsolved supports at the cost of + the extra solves, and which requires that every infinite parameter + have a support grid. """ mutable struct MBM{O, T} <: AbstractReformulationMethod optimizer::O default_M::T - M_sampler::Any + gp::Any + kappa::Float64 + budget::Float64 + min_solves::Int + detect_uniform_M::Bool # Constructor with optimizer (required) and optional default_M function MBM( - optimizer::O, default_M::T = 1e9; M_sampler = :auto + optimizer::O, + default_M::T = 1e9; + gp = nothing, + kappa::Real = 2.5, + budget::Real = 0.25, + min_solves::Int = 6, + detect_uniform_M::Bool = true ) where {O, T} - new{O, T}(optimizer, default_M, M_sampler) + kappa >= 0 || error("`kappa` must be nonnegative.") + 0 < budget <= 1 || error("`budget` must be in `(0, 1]`.") + min_solves >= 1 || error("`min_solves` must be at least 1.") + new{O, T}(optimizer, default_M, gp, Float64(kappa), + Float64(budget), min_solves, detect_uniform_M) end end @@ -397,7 +424,11 @@ mutable struct _MBM{O, T, M <: JuMP.AbstractModel} <: AbstractReformulationMetho optimizer::O M::Dict{LogicalVariableRef{M}, Any} default_M::T - M_sampler::Any + gp::Any + kappa::Float64 + budget::Float64 + min_solves::Int + detect_uniform_M::Bool subproblem_indicators::Vector{LogicalVariableRef{M}} # Cached submodels: indicator => GDPSubmodel. # Typed Any so extensions can store different types. @@ -408,7 +439,11 @@ mutable struct _MBM{O, T, M <: JuMP.AbstractModel} <: AbstractReformulationMetho method.optimizer, Dict{LogicalVariableRef{M}, Any}(), method.default_M, - method.M_sampler, + method.gp, + method.kappa, + method.budget, + method.min_solves, + method.detect_uniform_M, Vector{LogicalVariableRef{M}}(), Dict{LogicalVariableRef{M}, Any}() ) diff --git a/src/extension_api.jl b/src/extension_api.jl index 7c220b2..0e7e6aa 100644 --- a/src/extension_api.jl +++ b/src/extension_api.jl @@ -40,62 +40,31 @@ Y(t, x) function InfiniteLogical end """ - GPSampler(; kappa = 2.5, budget = 0.25, min_solves = 6, - kernel = nothing, detect_uniform_M = true) - -Creates a Gaussian-process M sampler for [`MBM`](@ref) on infinite -models. Instead of solving an M subproblem at every support of the -infinite parameters, the sampler solves a subset of the supports -selected by an upper-confidence-bound acquisition and fills the -remaining supports with the posterior upper confidence bound -`mean + kappa * sd`. The filled values are heuristic upper estimates -of the exact M values, not certificates; increase `kappa` for more -conservative estimates or use `MBM(...; M_sampler = :exact)` to solve -every support. This requires that InfiniteOpt and AbstractGPs be -imported first, in which case it is also the default M sampler (see -the `M_sampler` field of [`MBM`](@ref)). - -**Keyword Arguments** -- `kappa::Real`: Upper-confidence-bound multiplier used to select the - next support to solve and to fill unsolved supports (2.5). -- `budget::Real`: Fraction of the supports to solve exactly, in - `(0, 1]` (0.25). -- `min_solves::Int`: Minimum number of exactly solved supports (6). -- `kernel`: Covariance kernel for the GP fit. Defaults to a squared - exponential kernel whose lengthscale is selected by maximizing the - marginal likelihood; pass a `KernelFunctions` kernel to override. -- `detect_uniform_M::Bool`: If `true` (the default), M values that - agree at the first few supports are taken to be uniform and used - for every support. This is cheap and is what makes infinite - parameters without a support grid (e.g. dependent ones) workable, - but it assumes M does not vary elsewhere. Set it to `false` to - always fit the GP, which leaves the usual `kappa * sd` cushion on - the unsolved supports at the cost of the extra solves, and which - requires that every infinite parameter have a support grid. - -**Example** -```julia -julia> using DisjunctiveProgramming, InfiniteOpt, AbstractGPs, HiGHS - -julia> method = MBM(HiGHS.Optimizer, M_sampler = GPSampler(kappa = 4.0)) -``` -""" -function GPSampler end - -""" - sample_M_values(sampler, objectives, sub, method, support_grids) + sample_M_values(gp, objectives, sub, method, support_grids) Compute the MBM M values at the transcription supports of an infinite -model. `objectives` is the array of per-support objective expressions, -`sub` is the transcribed submodel wrapped as a `GDPSubmodel`, `method` -is the `_MBM` data, and `support_grids` is a function returning the +model. `gp` is the `gp` field of [`MBM`](@ref), `objectives` is the +array of per-support objective expressions, `sub` is the transcribed +submodel wrapped as a `GDPSubmodel`, `method` is the `_MBM` data (which +carries the sampling settings `kappa`, `budget`, `min_solves`, and +`detect_uniform_M`), and `support_grids` is a function returning the support vectors of the infinite parameters. It is a function because the supports are only well defined once M is known to vary over them, -so samplers that return early (or never need coordinates) must not -call it. Returns an array of M values shaped like -`objectives`, a scalar when M is uniform across the supports, or -`nothing` if an M subproblem is infeasible. Extensions implement -methods for their sampler types (e.g. [`GPSampler`](@ref)); the -`:exact` sampler solves every support. +so methods that return early (or never need coordinates) must not +call it. Returns an array of M values shaped like `objectives`, a +scalar when M is uniform across the supports, or `nothing` if an M +subproblem is infeasible. Extensions implement methods that dispatch +on `gp`: `nothing` solves an M subproblem at every support, while an +`AbstractGPs.AbstractGP` or a `KernelFunctions.Kernel` solves a +subset of the supports selected by an upper-confidence-bound +acquisition and fills the rest with the posterior upper confidence +bound `mean + kappa * sd`. The filled values are heuristic upper +estimates of the exact M values, not certificates. """ -function sample_M_values end +function sample_M_values(gp, objectives, sub, method, support_grids) + error("Unrecognized `gp` value `$(repr(gp))` for MBM on an " * + "infinite model. Use `nothing` to solve an M subproblem " * + "at every support, or load AbstractGPs and pass an " * + "`AbstractGPs.AbstractGP` or a `KernelFunctions.Kernel` " * + "to estimate M values with a Gaussian process.") +end diff --git a/src/mbm.jl b/src/mbm.jl index 4441d95..c5c17b4 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -88,7 +88,10 @@ function reformulate_disjunct_constraint( method::_MBM ) ref_cons = reformulate_disjunction(model, con, MBM( - method.optimizer, method.default_M, M_sampler = method.M_sampler)) + method.optimizer, method.default_M, gp = method.gp, + kappa = method.kappa, budget = method.budget, + min_solves = method.min_solves, + detect_uniform_M = method.detect_uniform_M)) new_ref_cons = Vector{JuMP.AbstractConstraint}() for ref_con in ref_cons append!(new_ref_cons, diff --git a/test/extensions/InfiniteDisjunctiveProgramming.jl b/test/extensions/InfiniteDisjunctiveProgramming.jl index 3526d62..710324c 100644 --- a/test/extensions/InfiniteDisjunctiveProgramming.jl +++ b/test/extensions/InfiniteDisjunctiveProgramming.jl @@ -376,7 +376,7 @@ function test_raw_M_infinite_scalar() @constraint(model, con, x >= 5, Disjunct(Y[1])) @constraint(model, con2, x <= 3, Disjunct(Y[2])) @disjunction(model, Y) - mbm = DP._MBM(MBM(HiGHS.Optimizer, M_sampler = :exact), model) + mbm = DP._MBM(MBM(HiGHS.Optimizer), model) sub = DP.copy_model_with_constraints( model, DP.DisjunctConstraintRef[con2], mbm) obj = DP.prepare_max_M_objective( @@ -399,7 +399,7 @@ function test_raw_M_infinite_param_function() @constraint(model, con, x <= f, Disjunct(Y[1])) @constraint(model, con2, x >= 0.5, Disjunct(Y[2])) @disjunction(model, Y) - mbm = DP._MBM(MBM(HiGHS.Optimizer, M_sampler = :exact), model) + mbm = DP._MBM(MBM(HiGHS.Optimizer), model) sub = DP.copy_model_with_constraints( model, DP.DisjunctConstraintRef[con2], mbm) obj = DP.prepare_max_M_objective( @@ -428,7 +428,7 @@ function test_raw_M_infinite_two_params() @constraint(model, con, x <= t + s, Disjunct(Y[1])) @constraint(model, con2, x >= 0.5, Disjunct(Y[2])) @disjunction(model, Y) - mbm = DP._MBM(MBM(HiGHS.Optimizer, M_sampler = :exact), model) + mbm = DP._MBM(MBM(HiGHS.Optimizer), model) sub = DP.copy_model_with_constraints( model, DP.DisjunctConstraintRef[con2], mbm) obj = DP.prepare_max_M_objective( @@ -452,14 +452,12 @@ function test_raw_M_infinite_dependent_params() @constraint(model, con, x >= 5, Disjunct(Y[1])) @constraint(model, con2, x <= 3, Disjunct(Y[2])) @disjunction(model, Y) - for sampler in (:exact, :auto) - mbm = DP._MBM(MBM(HiGHS.Optimizer, M_sampler = sampler), model) - sub = DP.copy_model_with_constraints( - model, DP.DisjunctConstraintRef[con2], mbm) - obj = DP.prepare_max_M_objective( - model, JuMP.constraint_object(con), sub) - @test DP.raw_M(sub, obj, mbm) == 5.0 - end + mbm = DP._MBM(MBM(HiGHS.Optimizer), model) + sub = DP.copy_model_with_constraints( + model, DP.DisjunctConstraintRef[con2], mbm) + obj = DP.prepare_max_M_objective( + model, JuMP.constraint_object(con), sub) + @test DP.raw_M(sub, obj, mbm) == 5.0 end # Piecewise-constant max-of-corners: returns the maximum value over diff --git a/test/extensions/InfiniteGPDisjunctiveProgramming.jl b/test/extensions/InfiniteGPDisjunctiveProgramming.jl index 796507b..107342f 100644 --- a/test/extensions/InfiniteGPDisjunctiveProgramming.jl +++ b/test/extensions/InfiniteGPDisjunctiveProgramming.jl @@ -1,40 +1,25 @@ -using InfiniteOpt, HiGHS, AbstractGPs, KernelFunctions +using InfiniteOpt, HiGHS, AbstractGPs import DisjunctiveProgramming as DP -# Helpers to access internal functions of the two extensions -const IGDP = Base.get_extension(DP, :InfiniteGPDisjunctiveProgramming) -const IDP = Base.get_extension(DP, :InfiniteDisjunctiveProgramming) - -function test_gp_sampler_creation() - sampler = GPSampler() - @test sampler isa IGDP.GPSampler - @test sampler.kappa == 2.5 - @test sampler.budget == 0.25 - @test sampler.min_solves == 6 - @test isnothing(sampler.kernel) - @test sampler.detect_uniform_M +function test_gp_mbm_kwargs() + method = MBM(HiGHS.Optimizer) + @test method.gp === nothing + @test method.kappa == 2.5 + @test method.budget == 0.25 + @test method.min_solves == 6 + @test method.detect_uniform_M kern = with_lengthscale(SqExponentialKernel(), 0.3) - sampler = GPSampler(kappa = 4.0, budget = 0.1, min_solves = 3, - kernel = kern, detect_uniform_M = false) - @test sampler.kappa == 4.0 - @test sampler.budget == 0.1 - @test sampler.min_solves == 3 - @test sampler.kernel === kern - @test !sampler.detect_uniform_M - @test_throws ErrorException GPSampler(kappa = -1) - @test_throws ErrorException GPSampler(budget = 0) - @test_throws ErrorException GPSampler(budget = 1.5) - @test_throws ErrorException GPSampler(min_solves = 0) -end - -function test_gp_sampler_resolution() - # the GP extension is loaded, so :auto resolves to a GPSampler - @test IDP._resolve_M_sampler(:auto) isa IGDP.GPSampler - @test IDP._resolve_M_sampler(:exact) === :exact - sampler = GPSampler(kappa = 3.0) - @test IDP._resolve_M_sampler(sampler) === sampler - @test MBM(HiGHS.Optimizer).M_sampler === :auto - @test MBM(HiGHS.Optimizer, M_sampler = :exact).M_sampler === :exact + method = MBM(HiGHS.Optimizer, gp = kern, kappa = 4.0, + budget = 0.1, min_solves = 3, detect_uniform_M = false) + @test method.gp === kern + @test method.kappa == 4.0 + @test method.budget == 0.1 + @test method.min_solves == 3 + @test !method.detect_uniform_M + @test_throws ErrorException MBM(HiGHS.Optimizer, kappa = -1) + @test_throws ErrorException MBM(HiGHS.Optimizer, budget = 0) + @test_throws ErrorException MBM(HiGHS.Optimizer, budget = 1.5) + @test_throws ErrorException MBM(HiGHS.Optimizer, min_solves = 0) end # Mirror of test_raw_M_infinite_scalar: uniform seed M values collapse @@ -47,7 +32,8 @@ function test_gp_raw_M_scalar() @constraint(model, con, x >= 5, Disjunct(Y[1])) @constraint(model, con2, x <= 3, Disjunct(Y[2])) @disjunction(model, Y) - mbm = DP._MBM(MBM(HiGHS.Optimizer), model) + mbm = DP._MBM( + MBM(HiGHS.Optimizer, gp = SqExponentialKernel()), model) sub = DP.copy_model_with_constraints( model, DP.DisjunctConstraintRef[con2], mbm) obj = DP.prepare_max_M_objective( @@ -59,7 +45,7 @@ end # sampler solves all of them exactly and must reproduce the exact # grid parameter function function test_gp_raw_M_matches_exact() - function pfunc_values(M_sampler, supports) + function pfunc_values(gp, supports) model = InfiniteGDPModel() @infinite_parameter(model, t ∈ [0, 1], supports = supports) @variable(model, 0 <= x <= 10, Infinite(t)) @@ -68,8 +54,7 @@ function test_gp_raw_M_matches_exact() @constraint(model, con, x <= f, Disjunct(Y[1])) @constraint(model, con2, x >= 0.5, Disjunct(Y[2])) @disjunction(model, Y) - mbm = DP._MBM( - MBM(HiGHS.Optimizer, M_sampler = M_sampler), model) + mbm = DP._MBM(MBM(HiGHS.Optimizer, gp = gp), model) sub = DP.copy_model_with_constraints( model, DP.DisjunctConstraintRef[con2], mbm) obj = DP.prepare_max_M_objective( @@ -79,11 +64,12 @@ function test_gp_raw_M_matches_exact() return [InfiniteOpt.raw_function(M)(t_val) for t_val in supports] end supports = [0.0, 0.25, 0.5, 0.75, 1.0] - exact_vals = pfunc_values(:exact, supports) - @test pfunc_values(GPSampler(), supports) == exact_vals - # a user kernel skips the lengthscale fit but solves the same supports - kern = with_lengthscale(SqExponentialKernel(), 0.2) - @test pfunc_values(GPSampler(kernel = kern), supports) == exact_vals + exact_vals = pfunc_values(nothing, supports) + # a kernel gets its lengthscale fit; a GP is used as given. Both + # solve the same supports here, so the M values match exactly + @test pfunc_values(SqExponentialKernel(), supports) == exact_vals + gp = GP(with_lengthscale(SqExponentialKernel(), 0.2)) + @test pfunc_values(gp, supports) == exact_vals end # an empty disjunct region makes the M subproblems infeasible; both @@ -103,16 +89,16 @@ function test_gp_infeasible_disjunct() @objective(model, Max, 𝔼(x, t)) return model end - for sampler in (:exact, GPSampler()) + for gp in (nothing, SqExponentialKernel()) model = build() @test_throws ErrorException optimize!(model, - gdp_method = MBM(HiGHS.Optimizer, M_sampler = sampler)) + gdp_method = MBM(HiGHS.Optimizer, gp = gp)) end end # optimum (10) needs M(t) >= 10 - 2t pointwise; the GP fill is heuristic function test_gp_mbm_solve_equivalence() - function solve_with(M_sampler) + function solve_with(method) model = InfiniteGDPModel(HiGHS.Optimizer) set_silent(model) @infinite_parameter(model, t ∈ [0, 1], num_supports = 20) @@ -123,20 +109,21 @@ function test_gp_mbm_solve_equivalence() @constraint(model, x >= 0.5, Disjunct(Y[2])) @disjunction(model, Y) @objective(model, Max, 𝔼(x, t)) - optimize!(model, - gdp_method = MBM(HiGHS.Optimizer, M_sampler = M_sampler)) + optimize!(model, gdp_method = method) @test termination_status(model) == MOI.OPTIMAL return objective_value(model) end - obj_exact = solve_with(:exact) - obj_auto = solve_with(:auto) - obj_gp = solve_with(GPSampler(kappa = 4.0, budget = 0.2)) + obj_exact = solve_with(MBM(HiGHS.Optimizer)) + obj_gp = solve_with( + MBM(HiGHS.Optimizer, gp = SqExponentialKernel())) + obj_tuned = solve_with(MBM(HiGHS.Optimizer, + gp = SqExponentialKernel(), kappa = 4.0, budget = 0.2)) @test obj_exact ≈ 10.0 atol = 1e-4 # over-M can't raise the optimum, under-M can only shave it a bit - @test obj_auto <= obj_exact + 1e-6 - @test obj_auto ≈ obj_exact atol = 1e-2 @test obj_gp <= obj_exact + 1e-6 @test obj_gp ≈ obj_exact atol = 1e-2 + @test obj_tuned <= obj_exact + 1e-6 + @test obj_tuned ≈ obj_exact atol = 1e-2 end # A periodic M must not read as uniform. With f(t) = 2|cos(2*pi*t)| @@ -145,7 +132,7 @@ end # which caps x at 8 and cuts the optimum from 10 down to 9. function test_gp_periodic_M_not_uniform() supports = [0.0, 0.25, 0.5, 0.75, 1.0] - function solve_with(M_sampler) + function solve_with(gp) model = InfiniteGDPModel(HiGHS.Optimizer) set_silent(model) @infinite_parameter(model, t ∈ [0, 1], supports = supports) @@ -156,12 +143,11 @@ function test_gp_periodic_M_not_uniform() @constraint(model, x >= 0.5, Disjunct(Y[2])) @disjunction(model, Y) @objective(model, Max, 𝔼(x, t)) - optimize!(model, gdp_method = MBM( - HiGHS.Optimizer, M_sampler = M_sampler)) + optimize!(model, gdp_method = MBM(HiGHS.Optimizer, gp = gp)) return objective_value(model) end - @test solve_with(:exact) ≈ 10.0 atol = 1e-6 - @test solve_with(GPSampler()) ≈ 10.0 atol = 1e-6 + @test solve_with(nothing) ≈ 10.0 atol = 1e-6 + @test solve_with(SqExponentialKernel()) ≈ 10.0 atol = 1e-6 end # With detection off the uniform M is not collapsed to a scalar: the @@ -177,7 +163,8 @@ function test_gp_detect_uniform_M_off() @constraint(model, con2, x <= 3, Disjunct(Y[2])) @disjunction(model, Y) mbm = DP._MBM(MBM(HiGHS.Optimizer, - M_sampler = GPSampler(detect_uniform_M = detect)), model) + gp = SqExponentialKernel(), + detect_uniform_M = detect), model) sub = DP.copy_model_with_constraints( model, DP.DisjunctConstraintRef[con2], mbm) obj = DP.prepare_max_M_objective( @@ -204,7 +191,8 @@ function test_gp_detect_uniform_M_off_dependent() @constraint(model, con2, x <= 3, Disjunct(Y[2])) @disjunction(model, Y) mbm = DP._MBM(MBM(HiGHS.Optimizer, - M_sampler = GPSampler(detect_uniform_M = false)), model) + gp = SqExponentialKernel(), + detect_uniform_M = false), model) sub = DP.copy_model_with_constraints( model, DP.DisjunctConstraintRef[con2], mbm) obj = DP.prepare_max_M_objective( @@ -212,7 +200,7 @@ function test_gp_detect_uniform_M_off_dependent() @test_throws ErrorException DP.raw_M(sub, obj, mbm) end -function test_gp_unknown_sampler_error() +function test_gp_unknown_gp_error() model = InfiniteGDPModel(HiGHS.Optimizer) set_silent(model) @infinite_parameter(model, t ∈ [0, 1], num_supports = 5) @@ -224,18 +212,17 @@ function test_gp_unknown_sampler_error() @disjunction(model, Y) @objective(model, Max, 𝔼(x, t)) @test_throws ErrorException optimize!(model, - gdp_method = MBM(HiGHS.Optimizer, M_sampler = :grid)) + gdp_method = MBM(HiGHS.Optimizer, gp = :grid)) end @testset "InfiniteGPDisjunctiveProgramming" begin - test_gp_sampler_creation() - test_gp_sampler_resolution() + test_gp_mbm_kwargs() test_gp_raw_M_scalar() test_gp_raw_M_matches_exact() test_gp_mbm_solve_equivalence() test_gp_periodic_M_not_uniform() test_gp_detect_uniform_M_off() test_gp_detect_uniform_M_off_dependent() - test_gp_unknown_sampler_error() + test_gp_unknown_gp_error() test_gp_infeasible_disjunct() end From 51b3fd98d6f60d73be471d714fc4e0921ff3b2ce Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Tue, 18 Aug 2026 14:48:58 -0400 Subject: [PATCH 07/10] AbstractGPsDisjunctiveProgramming --- Project.toml | 2 +- ...l => AbstractGPsDisjunctiveProgramming.jl} | 3 +- ...l => AbstractGPsDisjunctiveProgramming.jl} | 2 +- .../InfiniteDisjunctiveProgramming.jl | 31 +++++++++++++++++++ test/runtests.jl | 2 +- 5 files changed, 35 insertions(+), 5 deletions(-) rename ext/{InfiniteGPDisjunctiveProgramming.jl => AbstractGPsDisjunctiveProgramming.jl} (98%) rename test/extensions/{InfiniteGPDisjunctiveProgramming.jl => AbstractGPsDisjunctiveProgramming.jl} (99%) diff --git a/Project.toml b/Project.toml index 4fa94e8..868b4bc 100644 --- a/Project.toml +++ b/Project.toml @@ -12,8 +12,8 @@ AbstractGPs = "99985d1d-32ba-4be9-9821-2ec096f28918" InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57" [extensions] +AbstractGPsDisjunctiveProgramming = "AbstractGPs" InfiniteDisjunctiveProgramming = "InfiniteOpt" -InfiniteGPDisjunctiveProgramming = ["InfiniteOpt", "AbstractGPs"] [compat] AbstractGPs = "0.5" diff --git a/ext/InfiniteGPDisjunctiveProgramming.jl b/ext/AbstractGPsDisjunctiveProgramming.jl similarity index 98% rename from ext/InfiniteGPDisjunctiveProgramming.jl rename to ext/AbstractGPsDisjunctiveProgramming.jl index 1c46687..e8a37db 100644 --- a/ext/InfiniteGPDisjunctiveProgramming.jl +++ b/ext/AbstractGPsDisjunctiveProgramming.jl @@ -1,6 +1,5 @@ -module InfiniteGPDisjunctiveProgramming +module AbstractGPsDisjunctiveProgramming -import InfiniteOpt, JuMP import AbstractGPs import AbstractGPs.KernelFunctions import DisjunctiveProgramming as DP diff --git a/test/extensions/InfiniteGPDisjunctiveProgramming.jl b/test/extensions/AbstractGPsDisjunctiveProgramming.jl similarity index 99% rename from test/extensions/InfiniteGPDisjunctiveProgramming.jl rename to test/extensions/AbstractGPsDisjunctiveProgramming.jl index 107342f..2e70d03 100644 --- a/test/extensions/InfiniteGPDisjunctiveProgramming.jl +++ b/test/extensions/AbstractGPsDisjunctiveProgramming.jl @@ -215,7 +215,7 @@ function test_gp_unknown_gp_error() gdp_method = MBM(HiGHS.Optimizer, gp = :grid)) end -@testset "InfiniteGPDisjunctiveProgramming" begin +@testset "AbstractGPsDisjunctiveProgramming" begin test_gp_mbm_kwargs() test_gp_raw_M_scalar() test_gp_raw_M_matches_exact() diff --git a/test/extensions/InfiniteDisjunctiveProgramming.jl b/test/extensions/InfiniteDisjunctiveProgramming.jl index 710324c..6e96d8c 100644 --- a/test/extensions/InfiniteDisjunctiveProgramming.jl +++ b/test/extensions/InfiniteDisjunctiveProgramming.jl @@ -121,6 +121,36 @@ function test_requires_disaggregation() @test DP.requires_disaggregation(y) == true end +# Bound info for parameter refs in disjunct constraints: parameter +# functions report their support extrema, finite parameters a point, +# other parameters (-Inf, Inf); Hull/PSplit clamp the bounds to +# include 0 and error when a variable is missing bounds. +function test_parameter_bound_info() + model = InfiniteGDPModel() + @infinite_parameter(model, t in [0, 1], supports = [0.0, 0.5, 1.0]) + @infinite_parameter(model, s in [0, 1]) + @finite_parameter(model, p == 2.0) + @variable(model, 0 <= x <= 10, Infinite(t)) + @variable(model, y, Infinite(t)) + @parameter_function(model, pf == t -> 2t - 1) + @parameter_function(model, pf2 == (t, s) -> t + s) + @parameter_function(model, pf3 == s -> s) + @test DP.set_variable_bound_info(pf, BigM()) == (-1.0, 1.0) + @test DP.set_variable_bound_info(p, BigM()) == (2.0, 2.0) + @test DP.set_variable_bound_info(t, BigM()) == (-Inf, Inf) + @test DP.set_variable_bound_info(x, BigM()) == (0.0, 10.0) + @test DP.set_variable_bound_info(y, BigM()) == (-Inf, Inf) + # multi-parameter and support-less parameter functions fall back + @test DP.set_variable_bound_info(pf2, BigM()) == (-Inf, Inf) + @test DP.set_variable_bound_info(pf3, BigM()) == (-Inf, Inf) + # Hull and PSplit clamp the bounds to include 0 + @test DP.set_variable_bound_info(pf, Hull()) == (-1.0, 1.0) + @test DP.set_variable_bound_info(p, Hull()) == (0.0, 2.0) + @test DP.set_variable_bound_info(x, Hull()) == (0.0, 10.0) + @test DP.set_variable_bound_info(p, PSplit([[x]])) == (0.0, 2.0) + @test_throws ErrorException DP.set_variable_bound_info(y, Hull()) +end + function test_all_variables_infiniteopt() model = InfiniteGDPModel() @infinite_parameter(model, t ∈ [0, 1]) @@ -853,6 +883,7 @@ end test_is_parameter() test_is_parameter_concrete_dispatches() test_requires_disaggregation() + test_parameter_bound_info() test_variable_properties_infiniteopt() test_variable_properties_from_expr() test_variable_properties_from_quad_expr() diff --git a/test/runtests.jl b/test/runtests.jl index 569a12b..8034f6e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -24,5 +24,5 @@ include("constraints/disjunction.jl") include("print.jl") include("solve.jl") include("extensions/InfiniteDisjunctiveProgramming.jl") -include("extensions/InfiniteGPDisjunctiveProgramming.jl") +include("extensions/AbstractGPsDisjunctiveProgramming.jl") From 5f903607ebf1787d69c01cbd7ac1dcf4da1d3e22 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Tue, 18 Aug 2026 15:05:15 -0400 Subject: [PATCH 08/10] More keywords --- README.md | 2 +- ext/AbstractGPsDisjunctiveProgramming.jl | 30 ++++++------ src/datatypes.jl | 46 ++++++++++++++++-- src/extension_api.jl | 7 +-- src/mbm.jl | 4 +- .../AbstractGPsDisjunctiveProgramming.jl | 48 ++++++++++++++----- 6 files changed, 102 insertions(+), 35 deletions(-) diff --git a/README.md b/README.md index bbf80b4..5532b98 100644 --- a/README.md +++ b/README.md @@ -224,7 +224,7 @@ optimize!(model, gdp_method = Hull()) value(W) ``` -When the `MBM` reformulation is used on an infinite model, an M subproblem is solved at every support by default. Loading [AbstractGPs.jl](https://github.com/JuliaGaussianProcesses/AbstractGPs.jl) (`using AbstractGPs`) enables an additional extension that instead solves M at a subset of the supports and fills the rest with a conservative Gaussian-process estimate, which can substantially reduce the number of subproblem solves. To opt in, pass a GP or kernel via the `gp` keyword, e.g. `MBM(optimizer, gp = SqExponentialKernel())` (lengthscale selected by marginal likelihood) or `MBM(optimizer, gp = GP(with_lengthscale(Matern52Kernel(), 0.2)))` (used as given). The filled values are heuristic upper estimates rather than certificates; see the `MBM` docstring for the tuning knobs (`kappa`, `budget`, `min_solves`, `detect_uniform_M`). +When the `MBM` reformulation is used on an infinite model, an M subproblem is solved at every support by default. Loading [AbstractGPs.jl](https://github.com/JuliaGaussianProcesses/AbstractGPs.jl) (`using AbstractGPs`) enables an additional extension that instead solves M at a subset of the supports and fills the rest with a conservative Gaussian-process estimate, which can substantially reduce the number of subproblem solves. To opt in, pass a GP or kernel via the `gp` keyword, e.g. `MBM(optimizer, gp = SqExponentialKernel())` (lengthscale selected by marginal likelihood) or `MBM(optimizer, gp = GP(with_lengthscale(Matern52Kernel(), 0.2)))` (used as given). The filled values are heuristic upper estimates rather than certificates; see the `MBM` docstring for the tuning keywords (`kappa`, `budget`, `min_solves`, `detect_uniform_M`, `lengthscales`, `jitter`, `n_seeds`, `seeds`). ## Release Notes diff --git a/ext/AbstractGPsDisjunctiveProgramming.jl b/ext/AbstractGPsDisjunctiveProgramming.jl index e8a37db..0e2a611 100644 --- a/ext/AbstractGPsDisjunctiveProgramming.jl +++ b/ext/AbstractGPsDisjunctiveProgramming.jl @@ -7,9 +7,6 @@ import DisjunctiveProgramming as DP ################################################################################ # GP FITTING ################################################################################ -const _LENGTHSCALES = (0.05, 0.1, 0.2, 0.4, 0.8) -const _JITTER = 1e-8 - # Normalized to [0, 1]^d so one lengthscale works across dimensions function _support_coords(grids) idxs = CartesianIndices(length.(grids)) @@ -20,15 +17,15 @@ function _support_coords(grids) end # A user GP is used as the prior directly; a kernel gets its -# lengthscale selected by marginal likelihood -function _fit_posterior(gp::AbstractGPs.AbstractGP, X, y) - return AbstractGPs.posterior(gp(X, _JITTER), y) +# lengthscale selected by marginal likelihood over the candidates +function _fit_posterior(gp::AbstractGPs.AbstractGP, X, y, method) + return AbstractGPs.posterior(gp(X, method.jitter), y) end -function _fit_posterior(kernel::KernelFunctions.Kernel, X, y) +function _fit_posterior(kernel::KernelFunctions.Kernel, X, y, method) best_post, best_lp = nothing, -Inf - for ls in _LENGTHSCALES + for ls in method.lengthscales kern = KernelFunctions.with_lengthscale(kernel, ls) - fx = AbstractGPs.GP(kern)(X, _JITTER) + fx = AbstractGPs.GP(kern)(X, method.jitter) lp = AbstractGPs.logpdf(fx, y) if lp > best_lp best_post, best_lp = AbstractGPs.posterior(fx, y), lp @@ -37,14 +34,14 @@ function _fit_posterior(kernel::KernelFunctions.Kernel, X, y) return best_post end -function _mean_sd(gp, X, solved) +function _mean_sd(gp, X, solved, method) lis = collect(keys(solved)) y = [solved[li] for li in lis] ybar = sum(y) / length(y) # floored so near-equal solved values still cushion the filled ones ystd = max(sqrt(sum(abs2, y .- ybar) / max(length(y) - 1, 1)), 1e-2 * abs(ybar), 1e-8) - post = _fit_posterior(gp, X[lis], (y .- ybar) ./ ystd) + post = _fit_posterior(gp, X[lis], (y .- ybar) ./ ystd, method) mz = AbstractGPs.mean(post, X) vz = max.(AbstractGPs.var(post, X), 0.0) return mz .* ystd .+ ybar, sqrt.(vz) .* ystd @@ -70,9 +67,10 @@ function DP.sample_M_values( solved[li] = m return true end - # golden fractions; even spacing aliases with a periodic M - for s in unique([1, n, 1 + floor(Int, 0.618 * (n - 1)), - 1 + floor(Int, 0.382 * (n - 1))]) + # user-given seed fractions, or an evenly spaced grid of n_seeds + fracs = something(method.seeds, + range(0, 1, length = method.n_seeds)) + for s in unique(1 .+ round.(Int, fracs .* (n - 1))) solve_at(s) || return nothing end if method.detect_uniform_M @@ -84,7 +82,7 @@ function DP.sample_M_values( ceil(Int, method.budget * n), min(method.min_solves, n), n) X = _support_coords(support_grids()) while length(solved) < budget - ms, ss = _mean_sd(gp, X, solved) + ms, ss = _mean_sd(gp, X, solved, method) acq = ms .+ method.kappa .* ss for li in keys(solved) acq[li] = -Inf @@ -98,7 +96,7 @@ function DP.sample_M_values( end return M_vals end - ms, ss = _mean_sd(gp, X, solved) + ms, ss = _mean_sd(gp, X, solved, method) for (li, I) in enumerate(idxs) # exact M values are nonnegative M_vals[I] = get(solved, li, max(ms[li] + method.kappa * ss[li], 0.0)) end diff --git a/src/datatypes.jl b/src/datatypes.jl index 3196853..40812e4 100644 --- a/src/datatypes.jl +++ b/src/datatypes.jl @@ -387,11 +387,24 @@ A type for using the multiple big-M reformulation approach for disjunctive const `(0, 1]` (0.25). - `min_solves::Int`: Minimum number of exactly solved supports (6). - `detect_uniform_M::Bool`: If `true` (the default), M values that agree - at the first few supports are taken to be uniform and used for every + at the seed supports are taken to be uniform and used for every support. Set it to `false` to always fit the GP, which leaves the usual `kappa * sd` cushion on the unsolved supports at the cost of the extra solves, and which requires that every infinite parameter have a support grid. +- `lengthscales::Vector{Float64}`: Candidate lengthscales for the + marginal-likelihood kernel fit when `gp` is a kernel, relative to + support coordinates normalized to `[0, 1]`; ignored when `gp` is a + full Gaussian process ([0.05, 0.1, 0.2, 0.4, 0.8]). +- `jitter::Float64`: Observation-noise nugget added to the GP prior + when fitting (1e-8). +- `n_seeds::Int`: Number of supports solved before the first GP fit, + taken as an evenly spaced grid over the supports (4). +- `seeds::Any`: Optional vector of fractions in `[0, 1]` giving the + positions along the support grid to solve before the first GP fit, + overriding the evenly spaced grid (`nothing`). With + `detect_uniform_M`, evenly spaced seeds can read a periodic M as + uniform; pass unevenly spaced seeds to guard against that. """ mutable struct MBM{O, T} <: AbstractReformulationMethod optimizer::O @@ -401,6 +414,10 @@ mutable struct MBM{O, T} <: AbstractReformulationMethod budget::Float64 min_solves::Int detect_uniform_M::Bool + lengthscales::Vector{Float64} + jitter::Float64 + n_seeds::Int + seeds::Any # Constructor with optimizer (required) and optional default_M function MBM( @@ -410,13 +427,28 @@ mutable struct MBM{O, T} <: AbstractReformulationMethod kappa::Real = 2.5, budget::Real = 0.25, min_solves::Int = 6, - detect_uniform_M::Bool = true + detect_uniform_M::Bool = true, + lengthscales = [0.05, 0.1, 0.2, 0.4, 0.8], + jitter::Real = 1e-8, + n_seeds::Int = 4, + seeds = nothing ) where {O, T} kappa >= 0 || error("`kappa` must be nonnegative.") 0 < budget <= 1 || error("`budget` must be in `(0, 1]`.") min_solves >= 1 || error("`min_solves` must be at least 1.") + lengthscales = collect(Float64, lengthscales) + (!isempty(lengthscales) && all(>(0), lengthscales)) || + error("`lengthscales` must be positive and nonempty.") + jitter >= 0 || error("`jitter` must be nonnegative.") + n_seeds >= 2 || error("`n_seeds` must be at least 2.") + if seeds !== nothing + seeds = collect(Float64, seeds) + (!isempty(seeds) && all(f -> 0 <= f <= 1, seeds)) || + error("`seeds` must be fractions in `[0, 1]`.") + end new{O, T}(optimizer, default_M, gp, Float64(kappa), - Float64(budget), min_solves, detect_uniform_M) + Float64(budget), min_solves, detect_uniform_M, + lengthscales, Float64(jitter), n_seeds, seeds) end end @@ -429,6 +461,10 @@ mutable struct _MBM{O, T, M <: JuMP.AbstractModel} <: AbstractReformulationMetho budget::Float64 min_solves::Int detect_uniform_M::Bool + lengthscales::Vector{Float64} + jitter::Float64 + n_seeds::Int + seeds::Any subproblem_indicators::Vector{LogicalVariableRef{M}} # Cached submodels: indicator => GDPSubmodel. # Typed Any so extensions can store different types. @@ -444,6 +480,10 @@ mutable struct _MBM{O, T, M <: JuMP.AbstractModel} <: AbstractReformulationMetho method.budget, method.min_solves, method.detect_uniform_M, + method.lengthscales, + method.jitter, + method.n_seeds, + method.seeds, Vector{LogicalVariableRef{M}}(), Dict{LogicalVariableRef{M}, Any}() ) diff --git a/src/extension_api.jl b/src/extension_api.jl index 0e7e6aa..be7ae48 100644 --- a/src/extension_api.jl +++ b/src/extension_api.jl @@ -45,9 +45,10 @@ function InfiniteLogical end Compute the MBM M values at the transcription supports of an infinite model. `gp` is the `gp` field of [`MBM`](@ref), `objectives` is the array of per-support objective expressions, `sub` is the transcribed -submodel wrapped as a `GDPSubmodel`, `method` is the `_MBM` data (which -carries the sampling settings `kappa`, `budget`, `min_solves`, and -`detect_uniform_M`), and `support_grids` is a function returning the +submodel wrapped as a `GDPSubmodel`, `method` is the `_MBM` data +(which carries the sampling settings `kappa`, `budget`, `min_solves`, +`detect_uniform_M`, `lengthscales`, `jitter`, `n_seeds`, and `seeds`), +and `support_grids` is a function returning the support vectors of the infinite parameters. It is a function because the supports are only well defined once M is known to vary over them, so methods that return early (or never need coordinates) must not diff --git a/src/mbm.jl b/src/mbm.jl index c5c17b4..ef22f66 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -91,7 +91,9 @@ function reformulate_disjunct_constraint( method.optimizer, method.default_M, gp = method.gp, kappa = method.kappa, budget = method.budget, min_solves = method.min_solves, - detect_uniform_M = method.detect_uniform_M)) + detect_uniform_M = method.detect_uniform_M, + lengthscales = method.lengthscales, jitter = method.jitter, + n_seeds = method.n_seeds, seeds = method.seeds)) new_ref_cons = Vector{JuMP.AbstractConstraint}() for ref_con in ref_cons append!(new_ref_cons, diff --git a/test/extensions/AbstractGPsDisjunctiveProgramming.jl b/test/extensions/AbstractGPsDisjunctiveProgramming.jl index 2e70d03..dd9ffa0 100644 --- a/test/extensions/AbstractGPsDisjunctiveProgramming.jl +++ b/test/extensions/AbstractGPsDisjunctiveProgramming.jl @@ -8,18 +8,37 @@ function test_gp_mbm_kwargs() @test method.budget == 0.25 @test method.min_solves == 6 @test method.detect_uniform_M + @test method.lengthscales == [0.05, 0.1, 0.2, 0.4, 0.8] + @test method.jitter == 1e-8 + @test method.n_seeds == 4 + @test method.seeds === nothing kern = with_lengthscale(SqExponentialKernel(), 0.3) method = MBM(HiGHS.Optimizer, gp = kern, kappa = 4.0, - budget = 0.1, min_solves = 3, detect_uniform_M = false) + budget = 0.1, min_solves = 3, detect_uniform_M = false, + lengthscales = (0.1, 0.3), jitter = 1e-6, n_seeds = 6, + seeds = [0.0, 0.3, 1.0]) @test method.gp === kern @test method.kappa == 4.0 @test method.budget == 0.1 @test method.min_solves == 3 @test !method.detect_uniform_M + @test method.lengthscales == [0.1, 0.3] + @test method.jitter == 1e-6 + @test method.n_seeds == 6 + @test method.seeds == [0.0, 0.3, 1.0] @test_throws ErrorException MBM(HiGHS.Optimizer, kappa = -1) @test_throws ErrorException MBM(HiGHS.Optimizer, budget = 0) @test_throws ErrorException MBM(HiGHS.Optimizer, budget = 1.5) @test_throws ErrorException MBM(HiGHS.Optimizer, min_solves = 0) + @test_throws ErrorException MBM(HiGHS.Optimizer, + lengthscales = Float64[]) + @test_throws ErrorException MBM(HiGHS.Optimizer, + lengthscales = [-0.1]) + @test_throws ErrorException MBM(HiGHS.Optimizer, jitter = -1) + @test_throws ErrorException MBM(HiGHS.Optimizer, n_seeds = 1) + @test_throws ErrorException MBM(HiGHS.Optimizer, seeds = [1.5]) + @test_throws ErrorException MBM(HiGHS.Optimizer, + seeds = Float64[]) end # Mirror of test_raw_M_infinite_scalar: uniform seed M values collapse @@ -126,13 +145,15 @@ function test_gp_mbm_solve_equivalence() @test obj_tuned ≈ obj_exact atol = 1e-2 end -# A periodic M must not read as uniform. With f(t) = 2|cos(2*pi*t)| -# on these supports, M = 10 - f is 8 at supports 1, 3, 5 and 10 at -# supports 2, 4, so evenly spaced probes alias and collapse M to 8, -# which caps x at 8 and cuts the optimum from 10 down to 9. -function test_gp_periodic_M_not_uniform() +# Seed placement vs a periodic M. With f(t) = 2|cos(2*pi*t)| on these +# supports, M = 10 - f is 8 at supports 1, 3, 5 and 10 at supports +# 2, 4. Seeds that only hit the M = 8 supports alias the periodic M +# to uniform 8, which caps x at 8 and cuts the optimum from 10 down +# to 9; the default and denser seed grids see both values and stay +# exact. +function test_gp_periodic_M_seeds() supports = [0.0, 0.25, 0.5, 0.75, 1.0] - function solve_with(gp) + function solve_with(; kwargs...) model = InfiniteGDPModel(HiGHS.Optimizer) set_silent(model) @infinite_parameter(model, t ∈ [0, 1], supports = supports) @@ -143,11 +164,16 @@ function test_gp_periodic_M_not_uniform() @constraint(model, x >= 0.5, Disjunct(Y[2])) @disjunction(model, Y) @objective(model, Max, 𝔼(x, t)) - optimize!(model, gdp_method = MBM(HiGHS.Optimizer, gp = gp)) + optimize!(model, + gdp_method = MBM(HiGHS.Optimizer; kwargs...)) return objective_value(model) end - @test solve_with(nothing) ≈ 10.0 atol = 1e-6 - @test solve_with(SqExponentialKernel()) ≈ 10.0 atol = 1e-6 + @test solve_with() ≈ 10.0 atol = 1e-6 + @test solve_with(gp = SqExponentialKernel()) ≈ 10.0 atol = 1e-6 + @test solve_with(gp = SqExponentialKernel(), n_seeds = 5) ≈ + 10.0 atol = 1e-6 + @test solve_with(gp = SqExponentialKernel(), + seeds = [0.0, 0.5, 1.0]) ≈ 9.0 atol = 1e-6 end # With detection off the uniform M is not collapsed to a scalar: the @@ -220,7 +246,7 @@ end test_gp_raw_M_scalar() test_gp_raw_M_matches_exact() test_gp_mbm_solve_equivalence() - test_gp_periodic_M_not_uniform() + test_gp_periodic_M_seeds() test_gp_detect_uniform_M_off() test_gp_detect_uniform_M_off_dependent() test_gp_unknown_gp_error() From 3721c7b416445f38adfd47e0e942dee05587b673 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Wed, 19 Aug 2026 11:04:03 -0400 Subject: [PATCH 09/10] GPSampler config object --- README.md | 4 +- ext/AbstractGPsDisjunctiveProgramming.jl | 40 +++-- ext/InfiniteDisjunctiveProgramming.jl | 8 +- src/datatypes.jl | 151 ++++++++++-------- src/extension_api.jl | 33 ++-- src/mbm.jl | 8 +- .../AbstractGPsDisjunctiveProgramming.jl | 142 ++++++++-------- 7 files changed, 193 insertions(+), 193 deletions(-) diff --git a/README.md b/README.md index 5532b98..97c29a2 100644 --- a/README.md +++ b/README.md @@ -176,7 +176,7 @@ The following reformulation methods are currently supported: - `optimizer`: Optimizer to use when solving subproblems to determine M values. This is a required value. - `default_M`: Default big-M value to use if no big-M is specified for a logical variable (1e9). - - `gp`: Gaussian process (or kernel) used to estimate the M values across the supports of an infinite model. Default: `nothing`, which solves an M subproblem at every support. Pass an `AbstractGPs.AbstractGP` or a `KernelFunctions.Kernel` to instead solve a subset of the supports and fill the rest with a conservative GP estimate. Ignored for finite models, as are the tuning keywords `kappa`, `budget`, `min_solves`, and `detect_uniform_M` (see the `MBM` docstring). + - `sampler`: M-value sampler for infinite models. Default: `nothing`, which solves an M subproblem at every support. Pass a `GPSampler` to instead solve a subset of the supports and fill the rest with a conservative Gaussian-process estimate. Ignored for finite models. 5. [P-Split](https://arxiv.org/abs/2202.05198): This method reformulates each disjunct constraint into P constraints, each with a partitioned group defined by the user. This method requires that terms in the constraint be convex additively seperable with respect to each variable. The `PSplit` struct is created with the following required arguments: @@ -224,7 +224,7 @@ optimize!(model, gdp_method = Hull()) value(W) ``` -When the `MBM` reformulation is used on an infinite model, an M subproblem is solved at every support by default. Loading [AbstractGPs.jl](https://github.com/JuliaGaussianProcesses/AbstractGPs.jl) (`using AbstractGPs`) enables an additional extension that instead solves M at a subset of the supports and fills the rest with a conservative Gaussian-process estimate, which can substantially reduce the number of subproblem solves. To opt in, pass a GP or kernel via the `gp` keyword, e.g. `MBM(optimizer, gp = SqExponentialKernel())` (lengthscale selected by marginal likelihood) or `MBM(optimizer, gp = GP(with_lengthscale(Matern52Kernel(), 0.2)))` (used as given). The filled values are heuristic upper estimates rather than certificates; see the `MBM` docstring for the tuning keywords (`kappa`, `budget`, `min_solves`, `detect_uniform_M`, `lengthscales`, `jitter`, `n_seeds`, `seeds`). +When the `MBM` reformulation is used on an infinite model, an M subproblem is solved at every support by default. Loading [AbstractGPs.jl](https://github.com/JuliaGaussianProcesses/AbstractGPs.jl) (`using AbstractGPs`) enables an additional extension that instead solves M at a subset of the supports and fills the rest with a conservative Gaussian-process estimate, which can substantially reduce the number of subproblem solves. To opt in, pass a `GPSampler` via the `sampler` keyword, e.g. `MBM(optimizer, sampler = GPSampler())` (squared exponential kernel with the lengthscale selected by marginal likelihood), `GPSampler(Matern52Kernel())` (a custom kernel, lengthscale selected the same way), or `GPSampler(Matern52Kernel(), lengthscales = [0.2])` (lengthscale pinned). The filled values are heuristic upper estimates rather than certificates; see the `GPSampler` docstring for the tuning keywords (`kappa`, `budget`, `detect_uniform_M`, `lengthscales`, `jitter`, `seeds`). ## Release Notes diff --git a/ext/AbstractGPsDisjunctiveProgramming.jl b/ext/AbstractGPsDisjunctiveProgramming.jl index 0e2a611..09e3d8e 100644 --- a/ext/AbstractGPsDisjunctiveProgramming.jl +++ b/ext/AbstractGPsDisjunctiveProgramming.jl @@ -16,16 +16,12 @@ function _support_coords(grids) for I in vec(idxs)] end -# A user GP is used as the prior directly; a kernel gets its -# lengthscale selected by marginal likelihood over the candidates -function _fit_posterior(gp::AbstractGPs.AbstractGP, X, y, method) - return AbstractGPs.posterior(gp(X, method.jitter), y) -end -function _fit_posterior(kernel::KernelFunctions.Kernel, X, y, method) +# Lengthscale selected by marginal likelihood over the candidates +function _fit_posterior(kernel, X, y, sampler) best_post, best_lp = nothing, -Inf - for ls in method.lengthscales + for ls in sampler.lengthscales kern = KernelFunctions.with_lengthscale(kernel, ls) - fx = AbstractGPs.GP(kern)(X, method.jitter) + fx = AbstractGPs.GP(kern)(X, sampler.jitter) lp = AbstractGPs.logpdf(fx, y) if lp > best_lp best_post, best_lp = AbstractGPs.posterior(fx, y), lp @@ -34,14 +30,16 @@ function _fit_posterior(kernel::KernelFunctions.Kernel, X, y, method) return best_post end -function _mean_sd(gp, X, solved, method) +function _mean_sd(sampler, X, solved) lis = collect(keys(solved)) y = [solved[li] for li in lis] ybar = sum(y) / length(y) # floored so near-equal solved values still cushion the filled ones ystd = max(sqrt(sum(abs2, y .- ybar) / max(length(y) - 1, 1)), 1e-2 * abs(ybar), 1e-8) - post = _fit_posterior(gp, X[lis], (y .- ybar) ./ ystd, method) + kernel = something(sampler.kernel, + KernelFunctions.SqExponentialKernel()) + post = _fit_posterior(kernel, X[lis], (y .- ybar) ./ ystd, sampler) mz = AbstractGPs.mean(post, X) vz = max.(AbstractGPs.var(post, X), 0.0) return mz .* ystd .+ ybar, sqrt.(vz) .* ystd @@ -52,7 +50,7 @@ end ################################################################################ # Solve M at max-UCB selected supports, fill the rest with the bound function DP.sample_M_values( - gp::Union{AbstractGPs.AbstractGP, KernelFunctions.Kernel}, + sampler::DP.GPSampler, objectives::AbstractArray, sub::DP.GDPSubmodel, method::DP._MBM, @@ -67,23 +65,22 @@ function DP.sample_M_values( solved[li] = m return true end - # user-given seed fractions, or an evenly spaced grid of n_seeds - fracs = something(method.seeds, - range(0, 1, length = method.n_seeds)) + # an evenly spaced seed count, or user-given seed fractions + fracs = sampler.seeds isa Int ? + range(0, 1, length = sampler.seeds) : sampler.seeds for s in unique(1 .+ round.(Int, fracs .* (n - 1))) solve_at(s) || return nothing end - if method.detect_uniform_M + if sampler.detect_uniform_M # a uniform M needs no fit, and so no support grid either probes = collect(values(solved)) all(==(first(probes)), probes) && return first(probes) end - budget = clamp( - ceil(Int, method.budget * n), min(method.min_solves, n), n) + budget = min(ceil(Int, sampler.budget * n), n) X = _support_coords(support_grids()) while length(solved) < budget - ms, ss = _mean_sd(gp, X, solved, method) - acq = ms .+ method.kappa .* ss + ms, ss = _mean_sd(sampler, X, solved) + acq = ms .+ sampler.kappa .* ss for li in keys(solved) acq[li] = -Inf end @@ -96,9 +93,10 @@ function DP.sample_M_values( end return M_vals end - ms, ss = _mean_sd(gp, X, solved, method) + ms, ss = _mean_sd(sampler, X, solved) for (li, I) in enumerate(idxs) # exact M values are nonnegative - M_vals[I] = get(solved, li, max(ms[li] + method.kappa * ss[li], 0.0)) + M_vals[I] = get(solved, li, + max(ms[li] + sampler.kappa * ss[li], 0.0)) end return M_vals end diff --git a/ext/InfiniteDisjunctiveProgramming.jl b/ext/InfiniteDisjunctiveProgramming.jl index be572ec..2a492e6 100644 --- a/ext/InfiniteDisjunctiveProgramming.jl +++ b/ext/InfiniteDisjunctiveProgramming.jl @@ -333,7 +333,7 @@ end # Solve the M subproblem exactly at every support function DP.sample_M_values( - gp::Nothing, + sampler::Nothing, objectives::AbstractArray, sub::DP.GDPSubmodel, method::DP._MBM, @@ -349,7 +349,7 @@ function DP.sample_M_values( end # Transcribe mini_expr, compute the per-support M values with the -# method's gp, and aggregate to a scalar if uniform, else to a +# method's sampler, and aggregate to a scalar if uniform, else to a # parameter function on main. function DP.raw_M( sub::DP.GDPSubmodel{<:InfiniteOpt.InfiniteModel}, @@ -366,8 +366,8 @@ function DP.raw_M( transcribed = InfiniteOpt.transformation_model(sub.model) inner_sub = DP.GDPSubmodel(transcribed, JuMP.VariableRef[], Dict{JuMP.VariableRef, Vector{JuMP.VariableRef}}()) - M_vals = DP.sample_M_values(method.gp, objectives, inner_sub, - method, () -> _support_grids(sub, mini_expr)[2]) + M_vals = DP.sample_M_values(method.sampler, objectives, + inner_sub, method, () -> _support_grids(sub, mini_expr)[2]) M_vals === nothing && return nothing M_vals isa Number && return M_vals all(==(first(M_vals)), M_vals) && return first(M_vals) diff --git a/src/datatypes.jl b/src/datatypes.jl index 40812e4..ef82bce 100644 --- a/src/datatypes.jl +++ b/src/datatypes.jl @@ -375,80 +375,113 @@ A type for using the multiple big-M reformulation approach for disjunctive const **Fields** - `optimizer::O`: Optimizer to use when solving mini-models (required). - `default_M::T`: Default big-M value to use if no big-M is specified for a logical variable (1e9). -- `gp::Any`: Gaussian process (or kernel) used to estimate the M values - across the supports of an infinite model (`nothing`). `nothing` solves - an M subproblem at every support; an `AbstractGPs.AbstractGP` or a - `KernelFunctions.Kernel` solves a subset of the supports and fills the - rest with a conservative GP estimate (see [`sample_M_values`](@ref)). - Ignored for finite models, as are the remaining fields. -- `kappa::Float64`: Upper-confidence-bound multiplier used to select the - next support to solve and to fill unsolved supports (2.5). -- `budget::Float64`: Fraction of the supports to solve exactly, in - `(0, 1]` (0.25). -- `min_solves::Int`: Minimum number of exactly solved supports (6). -- `detect_uniform_M::Bool`: If `true` (the default), M values that agree - at the seed supports are taken to be uniform and used for every - support. Set it to `false` to always fit the GP, which leaves the - usual `kappa * sd` cushion on the unsolved supports at the cost of - the extra solves, and which requires that every infinite parameter - have a support grid. -- `lengthscales::Vector{Float64}`: Candidate lengthscales for the - marginal-likelihood kernel fit when `gp` is a kernel, relative to - support coordinates normalized to `[0, 1]`; ignored when `gp` is a - full Gaussian process ([0.05, 0.1, 0.2, 0.4, 0.8]). -- `jitter::Float64`: Observation-noise nugget added to the GP prior - when fitting (1e-8). -- `n_seeds::Int`: Number of supports solved before the first GP fit, - taken as an evenly spaced grid over the supports (4). -- `seeds::Any`: Optional vector of fractions in `[0, 1]` giving the - positions along the support grid to solve before the first GP fit, - overriding the evenly spaced grid (`nothing`). With - `detect_uniform_M`, evenly spaced seeds can read a periodic M as - uniform; pass unevenly spaced seeds to guard against that. +- `sampler::Any`: M-value sampler for infinite models (`nothing`). + `nothing` solves an M subproblem at every support; a + [`GPSampler`](@ref) solves a subset of the supports and fills the + rest with a conservative Gaussian-process estimate (see + [`sample_M_values`](@ref)). Ignored for finite models. """ mutable struct MBM{O, T} <: AbstractReformulationMethod optimizer::O default_M::T - gp::Any + sampler::Any + + # Constructor with optimizer (required) and optional default_M + function MBM( + optimizer::O, default_M::T = 1e9; sampler = nothing + ) where {O, T} + new{O, T}(optimizer, default_M, sampler) + end +end + +""" + GPSampler( + kernel = nothing; + kappa::Real = 2.5, + budget::Real = 0.25, + detect_uniform_M::Bool = true, + lengthscales = [0.05, 0.1, 0.2, 0.4, 0.8], + jitter::Real = 1e-8, + seeds = 4 + ) + +A Gaussian-process M sampler for the `sampler` field of [`MBM`](@ref) +on infinite models. Instead of solving an M subproblem at every +support, it solves the seed supports, then the supports selected by +an upper-confidence-bound acquisition until the budget is spent, and +fills the remaining supports with the posterior upper confidence +bound `mean + kappa * sd`. The filled values are heuristic upper +estimates of the exact M values, not certificates. Using it requires +that AbstractGPs be loaded. + +**Arguments** +- `kernel`: Covariance kernel for the GP fit; `nothing` (the + default) uses a squared exponential kernel. The lengthscale is + selected from `lengthscales` by marginal likelihood. +- `kappa::Real`: Upper-confidence-bound multiplier used to select the + next support to solve and to fill unsolved supports (2.5). +- `budget::Real`: Fraction of the supports to solve exactly, in + `(0, 1]` (0.25). The seed supports are always solved; + `budget = 1.0` solves every support. +- `detect_uniform_M::Bool`: If `true` (the default), M values that + agree at the seed supports are taken to be uniform and used for + every support. Set it to `false` to always fit the GP, which + leaves the usual `kappa * sd` cushion on the unsolved supports at + the cost of the extra solves, and which requires that every + infinite parameter have a support grid. +- `lengthscales`: Candidate lengthscales for the marginal-likelihood + kernel fit, relative to support coordinates normalized to + `[0, 1]` ([0.05, 0.1, 0.2, 0.4, 0.8]). Give a single candidate to + pin the lengthscale. +- `jitter::Real`: Observation-noise nugget added to the GP prior + when fitting (1e-8). +- `seeds`: Number of evenly spaced supports solved before the first + GP fit (4), or a vector of fractions in `[0, 1]` giving their + positions along the support grid. With `detect_uniform_M`, evenly + spaced seeds can read a periodic M as uniform; pass unevenly + spaced fractions to guard against that. + +**Example** +```julia +julia> using DisjunctiveProgramming, InfiniteOpt, AbstractGPs, HiGHS + +julia> method = MBM(HiGHS.Optimizer, sampler = GPSampler(kappa = 4.0)) +``` +""" +struct GPSampler + # Typed Any so base needs no AbstractGPs dependency + kernel::Any kappa::Float64 budget::Float64 - min_solves::Int detect_uniform_M::Bool lengthscales::Vector{Float64} jitter::Float64 - n_seeds::Int - seeds::Any + seeds::Union{Int, Vector{Float64}} - # Constructor with optimizer (required) and optional default_M - function MBM( - optimizer::O, - default_M::T = 1e9; - gp = nothing, + function GPSampler( + kernel = nothing; kappa::Real = 2.5, budget::Real = 0.25, - min_solves::Int = 6, detect_uniform_M::Bool = true, lengthscales = [0.05, 0.1, 0.2, 0.4, 0.8], jitter::Real = 1e-8, - n_seeds::Int = 4, - seeds = nothing - ) where {O, T} + seeds = 4 + ) kappa >= 0 || error("`kappa` must be nonnegative.") 0 < budget <= 1 || error("`budget` must be in `(0, 1]`.") - min_solves >= 1 || error("`min_solves` must be at least 1.") lengthscales = collect(Float64, lengthscales) (!isempty(lengthscales) && all(>(0), lengthscales)) || error("`lengthscales` must be positive and nonempty.") jitter >= 0 || error("`jitter` must be nonnegative.") - n_seeds >= 2 || error("`n_seeds` must be at least 2.") - if seeds !== nothing + if seeds isa Int + seeds >= 2 || error("`seeds` must be at least 2.") + else seeds = collect(Float64, seeds) (!isempty(seeds) && all(f -> 0 <= f <= 1, seeds)) || error("`seeds` must be fractions in `[0, 1]`.") end - new{O, T}(optimizer, default_M, gp, Float64(kappa), - Float64(budget), min_solves, detect_uniform_M, - lengthscales, Float64(jitter), n_seeds, seeds) + new(kernel, Float64(kappa), Float64(budget), + detect_uniform_M, lengthscales, Float64(jitter), seeds) end end @@ -456,15 +489,7 @@ mutable struct _MBM{O, T, M <: JuMP.AbstractModel} <: AbstractReformulationMetho optimizer::O M::Dict{LogicalVariableRef{M}, Any} default_M::T - gp::Any - kappa::Float64 - budget::Float64 - min_solves::Int - detect_uniform_M::Bool - lengthscales::Vector{Float64} - jitter::Float64 - n_seeds::Int - seeds::Any + sampler::Any subproblem_indicators::Vector{LogicalVariableRef{M}} # Cached submodels: indicator => GDPSubmodel. # Typed Any so extensions can store different types. @@ -475,15 +500,7 @@ mutable struct _MBM{O, T, M <: JuMP.AbstractModel} <: AbstractReformulationMetho method.optimizer, Dict{LogicalVariableRef{M}, Any}(), method.default_M, - method.gp, - method.kappa, - method.budget, - method.min_solves, - method.detect_uniform_M, - method.lengthscales, - method.jitter, - method.n_seeds, - method.seeds, + method.sampler, Vector{LogicalVariableRef{M}}(), Dict{LogicalVariableRef{M}, Any}() ) diff --git a/src/extension_api.jl b/src/extension_api.jl index be7ae48..aa1c517 100644 --- a/src/extension_api.jl +++ b/src/extension_api.jl @@ -40,32 +40,27 @@ Y(t, x) function InfiniteLogical end """ - sample_M_values(gp, objectives, sub, method, support_grids) + sample_M_values(sampler, objectives, sub, method, support_grids) Compute the MBM M values at the transcription supports of an infinite -model. `gp` is the `gp` field of [`MBM`](@ref), `objectives` is the -array of per-support objective expressions, `sub` is the transcribed -submodel wrapped as a `GDPSubmodel`, `method` is the `_MBM` data -(which carries the sampling settings `kappa`, `budget`, `min_solves`, -`detect_uniform_M`, `lengthscales`, `jitter`, `n_seeds`, and `seeds`), -and `support_grids` is a function returning the +model. `sampler` is the `sampler` field of [`MBM`](@ref), +`objectives` is the array of per-support objective expressions, `sub` +is the transcribed submodel wrapped as a `GDPSubmodel`, `method` is +the `_MBM` data, and `support_grids` is a function returning the support vectors of the infinite parameters. It is a function because the supports are only well defined once M is known to vary over them, so methods that return early (or never need coordinates) must not call it. Returns an array of M values shaped like `objectives`, a scalar when M is uniform across the supports, or `nothing` if an M subproblem is infeasible. Extensions implement methods that dispatch -on `gp`: `nothing` solves an M subproblem at every support, while an -`AbstractGPs.AbstractGP` or a `KernelFunctions.Kernel` solves a -subset of the supports selected by an upper-confidence-bound -acquisition and fills the rest with the posterior upper confidence -bound `mean + kappa * sd`. The filled values are heuristic upper -estimates of the exact M values, not certificates. +on `sampler`: `nothing` solves an M subproblem at every support, +while a [`GPSampler`](@ref) solves a subset of the supports and +fills the rest with a Gaussian-process upper confidence bound. """ -function sample_M_values(gp, objectives, sub, method, support_grids) - error("Unrecognized `gp` value `$(repr(gp))` for MBM on an " * - "infinite model. Use `nothing` to solve an M subproblem " * - "at every support, or load AbstractGPs and pass an " * - "`AbstractGPs.AbstractGP` or a `KernelFunctions.Kernel` " * - "to estimate M values with a Gaussian process.") +function sample_M_values(sampler, objectives, sub, method, support_grids) + error("Unrecognized `sampler` value `$(repr(sampler))` for MBM " * + "on an infinite model. Use `nothing` to solve an M " * + "subproblem at every support, or a `GPSampler` (with " * + "AbstractGPs loaded) to estimate M values with a " * + "Gaussian process.") end diff --git a/src/mbm.jl b/src/mbm.jl index ef22f66..8f74269 100644 --- a/src/mbm.jl +++ b/src/mbm.jl @@ -88,12 +88,8 @@ function reformulate_disjunct_constraint( method::_MBM ) ref_cons = reformulate_disjunction(model, con, MBM( - method.optimizer, method.default_M, gp = method.gp, - kappa = method.kappa, budget = method.budget, - min_solves = method.min_solves, - detect_uniform_M = method.detect_uniform_M, - lengthscales = method.lengthscales, jitter = method.jitter, - n_seeds = method.n_seeds, seeds = method.seeds)) + method.optimizer, method.default_M, + sampler = method.sampler)) new_ref_cons = Vector{JuMP.AbstractConstraint}() for ref_con in ref_cons append!(new_ref_cons, diff --git a/test/extensions/AbstractGPsDisjunctiveProgramming.jl b/test/extensions/AbstractGPsDisjunctiveProgramming.jl index dd9ffa0..2f0e818 100644 --- a/test/extensions/AbstractGPsDisjunctiveProgramming.jl +++ b/test/extensions/AbstractGPsDisjunctiveProgramming.jl @@ -1,44 +1,38 @@ using InfiniteOpt, HiGHS, AbstractGPs import DisjunctiveProgramming as DP -function test_gp_mbm_kwargs() - method = MBM(HiGHS.Optimizer) - @test method.gp === nothing - @test method.kappa == 2.5 - @test method.budget == 0.25 - @test method.min_solves == 6 - @test method.detect_uniform_M - @test method.lengthscales == [0.05, 0.1, 0.2, 0.4, 0.8] - @test method.jitter == 1e-8 - @test method.n_seeds == 4 - @test method.seeds === nothing - kern = with_lengthscale(SqExponentialKernel(), 0.3) - method = MBM(HiGHS.Optimizer, gp = kern, kappa = 4.0, - budget = 0.1, min_solves = 3, detect_uniform_M = false, - lengthscales = (0.1, 0.3), jitter = 1e-6, n_seeds = 6, - seeds = [0.0, 0.3, 1.0]) - @test method.gp === kern - @test method.kappa == 4.0 - @test method.budget == 0.1 - @test method.min_solves == 3 - @test !method.detect_uniform_M - @test method.lengthscales == [0.1, 0.3] - @test method.jitter == 1e-6 - @test method.n_seeds == 6 - @test method.seeds == [0.0, 0.3, 1.0] - @test_throws ErrorException MBM(HiGHS.Optimizer, kappa = -1) - @test_throws ErrorException MBM(HiGHS.Optimizer, budget = 0) - @test_throws ErrorException MBM(HiGHS.Optimizer, budget = 1.5) - @test_throws ErrorException MBM(HiGHS.Optimizer, min_solves = 0) - @test_throws ErrorException MBM(HiGHS.Optimizer, - lengthscales = Float64[]) - @test_throws ErrorException MBM(HiGHS.Optimizer, - lengthscales = [-0.1]) - @test_throws ErrorException MBM(HiGHS.Optimizer, jitter = -1) - @test_throws ErrorException MBM(HiGHS.Optimizer, n_seeds = 1) - @test_throws ErrorException MBM(HiGHS.Optimizer, seeds = [1.5]) - @test_throws ErrorException MBM(HiGHS.Optimizer, - seeds = Float64[]) +function test_gp_sampler_kwargs() + @test MBM(HiGHS.Optimizer).sampler === nothing + sampler = GPSampler() + @test sampler.kernel === nothing + @test sampler.kappa == 2.5 + @test sampler.budget == 0.25 + @test sampler.detect_uniform_M + @test sampler.lengthscales == [0.05, 0.1, 0.2, 0.4, 0.8] + @test sampler.jitter == 1e-8 + @test sampler.seeds == 4 + kern = SqExponentialKernel() + sampler = GPSampler(kern, kappa = 4.0, budget = 0.1, + detect_uniform_M = false, lengthscales = (0.1, 0.3), + jitter = 1e-6, seeds = [0.0, 0.3, 1.0]) + @test sampler.kernel === kern + @test sampler.kappa == 4.0 + @test sampler.budget == 0.1 + @test !sampler.detect_uniform_M + @test sampler.lengthscales == [0.1, 0.3] + @test sampler.jitter == 1e-6 + @test sampler.seeds == [0.0, 0.3, 1.0] + @test GPSampler(seeds = 6).seeds == 6 + @test MBM(HiGHS.Optimizer, sampler = sampler).sampler === sampler + @test_throws ErrorException GPSampler(kappa = -1) + @test_throws ErrorException GPSampler(budget = 0) + @test_throws ErrorException GPSampler(budget = 1.5) + @test_throws ErrorException GPSampler(lengthscales = Float64[]) + @test_throws ErrorException GPSampler(lengthscales = [-0.1]) + @test_throws ErrorException GPSampler(jitter = -1) + @test_throws ErrorException GPSampler(seeds = 1) + @test_throws ErrorException GPSampler(seeds = [1.5]) + @test_throws ErrorException GPSampler(seeds = Float64[]) end # Mirror of test_raw_M_infinite_scalar: uniform seed M values collapse @@ -52,7 +46,7 @@ function test_gp_raw_M_scalar() @constraint(model, con2, x <= 3, Disjunct(Y[2])) @disjunction(model, Y) mbm = DP._MBM( - MBM(HiGHS.Optimizer, gp = SqExponentialKernel()), model) + MBM(HiGHS.Optimizer, sampler = GPSampler()), model) sub = DP.copy_model_with_constraints( model, DP.DisjunctConstraintRef[con2], mbm) obj = DP.prepare_max_M_objective( @@ -60,11 +54,10 @@ function test_gp_raw_M_scalar() @test DP.raw_M(sub, obj, mbm) == 5.0 end -# With few supports the budget floor covers every support, so the GP -# sampler solves all of them exactly and must reproduce the exact -# grid parameter function +# With budget = 1.0 every support is solved exactly, so the GP +# sampler must reproduce the exact grid parameter function function test_gp_raw_M_matches_exact() - function pfunc_values(gp, supports) + function pfunc_values(sampler, supports) model = InfiniteGDPModel() @infinite_parameter(model, t ∈ [0, 1], supports = supports) @variable(model, 0 <= x <= 10, Infinite(t)) @@ -73,7 +66,8 @@ function test_gp_raw_M_matches_exact() @constraint(model, con, x <= f, Disjunct(Y[1])) @constraint(model, con2, x >= 0.5, Disjunct(Y[2])) @disjunction(model, Y) - mbm = DP._MBM(MBM(HiGHS.Optimizer, gp = gp), model) + mbm = DP._MBM( + MBM(HiGHS.Optimizer, sampler = sampler), model) sub = DP.copy_model_with_constraints( model, DP.DisjunctConstraintRef[con2], mbm) obj = DP.prepare_max_M_objective( @@ -84,11 +78,15 @@ function test_gp_raw_M_matches_exact() end supports = [0.0, 0.25, 0.5, 0.75, 1.0] exact_vals = pfunc_values(nothing, supports) - # a kernel gets its lengthscale fit; a GP is used as given. Both - # solve the same supports here, so the M values match exactly - @test pfunc_values(SqExponentialKernel(), supports) == exact_vals - gp = GP(with_lengthscale(SqExponentialKernel(), 0.2)) - @test pfunc_values(gp, supports) == exact_vals + # default kernel, user kernel, and pinned lengthscale all solve + # the same supports here, so the M values match exactly + @test pfunc_values(GPSampler(budget = 1.0), supports) == + exact_vals + @test pfunc_values( + GPSampler(SqExponentialKernel(), budget = 1.0), supports) == + exact_vals + @test pfunc_values(GPSampler(SqExponentialKernel(), + lengthscales = [0.2], budget = 1.0), supports) == exact_vals end # an empty disjunct region makes the M subproblems infeasible; both @@ -108,16 +106,16 @@ function test_gp_infeasible_disjunct() @objective(model, Max, 𝔼(x, t)) return model end - for gp in (nothing, SqExponentialKernel()) + for sampler in (nothing, GPSampler()) model = build() @test_throws ErrorException optimize!(model, - gdp_method = MBM(HiGHS.Optimizer, gp = gp)) + gdp_method = MBM(HiGHS.Optimizer, sampler = sampler)) end end # optimum (10) needs M(t) >= 10 - 2t pointwise; the GP fill is heuristic function test_gp_mbm_solve_equivalence() - function solve_with(method) + function solve_with(sampler) model = InfiniteGDPModel(HiGHS.Optimizer) set_silent(model) @infinite_parameter(model, t ∈ [0, 1], num_supports = 20) @@ -128,15 +126,14 @@ function test_gp_mbm_solve_equivalence() @constraint(model, x >= 0.5, Disjunct(Y[2])) @disjunction(model, Y) @objective(model, Max, 𝔼(x, t)) - optimize!(model, gdp_method = method) + optimize!(model, + gdp_method = MBM(HiGHS.Optimizer, sampler = sampler)) @test termination_status(model) == MOI.OPTIMAL return objective_value(model) end - obj_exact = solve_with(MBM(HiGHS.Optimizer)) - obj_gp = solve_with( - MBM(HiGHS.Optimizer, gp = SqExponentialKernel())) - obj_tuned = solve_with(MBM(HiGHS.Optimizer, - gp = SqExponentialKernel(), kappa = 4.0, budget = 0.2)) + obj_exact = solve_with(nothing) + obj_gp = solve_with(GPSampler()) + obj_tuned = solve_with(GPSampler(kappa = 4.0, budget = 0.2)) @test obj_exact ≈ 10.0 atol = 1e-4 # over-M can't raise the optimum, under-M can only shave it a bit @test obj_gp <= obj_exact + 1e-6 @@ -153,7 +150,7 @@ end # exact. function test_gp_periodic_M_seeds() supports = [0.0, 0.25, 0.5, 0.75, 1.0] - function solve_with(; kwargs...) + function solve_with(sampler) model = InfiniteGDPModel(HiGHS.Optimizer) set_silent(model) @infinite_parameter(model, t ∈ [0, 1], supports = supports) @@ -165,15 +162,14 @@ function test_gp_periodic_M_seeds() @disjunction(model, Y) @objective(model, Max, 𝔼(x, t)) optimize!(model, - gdp_method = MBM(HiGHS.Optimizer; kwargs...)) + gdp_method = MBM(HiGHS.Optimizer, sampler = sampler)) return objective_value(model) end - @test solve_with() ≈ 10.0 atol = 1e-6 - @test solve_with(gp = SqExponentialKernel()) ≈ 10.0 atol = 1e-6 - @test solve_with(gp = SqExponentialKernel(), n_seeds = 5) ≈ - 10.0 atol = 1e-6 - @test solve_with(gp = SqExponentialKernel(), - seeds = [0.0, 0.5, 1.0]) ≈ 9.0 atol = 1e-6 + @test solve_with(nothing) ≈ 10.0 atol = 1e-6 + @test solve_with(GPSampler()) ≈ 10.0 atol = 1e-6 + @test solve_with(GPSampler(seeds = 5)) ≈ 10.0 atol = 1e-6 + @test solve_with(GPSampler(seeds = [0.0, 0.5, 1.0])) ≈ + 9.0 atol = 1e-6 end # With detection off the uniform M is not collapsed to a scalar: the @@ -189,8 +185,7 @@ function test_gp_detect_uniform_M_off() @constraint(model, con2, x <= 3, Disjunct(Y[2])) @disjunction(model, Y) mbm = DP._MBM(MBM(HiGHS.Optimizer, - gp = SqExponentialKernel(), - detect_uniform_M = detect), model) + sampler = GPSampler(detect_uniform_M = detect)), model) sub = DP.copy_model_with_constraints( model, DP.DisjunctConstraintRef[con2], mbm) obj = DP.prepare_max_M_objective( @@ -217,8 +212,7 @@ function test_gp_detect_uniform_M_off_dependent() @constraint(model, con2, x <= 3, Disjunct(Y[2])) @disjunction(model, Y) mbm = DP._MBM(MBM(HiGHS.Optimizer, - gp = SqExponentialKernel(), - detect_uniform_M = false), model) + sampler = GPSampler(detect_uniform_M = false)), model) sub = DP.copy_model_with_constraints( model, DP.DisjunctConstraintRef[con2], mbm) obj = DP.prepare_max_M_objective( @@ -226,7 +220,7 @@ function test_gp_detect_uniform_M_off_dependent() @test_throws ErrorException DP.raw_M(sub, obj, mbm) end -function test_gp_unknown_gp_error() +function test_gp_unknown_sampler_error() model = InfiniteGDPModel(HiGHS.Optimizer) set_silent(model) @infinite_parameter(model, t ∈ [0, 1], num_supports = 5) @@ -238,17 +232,17 @@ function test_gp_unknown_gp_error() @disjunction(model, Y) @objective(model, Max, 𝔼(x, t)) @test_throws ErrorException optimize!(model, - gdp_method = MBM(HiGHS.Optimizer, gp = :grid)) + gdp_method = MBM(HiGHS.Optimizer, sampler = :grid)) end @testset "AbstractGPsDisjunctiveProgramming" begin - test_gp_mbm_kwargs() + test_gp_sampler_kwargs() test_gp_raw_M_scalar() test_gp_raw_M_matches_exact() test_gp_mbm_solve_equivalence() test_gp_periodic_M_seeds() test_gp_detect_uniform_M_off() test_gp_detect_uniform_M_off_dependent() - test_gp_unknown_gp_error() + test_gp_unknown_sampler_error() test_gp_infeasible_disjunct() end From 081e722edb40b3c0126edd289350916a1d2ae839 Mon Sep 17 00:00:00 2001 From: d227nguyen Date: Wed, 19 Aug 2026 13:39:58 -0400 Subject: [PATCH 10/10] datatypes and variable name updates --- ext/AbstractGPsDisjunctiveProgramming.jl | 107 +++++++++++++---------- ext/InfiniteDisjunctiveProgramming.jl | 33 ++++--- 2 files changed, 81 insertions(+), 59 deletions(-) diff --git a/ext/AbstractGPsDisjunctiveProgramming.jl b/ext/AbstractGPsDisjunctiveProgramming.jl index 09e3d8e..00af68c 100644 --- a/ext/AbstractGPsDisjunctiveProgramming.jl +++ b/ext/AbstractGPsDisjunctiveProgramming.jl @@ -8,41 +8,56 @@ import DisjunctiveProgramming as DP # GP FITTING ################################################################################ # Normalized to [0, 1]^d so one lengthscale works across dimensions -function _support_coords(grids) - idxs = CartesianIndices(length.(grids)) - los = [minimum(g) for g in grids] - rng = [max(maximum(g) - minimum(g), eps()) for g in grids] - return [[(grids[d][I[d]] - los[d]) / rng[d] for d in 1:length(grids)] - for I in vec(idxs)] +function _support_coords( + grids::Tuple{Vararg{Vector{Float64}}} + )::Vector{Vector{Float64}} + indices = CartesianIndices(length.(grids)) + mins = [minimum(g) for g in grids] + ranges = [max(maximum(g) - minimum(g), eps()) for g in grids] + return [[(grids[d][I[d]] - mins[d]) / ranges[d] + for d in 1:length(grids)] for I in vec(indices)] end # Lengthscale selected by marginal likelihood over the candidates -function _fit_posterior(kernel, X, y, sampler) - best_post, best_lp = nothing, -Inf - for ls in sampler.lengthscales - kern = KernelFunctions.with_lengthscale(kernel, ls) - fx = AbstractGPs.GP(kern)(X, sampler.jitter) - lp = AbstractGPs.logpdf(fx, y) - if lp > best_lp - best_post, best_lp = AbstractGPs.posterior(fx, y), lp +function _fit_posterior( + kernel::KernelFunctions.Kernel, + X::Vector{Vector{Float64}}, + y::Vector{Float64}, + sampler::DP.GPSampler + ) + best_posterior, best_log_prob = nothing, -Inf + for lengthscale in sampler.lengthscales + scaled_kernel = KernelFunctions.with_lengthscale( + kernel, lengthscale) + finite_gp = AbstractGPs.GP(scaled_kernel)(X, sampler.jitter) + log_prob = AbstractGPs.logpdf(finite_gp, y) + if log_prob > best_log_prob + best_posterior = AbstractGPs.posterior(finite_gp, y) + best_log_prob = log_prob end end - return best_post + return best_posterior end -function _mean_sd(sampler, X, solved) - lis = collect(keys(solved)) - y = [solved[li] for li in lis] - ybar = sum(y) / length(y) +function _mean_sd( + sampler::DP.GPSampler, + X::Vector{Vector{Float64}}, + solved::Dict{Int, Float64} + ) + solved_indices = collect(keys(solved)) + y = [solved[i] for i in solved_indices] + y_mean = sum(y) / length(y) # floored so near-equal solved values still cushion the filled ones - ystd = max(sqrt(sum(abs2, y .- ybar) / max(length(y) - 1, 1)), - 1e-2 * abs(ybar), 1e-8) + y_scale = max(sqrt(sum(abs2, y .- y_mean) / max(length(y) - 1, 1)), + 1e-2 * abs(y_mean), 1e-8) kernel = something(sampler.kernel, KernelFunctions.SqExponentialKernel()) - post = _fit_posterior(kernel, X[lis], (y .- ybar) ./ ystd, sampler) - mz = AbstractGPs.mean(post, X) - vz = max.(AbstractGPs.var(post, X), 0.0) - return mz .* ystd .+ ybar, sqrt.(vz) .* ystd + posterior = _fit_posterior( + kernel, X[solved_indices], (y .- y_mean) ./ y_scale, sampler) + posterior_mean = AbstractGPs.mean(posterior, X) + posterior_var = max.(AbstractGPs.var(posterior, X), 0.0) + return posterior_mean .* y_scale .+ y_mean, + sqrt.(posterior_var) .* y_scale end ################################################################################ @@ -54,22 +69,22 @@ function DP.sample_M_values( objectives::AbstractArray, sub::DP.GDPSubmodel, method::DP._MBM, - support_grids + support_grids::Function ) - idxs = collect(CartesianIndices(objectives)) - n = length(idxs) + indices = collect(CartesianIndices(objectives)) + n = length(indices) solved = Dict{Int, Float64}() - solve_at(li) = begin - m = DP.raw_M(sub, objectives[idxs[li]], method) - m === nothing && return false - solved[li] = m + solve_at(index::Int) = begin + M_val = DP.raw_M(sub, objectives[indices[index]], method) + M_val === nothing && return false + solved[index] = M_val return true end # an evenly spaced seed count, or user-given seed fractions - fracs = sampler.seeds isa Int ? + fractions = sampler.seeds isa Int ? range(0, 1, length = sampler.seeds) : sampler.seeds - for s in unique(1 .+ round.(Int, fracs .* (n - 1))) - solve_at(s) || return nothing + for index in unique(1 .+ round.(Int, fractions .* (n - 1))) + solve_at(index) || return nothing end if sampler.detect_uniform_M # a uniform M needs no fit, and so no support grid either @@ -79,24 +94,24 @@ function DP.sample_M_values( budget = min(ceil(Int, sampler.budget * n), n) X = _support_coords(support_grids()) while length(solved) < budget - ms, ss = _mean_sd(sampler, X, solved) - acq = ms .+ sampler.kappa .* ss - for li in keys(solved) - acq[li] = -Inf + means, sds = _mean_sd(sampler, X, solved) + acquisition = means .+ sampler.kappa .* sds + for index in keys(solved) + acquisition[index] = -Inf end - solve_at(argmax(acq)) || return nothing + solve_at(argmax(acquisition)) || return nothing end M_vals = Array{Float64}(undef, size(objectives)) if length(solved) == n # nothing left to estimate - for (li, I) in enumerate(idxs) - M_vals[I] = solved[li] + for (index, I) in enumerate(indices) + M_vals[I] = solved[index] end return M_vals end - ms, ss = _mean_sd(sampler, X, solved) - for (li, I) in enumerate(idxs) # exact M values are nonnegative - M_vals[I] = get(solved, li, - max(ms[li] + sampler.kappa * ss[li], 0.0)) + means, sds = _mean_sd(sampler, X, solved) + for (index, I) in enumerate(indices) # exact M values are nonnegative + M_vals[I] = get(solved, index, + max(means[index] + sampler.kappa * sds[index], 0.0)) end return M_vals end diff --git a/ext/InfiniteDisjunctiveProgramming.jl b/ext/InfiniteDisjunctiveProgramming.jl index 2a492e6..83aed38 100644 --- a/ext/InfiniteDisjunctiveProgramming.jl +++ b/ext/InfiniteDisjunctiveProgramming.jl @@ -41,20 +41,22 @@ end # Bound info for parameter refs in disjunct constraints: parameter # functions report their support extrema, finite parameters a point, # other parameters (-Inf, Inf). Returns nothing for variables. -function _parameter_bound_info(vref::InfiniteOpt.GeneralVariableRef) +function _parameter_bound_info( + vref::InfiniteOpt.GeneralVariableRef + )::Union{Nothing, Tuple{Float64, Float64}} _is_parameter(vref) || return nothing dvref = InfiniteOpt.dispatch_variable_ref(vref) if dvref isa InfiniteOpt.ParameterFunctionRef prefs = InfiniteOpt.parameter_list(dvref) length(prefs) == 1 || return (-Inf, Inf) - supps = InfiniteOpt.supports(first(prefs)) - isempty(supps) && return (-Inf, Inf) - f = InfiniteOpt.raw_function(dvref) - vals = [f(s) for s in supps] - return (minimum(vals), maximum(vals)) + supports = InfiniteOpt.supports(first(prefs)) + isempty(supports) && return (-Inf, Inf) + func = InfiniteOpt.raw_function(dvref) + func_values = [func(s) for s in supports] + return (minimum(func_values), maximum(func_values)) elseif dvref isa InfiniteOpt.FiniteParameterRef - v = InfiniteOpt.parameter_value(dvref) - return (v, v) + value = InfiniteOpt.parameter_value(dvref) + return (value, value) end return (-Inf, Inf) end @@ -69,7 +71,11 @@ function DP.set_variable_bound_info( end # Hull and PSplit require finite bounds that include 0 -function _clamped_bound_info(vref, info, method_name) +function _zero_inclusive_bounds( + vref::InfiniteOpt.GeneralVariableRef, + info::Union{Nothing, Tuple{Float64, Float64}}, + method_name::String + )::Tuple{Float64, Float64} info === nothing || return (min(0, info[1]), max(0, info[2])) if !JuMP.has_lower_bound(vref) || !JuMP.has_upper_bound(vref) error("Variable $vref must have both lower and upper " * @@ -82,14 +88,14 @@ end function DP.set_variable_bound_info( vref::InfiniteOpt.GeneralVariableRef, ::DP.Hull) - return _clamped_bound_info(vref, + return _zero_inclusive_bounds(vref, _parameter_bound_info(vref), "Hull") end function DP.set_variable_bound_info( vref::InfiniteOpt.GeneralVariableRef, ::Union{DP.PSplit, DP._PSplit}) - return _clamped_bound_info(vref, + return _zero_inclusive_bounds(vref, _parameter_bound_info(vref), "PSplit") end @@ -321,7 +327,8 @@ end # The infinite parameters of `mini_expr` and their supports, in the # ascending order of `parameter_refs`. Only defined when M varies over # the supports, so it is deferred until sample_M_values needs it. -function _support_grids(sub::DP.GDPSubmodel, mini_expr) +function _support_grids( + sub::DP.GDPSubmodel, mini_expr::JuMP.AbstractJuMPScalar) reverse_map = Dict(ws[1] => v for (v, ws) in sub.fwd_map) prefs = Tuple(get(reverse_map, p) do error("MBM cannot build a support grid over `$p`, which " * @@ -337,7 +344,7 @@ function DP.sample_M_values( objectives::AbstractArray, sub::DP.GDPSubmodel, method::DP._MBM, - support_grids + support_grids::Function ) M_vals = Array{Float64}(undef, size(objectives)) for I in eachindex(objectives)