Merge pull request #11 from Robbybp/matrix

Add `incidence_matrix` function

Author: Robbybp

Project.toml

@@ -8,6 +8,7 @@ BipartiteMatching = "79040ab4-24c8-4c92-950c-d48b5991a0f6"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
 julia = "1"

docs/make.jl

@@ -32,6 +32,7 @@ makedocs(
             "reference/identify_variables.md",
             "reference/incidence_graph.md",
             "reference/interface.md",
+            "reference/incidence_matrix.md",
         ],
     ],
 )

docs/src/reference/incidence_matrix.md

@@ -0,0 +1,9 @@
+# Incidence matrix
+
+```@meta
+CurrentModule = JuMPIn
+```
+
+```@docs
+incidence_matrix
+```

src/JuMPIn.jl

@@ -33,4 +33,8 @@ include("maximum_matching.jl")
 include("dulmage_mendelsohn.jl")
 include("interface.jl")
 
+# Methods to get incidence matrices (as SparseMatrixCSC) from
+# JuMP models or incidence graphs
+include("incidence_matrix.jl")
+
 end

src/incidence_matrix.jl

@@ -0,0 +1,92 @@
+#  ___________________________________________________________________________
+#
+#  JuMPIn.jl: JuMP Incidence Graph Analysis
+#  Copyright (c) 2023. Triad National Security, LLC. All rights reserved.
+#
+#  This program was produced under U.S. Government contract 89233218CNA000001
+#  for Los Alamos National Laboratory (LANL), which is operated by Triad
+#  National Security, LLC for the U.S. Department of Energy/National Nuclear
+#  Security Administration. All rights in the program are reserved by Triad
+#  National Security, LLC, and the U.S. Department of Energy/National Nuclear
+#  Security Administration. The Government is granted for itself and others
+#  acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license
+#  in this material to reproduce, prepare derivative works, distribute copies
+#  to the public, perform publicly and display publicly, and to permit others
+#  to do so.
+#
+#  This software is distributed under the 3-clause BSD license.
+#  ___________________________________________________________________________
+
+"""
+Utility functions for getting the incidence graph of JuMP constraints and
+variables.
+
+"""
+
+import JuMP
+import SparseArrays
+import Graphs
+
+import JuMPIn: get_bipartite_incidence_graph, IncidenceGraphInterface
+
+"""
+    incidence_matrix(constraints, variables)::SparseMatrixCSC
+
+Return a `SparseMatrixCSC` with entries corresponding to incidence
+between the provided constraints and variables.
+
+Rows correspond to constraints, columns correspond to variables, and
+matrix entries correspond to edges in the incidence graph. All matrix
+entries have values of `1.0`.
+
+"""
+function incidence_matrix(
+    constraints::Vector{<:JuMP.ConstraintRef},
+    variables::Vector{JuMP.VariableRef},
+)::SparseArrays.SparseMatrixCSC
+    graph, _, _ = get_bipartite_incidence_graph(constraints, variables)
+    A, B, E = graph
+    M = length(constraints)
+    row = Vector{Int64}()
+    col = Vector{Int64}()
+    val = Vector{Float64}()
+    for (i, j) in E
+        # NOTE: Here we exploit the graph's convention. This will need to
+        # change if we change how the graph is constructed.
+        push!(row, i)
+        push!(col, j - M)
+        push!(val, 1.0)
+    end
+    return SparseArrays.sparse(row, col, val)
+end
+
+"""
+    incidence_matrix(igraph::IncidenceGraphInterface)::SparseMatrixCSC
+
+Return the incidence matrix associated with the provided
+`IncidenceGraphInterface`.
+
+"""
+function incidence_matrix(
+    igraph::IncidenceGraphInterface
+)::SparseArrays.SparseMatrixCSC
+    row = Vector{Int64}()
+    col = Vector{Int64}()
+    val = Vector{Float64}()
+    M = length(igraph._con_node_map)
+    for e in Graphs.edges(igraph._graph)
+        i = Graphs.src(e)
+        j = Graphs.dst(e)
+        if typeof(igraph._nodes[i]) <: JuMP.VariableRef
+            i, j = j, i
+        end
+        # Now we know that j is a variable node
+        push!(row, i)
+        # NOTE: Here we exploit the graph's convention. This will need to
+        # change if we change how the graph is constructed, e.g. we add
+        # a node for the objective.
+        push!(col, j - M)
+        push!(val, 1.0)
+    end
+    return SparseArrays.sparse(row, col, val)
+end

test/incidence_matrix.jl

@@ -0,0 +1,79 @@
+#  ___________________________________________________________________________
+#
+#  JuMPIn.jl: JuMP Incidence Graph Analysis
+#  Copyright (c) 2023. Triad National Security, LLC. All rights reserved.
+#
+#  This program was produced under U.S. Government contract 89233218CNA000001
+#  for Los Alamos National Laboratory (LANL), which is operated by Triad
+#  National Security, LLC for the U.S. Department of Energy/National Nuclear
+#  Security Administration. All rights in the program are reserved by Triad
+#  National Security, LLC, and the U.S. Department of Energy/National Nuclear
+#  Security Administration. The Government is granted for itself and others
+#  acting on its behalf a nonexclusive, paid-up, irrevocable worldwide license
+#  in this material to reproduce, prepare derivative works, distribute copies
+#  to the public, perform publicly and display publicly, and to permit others
+#  to do so.
+#
+#  This software is distributed under the 3-clause BSD license.
+#  ___________________________________________________________________________
+
+module TestIncidenceMatrix
+
+import JuMP as jmp
+import MathOptInterface as moi
+import SparseArrays
+using Test: @test, @test_throws
+
+import JuMPIn as ji
+
+include("models.jl") # Models
+using .Models: make_degenerate_flow_model
+
+
+function test_incidence_matrix_from_constraints_and_variables()
+    m = jmp.Model()
+    @jmp.variable(m, x[1:3])
+    @jmp.constraint(m, eq1, 2*x[1] + 3*x[2] == 4)
+    @jmp.NLconstraint(m, eq2, 2*x[3]^1.5*x[2] == 1)
+    constraints = [eq1, eq2]
+    variables = [x[1], x[2], x[3]]
+    imat = ji.incidence_matrix(constraints, variables)
+    pred_row = [1, 1, 2, 2]
+    pred_col = [1, 2, 2, 3]
+    pred_val = [1.0, 1.0, 1.0, 1.0]
+    m = 2
+    n = 3
+    pred_mat = SparseArrays.sparse(pred_row, pred_col, pred_val, m, n)
+    @test imat == pred_mat
+end
+
+function test_incidence_matrix_from_incidence_graph()
+    m = jmp.Model()
+    @jmp.variable(m, x[1:3])
+    @jmp.constraint(m, eq1, 2*x[1] + 3*x[2] == 4)
+    @jmp.NLconstraint(m, eq2, 2*x[3]^1.5*x[2] == 1)
+    constraints = [eq1, eq2]
+    variables = [x[3], x[2], x[1]]
+    igraph = ji.IncidenceGraphInterface(constraints, variables)
+    imat = ji.incidence_matrix(igraph)
+    pred_row = [1, 1, 2, 2]
+    pred_col = [3, 2, 2, 1]
+    pred_val = [1.0, 1.0, 1.0, 1.0]
+    m = 2
+    n = 3
+    pred_mat = SparseArrays.sparse(pred_row, pred_col, pred_val, m, n)
+    @test imat == pred_mat
+end
+
+
+function runtests()
+    test_incidence_matrix_from_constraints_and_variables()
+    test_incidence_matrix_from_incidence_graph()
+    return
+end
+
+end # module TestIncidenceGraph
+
+if abspath(PROGRAM_FILE) == @__FILE__
+    TestIncidenceMatrix.runtests()
+end

test/runtests.jl

@@ -42,4 +42,9 @@ using Test
         TestInterface.runtests()
     end
 
+    @testset "IncidenceMatrix" begin
+        include("incidence_matrix.jl")
+        TestIncidenceMatrix.runtests()
+    end
+
 end