Merge branch 'master' into ar/VectorOfVariables

Author: andrewrosemberg

.github/workflows/ci.yml

@@ -17,7 +17,7 @@ jobs:
       fail-fast: false
       matrix:
         include:
-          - version: '1'
+          - version: 'lts'
             os: ubuntu-latest
             arch: x64
           - version: '1'

src/ConicProgram/ConicProgram.jl

@@ -461,15 +461,23 @@ function MOI.get(
     return MOI.get(model.model, attr, ci)
 end
 
+"""
+Method not supported for `DiffOpt.ConicProgram.Model` directly.
+However, a fallback is provided in `DiffOpt`.
+"""
 function MOI.get(::Model, ::DiffOpt.ForwardObjectiveSensitivity)
-    return error(
-        "ForwardObjectiveSensitivity is not implemented for the Conic Optimization backend",
+    return throw(
+        MOI.UnsupportedAttribute(DiffOpt.ForwardObjectiveSensitivity()),
     )
 end
 
+"""
+Method not supported for `DiffOpt.ConicProgram.Model` directly.
+However, a fallback is provided in `DiffOpt`.
+"""
 function MOI.set(::Model, ::DiffOpt.ReverseObjectiveSensitivity, val)
-    return error(
-        "ReverseObjectiveSensitivity is not implemented for the Conic Optimization backend",
+    return throw(
+        MOI.UnsupportedAttribute(DiffOpt.ReverseObjectiveSensitivity()),
     )
 end
 

src/NonLinearProgram/NonLinearProgram.jl

@@ -553,45 +553,46 @@ function DiffOpt.reverse_differentiate!(model::Model; tol = 1e-6)
         # Compute Jacobian
         Δs, df_dp = _compute_sensitivity(model; tol = tol)
         Δp = if !iszero(model.input_cache.dobj)
-            model.input_cache.dobj * df_dp
+            df_dp'model.input_cache.dobj
         else
-            num_primal = length(cache.primal_vars)
-            # Fetch primal sensitivities
-            Δx = zeros(num_primal)
-            for (i, var_idx) in enumerate(cache.primal_vars)
-                if haskey(model.input_cache.dx, var_idx)
-                    Δx[i] = model.input_cache.dx[var_idx]
-                end
+            zeros(length(cache.params))
+        end
+        num_primal = length(cache.primal_vars)
+        # Fetch primal sensitivities
+        Δx = zeros(num_primal)
+        for (i, var_idx) in enumerate(cache.primal_vars)
+            if haskey(model.input_cache.dx, var_idx)
+                Δx[i] = model.input_cache.dx[var_idx]
             end
-            # Fetch dual sensitivities
-            num_constraints = length(cache.cons)
-            num_up = length(cache.has_up)
-            num_low = length(cache.has_low)
-            Δdual = zeros(num_constraints + num_up + num_low)
-            for (i, ci) in enumerate(cache.cons)
-                idx = form.nlp_index_2_constraint[ci]
-                if haskey(model.input_cache.dy, idx)
-                    Δdual[i] = model.input_cache.dy[idx]
-                end
+        end
+        # Fetch dual sensitivities
+        num_constraints = length(cache.cons)
+        num_up = length(cache.has_up)
+        num_low = length(cache.has_low)
+        Δdual = zeros(num_constraints + num_up + num_low)
+        for (i, ci) in enumerate(cache.cons)
+            idx = form.nlp_index_2_constraint[ci]
+            if haskey(model.input_cache.dy, idx)
+                Δdual[i] = model.input_cache.dy[idx]
             end
-            for (i, var_idx) in enumerate(cache.primal_vars[cache.has_low])
-                idx = form.constraint_lower_bounds[var_idx.value]
-                if haskey(model.input_cache.dy, idx)
-                    Δdual[num_constraints+i] = model.input_cache.dy[idx]
-                end
+        end
+        for (i, var_idx) in enumerate(cache.primal_vars[cache.has_low])
+            idx = form.constraint_lower_bounds[var_idx.value]
+            if haskey(model.input_cache.dy, idx)
+                Δdual[num_constraints+i] = model.input_cache.dy[idx]
             end
-            for (i, var_idx) in enumerate(cache.primal_vars[cache.has_up])
-                idx = form.constraint_upper_bounds[var_idx.value]
-                if haskey(model.input_cache.dy, idx)
-                    Δdual[num_constraints+num_low+i] = model.input_cache.dy[idx]
-                end
+        end
+        for (i, var_idx) in enumerate(cache.primal_vars[cache.has_up])
+            idx = form.constraint_upper_bounds[var_idx.value]
+            if haskey(model.input_cache.dy, idx)
+                Δdual[num_constraints+num_low+i] = model.input_cache.dy[idx]
             end
-            # Extract Parameter sensitivities
-            Δw = zeros(size(Δs, 1))
-            Δw[1:num_primal] = Δx
-            Δw[cache.index_duals] = Δdual
-            Δp = Δs' * Δw
         end
+        # Extract Parameter sensitivities
+        Δw = zeros(size(Δs, 1))
+        Δw[1:num_primal] = Δx
+        Δw[cache.index_duals] = Δdual
+        Δp += Δs' * Δw
 
         Δp_dict = Dict{MOI.ConstraintIndex,Float64}(
             form.var2ci[var_idx] => Δp[form.var2param[var_idx].value]

src/NonLinearProgram/nlp_utilities.jl

@@ -489,9 +489,16 @@ function _compute_sensitivity(model::Model; tol = 1e-6)
     # Dual bounds upper
     ∂s[((num_w+num_cons+num_lower+1):end), :] *= -_sense_multiplier
 
-    # dual wrt parameter
+    grad = _compute_gradient(model)
+    # `grad` = [∇ₓf(x,p); ∇ₚf(x,p)] where `x` is the primal vars, `p` is the params,
+    #  and `f(x,p)` is the objective function. we extract the components
+    #  so we can form `∇ₚfᵒ(x,p) = ∇ₓf(x,p) * ∇ₚxᵒ(p) + ∇ₚf(x,p) * ∇ₚpᵒ(p)`
+    #  where `ᵒ` denotes "optimal". note that parameters are fixed, so
+    #  pᵒ(p) = p and ∇ₚpᵒ(p) = 𝐈ₚ.
     primal_idx = [i.value for i in model.cache.primal_vars]
-    df_dx = _compute_gradient(model)[primal_idx]
-    df_dp = df_dx'∂s[1:num_vars, :]
+    params_idx = [i.value for i in model.cache.params]
+    df_dx = grad[primal_idx]  # ∇ₓf(x,p)
+    df_dp_direct = grad[params_idx]  # ∇ₚf(x,p)
+    df_dp = df_dx'∂s[1:num_vars, :] + df_dp_direct'  # ∇ₚfᵒ(x,p) = ∇ₓf(x,p) * ∇ₚxᵒ(p) + ∇ₚf(x,p) * 𝐈ₚ
     return ∂s, df_dp
 end

src/QuadraticProgram/QuadraticProgram.jl

@@ -501,15 +501,23 @@ function MOI.set(model::Model, ::LinearAlgebraSolver, linear_solver)
     return model.linear_solver = linear_solver
 end
 
+"""
+Method not supported for `DiffOpt.QuadraticProgram.Model` directly.
+However, a fallback is provided in `DiffOpt`.
+"""
 function MOI.get(::Model, ::DiffOpt.ForwardObjectiveSensitivity)
-    return error(
-        "ForwardObjectiveSensitivity is not implemented for the Quadratic Optimization backend",
+    return throw(
+        MOI.UnsupportedAttribute(DiffOpt.ForwardObjectiveSensitivity()),
     )
 end
 
+"""
+Method not supported for `DiffOpt.QuadraticProgram.Model` directly.
+However, a fallback is provided in `DiffOpt`.
+"""
 function MOI.set(::Model, ::DiffOpt.ReverseObjectiveSensitivity, val)
-    return error(
-        "ReverseObjectiveSensitivity is not implemented for the Quadratic Optimization backend",
+    return throw(
+        MOI.UnsupportedAttribute(DiffOpt.ReverseObjectiveSensitivity()),
     )
 end
 

src/diff_opt.jl

@@ -31,6 +31,7 @@ Base.@kwdef mutable struct InputCache
         MOIDD.DoubleDict{MOI.VectorAffineFunction{Float64}}() # also includes G for QPs
     objective::Union{Nothing,MOI.AbstractScalarFunction} = nothing
     factorization::Union{Nothing,Function} = nothing
+    allow_objective_and_solution_input::Bool = false
 end
 
 function Base.empty!(cache::InputCache)
@@ -122,6 +123,8 @@ MOI.set(model, DiffOpt.NonLinearKKTJacobianFactorization(), factorization)
 """
 struct NonLinearKKTJacobianFactorization <: MOI.AbstractModelAttribute end
 
+struct AllowObjectiveAndSolutionInput <: MOI.AbstractModelAttribute end
+
 """
     ForwardConstraintFunction <: MOI.AbstractConstraintAttribute
 
@@ -440,6 +443,15 @@ function MOI.set(
     return
 end
 
+function MOI.set(
+    model::AbstractModel,
+    ::AllowObjectiveAndSolutionInput,
+    allow::Bool,
+)
+    model.input_cache.allow_objective_and_solution_input = allow
+    return
+end
+
 function MOI.set(
     model::AbstractModel,
     ::ReverseVariablePrimal,

src/jump_moi_overloads.jl

@@ -29,6 +29,14 @@ function MOI.set(
     return MOI.set(JuMP.backend(model), attr, factorization)
 end
 
+function MOI.set(
+    model::JuMP.Model,
+    attr::AllowObjectiveAndSolutionInput,
+    allow::Bool,
+)
+    return MOI.set(JuMP.backend(model), attr, allow)
+end
+
 function MOI.set(
     model::JuMP.Model,
     attr::ForwardObjectiveFunction,

src/jump_wrapper.jl

@@ -143,3 +143,21 @@ Get the value of a variable output sensitivity for forward mode.
 function get_forward_variable(model::JuMP.Model, variable::JuMP.VariableRef)
     return MOI.get(model, ForwardVariablePrimal(), variable)
 end
+
+"""
+    set_reverse_objective(model::JuMP.Model, value::Number)
+
+Set the value of the objective input sensitivity for reverse mode.
+"""
+function set_reverse_objective(model::JuMP.Model, value::Number)
+    return MOI.set(model, ReverseObjectiveSensitivity(), value)
+end
+
+"""
+    get_forward_objective(model::JuMP.Model)
+
+Get the value of the objective output sensitivity for forward mode.
+"""
+function get_forward_objective(model::JuMP.Model)
+    return MOI.get(model, ForwardObjectiveSensitivity())
+end