diff --git a/docs/make.jl b/docs/make.jl index 78432e32..cf0c4bf3 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -17,10 +17,7 @@ Documenter.makedocs(; ), sitename = "ParametricOptInterface.jl", authors = "Tomás Gutierrez, and contributors", - pages = [ - "Home" => "index.md", - "reference.md", - ], + pages = ["Home" => "index.md", "reference.md"], checkdocs = :none, ) diff --git a/docs/src/index.md b/docs/src/index.md index bbf3d486..2cc7668f 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -283,3 +283,29 @@ set_attribute(model, POI.ConstraintsInterpretation(), POI.BOUNDS_AND_CONSTRAINTS @constraint(model, x >= p) optimize!(model) ``` + +## Parameters multiplying quadratic terms + +POI supports parameters that multiply quadratic variable terms in objectives +**only**. This creates cubic polynomial expressions of the form `c * p * x * y` +where `c` is a number, `p` is a parameter, and `x` and `y` are variables. After +parameter substitution, the objective is quadratic instead of cubic. + +Note that the maximum degree is 3 (cubic), at least one factor in each cubic +term must be a parameter, and pure cubic variable terms (for example, +`x * y * z` with no parameters) are not supported. + +```@repl +using JuMP, HiGHS +import ParametricOptInterface as POI +model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) +set_silent(model) +@variable(model, 0 <= x <= 10) +@variable(model, p in Parameter(2)) +@objective(model, Min, p * x^2 - 3x) +optimize!(model) +value(x) # x = 3 / 2p = 0.75 +set_parameter_value(p, 3) +optimize!(model) +value(x) # x = 3 / 2p = 0.5 +``` diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 6352876d..f6c25882 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -436,7 +436,7 @@ function MOI.delete(model::Optimizer, v::MOI.VariableIndex) MOI.delete(model.optimizer, v) MOI.delete(model.original_objective_cache, v) # TODO - what happens if the variable was in a SAF that was converted to bounds? - # solution: do not allow if that is the case (requires going trhought the scalar affine cache) + # solution: do not allow if that is the case (requires going through the scalar affine cache) # TODO - deleting a variable also deletes constraints for (F, S) in MOI.Utilities.DoubleDicts.nonempty_outer_keys( model.constraint_outer_to_inner, @@ -1258,6 +1258,7 @@ end function _empty_objective_function_caches!(model::Optimizer{T}) where {T} model.affine_objective_cache = nothing model.quadratic_objective_cache = nothing + model.cubic_objective_cache = nothing model.original_objective_cache = MOI.Utilities.ObjectiveContainer{T}() return end diff --git a/src/ParametricOptInterface.jl b/src/ParametricOptInterface.jl index 8ebe027f..d341adb9 100644 --- a/src/ParametricOptInterface.jl +++ b/src/ParametricOptInterface.jl @@ -74,11 +74,19 @@ const VariableMap = MOI.Utilities.CleverDicts.CleverDict{ const DoubleDict{T} = MOI.Utilities.DoubleDicts.DoubleDict{T} const DoubleDictInner{F,S,T} = MOI.Utilities.DoubleDicts.DoubleDictInner{F,S,T} +# +# cubic functions helpers +# + +include("cubic_types.jl") +include("cubic_parser.jl") + # # parametric functions # include("parametric_functions.jl") +include("parametric_cubic_function.jl") """ Optimizer{T, OT <: MOI.ModelLike} <: MOI.AbstractOptimizer @@ -151,6 +159,7 @@ mutable struct Optimizer{T,OT<:MOI.ModelLike} <: MOI.AbstractOptimizer # Clever cache of data (at most one can be !== nothing) affine_objective_cache::Union{Nothing,ParametricAffineFunction{T}} quadratic_objective_cache::Union{Nothing,ParametricQuadraticFunction{T}} + cubic_objective_cache::Union{Nothing,ParametricCubicFunction{T}} original_objective_cache::MOI.Utilities.ObjectiveContainer{T} # Store parametric expressions for product of variables quadratic_objective_cache_product::Dict{ @@ -226,7 +235,7 @@ mutable struct Optimizer{T,OT<:MOI.ModelLike} <: MOI.AbstractOptimizer # objective nothing, nothing, - # nothing, + nothing, # cubic_objective_cache MOI.Utilities.ObjectiveContainer{T}(), Dict{ Tuple{MOI.VariableIndex,MOI.VariableIndex}, @@ -275,5 +284,6 @@ end include("duals.jl") include("update_parameters.jl") include("MOI_wrapper.jl") +include("cubic_objective.jl") end # module diff --git a/src/cubic_objective.jl b/src/cubic_objective.jl new file mode 100644 index 00000000..15bd0534 --- /dev/null +++ b/src/cubic_objective.jl @@ -0,0 +1,178 @@ +# Copyright (c) 2020: Tomás Gutierrez and contributors +# +# Use of this source code is governed by an MIT-style license that can be found +# in the LICENSE.md file or at https://opensource.org/licenses/MIT. + +function MOI.set( + model::Optimizer{T}, + ::MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}, + f::MOI.ScalarNonlinearFunction, +) where {T} + # 1. Attempt to parse as cubic + parsed = _parse_cubic_expression(f, T) + if parsed === nothing + error( + "ScalarNonlinearFunction must be a valid cubic polynomial with " * + "parameters multiplying at most quadratic variable terms. " * + "Non-polynomial operations or degree > 3 are not supported.", + ) + end + + # 2. Create ParametricCubicFunction + cubic_func = ParametricCubicFunction(parsed) + + # 3. Compute current function for inner optimizer + current = _current_function(cubic_func, model) + + # 4. Set current function on inner optimizer + try + MOI.set( + model.optimizer, + MOI.ObjectiveFunction{typeof(current)}(), + current, + ) + catch e + # rethrow the original error with the additional info of the objective function that caused it + error( + "Failed to set cubic objective function, f = $f, on inner " * + "optimizer. " * + "This may be due to unsupported features in the cubic " * + "expression. " * + "Original error: $(e.msg)", + ) + end + + # 5. Clear old caches + _empty_objective_function_caches!(model) + + # 6. Store new cache + model.cubic_objective_cache = cubic_func + + # 7. Store original for retrieval if option is enabled + if model.save_original_objective_and_constraints + MOI.set( + model.original_objective_cache, + MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}(), + f, + ) + end + + return nothing +end + +function MOI.get( + model::Optimizer, + attr::MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}, +) + if model.cubic_objective_cache === nothing + error("No ScalarNonlinearFunction objective is set") + end + if !model.save_original_objective_and_constraints + error( + "Cannot retrieve original objective: save_original_objective_and_constraints is false", + ) + end + return MOI.get(model.original_objective_cache, attr) +end + +function MOI.supports( + ::Optimizer, + ::MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}, +) + return true +end + +""" + _update_cubic_objective!(model::Optimizer{T}) where {T} + +Update the cubic objective after parameters have changed. +Uses incremental modifications (ScalarQuadraticCoefficientChange, ScalarCoefficientChange, +ScalarConstantChange) for efficiency when the solver supports them. +Falls back to rebuilding the full objective if incremental modifications are not supported. +""" +function _update_cubic_objective!(model::Optimizer{T}) where {T} + if model.cubic_objective_cache === nothing + return + end + pf = model.cubic_objective_cache + + # Check if any changes are needed by computing deltas + delta_constant = _delta_parametric_constant(model, pf) + delta_affine = _delta_parametric_affine_terms(model, pf) + delta_quadratic = _delta_parametric_quadratic_terms(model, pf) + + if iszero(delta_constant) && + isempty(delta_affine) && + isempty(delta_quadratic) + return # No changes needed + end + + _try_incremental_cubic_update!( + model, + pf, + delta_constant, + delta_affine, + delta_quadratic, + ) + + return nothing +end + +""" + _try_incremental_cubic_update!(model, pf, delta_constant, delta_affine, delta_quadratic) + +Apply incremental coefficient updates to the inner optimizer's objective. +""" +function _try_incremental_cubic_update!( + model::Optimizer{T}, + pf::ParametricCubicFunction{T}, + delta_constant::T, + delta_affine::Dict{MOI.VariableIndex,T}, + delta_quadratic::Dict{Tuple{MOI.VariableIndex,MOI.VariableIndex},T}, +) where {T} + # Get the current objective function type from the inner optimizer + F = MOI.get(model.optimizer, MOI.ObjectiveFunctionType()) + + # Compute full new values (not deltas) for robustness + # The delta was used to detect changes; now apply full new coefficients + new_quad_terms = _parametric_quadratic_terms(model, pf) + new_affine_terms = _parametric_affine_terms(model, pf) + new_constant = _parametric_constant(model, pf) + + # Apply quadratic coefficient changes + # MOI convention: + # - Off-diagonal (v1 != v2): coefficient C means C*v1*v2 (use as-is) + # - Diagonal (v1 == v2): coefficient C means (C/2)*v1^2 (multiply by 2) + for ((var1, var2), _) in delta_quadratic + new_coef = new_quad_terms[(var1, var2)] + # Apply MOI coefficient convention + moi_coef = var1 == var2 ? new_coef * 2 : new_coef + MOI.modify( + model.optimizer, + MOI.ObjectiveFunction{F}(), + MOI.ScalarQuadraticCoefficientChange(var1, var2, moi_coef), + ) + end + + # Apply affine coefficient changes (use full new coefficient) + for (var, _) in delta_affine + new_coef = new_affine_terms[var] + MOI.modify( + model.optimizer, + MOI.ObjectiveFunction{F}(), + MOI.ScalarCoefficientChange(var, new_coef), + ) + end + + # Apply constant change + if !iszero(delta_constant) + pf.current_constant = new_constant + MOI.modify( + model.optimizer, + MOI.ObjectiveFunction{F}(), + MOI.ScalarConstantChange(pf.current_constant), + ) + end + + return nothing +end diff --git a/src/cubic_parser.jl b/src/cubic_parser.jl new file mode 100644 index 00000000..e18f58dd --- /dev/null +++ b/src/cubic_parser.jl @@ -0,0 +1,492 @@ +# Copyright (c) 2020: Tomás Gutierrez and contributors +# +# Use of this source code is governed by an MIT-style license that can be found +# in the LICENSE.md file or at https://opensource.org/licenses/MIT. + +""" + _Monomial{T} + +Intermediate representation of a monomial during parsing. +""" +struct _Monomial{T} + coefficient::T + variables::Vector{MOI.VariableIndex} # includes both vars and params +end + +function _Monomial{T}(coefficient::T) where {T} + return _Monomial{T}(coefficient, MOI.VariableIndex[]) +end + +function _Monomial{T}(coefficient::T, var::MOI.VariableIndex) where {T} + return _Monomial{T}(coefficient, [var]) +end + +""" + _monomial_degree(m::_Monomial) -> Int + +Total degree of a monomial (number of variable/parameter factors). +""" +function _monomial_degree(m::_Monomial) + return length(m.variables) +end + +""" + _multiply_monomials(m1::_Monomial{T}, m2::_Monomial{T}) where {T} + +Multiply two monomials together. +""" +function _multiply_monomials(m1::_Monomial{T}, m2::_Monomial{T}) where {T} + return _Monomial{T}( + m1.coefficient * m2.coefficient, + vcat(m1.variables, m2.variables), + ) +end + +""" + _scale_monomial(m::_Monomial{T}, scalar::T) where {T} + +Scale a monomial by a scalar. +""" +function _scale_monomial(m::_Monomial{T}, scalar::T) where {T} + return _Monomial{T}(m.coefficient * scalar, copy(m.variables)) +end + +""" + _ParsedCubicExpression{T} + +Result of parsing a ScalarNonlinearFunction into cubic polynomial form. +""" +struct _ParsedCubicExpression{T} + pvv::Vector{_ScalarCubicTerm{T}} # Cubic terms with 1 parameter and 2 variables + ppv::Vector{_ScalarCubicTerm{T}} # Cubic terms with 2 parameters and 1 variable + ppp::Vector{_ScalarCubicTerm{T}} # Cubic terms with 3 parameters + + vv::Vector{MOI.ScalarQuadraticTerm{T}} # Quadratic terms with 2 variables + pv::Vector{MOI.ScalarQuadraticTerm{T}} # Quadratic terms with 1 parameter and 1 variable + pp::Vector{MOI.ScalarQuadraticTerm{T}} # Quadratic terms with 2 parameters + + v::Vector{MOI.ScalarAffineTerm{T}} # Affine terms with 1 variable + p::Vector{MOI.ScalarAffineTerm{T}} # Affine terms with 1 parameter + + constant::T # Constant term +end + +""" + _expand_to_monomials(arg, ::Type{T}) where {T} -> Union{Vector{_Monomial{T}}, Nothing} + +Expand an expression argument to a list of monomials. +Returns `nothing` if the expression is not a valid polynomial. +""" +function _expand_to_monomials(arg::Real, ::Type{T}) where {T} + return [_Monomial{T}(T(arg))] +end + +function _expand_to_monomials(arg::MOI.VariableIndex, ::Type{T}) where {T} + return [_Monomial{T}(one(T), arg)] +end + +function _expand_to_monomials( + arg::MOI.ScalarAffineFunction{T}, + ::Type{T}, +) where {T} + monomials = _Monomial{T}[] + for term in arg.terms + push!(monomials, _Monomial{T}(term.coefficient, term.variable)) + end + if !iszero(arg.constant) + push!(monomials, _Monomial{T}(arg.constant)) + end + return monomials +end + +function _expand_to_monomials( + arg::MOI.ScalarQuadraticFunction{T}, + ::Type{T}, +) where {T} + monomials = _Monomial{T}[] + # Quadratic terms + # MOI convention: + # - Off-diagonal (v1 != v2): coefficient C represents C*v1*v2 + # - Diagonal (v1 == v2): coefficient C represents (C/2)*v1^2 + for term in arg.quadratic_terms + coef = term.coefficient + if term.variable_1 == term.variable_2 + coef = coef / 2 # Diagonal: undo MOI's factor of 2 + end + # Off-diagonal: use coefficient as-is + push!(monomials, _Monomial{T}(coef, [term.variable_1, term.variable_2])) + end + # Affine terms + for term in arg.affine_terms + push!(monomials, _Monomial{T}(term.coefficient, term.variable)) + end + # Constant + if !iszero(arg.constant) + push!(monomials, _Monomial{T}(arg.constant)) + end + return monomials +end + +function _expand_to_monomials( + f::MOI.ScalarNonlinearFunction, + ::Type{T}, +) where {T} + head = f.head + args = f.args + + if head == :+ + return _expand_addition(args, T) + elseif head == :- + return _expand_subtraction(args, T) + elseif head == :* + return _expand_multiplication(args, T) + elseif head == :/ + return _expand_division(args, T) + elseif head == :^ + return _expand_power(args, T) + end + + return nothing +end + +""" + _expand_addition(args, ::Type{T}) where {T} + +Expand addition: collect monomials from all arguments. +""" +function _expand_addition(args, ::Type{T}) where {T} + result = _Monomial{T}[] + for arg in args + monomials = _expand_to_monomials(arg, T) + if monomials === nothing + return nothing + end + append!(result, monomials) + end + return result +end + +""" + _expand_subtraction(args, ::Type{T}) where {T} + +Expand subtraction: first arg positive, rest negative. +""" +function _expand_subtraction(args, ::Type{T}) where {T} + result = _Monomial{T}[] + + if length(args) == 1 + # Unary minus + monomials = _expand_to_monomials(args[1], T) + if monomials === nothing + return nothing + end + for m in monomials + push!(result, _scale_monomial(m, -one(T))) + end + else + # Binary subtraction + for (i, arg) in enumerate(args) + monomials = _expand_to_monomials(arg, T) + if monomials === nothing + return nothing + end + if i == 1 + append!(result, monomials) + else + for m in monomials + push!(result, _scale_monomial(m, -one(T))) + end + end + end + end + return result +end + +""" + _expand_multiplication(args, ::Type{T}) where {T} + +Expand multiplication: multiply all arguments together. +""" +function _expand_multiplication(args, ::Type{T}) where {T} + # Start with identity monomial + result = [_Monomial{T}(one(T))] + + for arg in args + monomials = _expand_to_monomials(arg, T) + if monomials === nothing + return nothing + end + + # Multiply each result monomial with each new monomial + new_result = _Monomial{T}[] + for m1 in result + for m2 in monomials + push!(new_result, _multiply_monomials(m1, m2)) + end + end + result = new_result + end + + return result +end + +""" + _expand_division(args, ::Type{T}) where {T} + +Expand division: multiply numerator by the inverse of denominator +""" +function _expand_division(args, ::Type{T}) where {T} + if length(args) != 2 + return nothing + end + + numerator = args[1] + denominator = args[2] + + # denominator must be a nonzero constant (no variables or parameters) + if !(denominator isa Real) || iszero(denominator) + return nothing + end + + return _expand_multiplication([one(T) / denominator, numerator], T) +end + +""" + _expand_power(args, ::Type{T}) where {T} + +Expand power: x^n becomes x*x*...*x (n times). +""" +function _expand_power(args, ::Type{T}) where {T} + if length(args) != 2 + return nothing + end + + base = args[1] + exponent = args[2] + + # Exponent must be a non-negative integer + if !(exponent isa Integer) || exponent < 0 + return nothing + end + + n = Int(exponent) + + if n == 0 + return [_Monomial{T}(one(T))] + end + + base_monomials = _expand_to_monomials(base, T) + if base_monomials === nothing + return nothing + end + + # x^n = x * x * ... * x (n times) + result = base_monomials + for _ in 2:n + new_result = _Monomial{T}[] + for m1 in result + for m2 in base_monomials + push!(new_result, _multiply_monomials(m1, m2)) + end + end + result = new_result + end + + return result +end + +""" + _combine_like_monomials(monomials::Vector{_Monomial{T}}) where {T} + +Combine like monomials (same variables, regardless of order). +""" +function _combine_like_monomials(monomials::Vector{_Monomial{T}}) where {T} + # Use a dict keyed by sorted variable tuple + combined = Dict{Vector{MOI.VariableIndex},T}() + + for m in monomials + # Sort variables for canonical key + key = sort(m.variables, by = v -> v.value) + combined[key] = get(combined, key, zero(T)) + m.coefficient + end + + result = _Monomial{T}[] + for (vars, coef) in combined + if !iszero(coef) + push!(result, _Monomial{T}(coef, vars)) + end + end + + return result +end + +""" + _classify_monomial(m::_Monomial) -> Symbol + +Classify a monomial by its structure. +""" +function _classify_monomial(m::_Monomial) + degree = _monomial_degree(m) + num_params = count(_is_parameter, m.variables) + + if degree == 0 + return :constant + elseif degree == 1 + return num_params == 1 ? :p : :v + elseif degree == 2 + if num_params == 0 + return :vv + elseif num_params == 1 + return :pv + else + return :pp + end + elseif degree == 3 + if num_params == 0 + return :vvv # Invalid - no parameter + elseif num_params == 1 + return :pvv + elseif num_params == 2 + return :ppv + else + return :ppp + end + else + return :invalid # Degree > 3 + end +end + +""" + _parse_cubic_expression(f::MOI.ScalarNonlinearFunction, ::Type{T}) where {T} -> Union{_ParsedCubicExpression{T}, Nothing} + +Parse a ScalarNonlinearFunction and return a _ParsedCubicExpression if it represents +a valid cubic polynomial (with parameters multiplying at most quadratic variable terms). + +Returns `nothing` if the expression: +- Contains non-polynomial operations (sin, exp, etc.) +- Has degree > 3 in any monomial +- Has a cubic term with no parameters (x*y*z) +""" +function _parse_cubic_expression( + f::MOI.ScalarNonlinearFunction, + ::Type{T}, +) where {T} + # Expand to monomials + monomials = _expand_to_monomials(f, T) + if monomials === nothing + return nothing + end + + # Combine like terms + monomials = _combine_like_monomials(monomials) + + # Classify and collect terms + cubic_terms = _ScalarCubicTerm{T}[] + + cubic_ppp = _ScalarCubicTerm{T}[] + cubic_ppv = _ScalarCubicTerm{T}[] + cubic_pvv = _ScalarCubicTerm{T}[] + + quadratic_pp = MOI.ScalarQuadraticTerm{T}[] + quadratic_pv = MOI.ScalarQuadraticTerm{T}[] + quadratic_vv = MOI.ScalarQuadraticTerm{T}[] + + affine_p = MOI.ScalarAffineTerm{T}[] + affine_v = MOI.ScalarAffineTerm{T}[] + + constant = zero(T) + + for m in monomials + classification = _classify_monomial(m) + + if classification == :invalid || classification == :vvv + return nothing # Invalid degree or no parameter in cubic + elseif classification == :constant + constant += m.coefficient + elseif classification == :v + push!( + affine_v, + MOI.ScalarAffineTerm{T}(m.coefficient, m.variables[1]), + ) + elseif classification == :p + push!( + affine_p, + MOI.ScalarAffineTerm{T}(m.coefficient, m.variables[1]), + ) + elseif classification == :pp + p1 = m.variables[1] + p2 = m.variables[2] + divisor = p1 == p2 ? T(2) : T(1) # Diagonal vs off-diagonal + push!( + quadratic_pp, + MOI.ScalarQuadraticTerm{T}(m.coefficient * divisor, p1, p2), + ) + elseif classification == :pv + # Convention: variable_1 = parameter, variable_2 = variable + # This matches the expectation in _parametric_affine_terms and + # _delta_parametric_affine_terms + is_param = _is_parameter(m.variables[1]) + p_idx_v = ifelse(is_param, m.variables[1], m.variables[2]) + v_idx_v = ifelse(is_param, m.variables[2], m.variables[1]) + push!( + quadratic_pv, + MOI.ScalarQuadraticTerm{T}(m.coefficient, p_idx_v, v_idx_v), + ) + elseif classification == :vv + v1 = m.variables[1] + v2 = m.variables[2] + divisor = v1 == v2 ? T(2) : T(1) # Diagonal vs off-diagonal + push!( + quadratic_vv, + MOI.ScalarQuadraticTerm{T}(m.coefficient * divisor, v1, v2), + ) + elseif classification == :ppp + push!( + cubic_ppp, + _make_cubic_term( + m.coefficient, + m.variables[1], + m.variables[2], + m.variables[3], + ), + ) + elseif classification == :ppv + push!( + cubic_ppv, + _make_cubic_term( + m.coefficient, + m.variables[1], + m.variables[2], + m.variables[3], + ), + ) + else # classification == :pvv + push!( + cubic_pvv, + _make_cubic_term( + m.coefficient, + m.variables[1], + m.variables[2], + m.variables[3], + ), + ) + end + end + + return _ParsedCubicExpression{T}( + cubic_pvv, + cubic_ppv, + cubic_ppp, + quadratic_vv, + quadratic_pv, + quadratic_pp, + affine_v, + affine_p, + constant, + ) +end + +# Convenience method with type inference +function _parse_cubic_expression(f::MOI.ScalarNonlinearFunction) + return _parse_cubic_expression(f, Float64) +end diff --git a/src/cubic_types.jl b/src/cubic_types.jl new file mode 100644 index 00000000..6afe667f --- /dev/null +++ b/src/cubic_types.jl @@ -0,0 +1,69 @@ +# Copyright (c) 2020: Tomás Gutierrez and contributors +# +# Use of this source code is governed by an MIT-style license that can be found +# in the LICENSE.md file or at https://opensource.org/licenses/MIT. + +""" + _ScalarCubicTerm{T} + +Represents a cubic term of the form `coefficient * index_1 * index_2 * index_3`. + +Each index can be either a variable (MOI.VariableIndex) or a parameter (encoded as +VariableIndex with value > PARAMETER_INDEX_THRESHOLD). + +# Fields +- `coefficient::T`: The numeric coefficient +- `index_1::MOI.VariableIndex`: First factor (parameter always) +- `index_2::MOI.VariableIndex`: Second factor (variable or parameter) +- `index_3::MOI.VariableIndex`: Third factor (variable or parameter) + +# Convention +Indices are stored in canonical order: +- Parameters come before variables +This ensures `2*p*x*y` and `2*x*p*y` produce the same term. +""" +struct _ScalarCubicTerm{T} + coefficient::T + index_1::MOI.VariableIndex + index_2::MOI.VariableIndex + index_3::MOI.VariableIndex +end + +""" + _normalize_cubic_indices(idx1, idx2, idx3) -> (idx1, idx2, idx3) + +Normalize cubic term indices to canonical order: +- Parameters come before variables +""" +function _normalize_cubic_indices( + idx1::MOI.VariableIndex, + idx2::MOI.VariableIndex, + idx3::MOI.VariableIndex, +) + params = MOI.VariableIndex[] + vars = MOI.VariableIndex[] + for idx in (idx1, idx2, idx3) + if _is_parameter(idx) + push!(params, idx) + else + push!(vars, idx) + end + end + all_indices = vcat(params, vars) + return all_indices[1], all_indices[2], all_indices[3] +end + +""" + _make_cubic_term(coefficient::T, idx1, idx2, idx3) where {T} + +Create a cubic term with normalized index order. +""" +function _make_cubic_term( + coefficient::T, + idx1::MOI.VariableIndex, + idx2::MOI.VariableIndex, + idx3::MOI.VariableIndex, +) where {T} + n1, n2, n3 = _normalize_cubic_indices(idx1, idx2, idx3) + return _ScalarCubicTerm{T}(coefficient, n1, n2, n3) +end diff --git a/src/parametric_cubic_function.jl b/src/parametric_cubic_function.jl new file mode 100644 index 00000000..99578e77 --- /dev/null +++ b/src/parametric_cubic_function.jl @@ -0,0 +1,469 @@ +# Copyright (c) 2020: Tomás Gutierrez and contributors +# +# Use of this source code is governed by an MIT-style license that can be found +# in the LICENSE.md file or at https://opensource.org/licenses/MIT. + +""" + ParametricCubicFunction{T} <: ParametricFunction{T} + +Represents a cubic function where parameters multiply up to quadratic variable terms. + +Supports the general form: + constant + Σ(affine) + Σ(quadratic) + Σ(cubic) + +After parameter substitution, cubic terms become: +- PVV (p*x*y) → quadratic term (c*p_val*x*y) +- PPV (p*q*x) → affine term (c*p_val*q_val*x) +- PPP (p*q*r) → constant (c*p_val*q_val*r_val) +""" +mutable struct ParametricCubicFunction{T} <: ParametricFunction{T} + # === Cubic terms (degree 3) - split by type like quadratic terms === + pvv::Vector{_ScalarCubicTerm{T}} # p*x*y → becomes quadratic + ppv::Vector{_ScalarCubicTerm{T}} # p*q*x → becomes affine + ppp::Vector{_ScalarCubicTerm{T}} # p*q*r → becomes constant + + # === Quadratic terms (degree 2) - same pattern as ParametricQuadraticFunction === + pv::Vector{MOI.ScalarQuadraticTerm{T}} # p*x → becomes affine + pp::Vector{MOI.ScalarQuadraticTerm{T}} # p*q → becomes constant + vv::Vector{MOI.ScalarQuadraticTerm{T}} # x*y → stays quadratic + + # === Affine terms (degree 1) === + p::Vector{MOI.ScalarAffineTerm{T}} # p → becomes constant + v::Vector{MOI.ScalarAffineTerm{T}} # x → stays affine + + # === Constant (degree 0) === + c::T + + # === Caches for efficient updates === + # Quadratic coefficients (from vv + pvv terms) - tracks current values in solver + quadratic_data::Dict{Tuple{MOI.VariableIndex,MOI.VariableIndex},T} + # Affine coefficients (from v + pv + ppv terms) + affine_data::Dict{MOI.VariableIndex,T} + # Affine coefficients not dependent on parameters + affine_data_np::Dict{MOI.VariableIndex,T} + # Current constant after parameter substitution + current_constant::T + # Set constant (for constraint handling, not used for objectives) + set_constant::T +end + +""" + ParametricCubicFunction(parsed::_ParsedCubicExpression{T}) where {T} + +Construct a ParametricCubicFunction from a _ParsedCubicExpression. +""" +function ParametricCubicFunction(parsed::_ParsedCubicExpression{T}) where {T} + + # Find variables related to parameters (from pv and ppv terms) + v_in_param_terms = Set{MOI.VariableIndex}() + for term in parsed.pv + push!(v_in_param_terms, term.variable_2) + end + for term in parsed.ppv + var = term.index_3 + push!(v_in_param_terms, var) + end + + # Split affine data + affine_data = Dict{MOI.VariableIndex,T}() + affine_data_np = Dict{MOI.VariableIndex,T}() + for term in parsed.v + if term.variable in v_in_param_terms + affine_data[term.variable] = + get(affine_data, term.variable, zero(T)) + term.coefficient + else + affine_data_np[term.variable] = + get(affine_data_np, term.variable, zero(T)) + term.coefficient + end + end + + # Find variable pairs related to parameters (from pvv terms) + var_pairs_in_param_terms = Set{Tuple{MOI.VariableIndex,MOI.VariableIndex}}() + for term in parsed.pvv + first_is_greater = term.index_2.value > term.index_3.value + v1 = ifelse(first_is_greater, term.index_3, term.index_2) + v2 = ifelse(first_is_greater, term.index_2, term.index_3) + push!(var_pairs_in_param_terms, (v1, v2)) + end + + # Initialize quadratic data + # Note: vv terms come from the parsed quadratic_func, which already has + # the MOI coefficient convention applied. We need to convert to internal form. + # MOI convention: + # - Off-diagonal (v1 != v2): coefficient C means C*v1*v2 (use as-is) + # - Diagonal (v1 == v2): coefficient C means (C/2)*v1^2 (divide by 2) + quadratic_data = Dict{Tuple{MOI.VariableIndex,MOI.VariableIndex},T}() + for term in parsed.vv + first_is_greater = term.variable_1.value > term.variable_2.value + v1 = ifelse(first_is_greater, term.variable_2, term.variable_1) + v2 = ifelse(first_is_greater, term.variable_1, term.variable_2) + coef = term.coefficient + if term.variable_1 == term.variable_2 + coef = coef / 2 # Diagonal: undo MOI's factor + end + quadratic_data[(v1, v2)] = get(quadratic_data, (v1, v2), zero(T)) + coef + end + # Add entries for pvv terms (will be updated with parameter values later) + for pair in var_pairs_in_param_terms + if !haskey(quadratic_data, pair) + quadratic_data[pair] = zero(T) + end + end + + return ParametricCubicFunction{T}( + parsed.pvv, + parsed.ppv, + parsed.ppp, + parsed.pv, + parsed.pp, + parsed.vv, + parsed.p, + parsed.v, + parsed.constant, + quadratic_data, + affine_data, + affine_data_np, + zero(T), # current_constant (computed later) + zero(T), # set_constant + ) +end + +# Accessors for cubic terms by type (direct field access) +_cubic_pvv_terms(f::ParametricCubicFunction) = f.pvv +_cubic_ppv_terms(f::ParametricCubicFunction) = f.ppv +_cubic_ppp_terms(f::ParametricCubicFunction) = f.ppp + +""" + _effective_param_value(model, pi::ParameterIndex) + +Get the effective parameter value: updated value if available, otherwise current value. +""" +function _effective_param_value(model, pi::ParameterIndex) + if haskey(model.updated_parameters, pi) && + !isnan(model.updated_parameters[pi]) + return model.updated_parameters[pi] + end + return model.parameters[pi] +end + +""" + _parametric_constant(model, f::ParametricCubicFunction{T}) where {T} + +Compute the constant term after parameter substitution. +Includes contributions from: c + p terms + pp terms + ppp cubic terms +""" +function _parametric_constant(model, f::ParametricCubicFunction{T}) where {T} + constant = f.c + + # From affine parameter terms (p) + for term in f.p + constant += + term.coefficient * + _effective_param_value(model, p_idx(term.variable)) + end + + # From quadratic parameter-parameter terms (pp) + # MOI convention: diagonal C means C/2*p^2, off-diagonal C means C*p1*p2 + for term in f.pp + divisor = term.variable_1 == term.variable_2 ? 2 : 1 + constant += + (term.coefficient / divisor) * + _effective_param_value(model, p_idx(term.variable_1)) * + _effective_param_value(model, p_idx(term.variable_2)) + end + + # From cubic ppp terms (all 3 indices are parameters) + for term in _cubic_ppp_terms(f) + p1 = term.index_1 + p2 = term.index_2 + p3 = term.index_3 + constant += + term.coefficient * + _effective_param_value(model, p_idx(p1)) * + _effective_param_value(model, p_idx(p2)) * + _effective_param_value(model, p_idx(p3)) + end + + return constant +end + +""" + _parametric_affine_terms(model, f::ParametricCubicFunction{T}) where {T} + +Compute affine coefficients after parameter substitution. +Includes contributions from: v terms + pv terms + ppv cubic terms +""" +function _parametric_affine_terms( + model, + f::ParametricCubicFunction{T}, +) where {T} + # Start with non-parametric terms + terms_dict = copy(f.affine_data) + + # Add contributions from pv terms (parameter * variable) + # These are always off-diagonal (p != v), so coefficient is used as-is + for term in f.pv + var = term.variable_2 + coef = term.coefficient + p_val = _effective_param_value(model, p_idx(term.variable_1)) + terms_dict[var] = get(terms_dict, var, zero(T)) + coef * p_val + end + + # Add contributions from ppv cubic terms + for term in _cubic_ppv_terms(f) + var = term.index_3 + p1_val = _effective_param_value(model, p_idx(term.index_1)) + p2_val = _effective_param_value(model, p_idx(term.index_2)) + terms_dict[var] = + get(terms_dict, var, zero(T)) + term.coefficient * p1_val * p2_val + end + + return terms_dict +end + +""" + _parametric_quadratic_terms(model, f::ParametricCubicFunction{T}) where {T} + +Compute quadratic coefficients after parameter substitution. +Includes contributions from: vv terms + pvv terms +""" +function _parametric_quadratic_terms( + model, + f::ParametricCubicFunction{T}, +) where {T} + # Start with vv terms + terms_dict = copy(f.quadratic_data) + + # Add contributions from pvv cubic terms + for term in _cubic_pvv_terms(f) + p = term.index_1 + first_is_greater = term.index_2.value > term.index_3.value + v1 = ifelse(first_is_greater, term.index_3, term.index_2) + v2 = ifelse(first_is_greater, term.index_2, term.index_3) + var_pair = (v1, v2) + p_val = _effective_param_value(model, p_idx(p)) + terms_dict[var_pair] = + get(terms_dict, var_pair, zero(T)) + term.coefficient * p_val + end + + return terms_dict +end + +""" + _current_function(f::ParametricCubicFunction{T}, model) where {T} + +Evaluate the cubic function with current parameter values and return +the appropriate MOI function type. +""" +function _current_function(f::ParametricCubicFunction{T}, model) where {T} + # Get current values + quad_data = _parametric_quadratic_terms(model, f) + affine_data = _parametric_affine_terms(model, f) + constant = _parametric_constant(model, f) + + # Build quadratic terms + # MOI convention: + # - Off-diagonal (v1 != v2): coefficient C means C*v1*v2 (use as-is) + # - Diagonal (v1 == v2): coefficient C means (C/2)*v1^2 (multiply by 2) + quadratic_terms = MOI.ScalarQuadraticTerm{T}[] + for ((v1, v2), coef) in quad_data + if !iszero(coef) + moi_coef = v1 == v2 ? coef * 2 : coef + push!(quadratic_terms, MOI.ScalarQuadraticTerm{T}(moi_coef, v1, v2)) + end + end + + # Build affine terms + affine_terms = MOI.ScalarAffineTerm{T}[] + for (v, coef) in affine_data + if !iszero(coef) + push!(affine_terms, MOI.ScalarAffineTerm{T}(coef, v)) + end + end + # Add non-parametric affine terms + for (v, coef) in f.affine_data_np + if !iszero(coef) + push!(affine_terms, MOI.ScalarAffineTerm{T}(coef, v)) + end + end + + # Note: We don't update f.affine_data or f.quadratic_data here. + # These store the BASE coefficients (from v and vv terms) and must remain unchanged. + # current_constant is the only cache we update for reference. + f.current_constant = constant + + # Always return a ScalarQuadraticFunction, even if it has no quadratic terms. + return MOI.ScalarQuadraticFunction{T}( + quadratic_terms, + affine_terms, + constant, + ) +end + +# === Delta functions for efficient updates === + +""" + _delta_parametric_constant(model, f::ParametricCubicFunction{T}) where {T} + +Compute the CHANGE in constant when parameters are updated. +""" +function _delta_parametric_constant( + model, + f::ParametricCubicFunction{T}, +) where {T} + delta = zero(T) + + # From p terms + for term in f.p + p_i = p_idx(term.variable) + if haskey(model.updated_parameters, p_i) && + !isnan(model.updated_parameters[p_i]) + old_val = model.parameters[p_i] + new_val = model.updated_parameters[p_i] + delta += term.coefficient * (new_val - old_val) + end + end + + # From pp terms + for term in f.pp + pi1 = p_idx(term.variable_1) + pi2 = p_idx(term.variable_2) + updated1 = + haskey(model.updated_parameters, pi1) && + !isnan(model.updated_parameters[pi1]) + updated2 = + haskey(model.updated_parameters, pi2) && + !isnan(model.updated_parameters[pi2]) + + if updated1 || updated2 + divisor = term.variable_1 == term.variable_2 ? 2 : 1 + old_val = + (term.coefficient / divisor) * + model.parameters[pi1] * + model.parameters[pi2] + new_p1 = + updated1 ? model.updated_parameters[pi1] : model.parameters[pi1] + new_p2 = + updated2 ? model.updated_parameters[pi2] : model.parameters[pi2] + new_val = (term.coefficient / divisor) * new_p1 * new_p2 + delta += new_val - old_val + end + end + + # From ppp cubic terms + for term in _cubic_ppp_terms(f) + pi1 = p_idx(term.index_1) + pi2 = p_idx(term.index_2) + pi3 = p_idx(term.index_3) + updated1 = + haskey(model.updated_parameters, pi1) && + !isnan(model.updated_parameters[pi1]) + updated2 = + haskey(model.updated_parameters, pi2) && + !isnan(model.updated_parameters[pi2]) + updated3 = + haskey(model.updated_parameters, pi3) && + !isnan(model.updated_parameters[pi3]) + + if updated1 || updated2 || updated3 + old_val = + term.coefficient * + model.parameters[pi1] * + model.parameters[pi2] * + model.parameters[pi3] + new_p1 = + updated1 ? model.updated_parameters[pi1] : model.parameters[pi1] + new_p2 = + updated2 ? model.updated_parameters[pi2] : model.parameters[pi2] + new_p3 = + updated3 ? model.updated_parameters[pi3] : model.parameters[pi3] + new_val = term.coefficient * new_p1 * new_p2 * new_p3 + delta += new_val - old_val + end + end + + return delta +end + +""" + _delta_parametric_affine_terms(model, f::ParametricCubicFunction{T}) where {T} + +Compute the CHANGE in affine coefficients when parameters are updated. +""" +function _delta_parametric_affine_terms( + model, + f::ParametricCubicFunction{T}, +) where {T} + delta_dict = Dict{MOI.VariableIndex,T}() + + # From pv terms (parameter * variable, always off-diagonal) + for term in f.pv + p_i = p_idx(term.variable_1) + if haskey(model.updated_parameters, p_i) && + !isnan(model.updated_parameters[p_i]) + var = term.variable_2 + coef = term.coefficient # Off-diagonal: use as-is + old_val = model.parameters[p_i] + new_val = model.updated_parameters[p_i] + delta_dict[var] = + get(delta_dict, var, zero(T)) + coef * (new_val - old_val) + end + end + + # From ppv cubic terms + for term in _cubic_ppv_terms(f) + var = term.index_3 + pi1 = p_idx(term.index_1) + pi2 = p_idx(term.index_2) + updated1 = + haskey(model.updated_parameters, pi1) && + !isnan(model.updated_parameters[pi1]) + updated2 = + haskey(model.updated_parameters, pi2) && + !isnan(model.updated_parameters[pi2]) + + if updated1 || updated2 + old_val = + term.coefficient * model.parameters[pi1] * model.parameters[pi2] + new_p1 = + updated1 ? model.updated_parameters[pi1] : model.parameters[pi1] + new_p2 = + updated2 ? model.updated_parameters[pi2] : model.parameters[pi2] + new_val = term.coefficient * new_p1 * new_p2 + delta_dict[var] = + get(delta_dict, var, zero(T)) + (new_val - old_val) + end + end + + return delta_dict +end + +""" + _delta_parametric_quadratic_terms(model, f::ParametricCubicFunction{T}) where {T} + +Compute the CHANGE in quadratic coefficients when parameters are updated. +""" +function _delta_parametric_quadratic_terms( + model, + f::ParametricCubicFunction{T}, +) where {T} + delta_dict = Dict{Tuple{MOI.VariableIndex,MOI.VariableIndex},T}() + + for term in _cubic_pvv_terms(f) + p_i = p_idx(term.index_1) + first_is_greater = term.index_2.value > term.index_3.value + v1 = ifelse(first_is_greater, term.index_3, term.index_2) + v2 = ifelse(first_is_greater, term.index_2, term.index_3) + var_pair = (v1, v2) + + if haskey(model.updated_parameters, p_i) && + !isnan(model.updated_parameters[p_i]) + old_val = model.parameters[p_i] + new_val = model.updated_parameters[p_i] + delta = term.coefficient * (new_val - old_val) + delta_dict[var_pair] = get(delta_dict, var_pair, zero(T)) + delta + end + end + + return delta_dict +end diff --git a/src/update_parameters.jl b/src/update_parameters.jl index 4f4b08b7..0a365cd6 100644 --- a/src/update_parameters.jl +++ b/src/update_parameters.jl @@ -296,6 +296,7 @@ function update_parameters!(model::Optimizer) _update_vector_quadratic_constraints!(model) _update_affine_objective!(model) _update_quadratic_objective!(model) + _update_cubic_objective!(model) # Update parameters and put NaN to indicate that the parameter has been # updated diff --git a/test/test_MathOptInterface.jl b/test/test_MathOptInterface.jl index 268e1b5f..de82313c 100644 --- a/test/test_MathOptInterface.jl +++ b/test/test_MathOptInterface.jl @@ -297,14 +297,27 @@ function test_moi_highs() ) MOI.set(model, MOI.Silent(), true) MOI.set(model, MOI.RawOptimizerAttribute("presolve"), "off") - MOI.Test.runtests(model, MOI.Test.Config(; atol = 1e-7); exclude = []) + # Exclude nonlinear tests that use functions outside POI's cubic polynomial support + # (general nonlinear functions like sin/cos/exp, or degree > 3 polynomials) + # POI only supports ScalarNonlinearFunction when it's a cubic polynomial with + # parameters multiplying at most quadratic variable terms + nonlinear_excludes = ["test_nonlinear_duals", "test_nonlinear_expression_"] + MOI.Test.runtests( + model, + MOI.Test.Config(; atol = 1e-7); + exclude = nonlinear_excludes, + ) model = POI.Optimizer( MOI.instantiate(HiGHS.Optimizer; with_bridge_type = Float64), ) MOI.set(model, MOI.Silent(), true) MOI.set(model, MOI.RawOptimizerAttribute("presolve"), "off") - MOI.Test.runtests(model, MOI.Test.Config(; atol = 1e-7); exclude = []) + MOI.Test.runtests( + model, + MOI.Test.Config(; atol = 1e-7); + exclude = nonlinear_excludes, + ) return end @@ -351,6 +364,9 @@ function test_moi_ipopt() # - CachingOptimizer does not throw if optimizer not attached "test_model_copy_to_UnsupportedAttribute", "test_model_copy_to_UnsupportedConstraint", + # - POI only supports cubic polynomial ScalarNonlinearFunction + "test_nonlinear_duals", + "test_nonlinear_expression_", ], ) return @@ -2106,7 +2122,7 @@ MOI.Utilities.@model( (MOI.ScalarAffineFunction,), (), () -); +) MOI.Utilities.@model( Model185_2, @@ -2118,7 +2134,7 @@ MOI.Utilities.@model( (), (), (MOI.VectorAffineFunction,) -); +) function test_issue_185() inner = Model185{Float64}() diff --git a/test/test_cubic.jl b/test/test_cubic.jl new file mode 100644 index 00000000..4da52914 --- /dev/null +++ b/test/test_cubic.jl @@ -0,0 +1,1628 @@ +# Copyright (c) 2020: Tomás Gutierrez and contributors +# +# Use of this source code is governed by an MIT-style license that can be found +# in the LICENSE.md file or at https://opensource.org/licenses/MIT. + +module TestCubic + +using Test +using JuMP + +import HiGHS +import ParametricOptInterface as POI +import MathOptInterface as MOI + +const ATOL = 1e-4 + +function runtests() + for name in names(@__MODULE__; all = true) + if startswith("$name", "test_") + @testset "$(name)" begin + getfield(@__MODULE__, name)() + end + end + end + return +end + +# ============================================================================ +# Helper mock optimizer +# ============================================================================ + +""" +Mock optimizer that rejects ScalarQuadraticFunction objectives. +""" +struct NoQuadObjModel <: MOI.ModelLike + inner::MOI.Utilities.Model{Float64} +end + +NoQuadObjModel() = NoQuadObjModel(MOI.Utilities.Model{Float64}()) + +MOI.add_variable(m::NoQuadObjModel) = MOI.add_variable(m.inner) +MOI.is_empty(m::NoQuadObjModel) = MOI.is_empty(m.inner) +MOI.empty!(m::NoQuadObjModel) = MOI.empty!(m.inner) +function MOI.is_valid(m::NoQuadObjModel, vi::MOI.VariableIndex) + return MOI.is_valid(m.inner, vi) +end + +function MOI.set( + ::NoQuadObjModel, + ::MOI.ObjectiveFunction{<:MOI.ScalarQuadraticFunction}, + ::MOI.ScalarQuadraticFunction, +) + return error("Quadratic objectives not supported") +end + +function MOI.set(m::NoQuadObjModel, attr::MOI.ObjectiveSense, v) + return MOI.set(m.inner, attr, v) +end + +# ============================================================================ +# Cubic Types +# ============================================================================ + +function test_normalize_cubic_indices_mixed() + x = MOI.VariableIndex(2) + y = MOI.VariableIndex(1) + p = POI.v_idx(POI.ParameterIndex(1)) + + n1, n2, n3 = POI._normalize_cubic_indices(x, p, y) + # Parameter should come first + @test POI._is_parameter(n1) + @test !POI._is_parameter(n2) + @test !POI._is_parameter(n3) + return +end + +function test_make_cubic_term_normalization() + x = MOI.VariableIndex(2) + y = MOI.VariableIndex(1) + p = POI.v_idx(POI.ParameterIndex(1)) + + term = POI._make_cubic_term(3.0, x, p, y) + @test term.coefficient == 3.0 + # index_1 should be parameter + @test POI._is_parameter(term.index_1) + return +end + +# ============================================================================ +# Parser - Valid cubic terms (pvv) +# ============================================================================ + +function test_cubic_parse_single_pvv_term() + x = MOI.VariableIndex(1) + y = MOI.VariableIndex(2) + p = POI.v_idx(POI.ParameterIndex(1)) + + # 2 * x * y * p + f = MOI.ScalarNonlinearFunction(:*, Any[2.0, x, y, p]) + result = POI._parse_cubic_expression(f, Float64) + + @test result !== nothing + pvv = result.pvv + @test length(pvv) == 1 + @test pvv[1].coefficient == 2.0 + @test pvv[1].index_1 == p + @test pvv[1].index_2 == x + @test pvv[1].index_3 == y + return +end + +function test_cubic_parse_squared_variable() + x = MOI.VariableIndex(1) + p = POI.v_idx(POI.ParameterIndex(1)) + + # 3 * p * x^2 using power operator + f = MOI.ScalarNonlinearFunction( + :*, + Any[3.0, p, MOI.ScalarNonlinearFunction(:^, Any[x, 2])], + ) + result = POI._parse_cubic_expression(f, Float64) + + @test result !== nothing + pvv = result.pvv + @test length(pvv) == 1 + @test pvv[1].coefficient == 3.0 + # Check that both variables are x (squared variable) + v1 = pvv[1].index_2 + v2 = pvv[1].index_3 + @test v1 == x + @test v2 == x + return +end + +function test_cubic_parse_parenthesis_variations() + x = MOI.VariableIndex(1) + y = MOI.VariableIndex(2) + p = POI.v_idx(POI.ParameterIndex(1)) + + # Flat: 2 * x * y * p + f1 = MOI.ScalarNonlinearFunction(:*, Any[2.0, x, y, p]) + + # Left-associative: ((2*x)*y)*p + f2 = MOI.ScalarNonlinearFunction( + :*, + Any[ + MOI.ScalarNonlinearFunction( + :*, + Any[MOI.ScalarNonlinearFunction(:*, Any[2.0, x]), y], + ), + p, + ], + ) + + # Grouped: (2*p) * (x*y) + f3 = MOI.ScalarNonlinearFunction( + :*, + Any[ + MOI.ScalarNonlinearFunction(:*, Any[2.0, p]), + MOI.ScalarNonlinearFunction(:*, Any[x, y]), + ], + ) + + r1 = POI._parse_cubic_expression(f1, Float64) + r2 = POI._parse_cubic_expression(f2, Float64) + r3 = POI._parse_cubic_expression(f3, Float64) + + # All should parse to equivalent results + for r in [r1, r2, r3] + @test r !== nothing + pvv = r.pvv + @test length(pvv) == 1 + @test pvv[1].coefficient == 2.0 + end + return +end + +function test_cubic_parse_multiple_pvv_different_vars() + x = MOI.VariableIndex(1) + y = MOI.VariableIndex(2) + z = MOI.VariableIndex(3) + p = POI.v_idx(POI.ParameterIndex(1)) + + # p*x*y + 2*p*x*z (two different pvv terms) + f = MOI.ScalarNonlinearFunction( + :+, + Any[ + MOI.ScalarNonlinearFunction(:*, Any[1.0, p, x, y]), + MOI.ScalarNonlinearFunction(:*, Any[2.0, p, x, z]), + ], + ) + result = POI._parse_cubic_expression(f, Float64) + @test result !== nothing + @test length(result.pvv) == 2 + return +end + +# ============================================================================ +# Parser - Valid cubic terms (ppv, ppp) +# ============================================================================ + +function test_cubic_parse_ppv_term() + x = MOI.VariableIndex(1) + p = POI.v_idx(POI.ParameterIndex(1)) + q = POI.v_idx(POI.ParameterIndex(2)) + + # 2 * p * q * x + f = MOI.ScalarNonlinearFunction(:*, Any[2.0, p, q, x]) + result = POI._parse_cubic_expression(f, Float64) + + @test result !== nothing + ppv = result.ppv + @test length(ppv) == 1 + @test ppv[1].coefficient == 2.0 + return +end + +function test_cubic_parse_ppp_term() + p = POI.v_idx(POI.ParameterIndex(1)) + q = POI.v_idx(POI.ParameterIndex(2)) + r = POI.v_idx(POI.ParameterIndex(3)) + + # 3 * p * q * r + f = MOI.ScalarNonlinearFunction(:*, Any[3.0, p, q, r]) + result = POI._parse_cubic_expression(f, Float64) + + @test result !== nothing + ppp = result.ppp + @test length(ppp) == 1 + @test ppp[1].coefficient == 3.0 + return +end + +# ============================================================================ +# Parser - Valid quadratic terms (pp, pv, vv) +# ============================================================================ + +function test_cubic_parse_pp_terms() + p = POI.v_idx(POI.ParameterIndex(1)) + q = POI.v_idx(POI.ParameterIndex(2)) + + # p * q (quadratic in parameters only) + f = MOI.ScalarNonlinearFunction(:*, Any[3.0, p, q]) + result = POI._parse_cubic_expression(f, Float64) + @test result !== nothing + @test length(result.pp) == 1 + @test result.pp[1].coefficient == 3.0 + return +end + +function test_cubic_parse_pp_same_parameter() + p = POI.v_idx(POI.ParameterIndex(1)) + + # p^2 (diagonal quadratic in parameters) + f = MOI.ScalarNonlinearFunction(:^, Any[p, 2]) + result = POI._parse_cubic_expression(f, Float64) + @test result !== nothing + @test length(result.pp) == 1 + return +end + +function test_cubic_parse_pv_terms() + x = MOI.VariableIndex(1) + p = POI.v_idx(POI.ParameterIndex(1)) + + # 4 * p * x (quadratic with one parameter and one variable) + f = MOI.ScalarNonlinearFunction(:*, Any[4.0, p, x]) + result = POI._parse_cubic_expression(f, Float64) + @test result !== nothing + @test length(result.pv) == 1 + @test result.pv[1].coefficient == 4.0 + return +end + +function test_cubic_parse_pv_param_first() + # Test the branch where m.variables[1] is already the parameter + p = POI.v_idx(POI.ParameterIndex(1)) + x = MOI.VariableIndex(1) + + # The monomial from _combine_like_monomials sorts by value, so the + # variable (lower value) comes first. We construct a monomial explicitly + # where the parameter is first in the raw expression to trigger both + # branches of the pv classification. + # With flat mult: p * x - monomials get combined with sorted key, + # so let's just verify both orderings work. + f1 = MOI.ScalarNonlinearFunction(:*, Any[p, x]) + f2 = MOI.ScalarNonlinearFunction(:*, Any[x, p]) + r1 = POI._parse_cubic_expression(f1, Float64) + r2 = POI._parse_cubic_expression(f2, Float64) + + for r in [r1, r2] + @test r !== nothing + @test length(r.pv) == 1 + # Convention: variable_1 = parameter, variable_2 = variable + @test POI._is_parameter(r.pv[1].variable_1) + @test !POI._is_parameter(r.pv[1].variable_2) + end + return +end + +# ============================================================================ +# Parser - Input type handling (SAF, SQF) +# ============================================================================ + +function test_parse_nonlinear_with_saf() + saf = MOI.ScalarAffineFunction( + [MOI.ScalarAffineTerm(1.0, MOI.VariableIndex(1))], + 1.3, + ) + f = MOI.ScalarNonlinearFunction(:+, Any[saf, 1.0]) + result = POI._parse_cubic_expression(f, Float64) + @test result.constant == 2.3 + return +end + +function test_cubic_parse_quadratic_function_input() + x = MOI.VariableIndex(1) + y = MOI.VariableIndex(2) + p = POI.v_idx(POI.ParameterIndex(1)) + + # ScalarQuadraticFunction as argument inside a nonlinear expression + sqf = MOI.ScalarQuadraticFunction( + [MOI.ScalarQuadraticTerm(2.0, x, y)], # off-diagonal: 2*x*y + [MOI.ScalarAffineTerm(3.0, x)], + 1.5, + ) + f = MOI.ScalarNonlinearFunction(:+, Any[sqf, p]) + result = POI._parse_cubic_expression(f, Float64) + @test result !== nothing + @test length(result.vv) == 1 + @test result.vv[1].coefficient == 2.0 + @test length(result.v) == 1 + @test result.v[1].coefficient == 3.0 + @test length(result.p) == 1 + @test result.constant == 1.5 + return +end + +function test_cubic_parse_quadratic_function_diagonal() + x = MOI.VariableIndex(1) + + # ScalarQuadraticFunction with diagonal term: x^2 + # MOI convention: diagonal coefficient C means (C/2)*x^2, so C=2 means x^2 + sqf = MOI.ScalarQuadraticFunction( + [MOI.ScalarQuadraticTerm(2.0, x, x)], # diagonal: (2/2)*x^2 = x^2 + MOI.ScalarAffineTerm{Float64}[], + 0.0, + ) + f = MOI.ScalarNonlinearFunction(:+, Any[sqf, 0.0]) + result = POI._parse_cubic_expression(f, Float64) + @test result !== nothing + @test length(result.vv) == 1 + # After parsing: coef should have MOI convention applied + return +end + +function test_cubic_parse_convenience_method() + x = MOI.VariableIndex(1) + p = POI.v_idx(POI.ParameterIndex(1)) + + # Test the convenience method without specifying type + f = MOI.ScalarNonlinearFunction(:*, Any[2.0, p, x, x]) + result = POI._parse_cubic_expression(f) + @test result !== nothing + @test length(result.pvv) == 1 + @test result.pvv[1].coefficient == 2.0 + return +end + +# ============================================================================ +# Parser - Operations (addition, subtraction, combination) +# ============================================================================ + +function test_cubic_parse_subtraction() + x = MOI.VariableIndex(1) + y = MOI.VariableIndex(2) + p = POI.v_idx(POI.ParameterIndex(1)) + + # x*y*p - 2*x (one cubic, one affine with negative coef) + f = MOI.ScalarNonlinearFunction( + :-, + Any[ + MOI.ScalarNonlinearFunction(:*, Any[x, y, p]), + MOI.ScalarNonlinearFunction(:*, Any[2.0, x]), + ], + ) + result = POI._parse_cubic_expression(f, Float64) + + @test result !== nothing + pvv = result.pvv + @test length(pvv) == 1 + # Check affine term via quadratic_func + affine = result.v + @test length(affine) == 1 + @test affine[1].coefficient == -2.0 + return +end + +function test_cubic_parse_subtraction_multiple_args() + x = MOI.VariableIndex(1) + y = MOI.VariableIndex(2) + p = POI.v_idx(POI.ParameterIndex(1)) + + # p*x*y - x - 1 (binary subtraction, second operand is sum) + f = MOI.ScalarNonlinearFunction( + :-, + Any[ + MOI.ScalarNonlinearFunction(:*, Any[p, x, y]), + MOI.ScalarNonlinearFunction(:+, Any[x, 1.0]), + ], + ) + result = POI._parse_cubic_expression(f, Float64) + @test result !== nothing + @test length(result.pvv) == 1 + @test length(result.v) == 1 + @test result.v[1].coefficient == -1.0 + @test result.constant == -1.0 + return +end + +function test_cubic_parse_unary_minus() + x = MOI.VariableIndex(1) + y = MOI.VariableIndex(2) + p = POI.v_idx(POI.ParameterIndex(1)) + + # -x*y*p (negation of cubic term) + f = MOI.ScalarNonlinearFunction( + :-, + Any[MOI.ScalarNonlinearFunction(:*, Any[x, y, p])], + ) + result = POI._parse_cubic_expression(f, Float64) + + @test result !== nothing + pvv = result.pvv + @test length(pvv) == 1 + @test pvv[1].coefficient == -1.0 + return +end + +function test_cubic_parse_term_combination() + x = MOI.VariableIndex(1) + y = MOI.VariableIndex(2) + p = POI.v_idx(POI.ParameterIndex(1)) + + # x*y*p + 2*x*y*p = 3*x*y*p (should combine into single term) + f = MOI.ScalarNonlinearFunction( + :+, + Any[ + MOI.ScalarNonlinearFunction(:*, Any[x, y, p]), + MOI.ScalarNonlinearFunction(:*, Any[2.0, x, y, p]), + ], + ) + result = POI._parse_cubic_expression(f, Float64) + + @test result !== nothing + pvv = result.pvv + @test length(pvv) == 1 # combined into single term + @test pvv[1].coefficient == 3.0 + return +end + +function test_cubic_parse_zero_coefficient_elimination() + x = MOI.VariableIndex(1) + y = MOI.VariableIndex(2) + p = POI.v_idx(POI.ParameterIndex(1)) + + # p*x*y - p*x*y = 0 (coefficients cancel) + f = MOI.ScalarNonlinearFunction( + :-, + Any[ + MOI.ScalarNonlinearFunction(:*, Any[1.0, p, x, y]), + MOI.ScalarNonlinearFunction(:*, Any[1.0, p, x, y]), + ], + ) + result = POI._parse_cubic_expression(f, Float64) + @test result !== nothing + @test isempty(result.pvv) + return +end + +# ============================================================================ +# Parser - Mixed degrees +# ============================================================================ + +function test_cubic_parse_mixed_all_degrees() + x = MOI.VariableIndex(1) + y = MOI.VariableIndex(2) + p = POI.v_idx(POI.ParameterIndex(1)) + q = POI.v_idx(POI.ParameterIndex(2)) + + # 2*p*x*y + 3*p*q*x + p*q*p + x^2 + 5*p*x + 7*x + 2*p + 10 + f = MOI.ScalarNonlinearFunction( + :+, + Any[ + MOI.ScalarNonlinearFunction(:*, Any[2.0, p, x, y]), # pvv + MOI.ScalarNonlinearFunction(:*, Any[3.0, p, q, x]), # ppv + MOI.ScalarNonlinearFunction(:*, Any[1.0, p, q, p]), # ppp + MOI.ScalarNonlinearFunction(:^, Any[x, 2]), # vv + MOI.ScalarNonlinearFunction(:*, Any[5.0, p, x]), # pv + MOI.ScalarNonlinearFunction(:*, Any[7.0, x]), # v + MOI.ScalarNonlinearFunction(:*, Any[2.0, p]), # p + 10.0, # constant + ], + ) + result = POI._parse_cubic_expression(f, Float64) + @test result !== nothing + @test length(result.pvv) == 1 + @test result.pvv[1].coefficient == 2.0 + @test length(result.ppv) == 1 + @test result.ppv[1].coefficient == 3.0 + @test length(result.ppp) == 1 + @test length(result.vv) == 1 + @test length(result.pv) == 1 + @test length(result.v) == 1 + @test result.v[1].coefficient == 7.0 + @test length(result.p) == 1 + @test result.p[1].coefficient == 2.0 + @test result.constant == 10.0 + return +end + +# ============================================================================ +# Parser - Invalid expressions (rejection) +# ============================================================================ + +function test_cubic_parse_non_polynomial_rejected() + x = MOI.VariableIndex(1) + p = POI.v_idx(POI.ParameterIndex(1)) + + # sin(x) * p - should be rejected + f = MOI.ScalarNonlinearFunction( + :*, + Any[MOI.ScalarNonlinearFunction(:sin, Any[x]), p], + ) + result = POI._parse_cubic_expression(f, Float64) + + @test result === nothing + + # sin(x) * p - should be rejected + f = MOI.ScalarNonlinearFunction( + :+, + Any[MOI.ScalarNonlinearFunction(:sin, Any[x]), p], + ) + result = POI._parse_cubic_expression(f, Float64) + + @test result === nothing + return +end + +function test_cubic_parse_invalid_degree_4() + x = MOI.VariableIndex(1) + y = MOI.VariableIndex(2) + z = MOI.VariableIndex(3) + p = POI.v_idx(POI.ParameterIndex(1)) + + # x * y * z * p (degree 4) should return nothing + f = MOI.ScalarNonlinearFunction(:*, Any[x, y, z, p]) + result = POI._parse_cubic_expression(f, Float64) + + @test result === nothing + return +end + +function test_cubic_parse_three_vars_no_param() + x = MOI.VariableIndex(1) + y = MOI.VariableIndex(2) + z = MOI.VariableIndex(3) + + # x * y * z (3 variables, 0 parameters) should be rejected + f = MOI.ScalarNonlinearFunction(:*, Any[x, y, z]) + result = POI._parse_cubic_expression(f, Float64) + + @test result === nothing + return +end + +function test_cubic_parse_unary_minus_invalid() + x = MOI.VariableIndex(1) + + # -(sin(x)) should propagate nothing from unary minus + f = MOI.ScalarNonlinearFunction( + :-, + Any[MOI.ScalarNonlinearFunction(:sin, Any[x])], + ) + result = POI._parse_cubic_expression(f, Float64) + @test result === nothing + return +end + +function test_cubic_parse_binary_subtraction_invalid_rhs() + x = MOI.VariableIndex(1) + + # x - sin(x) should propagate nothing from binary subtraction + f = MOI.ScalarNonlinearFunction( + :-, + Any[x, MOI.ScalarNonlinearFunction(:sin, Any[x])], + ) + result = POI._parse_cubic_expression(f, Float64) + @test result === nothing + return +end + +function test_cubic_parse_power_of_invalid_base() + x = MOI.VariableIndex(1) + + # sin(x)^2 should propagate nothing from power + f = MOI.ScalarNonlinearFunction( + :^, + Any[MOI.ScalarNonlinearFunction(:sin, Any[x]), 2], + ) + result = POI._parse_cubic_expression(f, Float64) + @test result === nothing + return +end + +# ============================================================================ +# Parser - Power operator edge cases +# ============================================================================ + +function test_cubic_parse_power_zero_exponent() + x = MOI.VariableIndex(1) + + # x^0 = 1 (constant) + f = MOI.ScalarNonlinearFunction( + :+, + Any[MOI.ScalarNonlinearFunction(:^, Any[x, 0]), 5.0], + ) + result = POI._parse_cubic_expression(f, Float64) + @test result !== nothing + @test result.constant == 6.0 + return +end + +function test_cubic_parse_power_negative_exponent() + x = MOI.VariableIndex(1) + + # x^(-1) should be rejected + f = MOI.ScalarNonlinearFunction(:^, Any[x, -1]) + result = POI._parse_cubic_expression(f, Float64) + @test result === nothing + return +end + +function test_cubic_parse_power_non_integer_exponent() + x = MOI.VariableIndex(1) + + # x^1.5 should be rejected + f = MOI.ScalarNonlinearFunction(:^, Any[x, 1.5]) + result = POI._parse_cubic_expression(f, Float64) + @test result === nothing + return +end + +function test_cubic_parse_power_wrong_arity() + x = MOI.VariableIndex(1) + + # Power with 1 or 3 args should be rejected + f = MOI.ScalarNonlinearFunction(:^, Any[x]) + result = POI._parse_cubic_expression(f, Float64) + @test result === nothing + + f = MOI.ScalarNonlinearFunction(:^, Any[x, 2, 3]) + result = POI._parse_cubic_expression(f, Float64) + @test result === nothing + return +end + +# ============================================================================ +# Parser - Division edge cases +# ============================================================================ + +function test_cubic_parse_division_by_zero() + x = MOI.VariableIndex(1) + + # x / 0 should be rejected + f = MOI.ScalarNonlinearFunction(:/, Any[x, 0.0]) + result = POI._parse_cubic_expression(f, Float64) + @test result === nothing + return +end + +function test_cubic_parse_division_by_variable() + x = MOI.VariableIndex(1) + y = MOI.VariableIndex(2) + + # x / y should be rejected (variable denominator) + f = MOI.ScalarNonlinearFunction(:/, Any[x, y]) + result = POI._parse_cubic_expression(f, Float64) + @test result === nothing + return +end + +function test_cubic_parse_division_wrong_arity() + x = MOI.VariableIndex(1) + + # Division with 1 or 3 args should be rejected + f = MOI.ScalarNonlinearFunction(:/, Any[x]) + result = POI._parse_cubic_expression(f, Float64) + @test result === nothing + + f = MOI.ScalarNonlinearFunction(:/, Any[x, 2.0, 3.0]) + result = POI._parse_cubic_expression(f, Float64) + @test result === nothing + return +end + +# ============================================================================ +# ParametricCubicFunction - Constructor +# ============================================================================ + +function test_parametric_cubic_function_constructor_with_ppv() + x = MOI.VariableIndex(1) + p = POI.v_idx(POI.ParameterIndex(1)) + q = POI.v_idx(POI.ParameterIndex(2)) + + # p*q*x + 2*x (ppv term + affine term on same variable) + f = MOI.ScalarNonlinearFunction( + :+, + Any[ + MOI.ScalarNonlinearFunction(:*, Any[1.0, p, q, x]), + MOI.ScalarNonlinearFunction(:*, Any[2.0, x]), + ], + ) + parsed = POI._parse_cubic_expression(f, Float64) + @test parsed !== nothing + pcf = POI.ParametricCubicFunction(parsed) + @test length(pcf.ppv) == 1 + # x should be in affine_data since it's shared with ppv + @test haskey(pcf.affine_data, x) + @test pcf.affine_data[x] == 2.0 + return +end + +function test_parametric_cubic_function_constructor_np_affine() + x = MOI.VariableIndex(1) + y = MOI.VariableIndex(2) + p = POI.v_idx(POI.ParameterIndex(1)) + + # p*x^2 + 3*y (pvv term on x, non-parametric affine on y) + f = MOI.ScalarNonlinearFunction( + :+, + Any[ + MOI.ScalarNonlinearFunction(:*, Any[1.0, p, x, x]), + MOI.ScalarNonlinearFunction(:*, Any[3.0, y]), + ], + ) + parsed = POI._parse_cubic_expression(f, Float64) + @test parsed !== nothing + pcf = POI.ParametricCubicFunction(parsed) + # y should be in affine_data_np since it's not related to any pv or ppv term + @test haskey(pcf.affine_data_np, y) + @test pcf.affine_data_np[y] == 3.0 + return +end + +# ============================================================================ +# MOI Objective Interface +# ============================================================================ + +function test_update_cubic_objective_no_cache() + model = POI.Optimizer(HiGHS.Optimizer()) + # No cubic objective set, cache is nothing — should return early + POI._update_cubic_objective!(model) + return +end + +function test_cubic_objective_supports() + model = POI.Optimizer(HiGHS.Optimizer()) + @test MOI.supports( + model, + MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}(), + ) + return +end + +function test_cubic_objective_set_invalid_expression() + model = POI.Optimizer(HiGHS.Optimizer()) + + x = MOI.add_variable(model) + y = MOI.add_variable(model) + z = MOI.add_variable(model) + + # x*y*z (3 variables, no parameter) should error + f = MOI.ScalarNonlinearFunction(:*, Any[x, y, z]) + @test_throws ErrorException MOI.set( + model, + MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}(), + f, + ) + return +end + +function test_cubic_objective_set_error_on_inner_optimizer() + mock = NoQuadObjModel() + model = POI.Optimizer(mock) + + x = MOI.add_variable(model) + p, _ = MOI.add_constrained_variable(model, MOI.Parameter(2.0)) + p_v = POI.v_idx(POI.p_idx(p)) + + # p * x^2 -- parsed successfully but inner optimizer rejects SQF + f = MOI.ScalarNonlinearFunction(:*, Any[1.0, p_v, x, x]) + err = try + MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}(), f) + nothing + catch e + e + end + @test err isa ErrorException + @test occursin("Failed to set cubic objective function", err.msg) + @test occursin("Quadratic objectives not supported", err.msg) + return +end + +function test_cubic_objective_get_no_cache() + model = POI.Optimizer(HiGHS.Optimizer()) + + @test_throws ErrorException MOI.get( + model, + MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}(), + ) + return +end + +function test_cubic_objective_get_save_original() + model = POI.Optimizer( + HiGHS.Optimizer(); + save_original_objective_and_constraints = true, + ) + MOI.set(model, MOI.Silent(), true) + + x = MOI.add_variable(model) + p, _ = MOI.add_constrained_variable(model, MOI.Parameter(2.0)) + + # Set a cubic objective: p*x^2 + p_v = POI.v_idx(POI.p_idx(p)) + f = MOI.ScalarNonlinearFunction(:*, Any[1.0, p_v, x, x]) + MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}(), f) + + # Should be able to retrieve + retrieved = + MOI.get(model, MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}()) + @test retrieved isa MOI.ScalarNonlinearFunction + return +end + +function test_cubic_objective_get_no_save() + model = POI.Optimizer( + HiGHS.Optimizer(); + save_original_objective_and_constraints = false, + ) + MOI.set(model, MOI.Silent(), true) + + x = MOI.add_variable(model) + p, _ = MOI.add_constrained_variable(model, MOI.Parameter(2.0)) + + # Set a cubic objective: p*x^2 + p_v = POI.v_idx(POI.p_idx(p)) + f = MOI.ScalarNonlinearFunction(:*, Any[1.0, p_v, x, x]) + MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}(), f) + + # Should error since save_original is false + @test_throws ErrorException MOI.get( + model, + MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}(), + ) + return +end + +# ============================================================================ +# JuMP Integration - Basic pvv +# ============================================================================ + +function test_jump_cubic_pvv_basic() + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, 0 <= x <= 10) + @variable(model, -1 <= y <= 10) + @variable(model, p in MOI.Parameter(1.0)) + + # Minimize: x^2 + p*x*y + y^2 - 3x + @constraint(model, x + y >= 0) + @objective(model, Min, x^2 + p * x * y + y^2 - 3 * x) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ -3.0 atol = ATOL + @test value(x) ≈ 2.0 atol = ATOL + @test value(y) ≈ -1.0 atol = ATOL + + # Change p to 0.5 + set_parameter_value(p, 0.5) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ -2.4 atol = ATOL + @test value(x) ≈ 1.6 atol = ATOL + @test value(y) ≈ -0.4 atol = ATOL + + # Change p to 0 (removes cross term) + set_parameter_value(p, 0.0) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ -9 / 4 atol = ATOL + @test value(x) ≈ 3 / 2 atol = ATOL + @test value(y) ≈ 0.0 atol = ATOL + return +end + +function test_jump_cubic_pvv_same() + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, 0 <= x <= 10) + @variable(model, p in MOI.Parameter(1.0)) + + # Minimize: p*x^2 - 3x + @constraint(model, x >= 0) + @objective(model, Min, p * x^2 - 3 * x) + + # p=1: optimal at x=3/2, obj=-9/4 + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ -9 / 4 atol = ATOL + @test value(x) ≈ 3 / 2 atol = ATOL + + # p=0.5: optimal at x=3, obj=-9/2 + set_parameter_value(p, 0.5) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ -9 / 2 atol = ATOL + @test value(x) ≈ 3 atol = ATOL + + # p=0: linear, optimal at x=10, obj=-30 + set_parameter_value(p, 0.0) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ -30.0 atol = ATOL + @test value(x) ≈ 10.0 atol = ATOL + return +end + +function test_jump_cubic_pvv_with_constant_and_affine() + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, 0 <= x <= 10) + @variable(model, 0 <= y <= 10) + @variable(model, p in MOI.Parameter(2.0)) + + # Minimize: p*x*y + x^2 + y^2 - 6x - 6y + 10 + # With p=2: (x+y)^2 - 6(x+y) + 10 = (s-3)^2 + 1, optimal at s=3, obj=1 + @constraint(model, x >= 0) + @constraint(model, y >= 0) + @objective(model, Min, p * x * y + x^2 + y^2 - 6 * x - 6 * y + 10) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ 1.0 atol = ATOL + @test value(x) + value(y) ≈ 3.0 atol = ATOL + + # p=0: x^2 + y^2 - 6x - 6y + 10, optimal at x=3, y=3, obj=-8 + set_parameter_value(p, 0.0) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ -8.0 atol = ATOL + @test value(x) ≈ 3.0 atol = ATOL + @test value(y) ≈ 3.0 atol = ATOL + return +end + +function test_jump_cubic_pvv_multiple_cross_terms() + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, 0 <= x <= 10) + @variable(model, 0 <= y <= 10) + @variable(model, 0 <= z <= 10) + @variable(model, p in MOI.Parameter(1.0)) + + # Minimize: x^2 + y^2 + z^2 + p*x*y + p*x*z - 6x - 4y - 2z + # With p=0: separable quadratics, x=3,y=2,z=1, obj=-14 + @constraint(model, x >= 0) + @constraint(model, y >= 0) + @constraint(model, z >= 0) + @objective( + model, + Min, + x^2 + y^2 + z^2 + p * x * y + p * x * z - 6 * x - 4 * y - 2 * z, + ) + + set_parameter_value(p, 0.0) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ -14.0 atol = ATOL + + # With p=1: coupled system + set_parameter_value(p, 1.0) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + x1 = value(x) + y1 = value(y) + z1 = value(z) + obj = objective_value(model) + expected = x1^2 + y1^2 + z1^2 + x1 * y1 + x1 * z1 - 6x1 - 4y1 - 2z1 + @test obj ≈ expected atol = ATOL + return +end + +function test_jump_cubic_pvv_reverse_var_order() + # Exercise variable ordering by having pvv with y*x where y.index > x.index + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, 0 <= x <= 10) + @variable(model, -1 <= y <= 10) + @variable(model, p in MOI.Parameter(2.0)) + + # Minimize: x^2 + y^2 + p*y*x - 3x (y comes before x in term) + @constraint(model, x + y >= 0) + @objective(model, Min, x^2 + y^2 + p * y * x - 3 * x) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + x1 = value(x) + y1 = value(y) + @test objective_value(model) ≈ x1^2 + y1^2 + 2 * y1 * x1 - 3 * x1 atol = + ATOL + + # Change p=0 + set_parameter_value(p, 0.0) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test value(x) ≈ 3 / 2 atol = ATOL + @test value(y) ≈ 0.0 atol = ATOL + return +end + +# ============================================================================ +# JuMP Integration - Basic ppv +# ============================================================================ + +function test_jump_cubic_ppv_basic() + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, x >= 0) + @variable(model, p in MOI.Parameter(2.0)) + @variable(model, q in MOI.Parameter(3.0)) + + # Minimize: x + p*q*x = x*(1 + p*q) + # With p=2, q=3: 7x, optimal at x=1, obj=7 + @constraint(model, x >= 1) + @objective(model, Min, x + p * q * x) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ 7.0 atol = ATOL + + # p=1, q=1: 2x, obj=2 + set_parameter_value(p, 1.0) + set_parameter_value(q, 1.0) + optimize!(model) + @test objective_value(model) ≈ 2.0 atol = ATOL +end + +function test_jump_cubic_ppv_same() + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, x >= 0) + @variable(model, p in MOI.Parameter(2.0)) + + # Minimize: x + p^2*x = x*(1 + p^2) + # With p=2: 5x, optimal at x=1, obj=5 + @constraint(model, x >= 1) + @objective(model, Min, x + p * p * x) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ 5.0 atol = ATOL + + # p=1: 2x, obj=2 + set_parameter_value(p, 1.0) + optimize!(model) + @test objective_value(model) ≈ 2.0 atol = ATOL +end + +# ============================================================================ +# JuMP Integration - Basic ppp +# ============================================================================ + +function test_jump_cubic_ppp_basic() + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, x >= 0) + @variable(model, p in MOI.Parameter(2.0)) + @variable(model, q in MOI.Parameter(3.0)) + @variable(model, r in MOI.Parameter(4.0)) + + # Minimize: x + p*q*r + # With p=2, q=3, r=4: x + 24, optimal at x=1, obj=25 + @constraint(model, x >= 1) + @objective(model, Min, x + p * q * r) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ 25.0 atol = ATOL + + # p=1, q=1, r=1: x + 1, obj=2 + set_parameter_value(p, 1.0) + set_parameter_value(q, 1.0) + set_parameter_value(r, 1.0) + optimize!(model) + @test objective_value(model) ≈ 2.0 atol = ATOL + return +end + +function test_jump_cubic_ppp_same() + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, x >= 0) + @variable(model, p in MOI.Parameter(2.0)) + + # Minimize: x + p^3 + # With p=2: x + 8, optimal at x=1, obj=9 + @constraint(model, x >= 1) + @objective(model, Min, x + p * p * p) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ 9.0 atol = ATOL + + # p=1: x + 1, obj=2 + set_parameter_value(p, 1.0) + optimize!(model) + @test objective_value(model) ≈ 2.0 atol = ATOL + return +end + +# ============================================================================ +# JuMP Integration - Specific term types in objective (pp, p affine) +# ============================================================================ + +function test_jump_cubic_with_pp_terms() + # pp terms alongside a cubic term (forces ScalarNonlinearFunction path) + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, 0 <= x <= 10) + @variable(model, p in MOI.Parameter(2.0)) + @variable(model, q in MOI.Parameter(3.0)) + + # p*x^2 + p*q + x (pvv + pp + v) + # With p=2, q=3: 2x^2 + 6 + x, at x=0: obj=6 + @constraint(model, x >= 0) + @objective(model, Min, p * x^2 + p * q + x) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ 6.0 atol = ATOL + @test value(x) ≈ 0.0 atol = ATOL + + # p=1, q=1: x^2 + 1 + x, at x=0: obj=1 + set_parameter_value(p, 1.0) + set_parameter_value(q, 1.0) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ 1.0 atol = ATOL + @test value(x) ≈ 0.0 atol = ATOL + return +end + +function test_jump_cubic_with_pp_same_param() + # p^2 term alongside a cubic term + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, 0 <= x <= 10) + @variable(model, p in MOI.Parameter(3.0)) + + # p*x^2 + p*p + x (pvv + pp + v) + # With p=3: 3x^2 + 9 + x, at x=0: obj=9 + @constraint(model, x >= 0) + @objective(model, Min, p * x^2 + p * p + x) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ 9.0 atol = ATOL + @test value(x) ≈ 0.0 atol = ATOL + + # p=0: x, at x=0: obj=0 + set_parameter_value(p, 0.0) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ 0.0 atol = ATOL + @test value(x) ≈ 0.0 atol = ATOL + return +end + +function test_jump_cubic_p_affine_in_cubic_expr() + # p affine terms alongside pvv (forces cubic path) + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, 0 <= x <= 10) + @variable(model, p in MOI.Parameter(5.0)) + + # p*x^2 + 3*p + x (pvv + p_affine + v) + # With p=5: 5x^2 + 15 + x, at x=0: obj=15 + @constraint(model, x >= 0) + @objective(model, Min, p * x^2 + 3 * p + x) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ 15.0 atol = ATOL + @test value(x) ≈ 0.0 atol = ATOL + + # p=0: x, at x=0: obj=0 + set_parameter_value(p, 0.0) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ 0.0 atol = ATOL + return +end + +# ============================================================================ +# JuMP Integration - Mixed term types +# ============================================================================ + +function test_jump_cubic_mixed_pvv_ppv_ppp() + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, 0 <= x <= 10) + @variable(model, 0 <= y <= 10) + @variable(model, p in MOI.Parameter(1.0)) + @variable(model, q in MOI.Parameter(2.0)) + + # Minimize: p*x*y + p*q*x + p*q*p + x^2 + y^2 + # With p=1, q=2: x*y + 2*x + 2 + x^2 + y^2 + @constraint(model, x >= 0) + @constraint(model, y >= 0) + @objective(model, Min, p * x * y + p * q * x + p * q * p + x^2 + y^2) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + obj1 = objective_value(model) + x1 = value(x) + y1 = value(y) + @test obj1 ≈ 1.0 * x1 * y1 + 2.0 * x1 + 2.0 + x1^2 + y1^2 atol = ATOL + + # p=2, q=3: 2*x*y + 6*x + 12 + x^2 + y^2 + set_parameter_value(p, 2.0) + set_parameter_value(q, 3.0) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + obj2 = objective_value(model) + x2 = value(x) + y2 = value(y) + @test obj2 ≈ 2.0 * x2 * y2 + 6.0 * x2 + 12.0 + x2^2 + y2^2 atol = ATOL + return +end + +function test_jump_cubic_all_term_types_combined() + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, 0 <= x <= 10) + @variable(model, 0 <= y <= 10) + @variable(model, p in MOI.Parameter(1.0)) + @variable(model, q in MOI.Parameter(1.0)) + + # f = p*x*y + p*q*x + p*q*p + x^2 + y^2 + p*x + 2*p + 3*x + 5 + # With p=1,q=1: x^2 + y^2 + x*y + 5x + 8 + @constraint(model, x >= 0) + @constraint(model, y >= 0) + @objective( + model, + Min, + p * x * y + + p * q * x + + p * q * p + + x^2 + + y^2 + + p * x + + 2 * p + + 3 * x + + 5, + ) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + x1 = value(x) + y1 = value(y) + obj1 = objective_value(model) + expected1 = x1^2 + y1^2 + x1 * y1 + 5 * x1 + 8 + @test obj1 ≈ expected1 atol = ATOL + + # p=2, q=0.5: x^2 + y^2 + 2*x*y + 6*x + 11 + set_parameter_value(p, 2.0) + set_parameter_value(q, 0.5) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + x2 = value(x) + y2 = value(y) + obj2 = objective_value(model) + expected2 = x2^2 + y2^2 + 2 * x2 * y2 + 6 * x2 + 11 + @test obj2 ≈ expected2 atol = ATOL + return +end + +# ============================================================================ +# JuMP Integration - Parameter updates +# ============================================================================ + +function test_jump_cubic_parameter_initially_zero() + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, 0 <= x <= 10) + @variable(model, -1 <= y <= 10) + @variable(model, p in MOI.Parameter(0.0)) + + @constraint(model, x + y >= 0) + @objective(model, Min, x^2 + p * x * y + y^2 - 3 * x) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ -9 / 4 atol = ATOL + @test value(x) ≈ 3 / 2 atol = ATOL + @test value(y) ≈ 0.0 atol = ATOL + + # Change p to 0.5 + set_parameter_value(p, 0.5) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ -2.4 atol = ATOL + @test value(x) ≈ 1.6 atol = ATOL + @test value(y) ≈ -0.4 atol = ATOL + + return +end + +function test_jump_cubic_multiple_parameter_updates() + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, 0 <= x <= 10) + @variable(model, p in MOI.Parameter(1.0)) + + @constraint(model, x >= 0) + @objective(model, Min, p * x^2 - 4 * x) + + # p=1: x^2 - 4x, optimal at x=2, obj=-4 + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test value(x) ≈ 2.0 atol = ATOL + @test objective_value(model) ≈ -4.0 atol = ATOL + + # p=2: 2x^2 - 4x, optimal at x=1, obj=-2 + set_parameter_value(p, 2.0) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test value(x) ≈ 1.0 atol = ATOL + @test objective_value(model) ≈ -2.0 atol = ATOL + + # p=4: 4x^2 - 4x, optimal at x=0.5, obj=-1 + set_parameter_value(p, 4.0) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test value(x) ≈ 0.5 atol = ATOL + @test objective_value(model) ≈ -1.0 atol = ATOL + + # p=0.5: 0.5x^2 - 4x, optimal at x=4, obj=-8 + set_parameter_value(p, 0.5) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test value(x) ≈ 4.0 atol = ATOL + @test objective_value(model) ≈ -8.0 atol = ATOL + return +end + +function test_jump_cubic_partial_parameter_update() + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, x >= 0) + @variable(model, p in MOI.Parameter(2.0)) + @variable(model, q in MOI.Parameter(3.0)) + + # Minimize: x + p*q*x = x*(1 + p*q) + # With p=2, q=3: 7x, obj=7 + @constraint(model, x >= 1) + @objective(model, Min, x + p * q * x) + + optimize!(model) + @test objective_value(model) ≈ 7.0 atol = ATOL + + # Only update p to 0, keep q=3: x, obj=1 + set_parameter_value(p, 0.0) + optimize!(model) + @test objective_value(model) ≈ 1.0 atol = ATOL + return +end + +function test_jump_cubic_pvv_update_no_change() + # Test the case where parameter update results in no actual change + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, 0 <= x <= 10) + @variable(model, p in MOI.Parameter(1.0)) + + @constraint(model, x >= 0) + @objective(model, Min, p * x^2 - 2 * x) + + # First solve + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + obj1 = objective_value(model) + + # "Update" parameter to same value + set_parameter_value(p, 1.0) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ obj1 atol = ATOL + return +end + +# ============================================================================ +# JuMP Integration - Negative parameters +# ============================================================================ + +function test_jump_cubic_negative_parameters() + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, 0 <= x <= 10) + @variable(model, p in MOI.Parameter(-1.0)) + + # x^2 + p*x^2 - x = (1+p)*x^2 - x + # With p=-0.5: 0.5x^2 - x, optimal at x=1, obj=-0.5 + @constraint(model, x >= 0) + @objective(model, Min, x^2 + p * x^2 - x) + + set_parameter_value(p, -0.5) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ -0.5 atol = ATOL + @test value(x) ≈ 1.0 atol = ATOL + + # p=0: x^2 - x, optimal at x=0.5, obj=-0.25 + set_parameter_value(p, 0.0) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ -0.25 atol = ATOL + @test value(x) ≈ 0.5 atol = ATOL + return +end + +function test_jump_cubic_ppv_negative_params() + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, x >= 0) + @variable(model, p in MOI.Parameter(-1.0)) + @variable(model, q in MOI.Parameter(2.0)) + + # Minimize: x + p*q*x = x*(1 + p*q) = x*(1-2) = -x + # Subject to: 0 <= x <= 5 + @constraint(model, x <= 5) + @objective(model, Min, x + p * q * x) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + # f(x) = x + (-1)(2)(x) = -x, min at x=5, obj=-5 + @test objective_value(model) ≈ -5.0 atol = ATOL + @test value(x) ≈ 5.0 atol = ATOL + + # p=1, q=1: 2x, min at x=0 + set_parameter_value(p, 1.0) + set_parameter_value(q, 1.0) + optimize!(model) + @test objective_value(model) ≈ 0.0 atol = ATOL + @test value(x) ≈ 0.0 atol = ATOL + return +end + +function test_jump_cubic_ppp_negative_params() + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, x >= 0) + @variable(model, p in MOI.Parameter(-2.0)) + @variable(model, q in MOI.Parameter(3.0)) + @variable(model, r in MOI.Parameter(1.0)) + + # Minimize: x + p*q*r = x + (-2)(3)(1) = x - 6 + @constraint(model, x >= 1) + @objective(model, Min, x + p * q * r) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ -5.0 atol = ATOL + @test value(x) ≈ 1.0 atol = ATOL + + # p=1, q=1, r=1: x + 1, obj=2 + set_parameter_value(p, 1.0) + set_parameter_value(q, 1.0) + set_parameter_value(r, 1.0) + optimize!(model) + @test objective_value(model) ≈ 2.0 atol = ATOL + return +end + +# ============================================================================ +# JuMP Integration - Division and direct model +# ============================================================================ + +function test_jump_cubic_parameter_division_by_constant() + model = direct_model(POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, 0 <= x <= 10) + @variable(model, 0 <= y <= 10) + @variable(model, p in MOI.Parameter(0.0)) + + @constraint(model, x + y >= 0) + @objective(model, Min, x^2 + p * x * y / 1 + y^2 - 3 * x) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ -9 / 4 atol = ATOL + @test value(x) ≈ 3 / 2 atol = ATOL + @test value(y) ≈ 0.0 atol = ATOL + + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, 0 <= x <= 10) + @variable(model, -1 <= y <= 10) + @variable(model, p in MOI.Parameter(1.0)) + + # Minimize: x^2 + 0.5*x*y + y^2 - 3x + @constraint(model, x + y >= 0) + @objective(model, Min, x^2 + p * x * y / 2 + y^2 - 3 * x) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ -2.4 atol = ATOL + @test value(x) ≈ 1.6 atol = ATOL + @test value(y) ≈ -0.4 atol = ATOL + + return +end + +function test_jump_cubic_division_in_ppv() + model = Model(() -> POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, x >= 0) + @variable(model, p in MOI.Parameter(4.0)) + @variable(model, q in MOI.Parameter(6.0)) + + # Minimize: x + p*q*x/2 = x*(1 + p*q/2) + # With p=4, q=6: x*(1+12)=13x + @constraint(model, x >= 1) + @objective(model, Min, x + p * q * x / 2) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ 13.0 atol = ATOL + return +end + +function test_jump_cubic_direct_model_ppv() + model = direct_model(POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, x >= 0) + @variable(model, p in MOI.Parameter(3.0)) + @variable(model, q in MOI.Parameter(2.0)) + + # Minimize: x + p*q*x = x*(1 + 6) = 7x + @constraint(model, x >= 1) + @objective(model, Min, x + p * q * x) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ 7.0 atol = ATOL + + # p=0: x, obj=1 + set_parameter_value(p, 0.0) + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ 1.0 atol = ATOL + return +end + +function test_jump_cubic_direct_model_ppp() + model = direct_model(POI.Optimizer(HiGHS.Optimizer())) + set_silent(model) + + @variable(model, x >= 0) + @variable(model, p in MOI.Parameter(2.0)) + @variable(model, q in MOI.Parameter(3.0)) + @variable(model, r in MOI.Parameter(4.0)) + + # Minimize: x + p*q*r = x + 24 + @constraint(model, x >= 1) + @objective(model, Min, x + p * q * r) + + optimize!(model) + @test termination_status(model) in (OPTIMAL, LOCALLY_SOLVED) + @test objective_value(model) ≈ 25.0 atol = ATOL + + set_parameter_value(p, 0.0) + optimize!(model) + @test objective_value(model) ≈ 1.0 atol = ATOL + return +end + +end # module + +TestCubic.runtests()