Add back qp_test

src/ConicProgram/ConicProgram.jl

@@ -398,9 +398,12 @@ function DiffOpt._get_db(
     i = MOI.Utilities.rows(model.model.constraints, ci) # vector
     # i = ci.value
     n = length(model.x) # columns in A
+    # Since `b` in https://arxiv.org/pdf/1904.09043.pdf is the constant in the right-hand side and
+    # `b` in MOI is the constant on the left-hand side, we have the opposite sign here
     # db = - dQ[n+1:n+m, end] + dQ[end, n+1:n+m]'
     g = model.back_grad_cache.g
     πz = model.back_grad_cache.πz
+    # `g[end] * πz[n .+ i] - πz[end] * g[n .+ i]`
     return DiffOpt.lazy_combination(-, πz, g, length(g), n .+ i)
 end
 
@@ -415,7 +418,8 @@ function DiffOpt._get_dA(
     # dA = - dQ[1:n, n+1:n+m]' + dQ[n+1:n+m, 1:n]
     g = model.back_grad_cache.g
     πz = model.back_grad_cache.πz
-    return DiffOpt.lazy_combination(-, g, πz, i, n .+ (1:n))
+    #return DiffOpt.lazy_combination(-, g, πz, n .+ i, 1:n)
+    return g[n.+i] * πz[1:n]' - πz[n.+i] * g[1:n]'
 end
 
 function MOI.get(

src/bridges.jl

@@ -17,6 +17,15 @@ function MOI.set(
     )
 end
 
+function MOI.get(
+    model::MOI.ModelLike,
+    attr::ReverseConstraintFunction,
+    bridge::MOI.Bridges.Constraint.VectorizeBridge,
+)
+    return MOI.Utilities.eachscalar(
+        MOI.get(model, attr, bridge.vector_constraint),
+    )[1]
+end
 function MOI.get(
     model::MOI.ModelLike,
     attr::DiffOpt.ReverseConstraintFunction,

src/jump_moi_overloads.jl

@@ -131,82 +131,6 @@ function Base.isapprox(
     return isapprox(standard_form(func1), standard_form(func2); kws...)
 end
 
-# In the future, we could replace by https://github.com/jump-dev/MathOptInterface.jl/pull/1238
-"""
-    VectorScalarAffineFunction{T, VT} <: MOI.AbstractScalarFunction
-
-Represents the function `x ⋅ terms + constant`
-as an `MOI.AbstractScalarFunction` where `x[i] = MOI.VariableIndex(i)`.
-Use [`standard_form`](@ref) to convert it to a `MOI.ScalarAffineFunction{T}`.
-"""
-struct VectorScalarAffineFunction{T,VT} <: MOI.AbstractScalarFunction
-    terms::VT
-    constant::T
-end
-
-MOI.constant(func::VectorScalarAffineFunction) = func.constant
-
-function JuMP.coefficient(
-    func::VectorScalarAffineFunction,
-    vi::MOI.VariableIndex,
-)
-    return func.terms[vi.value]
-end
-
-function Base.convert(
-    ::Type{MOI.ScalarAffineFunction{T}},
-    func::VectorScalarAffineFunction,
-) where {T}
-    return MOI.ScalarAffineFunction{T}(
-        # TODO we should do better if the vector is a `SparseVector`, I think
-        #      I have some code working for both vector types in Polyhedra.jl
-        MOI.ScalarAffineTerm{T}[
-            MOI.ScalarAffineTerm{T}(func.terms[i], MOI.VariableIndex(i)) for
-            i in eachindex(func.terms) if !iszero(func.terms[i])
-        ],
-        func.constant,
-    )
-end
-
-function standard_form(func::VectorScalarAffineFunction{T}) where {T}
-    return convert(MOI.ScalarAffineFunction{T}, func)
-end
-
-function MOI.Utilities.operate(
-    ::typeof(-),
-    ::Type{T},
-    func::VectorScalarAffineFunction{T},
-) where {T}
-    return VectorScalarAffineFunction(
-        LazyArrays.ApplyArray(-, func.terms),
-        -func.constant,
-    )
-end
-
-"""
-    struct MatrixScalarQuadraticFunction{T, VT, MT} <: MOI.AbstractScalarFunction
-        affine::VectorScalarAffineFunction{T,VT}
-        terms::MT
-    end
-
-Represents the function `x' * terms * x / 2 + affine` as an
-`MOI.AbstractScalarFunction` where `x[i] = MOI.VariableIndex(i)`.
-Use [`standard_form`](@ref) to convert it to a `MOI.ScalarQuadraticFunction{T}`.
-"""
-struct MatrixScalarQuadraticFunction{T,VT,MT} <: MOI.AbstractScalarFunction
-    affine::VectorScalarAffineFunction{T,VT}
-    terms::MT
-end
-
-MOI.constant(func::MatrixScalarQuadraticFunction) = MOI.constant(func.affine)
-
-function JuMP.coefficient(
-    func::MatrixScalarQuadraticFunction,
-    vi::MOI.VariableIndex,
-)
-    return JuMP.coefficient(func.affine, vi)
-end
-
 function quad_sym_half(
     func::MatrixScalarQuadraticFunction,
     vi1::MOI.VariableIndex,
@@ -250,50 +174,6 @@ function standard_form(func::MatrixScalarQuadraticFunction{T}) where {T}
     return convert(MOI.ScalarQuadraticFunction{T}, func)
 end
 
-"""
-    MatrixVectorAffineFunction{T, VT} <: MOI.AbstractVectorFunction
-
-Represents the function `terms * x + constant`
-as an `MOI.AbstractVectorFunction` where `x[i] = MOI.VariableIndex(i)`.
-Use [`standard_form`](@ref) to convert it to a `MOI.VectorAffineFunction{T}`.
-"""
-struct MatrixVectorAffineFunction{AT,VT} <: MOI.AbstractVectorFunction
-    terms::AT
-    constants::VT
-end
-
-MOI.constant(func::MatrixVectorAffineFunction) = func.constants
-
-function Base.convert(
-    ::Type{MOI.VectorAffineFunction{T}},
-    func::MatrixVectorAffineFunction,
-) where {T}
-    return MOI.VectorAffineFunction{T}(
-        MOI.VectorAffineTerm{T}[
-            # TODO we should do better if the matrix is a `SparseMatrixCSC`
-            MOI.VectorAffineTerm(
-                i,
-                MOI.ScalarAffineTerm{T}(func.terms[i, j], MOI.VariableIndex(j)),
-            ) for i in 1:size(func.terms, 1) for
-            j in 1:size(func.terms, 2) if !iszero(func.terms[i, j])
-        ],
-        func.constants,
-    )
-end
-
-function standard_form(func::MatrixVectorAffineFunction{T}) where {T}
-    return convert(MOI.VectorAffineFunction{T}, func)
-end
-
-# Only used for testing at the moment so performance is not critical so
-# converting to standard form is ok
-function MOI.Utilities.isapprox_zero(
-    func::Union{VectorScalarAffineFunction,MatrixScalarQuadraticFunction},
-    tol,
-)
-    return MOI.Utilities.isapprox_zero(standard_form(func), tol)
-end
-
 """
     IndexMappedFunction{F<:MOI.AbstractFunction} <: AbstractLazyScalarFunction
 

src/utils.jl

@@ -152,6 +152,127 @@ function sparse_array_representation(
     )
 end
 
+# In the future, we could replace by https://github.com/jump-dev/MathOptInterface.jl/pull/1238
+"""
+    VectorScalarAffineFunction{T, VT} <: MOI.AbstractScalarFunction
+
+Represents the function `x ⋅ terms + constant`
+as an `MOI.AbstractScalarFunction` where `x[i] = MOI.VariableIndex(i)`.
+Use [`standard_form`](@ref) to convert it to a `MOI.ScalarAffineFunction{T}`.
+"""
+struct VectorScalarAffineFunction{T,VT} <: MOI.AbstractScalarFunction
+    terms::VT
+    constant::T
+end
+MOI.constant(func::VectorScalarAffineFunction) = func.constant
+function JuMP.coefficient(
+    func::VectorScalarAffineFunction,
+    vi::MOI.VariableIndex,
+)
+    return func.terms[vi.value]
+end
+function Base.convert(
+    ::Type{MOI.ScalarAffineFunction{T}},
+    func::VectorScalarAffineFunction,
+) where {T}
+    return MOI.ScalarAffineFunction{T}(
+        # TODO we should do better if the vector is a `SparseVector`, I think
+        #      I have some code working for both vector types in Polyhedra.jl
+        MOI.ScalarAffineTerm{T}[
+            MOI.ScalarAffineTerm{T}(func.terms[i], MOI.VariableIndex(i)) for
+            i in eachindex(func.terms) if !iszero(func.terms[i])
+        ],
+        func.constant,
+    )
+end
+function standard_form(func::VectorScalarAffineFunction{T}) where {T}
+    return convert(MOI.ScalarAffineFunction{T}, func)
+end
+
+function MOI.Utilities.operate(
+    ::typeof(-),
+    ::Type{T},
+    func::VectorScalarAffineFunction{T},
+) where {T}
+    return VectorScalarAffineFunction(
+        LazyArrays.ApplyArray(-, func.terms),
+        -func.constant,
+    )
+end
+
+"""
+    struct MatrixScalarQuadraticFunction{T, VT, MT} <: MOI.AbstractScalarFunction
+        affine::VectorScalarAffineFunction{T,VT}
+        terms::MT
+    end
+
+Represents the function `x' * terms * x / 2 + affine` as an
+`MOI.AbstractScalarFunction` where `x[i] = MOI.VariableIndex(i)`.
+Use [`standard_form`](@ref) to convert it to a `MOI.ScalarQuadraticFunction{T}`.
+"""
+struct MatrixScalarQuadraticFunction{T,VT,MT} <: MOI.AbstractScalarFunction
+    affine::VectorScalarAffineFunction{T,VT}
+    terms::MT
+end
+MOI.constant(func::MatrixScalarQuadraticFunction) = MOI.constant(func.affine)
+function JuMP.coefficient(
+    func::MatrixScalarQuadraticFunction,
+    vi::MOI.VariableIndex,
+)
+    return JuMP.coefficient(func.affine, vi)
+end
+
+"""
+    MatrixVectorAffineFunction{T, VT} <: MOI.AbstractVectorFunction
+
+Represents the function `terms * x + constant`
+as an `MOI.AbstractVectorFunction` where `x[i] = MOI.VariableIndex(i)`.
+Use [`standard_form`](@ref) to convert it to a `MOI.VectorAffineFunction{T}`.
+"""
+struct MatrixVectorAffineFunction{AT,VT} <: MOI.AbstractVectorFunction
+    terms::AT
+    constants::VT
+end
+MOI.constant(func::MatrixVectorAffineFunction) = func.constants
+function Base.convert(
+    ::Type{MOI.VectorAffineFunction{T}},
+    func::MatrixVectorAffineFunction,
+) where {T}
+    return MOI.VectorAffineFunction{T}(
+        MOI.VectorAffineTerm{T}[
+            # TODO we should do better if the matrix is a `SparseMatrixCSC`
+            MOI.VectorAffineTerm(
+                i,
+                MOI.ScalarAffineTerm{T}(func.terms[i, j], MOI.VariableIndex(j)),
+            ) for i in 1:size(func.terms, 1) for
+            j in 1:size(func.terms, 2) if !iszero(func.terms[i, j])
+        ],
+        func.constants,
+    )
+end
+function standard_form(func::MatrixVectorAffineFunction{T}) where {T}
+    return convert(MOI.VectorAffineFunction{T}, func)
+end
+
+# Only used for testing at the moment so performance is not critical so
+# converting to standard form is ok
+function MOIU.isapprox_zero(
+    func::Union{VectorScalarAffineFunction,MatrixScalarQuadraticFunction},
+    tol,
+)
+    return MOIU.isapprox_zero(standard_form(func), tol)
+end
+
+_scalar(::Type{<:MatrixVectorAffineFunction}) = VectorScalarAffineFunction
+_scalar(::Type{<:SparseVectorAffineFunction}) = SparseScalarAffineFunction
+
+function Base.getindex(
+    it::MOI.Utilities.ScalarFunctionIterator{F},
+    output_index::Integer,
+) where {F<:Union{MatrixVectorAffineFunction,SparseVectorAffineFunction}}
+    return _scalar(F)(it.f.terms[output_index, :], it.f.constants[output_index])
+end
+
 function _index_map_to_oneto!(index_map, v::MOI.VariableIndex)
     if !haskey(index_map, v)
         n = length(index_map.var_map)

test/linear_program.jl

@@ -12,8 +12,8 @@ import Ipopt
 import MathOptInterface as MOI
 import SCS
 
-const ATOL = 2e-4
-const RTOL = 2e-4
+const ATOL = 1e-2
+const RTOL = 1e-2
 
 function runtests()
     for name in names(@__MODULE__; all = true)

test/quadratic_program.jl

@@ -322,6 +322,7 @@ function test_differentiating_non_trivial_convex_qp_moi()
     dbb = vec(grads_actual[6])
     qp_test(
         Ipopt.Optimizer,
+        DiffOpt.ConicProgram.Model,
         true,
         true,
         true;