Algorithm interfaces

IncidenceGraphInterface(
    model;
    include_inequality = false,
    include_active_inequalities = false,
    tolerance = 0.0,
)

A bipartite incidence graph of JuMP constraints and variables.

This is the primary data type accepted by the algorithms implemented in the remainder of this module. This type can be instantiated with a JuMP model or a tuple of (graph, con_node_map, var_node_map), as returned by get_bipartite_incidence_graph. If a model is provided, optional arguments are the same as those provided to get_bipartite_incidence_graph.

Note that the fields of this struct are private, and may change behavior in a future release without warning.

Example using only equality constraints

using JuMP
import MathProgIncidence
m = Model()
@variable(m, v[1:3] >= 0)
@constraint(m, eq_1, v[1] + v[3]^2 == 1.0)
@NLconstraint(m, eq_2, v[1]*v[2]^1.5 == 2.0)
graph = MathProgIncidence.IncidenceGraphInterface(m)

Example including active inequality constraints

using JuMP
import Ipopt
import MathProgIncidence
m = Model(Ipopt.Optimizer)
@variable(m, v[1:3] >= 0)
@NLconstraint(m, eq_1, v[1]*v[2]^1.5 == 2.0)
@constraint(m, ineq_1, v[1] + v[2] + v[3] >= 7)
@objective(m, Min, v[1]^2 + 2*v[2]^2 + 3*v[3]^2)
optimize!(m)
graph = MathProgIncidence.IncidenceGraphInterface(
    m, include_active_inequalities = true, tolerance = 1e-6
)

get_adjacent(
    igraph::IncidenceGraphInterface,
    constraint::JuMP.ConstriantRef,
)::Vector{JuMP.VariableRef}

Return the variables adjacent to a constraint in an incidence graph.

get_adjacent(
    igraph::IncidenceGraphInterface,
    variable::JuMP.VariableRef,
)::Vector{JuMP.ConstraintRef}

Return the constraints adjacent to a variable in an incidence graph.

Example

julia> using JuMP

julia> import MathProgIncidence

julia> m = Model();

julia> @variable(m, v[1:3]);

julia> @constraint(m, eq_1, v[1] + v[3] == 1);

julia> @NLconstraint(m, eq_2, v[1]*v[2]^3 == 2);

julia> igraph = MathProgIncidence.IncidenceGraphInterface(m);

julia> adj_cons = MathProgIncidence.get_adjacent(igraph, v[1]);

julia> display(adj_cons)
2-element Vector{ConstraintRef}:
 eq_1 : v[1] + v[3] = 1.0
 v[1] * v[2] ^ 3.0 - 2.0 = 0

maximum_matching(igraph::IncidenceGraphInterface)::Dict

Compute a maximum matching of variables and constraints in the incidence graph.

The returned Dict maps JuMP ConstraintRefs to their matched VariableRefs.

Example

julia> using JuMP

julia> import MathProgIncidence

julia> m = Model();

julia> @variable(m, v[1:3]);

julia> @constraint(m, eq_1, v[1] + v[3] == 1);

julia> @NLconstraint(m, eq_2, v[1]*v[2]^3 == 2);

julia> igraph = MathProgIncidence.IncidenceGraphInterface(m);

julia> matching = MathProgIncidence.maximum_matching(igraph);

julia> display(matching)
Dict{ConstraintRef, VariableRef} with 2 entries:
  v[1] * v[2] ^ 3.0 - 2.0 = 0 => v[2]
  eq_1 : v[1] + v[3] = 1.0 => v[1]

dulmage_mendelsohn(igraph::IncidenceGraphInterface)

Return the Dulmage-Mendelsohn partition of variables and constraints in an incidence graph.

The returned type is a Tuple of two NamedTuples, (con_dmp, var_dmp). These NamedTuples have the following fields:

con_dmp:

• underconstrained – The constraints matched with variables that can possibly be unmatched in a maximum cardinality matching
• square – The constraints that cannot possibly be unmatched in a maximum matching
• overconstrained – The constraints that are matched, but can possibly be unmatched in a maximum matching
• unmatched – The constraints that are not matched in the maximum matching that was found

var_dmp:

• unmatched – The variables that are not matched in the maximum matching that was found
• underconstrained – The variables that can possibly be unmatched in a maximum matching
• square – The variables that cannot possibly be unmatched in a maximum matching
• overconstrained – The variables matched with constraints that can possibly be unmatched in a maximum cardinality matching

The Dulmage-Mendelsohn partition groups nodes in a bipartite graph into three unique subsets. In the application to constraints and variables, these may be thought of as:

• the "overconstrained subsystem", which has more constraints than variables,
• the "underconstrained subsystem", which has more variables than constraints,
• and the "square subsystem", which has the same number of variables as constraints

In the NamedTuples returned by this function, the constraints in the overconstrained subsystem are split into overconstrained and unmatched, while the variables in the underconstrained subsystem are split into underconstrained and unmatched. This is because it is useful to explicitly check whether there are any unmatched variables and constraints, and also useful to recover the maximum matching by zip-ing corresponding variables and constraints.

Example

julia> using JuMP

julia> import MathProgIncidence

julia> m = Model();

julia> @variable(m, v[1:4]);

julia> @constraint(m, eq_1, v[1] + v[3] == 1);

julia> @NLconstraint(m, eq_2, v[1]*v[2]^3 == 2);

julia> @constraint(m, eq_3, v[4]^2 == 3);

julia> igraph = MathProgIncidence.IncidenceGraphInterface(m);

julia> con_dmp, var_dmp = MathProgIncidence.dulmage_mendelsohn(igraph);

julia> # Assert that there are no unmatched constraints

julia> @assert isempty(con_dmp.unmatched);

julia> display(var_dmp.unmatched)
1-element Vector{VariableRef}:
 v[3]

julia> display(var_dmp.underconstrained)
2-element Vector{VariableRef}:
 v[1]
 v[2]

julia> display(con_dmp.underconstrained)
2-element Vector{ConstraintRef}:
 eq_1 : v[1] + v[3] = 1.0
 v[1] * v[2] ^ 3.0 - 2.0 = 0
 
julia> display(var_dmp.square)
1-element Vector{VariableRef}:
 v[4]

julia> display(con_dmp.square)
1-element Vector{ConstraintRef}:
 eq_3 : v[4]² = 3.0

julia> # As there are no unmatched constraints, the overconstrained subsystem is empty

DulmageMendelsohnDecomposition <: NamedTuple

The return type of dulmage_mendelsohn. Fieldnames are con and var, which each contain the fieldnames documented by dulmage_mendelsohn.

connected_components(igraph::IncidenceGraphInterface)

Return the connected components of a bipartite incidence graph of constraints and variables.

The connected components are returned as a NamedTuple with fields con and var, each a vector-of-vectors containing the constraints/rows and variables/columns in each connected component. Note that the input graph is undirected, so there is no distinction between strongly and weakly connected components.

Example

julia> using JuMP

julia> import MathProgIncidence

julia> m = Model();

julia> @variable(m, x[1:2] >= 0);

julia> @constraint(m, eq1, x[1] == 1);

julia> @constraint(m, eq2, x[2]^2 == 2);

julia> igraph = MathProgIncidence.IncidenceGraphInterface(m);

julia> cc = MathProgIncidence.connected_components(igraph);

julia> cc.con
2-element Vector{Vector{ConstraintRef}}:
 [eq1 : x[1] = 1]
 [eq2 : x[2]² = 2]

julia> cc.var
2-element Vector{Vector{VariableRef}}:
 [x[1]]
 [x[2]]

connected_components(constraints, variables)

Return the connected components of a bipartite incidence graph of constraints and variables.

The method that accepts constraints and variables directly is convenient for working with the output of the Dulmage-Mendelsohn partition. It is often used to decompose and help debug the over and under-constrained subsystems.

Example

julia> using JuMP

julia> import MathProgIncidence

julia> m = Model();

julia> @variable(m, x[1:4] >= 0);

julia> @constraint(m, eq1, x[1] + x[3] == 7);

julia> @constraint(m, eq2, x[2]^2 + x[4]^2 == 1);

julia> igraph = MathProgIncidence.IncidenceGraphInterface(m);

julia> con_dmp, var_dmp = MathProgIncidence.dulmage_mendelsohn(igraph);

julia> uc_con = con_dmp.underconstrained;

julia> uc_var = [var_dmp.unmatched..., var_dmp.underconstrained...];

julia> cc = MathProgIncidence.connected_components(uc_con, uc_var);

julia> cc.con
2-element Vector{Vector{ConstraintRef}}:
 [eq1 : x[1] + x[3] = 7]
 [eq2 : x[2]² + x[4]² = 1]

julia> cc.var
2-element Vector{Vector{VariableRef}}:
 [x[3], x[1]]
 [x[4], x[2]]

ConnectedComponentDecomposition <: NamedTuple

The return type of connected_components. Fields are rows and columns, which each contain a vector of vectors of either constraints/variables or integers.

block_triangularize(igraph::IncidenceGraphInterface)::Vector{Tuple{Vector, Vector}}

Return an ordered partition of constraints and variables that puts the incidence matrix into block-lower triangular form.

The provided constraints and variables must be well-determined. That is, a maximum matching must cover all constraints and variables provided.

The return type is a vector of tuples of vectors of constraints and variables.

The subsets of the partition correspond to diagonal blocks of the incidence matrix. These subsets are strongly connected components of a directed version of the bipartite incidence graph and have the strong Hall property (are irreducible). The block-triangular ordering of these subsets corresponds to a topological order of the directed acyclic graph of strongly connected components.

Re-ordering the constraints and variables into the order given will yield an incidence matrix that is block-lower triangular. As with all incidence matrices in this package, rows correspond to constraints and columns correspond to variables. This means that the subsets are provided in a topological order of the precedence graph, where each subset only uses the variables of those that precede it (as well as its own variables).

Example

julia> using JuMP; import MathProgIncidence as MPIN

julia> m = Model(); @variable(m, x[1:3] >= 0);

julia> c1 = @constraint(m, x[1] + x[2]^2 + x[3]^3 == 4);

julia> c2 = @constraint(m, x[1] + x[3] == 1);

julia> c3 = @constraint(m, x[1] * x[3]^0.5 == 3);

julia> igraph = MPIN.IncidenceGraphInterface(m);

julia> blocks = MPIN.block_triangularize(igraph);

julia> length(blocks)
2

julia> for (i, (cb, vb)) in enumerate(blocks)
           println()
           println("Block $i")
           println("-------")
           for (c, v) in zip(cb, vb) # cb and vb have the same length
               println("$v, $c")
           end
       end

Block 1
-------
x[1], (x[1] * (x[3] ^ 0.5)) - 3.0 = 0
x[3], x[1] + x[3] = 1

Block 2
-------
x[2], ((x[2]² + x[1]) + (x[3] ^ 3.0)) - 4.0 = 0


julia> corder = JuMP.ConstraintRef[]; vorder = JuMP.VariableRef[];

julia> for (cb, vb) in blocks
           append!(corder, cb)
           append!(vorder, vb)
       end

julia> MPIN.incidence_matrix(corder, vorder)
3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 7 stored entries:
 1.0  1.0   ⋅
 1.0  1.0   ⋅
 1.0  1.0  1.0

block_triangularize(constraints, variables)::Vector{Tuple{Vector, Vector}}

Return an ordered partition of constraints and variables that puts the incidence matrix into block-lower triangular form.

This method accepts vectors of constraints and variables and is useful for performing the block-triangular decomposition on the well-constrained subsystem from the Dulmage-Mendelsohn decomposition.

The return type is a vector of Subsystems.

Example

julia> using JuMP; import MathProgIncidence as MPIN

julia> m = Model(); @variable(m, x[1:5] >= 0);

julia> c1 = @constraint(m, x[1] + x[2]^2 + x[3]^3 + x[4]^4 +x[5]^5 == 5);

julia> c2 = @constraint(m, x[1] == 1);

julia> c3 = @constraint(m, x[1] * x[3]^0.5 == 3);

julia> c4 = @constraint(m, x[1] - x[2] + x[3] - x[4] == 10);

julia> igraph = MPIN.IncidenceGraphInterface(m);

julia> cdmp, vdmp = MPIN.dulmage_mendelsohn(igraph);

julia> blocks = MPIN.block_triangularize(cdmp.square, vdmp.square);

julia> length(blocks)
2

julia> for (i, (cb, vb)) in enumerate(blocks)
           println("Block $i")
           println("-------")
           for (c, v) in zip(cb, vb) # cb and vb have the same length
               println("$v, $c")
           end
           println()
       end
Block 1
-------
x[1], x[1] = 1

Block 2
-------
x[3], (x[1] * (x[3] ^ 0.5)) - 3.0 = 0

Subsystem <: NamedTuple

A set of constraints and variables (or row and column indices). Fieldnames are con and var, each of which contains a vector of constraints/variables or integers.

limited_bfs(
    igraph::IncidenceGraphInterface, root; depth = 1
)::IncidenceSubtree

Perform a limited-depth breadth-first search (BFS) starting from the root variable, constraint, or expression. The resulting BFS tree is returned as an IncidenceSubtree, which can be displayed on the console.

Example

julia> using JuMP; import MathProgIncidence as MPIN

julia> model = Model(); @variable(model, x[1:5] >= 0);

julia> c1 = @constraint(model, x[1] + x[2]^2 + x[3]^3 + x[4]^4 +x[5]^5 == 5);

julia> c2 = @constraint(model, x[1] == 1);

julia> c3 = @constraint(model, x[1] * x[3]^0.5 == 3);

julia> c4 = @constraint(model, x[1] - x[2] + x[3] - x[4] == 10);

julia> igraph = MPIN.IncidenceGraphInterface(model)

julia> MPIN.limited_bfs(igraph, x[3])
x[3]
├ ((x[2]² + x[1]) + (x[5] ^ 5) + (x[4] ^ 4) + (x[3] ^ 3)) - 5 = 0
├ (x[1] * (x[3] ^ 0.5)) - 3 = 0
└ x[1] - x[2] + x[3] - x[4] = 10

julia> MPIN.limited_bfs(igraph, x[3]; depth = 2)
x[3]
├ ((x[2]² + x[1]) + (x[5] ^ 5) + (x[4] ^ 4) + (x[3] ^ 3)) - 5 = 0
│ ├ x[2]
│ ├ x[1]
│ ├ x[5]
│ └ x[4]
├ (x[1] * (x[3] ^ 0.5)) - 3 = 0
└ x[1] - x[2] + x[3] - x[4] = 10

limited_bfs(root; depth = 1, kwds...)

Convenience methods for limited_bfs that only require the root as a JuMP variable, constraint, or expression. The model is automatically inferred.

IncidenceSubtree

The return type of limited_bfs. It contains a vector of nodes (JuMP variables, constraints, or expressions) and Graphs.jl DiGraph containing the edges.