Skip to content

Commit

Permalink
Add lp_matrix_data (#3573)
Browse files Browse the repository at this point in the history
  • Loading branch information
odow authored Nov 21, 2023
1 parent 3f59da9 commit 3608379
Show file tree
Hide file tree
Showing 5 changed files with 367 additions and 102 deletions.
66 changes: 66 additions & 0 deletions docs/src/manual/models.md
Original file line number Diff line number Diff line change
Expand Up @@ -439,6 +439,72 @@ optimize!(model)
reduced_cost(x) # Works
```

## Get the matrix representation

Use [`lp_matrix_data`](@ref) to return a data structure that represents the
matrix form of a linear program.

```jldoctest
julia> begin
model = Model()
@variable(model, x >= 1)
@variable(model, 2 <= y)
@variable(model, 3 <= z <= 4)
@constraint(model, x == 5)
@constraint(model, 2x + 3y <= 6)
@constraint(model, -4y >= 5z + 7)
@constraint(model, -1 <= x + y <= 2)
@objective(model, Max, 1 + 2x)
end;
julia> data = lp_matrix_data(model);
julia> data.A
4×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 7 stored entries:
1.0 ⋅ ⋅
⋅ -4.0 -5.0
2.0 3.0 ⋅
1.0 1.0 ⋅
julia> data.b_lower
4-element Vector{Float64}:
5.0
7.0
-Inf
-1.0
julia> data.b_upper
4-element Vector{Float64}:
5.0
Inf
6.0
2.0
julia> data.x_lower
3-element Vector{Float64}:
1.0
2.0
3.0
julia> data.x_upper
3-element Vector{Float64}:
Inf
Inf
4.0
julia> data.c
3-element Vector{Float64}:
2.0
0.0
0.0
julia> data.c_offset
1.0
julia> data.sense
MAX_SENSE::OptimizationSense = 1
```

## Backends

!!! info
Expand Down
1 change: 1 addition & 0 deletions src/JuMP.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1091,6 +1091,7 @@ include("complement.jl")
include("copy.jl")
include("feasibility_checker.jl")
include("file_formats.jl")
include("lp_matrix_data.jl")
include("lp_sensitivity2.jl")
include("indicator.jl")
include("reified.jl")
Expand Down
192 changes: 192 additions & 0 deletions src/lp_matrix_data.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,192 @@
# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at http://mozilla.org/MPL/2.0/.

"""
LPMatrixData{T}
The struct returned by [`lp_matrix_data`](@ref). See [`lp_matrix_data`](@ref)
for a description of the public fields.
"""
struct LPMatrixData{T}
A::SparseArrays.SparseMatrixCSC{T,Int}
b_lower::Vector{T}
b_upper::Vector{T}
x_lower::Vector{T}
x_upper::Vector{T}
c::Vector{T}
c_offset::T
sense::MOI.OptimizationSense
variables::Vector{GenericVariableRef{T}}
affine_constraints::Vector{ConstraintRef}
variable_constraints::Vector{ConstraintRef}
end

"""
lp_matrix_data(model::GenericModel{T})
Given a JuMP model of a linear program, return an [`LPMatrixData{T}`](@ref)
struct storing data for an equivalent linear program in the form:
```math
\\begin{aligned}
\\min & c^\\top x + c_0\\
& b_l \\le A x \\le b_u \\
& x_l \\le x \\le x_u
\\end{aligned}
```
## Fields
The struct returned by [`lp_matrix_data`](@ref) has the fields:
* `A::SparseArrays.SparseMatrixCSC{T,Int}`: the constraint matrix in sparse
matrix form.
* `b_lower::Vector{T}`: the dense vector of row lower bounds. If missing, the
value of `typemin(T)` is used.
* `b_upper::Vector{T}`: the dense vector of row upper bounds. If missing, the
value of `typemax(T)` is used.
* `x_lower::Vector{T}`: the dense vector of variable lower bounds. If missing,
the value of `typemin(T)` is used.
* `x_upper::Vector{T}`: the dense vector of variable upper bounds. If missing,
the value of `typemax(T)` is used.
* `c::Vector{T}`: the dense vector of linear objective coefficiennts
* `c_offset::T`: the constant term in the objective function.
* `sense::MOI.OptimizationSense`: the objective sense of the model.
* `variables::Vector{GenericVariableRef{T}}`: a vector of [`GenericVariableRef`](@ref),
corresponding to order of the columns in the matrix form.
* `affine_constraints::Vector{ConstraintRef}`: a vector of [`ConstraintRef`](@ref),
corresponding to the order of rows in the matrix form.
## Limitations
The models supported by [`lp_matrix_data`](@ref) are intentionally limited
to linear programs.
If your model has integrality, use [`relax_integrality`](@ref) to remove integer
restrictions before calling [`lp_matrix_data`](@ref).
"""
function lp_matrix_data(model::GenericModel{T}) where {T}
variables = all_variables(model)
columns = Dict(var => i for (i, var) in enumerate(variables))
n = length(columns)
cache = (;
x_l = fill(typemin(T), n),
x_u = fill(typemax(T), n),
c = zeros(T, n),
c_offset = Ref{T}(zero(T)),
b_l = T[],
b_u = T[],
I = Int[],
J = Int[],
V = T[],
variable_to_column = columns,
bound_constraints = ConstraintRef[],
affine_constraints = ConstraintRef[],
)
for (F, S) in list_of_constraint_types(model)
_fill_standard_form(model, F, S, cache)
end
_fill_standard_form(model, objective_function_type(model), cache)
return LPMatrixData(
SparseArrays.sparse(cache.I, cache.J, cache.V, length(cache.b_l), n),
cache.b_l,
cache.b_u,
cache.x_l,
cache.x_u,
cache.c,
cache.c_offset[],
MOI.get(model, MOI.ObjectiveSense()),
variables,
cache.affine_constraints,
cache.bound_constraints,
)
end

_bounds(s::MOI.LessThan{T}) where {T} = (typemin(T), s.upper)
_bounds(s::MOI.GreaterThan{T}) where {T} = (s.lower, typemax(T))
_bounds(s::MOI.EqualTo{T}) where {T} = (s.value, s.value)
_bounds(s::MOI.Interval{T}) where {T} = (s.lower, s.upper)

function _fill_standard_form(
model::GenericModel{T},
::Type{GenericVariableRef{T}},
::Type{S},
cache::Any,
) where {
T,
S<:Union{MOI.LessThan{T},MOI.GreaterThan{T},MOI.EqualTo{T},MOI.Interval{T}},
}
for c in all_constraints(model, GenericVariableRef{T}, S)
push!(cache.bound_constraints, c)
c_obj = constraint_object(c)
i = cache.variable_to_column[c_obj.func]
l, u = _bounds(c_obj.set)
cache.x_l[i] = max(cache.x_l[i], l)
cache.x_u[i] = min(cache.x_u[i], u)
end
return
end

function _fill_standard_form(
model::GenericModel{T},
::Type{F},
::Type{S},
cache::Any,
) where {
T,
F<:GenericAffExpr{T,GenericVariableRef{T}},
S<:Union{MOI.LessThan{T},MOI.GreaterThan{T},MOI.EqualTo{T},MOI.Interval{T}},
}
for c in all_constraints(model, F, S)
push!(cache.affine_constraints, c)
c_obj = constraint_object(c)
@assert iszero(c_obj.func.constant)
row = length(cache.b_l) + 1
l, u = _bounds(c_obj.set)
push!(cache.b_l, l)
push!(cache.b_u, u)
for (x, coef) in c_obj.func.terms
push!(cache.I, row)
push!(cache.J, cache.variable_to_column[x])
push!(cache.V, coef)
end
end
return
end

function _fill_standard_form(
::GenericModel{T},
::Type{F},
::Type{S},
::Any,
) where {T,F,S}
return error("Unsupported constraint type in `lp_matrix_data`: $F -in- $S")
end

function _fill_standard_form(
model::GenericModel{T},
::Type{GenericVariableRef{T}},
cache::Any,
) where {T}
cache.c[cache.variable_to_column[objective_function(model)]] = one(T)
cache.c_offset[] = zero(T)
return
end

function _fill_standard_form(
model::GenericModel{T},
::Type{GenericAffExpr{T,GenericVariableRef{T}}},
cache::Any,
) where {T}
f = objective_function(model)
for (k, v) in f.terms
cache.c[cache.variable_to_column[k]] += v
end
cache.c_offset[] = f.constant
return
end

function _fill_standard_form(::GenericModel{T}, ::Type{F}, ::Any) where {T,F}
return error("Unsupported objective type in `lp_matrix_data`: $F")
end
111 changes: 9 additions & 102 deletions src/lp_sensitivity2.jl
Original file line number Diff line number Diff line change
Expand Up @@ -325,114 +325,21 @@ end
"""
_standard_form_matrix(model::GenericModel)
Given a problem:
r_l <= Ax <= r_u
c_l <= x <= c_u
Return the standard form:
[A -I] [x, y] = 0
[c_l, r_l] <= [x, y] <= [c_u, r_u]
`columns` maps the variable references to column indices.
See [`lp_matrix_data`](@ref) instead.
"""
function _standard_form_matrix(model::GenericModel{T}) where {T}
columns = Dict(var => i for (i, var) in enumerate(all_variables(model)))
n = length(columns)
c_l, c_u = fill(typemin(T), n), fill(typemax(T), n)
r_l, r_u = T[], T[]
I, J, V = Int[], Int[], T[]
bound_constraints = ConstraintRef[]
affine_constraints = ConstraintRef[]
for (F, S) in list_of_constraint_types(model)
_fill_standard_form(
model,
columns,
bound_constraints,
affine_constraints,
F,
S,
c_l,
c_u,
r_l,
r_u,
I,
J,
V,
)
end
matrix = lp_matrix_data(model)
I = SparseArrays.spdiagm(fill(-one(T), length(matrix.affine_constraints)))
return (
columns = columns,
lower = vcat(c_l, r_l),
upper = vcat(c_u, r_u),
A = SparseArrays.sparse(I, J, V, length(r_l), n + length(r_l)),
bounds = bound_constraints,
constraints = affine_constraints,
columns = Dict(x => i for (i, x) in enumerate(matrix.variables)),
lower = vcat(matrix.x_lower, matrix.b_lower),
upper = vcat(matrix.x_upper, matrix.b_upper),
A = hcat(matrix.A, I),
bounds = matrix.variable_constraints,
constraints = matrix.affine_constraints,
)
end

function _fill_standard_form(
model::GenericModel{T},
x::Dict{GenericVariableRef{T},Int},
bound_constraints::Vector{ConstraintRef},
::Vector{ConstraintRef},
F::Type{GenericVariableRef{T}},
S::Type,
c_l::Vector{T},
c_u::Vector{T},
::Vector{T},
::Vector{T},
::Vector{Int},
::Vector{Int},
::Vector{T},
) where {T}
for c in all_constraints(model, F, S)
push!(bound_constraints, c)
c_obj = constraint_object(c)
i = x[c_obj.func]
set = MOI.Interval(c_obj.set)
c_l[i] = max(c_l[i], set.lower)
c_u[i] = min(c_u[i], set.upper)
end
return
end

function _fill_standard_form(
model::GenericModel{T},
x::Dict{GenericVariableRef{T},Int},
::Vector{ConstraintRef},
affine_constraints::Vector{ConstraintRef},
F::Type{<:GenericAffExpr},
S::Type,
::Vector{T},
::Vector{T},
r_l::Vector{T},
r_u::Vector{T},
I::Vector{Int},
J::Vector{Int},
V::Vector{T},
) where {T}
for c in all_constraints(model, F, S)
push!(affine_constraints, c)
c_obj = constraint_object(c)
@assert iszero(c_obj.func.constant)
row = length(r_l) + 1
set = MOI.Interval(c_obj.set)
push!(r_l, set.lower)
push!(r_u, set.upper)
for (var, coef) in c_obj.func.terms
push!(I, row)
push!(J, x[var])
push!(V, coef)
end
push!(I, row)
push!(J, length(x) + row)
push!(V, -one(T))
end
return
end

_convert_nonbasic_status(::MOI.LessThan) = MOI.NONBASIC_AT_UPPER
_convert_nonbasic_status(::MOI.GreaterThan) = MOI.NONBASIC_AT_LOWER
_convert_nonbasic_status(::Any) = MOI.NONBASIC
Expand Down
Loading

0 comments on commit 3608379

Please sign in to comment.