From 018deed49cbc3a7537da9cd9f7b658be3d6fd775 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Fri, 27 Dec 2024 20:25:04 -0300 Subject: [PATCH 01/26] [WIP] Integrate with POI to improve UX --- Project.toml | 1 + src/DiffOpt.jl | 12 + src/diff_opt.jl | 11 + src/jump_moi_overloads.jl | 15 + src/moi_wrapper.jl | 77 ++++- src/parameters.jl | 629 ++++++++++++++++++++++++++++++++++++++ test/parameters.jl | 468 ++++++++++++++++++++++++++++ test/runtests.jl | 3 + 8 files changed, 1206 insertions(+), 10 deletions(-) create mode 100644 src/parameters.jl create mode 100644 test/parameters.jl diff --git a/Project.toml b/Project.toml index d73a81b42..2cc3a7687 100644 --- a/Project.toml +++ b/Project.toml @@ -12,6 +12,7 @@ LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" MathOptSetDistances = "3b969827-a86c-476c-9527-bb6f1a8fbad5" +ParametricOptInterface = "0ce4ce61-57bf-432b-a095-efac525d185e" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] diff --git a/src/DiffOpt.jl b/src/DiffOpt.jl index 9ab20101d..57c59fb5e 100644 --- a/src/DiffOpt.jl +++ b/src/DiffOpt.jl @@ -12,6 +12,7 @@ import LazyArrays import LinearAlgebra import MathOptInterface as MOI import MathOptSetDistances as MOSD +import ParametricOptInterface as POI import SparseArrays include("utils.jl") @@ -19,6 +20,7 @@ include("product_of_sets.jl") include("diff_opt.jl") include("moi_wrapper.jl") include("jump_moi_overloads.jl") +include("parameters.jl") include("copy_dual.jl") include("bridges.jl") @@ -38,6 +40,16 @@ function add_all_model_constructors(model) return end +function add_conic_model_constructor(model) + add_model_constructor(model, ConicProgram.Model) + return +end + +function add_quadratic_model_constructor(model) + add_model_constructor(model, QuadraticProgram.Model) + return +end + export diff_optimizer end # module diff --git a/src/diff_opt.jl b/src/diff_opt.jl index 77ba7d968..b76116d17 100644 --- a/src/diff_opt.jl +++ b/src/diff_opt.jl @@ -58,6 +58,17 @@ The output solution differentials can be queried with the attribute """ function forward_differentiate! end +""" + empty_input_sensitivities!(model::MOI.ModelLike) + +Empty the input sensitivities of the model. +Sets to zero all the sensitivities set by the user with method such as: +- `MOI.set(model, DiffOpt.ReverseVariablePrimal(), variable_index, value)` +- `MOI.set(model, DiffOpt.ForwardObjectiveFunction(), expression)` +- `MOI.set(model, DiffOpt.ForwardConstraintFunction(), index, expression)` +""" +function empty_input_sensitivities! end + """ ForwardObjectiveFunction <: MOI.AbstractModelAttribute diff --git a/src/jump_moi_overloads.jl b/src/jump_moi_overloads.jl index 9459f840f..6f8c29ae8 100644 --- a/src/jump_moi_overloads.jl +++ b/src/jump_moi_overloads.jl @@ -307,6 +307,11 @@ function forward_differentiate!(model::JuMP.Model) return forward_differentiate!(JuMP.backend(model)) end +function empty_input_sensitivities!(model::JuMP.Model) + empty_input_sensitivities!(JuMP.backend(model)) + return +end + # MOI.Utilities function reverse_differentiate!(model::MOI.Utilities.CachingOptimizer) @@ -317,6 +322,11 @@ function forward_differentiate!(model::MOI.Utilities.CachingOptimizer) return forward_differentiate!(model.optimizer) end +function empty_input_sensitivities!(model::MOI.Utilities.CachingOptimizer) + empty_input_sensitivities!(model.optimizer) + return +end + # MOIB function reverse_differentiate!(model::MOI.Bridges.AbstractBridgeOptimizer) @@ -326,3 +336,8 @@ end function forward_differentiate!(model::MOI.Bridges.AbstractBridgeOptimizer) return forward_differentiate!(model.model) end + +function empty_input_sensitivities!(model::MOI.Bridges.AbstractBridgeOptimizer) + empty_input_sensitivities!(model.model) + return +end diff --git a/src/moi_wrapper.jl b/src/moi_wrapper.jl index 3570d94d5..d577765a7 100644 --- a/src/moi_wrapper.jl +++ b/src/moi_wrapper.jl @@ -4,7 +4,37 @@ # in the LICENSE.md file or at https://opensource.org/licenses/MIT. """ - diff_optimizer(optimizer_constructor)::Optimizer + DiffMethod + +An enum to define the differentiation method. + +## Values + +Possible values are: + + * [`DIFF_AUTOMATIC`]: Automatic differentiation: tries all differentiation methods and chooses the first that works. This can be slower than manually choosing a method. + * [`DIFF_CONIC`]: Conic optimization based differentiation: works for conic programs or programs that can be transformed into conic programs by JuMP. + * [`DIFF_QUADRATIC`]: Quadratic optimization based differentiation: works for quadratic programs or programs that can be transformed into conic programs by JuMP. +""" +@enum(DiffMethod, DIFF_AUTOMATIC, DIFF_CONIC, DIFF_QUADRATIC) + +@doc( + "Automatic differentiation: tries all differentiation methods and chooses the first that works. This can be slower than manually choosing a method.", + DIFF_AUTOMATIC +) + +@doc( + "Conic optimization based differentiation: works for conic programs or programs that can be transformed into conic programs by JuMP.", + DIFF_CONIC +) + +@doc( + "Quadratic optimization based differentiation: works for quadratic programs or programs that can be transformed into conic programs by JuMP.", + DIFF_QUADRATIC +) + +""" + diff_optimizer(optimizer_constructor) Creates a `DiffOpt.Optimizer`, which is an MOI layer with an internal optimizer and other utility methods. Results (primal, dual and slack values) are obtained @@ -21,19 +51,33 @@ julia> x = model.add_variable(model) julia> model.add_constraint(model, ...) ``` """ -function diff_optimizer(optimizer_constructor)::Optimizer +function diff_optimizer( + optimizer_constructor; + method = DIFF_AUTOMATIC, + with_parametric_opt_interface::Bool = false, + with_bridge_type = Float64, + with_cache::Bool = true, +) optimizer = - MOI.instantiate(optimizer_constructor; with_bridge_type = Float64) + MOI.instantiate(optimizer_constructor; with_bridge_type = with_bridge_type) # When we do `MOI.copy_to(diff, optimizer)` we need to efficiently `MOI.get` # the model information from `optimizer`. However, 1) `optimizer` may not # implement some getters or it may be inefficient and 2) the getters may be # unimplemented or inefficient through some bridges. # For this reason we add a cache layer, the same cache JuMP adds. - caching_opt = MOI.Utilities.CachingOptimizer( - MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), - optimizer, - ) - return Optimizer(caching_opt) + caching_opt = if with_cache + MOI.Utilities.CachingOptimizer( + MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), + optimizer, + ) + else + optimizer + end + if with_parametric_opt_interface + return POI.Optimizer(Optimizer(caching_opt, method = method)) + else + return Optimizer(caching_opt, method = method) + end end mutable struct Optimizer{OT<:MOI.ModelLike} <: MOI.AbstractOptimizer @@ -49,10 +93,18 @@ mutable struct Optimizer{OT<:MOI.ModelLike} <: MOI.AbstractOptimizer # sensitivity input cache using MOI like sparse format input_cache::InputCache - function Optimizer(optimizer::OT) where {OT<:MOI.ModelLike} + function Optimizer(optimizer::OT; method = DIFF_AUTOMATIC) where {OT<:MOI.ModelLike} output = new{OT}(optimizer, Any[], nothing, nothing, nothing, InputCache()) - add_all_model_constructors(output) + if method == DIFF_CONIC + add_conic_model_constructor(output) + elseif method == DIFF_QUADRATIC + add_quadratic_model_constructor(output) + elseif method == DIFF_AUTOMATIC + add_all_model_constructors(output) + else + add_model_constructor(output, method) + end return output end end @@ -552,6 +604,11 @@ function forward_differentiate!(model::Optimizer) return forward_differentiate!(diff) end +function empty_input_sensitivities!(model::Optimizer) + empty!(model.input_cache) + return +end + function _instantiate_with_bridges(model_constructor) model = MOI.Bridges.LazyBridgeOptimizer(MOI.instantiate(model_constructor)) # We don't add any variable bridge here because: diff --git a/src/parameters.jl b/src/parameters.jl new file mode 100644 index 000000000..d35276fbf --- /dev/null +++ b/src/parameters.jl @@ -0,0 +1,629 @@ +# Copyright (c) 2020: Akshay Sharma 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. + +# block other methods + +MOI.supports(::POI.Optimizer, ::ForwardObjectiveFunction) = false + +function MOI.set( + ::POI.Optimizer, + ::ForwardObjectiveFunction, + _, +) + error( + "Forward objective function is not supported when " * + "`with_parametric_opt_interface` is set to `true` in " * + "`diff_optimizer`." * + "Use parameters to set the forward sensitivity.", + ) +end + +MOI.supports(::POI.Optimizer, ::ForwardConstraintFunction) = false + +function MOI.set( + ::POI.Optimizer, + ::ForwardConstraintFunction, + ::MOI.ConstraintIndex, + _, +) + error( + "Forward constraint function is not supported when " * + "`with_parametric_opt_interface` is set to `true` in " * + "`diff_optimizer`." * + "Use parameters to set the forward sensitivity.", + ) +end + +MOI.supports(::POI.Optimizer, ::ReverseObjectiveFunction) = false + +function MOI.get( + ::POI.Optimizer, + ::ReverseObjectiveFunction, +) + error( + "Reverse objective function is not supported when " * + "`with_parametric_opt_interface` is set to `true` in " * + "`diff_optimizer`." * + "Use parameters to get the reverse sensitivity.", + ) +end + +MOI.supports(::POI.Optimizer, ::ReverseConstraintFunction) = false + +function MOI.get( + ::POI.Optimizer, + ::ReverseConstraintFunction, + ::MOI.ConstraintIndex, +) + error( + "Reverse constraint function is not supported when " * + "`with_parametric_opt_interface` is set to `true` in " * + "`diff_optimizer`." * + "Use parameters to get the reverse sensitivity.", + ) +end + +# functions to be used with ParametricOptInterface.jl + +struct ForwardConstraintSet <: MOI.AbstractConstraintAttribute end + +struct ReverseConstraintSet <: MOI.AbstractConstraintAttribute end + +mutable struct SensitivityData{T} + parameter_input_forward::Dict{MOI.VariableIndex,T} + parameter_output_backward::Dict{MOI.VariableIndex,T} +end + +function SensitivityData{T}() where {T} + return SensitivityData{T}(Dict{MOI.VariableIndex,T}(), Dict{MOI.VariableIndex,T}()) +end + +function _get_sensitivity_data(model::POI.Optimizer{T})::SensitivityData{T} where {T} + _initialize_sensitivity_data!(model) + return model.ext[:_sensitivity_data]::SensitivityData{T} +end + +function _initialize_sensitivity_data!(model::POI.Optimizer{T}) where {T} + if !haskey(model.ext, :_sensitivity_data) + model.ext[:_sensitivity_data] = SensitivityData{T}() + end + return +end + +const DoubleDictInner = MOI.Utilities.DoubleDicts.DoubleDictInner + +# forward mode + +function _constraint_set_forward!( + model::POI.Optimizer{T}, + affine_constraint_cache_dict, + ::Type{P}, +) where {T, P <: POI.ParametricAffineFunction} + sensitivity_data = _get_sensitivity_data(model) + for (inner_ci, pf) in affine_constraint_cache_dict + cte = zero(T) + terms = MOI.ScalarAffineTerm{T}[] + sizehint!(terms, 0) + for term in POI.affine_parameter_terms(pf) + p = term.variable + sensitivity = get(sensitivity_data.parameter_input_forward, p, 0.0) + cte += sensitivity * term.coefficient + end + if !iszero(cte) + MOI.set( + model.optimizer, + ForwardConstraintFunction(), + inner_ci, + MOI.ScalarAffineFunction{T}(terms, cte), + ) + end + end + return +end + +function _constraint_set_forward!( + model::POI.Optimizer{T}, + vector_affine_constraint_cache_dict, + ::Type{P}, +) where {T, P <: POI.ParametricVectorAffineFunction} + sensitivity_data = _get_sensitivity_data(model) + for (inner_ci, pf) in vector_affine_constraint_cache_dict + cte = zeros(T, length(pf.c)) + terms = MOI.VectorAffineTerm{T}[] + sizehint!(terms, 0) + for term in POI.vector_affine_parameter_terms(pf) + p = term.scalar_term.variable + sensitivity = get(sensitivity_data.parameter_input_forward, p, 0.0) + cte[term.output_index] += sensitivity * term.scalar_term.coefficient + end + if !iszero(cte) + MOI.set( + model.optimizer, + ForwardConstraintFunction(), + inner_ci, + MOI.VectorAffineFunction{T}(terms, cte), + ) + end + end + return +end + +function _constraint_set_forward!( + model::POI.Optimizer{T}, + quadratic_constraint_cache_dict, + ::Type{P}, +) where {T, P <: POI.ParametricQuadraticFunction} + sensitivity_data = _get_sensitivity_data(model) + for (inner_ci, pf) in quadratic_constraint_cache_dict + cte = zero(T) + terms = MOI.ScalarAffineTerm{T}[] + for term in POI.affine_parameter_terms(pf) + p = term.variable + sensitivity = get(sensitivity_data.parameter_input_forward, p, 0.0) + cte += sensitivity * term.coefficient + end + for term in POI.quadratic_parameter_parameter_terms(pf) + p_1 = term.variable_1 + p_2 = term.variable_2 + sensitivity_1 = get(sensitivity_data.parameter_input_forward, p_1, 0.0) + sensitivity_2 = get(sensitivity_data.parameter_input_forward, p_2, 0.0) + cte += + sensitivity_1 * term.coefficient * + MOI.get(model, MOI.VariablePrimal(), p_2) / + ifelse(term.variable_1 === term.variable_2, 2, 1) + cte += + sensitivity_2 * term.coefficient * + MOI.get(model, MOI.VariablePrimal(), p_1) / + ifelse(term.variable_1 === term.variable_2, 2, 1) + end + sizehint!(terms, length(POI.quadratic_parameter_variable_terms(pf))) + for term in POI.quadratic_parameter_variable_terms(pf) + p = term.variable_1 + sensitivity = get(sensitivity_data.parameter_input_forward, p, NaN) + if !isnan(sensitivity) + push!( + terms, + MOI.ScalarAffineTerm{T}( + sensitivity * term.coefficient, + term.variable_2, + ), + ) + end + end + if !iszero(cte) || !isempty(terms) + MOI.set( + model.optimizer, + ForwardConstraintFunction(), + inner_ci, + MOI.ScalarAffineFunction{T}(terms, cte), + ) + end + end + return +end + +function _affine_objective_set_forward!(model::POI.Optimizer{T}) where {T} + cte = zero(T) + terms = MOI.ScalarAffineTerm{T}[] + pf = model.affine_objective_cache + sizehint!(terms, 0) + sensitivity_data = _get_sensitivity_data(model) + for term in POI.affine_parameter_terms(pf) + p = term.variable + sensitivity = get(sensitivity_data.parameter_input_forward, p, 0.0) + cte += sensitivity * term.coefficient + end + if !iszero(cte) + MOI.set( + model.optimizer, + ForwardObjectiveFunction(), + MOI.ScalarAffineFunction{T}(terms, cte), + ) + end + return +end + +function _quadratic_objective_set_forward!(model::POI.Optimizer{T}) where {T} + cte = zero(T) + pf = MOI.get(model, POI.ParametricObjectiveFunction{POI.ParametricQuadraticFunction{T}}()) + sensitivity_data = _get_sensitivity_data(model) + for term in POI.affine_parameter_terms(pf) + p = term.variable + sensitivity = get(sensitivity_data.parameter_input_forward, p, 0.0) + cte += sensitivity * term.coefficient + end + for term in POI.quadratic_parameter_parameter_terms(pf) + p_1 = term.variable_1 + p_2 = term.variable_2 + sensitivity_1 = get(sensitivity_data.parameter_input_forward, p_1, 0.0) + sensitivity_2 = get(sensitivity_data.parameter_input_forward, p_2, 0.0) + cte += + sensitivity_1 * term.coefficient * + MOI.get(model, MOI.VariablePrimal(), p_2) / + ifelse(term.variable_1 === term.variable_2, 2, 1) + cte += + sensitivity_2 * term.coefficient + MOI.get(model, MOI.VariablePrimal(), p_1) / + ifelse(term.variable_1 === term.variable_2, 2, 1) + end + terms = MOI.ScalarAffineTerm{T}[] + sizehint!(terms, length(POI.quadratic_parameter_variable_terms(pf))) + for term in POI.quadratic_parameter_variable_terms(pf) + p = term.variable_1 + sensitivity = get(sensitivity_data.parameter_input_forward, p, NaN) + if !isnan(sensitivity) + push!( + terms, + MOI.ScalarAffineTerm{T}( + sensitivity * term.coefficient, + term.variable_2, + ), + ) + end + end + if !iszero(cte) || !isempty(terms) + MOI.set( + model.optimizer, + ForwardObjectiveFunction(), + MOI.ScalarAffineFunction{T}(terms, cte), + ) + end + return +end + +function empty_input_sensitivities!(model::POI.Optimizer) + empty_input_sensitivities!(model.optimizer) + return +end + +function forward_differentiate!(model::POI.Optimizer{T}) where {T} + empty_input_sensitivities!(model) + ctr_types = MOI.get(model, POI.ListOfParametricConstraintTypesPresent()) + for (F, S, P) in ctr_types + dict = MOI.get(model, POI.DictOfParametricConstraintIndicesAndFunctions{F, S, P}()) + _constraint_set_forward!(model, dict, P) + end + obj_type = MOI.get(model, POI.ParametricObjectiveType()) + if obj_type <: POI.ParametricAffineFunction + _affine_objective_set_forward!(model) + elseif obj_type <: POI.ParametricQuadraticFunction + _quadratic_objective_set_forward!(model) + end + forward_differentiate!(model.optimizer) + return +end + +# function MOI.set( +# model::POI.Optimizer, +# ::ForwardParameter, +# variable::MOI.VariableIndex, +# value::Number, +# ) +# if _is_variable(model, variable) +# error("Trying to set a forward parameter sensitivity for a variable") +# end +# parameter = variable +# sensitivity_data = _get_sensitivity_data(model) +# sensitivity_data.parameter_input_forward[parameter] = value +# return +# end + +function MOI.set( + model::POI.Optimizer, + ::ForwardConstraintSet, + ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Parameter{T}}, + value::Number, +) where {T} + variable = MOI.VariableIndex(ci.value) + if _is_variable(model, variable) + error("Trying to set a forward parameter sensitivity for a variable") + end + sensitivity_data = _get_sensitivity_data(model) + sensitivity_data.parameter_input_forward[variable] = value + return +end + +function MOI.get( + model::POI.Optimizer, + attr::ForwardVariablePrimal, + variable::MOI.VariableIndex, +) + if _is_parameter(model, variable) + error("Trying to get a forward variable sensitivity for a parameter") + end + return MOI.get(model.optimizer, attr, model.variables[variable]) +end + +# reverse mode + +function _constraint_get_reverse!( + model::POI.Optimizer{T}, + affine_constraint_cache_dict, + ::Type{P}, +) where {T, P <: POI.ParametricAffineFunction} + sensitivity_data = _get_sensitivity_data(model) + for (inner_ci, pf) in affine_constraint_cache_dict + terms = POI.affine_parameter_terms(pf) + if isempty(terms) + continue + end + grad_pf_cte = MOI.constant( + MOI.get( + model.optimizer, + ReverseConstraintFunction(), + inner_ci, + ), + ) + for term in terms + p = term.variable + value = get!(sensitivity_data.parameter_output_backward, p, 0.0) + sensitivity_data.parameter_output_backward[p] = + value + term.coefficient * grad_pf_cte + end + end + return +end + +function _constraint_get_reverse!( + model::POI.Optimizer{T}, + vector_affine_constraint_cache_dict, + ::Type{P}, +) where {T, P <: POI.ParametricVectorAffineFunction} + sensitivity_data = _get_sensitivity_data(model) + for (inner_ci, pf) in vector_affine_constraint_cache_dict + terms = POI.vector_affine_parameter_terms(pf) + if isempty(terms) + continue + end + grad_pf_cte = MOI.constant( + MOI.get( + model.optimizer, + ReverseConstraintFunction(), + inner_ci, + ), + ) + for term in terms + p = term.scalar_term.variable + value = get!(sensitivity_data.parameter_output_backward, p, 0.0) + sensitivity_data.parameter_output_backward[p] = + value + + term.scalar_term.coefficient * grad_pf_cte[term.output_index] + end + end + return +end + +function _constraint_get_reverse!( + model::POI.Optimizer{T}, + quadratic_constraint_cache_dict, + ::Type{P}, +) where {T, P <: POI.ParametricQuadraticFunction} + sensitivity_data = _get_sensitivity_data(model) + for (inner_ci, pf) in quadratic_constraint_cache_dict + p_terms = POI.affine_parameter_terms(pf) + pp_terms = POI.quadratic_parameter_parameter_terms(pf) + pv_terms = POI.quadratic_parameter_variable_terms(pf) + if isempty(p_terms) && isempty(pp_terms) && isempty(pv_terms) + continue + end + grad_pf = MOI.get( + model.optimizer, + ReverseConstraintFunction(), + inner_ci, + ) + grad_pf_cte = MOI.constant(grad_pf) + for term in p_terms + p = term.variable + value = get!(sensitivity_data.parameter_output_backward, p, 0.0) + sensitivity_data.parameter_output_backward[p] = + value + term.coefficient * grad_pf_cte + end + for term in pp_terms + p_1 = term.variable_1 + p_2 = term.variable_2 + value_1 = get!(sensitivity_data.parameter_output_backward, p_1, 0.0) + value_2 = get!(sensitivity_data.parameter_output_backward, p_2, 0.0) + # TODO: why there is no factor of 2 here???? + # ANS: probably because it was SET + sensitivity_data.parameter_output_backward[p_1] = + value_1 + + term.coefficient * grad_pf_cte * + MOI.get(model, MOI.VariablePrimal(), p_2) / + ifelse(term.variable_1 === term.variable_2, 1, 1) + sensitivity_data.parameter_output_backward[p_2] = + value_2 + + term.coefficient * grad_pf_cte * + MOI.get(model, MOI.VariablePrimal(), p_1) / + ifelse(term.variable_1 === term.variable_2, 1, 1) + end + for term in pv_terms + p = term.variable_1 + v = term.variable_2 # check if inner or outer (should be inner) + value = get!(sensitivity_data.parameter_output_backward, p, 0.0) + sensitivity_data.parameter_output_backward[p] = + value + term.coefficient * JuMP.coefficient(grad_pf, v) # * fixed value of the parameter ? + end + end + return +end + +function _affine_objective_get_reverse!(model::POI.Optimizer{T}) where {T} + pf = MOI.get(model, POI.ParametricObjectiveFunction{POI.ParametricAffineFunction{T}}()) + terms = POI.affine_parameter_terms(pf) + if isempty(terms) + return + end + sensitivity_data = _get_sensitivity_data(model) + grad_pf = MOI.get(model.optimizer, ReverseObjectiveFunction()) + grad_pf_cte = MOI.constant(grad_pf) + for term in terms + p = term.variable + value = get!(sensitivity_data.parameter_output_backward, p, 0.0) + sensitivity_data.parameter_output_backward[p] = + value + term.coefficient * grad_pf_cte + end + return +end +function _quadratic_objective_get_reverse!(model::POI.Optimizer{T}) where {T} + pf = MOI.get(model, POI.ParametricObjectiveFunction{POI.ParametricQuadraticFunction{T}}()) + p_terms = POI.affine_parameter_terms(pf) + pp_terms = POI.quadratic_parameter_parameter_terms(pf) + pv_terms = POI.quadratic_parameter_variable_terms(pf) + if isempty(p_terms) && isempty(pp_terms) && isempty(pv_terms) + return + end + sensitivity_data = _get_sensitivity_data(model) + grad_pf = MOI.get(model.optimizer, ReverseObjectiveFunction()) + grad_pf_cte = MOI.constant(grad_pf) + for term in p_terms + p = term.variable + value = get!(sensitivity_data.parameter_output_backward, p, 0.0) + sensitivity_data.parameter_output_backward[p] = + value + term.coefficient * grad_pf_cte + end + for term in pp_terms + p_1 = term.variable_1 + p_2 = term.variable_2 + value_1 = get!(sensitivity_data.parameter_output_backward, p_1, 0.0) + value_2 = get!(sensitivity_data.parameter_output_backward, p_2, 0.0) + sensitivity_data.parameter_output_backward[p_1] = + value_1 + + term.coefficient * grad_pf_cte * + MOI.get(model, MOI.VariablePrimal(), p_2) / + ifelse(term.variable_1 === term.variable_2, 2, 1) + sensitivity_data.parameter_output_backward[p_2] = + value_2 + + term.coefficient * grad_pf_cte * + MOI.get(model, MOI.VariablePrimal(), p_1) / + ifelse(term.variable_1 === term.variable_2, 2, 1) + end + for term in pv_terms + p = term.variable_1 + v = term.variable_2 # check if inner or outer (should be inner) + value = get!(sensitivity_data.parameter_output_backward, p, 0.0) + sensitivity_data.parameter_output_backward[p] = + value + term.coefficient * JuMP.coefficient(grad_pf, v) # * fixed value of the parameter ? + end + return +end + +function reverse_differentiate!(model::POI.Optimizer) + reverse_differentiate!(model.optimizer) + sensitivity_data = _get_sensitivity_data(model) + empty!(sensitivity_data.parameter_output_backward) + sizehint!(sensitivity_data.parameter_output_backward, length(model.parameters)) + ctr_types = MOI.get(model, POI.ListOfParametricConstraintTypesPresent()) + for (F, S, P) in ctr_types + dict = MOI.get(model, POI.DictOfParametricConstraintIndicesAndFunctions{F, S, P}()) + _constraint_get_reverse!(model, dict, P) + end + obj_type = MOI.get(model, POI.ParametricObjectiveType()) + if obj_type <: POI.ParametricAffineFunction + _affine_objective_get_reverse!(model) + elseif obj_type <: POI.ParametricQuadraticFunction + _quadratic_objective_get_reverse!(model) + end + return +end + +function _is_parameter(model::POI.Optimizer{T}, variable::MOI.VariableIndex) where {T} + MOI.is_valid(model, MOI.ConstraintIndex{MOI.VariableIndex,MOI.Parameter{T}}(variable.value)) +end + +function _is_variable(model::POI.Optimizer{T}, variable::MOI.VariableIndex) where {T} + MOI.is_valid(model, variable) && !MOI.is_valid(model, MOI.ConstraintIndex{MOI.VariableIndex,MOI.Parameter{T}}(variable.value)) +end + +function MOI.set( + model::POI.Optimizer, + attr::ReverseVariablePrimal, + variable::MOI.VariableIndex, + value::Number, +) + if _is_parameter(model, variable) + error("Trying to set a backward variable sensitivity for a parameter") + end + MOI.set(model.optimizer, attr, variable, value) + return +end + +# MOI.is_set_by_optimize(::ReverseParameter) = true + +# function MOI.get( +# model::POI.Optimizer, +# ::ReverseParameter, +# variable::MOI.VariableIndex, +# ) +# if _is_variable(model, variable) +# error("Trying to get a backward parameter sensitivity for a variable") +# end +# p = variable +# sensitivity_data = _get_sensitivity_data(model) +# return get(sensitivity_data.parameter_output_backward, p, 0.0) +# end + +MOI.is_set_by_optimize(::ReverseConstraintSet) = true + +function MOI.get( + model::POI.Optimizer, + ::ReverseConstraintSet, + ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Parameter{T}}, +) where {T} + variable = MOI.VariableIndex(ci.value) + if _is_variable(model, variable) + error("Trying to get a backward parameter sensitivity for a variable") + end + sensitivity_data = _get_sensitivity_data(model) + return get(sensitivity_data.parameter_output_backward, variable, 0.0) +end + +# extras to handle model_dirty + +function MOI.get( + model::JuMP.Model, + attr::ReverseConstraintSet, + var_ref::JuMP.ConstraintRef, +) + JuMP.check_belongs_to_model(var_ref, model) + return _moi_get_result(JuMP.backend(model), attr, JuMP.index(var_ref)) +end + +""" +""" +function set_forward_parameter( + model::JuMP.Model, + variable::JuMP.VariableRef, + value::Number, +) + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(variable), value) +end + +""" +""" +function get_reverse_parameter( + model::JuMP.Model, + variable::JuMP.VariableRef, +) + return MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(variable)) +end + +""" +""" +function set_reverse_variable( + model::JuMP.Model, + variable::JuMP.VariableRef, + value::Number, +) + MOI.set(model, DiffOpt.ReverseVariablePrimal(), variable, value) +end + +""" +""" +function get_forward_variable( + model::JuMP.Model, + variable::JuMP.VariableRef, +) + return MOI.get(model, DiffOpt.ForwardVariablePrimal(), variable) +end diff --git a/test/parameters.jl b/test/parameters.jl new file mode 100644 index 000000000..af2a32292 --- /dev/null +++ b/test/parameters.jl @@ -0,0 +1,468 @@ +# Copyright (c) 2020: Akshay Sharma 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. + +# using Revise + +module TestParameters + +using Test +using JuMP +import DiffOpt +import MathOptInterface as MOI +import HiGHS +import SCS + +function runtests() + for name in names(@__MODULE__; all = true) + if startswith("$name", "test_") + @testset "$(name)" begin + getfield(@__MODULE__, name)() + end + end + end + return +end + +function test_diff_rhs() + model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer, with_parametric_opt_interface = true)) + set_silent(model) + @variable(model, x) + @variable(model, p in Parameter(3.0)) + @constraint(model, cons, x >= 3 * p) + @objective(model, Min, 2x) + optimize!(model) + @test value(x) ≈ 9 + # the function is + # x(p) = 3p, hence x'(p) = 3 + # differentiate w.r.t. p + # MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), 1) + DiffOpt.set_forward_parameter(model, p, 1) + DiffOpt.forward_differentiate!(model) + # @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 3 + @test DiffOpt.get_forward_variable(model, x) ≈ 3 + # again with different "direction" + # MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), 2) + DiffOpt.set_forward_parameter(model, p, 2) + DiffOpt.forward_differentiate!(model) + # @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 6 + @test DiffOpt.get_forward_variable(model, x) ≈ 6 + # + set_parameter_value(p, 2.0) + optimize!(model) + @test value(x) ≈ 6 + # differentiate w.r.t. p + # MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), 1) + DiffOpt.set_forward_parameter(model, p, 1) + DiffOpt.forward_differentiate!(model) + # @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 3 + @test DiffOpt.get_forward_variable(model, x) ≈ 3 + # again with different "direction" + # MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), 2) + DiffOpt.set_forward_parameter(model, p, 2) + DiffOpt.forward_differentiate!(model) + # @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 6 + @test DiffOpt.get_forward_variable(model, x) ≈ 6 + # + # test reverse mode + # + # MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, 1) + DiffOpt.set_reverse_variable(model, x, 1) + DiffOpt.reverse_differentiate!(model) + # @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) ≈ 3 + @test DiffOpt.get_reverse_parameter(model, p) ≈ 3 + # again with different "direction" + # MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, 2) + DiffOpt.set_reverse_variable(model, x, 2) + DiffOpt.reverse_differentiate!(model) + # @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) ≈ 6 + @test DiffOpt.get_reverse_parameter(model, p) ≈ 6 + return +end + +function test_diff_vector_rhs() + model = direct_model(DiffOpt.diff_optimizer(SCS.Optimizer, with_parametric_opt_interface = true)) + set_silent(model) + @variable(model, x) + @variable(model, p in Parameter(3.0)) + @constraint(model, cons, [x - 3 * p] in MOI.Zeros(1)) + + # FIXME + @constraint(model, fake_soc, [0, 0, 0] in SecondOrderCone()) + + @objective(model, Min, 2x) + optimize!(model) + @test isapprox(value(x), 9, atol = 1e-3) + # the function is + # x(p) = 3p, hence x'(p) = 3 + # differentiate w.r.t. p + for p_val in 0:3 + set_parameter_value(p, p_val) + optimize!(model) + @test isapprox( + value(x), + 3 * p_val, + atol = 1e-3, + ) + for direction in 0:3 + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction) + DiffOpt.forward_differentiate!(model) + @test isapprox( + MOI.get(model, DiffOpt.ForwardVariablePrimal(), x), + direction * 3, + atol = 1e-3, + ) + # reverse mode + MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, direction) + DiffOpt.reverse_differentiate!(model) + @test isapprox( + MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)), + direction * 3, + atol = 1e-3, + ) + end + end + return +end + +function test_affine_changes() + model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer, with_parametric_opt_interface = true)) + direct_model + set_silent(model) + p_val = 3.0 + pc_val = 1.0 + @variable(model, x) + @variable(model, p in Parameter(p_val)) + @variable(model, pc in Parameter(pc_val)) + @constraint(model, cons, pc * x >= 3 * p) + @objective(model, Min, 2x) + optimize!(model) + @test value(x) ≈ 3 * p_val / pc_val + # the function is + # x(p, pc) = 3p / pc, hence dx/dp = 3 / pc, dx/dpc = -3p / pc^2 + # differentiate w.r.t. p + for direction_p in 1:1#2 + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) + DiffOpt.forward_differentiate!(model) + @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ + direction_p * 3 / pc_val + end + # update p + p_val = 2.0 + set_parameter_value(p, p_val) + optimize!(model) + @test value(x) ≈ 3 * p_val / pc_val + # differentiate w.r.t. p + for direction_p in 1:2 + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) + DiffOpt.forward_differentiate!(model) + @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ + direction_p * 3 / pc_val + end + # differentiate w.r.t. pc + # stop differentiating with respect to p + direction_p = 0.0 + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) + for direction_pc in 1:2 + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(pc), direction_pc) + DiffOpt.forward_differentiate!(model) + @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ + -direction_pc * 3 * p_val / pc_val^2 + end + # update pc + pc_val = 2.0 + set_parameter_value(pc, pc_val) + optimize!(model) + @test value(x) ≈ 3 * p_val / pc_val + for direction_pc in 1:1#2 + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(pc), direction_pc) + DiffOpt.forward_differentiate!(model) + @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ + -direction_pc * 3 * p_val / pc_val^2 + end + # test combines directions + for direction_pc in 1:2, direction_p in 1:2 + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(pc), direction_pc) + DiffOpt.forward_differentiate!(model) + @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ + -direction_pc * 3 * p_val / pc_val^2 + direction_p * 3 / pc_val + end + return +end + +function test_affine_changes_compact() + model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer, with_parametric_opt_interface = true)) + set_silent(model) + p_val = 3.0 + pc_val = 1.0 + @variable(model, x) + @variable(model, p in Parameter(p_val)) + @variable(model, pc in Parameter(pc_val)) + @constraint(model, cons, pc * x >= 3 * p) + @objective(model, Min, 2x) + # the function is + # x(p, pc) = 3p / pc, hence dx/dp = 3 / pc, dx/dpc = -3p / pc^2 + for p_val in 1:3, pc_val in 1:3 + set_parameter_value(p, p_val) + set_parameter_value(pc, pc_val) + optimize!(model) + @test value(x) ≈ 3 * p_val / pc_val + for direction_pc in 0:2, direction_p in 0:2 + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(pc), direction_pc) + DiffOpt.forward_differentiate!(model) + @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ + -direction_pc * 3 * p_val / pc_val^2 + + direction_p * 3 / pc_val + end + for direction_x in 0:2 + MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, direction_x) + DiffOpt.reverse_differentiate!(model) + @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) ≈ + direction_x * 3 / pc_val + @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(pc)) ≈ + -direction_x * 3 * p_val / pc_val^2 + end + end + return +end + +function test_quadratic_rhs_changes() + model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer, with_parametric_opt_interface = true)) + set_silent(model) + p_val = 2.0 + q_val = 2.0 + r_val = 2.0 + s_val = 2.0 + t_val = 2.0 + @variable(model, x) + @variable(model, p in Parameter(p_val)) + @variable(model, q in Parameter(q_val)) + @variable(model, r in Parameter(r_val)) + @variable(model, s in Parameter(s_val)) + @variable(model, t in Parameter(t_val)) + @constraint(model, cons, 11 * t * x >= 1 + 3 * p * q + 5 * r^2 + 7 * s) + @objective(model, Min, 2x) + # the function is + # x(p, q, r, s, t) = (1 + 3pq + 5r^2 + 7s) / (11t) + # hence + # dx/dp = 3q / (11t) + # dx/dq = 3p / (11t) + # dx/dr = 10r / (11t) + # dx/ds = 7 / (11t) + # dx/dt = - (1 + 3pq + 5r^2 + 7s) / (11t^2) + optimize!(model) + @test value(x) ≈ + (1 + 3 * p_val * q_val + 5 * r_val^2 + 7 * s_val) / (11 * t_val) + for p_val in 2:3, q_val in 2:3, r_val in 2:3, s_val in 2:3, t_val in 2:3 + set_parameter_value(p, p_val) + set_parameter_value(q, q_val) + set_parameter_value(r, r_val) + set_parameter_value(s, s_val) + set_parameter_value(t, t_val) + optimize!(model) + @test value(x) ≈ + (1 + 3 * p_val * q_val + 5 * r_val^2 + 7 * s_val) / (11 * t_val) + for dir_p in 0:2, dir_q in 0:2, dir_r in 0:2, dir_s in 0:2, dir_t in 0:2 + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), dir_p) + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(q), dir_q) + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(r), dir_r) + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(s), dir_s) + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(t), dir_t) + DiffOpt.forward_differentiate!(model) + @test isapprox( + MOI.get(model, DiffOpt.ForwardVariablePrimal(), x), + dir_p * 3 * q_val / (11 * t_val) + + dir_q * 3 * p_val / (11 * t_val) + + dir_r * 10 * r_val / (11 * t_val) + + dir_s * 7 / (11 * t_val) + + dir_t * ( + -(1 + 3 * p_val * q_val + 5 * r_val^2 + 7 * s_val) / + (11 * t_val^2) + ), + atol = 1e-10, + ) + end + for dir_x in 0:3 + MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, dir_x) + DiffOpt.reverse_differentiate!(model) + @test isapprox( + MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)), + dir_x * 3 * q_val / (11 * t_val), + atol = 1e-10, + ) + @test isapprox( + MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(q)), + dir_x * 3 * p_val / (11 * t_val), + atol = 1e-10, + ) + @test isapprox( + MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(r)), + dir_x * 10 * r_val / (11 * t_val), + atol = 1e-10, + ) + @test isapprox( + MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(s)), + dir_x * 7 / (11 * t_val), + atol = 1e-10, + ) + @test isapprox( + MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(t)), + dir_x * ( + -(1 + 3 * p_val * q_val + 5 * r_val^2 + 7 * s_val) / + (11 * t_val^2) + ), + atol = 1e-10, + ) + end + end + return +end + +function test_affine_changes_compact_max() + model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer, with_parametric_opt_interface = true)) + set_silent(model) + p_val = 3.0 + pc_val = 1.0 + @variable(model, x) + @variable(model, p in Parameter(p_val)) + @variable(model, pc in Parameter(pc_val)) + @constraint(model, cons, pc * x >= 3 * p) + @objective(model, Max, -2x) + # the function is + # x(p, pc) = 3p / pc, hence dx/dp = 3 / pc, dx/dpc = -3p / pc^2 + for p_val in 1:3, pc_val in 1:3 + set_parameter_value(p, p_val) + set_parameter_value(pc, pc_val) + optimize!(model) + @test value(x) ≈ 3 * p_val / pc_val + for direction_pc in 0:2, direction_p in 0:2 + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(pc), direction_pc) + DiffOpt.forward_differentiate!(model) + @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ + -direction_pc * 3 * p_val / pc_val^2 + + direction_p * 3 / pc_val + end + end + return +end + +function test_diff_affine_objective() + model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer, with_parametric_opt_interface = true)) + set_silent(model) + p_val = 3.0 + @variable(model, x) + @variable(model, p in Parameter(p_val)) + @constraint(model, cons, x >= 3) + @objective(model, Min, 2x + 3p) + # x(p, pc) = 3, hence dx/dp = 0 + for p_val in 1:2 + set_parameter_value(p, p_val) + optimize!(model) + @test value(x) ≈ 3 + for direction_p in 0:2 + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) + DiffOpt.forward_differentiate!(model) + @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 0.0 + # reverse mode + MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, direction_p) + DiffOpt.reverse_differentiate!(model) + @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) ≈ 0.0 + end + end + return +end + +function test_diff_quadratic_objective() + model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer, with_parametric_opt_interface = true)) + set_silent(model) + p_val = 3.0 + @variable(model, x) + @variable(model, p in Parameter(p_val)) + @constraint(model, cons, x >= 3) + @objective(model, Min, p * x) + # x(p, pc) = 3, hence dx/dp = 0 + for p_val in 1:2 + set_parameter_value(p, p_val) + optimize!(model) + @test value(x) ≈ 3 + for direction_p in 0:2 + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) + DiffOpt.forward_differentiate!(model) + @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 0.0 + # reverse mode + MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, direction_p) + DiffOpt.reverse_differentiate!(model) + @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) ≈ 0.0 + end + end + return +end + +function test_quadratic_objective_qp() + model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer, with_parametric_opt_interface = true)) + set_silent(model) + p_val = 3.0 + @variable(model, x) + @variable(model, p in Parameter(p_val)) + @constraint(model, cons, x >= -10) + @objective(model, Min, 3 * p * x + x * x + 5 * p + 7 * p^2) + # 2x + 3p = 0, hence x = -3p/2 + # hence dx/dp = -3/2 + for p_val in 3:3 + set_parameter_value(p, p_val) + optimize!(model) + @test value(x) ≈ -3p_val / 2 atol = 1e-4 + for direction_p in 0:2 + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) + DiffOpt.forward_differentiate!(model) + @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ + direction_p * (-3 / 2) atol = 1e-4 + # reverse mode + MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, direction_p) + DiffOpt.reverse_differentiate!(model) + @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) ≈ + direction_p * (-3 / 2) atol = 1e-4 + end + end + return +end + +function test_diff_errors() + model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer, with_parametric_opt_interface = true)) + set_silent(model) + @variable(model, x) + @variable(model, p in Parameter(3.0)) + @constraint(model, cons, x >= 3 * p) + @objective(model, Min, 2x) + optimize!(model) + @test value(x) ≈ 9 + + @test_throws ErrorException MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(x), 1) + @test_throws ErrorException MOI.set( + model, + DiffOpt.ReverseVariablePrimal(), + p, + 1, + ) + @test_throws ErrorException MOI.get( + model, + DiffOpt.ForwardVariablePrimal(), + p, + ) + @test_throws ErrorException MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(x)) + + @test_throws ErrorException MOI.set(model, DiffOpt.ForwardObjectiveFunction(), 3 * x) + @test_throws ErrorException MOI.set(model, DiffOpt.ForwardConstraintFunction(), cons, 1 + 7 * x) + @test_throws ErrorException MOI.get(model, DiffOpt.ReverseObjectiveFunction()) + @test_throws ErrorException MOI.get(model, DiffOpt.ReverseConstraintFunction(), cons) + + return +end + +end # module + +TestParameters.runtests() diff --git a/test/runtests.jl b/test/runtests.jl index 70844576e..626d5eba4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,9 @@ # 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. +# TODO: remove this +Pkg.add(url="https://github.com/jump-dev/ParametricOptInterface.jl", rev="jg/diff2") + using Test @testset "$file" for file in readdir(@__DIR__) From 381581946c391a1115f4fc7c825d4ef54fe5ee53 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Sat, 28 Dec 2024 03:48:34 -0300 Subject: [PATCH 02/26] add missing import --- test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/runtests.jl b/test/runtests.jl index 626d5eba4..704250644 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,6 +4,7 @@ # in the LICENSE.md file or at https://opensource.org/licenses/MIT. # TODO: remove this +import Pkg Pkg.add(url="https://github.com/jump-dev/ParametricOptInterface.jl", rev="jg/diff2") using Test From a1e212c0ff7f466526ab1e9cd8203784ed9e9561 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Sat, 28 Dec 2024 04:15:43 -0300 Subject: [PATCH 03/26] temp change to proj toml --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 2cc3a7687..864cbfe95 100644 --- a/Project.toml +++ b/Project.toml @@ -14,6 +14,7 @@ MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" MathOptSetDistances = "3b969827-a86c-476c-9527-bb6f1a8fbad5" ParametricOptInterface = "0ce4ce61-57bf-432b-a095-efac525d185e" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" [compat] BlockDiagonals = "0.1" From a636915608562b924afe520f1915136d74ba398e Mon Sep 17 00:00:00 2001 From: joaquimg Date: Sat, 28 Dec 2024 04:18:59 -0300 Subject: [PATCH 04/26] format --- src/moi_wrapper.jl | 15 ++- src/parameters.jl | 169 ++++++++++++++++----------- test/parameters.jl | 281 +++++++++++++++++++++++++++++++++++++-------- test/runtests.jl | 5 +- 4 files changed, 345 insertions(+), 125 deletions(-) diff --git a/src/moi_wrapper.jl b/src/moi_wrapper.jl index d577765a7..1a40b26f9 100644 --- a/src/moi_wrapper.jl +++ b/src/moi_wrapper.jl @@ -58,8 +58,10 @@ function diff_optimizer( with_bridge_type = Float64, with_cache::Bool = true, ) - optimizer = - MOI.instantiate(optimizer_constructor; with_bridge_type = with_bridge_type) + optimizer = MOI.instantiate( + optimizer_constructor; + with_bridge_type = with_bridge_type, + ) # When we do `MOI.copy_to(diff, optimizer)` we need to efficiently `MOI.get` # the model information from `optimizer`. However, 1) `optimizer` may not # implement some getters or it may be inefficient and 2) the getters may be @@ -74,9 +76,9 @@ function diff_optimizer( optimizer end if with_parametric_opt_interface - return POI.Optimizer(Optimizer(caching_opt, method = method)) + return POI.Optimizer(Optimizer(caching_opt; method = method)) else - return Optimizer(caching_opt, method = method) + return Optimizer(caching_opt; method = method) end end @@ -93,7 +95,10 @@ mutable struct Optimizer{OT<:MOI.ModelLike} <: MOI.AbstractOptimizer # sensitivity input cache using MOI like sparse format input_cache::InputCache - function Optimizer(optimizer::OT; method = DIFF_AUTOMATIC) where {OT<:MOI.ModelLike} + function Optimizer( + optimizer::OT; + method = DIFF_AUTOMATIC, + ) where {OT<:MOI.ModelLike} output = new{OT}(optimizer, Any[], nothing, nothing, nothing, InputCache()) if method == DIFF_CONIC diff --git a/src/parameters.jl b/src/parameters.jl index d35276fbf..47ee260c6 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -7,12 +7,8 @@ MOI.supports(::POI.Optimizer, ::ForwardObjectiveFunction) = false -function MOI.set( - ::POI.Optimizer, - ::ForwardObjectiveFunction, - _, -) - error( +function MOI.set(::POI.Optimizer, ::ForwardObjectiveFunction, _) + return error( "Forward objective function is not supported when " * "`with_parametric_opt_interface` is set to `true` in " * "`diff_optimizer`." * @@ -28,7 +24,7 @@ function MOI.set( ::MOI.ConstraintIndex, _, ) - error( + return error( "Forward constraint function is not supported when " * "`with_parametric_opt_interface` is set to `true` in " * "`diff_optimizer`." * @@ -38,11 +34,8 @@ end MOI.supports(::POI.Optimizer, ::ReverseObjectiveFunction) = false -function MOI.get( - ::POI.Optimizer, - ::ReverseObjectiveFunction, -) - error( +function MOI.get(::POI.Optimizer, ::ReverseObjectiveFunction) + return error( "Reverse objective function is not supported when " * "`with_parametric_opt_interface` is set to `true` in " * "`diff_optimizer`." * @@ -57,7 +50,7 @@ function MOI.get( ::ReverseConstraintFunction, ::MOI.ConstraintIndex, ) - error( + return error( "Reverse constraint function is not supported when " * "`with_parametric_opt_interface` is set to `true` in " * "`diff_optimizer`." * @@ -77,10 +70,15 @@ mutable struct SensitivityData{T} end function SensitivityData{T}() where {T} - return SensitivityData{T}(Dict{MOI.VariableIndex,T}(), Dict{MOI.VariableIndex,T}()) + return SensitivityData{T}( + Dict{MOI.VariableIndex,T}(), + Dict{MOI.VariableIndex,T}(), + ) end -function _get_sensitivity_data(model::POI.Optimizer{T})::SensitivityData{T} where {T} +function _get_sensitivity_data( + model::POI.Optimizer{T}, +)::SensitivityData{T} where {T} _initialize_sensitivity_data!(model) return model.ext[:_sensitivity_data]::SensitivityData{T} end @@ -100,7 +98,7 @@ function _constraint_set_forward!( model::POI.Optimizer{T}, affine_constraint_cache_dict, ::Type{P}, -) where {T, P <: POI.ParametricAffineFunction} +) where {T,P<:POI.ParametricAffineFunction} sensitivity_data = _get_sensitivity_data(model) for (inner_ci, pf) in affine_constraint_cache_dict cte = zero(T) @@ -127,7 +125,7 @@ function _constraint_set_forward!( model::POI.Optimizer{T}, vector_affine_constraint_cache_dict, ::Type{P}, -) where {T, P <: POI.ParametricVectorAffineFunction} +) where {T,P<:POI.ParametricVectorAffineFunction} sensitivity_data = _get_sensitivity_data(model) for (inner_ci, pf) in vector_affine_constraint_cache_dict cte = zeros(T, length(pf.c)) @@ -154,7 +152,7 @@ function _constraint_set_forward!( model::POI.Optimizer{T}, quadratic_constraint_cache_dict, ::Type{P}, -) where {T, P <: POI.ParametricQuadraticFunction} +) where {T,P<:POI.ParametricQuadraticFunction} sensitivity_data = _get_sensitivity_data(model) for (inner_ci, pf) in quadratic_constraint_cache_dict cte = zero(T) @@ -167,14 +165,18 @@ function _constraint_set_forward!( for term in POI.quadratic_parameter_parameter_terms(pf) p_1 = term.variable_1 p_2 = term.variable_2 - sensitivity_1 = get(sensitivity_data.parameter_input_forward, p_1, 0.0) - sensitivity_2 = get(sensitivity_data.parameter_input_forward, p_2, 0.0) + sensitivity_1 = + get(sensitivity_data.parameter_input_forward, p_1, 0.0) + sensitivity_2 = + get(sensitivity_data.parameter_input_forward, p_2, 0.0) cte += - sensitivity_1 * term.coefficient * + sensitivity_1 * + term.coefficient * MOI.get(model, MOI.VariablePrimal(), p_2) / ifelse(term.variable_1 === term.variable_2, 2, 1) cte += - sensitivity_2 * term.coefficient * + sensitivity_2 * + term.coefficient * MOI.get(model, MOI.VariablePrimal(), p_1) / ifelse(term.variable_1 === term.variable_2, 2, 1) end @@ -227,7 +229,10 @@ end function _quadratic_objective_set_forward!(model::POI.Optimizer{T}) where {T} cte = zero(T) - pf = MOI.get(model, POI.ParametricObjectiveFunction{POI.ParametricQuadraticFunction{T}}()) + pf = MOI.get( + model, + POI.ParametricObjectiveFunction{POI.ParametricQuadraticFunction{T}}(), + ) sensitivity_data = _get_sensitivity_data(model) for term in POI.affine_parameter_terms(pf) p = term.variable @@ -240,13 +245,13 @@ function _quadratic_objective_set_forward!(model::POI.Optimizer{T}) where {T} sensitivity_1 = get(sensitivity_data.parameter_input_forward, p_1, 0.0) sensitivity_2 = get(sensitivity_data.parameter_input_forward, p_2, 0.0) cte += - sensitivity_1 * term.coefficient * + sensitivity_1 * + term.coefficient * MOI.get(model, MOI.VariablePrimal(), p_2) / ifelse(term.variable_1 === term.variable_2, 2, 1) - cte += - sensitivity_2 * term.coefficient - MOI.get(model, MOI.VariablePrimal(), p_1) / - ifelse(term.variable_1 === term.variable_2, 2, 1) + cte += sensitivity_2 * term.coefficient + MOI.get(model, MOI.VariablePrimal(), p_1) / + ifelse(term.variable_1 === term.variable_2, 2, 1) end terms = MOI.ScalarAffineTerm{T}[] sizehint!(terms, length(POI.quadratic_parameter_variable_terms(pf))) @@ -282,7 +287,10 @@ function forward_differentiate!(model::POI.Optimizer{T}) where {T} empty_input_sensitivities!(model) ctr_types = MOI.get(model, POI.ListOfParametricConstraintTypesPresent()) for (F, S, P) in ctr_types - dict = MOI.get(model, POI.DictOfParametricConstraintIndicesAndFunctions{F, S, P}()) + dict = MOI.get( + model, + POI.DictOfParametricConstraintIndicesAndFunctions{F,S,P}(), + ) _constraint_set_forward!(model, dict, P) end obj_type = MOI.get(model, POI.ParametricObjectiveType()) @@ -342,7 +350,7 @@ function _constraint_get_reverse!( model::POI.Optimizer{T}, affine_constraint_cache_dict, ::Type{P}, -) where {T, P <: POI.ParametricAffineFunction} +) where {T,P<:POI.ParametricAffineFunction} sensitivity_data = _get_sensitivity_data(model) for (inner_ci, pf) in affine_constraint_cache_dict terms = POI.affine_parameter_terms(pf) @@ -350,11 +358,7 @@ function _constraint_get_reverse!( continue end grad_pf_cte = MOI.constant( - MOI.get( - model.optimizer, - ReverseConstraintFunction(), - inner_ci, - ), + MOI.get(model.optimizer, ReverseConstraintFunction(), inner_ci), ) for term in terms p = term.variable @@ -370,7 +374,7 @@ function _constraint_get_reverse!( model::POI.Optimizer{T}, vector_affine_constraint_cache_dict, ::Type{P}, -) where {T, P <: POI.ParametricVectorAffineFunction} +) where {T,P<:POI.ParametricVectorAffineFunction} sensitivity_data = _get_sensitivity_data(model) for (inner_ci, pf) in vector_affine_constraint_cache_dict terms = POI.vector_affine_parameter_terms(pf) @@ -378,11 +382,7 @@ function _constraint_get_reverse!( continue end grad_pf_cte = MOI.constant( - MOI.get( - model.optimizer, - ReverseConstraintFunction(), - inner_ci, - ), + MOI.get(model.optimizer, ReverseConstraintFunction(), inner_ci), ) for term in terms p = term.scalar_term.variable @@ -399,7 +399,7 @@ function _constraint_get_reverse!( model::POI.Optimizer{T}, quadratic_constraint_cache_dict, ::Type{P}, -) where {T, P <: POI.ParametricQuadraticFunction} +) where {T,P<:POI.ParametricQuadraticFunction} sensitivity_data = _get_sensitivity_data(model) for (inner_ci, pf) in quadratic_constraint_cache_dict p_terms = POI.affine_parameter_terms(pf) @@ -408,11 +408,8 @@ function _constraint_get_reverse!( if isempty(p_terms) && isempty(pp_terms) && isempty(pv_terms) continue end - grad_pf = MOI.get( - model.optimizer, - ReverseConstraintFunction(), - inner_ci, - ) + grad_pf = + MOI.get(model.optimizer, ReverseConstraintFunction(), inner_ci) grad_pf_cte = MOI.constant(grad_pf) for term in p_terms p = term.variable @@ -429,12 +426,14 @@ function _constraint_get_reverse!( # ANS: probably because it was SET sensitivity_data.parameter_output_backward[p_1] = value_1 + - term.coefficient * grad_pf_cte * + term.coefficient * + grad_pf_cte * MOI.get(model, MOI.VariablePrimal(), p_2) / ifelse(term.variable_1 === term.variable_2, 1, 1) sensitivity_data.parameter_output_backward[p_2] = value_2 + - term.coefficient * grad_pf_cte * + term.coefficient * + grad_pf_cte * MOI.get(model, MOI.VariablePrimal(), p_1) / ifelse(term.variable_1 === term.variable_2, 1, 1) end @@ -450,7 +449,10 @@ function _constraint_get_reverse!( end function _affine_objective_get_reverse!(model::POI.Optimizer{T}) where {T} - pf = MOI.get(model, POI.ParametricObjectiveFunction{POI.ParametricAffineFunction{T}}()) + pf = MOI.get( + model, + POI.ParametricObjectiveFunction{POI.ParametricAffineFunction{T}}(), + ) terms = POI.affine_parameter_terms(pf) if isempty(terms) return @@ -467,7 +469,10 @@ function _affine_objective_get_reverse!(model::POI.Optimizer{T}) where {T} return end function _quadratic_objective_get_reverse!(model::POI.Optimizer{T}) where {T} - pf = MOI.get(model, POI.ParametricObjectiveFunction{POI.ParametricQuadraticFunction{T}}()) + pf = MOI.get( + model, + POI.ParametricObjectiveFunction{POI.ParametricQuadraticFunction{T}}(), + ) p_terms = POI.affine_parameter_terms(pf) pp_terms = POI.quadratic_parameter_parameter_terms(pf) pv_terms = POI.quadratic_parameter_variable_terms(pf) @@ -490,12 +495,14 @@ function _quadratic_objective_get_reverse!(model::POI.Optimizer{T}) where {T} value_2 = get!(sensitivity_data.parameter_output_backward, p_2, 0.0) sensitivity_data.parameter_output_backward[p_1] = value_1 + - term.coefficient * grad_pf_cte * + term.coefficient * + grad_pf_cte * MOI.get(model, MOI.VariablePrimal(), p_2) / ifelse(term.variable_1 === term.variable_2, 2, 1) sensitivity_data.parameter_output_backward[p_2] = value_2 + - term.coefficient * grad_pf_cte * + term.coefficient * + grad_pf_cte * MOI.get(model, MOI.VariablePrimal(), p_1) / ifelse(term.variable_1 === term.variable_2, 2, 1) end @@ -513,10 +520,16 @@ function reverse_differentiate!(model::POI.Optimizer) reverse_differentiate!(model.optimizer) sensitivity_data = _get_sensitivity_data(model) empty!(sensitivity_data.parameter_output_backward) - sizehint!(sensitivity_data.parameter_output_backward, length(model.parameters)) + sizehint!( + sensitivity_data.parameter_output_backward, + length(model.parameters), + ) ctr_types = MOI.get(model, POI.ListOfParametricConstraintTypesPresent()) for (F, S, P) in ctr_types - dict = MOI.get(model, POI.DictOfParametricConstraintIndicesAndFunctions{F, S, P}()) + dict = MOI.get( + model, + POI.DictOfParametricConstraintIndicesAndFunctions{F,S,P}(), + ) _constraint_get_reverse!(model, dict, P) end obj_type = MOI.get(model, POI.ParametricObjectiveType()) @@ -528,12 +541,25 @@ function reverse_differentiate!(model::POI.Optimizer) return end -function _is_parameter(model::POI.Optimizer{T}, variable::MOI.VariableIndex) where {T} - MOI.is_valid(model, MOI.ConstraintIndex{MOI.VariableIndex,MOI.Parameter{T}}(variable.value)) +function _is_parameter( + model::POI.Optimizer{T}, + variable::MOI.VariableIndex, +) where {T} + return MOI.is_valid( + model, + MOI.ConstraintIndex{MOI.VariableIndex,MOI.Parameter{T}}(variable.value), + ) end -function _is_variable(model::POI.Optimizer{T}, variable::MOI.VariableIndex) where {T} - MOI.is_valid(model, variable) && !MOI.is_valid(model, MOI.ConstraintIndex{MOI.VariableIndex,MOI.Parameter{T}}(variable.value)) +function _is_variable( + model::POI.Optimizer{T}, + variable::MOI.VariableIndex, +) where {T} + return MOI.is_valid(model, variable) && + !MOI.is_valid( + model, + MOI.ConstraintIndex{MOI.VariableIndex,MOI.Parameter{T}}(variable.value), + ) end function MOI.set( @@ -597,16 +623,22 @@ function set_forward_parameter( variable::JuMP.VariableRef, value::Number, ) - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(variable), value) + return MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(variable), + value, + ) end """ """ -function get_reverse_parameter( - model::JuMP.Model, - variable::JuMP.VariableRef, -) - return MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(variable)) +function get_reverse_parameter(model::JuMP.Model, variable::JuMP.VariableRef) + return MOI.get( + model, + DiffOpt.ReverseConstraintSet(), + ParameterRef(variable), + ) end """ @@ -616,14 +648,11 @@ function set_reverse_variable( variable::JuMP.VariableRef, value::Number, ) - MOI.set(model, DiffOpt.ReverseVariablePrimal(), variable, value) + return MOI.set(model, DiffOpt.ReverseVariablePrimal(), variable, value) end """ """ -function get_forward_variable( - model::JuMP.Model, - variable::JuMP.VariableRef, -) +function get_forward_variable(model::JuMP.Model, variable::JuMP.VariableRef) return MOI.get(model, DiffOpt.ForwardVariablePrimal(), variable) end diff --git a/test/parameters.jl b/test/parameters.jl index af2a32292..382b85802 100644 --- a/test/parameters.jl +++ b/test/parameters.jl @@ -26,7 +26,12 @@ function runtests() end function test_diff_rhs() - model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer, with_parametric_opt_interface = true)) + model = Model( + () -> DiffOpt.diff_optimizer( + HiGHS.Optimizer; + with_parametric_opt_interface = true, + ), + ) set_silent(model) @variable(model, x) @variable(model, p in Parameter(3.0)) @@ -82,7 +87,12 @@ function test_diff_rhs() end function test_diff_vector_rhs() - model = direct_model(DiffOpt.diff_optimizer(SCS.Optimizer, with_parametric_opt_interface = true)) + model = direct_model( + DiffOpt.diff_optimizer( + SCS.Optimizer; + with_parametric_opt_interface = true, + ), + ) set_silent(model) @variable(model, x) @variable(model, p in Parameter(3.0)) @@ -100,13 +110,14 @@ function test_diff_vector_rhs() for p_val in 0:3 set_parameter_value(p, p_val) optimize!(model) - @test isapprox( - value(x), - 3 * p_val, - atol = 1e-3, - ) + @test isapprox(value(x), 3 * p_val, atol = 1e-3) for direction in 0:3 - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + direction, + ) DiffOpt.forward_differentiate!(model) @test isapprox( MOI.get(model, DiffOpt.ForwardVariablePrimal(), x), @@ -127,7 +138,12 @@ function test_diff_vector_rhs() end function test_affine_changes() - model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer, with_parametric_opt_interface = true)) + model = Model( + () -> DiffOpt.diff_optimizer( + HiGHS.Optimizer; + with_parametric_opt_interface = true, + ), + ) direct_model set_silent(model) p_val = 3.0 @@ -143,7 +159,12 @@ function test_affine_changes() # x(p, pc) = 3p / pc, hence dx/dp = 3 / pc, dx/dpc = -3p / pc^2 # differentiate w.r.t. p for direction_p in 1:1#2 - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + direction_p, + ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ direction_p * 3 / pc_val @@ -155,7 +176,12 @@ function test_affine_changes() @test value(x) ≈ 3 * p_val / pc_val # differentiate w.r.t. p for direction_p in 1:2 - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + direction_p, + ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ direction_p * 3 / pc_val @@ -165,7 +191,12 @@ function test_affine_changes() direction_p = 0.0 MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) for direction_pc in 1:2 - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(pc), direction_pc) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(pc), + direction_pc, + ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ -direction_pc * 3 * p_val / pc_val^2 @@ -176,15 +207,30 @@ function test_affine_changes() optimize!(model) @test value(x) ≈ 3 * p_val / pc_val for direction_pc in 1:1#2 - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(pc), direction_pc) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(pc), + direction_pc, + ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ -direction_pc * 3 * p_val / pc_val^2 end # test combines directions for direction_pc in 1:2, direction_p in 1:2 - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(pc), direction_pc) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + direction_p, + ) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(pc), + direction_pc, + ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ -direction_pc * 3 * p_val / pc_val^2 + direction_p * 3 / pc_val @@ -193,7 +239,12 @@ function test_affine_changes() end function test_affine_changes_compact() - model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer, with_parametric_opt_interface = true)) + model = Model( + () -> DiffOpt.diff_optimizer( + HiGHS.Optimizer; + with_parametric_opt_interface = true, + ), + ) set_silent(model) p_val = 3.0 pc_val = 1.0 @@ -210,8 +261,18 @@ function test_affine_changes_compact() optimize!(model) @test value(x) ≈ 3 * p_val / pc_val for direction_pc in 0:2, direction_p in 0:2 - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(pc), direction_pc) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + direction_p, + ) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(pc), + direction_pc, + ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ -direction_pc * 3 * p_val / pc_val^2 + @@ -220,17 +281,28 @@ function test_affine_changes_compact() for direction_x in 0:2 MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, direction_x) DiffOpt.reverse_differentiate!(model) - @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) ≈ - direction_x * 3 / pc_val - @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(pc)) ≈ - -direction_x * 3 * p_val / pc_val^2 + @test MOI.get( + model, + DiffOpt.ReverseConstraintSet(), + ParameterRef(p), + ) ≈ direction_x * 3 / pc_val + @test MOI.get( + model, + DiffOpt.ReverseConstraintSet(), + ParameterRef(pc), + ) ≈ -direction_x * 3 * p_val / pc_val^2 end end return end function test_quadratic_rhs_changes() - model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer, with_parametric_opt_interface = true)) + model = Model( + () -> DiffOpt.diff_optimizer( + HiGHS.Optimizer; + with_parametric_opt_interface = true, + ), + ) set_silent(model) p_val = 2.0 q_val = 2.0 @@ -266,11 +338,36 @@ function test_quadratic_rhs_changes() @test value(x) ≈ (1 + 3 * p_val * q_val + 5 * r_val^2 + 7 * s_val) / (11 * t_val) for dir_p in 0:2, dir_q in 0:2, dir_r in 0:2, dir_s in 0:2, dir_t in 0:2 - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), dir_p) - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(q), dir_q) - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(r), dir_r) - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(s), dir_s) - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(t), dir_t) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + dir_p, + ) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(q), + dir_q, + ) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(r), + dir_r, + ) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(s), + dir_s, + ) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(t), + dir_t, + ) DiffOpt.forward_differentiate!(model) @test isapprox( MOI.get(model, DiffOpt.ForwardVariablePrimal(), x), @@ -322,7 +419,12 @@ function test_quadratic_rhs_changes() end function test_affine_changes_compact_max() - model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer, with_parametric_opt_interface = true)) + model = Model( + () -> DiffOpt.diff_optimizer( + HiGHS.Optimizer; + with_parametric_opt_interface = true, + ), + ) set_silent(model) p_val = 3.0 pc_val = 1.0 @@ -339,8 +441,18 @@ function test_affine_changes_compact_max() optimize!(model) @test value(x) ≈ 3 * p_val / pc_val for direction_pc in 0:2, direction_p in 0:2 - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(pc), direction_pc) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + direction_p, + ) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(pc), + direction_pc, + ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ -direction_pc * 3 * p_val / pc_val^2 + @@ -351,7 +463,12 @@ function test_affine_changes_compact_max() end function test_diff_affine_objective() - model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer, with_parametric_opt_interface = true)) + model = Model( + () -> DiffOpt.diff_optimizer( + HiGHS.Optimizer; + with_parametric_opt_interface = true, + ), + ) set_silent(model) p_val = 3.0 @variable(model, x) @@ -364,20 +481,34 @@ function test_diff_affine_objective() optimize!(model) @test value(x) ≈ 3 for direction_p in 0:2 - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + direction_p, + ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 0.0 # reverse mode MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, direction_p) DiffOpt.reverse_differentiate!(model) - @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) ≈ 0.0 + @test MOI.get( + model, + DiffOpt.ReverseConstraintSet(), + ParameterRef(p), + ) ≈ 0.0 end end return end function test_diff_quadratic_objective() - model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer, with_parametric_opt_interface = true)) + model = Model( + () -> DiffOpt.diff_optimizer( + HiGHS.Optimizer; + with_parametric_opt_interface = true, + ), + ) set_silent(model) p_val = 3.0 @variable(model, x) @@ -390,20 +521,34 @@ function test_diff_quadratic_objective() optimize!(model) @test value(x) ≈ 3 for direction_p in 0:2 - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + direction_p, + ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 0.0 # reverse mode MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, direction_p) DiffOpt.reverse_differentiate!(model) - @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) ≈ 0.0 + @test MOI.get( + model, + DiffOpt.ReverseConstraintSet(), + ParameterRef(p), + ) ≈ 0.0 end end return end function test_quadratic_objective_qp() - model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer, with_parametric_opt_interface = true)) + model = Model( + () -> DiffOpt.diff_optimizer( + HiGHS.Optimizer; + with_parametric_opt_interface = true, + ), + ) set_silent(model) p_val = 3.0 @variable(model, x) @@ -415,24 +560,37 @@ function test_quadratic_objective_qp() for p_val in 3:3 set_parameter_value(p, p_val) optimize!(model) - @test value(x) ≈ -3p_val / 2 atol = 1e-4 + @test value(x) ≈ -3p_val / 2 atol = 1e-4 for direction_p in 0:2 - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + direction_p, + ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ direction_p * (-3 / 2) atol = 1e-4 # reverse mode MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, direction_p) DiffOpt.reverse_differentiate!(model) - @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) ≈ - direction_p * (-3 / 2) atol = 1e-4 + @test MOI.get( + model, + DiffOpt.ReverseConstraintSet(), + ParameterRef(p), + ) ≈ direction_p * (-3 / 2) atol = 1e-4 end end return end function test_diff_errors() - model = Model(() -> DiffOpt.diff_optimizer(HiGHS.Optimizer, with_parametric_opt_interface = true)) + model = Model( + () -> DiffOpt.diff_optimizer( + HiGHS.Optimizer; + with_parametric_opt_interface = true, + ), + ) set_silent(model) @variable(model, x) @variable(model, p in Parameter(3.0)) @@ -441,7 +599,12 @@ function test_diff_errors() optimize!(model) @test value(x) ≈ 9 - @test_throws ErrorException MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(x), 1) + @test_throws ErrorException MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(x), + 1, + ) @test_throws ErrorException MOI.set( model, DiffOpt.ReverseVariablePrimal(), @@ -453,12 +616,32 @@ function test_diff_errors() DiffOpt.ForwardVariablePrimal(), p, ) - @test_throws ErrorException MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(x)) + @test_throws ErrorException MOI.get( + model, + DiffOpt.ReverseConstraintSet(), + ParameterRef(x), + ) - @test_throws ErrorException MOI.set(model, DiffOpt.ForwardObjectiveFunction(), 3 * x) - @test_throws ErrorException MOI.set(model, DiffOpt.ForwardConstraintFunction(), cons, 1 + 7 * x) - @test_throws ErrorException MOI.get(model, DiffOpt.ReverseObjectiveFunction()) - @test_throws ErrorException MOI.get(model, DiffOpt.ReverseConstraintFunction(), cons) + @test_throws ErrorException MOI.set( + model, + DiffOpt.ForwardObjectiveFunction(), + 3 * x, + ) + @test_throws ErrorException MOI.set( + model, + DiffOpt.ForwardConstraintFunction(), + cons, + 1 + 7 * x, + ) + @test_throws ErrorException MOI.get( + model, + DiffOpt.ReverseObjectiveFunction(), + ) + @test_throws ErrorException MOI.get( + model, + DiffOpt.ReverseConstraintFunction(), + cons, + ) return end diff --git a/test/runtests.jl b/test/runtests.jl index 704250644..95c8d9949 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,7 +5,10 @@ # TODO: remove this import Pkg -Pkg.add(url="https://github.com/jump-dev/ParametricOptInterface.jl", rev="jg/diff2") +Pkg.add(; + url = "https://github.com/jump-dev/ParametricOptInterface.jl", + rev = "jg/diff2", +) using Test From 4448a46e5f59605c7231e955b7c7802a9602aba6 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Sat, 28 Dec 2024 04:49:59 -0300 Subject: [PATCH 05/26] simplify method setting to sue model constructor --- src/DiffOpt.jl | 10 ---------- src/moi_wrapper.jl | 45 +++++---------------------------------------- 2 files changed, 5 insertions(+), 50 deletions(-) diff --git a/src/DiffOpt.jl b/src/DiffOpt.jl index 57c59fb5e..f877ddc5f 100644 --- a/src/DiffOpt.jl +++ b/src/DiffOpt.jl @@ -40,16 +40,6 @@ function add_all_model_constructors(model) return end -function add_conic_model_constructor(model) - add_model_constructor(model, ConicProgram.Model) - return -end - -function add_quadratic_model_constructor(model) - add_model_constructor(model, QuadraticProgram.Model) - return -end - export diff_optimizer end # module diff --git a/src/moi_wrapper.jl b/src/moi_wrapper.jl index 1a40b26f9..9f15a20c8 100644 --- a/src/moi_wrapper.jl +++ b/src/moi_wrapper.jl @@ -3,36 +3,6 @@ # 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. -""" - DiffMethod - -An enum to define the differentiation method. - -## Values - -Possible values are: - - * [`DIFF_AUTOMATIC`]: Automatic differentiation: tries all differentiation methods and chooses the first that works. This can be slower than manually choosing a method. - * [`DIFF_CONIC`]: Conic optimization based differentiation: works for conic programs or programs that can be transformed into conic programs by JuMP. - * [`DIFF_QUADRATIC`]: Quadratic optimization based differentiation: works for quadratic programs or programs that can be transformed into conic programs by JuMP. -""" -@enum(DiffMethod, DIFF_AUTOMATIC, DIFF_CONIC, DIFF_QUADRATIC) - -@doc( - "Automatic differentiation: tries all differentiation methods and chooses the first that works. This can be slower than manually choosing a method.", - DIFF_AUTOMATIC -) - -@doc( - "Conic optimization based differentiation: works for conic programs or programs that can be transformed into conic programs by JuMP.", - DIFF_CONIC -) - -@doc( - "Quadratic optimization based differentiation: works for quadratic programs or programs that can be transformed into conic programs by JuMP.", - DIFF_QUADRATIC -) - """ diff_optimizer(optimizer_constructor) @@ -53,7 +23,7 @@ julia> model.add_constraint(model, ...) """ function diff_optimizer( optimizer_constructor; - method = DIFF_AUTOMATIC, + method = nothing, with_parametric_opt_interface::Bool = false, with_bridge_type = Float64, with_cache::Bool = true, @@ -97,18 +67,13 @@ mutable struct Optimizer{OT<:MOI.ModelLike} <: MOI.AbstractOptimizer function Optimizer( optimizer::OT; - method = DIFF_AUTOMATIC, + method = nothing, ) where {OT<:MOI.ModelLike} output = new{OT}(optimizer, Any[], nothing, nothing, nothing, InputCache()) - if method == DIFF_CONIC - add_conic_model_constructor(output) - elseif method == DIFF_QUADRATIC - add_quadratic_model_constructor(output) - elseif method == DIFF_AUTOMATIC - add_all_model_constructors(output) - else - add_model_constructor(output, method) + add_all_model_constructors(output) + if method !== nothing + output.model_constructor = method end return output end From 1da3d2aa4f18a365945bb4d427a8a8467af699e5 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Sat, 28 Dec 2024 05:51:07 -0300 Subject: [PATCH 06/26] add possible fix to scalarize bridge error --- src/bridges.jl | 21 +++++++++++++++++++++ src/utils.jl | 17 +++++++++++++++++ test/parameters.jl | 24 ++++++++++++++++++++++++ 3 files changed, 62 insertions(+) diff --git a/src/bridges.jl b/src/bridges.jl index 02299d65a..b928a6123 100644 --- a/src/bridges.jl +++ b/src/bridges.jl @@ -43,6 +43,27 @@ function MOI.get( MOI.get(model, attr, bridge.vector_constraint), )[1] end + +function MOI.set( + model::MOI.ModelLike, + attr::ForwardConstraintFunction, + bridge::MOI.Bridges.Constraint.ScalarizeBridge{T}, + value, +) where {T} + MOI.set.(model, attr, bridge.scalar_constraints, value) + return +end + +function MOI.get( + model::MOI.ModelLike, + attr::ReverseConstraintFunction, + bridge::MOI.Bridges.Constraint.ScalarizeBridge, +) + return _vectorize( + MOI.get.(model, attr, bridge.scalar_constraints), + ) +end + function MOI.get( model::MOI.ModelLike, attr::DiffOpt.ReverseConstraintFunction, diff --git a/src/utils.jl b/src/utils.jl index 7e0b738c0..0a24f441e 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -3,6 +3,8 @@ # 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. +# Helpers for objective function + # Representation of MOI functions using SparseArrays # Might be able to replace in the future by a function in MOI, see # https://github.com/jump-dev/MathOptInterface.jl/pull/1238 @@ -152,6 +154,8 @@ function sparse_array_representation( ) end +# Helpers for scalar constraints + # In the future, we could replace by https://github.com/jump-dev/MathOptInterface.jl/pull/1238 """ VectorScalarAffineFunction{T, VT} <: MOI.AbstractScalarFunction @@ -228,6 +232,8 @@ function JuMP.coefficient( return JuMP.coefficient(func.affine, vi) end +# Helpers for Vector constraints + """ MatrixVectorAffineFunction{T, VT} <: MOI.AbstractVectorFunction @@ -260,6 +266,17 @@ function standard_form(func::MatrixVectorAffineFunction{T}) where {T} return convert(MOI.VectorAffineFunction{T}, func) end +_get_terms(f::VectorScalarAffineFunction) = f.terms +_get_constant(f::VectorScalarAffineFunction) = f.constant +function _vectorize(vec::Vector{T}) where {T<:VectorScalarAffineFunction} + terms = LazyArrays.@~ LazyArrays.ApplyArray(hcat, _get_terms.(vec)...)' + constants = LazyArrays.ApplyArray(vcat, _get_constant.(vec)...) + AT = typeof(terms) + VT = typeof(constants) + ret = MatrixVectorAffineFunction{AT,VT}(terms, constants) + return ret +end + # Only used for testing at the moment so performance is not critical so # converting to standard form is ok function MOIU.isapprox_zero( diff --git a/test/parameters.jl b/test/parameters.jl index 382b85802..0a54baabe 100644 --- a/test/parameters.jl +++ b/test/parameters.jl @@ -646,6 +646,30 @@ function test_diff_errors() return end +function test_scalarize_bridge() + using JuMP, DiffOpt, HiGHS + + b = [1.0, 2.0] + + m = Model( + () -> DiffOpt.diff_optimizer( + HiGHS.Optimizer; + with_parametric_opt_interface = true, + ), + ) + @variable(m, x[1:2] >= 0) + @variable(m, c[1:2] in MOI.Parameter.(b)) + @constraint(m, con, x <= c) + @objective(m, Max, sum(x)) + optimize!(m) + + MOI.set(m, DiffOpt.ReverseVariablePrimal(), m[:x][1], 1.0) + DiffOpt.reverse_differentiate!(m) + @test MOI.get.(m, DiffOpt.ReverseConstraintSet(), ParameterRef.(c)) == [1, 0] + + return +end + end # module TestParameters.runtests() From 4498738c578b2d44c1b60947b549cd1ea43d82f8 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Sat, 28 Dec 2024 06:08:08 -0300 Subject: [PATCH 07/26] add pkg to project --- test/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/test/Project.toml b/test/Project.toml index fdb058f22..bd7cd5af4 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -16,6 +16,7 @@ SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" [compat] HiGHS = "1" From e9caaa4f1f938d454ced0b39929e43b969d716b0 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Sat, 28 Dec 2024 06:09:09 -0300 Subject: [PATCH 08/26] format --- src/bridges.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/bridges.jl b/src/bridges.jl index b928a6123..ef1bab80a 100644 --- a/src/bridges.jl +++ b/src/bridges.jl @@ -59,9 +59,7 @@ function MOI.get( attr::ReverseConstraintFunction, bridge::MOI.Bridges.Constraint.ScalarizeBridge, ) - return _vectorize( - MOI.get.(model, attr, bridge.scalar_constraints), - ) + return _vectorize(MOI.get.(model, attr, bridge.scalar_constraints)) end function MOI.get( From 43613c9bb552ac3c4839667c33b7fe3a83a228a9 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Sat, 11 Jan 2025 22:14:03 -0300 Subject: [PATCH 09/26] improvements --- README.md | 75 +++++++++++++++++++++++++++++++++++++++++++++- docs/src/manual.md | 24 ++++++++------- src/DiffOpt.jl | 5 ++++ src/parameters.jl | 61 ++----------------------------------- test/parameters.jl | 25 ---------------- test/utils.jl | 14 ++++----- 6 files changed, 102 insertions(+), 102 deletions(-) diff --git a/README.md b/README.md index 36541bab2..1a963244d 100644 --- a/README.md +++ b/README.md @@ -31,7 +31,80 @@ examples, tutorials, and an API reference. ## Use with JuMP -Use DiffOpt with JuMP by following this brief example: +### DiffOpt-JuMP API with `Parameters` + +```julia +using JuMP, DiffOpt, HiGHS + +# model = Model( +# () -> DiffOpt.diff_optimizer( +# HiGHS.Optimizer; +# with_parametric_opt_interface = true, +# ), +# ) +model = DiffOpt.diff_model( + HiGHS.Optimizer; + with_parametric_opt_interface = true, +) +set_silent(model) + +p_val = 4.0 +pc_val = 2.0 +@variable(model, x) +@variable(model, p in Parameter(p_val)) +@variable(model, pc in Parameter(pc_val)) +@constraint(model, cons, pc * x >= 3 * p) #??? InvalidConstraintRef TODO +@objective(model, Min, 2x) +optimize!(model) +@show value(x) == 3 * p_val / pc_val + +# the function is +# x(p, pc) = 3p / pc +# hence, +# dx/dp = 3 / pc +# dx/dpc = -3p / pc^2 + +# First, try forward mode AD + +# differentiate w.r.t. p +direction_p = 3.0 +DiffOpt.set_forward_parameter(model, p, direction_p) +DiffOpt.forward_differentiate!(model) +@show DiffOpt.get_forward_variable(model, x) == direction_p * 3 / pc_val + +# update p and pc +p_val = 2.0 +pc_val = 6.0 +set_parameter_value(p, p_val) +set_parameter_value(pc, pc_val) +# re-optimize +optimize!(model) +# check solution +@show value(x) ≈ 3 * p_val / pc_val + +# stop differentiating with respect to p +DiffOpt.empty_input_sensitivities!(model) +# differentiate w.r.t. pc +direction_pc = 10.0 +DiffOpt.set_forward_parameter(model, pc, direction_pc) +DiffOpt.forward_differentiate!(model) +@show abs(DiffOpt.get_forward_variable(model, x) - + -direction_pc * 3 * p_val / pc_val^2) < 1e-5 + +# always a good practice to clear previously set sensitivities +DiffOpt.empty_input_sensitivities!(model) +# Now, reverse model AD +direction_x = 10.0 +DiffOpt.set_reverse_variable(model, x, direction_x) +DiffOpt.reverse_differentiate!(model) +@show DiffOpt.get_reverse_parameter(model, p) == direction_x * 3 / pc_val +@show abs(DiffOpt.get_reverse_parameter(model, pc) - + -direction_x * 3 * p_val / pc_val^2) < 1e-5 +``` + +### Low level DiffOpt-JuMP API: + +A brief example: ```julia using JuMP, DiffOpt, HiGHS diff --git a/docs/src/manual.md b/docs/src/manual.md index 57f007219..12c9be0a2 100644 --- a/docs/src/manual.md +++ b/docs/src/manual.md @@ -4,9 +4,9 @@ As of now, this package only works for optimization models that can be written either in convex conic form or convex quadratic form. -## Supported objectives & constraints - scheme 1 +## Supported objectives & constraints - `QuadraticProgram` backend -For `QPTH`/`OPTNET` style backend, the package supports following `Function-in-Set` constraints: +For `QuadraticProgram` backend, the package supports following `Function-in-Set` constraints: | MOI Function | MOI Set | |:-------|:---------------| @@ -26,9 +26,9 @@ and the following objective types: | `ScalarQuadraticFunction` | -## Supported objectives & constraints - scheme 2 +## Supported objectives & constraints - `ConicProgram` backend -For `DiffCP`/`CVXPY` style backend, the package supports following `Function-in-Set` constraints: +For the `ConicProgram` backend backend, the package supports following `Function-in-Set` constraints: | MOI Function | MOI Set | |:-------|:---------------| @@ -50,18 +50,22 @@ and the following objective types: | `VariableIndex` | | `ScalarAffineFunction` | +Other conic sets such as `RotatedSecondOrderCone` and `PositiveSemidefiniteConeSquare` are supported through bridges. -## Creating a differentiable optimizer + +## Creating a differentiable MOI optimizer You can create a differentiable optimizer over an existing MOI solver by using the `diff_optimizer` utility. ```@docs diff_optimizer ``` -## Adding new sets and constraints +## Creating a differentiable JuMP model -The DiffOpt `Optimizer` behaves similarly to other MOI Optimizers -and implements the `MOI.AbstractOptimizer` API. +You initialize a differentiable JuMP model by using the `diff_model` utility. +```@docs +diff_model +``` ## Projections on cone sets @@ -104,6 +108,4 @@ In the light of above, DiffOpt differentiates program variables ``x``, ``s``, `` - OptNet: Differentiable Optimization as a Layer in Neural Networks ### Backward Pass vector -One possible point of confusion in finding Jacobians is the role of the backward pass vector - above eqn (7), *OptNet: Differentiable Optimization as a Layer in Neural Networks*. While differentiating convex programs, it is often the case that we don't want to find the actual derivatives, rather we might be interested in computing the product of Jacobians with a *backward pass vector*, often used in backprop in machine learning/automatic differentiation. This is what happens in scheme 1 of `DiffOpt` backend. - -But, for the conic system (scheme 2), we provide perturbations in conic data (`dA`, `db`, `dc`) to compute pertubations (`dx`, `dy`, `dz`) in input variables. Unlike the quadratic case, these perturbations are actual derivatives, not the product with a backward pass vector. This is an important distinction between the two schemes of differential optimization. +One possible point of confusion in finding Jacobians is the role of the backward pass vector - above eqn (7), *OptNet: Differentiable Optimization as a Layer in Neural Networks*. While differentiating convex programs, it is often the case that we don't want to find the actual derivatives, rather we might be interested in computing the product of Jacobians with a *backward pass vector*, often used in backpropagation in machine learning/automatic differentiation. This is what happens in `DiffOpt` backends. diff --git a/src/DiffOpt.jl b/src/DiffOpt.jl index f877ddc5f..628932e54 100644 --- a/src/DiffOpt.jl +++ b/src/DiffOpt.jl @@ -19,6 +19,7 @@ include("utils.jl") include("product_of_sets.jl") include("diff_opt.jl") include("moi_wrapper.jl") +include("jump_wrapper.jl") include("jump_moi_overloads.jl") include("parameters.jl") @@ -42,4 +43,8 @@ end export diff_optimizer + +# TODO +# add precompilation statements + end # module diff --git a/src/parameters.jl b/src/parameters.jl index 47ee260c6..757e98b38 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -278,13 +278,14 @@ function _quadratic_objective_set_forward!(model::POI.Optimizer{T}) where {T} return end -function empty_input_sensitivities!(model::POI.Optimizer) +function empty_input_sensitivities!(model::POI.Optimizer{T}) where {T} empty_input_sensitivities!(model.optimizer) + model.ext[:_sensitivity_data] = SensitivityData{T}() return end function forward_differentiate!(model::POI.Optimizer{T}) where {T} - empty_input_sensitivities!(model) + empty_input_sensitivities!(model.optimizer) ctr_types = MOI.get(model, POI.ListOfParametricConstraintTypesPresent()) for (F, S, P) in ctr_types dict = MOI.get( @@ -575,21 +576,6 @@ function MOI.set( return end -# MOI.is_set_by_optimize(::ReverseParameter) = true - -# function MOI.get( -# model::POI.Optimizer, -# ::ReverseParameter, -# variable::MOI.VariableIndex, -# ) -# if _is_variable(model, variable) -# error("Trying to get a backward parameter sensitivity for a variable") -# end -# p = variable -# sensitivity_data = _get_sensitivity_data(model) -# return get(sensitivity_data.parameter_output_backward, p, 0.0) -# end - MOI.is_set_by_optimize(::ReverseConstraintSet) = true function MOI.get( @@ -615,44 +601,3 @@ function MOI.get( JuMP.check_belongs_to_model(var_ref, model) return _moi_get_result(JuMP.backend(model), attr, JuMP.index(var_ref)) end - -""" -""" -function set_forward_parameter( - model::JuMP.Model, - variable::JuMP.VariableRef, - value::Number, -) - return MOI.set( - model, - DiffOpt.ForwardConstraintSet(), - ParameterRef(variable), - value, - ) -end - -""" -""" -function get_reverse_parameter(model::JuMP.Model, variable::JuMP.VariableRef) - return MOI.get( - model, - DiffOpt.ReverseConstraintSet(), - ParameterRef(variable), - ) -end - -""" -""" -function set_reverse_variable( - model::JuMP.Model, - variable::JuMP.VariableRef, - value::Number, -) - return MOI.set(model, DiffOpt.ReverseVariablePrimal(), variable, value) -end - -""" -""" -function get_forward_variable(model::JuMP.Model, variable::JuMP.VariableRef) - return MOI.get(model, DiffOpt.ForwardVariablePrimal(), variable) -end diff --git a/test/parameters.jl b/test/parameters.jl index 0a54baabe..2802bd39e 100644 --- a/test/parameters.jl +++ b/test/parameters.jl @@ -144,7 +144,6 @@ function test_affine_changes() with_parametric_opt_interface = true, ), ) - direct_model set_silent(model) p_val = 3.0 pc_val = 1.0 @@ -646,30 +645,6 @@ function test_diff_errors() return end -function test_scalarize_bridge() - using JuMP, DiffOpt, HiGHS - - b = [1.0, 2.0] - - m = Model( - () -> DiffOpt.diff_optimizer( - HiGHS.Optimizer; - with_parametric_opt_interface = true, - ), - ) - @variable(m, x[1:2] >= 0) - @variable(m, c[1:2] in MOI.Parameter.(b)) - @constraint(m, con, x <= c) - @objective(m, Max, sum(x)) - optimize!(m) - - MOI.set(m, DiffOpt.ReverseVariablePrimal(), m[:x][1], 1.0) - DiffOpt.reverse_differentiate!(m) - @test MOI.get.(m, DiffOpt.ReverseConstraintSet(), ParameterRef.(c)) == [1, 0] - - return -end - end # module TestParameters.runtests() diff --git a/test/utils.jl b/test/utils.jl index 10e37eacb..49cfac385 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -36,7 +36,7 @@ PMLR, 2017. https://arxiv.org/pdf/1703.00443.pdf """ function qp_test( solver, - diff_model, + _diff_model, lt::Bool, set_zero::Bool, canonicalize::Bool; @@ -88,7 +88,7 @@ function qp_test( atol = ATOL, rtol = RTOL, ) - is_conic_qp = !all(iszero, Q) && diff_model == DiffOpt.ConicProgram.Model + is_conic_qp = !all(iszero, Q) && _diff_model == DiffOpt.ConicProgram.Model n = length(q) @assert n == LinearAlgebra.checksquare(Q) @assert n == size(A, 2) @@ -162,7 +162,7 @@ function qp_test( end @_test(convert(Vector{Float64}, _λ), λ) - MOI.set(model, DiffOpt.ModelConstructor(), diff_model) + MOI.set(model, DiffOpt.ModelConstructor(), _diff_model) #dobjb = v' * (dQb / 2.0) * v + dqb' * v # TODO, it should .- @@ -344,7 +344,7 @@ function qp_test( return end -function qp_test(solver, diff_model; kws...) +function qp_test(solver, _diff_model; kws...) @testset "With $(lt ? "LessThan" : "GreaterThan") constraints" for lt in [ true, #false, @@ -359,7 +359,7 @@ function qp_test(solver, diff_model; kws...) true, #false, ] - qp_test(solver, diff_model, lt, set_zero, canonicalize; kws...) + qp_test(solver, _diff_model, lt, set_zero, canonicalize; kws...) end end end @@ -367,11 +367,11 @@ function qp_test(solver, diff_model; kws...) end function qp_test(solver; kws...) - @testset "With $diff_model" for diff_model in [ + @testset "With $_diff_model" for _diff_model in [ DiffOpt.QuadraticProgram.Model, DiffOpt.ConicProgram.Model, ] - qp_test(solver, diff_model; kws...) + qp_test(solver, _diff_model; kws...) end return end From e66502bd832ce82e14fac2ff910a262dbcfb6dc7 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Sat, 11 Jan 2025 22:14:37 -0300 Subject: [PATCH 10/26] remove jump wrapper --- src/DiffOpt.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/DiffOpt.jl b/src/DiffOpt.jl index 628932e54..edb0acdb2 100644 --- a/src/DiffOpt.jl +++ b/src/DiffOpt.jl @@ -19,7 +19,6 @@ include("utils.jl") include("product_of_sets.jl") include("diff_opt.jl") include("moi_wrapper.jl") -include("jump_wrapper.jl") include("jump_moi_overloads.jl") include("parameters.jl") From 526bce9dfeb066a196be3ba347baa1c620a5ce0a Mon Sep 17 00:00:00 2001 From: joaquimg Date: Sat, 11 Jan 2025 17:17:15 -0800 Subject: [PATCH 11/26] clean tests --- Project.toml | 1 + test/runtests.jl | 7 ------- 2 files changed, 1 insertion(+), 7 deletions(-) diff --git a/Project.toml b/Project.toml index 651fe308e..c66179a25 100644 --- a/Project.toml +++ b/Project.toml @@ -24,4 +24,5 @@ JuMP = "1" LazyArrays = "0.21, 0.22, 1" MathOptInterface = "1.18" MathOptSetDistances = "0.2.7" +ParametricOptInterface = "0.9.0" julia = "1.6" diff --git a/test/runtests.jl b/test/runtests.jl index 95c8d9949..70844576e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,13 +3,6 @@ # 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. -# TODO: remove this -import Pkg -Pkg.add(; - url = "https://github.com/jump-dev/ParametricOptInterface.jl", - rev = "jg/diff2", -) - using Test @testset "$file" for file in readdir(@__DIR__) From db73bb98a977ae061c8ff3fcca17978d826a2292 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Sat, 11 Jan 2025 17:31:43 -0800 Subject: [PATCH 12/26] fix readme --- README.md | 28 ++++++++++++---------------- 1 file changed, 12 insertions(+), 16 deletions(-) diff --git a/README.md b/README.md index e600ae132..b14f13be9 100644 --- a/README.md +++ b/README.md @@ -36,15 +36,11 @@ examples, tutorials, and an API reference. ```julia using JuMP, DiffOpt, HiGHS -# model = Model( -# () -> DiffOpt.diff_optimizer( -# HiGHS.Optimizer; -# with_parametric_opt_interface = true, -# ), -# ) -model = DiffOpt.diff_model( - HiGHS.Optimizer; - with_parametric_opt_interface = true, +model = Model( + () -> DiffOpt.diff_optimizer( + HiGHS.Optimizer; + with_parametric_opt_interface = true, + ), ) set_silent(model) @@ -68,9 +64,9 @@ optimize!(model) # differentiate w.r.t. p direction_p = 3.0 -DiffOpt.set_forward_parameter(model, p, direction_p) +MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) DiffOpt.forward_differentiate!(model) -@show DiffOpt.get_forward_variable(model, x) == direction_p * 3 / pc_val +@show MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) == direction_p * 3 / pc_val # update p and pc p_val = 2.0 @@ -86,19 +82,19 @@ optimize!(model) DiffOpt.empty_input_sensitivities!(model) # differentiate w.r.t. pc direction_pc = 10.0 -DiffOpt.set_forward_parameter(model, pc, direction_pc) +MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(pc), direction_pc) DiffOpt.forward_differentiate!(model) -@show abs(DiffOpt.get_forward_variable(model, x) - +@show abs(MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) - -direction_pc * 3 * p_val / pc_val^2) < 1e-5 # always a good practice to clear previously set sensitivities DiffOpt.empty_input_sensitivities!(model) # Now, reverse model AD direction_x = 10.0 -DiffOpt.set_reverse_variable(model, x, direction_x) +MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, direction_x) DiffOpt.reverse_differentiate!(model) -@show DiffOpt.get_reverse_parameter(model, p) == direction_x * 3 / pc_val -@show abs(DiffOpt.get_reverse_parameter(model, pc) - +@show MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) == direction_x * 3 / pc_val +@show abs(MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(pc)) - -direction_x * 3 * p_val / pc_val^2) < 1e-5 ``` From d5bd8033a9d8314e6322ff9e511c44256143b76b Mon Sep 17 00:00:00 2001 From: joaquimg Date: Sat, 11 Jan 2025 17:53:18 -0800 Subject: [PATCH 13/26] use intermediary API --- test/parameters.jl | 36 ++++++++++++------------------------ 1 file changed, 12 insertions(+), 24 deletions(-) diff --git a/test/parameters.jl b/test/parameters.jl index 2802bd39e..9c901d579 100644 --- a/test/parameters.jl +++ b/test/parameters.jl @@ -42,47 +42,35 @@ function test_diff_rhs() # the function is # x(p) = 3p, hence x'(p) = 3 # differentiate w.r.t. p - # MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), 1) - DiffOpt.set_forward_parameter(model, p, 1) + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), 1) DiffOpt.forward_differentiate!(model) - # @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 3 - @test DiffOpt.get_forward_variable(model, x) ≈ 3 + @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 3 # again with different "direction" - # MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), 2) - DiffOpt.set_forward_parameter(model, p, 2) + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), 2) DiffOpt.forward_differentiate!(model) - # @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 6 - @test DiffOpt.get_forward_variable(model, x) ≈ 6 + @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 6 # set_parameter_value(p, 2.0) optimize!(model) @test value(x) ≈ 6 # differentiate w.r.t. p - # MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), 1) - DiffOpt.set_forward_parameter(model, p, 1) + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), 1) DiffOpt.forward_differentiate!(model) - # @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 3 - @test DiffOpt.get_forward_variable(model, x) ≈ 3 + @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 3 # again with different "direction" - # MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), 2) - DiffOpt.set_forward_parameter(model, p, 2) + MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), 2) DiffOpt.forward_differentiate!(model) - # @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 6 - @test DiffOpt.get_forward_variable(model, x) ≈ 6 + @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 6 # # test reverse mode # - # MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, 1) - DiffOpt.set_reverse_variable(model, x, 1) + MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, 1) DiffOpt.reverse_differentiate!(model) - # @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) ≈ 3 - @test DiffOpt.get_reverse_parameter(model, p) ≈ 3 + @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) ≈ 3 # again with different "direction" - # MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, 2) - DiffOpt.set_reverse_variable(model, x, 2) + MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, 2) DiffOpt.reverse_differentiate!(model) - # @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) ≈ 6 - @test DiffOpt.get_reverse_parameter(model, p) ≈ 6 + @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) ≈ 6 return end From 3f995857dd9d55b5ed48332da4c1c4ee8fa5dfdc Mon Sep 17 00:00:00 2001 From: joaquimg Date: Sat, 11 Jan 2025 18:54:29 -0800 Subject: [PATCH 14/26] format --- src/DiffOpt.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/DiffOpt.jl b/src/DiffOpt.jl index edb0acdb2..f09d26082 100644 --- a/src/DiffOpt.jl +++ b/src/DiffOpt.jl @@ -42,7 +42,6 @@ end export diff_optimizer - # TODO # add precompilation statements From 24347f95ce080e332496a8d66f7dbda03e540127 Mon Sep 17 00:00:00 2001 From: Joaquim Dias Garcia Date: Mon, 13 Jan 2025 09:49:18 -0800 Subject: [PATCH 15/26] Apply suggestions from code review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Benoît Legat --- Project.toml | 1 - docs/src/manual.md | 2 +- src/bridges.jl | 4 ++-- src/moi_wrapper.jl | 2 +- test/Project.toml | 1 - 5 files changed, 4 insertions(+), 6 deletions(-) diff --git a/Project.toml b/Project.toml index c66179a25..32b65d1e1 100644 --- a/Project.toml +++ b/Project.toml @@ -14,7 +14,6 @@ MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" MathOptSetDistances = "3b969827-a86c-476c-9527-bb6f1a8fbad5" ParametricOptInterface = "0ce4ce61-57bf-432b-a095-efac525d185e" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" [compat] BlockDiagonals = "0.1" diff --git a/docs/src/manual.md b/docs/src/manual.md index 12c9be0a2..52edd6112 100644 --- a/docs/src/manual.md +++ b/docs/src/manual.md @@ -28,7 +28,7 @@ and the following objective types: ## Supported objectives & constraints - `ConicProgram` backend -For the `ConicProgram` backend backend, the package supports following `Function-in-Set` constraints: +For the `ConicProgram` backend, the package supports following `Function-in-Set` constraints: | MOI Function | MOI Set | |:-------|:---------------| diff --git a/src/bridges.jl b/src/bridges.jl index ef1bab80a..2b0eb46b4 100644 --- a/src/bridges.jl +++ b/src/bridges.jl @@ -47,9 +47,9 @@ end function MOI.set( model::MOI.ModelLike, attr::ForwardConstraintFunction, - bridge::MOI.Bridges.Constraint.ScalarizeBridge{T}, + bridge::MOI.Bridges.Constraint.ScalarizeBridge, value, -) where {T} +) MOI.set.(model, attr, bridge.scalar_constraints, value) return end diff --git a/src/moi_wrapper.jl b/src/moi_wrapper.jl index 9f15a20c8..f47a83f94 100644 --- a/src/moi_wrapper.jl +++ b/src/moi_wrapper.jl @@ -39,7 +39,7 @@ function diff_optimizer( # For this reason we add a cache layer, the same cache JuMP adds. caching_opt = if with_cache MOI.Utilities.CachingOptimizer( - MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), + MOI.Utilities.UniversalFallback(MOI.Utilities.Model{with_bridge_type}()), optimizer, ) else diff --git a/test/Project.toml b/test/Project.toml index bd7cd5af4..fdb058f22 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -16,7 +16,6 @@ SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" [compat] HiGHS = "1" From 91865cf1483b410d642325749489a332c1cfd3b0 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Mon, 13 Jan 2025 18:58:52 -0800 Subject: [PATCH 16/26] add suggestion --- Project.toml | 2 +- src/moi_wrapper.jl | 4 +++- src/parameters.jl | 27 ++++++--------------------- 3 files changed, 10 insertions(+), 23 deletions(-) diff --git a/Project.toml b/Project.toml index 32b65d1e1..42b60eeb1 100644 --- a/Project.toml +++ b/Project.toml @@ -22,6 +22,6 @@ IterativeSolvers = "0.9" JuMP = "1" LazyArrays = "0.21, 0.22, 1" MathOptInterface = "1.18" -MathOptSetDistances = "0.2.7" +MathOptSetDistances = "0.2.9" ParametricOptInterface = "0.9.0" julia = "1.6" diff --git a/src/moi_wrapper.jl b/src/moi_wrapper.jl index f47a83f94..2f3038862 100644 --- a/src/moi_wrapper.jl +++ b/src/moi_wrapper.jl @@ -39,7 +39,9 @@ function diff_optimizer( # For this reason we add a cache layer, the same cache JuMP adds. caching_opt = if with_cache MOI.Utilities.CachingOptimizer( - MOI.Utilities.UniversalFallback(MOI.Utilities.Model{with_bridge_type}()), + MOI.Utilities.UniversalFallback( + MOI.Utilities.Model{with_bridge_type}(), + ), optimizer, ) else diff --git a/src/parameters.jl b/src/parameters.jl index 757e98b38..5deebb030 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -76,22 +76,22 @@ function SensitivityData{T}() where {T} ) end +const _SENSITIVITY_DATA = :_sensitivity_data + function _get_sensitivity_data( model::POI.Optimizer{T}, )::SensitivityData{T} where {T} _initialize_sensitivity_data!(model) - return model.ext[:_sensitivity_data]::SensitivityData{T} + return model.ext[_SENSITIVITY_DATA]::SensitivityData{T} end function _initialize_sensitivity_data!(model::POI.Optimizer{T}) where {T} - if !haskey(model.ext, :_sensitivity_data) - model.ext[:_sensitivity_data] = SensitivityData{T}() + if !haskey(model.ext, _SENSITIVITY_DATA) + model.ext[_SENSITIVITY_DATA] = SensitivityData{T}() end return end -const DoubleDictInner = MOI.Utilities.DoubleDicts.DoubleDictInner - # forward mode function _constraint_set_forward!( @@ -280,7 +280,7 @@ end function empty_input_sensitivities!(model::POI.Optimizer{T}) where {T} empty_input_sensitivities!(model.optimizer) - model.ext[:_sensitivity_data] = SensitivityData{T}() + model.ext[_SENSITIVITY_DATA] = SensitivityData{T}() return end @@ -304,21 +304,6 @@ function forward_differentiate!(model::POI.Optimizer{T}) where {T} return end -# function MOI.set( -# model::POI.Optimizer, -# ::ForwardParameter, -# variable::MOI.VariableIndex, -# value::Number, -# ) -# if _is_variable(model, variable) -# error("Trying to set a forward parameter sensitivity for a variable") -# end -# parameter = variable -# sensitivity_data = _get_sensitivity_data(model) -# sensitivity_data.parameter_input_forward[parameter] = value -# return -# end - function MOI.set( model::POI.Optimizer, ::ForwardConstraintSet, From 9ec668b03c4b42c32f176d571654c35dee7a0174 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Mon, 13 Jan 2025 20:42:39 -0800 Subject: [PATCH 17/26] use Parameter set --- README.md | 8 +-- src/DiffOpt.jl | 2 +- src/jump_moi_overloads.jl | 31 +++++++++ src/parameters.jl | 19 ++---- test/parameters.jl | 139 ++++++++++++++++++++++++-------------- 5 files changed, 129 insertions(+), 70 deletions(-) diff --git a/README.md b/README.md index b14f13be9..51a8cacaf 100644 --- a/README.md +++ b/README.md @@ -64,7 +64,7 @@ optimize!(model) # differentiate w.r.t. p direction_p = 3.0 -MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) +MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), Parameter(direction_p)) DiffOpt.forward_differentiate!(model) @show MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) == direction_p * 3 / pc_val @@ -82,7 +82,7 @@ optimize!(model) DiffOpt.empty_input_sensitivities!(model) # differentiate w.r.t. pc direction_pc = 10.0 -MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(pc), direction_pc) +MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(pc), Parameter(direction_pc)) DiffOpt.forward_differentiate!(model) @show abs(MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) - -direction_pc * 3 * p_val / pc_val^2) < 1e-5 @@ -93,8 +93,8 @@ DiffOpt.empty_input_sensitivities!(model) direction_x = 10.0 MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, direction_x) DiffOpt.reverse_differentiate!(model) -@show MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) == direction_x * 3 / pc_val -@show abs(MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(pc)) - +@show MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) == MOI.Parameter(direction_x * 3 / pc_val) +@show abs(MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(pc)).value - -direction_x * 3 * p_val / pc_val^2) < 1e-5 ``` diff --git a/src/DiffOpt.jl b/src/DiffOpt.jl index f09d26082..71dd1a005 100644 --- a/src/DiffOpt.jl +++ b/src/DiffOpt.jl @@ -19,8 +19,8 @@ include("utils.jl") include("product_of_sets.jl") include("diff_opt.jl") include("moi_wrapper.jl") -include("jump_moi_overloads.jl") include("parameters.jl") +include("jump_moi_overloads.jl") include("copy_dual.jl") include("bridges.jl") diff --git a/src/jump_moi_overloads.jl b/src/jump_moi_overloads.jl index 6f8c29ae8..32ab6dad8 100644 --- a/src/jump_moi_overloads.jl +++ b/src/jump_moi_overloads.jl @@ -80,6 +80,37 @@ function MOI.get( return _moi_get_result(JuMP.backend(model), attr, JuMP.index(var_ref)) end +# extras to handle model_dirty + +function MOI.get( + model::JuMP.Model, + attr::ReverseConstraintSet, + var_ref::JuMP.ConstraintRef, +) + JuMP.check_belongs_to_model(var_ref, model) + return _moi_get_result(JuMP.backend(model), attr, JuMP.index(var_ref)) +end + +function MOI.set( + model::JuMP.Model, + attr::ForwardConstraintSet, + con_ref::JuMP.ConstraintRef, + set::MOI.AbstractScalarSet, +) + JuMP.check_belongs_to_model(con_ref, model) + return MOI.set(JuMP.backend(model), attr, JuMP.index(con_ref), set) +end + +function MOI.set( + model::JuMP.Model, + attr::ForwardConstraintSet, + con_ref::JuMP.ConstraintRef, + set::JuMP.AbstractScalarSet, +) + JuMP.check_belongs_to_model(con_ref, model) + return MOI.set(JuMP.backend(model), attr, JuMP.index(con_ref), JuMP.moi_set(set)) +end + """ abstract type AbstractLazyScalarFunction <: MOI.AbstractScalarFunction end diff --git a/src/parameters.jl b/src/parameters.jl index 5deebb030..65798949c 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -308,14 +308,14 @@ function MOI.set( model::POI.Optimizer, ::ForwardConstraintSet, ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Parameter{T}}, - value::Number, + set::MOI.Parameter, ) where {T} variable = MOI.VariableIndex(ci.value) if _is_variable(model, variable) error("Trying to set a forward parameter sensitivity for a variable") end sensitivity_data = _get_sensitivity_data(model) - sensitivity_data.parameter_input_forward[variable] = value + sensitivity_data.parameter_input_forward[variable] = set.value return end @@ -573,16 +573,7 @@ function MOI.get( error("Trying to get a backward parameter sensitivity for a variable") end sensitivity_data = _get_sensitivity_data(model) - return get(sensitivity_data.parameter_output_backward, variable, 0.0) -end - -# extras to handle model_dirty - -function MOI.get( - model::JuMP.Model, - attr::ReverseConstraintSet, - var_ref::JuMP.ConstraintRef, -) - JuMP.check_belongs_to_model(var_ref, model) - return _moi_get_result(JuMP.backend(model), attr, JuMP.index(var_ref)) + return MOI.Parameter{T}( + get(sensitivity_data.parameter_output_backward, variable, 0.0), + ) end diff --git a/test/parameters.jl b/test/parameters.jl index 9c901d579..8d8794303 100644 --- a/test/parameters.jl +++ b/test/parameters.jl @@ -14,6 +14,9 @@ import MathOptInterface as MOI import HiGHS import SCS +Base.isapprox(x::MOI.Parameter, y::MOI.Parameter; atol = 1e-10) = + isapprox(x.value, y.value, atol = atol) + function runtests() for name in names(@__MODULE__; all = true) if startswith("$name", "test_") @@ -42,11 +45,21 @@ function test_diff_rhs() # the function is # x(p) = 3p, hence x'(p) = 3 # differentiate w.r.t. p - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), 1) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + Parameter(1.0), + ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 3 # again with different "direction" - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), 2) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + Parameter(2.0), + ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 6 # @@ -54,11 +67,21 @@ function test_diff_rhs() optimize!(model) @test value(x) ≈ 6 # differentiate w.r.t. p - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), 1) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + Parameter(1.0), + ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 3 # again with different "direction" - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), 2) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + Parameter(2.0), + ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 6 # @@ -66,11 +89,13 @@ function test_diff_rhs() # MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, 1) DiffOpt.reverse_differentiate!(model) - @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) ≈ 3 + @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) ≈ + MOI.Parameter(3.0) # again with different "direction" MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, 2) DiffOpt.reverse_differentiate!(model) - @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) ≈ 6 + @test MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)) ≈ + MOI.Parameter(6.0) return end @@ -99,12 +124,12 @@ function test_diff_vector_rhs() set_parameter_value(p, p_val) optimize!(model) @test isapprox(value(x), 3 * p_val, atol = 1e-3) - for direction in 0:3 + for direction in 0.0:3.0 MOI.set( model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), - direction, + Parameter(direction), ) DiffOpt.forward_differentiate!(model) @test isapprox( @@ -117,7 +142,7 @@ function test_diff_vector_rhs() DiffOpt.reverse_differentiate!(model) @test isapprox( MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)), - direction * 3, + MOI.Parameter(direction * 3), atol = 1e-3, ) end @@ -145,12 +170,12 @@ function test_affine_changes() # the function is # x(p, pc) = 3p / pc, hence dx/dp = 3 / pc, dx/dpc = -3p / pc^2 # differentiate w.r.t. p - for direction_p in 1:1#2 + for direction_p in 1.0:2.0 MOI.set( model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), - direction_p, + Parameter(direction_p), ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ @@ -162,12 +187,12 @@ function test_affine_changes() optimize!(model) @test value(x) ≈ 3 * p_val / pc_val # differentiate w.r.t. p - for direction_p in 1:2 + for direction_p in 1.0:2.0 MOI.set( model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), - direction_p, + Parameter(direction_p), ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ @@ -176,13 +201,18 @@ function test_affine_changes() # differentiate w.r.t. pc # stop differentiating with respect to p direction_p = 0.0 - MOI.set(model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), direction_p) - for direction_pc in 1:2 + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + Parameter(direction_p), + ) + for direction_pc in 1.0:2.0 MOI.set( model, DiffOpt.ForwardConstraintSet(), ParameterRef(pc), - direction_pc, + Parameter(direction_pc), ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ @@ -193,12 +223,12 @@ function test_affine_changes() set_parameter_value(pc, pc_val) optimize!(model) @test value(x) ≈ 3 * p_val / pc_val - for direction_pc in 1:1#2 + for direction_pc in 1.0:2.0 MOI.set( model, DiffOpt.ForwardConstraintSet(), ParameterRef(pc), - direction_pc, + Parameter(direction_pc), ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ @@ -210,13 +240,13 @@ function test_affine_changes() model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), - direction_p, + Parameter(direction_p), ) MOI.set( model, DiffOpt.ForwardConstraintSet(), ParameterRef(pc), - direction_pc, + Parameter(direction_pc), ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ @@ -247,18 +277,18 @@ function test_affine_changes_compact() set_parameter_value(pc, pc_val) optimize!(model) @test value(x) ≈ 3 * p_val / pc_val - for direction_pc in 0:2, direction_p in 0:2 + for direction_pc in 0.0:2.0, direction_p in 0.0:2.0 MOI.set( model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), - direction_p, + Parameter(direction_p), ) MOI.set( model, DiffOpt.ForwardConstraintSet(), ParameterRef(pc), - direction_pc, + Parameter(direction_pc), ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ @@ -272,12 +302,12 @@ function test_affine_changes_compact() model, DiffOpt.ReverseConstraintSet(), ParameterRef(p), - ) ≈ direction_x * 3 / pc_val + ) ≈ MOI.Parameter(direction_x * 3 / pc_val) @test MOI.get( model, DiffOpt.ReverseConstraintSet(), ParameterRef(pc), - ) ≈ -direction_x * 3 * p_val / pc_val^2 + ) ≈ MOI.Parameter(-direction_x * 3 * p_val / pc_val^2) end end return @@ -324,36 +354,41 @@ function test_quadratic_rhs_changes() optimize!(model) @test value(x) ≈ (1 + 3 * p_val * q_val + 5 * r_val^2 + 7 * s_val) / (11 * t_val) - for dir_p in 0:2, dir_q in 0:2, dir_r in 0:2, dir_s in 0:2, dir_t in 0:2 + for dir_p in 0.0:2.0, + dir_q in 0.0:2.0, + dir_r in 0.0:2.0, + dir_s in 0.0:2.0, + dir_t in 0.0:2.0 + MOI.set( model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), - dir_p, + Parameter(dir_p), ) MOI.set( model, DiffOpt.ForwardConstraintSet(), ParameterRef(q), - dir_q, + Parameter(dir_q), ) MOI.set( model, DiffOpt.ForwardConstraintSet(), ParameterRef(r), - dir_r, + Parameter(dir_r), ) MOI.set( model, DiffOpt.ForwardConstraintSet(), ParameterRef(s), - dir_s, + Parameter(dir_s), ) MOI.set( model, DiffOpt.ForwardConstraintSet(), ParameterRef(t), - dir_t, + Parameter(dir_t), ) DiffOpt.forward_differentiate!(model) @test isapprox( @@ -374,29 +409,31 @@ function test_quadratic_rhs_changes() DiffOpt.reverse_differentiate!(model) @test isapprox( MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(p)), - dir_x * 3 * q_val / (11 * t_val), + MOI.Parameter(dir_x * 3 * q_val / (11 * t_val)), atol = 1e-10, ) @test isapprox( MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(q)), - dir_x * 3 * p_val / (11 * t_val), + MOI.Parameter(dir_x * 3 * p_val / (11 * t_val)), atol = 1e-10, ) @test isapprox( MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(r)), - dir_x * 10 * r_val / (11 * t_val), + MOI.Parameter(dir_x * 10 * r_val / (11 * t_val)), atol = 1e-10, ) @test isapprox( MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(s)), - dir_x * 7 / (11 * t_val), + MOI.Parameter(dir_x * 7 / (11 * t_val)), atol = 1e-10, ) @test isapprox( MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(t)), - dir_x * ( - -(1 + 3 * p_val * q_val + 5 * r_val^2 + 7 * s_val) / - (11 * t_val^2) + MOI.Parameter( + dir_x * ( + -(1 + 3 * p_val * q_val + 5 * r_val^2 + 7 * s_val) / + (11 * t_val^2) + ), ), atol = 1e-10, ) @@ -427,18 +464,18 @@ function test_affine_changes_compact_max() set_parameter_value(pc, pc_val) optimize!(model) @test value(x) ≈ 3 * p_val / pc_val - for direction_pc in 0:2, direction_p in 0:2 + for direction_pc in 0.0:2.0, direction_p in 0.0:2.0 MOI.set( model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), - direction_p, + Parameter(direction_p), ) MOI.set( model, DiffOpt.ForwardConstraintSet(), ParameterRef(pc), - direction_pc, + Parameter(direction_pc), ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ @@ -467,12 +504,12 @@ function test_diff_affine_objective() set_parameter_value(p, p_val) optimize!(model) @test value(x) ≈ 3 - for direction_p in 0:2 + for direction_p in 0.0:2.0 MOI.set( model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), - direction_p, + Parameter(direction_p), ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 0.0 @@ -483,7 +520,7 @@ function test_diff_affine_objective() model, DiffOpt.ReverseConstraintSet(), ParameterRef(p), - ) ≈ 0.0 + ) ≈ MOI.Parameter(0.0) end end return @@ -507,12 +544,12 @@ function test_diff_quadratic_objective() set_parameter_value(p, p_val) optimize!(model) @test value(x) ≈ 3 - for direction_p in 0:2 + for direction_p in 0.0:2.0 MOI.set( model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), - direction_p, + Parameter(direction_p), ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 0.0 @@ -523,7 +560,7 @@ function test_diff_quadratic_objective() model, DiffOpt.ReverseConstraintSet(), ParameterRef(p), - ) ≈ 0.0 + ) ≈ MOI.Parameter(0.0) end end return @@ -548,12 +585,12 @@ function test_quadratic_objective_qp() set_parameter_value(p, p_val) optimize!(model) @test value(x) ≈ -3p_val / 2 atol = 1e-4 - for direction_p in 0:2 + for direction_p in 0.0:2.0 MOI.set( model, DiffOpt.ForwardConstraintSet(), ParameterRef(p), - direction_p, + Parameter(direction_p), ) DiffOpt.forward_differentiate!(model) @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ @@ -565,7 +602,7 @@ function test_quadratic_objective_qp() model, DiffOpt.ReverseConstraintSet(), ParameterRef(p), - ) ≈ direction_p * (-3 / 2) atol = 1e-4 + ) ≈ MOI.Parameter(direction_p * (-3 / 2)) atol = 1e-4 end end return @@ -590,7 +627,7 @@ function test_diff_errors() model, DiffOpt.ForwardConstraintSet(), ParameterRef(x), - 1, + Parameter(1.0), ) @test_throws ErrorException MOI.set( model, From 4cddbe068e9487dbaea56cee0f7ed2a95c026833 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Mon, 13 Jan 2025 20:44:30 -0800 Subject: [PATCH 18/26] todo was fixed --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 51a8cacaf..2af524462 100644 --- a/README.md +++ b/README.md @@ -49,7 +49,7 @@ pc_val = 2.0 @variable(model, x) @variable(model, p in Parameter(p_val)) @variable(model, pc in Parameter(pc_val)) -@constraint(model, cons, pc * x >= 3 * p) #??? InvalidConstraintRef TODO +@constraint(model, cons, pc * x >= 3 * p) @objective(model, Min, 2x) optimize!(model) @show value(x) == 3 * p_val / pc_val From a0ed45f32344168649ed83307ecec1511ef3fb71 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Mon, 13 Jan 2025 20:49:03 -0800 Subject: [PATCH 19/26] format --- src/jump_moi_overloads.jl | 7 ++++++- test/parameters.jl | 5 +++-- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/src/jump_moi_overloads.jl b/src/jump_moi_overloads.jl index 32ab6dad8..13f9ee582 100644 --- a/src/jump_moi_overloads.jl +++ b/src/jump_moi_overloads.jl @@ -108,7 +108,12 @@ function MOI.set( set::JuMP.AbstractScalarSet, ) JuMP.check_belongs_to_model(con_ref, model) - return MOI.set(JuMP.backend(model), attr, JuMP.index(con_ref), JuMP.moi_set(set)) + return MOI.set( + JuMP.backend(model), + attr, + JuMP.index(con_ref), + JuMP.moi_set(set), + ) end """ diff --git a/test/parameters.jl b/test/parameters.jl index 8d8794303..294fa1e8a 100644 --- a/test/parameters.jl +++ b/test/parameters.jl @@ -14,8 +14,9 @@ import MathOptInterface as MOI import HiGHS import SCS -Base.isapprox(x::MOI.Parameter, y::MOI.Parameter; atol = 1e-10) = - isapprox(x.value, y.value, atol = atol) +function Base.isapprox(x::MOI.Parameter, y::MOI.Parameter; atol = 1e-10) + return isapprox(x.value, y.value; atol = atol) +end function runtests() for name in names(@__MODULE__; all = true) From ede336456f49f22464af6113c7ba727b2467f4b9 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Mon, 13 Jan 2025 21:20:40 -0800 Subject: [PATCH 20/26] update docs for newer Flux --- docs/src/examples/custom-relu.jl | 23 ++++++++--------------- docs/src/examples/polyhedral_project.jl | 23 ++++++++--------------- 2 files changed, 16 insertions(+), 30 deletions(-) diff --git a/docs/src/examples/custom-relu.jl b/docs/src/examples/custom-relu.jl index 13e795e5c..e6b440f28 100644 --- a/docs/src/examples/custom-relu.jl +++ b/docs/src/examples/custom-relu.jl @@ -32,7 +32,7 @@ function matrix_relu( @variable(model, x[1:layer_size, 1:batch_size] >= 0) @objective(model, Min, x[:]'x[:] - 2y[:]'x[:]) optimize!(model) - return value.(x) + return Float32.(value.(x)) end # Define the reverse differentiation rule, for the function we defined above. @@ -76,13 +76,13 @@ m = Flux.Chain( N = 1000 # batch size ## Preprocessing train data -imgs = MLDatasets.MNIST.traintensor(1:N) -labels = MLDatasets.MNIST.trainlabels(1:N) +imgs = MLDatasets.MNIST(split=:train).features[:,:,1:N] +labels = MLDatasets.MNIST(split=:train).targets[1:N] train_X = float.(reshape(imgs, size(imgs, 1) * size(imgs, 2), N)) # stack images train_Y = Flux.onehotbatch(labels, 0:9); ## Preprocessing test data -test_imgs = MLDatasets.MNIST.testtensor(1:N) -test_labels = MLDatasets.MNIST.testlabels(1:N) +test_imgs = MLDatasets.MNIST(split=:test).features[:,:,1:N] +test_labels = MLDatasets.MNIST(split=:test).targets[1:N]; test_X = float.(reshape(test_imgs, size(test_imgs, 1) * size(test_imgs, 2), N)) test_Y = Flux.onehotbatch(test_labels, 0:9); @@ -97,19 +97,12 @@ dataset = repeated((train_X, train_Y), epochs); # ## Network training # training loss function, Flux optimizer -custom_loss(x, y) = Flux.crossentropy(m(x), y) -opt = Flux.Adam() -evalcb = () -> @show(custom_loss(train_X, train_Y)) +custom_loss(m, x, y) = Flux.crossentropy(m(x), y) +opt = Flux.setup(Flux.Adam(), m) # Train to optimize network parameters -@time Flux.train!( - custom_loss, - Flux.params(m), - dataset, - opt, - cb = Flux.throttle(evalcb, 5), -); +@time Flux.train!(custom_loss, m, dataset, opt); # Although our custom implementation takes time, it is able to reach similar # accuracy as the usual ReLU function implementation. diff --git a/docs/src/examples/polyhedral_project.jl b/docs/src/examples/polyhedral_project.jl index e1013b94d..dc9acdf20 100644 --- a/docs/src/examples/polyhedral_project.jl +++ b/docs/src/examples/polyhedral_project.jl @@ -54,7 +54,7 @@ function (polytope::Polytope{N})( ) @objective(model, Min, dot(x - y, x - y)) optimize!(model) - return JuMP.value.(x) + return Float32.(JuMP.value.(x)) end # The `@functor` macro from Flux implements auxiliary functions for collecting the parameters of @@ -122,13 +122,13 @@ m = Flux.Chain( M = 500 # batch size ## Preprocessing train data -imgs = MLDatasets.MNIST.traintensor(1:M) -labels = MLDatasets.MNIST.trainlabels(1:M); +imgs = MLDatasets.MNIST(split=:train).features[:,:,1:M] +labels = MLDatasets.MNIST(split=:train).targets[1:M] train_X = float.(reshape(imgs, size(imgs, 1) * size(imgs, 2), M)) # stack images train_Y = Flux.onehotbatch(labels, 0:9); ## Preprocessing test data -test_imgs = MLDatasets.MNIST.testtensor(1:M) -test_labels = MLDatasets.MNIST.testlabels(1:M) +test_imgs = MLDatasets.MNIST(split=:test).features[:,:,1:M] +test_labels = MLDatasets.MNIST(split=:test).targets[1:M] test_X = float.(reshape(test_imgs, size(test_imgs, 1) * size(test_imgs, 2), M)) test_Y = Flux.onehotbatch(test_labels, 0:9); @@ -142,19 +142,12 @@ dataset = repeated((train_X, train_Y), epochs); # ## Network training # training loss function, Flux optimizer -custom_loss(x, y) = Flux.crossentropy(m(x), y) -opt = Flux.ADAM() -evalcb = () -> @show(custom_loss(train_X, train_Y)) +custom_loss(m, x, y) = Flux.crossentropy(m(x), y) +opt = Flux.setup(Flux.Adam(), m) # Train to optimize network parameters -@time Flux.train!( - custom_loss, - Flux.params(m), - dataset, - opt, - cb = Flux.throttle(evalcb, 5), -); +@time Flux.train!(custom_loss, m, dataset, opt); # Although our custom implementation takes time, it is able to reach similar # accuracy as the usual ReLU function implementation. From fe2446c271127c279f8bdb7bb83b86c79d602b89 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Mon, 13 Jan 2025 21:22:47 -0800 Subject: [PATCH 21/26] format --- docs/src/examples/custom-relu.jl | 8 ++++---- docs/src/examples/polyhedral_project.jl | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/docs/src/examples/custom-relu.jl b/docs/src/examples/custom-relu.jl index e6b440f28..ee3477c30 100644 --- a/docs/src/examples/custom-relu.jl +++ b/docs/src/examples/custom-relu.jl @@ -76,13 +76,13 @@ m = Flux.Chain( N = 1000 # batch size ## Preprocessing train data -imgs = MLDatasets.MNIST(split=:train).features[:,:,1:N] -labels = MLDatasets.MNIST(split=:train).targets[1:N] +imgs = MLDatasets.MNIST(; split = :train).features[:, :, 1:N] +labels = MLDatasets.MNIST(; split = :train).targets[1:N] train_X = float.(reshape(imgs, size(imgs, 1) * size(imgs, 2), N)) # stack images train_Y = Flux.onehotbatch(labels, 0:9); ## Preprocessing test data -test_imgs = MLDatasets.MNIST(split=:test).features[:,:,1:N] -test_labels = MLDatasets.MNIST(split=:test).targets[1:N]; +test_imgs = MLDatasets.MNIST(; split = :test).features[:, :, 1:N] +test_labels = MLDatasets.MNIST(; split = :test).targets[1:N]; test_X = float.(reshape(test_imgs, size(test_imgs, 1) * size(test_imgs, 2), N)) test_Y = Flux.onehotbatch(test_labels, 0:9); diff --git a/docs/src/examples/polyhedral_project.jl b/docs/src/examples/polyhedral_project.jl index dc9acdf20..b0f719db8 100644 --- a/docs/src/examples/polyhedral_project.jl +++ b/docs/src/examples/polyhedral_project.jl @@ -122,13 +122,13 @@ m = Flux.Chain( M = 500 # batch size ## Preprocessing train data -imgs = MLDatasets.MNIST(split=:train).features[:,:,1:M] -labels = MLDatasets.MNIST(split=:train).targets[1:M] +imgs = MLDatasets.MNIST(; split = :train).features[:, :, 1:M] +labels = MLDatasets.MNIST(; split = :train).targets[1:M] train_X = float.(reshape(imgs, size(imgs, 1) * size(imgs, 2), M)) # stack images train_Y = Flux.onehotbatch(labels, 0:9); ## Preprocessing test data -test_imgs = MLDatasets.MNIST(split=:test).features[:,:,1:M] -test_labels = MLDatasets.MNIST(split=:test).targets[1:M] +test_imgs = MLDatasets.MNIST(; split = :test).features[:, :, 1:M] +test_labels = MLDatasets.MNIST(; split = :test).targets[1:M] test_X = float.(reshape(test_imgs, size(test_imgs, 1) * size(test_imgs, 2), M)) test_Y = Flux.onehotbatch(test_labels, 0:9); From bddb5ef3122c17bee98346754fdb92583ad0c46a Mon Sep 17 00:00:00 2001 From: joaquimg Date: Mon, 13 Jan 2025 21:32:13 -0800 Subject: [PATCH 22/26] kwargs --- src/moi_wrapper.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/moi_wrapper.jl b/src/moi_wrapper.jl index 2f3038862..6493f6831 100644 --- a/src/moi_wrapper.jl +++ b/src/moi_wrapper.jl @@ -28,10 +28,7 @@ function diff_optimizer( with_bridge_type = Float64, with_cache::Bool = true, ) - optimizer = MOI.instantiate( - optimizer_constructor; - with_bridge_type = with_bridge_type, - ) + optimizer = MOI.instantiate(optimizer_constructor; with_bridge_type) # When we do `MOI.copy_to(diff, optimizer)` we need to efficiently `MOI.get` # the model information from `optimizer`. However, 1) `optimizer` may not # implement some getters or it may be inefficient and 2) the getters may be From a98f54897e428e22d142b74e31357f95d207f717 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Mon, 13 Jan 2025 21:55:19 -0800 Subject: [PATCH 23/26] remove diff model --- docs/src/manual.md | 7 ------- 1 file changed, 7 deletions(-) diff --git a/docs/src/manual.md b/docs/src/manual.md index 52edd6112..614c4be8b 100644 --- a/docs/src/manual.md +++ b/docs/src/manual.md @@ -60,13 +60,6 @@ You can create a differentiable optimizer over an existing MOI solver by using t diff_optimizer ``` -## Creating a differentiable JuMP model - -You initialize a differentiable JuMP model by using the `diff_model` utility. -```@docs -diff_model -``` - ## Projections on cone sets DiffOpt requires taking projections and finding projection gradients of vectors while computing the jacobians. For this purpose, we use [MathOptSetDistances.jl](https://github.com/matbesancon/MathOptSetDistances.jl), which is a dedicated package for computing set distances, projections and projection gradients. From e713f25dfe59e14f3f668a9a1e63c4987bc993e4 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Fri, 31 Jan 2025 12:50:38 -0800 Subject: [PATCH 24/26] suggestions --- src/jump_moi_overloads.jl | 20 +++++++++++--------- src/moi_wrapper.jl | 28 +++++++++++++++++----------- 2 files changed, 28 insertions(+), 20 deletions(-) diff --git a/src/jump_moi_overloads.jl b/src/jump_moi_overloads.jl index 13f9ee582..851cac963 100644 --- a/src/jump_moi_overloads.jl +++ b/src/jump_moi_overloads.jl @@ -3,6 +3,15 @@ # 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. +# FIXME +# Some function in this file are overloads to skip JuMP dirty state. +# Workaround for https://github.com/jump-dev/JuMP.jl/issues/2797 +# This workaround is necessary because once some attributes are set the JuMP +# model changes to a dirty state, then getting some attributes is blocked. +# However, getting and setting forward and backward sensitivities is +# done after the model is optimized, so we add function to bypass the +# dirty state. + function MOI.set( model::JuMP.Model, attr::ForwardObjectiveFunction, @@ -54,7 +63,7 @@ function MOI.get( return JuMP.jump_function(model, moi_func) end -# FIXME Workaround for https://github.com/jump-dev/JuMP.jl/issues/2797 +# see FIXME comment in the top of the file function _moi_get_result(model::MOI.ModelLike, args...) if MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED throw(OptimizeNotCalled()) @@ -80,8 +89,6 @@ function MOI.get( return _moi_get_result(JuMP.backend(model), attr, JuMP.index(var_ref)) end -# extras to handle model_dirty - function MOI.get( model::JuMP.Model, attr::ReverseConstraintSet, @@ -108,12 +115,7 @@ function MOI.set( set::JuMP.AbstractScalarSet, ) JuMP.check_belongs_to_model(con_ref, model) - return MOI.set( - JuMP.backend(model), - attr, - JuMP.index(con_ref), - JuMP.moi_set(set), - ) + return MOI.set(model, attr, con_ref, JuMP.moi_set(set)) end """ diff --git a/src/moi_wrapper.jl b/src/moi_wrapper.jl index 6493f6831..c75d9e013 100644 --- a/src/moi_wrapper.jl +++ b/src/moi_wrapper.jl @@ -17,13 +17,13 @@ One define a differentiable model by using any solver of choice. Example: julia> import DiffOpt, HiGHS julia> model = DiffOpt.diff_optimizer(HiGHS.Optimizer) +julia> set_attribute(model, DiffOpt.ModelConstructor, DiffOpt.QuadraticProgram.Model) # optional selection of diff method julia> x = model.add_variable(model) julia> model.add_constraint(model, ...) ``` """ function diff_optimizer( optimizer_constructor; - method = nothing, with_parametric_opt_interface::Bool = false, with_bridge_type = Float64, with_cache::Bool = true, @@ -45,9 +45,9 @@ function diff_optimizer( optimizer end if with_parametric_opt_interface - return POI.Optimizer(Optimizer(caching_opt; method = method)) + return POI.Optimizer(Optimizer(caching_opt)) else - return Optimizer(caching_opt; method = method) + return Optimizer(caching_opt) end end @@ -64,16 +64,10 @@ mutable struct Optimizer{OT<:MOI.ModelLike} <: MOI.AbstractOptimizer # sensitivity input cache using MOI like sparse format input_cache::InputCache - function Optimizer( - optimizer::OT; - method = nothing, - ) where {OT<:MOI.ModelLike} + function Optimizer(optimizer::OT) where {OT<:MOI.ModelLike} output = new{OT}(optimizer, Any[], nothing, nothing, nothing, InputCache()) add_all_model_constructors(output) - if method !== nothing - output.model_constructor = method - end return output end end @@ -497,6 +491,14 @@ end Determines which subtype of [`DiffOpt.AbstractModel`](@ref) to use for differentiation. When set to `nothing`, the first one out of `model.model_constructors` that support the problem is used. + +Examples: + +```julia +julia> MOI.set(model, DiffOpt.ModelConstructor(), DiffOpt.QuadraticProgram.Model) + +julia> MOI.set(model, DiffOpt.ModelConstructor(), DiffOpt.ConicProgram.Model) +``` """ struct ModelConstructor <: MOI.AbstractOptimizerAttribute end @@ -616,7 +618,11 @@ function _diff(model::Optimizer) end if isnothing(model.diff) error( - "No differentiation model supports the problem. If you believe it should be supported, say by `DiffOpt.QuadraticProgram.Model`, use `MOI.set(model, DiffOpt.ModelConstructor, DiffOpt.QuadraticProgram.Model)` and try again to see an error indicating why it is not supported.", + "No differentiation model supports the problem. If you " * + "believe it should be supported, say by " * + "`DiffOpt.QuadraticProgram.Model`, use " * + "`MOI.set(model, DiffOpt.ModelConstructor, DiffOpt.QuadraticProgram.Model)`" * + "and try again to see an error indicating why it is not supported." ) end else From e8099a1fbdf278e9c252bc207c446110da2a2d66 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Fri, 31 Jan 2025 13:15:54 -0800 Subject: [PATCH 25/26] format --- src/moi_wrapper.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/moi_wrapper.jl b/src/moi_wrapper.jl index c75d9e013..e3e00005b 100644 --- a/src/moi_wrapper.jl +++ b/src/moi_wrapper.jl @@ -622,7 +622,7 @@ function _diff(model::Optimizer) "believe it should be supported, say by " * "`DiffOpt.QuadraticProgram.Model`, use " * "`MOI.set(model, DiffOpt.ModelConstructor, DiffOpt.QuadraticProgram.Model)`" * - "and try again to see an error indicating why it is not supported." + "and try again to see an error indicating why it is not supported.", ) end else From 35cd234dc7716ae333451b41a38e1747155157b4 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Fri, 31 Jan 2025 14:24:29 -0800 Subject: [PATCH 26/26] fix examples --- docs/src/examples/custom-relu.jl | 6 +++--- docs/src/examples/polyhedral_project.jl | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/src/examples/custom-relu.jl b/docs/src/examples/custom-relu.jl index ee3477c30..ce68ec5a5 100644 --- a/docs/src/examples/custom-relu.jl +++ b/docs/src/examples/custom-relu.jl @@ -42,9 +42,9 @@ function ChainRulesCore.rrule(::typeof(matrix_relu), y::Matrix{T}) where {T} function pullback_matrix_relu(dl_dx) ## some value from the backpropagation (e.g., loss) is denoted by `l` ## so `dl_dy` is the derivative of `l` wrt `y` - x = model[:x] # load decision variable `x` into scope - dl_dy = zeros(T, size(dl_dx)) - dl_dq = zeros(T, size(dl_dx)) + x = model[:x]::Matrix{JuMP.VariableRef} # load decision variable `x` into scope + dl_dy = zeros(T, size(x)) + dl_dq = zeros(T, size(x)) ## set sensitivities MOI.set.(model, DiffOpt.ReverseVariablePrimal(), x[:], dl_dx[:]) ## compute grad diff --git a/docs/src/examples/polyhedral_project.jl b/docs/src/examples/polyhedral_project.jl index b0f719db8..ae2421cd4 100644 --- a/docs/src/examples/polyhedral_project.jl +++ b/docs/src/examples/polyhedral_project.jl @@ -75,12 +75,12 @@ function ChainRulesCore.rrule( model = direct_model(DiffOpt.diff_optimizer(Ipopt.Optimizer)) xv = polytope(y; model = model) function pullback_matrix_projection(dl_dx) - layer_size, batch_size = size(dl_dx) dl_dx = ChainRulesCore.unthunk(dl_dx) ## `dl_dy` is the derivative of `l` wrt `y` - x = model[:x] + x = model[:x]::Matrix{JuMP.VariableRef} + layer_size, batch_size = size(x) ## grad wrt input parameters - dl_dy = zeros(size(dl_dx)) + dl_dy = zeros(size(x)) ## grad wrt layer parameters dl_dw = zero.(polytope.w) dl_db = zero(polytope.b)