From 27a067945c44142471a4a0a226f106ab229e281a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sun, 12 Jan 2025 12:36:39 +0100 Subject: [PATCH 1/4] Add error for missing starting value --- src/ConicProgram/ConicProgram.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/ConicProgram/ConicProgram.jl b/src/ConicProgram/ConicProgram.jl index 69df7e09..5b5a73bd 100644 --- a/src/ConicProgram/ConicProgram.jl +++ b/src/ConicProgram/ConicProgram.jl @@ -183,6 +183,14 @@ function _gradient_cache(model::Model) ) b = model.model.constraints.constants + if any(isnan, model.y) || length(model.y) < length(b) + error("Some constraints are missing a value for the `ConstraintDualStart` attribute.") + end + + if any(isnan, model.s) || length(model.s) < length(b) + error("Some constraints are missing a value for the `ConstraintPrimalStart` attribute.") + end + if MOI.get(model, MOI.ObjectiveSense()) == MOI.FEASIBILITY_SENSE c = SparseArrays.spzeros(size(A, 2)) else From 6dd82d6155f92c94cf3f9908baed79fca5eee1dc Mon Sep 17 00:00:00 2001 From: Joaquim Date: Fri, 31 Jan 2025 15:58:49 -0800 Subject: [PATCH 2/4] Integrate with POI to improve UX (#262) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * [WIP] Integrate with POI to improve UX * add missing import * temp change to proj toml * format * simplify method setting to sue model constructor * add possible fix to scalarize bridge error * add pkg to project * format * improvements * remove jump wrapper * clean tests * fix readme * use intermediary API * format * Apply suggestions from code review Co-authored-by: Benoît Legat * add suggestion * use Parameter set * todo was fixed * format * update docs for newer Flux * format * kwargs * remove diff model * suggestions * format * fix examples --------- Co-authored-by: Benoît Legat --- Project.toml | 4 +- README.md | 71 ++- docs/src/examples/custom-relu.jl | 29 +- docs/src/examples/polyhedral_project.jl | 29 +- docs/src/manual.md | 21 +- src/DiffOpt.jl | 5 + src/bridges.jl | 19 + src/diff_opt.jl | 11 + src/jump_moi_overloads.jl | 55 +- src/moi_wrapper.jl | 52 +- src/parameters.jl | 579 ++++++++++++++++++++ src/utils.jl | 17 + test/parameters.jl | 676 ++++++++++++++++++++++++ test/utils.jl | 14 +- 14 files changed, 1513 insertions(+), 69 deletions(-) create mode 100644 src/parameters.jl create mode 100644 test/parameters.jl diff --git a/Project.toml b/Project.toml index a15809f0..42b60eeb 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] @@ -21,5 +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/README.md b/README.md index a90b9a9b..2af52446 100644 --- a/README.md +++ b/README.md @@ -31,7 +31,76 @@ 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, + ), +) +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) +@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 +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 + +# 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 +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 + +# always a good practice to clear previously set sensitivities +DiffOpt.empty_input_sensitivities!(model) +# Now, reverse model AD +direction_x = 10.0 +MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, direction_x) +DiffOpt.reverse_differentiate!(model) +@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 +``` + +### Low level DiffOpt-JuMP API: + +A brief example: ```julia using JuMP, DiffOpt, HiGHS diff --git a/docs/src/examples/custom-relu.jl b/docs/src/examples/custom-relu.jl index 13e795e5..ce68ec5a 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. @@ -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 @@ -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 e1013b94..ae2421cd 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 @@ -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) @@ -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. diff --git a/docs/src/manual.md b/docs/src/manual.md index 57f00721..614c4be8 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, the package supports following `Function-in-Set` constraints: | MOI Function | MOI Set | |:-------|:---------------| @@ -50,19 +50,16 @@ 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 - -The DiffOpt `Optimizer` behaves similarly to other MOI Optimizers -and implements the `MOI.AbstractOptimizer` API. - ## 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. @@ -104,6 +101,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 9ab20101..71dd1a00 100644 --- a/src/DiffOpt.jl +++ b/src/DiffOpt.jl @@ -12,12 +12,14 @@ import LazyArrays import LinearAlgebra import MathOptInterface as MOI import MathOptSetDistances as MOSD +import ParametricOptInterface as POI import SparseArrays include("utils.jl") include("product_of_sets.jl") include("diff_opt.jl") include("moi_wrapper.jl") +include("parameters.jl") include("jump_moi_overloads.jl") include("copy_dual.jl") @@ -40,4 +42,7 @@ end export diff_optimizer +# TODO +# add precompilation statements + end # module diff --git a/src/bridges.jl b/src/bridges.jl index 02299d65..2b0eb46b 100644 --- a/src/bridges.jl +++ b/src/bridges.jl @@ -43,6 +43,25 @@ 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, + value, +) + 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/diff_opt.jl b/src/diff_opt.jl index 22781eec..c91ab9e8 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 9459f840..851cac96 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,6 +89,35 @@ function MOI.get( return _moi_get_result(JuMP.backend(model), attr, JuMP.index(var_ref)) end +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(model, attr, con_ref, JuMP.moi_set(set)) +end + """ abstract type AbstractLazyScalarFunction <: MOI.AbstractScalarFunction end @@ -307,6 +345,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 +360,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 +374,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 3570d94d..e3e00005 100644 --- a/src/moi_wrapper.jl +++ b/src/moi_wrapper.jl @@ -4,7 +4,7 @@ # in the LICENSE.md file or at https://opensource.org/licenses/MIT. """ - diff_optimizer(optimizer_constructor)::Optimizer + 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 @@ -17,23 +17,38 @@ 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)::Optimizer - optimizer = - MOI.instantiate(optimizer_constructor; with_bridge_type = Float64) +function diff_optimizer( + optimizer_constructor; + with_parametric_opt_interface::Bool = false, + with_bridge_type = Float64, + with_cache::Bool = true, +) + 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 # 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{with_bridge_type}(), + ), + optimizer, + ) + else + optimizer + end + if with_parametric_opt_interface + return POI.Optimizer(Optimizer(caching_opt)) + else + return Optimizer(caching_opt) + end end mutable struct Optimizer{OT<:MOI.ModelLike} <: MOI.AbstractOptimizer @@ -476,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 @@ -552,6 +575,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: @@ -590,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 diff --git a/src/parameters.jl b/src/parameters.jl new file mode 100644 index 00000000..65798949 --- /dev/null +++ b/src/parameters.jl @@ -0,0 +1,579 @@ +# 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, _) + return 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, + _, +) + return 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) + return 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, +) + return 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 + +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} +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 + +# 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{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.optimizer) + 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, + ::ForwardConstraintSet, + ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Parameter{T}}, + 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] = set.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} + 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} + return 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(::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 MOI.Parameter{T}( + get(sensitivity_data.parameter_output_backward, variable, 0.0), + ) +end diff --git a/src/utils.jl b/src/utils.jl index dad12f7a..d6764c1f 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 @@ -261,6 +267,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 new file mode 100644 index 00000000..294fa1e8 --- /dev/null +++ b/test/parameters.jl @@ -0,0 +1,676 @@ +# 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 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) + 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), + 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), + Parameter(2.0), + ) + DiffOpt.forward_differentiate!(model) + @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), + 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), + Parameter(2.0), + ) + DiffOpt.forward_differentiate!(model) + @test MOI.get(model, DiffOpt.ForwardVariablePrimal(), x) ≈ 6 + # + # test reverse mode + # + MOI.set(model, DiffOpt.ReverseVariablePrimal(), x, 1) + DiffOpt.reverse_differentiate!(model) + @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)) ≈ + MOI.Parameter(6.0) + 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.0:3.0 + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + Parameter(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)), + MOI.Parameter(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, + ), + ) + 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.0:2.0 + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + Parameter(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.0:2.0 + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + Parameter(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), + Parameter(direction_p), + ) + for direction_pc in 1.0:2.0 + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(pc), + Parameter(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.0:2.0 + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(pc), + Parameter(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), + Parameter(direction_p), + ) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(pc), + Parameter(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.0:2.0, direction_p in 0.0:2.0 + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + Parameter(direction_p), + ) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(pc), + Parameter(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), + ) ≈ MOI.Parameter(direction_x * 3 / pc_val) + @test MOI.get( + model, + DiffOpt.ReverseConstraintSet(), + ParameterRef(pc), + ) ≈ MOI.Parameter(-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.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), + Parameter(dir_p), + ) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(q), + Parameter(dir_q), + ) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(r), + Parameter(dir_r), + ) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(s), + Parameter(dir_s), + ) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(t), + Parameter(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)), + MOI.Parameter(dir_x * 3 * q_val / (11 * t_val)), + atol = 1e-10, + ) + @test isapprox( + MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(q)), + MOI.Parameter(dir_x * 3 * p_val / (11 * t_val)), + atol = 1e-10, + ) + @test isapprox( + MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(r)), + MOI.Parameter(dir_x * 10 * r_val / (11 * t_val)), + atol = 1e-10, + ) + @test isapprox( + MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(s)), + MOI.Parameter(dir_x * 7 / (11 * t_val)), + atol = 1e-10, + ) + @test isapprox( + MOI.get(model, DiffOpt.ReverseConstraintSet(), ParameterRef(t)), + MOI.Parameter( + 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.0:2.0, direction_p in 0.0:2.0 + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + Parameter(direction_p), + ) + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(pc), + Parameter(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.0:2.0 + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + Parameter(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), + ) ≈ MOI.Parameter(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.0:2.0 + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + Parameter(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), + ) ≈ MOI.Parameter(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.0:2.0 + MOI.set( + model, + DiffOpt.ForwardConstraintSet(), + ParameterRef(p), + Parameter(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), + ) ≈ MOI.Parameter(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), + Parameter(1.0), + ) + @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/utils.jl b/test/utils.jl index 10e37eac..49cfac38 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 76fb4b96a22d613a295b71c58ae650b47244000f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sun, 12 Jan 2025 12:36:39 +0100 Subject: [PATCH 3/4] Add error for missing starting value --- src/ConicProgram/ConicProgram.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/ConicProgram/ConicProgram.jl b/src/ConicProgram/ConicProgram.jl index 69df7e09..5b5a73bd 100644 --- a/src/ConicProgram/ConicProgram.jl +++ b/src/ConicProgram/ConicProgram.jl @@ -183,6 +183,14 @@ function _gradient_cache(model::Model) ) b = model.model.constraints.constants + if any(isnan, model.y) || length(model.y) < length(b) + error("Some constraints are missing a value for the `ConstraintDualStart` attribute.") + end + + if any(isnan, model.s) || length(model.s) < length(b) + error("Some constraints are missing a value for the `ConstraintPrimalStart` attribute.") + end + if MOI.get(model, MOI.ObjectiveSense()) == MOI.FEASIBILITY_SENSE c = SparseArrays.spzeros(size(A, 2)) else From 72b50c534d8d4b03a6b3e2139bcdbcfc4a98ba93 Mon Sep 17 00:00:00 2001 From: joaquimg Date: Fri, 31 Jan 2025 17:24:44 -0800 Subject: [PATCH 4/4] format --- src/ConicProgram/ConicProgram.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/ConicProgram/ConicProgram.jl b/src/ConicProgram/ConicProgram.jl index 5b5a73bd..89a79183 100644 --- a/src/ConicProgram/ConicProgram.jl +++ b/src/ConicProgram/ConicProgram.jl @@ -184,11 +184,15 @@ function _gradient_cache(model::Model) b = model.model.constraints.constants if any(isnan, model.y) || length(model.y) < length(b) - error("Some constraints are missing a value for the `ConstraintDualStart` attribute.") + error( + "Some constraints are missing a value for the `ConstraintDualStart` attribute.", + ) end if any(isnan, model.s) || length(model.s) < length(b) - error("Some constraints are missing a value for the `ConstraintPrimalStart` attribute.") + error( + "Some constraints are missing a value for the `ConstraintPrimalStart` attribute.", + ) end if MOI.get(model, MOI.ObjectiveSense()) == MOI.FEASIBILITY_SENSE