From d924cb3a5b9b1ccfbeb8ebb13ba9b47103fbca54 Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Thu, 30 Oct 2025 23:39:17 -0400 Subject: [PATCH 1/7] start implementation of limited BFS --- src/interface.jl | 71 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) diff --git a/src/interface.jl b/src/interface.jl index 2b4662f..ba604e7 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -713,3 +713,74 @@ end block_triangularize(model::JuMP.Model) = block_triangularize(IncidenceGraphInterface(model)) block_triangularize(matrix::SparseMatrixCSC) = block_triangularize(IncidenceGraphInterface(matrix)) block_triangularize(matrix::Matrix) = block_triangularize(IncidenceGraphInterface(matrix)) + +const NodeType = Union{JuMP.VariableRef,JuMP.ConstraintRef,Int} + +struct IncidenceSubtree{T} + _nodes::Vector{NodeType} + _dag::Graphs.DiGraph{T} +end + +function _collect_lines!(lines::Vector, tree::IncidenceSubtree; node = 1, indent = "") + indent = (indent == "" ? "|- " : join(["| ", indent])) + for i in Graphs.neighbors(tree._dag, node) + child = tree._nodes[i] + push!(lines, "$indent$child") + _collect_lines!(lines, tree; node = i, indent) + end + return +end + +function Base.show(io::IO, tree::IncidenceSubtree) + root = tree._nodes[1] + lines = ["$root"] + _collect_lines!(lines, tree) + for line in lines + println(io, line) + end + return +end + +function limited_bfs( + igraph::IncidenceGraphInterface, + root::NodeType; + depth::Int = 1, +) + root = if root in keys(igraph._var_node_map) + igraph._var_node_map[root] + elseif root in keys(igraph._con_node_map) + igraph._con_node_map[root] + else + error("root is not a node in this graph") + end + # What is a natural data structure to return as a BFS tree? + # - An adjacency list? + # - First I just need assemble the set of vertices traversed + g = igraph._graph + nodes = NodeType[root] + idx = 1 + # I either have to use a length-n array or use a Set/Dict + parents = fill(-1, length(igraph._nodes)) + parents[root] = 0 + # FIXME: This explores at most the first depth nodes + for i in 1:depth + node = nodes[idx] + parent_idx = idx + index_updated = false + for other in Graphs.neighbors(g, node) + if parents[other] == -1 + parents[other] = parent_idx + push!(nodes, other) + if !index_updated + idx += 1 + index_updated = true + end + end + end + end + parents = map(n -> parents[n], nodes) + digraph = Graphs.tree(parents) + subtree_nodes = convert(Vector{NodeType}, map(n -> igraph._nodes[n], nodes)) + tree = IncidenceSubtree(subtree_nodes, digraph) + return tree +end From f147fc08fb90817854e0a6cf1381d1a91bb26bfa Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Fri, 31 Oct 2025 10:29:00 -0400 Subject: [PATCH 2/7] fix limited bfs implementation --- src/interface.jl | 46 +++++++++++++++++++++++++++++++--------------- 1 file changed, 31 insertions(+), 15 deletions(-) diff --git a/src/interface.jl b/src/interface.jl index ba604e7..08f5873 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -753,30 +753,46 @@ function limited_bfs( else error("root is not a node in this graph") end - # What is a natural data structure to return as a BFS tree? - # - An adjacency list? - # - First I just need assemble the set of vertices traversed g = igraph._graph - nodes = NodeType[root] + # This is the queue + nodes = Int[root] idx = 1 - # I either have to use a length-n array or use a Set/Dict - parents = fill(-1, length(igraph._nodes)) - parents[root] = 0 - # FIXME: This explores at most the first depth nodes - for i in 1:depth + # I use a dict rather than an array here to avoid constructing an O(n) + # array for every (potentially) small BFS in a (potentially) large graph + parents = Dict(root => 0) + # We differ from traditional BFS by imposing a depth limit. We store the index + # (in the queue) corresponding to the start of the current "level". + maxdepth = depth + current_depth = 0 + levelstart = 1 + while idx <= length(nodes) node = nodes[idx] parent_idx = idx - index_updated = false + idx += 1 for other in Graphs.neighbors(g, node) - if parents[other] == -1 + if other ∉ keys(parents) + if parent_idx >= levelstart + # When we first encounter the child of a node in the current level, + # we advance (a) the depth and (b) the level start index. + # The level starts at the child node we are about to add. We do the + # checking here so we can break before adding the child if doing + # so would violate the depth limit. + levelstart = length(nodes) + 1 + current_depth += 1 + if current_depth > maxdepth + # Break before adding this child to the queue + break + end + end parents[other] = parent_idx push!(nodes, other) - if !index_updated - idx += 1 - index_updated = true - end end end + # Once we exceed the depth limit, we must stop processing children + # of future nodes. + if current_depth > maxdepth + break + end end parents = map(n -> parents[n], nodes) digraph = Graphs.tree(parents) From fa4cbef6949a15f4336886c58af050601fa615c5 Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Fri, 31 Oct 2025 13:10:16 -0400 Subject: [PATCH 3/7] move lower-level BFS implementation to its own file --- src/MathProgIncidence.jl | 1 + src/bfs.jl | 66 +++++++++++++++++++++++++ src/interface.jl | 101 +++++++++++++++++++++------------------ 3 files changed, 122 insertions(+), 46 deletions(-) create mode 100644 src/bfs.jl diff --git a/src/MathProgIncidence.jl b/src/MathProgIncidence.jl index 350fdf5..6fb2355 100644 --- a/src/MathProgIncidence.jl +++ b/src/MathProgIncidence.jl @@ -32,6 +32,7 @@ include("incidence_graph.jl") include("maximum_matching.jl") include("dulmage_mendelsohn.jl") include("block_triangularize.jl") +include("bfs.jl") include("interface.jl") # Methods to get incidence matrices (as SparseMatrixCSC) from diff --git a/src/bfs.jl b/src/bfs.jl new file mode 100644 index 0000000..cb3a2b0 --- /dev/null +++ b/src/bfs.jl @@ -0,0 +1,66 @@ +# ___________________________________________________________________________ +# +# MathProgIncidence.jl: Math Programming 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. +# ___________________________________________________________________________ + +import Graphs + +function _limited_bfs(adjacency_list::Vector{Vector{Int}}, root::Int; depth::Int = 1) + # This is the queue + nodes = Int[root] + idx = 1 + # I use a dict rather than an array here to avoid constructing an O(n) + # array for every (potentially) small BFS in a (potentially) large graph + parents = Dict(root => 0) + # We differ from traditional BFS by imposing a depth limit. We store the index + # (in the queue) corresponding to the start of the current "level". + maxdepth = depth + current_depth = 0 + levelstart = 1 + while idx <= length(nodes) + node = nodes[idx] + parent_idx = idx + idx += 1 + for other in adjacency_list[node] + if other ∉ keys(parents) + if parent_idx >= levelstart + # When we first encounter the child of a node in the current level, + # we advance (a) the depth and (b) the level start index. + # The level starts at the child node we are about to add. We do the + # checking here so we can break before adding the child if doing + # so would violate the depth limit. + levelstart = length(nodes) + 1 + current_depth += 1 + if current_depth > maxdepth + # Break before adding this child to the queue + break + end + end + parents[other] = parent_idx + push!(nodes, other) + end + end + # Once we exceed the depth limit, we must stop processing children + # of future nodes. + if current_depth > maxdepth + break + end + end + parents = map(n -> parents[n], nodes) + dag = Graphs.tree(parents) # Returns a DiGraph + return nodes, dag +end diff --git a/src/interface.jl b/src/interface.jl index 08f5873..2eb0d58 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -714,7 +714,16 @@ block_triangularize(model::JuMP.Model) = block_triangularize(IncidenceGraphInter block_triangularize(matrix::SparseMatrixCSC) = block_triangularize(IncidenceGraphInterface(matrix)) block_triangularize(matrix::Matrix) = block_triangularize(IncidenceGraphInterface(matrix)) -const NodeType = Union{JuMP.VariableRef,JuMP.ConstraintRef,Int} +const NodeType = Union{ + JuMP.VariableRef, + JuMP.ConstraintRef, + Int, + # This is starting to get a little excessive. + # Should I just use Vector{Any}? + JuMP.AffExpr, + JuMP.QuadExpr, + JuMP.NonlinearExpr, +} struct IncidenceSubtree{T} _nodes::Vector{NodeType} @@ -753,50 +762,50 @@ function limited_bfs( else error("root is not a node in this graph") end - g = igraph._graph - # This is the queue - nodes = Int[root] - idx = 1 - # I use a dict rather than an array here to avoid constructing an O(n) - # array for every (potentially) small BFS in a (potentially) large graph - parents = Dict(root => 0) - # We differ from traditional BFS by imposing a depth limit. We store the index - # (in the queue) corresponding to the start of the current "level". - maxdepth = depth - current_depth = 0 - levelstart = 1 - while idx <= length(nodes) - node = nodes[idx] - parent_idx = idx - idx += 1 - for other in Graphs.neighbors(g, node) - if other ∉ keys(parents) - if parent_idx >= levelstart - # When we first encounter the child of a node in the current level, - # we advance (a) the depth and (b) the level start index. - # The level starts at the child node we are about to add. We do the - # checking here so we can break before adding the child if doing - # so would violate the depth limit. - levelstart = length(nodes) + 1 - current_depth += 1 - if current_depth > maxdepth - # Break before adding this child to the queue - break - end - end - parents[other] = parent_idx - push!(nodes, other) - end - end - # Once we exceed the depth limit, we must stop processing children - # of future nodes. - if current_depth > maxdepth - break - end - end - parents = map(n -> parents[n], nodes) - digraph = Graphs.tree(parents) - subtree_nodes = convert(Vector{NodeType}, map(n -> igraph._nodes[n], nodes)) - tree = IncidenceSubtree(subtree_nodes, digraph) + adjlist = Graphs.SimpleGraphs.adj(igraph._graph) + nodes, dag = _limited_bfs(adjlist, root; depth) + nodes = convert(Vector{NodeType}, map(n -> igraph._nodes[n], nodes)) + tree = IncidenceSubtree(nodes, dag) return tree end + +function limited_bfs(root::Union{JuMP.VariableRef,<:JuMP.ConstraintRef}; depth::Int = 1) + igraph = IncidenceGraphInterface(root.model) + return limited_bfs(igraph, root; depth) +end + +function limited_bfs( + igraph::IncidenceGraphInterface, + root::Union{ + JuMP.AffExpr, + JuMP.QuadExpr, + JuMP.NonlinearExpr, + }; + depth::Int = 1, +) + nnodes = length(igraph._nodes) + # Here, I assume the nodes are contiguous integers + expr_node = nnodes + 1 + # Get the expression's neighbors. These must be variables. + expr_neighbors = identify_unique_variables(JuMP.moi_function(root)) + # FIXME: I should just have an identify_unique_variables method that returns + # VariableRef from expressions. This is just a hack in the meantime. + model = igraph._nodes[1].model + # Convert VariableIndex to VariableRef + expr_neighbors = map(i -> JuMP.VariableRef(model, i), expr_neighbors) + expr_neighbors = map(v -> igraph._var_node_map[v], expr_neighbors) + g = igraph._graph + adjlist = map(n -> Graphs.neighbors(g, n), 1:nnodes) + # Add the new node's edges to our adjacency list + push!(adjlist, expr_neighbors) + subtree_nodes, dag = _limited_bfs(adjlist, expr_node; depth) + # We know that root (or its corresponding integer, expr_node) is the first + # node in subtree_nodes. We need to use this because we can't look it up + # in igraph._nodes. + jump_subtree_nodes = NodeType[root] + append!( + jump_subtree_nodes, + map(n -> igraph._nodes[n], subtree_nodes[2:end]), + ) + return IncidenceSubtree(jump_subtree_nodes, dag) +end From 2dc5eebe6dddcec83064fef36c040ef418d31736 Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Fri, 31 Oct 2025 13:14:26 -0400 Subject: [PATCH 4/7] typos --- src/identify_variables.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/identify_variables.jl b/src/identify_variables.jl index 1911b9e..db56482 100644 --- a/src/identify_variables.jl +++ b/src/identify_variables.jl @@ -199,7 +199,7 @@ function identify_unique_variables( end """ - identify_unique_variables(fcn)::Vector{JuMP.VariableIndex} + identify_unique_variables(fcn)::Vector{MOI.VariableIndex} Return the variables that appear in the provided MathOptInterface function. @@ -224,7 +224,7 @@ function identify_unique_variables( return _filter_duplicates(variables) end -# This method is used to handle function-in-set constraints. +# This method is used to handle variable-in-set constraints. function identify_unique_variables( var::MOI.VariableIndex )::Vector{MOI.VariableIndex} From e206ca0087f813d880d8f51a89c15933f4d0a6f7 Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Fri, 31 Oct 2025 14:15:23 -0400 Subject: [PATCH 5/7] limited_bfs methods for expressions --- src/identify_variables.jl | 33 ++++++++++++++++++++++++++ src/interface.jl | 50 +++++++++++++++++++++++---------------- 2 files changed, 62 insertions(+), 21 deletions(-) diff --git a/src/identify_variables.jl b/src/identify_variables.jl index db56482..1e77125 100644 --- a/src/identify_variables.jl +++ b/src/identify_variables.jl @@ -248,6 +248,39 @@ function identify_unique_variables( )) end +function identify_unique_variables(expr::JuMP.AffExpr)::Vector{JuMP.VariableRef} + return _filter_duplicates(collect(keys(expr.terms))) +end + +function identify_unique_variables(expr::JuMP.QuadExpr)::Vector{JuMP.VariableRef} + vars = collect(keys(expr.aff.terms)) + for term in keys(expr.terms) + push!(vars, term.a) + push!(vars, term.b) + end + return _filter_duplicates(vars) +end + +function identify_unique_variables(expr::JuMP.VariableRef)::Vector{JuMP.VariableRef} + return [expr] +end + +function identify_unique_variables(expr::T)::Vector{JuMP.VariableRef} where {T<:Real} + return [] +end + +function identify_unique_variables(expr::JuMP.NonlinearExpr)::Vector{JuMP.VariableRef} + vars = JuMP.VariableRef[] + for arg in expr.args + # NOTE: This could be made more efficient by avoiding the construction + # of many small vectors. However, this method is probably not critical + # for overall performance as of 2025-10-31. + append!(vars, identify_unique_variables(arg)) + end + return _filter_duplicates(vars) +end + + """ _identify_variables(fcn)::Vector{JuMP.VariableIndex} diff --git a/src/interface.jl b/src/interface.jl index 2eb0d58..7d2fb52 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -714,7 +714,7 @@ block_triangularize(model::JuMP.Model) = block_triangularize(IncidenceGraphInter block_triangularize(matrix::SparseMatrixCSC) = block_triangularize(IncidenceGraphInterface(matrix)) block_triangularize(matrix::Matrix) = block_triangularize(IncidenceGraphInterface(matrix)) -const NodeType = Union{ +const SubtreeNodeType = Union{ JuMP.VariableRef, JuMP.ConstraintRef, Int, @@ -726,7 +726,7 @@ const NodeType = Union{ } struct IncidenceSubtree{T} - _nodes::Vector{NodeType} + _nodes::Vector{SubtreeNodeType} _dag::Graphs.DiGraph{T} end @@ -752,7 +752,7 @@ end function limited_bfs( igraph::IncidenceGraphInterface, - root::NodeType; + root::Union{JuMP.VariableRef,<:JuMP.ConstraintRef,Int}; depth::Int = 1, ) root = if root in keys(igraph._var_node_map) @@ -764,35 +764,29 @@ function limited_bfs( end adjlist = Graphs.SimpleGraphs.adj(igraph._graph) nodes, dag = _limited_bfs(adjlist, root; depth) - nodes = convert(Vector{NodeType}, map(n -> igraph._nodes[n], nodes)) - tree = IncidenceSubtree(nodes, dag) - return tree + nodes = convert(Vector{SubtreeNodeType}, map(n -> igraph._nodes[n], nodes)) + return IncidenceSubtree(nodes, dag) end -function limited_bfs(root::Union{JuMP.VariableRef,<:JuMP.ConstraintRef}; depth::Int = 1) - igraph = IncidenceGraphInterface(root.model) +function limited_bfs( + root::Union{JuMP.VariableRef,<:JuMP.ConstraintRef}; + depth::Int = 1, + kwds..., +) + igraph = IncidenceGraphInterface(root.model; kwds...) return limited_bfs(igraph, root; depth) end function limited_bfs( igraph::IncidenceGraphInterface, - root::Union{ - JuMP.AffExpr, - JuMP.QuadExpr, - JuMP.NonlinearExpr, - }; + root::Union{JuMP.AffExpr,JuMP.QuadExpr,JuMP.NonlinearExpr}; depth::Int = 1, ) nnodes = length(igraph._nodes) # Here, I assume the nodes are contiguous integers expr_node = nnodes + 1 - # Get the expression's neighbors. These must be variables. - expr_neighbors = identify_unique_variables(JuMP.moi_function(root)) - # FIXME: I should just have an identify_unique_variables method that returns - # VariableRef from expressions. This is just a hack in the meantime. - model = igraph._nodes[1].model - # Convert VariableIndex to VariableRef - expr_neighbors = map(i -> JuMP.VariableRef(model, i), expr_neighbors) + # Get the expression's neighbors. These must be variables + expr_neighbors = identify_unique_variables(root) expr_neighbors = map(v -> igraph._var_node_map[v], expr_neighbors) g = igraph._graph adjlist = map(n -> Graphs.neighbors(g, n), 1:nnodes) @@ -802,10 +796,24 @@ function limited_bfs( # We know that root (or its corresponding integer, expr_node) is the first # node in subtree_nodes. We need to use this because we can't look it up # in igraph._nodes. - jump_subtree_nodes = NodeType[root] + jump_subtree_nodes = SubtreeNodeType[root] append!( jump_subtree_nodes, map(n -> igraph._nodes[n], subtree_nodes[2:end]), ) return IncidenceSubtree(jump_subtree_nodes, dag) end + +function limited_bfs( + root::Union{JuMP.AffExpr, JuMP.QuadExpr, JuMP.NonlinearExpr}; + depth::Int = 1, + kwds..., +) + expr_neighbors = identify_unique_variables(root) + if isempty(expr_neighbors) + return IncidenceSubtree(SubtreeNodeType[root], Graphs.DiGraph(1)) + else + igraph = IncidenceGraphInterface(expr_neighbors[1].model; kwds...) + return limited_bfs(igraph, root; depth) + end +end From c020965392d64af51d433e3be6be5fdef173589b Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Fri, 31 Oct 2025 14:30:56 -0400 Subject: [PATCH 6/7] slighly nicer printing --- src/interface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/interface.jl b/src/interface.jl index 7d2fb52..d0cb7f2 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -731,7 +731,7 @@ struct IncidenceSubtree{T} end function _collect_lines!(lines::Vector, tree::IncidenceSubtree; node = 1, indent = "") - indent = (indent == "" ? "|- " : join(["| ", indent])) + indent = (indent == "" ? "├ " : join(["│ ", indent])) for i in Graphs.neighbors(tree._dag, node) child = tree._nodes[i] push!(lines, "$indent$child") From 7332ff35fa978b4a29f278849e844303f43ef2b2 Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Fri, 31 Oct 2025 14:46:12 -0400 Subject: [PATCH 7/7] maybe slightly better printing --- src/interface.jl | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/interface.jl b/src/interface.jl index d0cb7f2..ae2e87e 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -732,9 +732,16 @@ end function _collect_lines!(lines::Vector, tree::IncidenceSubtree; node = 1, indent = "") indent = (indent == "" ? "├ " : join(["│ ", indent])) - for i in Graphs.neighbors(tree._dag, node) + orig_indent = indent + neighbors = Graphs.neighbors(tree._dag, node) + for i in neighbors child = tree._nodes[i] + if i == last(neighbors) + indent = replace(indent, "├" => "└") + end push!(lines, "$indent$child") + # This could be handled better... + indent = orig_indent _collect_lines!(lines, tree; node = i, indent) end return