From 85e6a0fae2af1962703dee066ee0362ebcde5066 Mon Sep 17 00:00:00 2001 From: robbybp Date: Mon, 30 Jan 2023 22:41:41 -0700 Subject: [PATCH 01/25] run black on incidence_analysis --- pyomo/contrib/incidence_analysis/__init__.py | 6 +- .../common/dulmage_mendelsohn.py | 31 ++-- pyomo/contrib/incidence_analysis/connected.py | 8 +- .../incidence_analysis/dulmage_mendelsohn.py | 40 ++--- pyomo/contrib/incidence_analysis/interface.py | 137 +++++++++--------- pyomo/contrib/incidence_analysis/matching.py | 13 +- .../incidence_analysis/triangularize.py | 23 ++- pyomo/contrib/incidence_analysis/util.py | 59 ++++---- 8 files changed, 158 insertions(+), 159 deletions(-) diff --git a/pyomo/contrib/incidence_analysis/__init__.py b/pyomo/contrib/incidence_analysis/__init__.py index 178bbc96db8..364046d51da 100644 --- a/pyomo/contrib/incidence_analysis/__init__.py +++ b/pyomo/contrib/incidence_analysis/__init__.py @@ -2,6 +2,6 @@ from .matching import maximum_matching from .interface import IncidenceGraphInterface from .util import ( - generate_strongly_connected_components, - solve_strongly_connected_components, - ) + generate_strongly_connected_components, + solve_strongly_connected_components, +) diff --git a/pyomo/contrib/incidence_analysis/common/dulmage_mendelsohn.py b/pyomo/contrib/incidence_analysis/common/dulmage_mendelsohn.py index 347538659d4..b861cce356a 100644 --- a/pyomo/contrib/incidence_analysis/common/dulmage_mendelsohn.py +++ b/pyomo/contrib/incidence_analysis/common/dulmage_mendelsohn.py @@ -27,8 +27,7 @@ def _get_projected_digraph(bg, matching, top_nodes): - """ - """ + """ """ digraph = DiGraph() digraph.add_nodes_from(top_nodes) for n in top_nodes: @@ -45,8 +44,7 @@ def _get_projected_digraph(bg, matching, top_nodes): def _get_reachable_from(digraph, sources): - """ - """ + """ """ _filter = set() reachable = [] for node in sources: @@ -96,14 +94,17 @@ def dulmage_mendelsohn(bg, top_nodes=None, matching=None): t_other = [t for t in top_nodes if t not in _filter] b_other = [b for b in bot_nodes if b not in _filter] - return (( - t_unmatched, - t_reachable, - t_matched_with_reachable, - t_other, - ), ( - b_unmatched, - b_reachable, - b_matched_with_reachable, - b_other, - )) + return ( + ( + t_unmatched, + t_reachable, + t_matched_with_reachable, + t_other, + ), + ( + b_unmatched, + b_reachable, + b_matched_with_reachable, + b_other, + ), + ) diff --git a/pyomo/contrib/incidence_analysis/connected.py b/pyomo/contrib/incidence_analysis/connected.py index f268e292991..67b7e0b7392 100644 --- a/pyomo/contrib/incidence_analysis/connected.py +++ b/pyomo/contrib/incidence_analysis/connected.py @@ -13,8 +13,7 @@ def get_independent_submatrices(matrix): - """ - """ + """ """ nxc = nx.algorithms.components nxb = nx.algorithms.bipartite from_biadjacency_matrix = nxb.matrix.from_biadjacency_matrix @@ -27,11 +26,10 @@ def get_independent_submatrices(matrix): # nodes have values in [N, N+M-1]. # We could also check the "bipartite" attribute of each node... row_blocks = [ - sorted([node for node in comp if node < N]) - for comp in connected_components + sorted([node for node in comp if node < N]) for comp in connected_components ] col_blocks = [ - sorted([node-N for node in comp if node >= N]) + sorted([node - N for node in comp if node >= N]) for comp in connected_components ] return row_blocks, col_blocks diff --git a/pyomo/contrib/incidence_analysis/dulmage_mendelsohn.py b/pyomo/contrib/incidence_analysis/dulmage_mendelsohn.py index 761b6ee69c2..64d77c64266 100644 --- a/pyomo/contrib/incidence_analysis/dulmage_mendelsohn.py +++ b/pyomo/contrib/incidence_analysis/dulmage_mendelsohn.py @@ -13,20 +13,22 @@ from pyomo.common.dependencies import networkx as nx from pyomo.contrib.incidence_analysis.common.dulmage_mendelsohn import ( dulmage_mendelsohn as dm_nx, - ) +) + """ This module imports the general Dulmage-Mendelsohn-on-a-graph function from "common" and implements an interface for coo_matrix-like objects. """ RowPartition = namedtuple( - "RowPartition", - ["unmatched", "overconstrained", "underconstrained", "square"], - ) + "RowPartition", + ["unmatched", "overconstrained", "underconstrained", "square"], +) ColPartition = namedtuple( - "ColPartition", - ["unmatched", "underconstrained", "overconstrained", "square"], - ) + "ColPartition", + ["unmatched", "underconstrained", "overconstrained", "square"], +) + def dulmage_mendelsohn(matrix_or_graph, top_nodes=None, matching=None): """ @@ -45,9 +47,9 @@ def dulmage_mendelsohn(matrix_or_graph, top_nodes=None, matching=None): graph = matrix_or_graph if top_nodes is None: raise ValueError( - "top_nodes must be specified if a graph is provided," - "\notherwise the result is ambiguous." - ) + "top_nodes must be specified if a graph is provided," + "\notherwise the result is ambiguous." + ) partition = dm_nx(graph, top_nodes=top_nodes, matching=matching) # RowPartition and ColPartition do not make sense for a general graph. # However, here we assume that this graph comes from a Pyomo model, @@ -74,17 +76,17 @@ def dulmage_mendelsohn(matrix_or_graph, top_nodes=None, matching=None): # Matrix rows have bipartite=0, columns have bipartite=1 bg = from_biadjacency_matrix(matrix) row_partition, col_partition = dm_nx( - bg, - top_nodes=list(range(M)), - matching=matching, - ) + bg, + top_nodes=list(range(M)), + matching=matching, + ) partition = ( - row_partition, - tuple([n - M for n in subset] for subset in col_partition) - # Column nodes have values in [M, M+N-1]. Apply the offset - # to get values corresponding to indices in user's matrix. - ) + row_partition, + tuple([n - M for n in subset] for subset in col_partition) + # Column nodes have values in [M, M+N-1]. Apply the offset + # to get values corresponding to indices in user's matrix. + ) partition = (RowPartition(*partition[0]), ColPartition(*partition[1])) return partition diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py index edb4f370424..8ccd88c9c2a 100644 --- a/pyomo/contrib/incidence_analysis/interface.py +++ b/pyomo/contrib/incidence_analysis/interface.py @@ -29,12 +29,13 @@ block_triangularize, get_diagonal_blocks, get_blocks_from_maps, - ) +) from pyomo.contrib.incidence_analysis.dulmage_mendelsohn import ( dulmage_mendelsohn, RowPartition, ColPartition, - ) +) + if scipy_available: from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP import scipy as sp @@ -50,10 +51,9 @@ def _check_unindexed(complist): for comp in complist: if comp.is_indexed(): raise ValueError( - "Variables and constraints must be unindexed " - "ComponentData objects. Got %s, which is indexed." - % comp.name - ) + "Variables and constraints must be unindexed " + "ComponentData objects. Got %s, which is indexed." % comp.name + ) def get_incidence_graph(variables, constraints, include_fixed=True): @@ -72,14 +72,14 @@ def get_incidence_graph(variables, constraints, include_fixed=True): Returns: -------- NetworkX Graph - + """ - _check_unindexed(variables+constraints) + _check_unindexed(variables + constraints) N, M = len(variables), len(constraints) graph = nx.Graph() graph.add_nodes_from(range(M), bipartite=0) - graph.add_nodes_from(range(M, M+N), bipartite=1) - var_node_map = ComponentMap((v, M+i) for i, v in enumerate(variables)) + graph.add_nodes_from(range(M, M + N), bipartite=1) + var_node_map = ComponentMap((v, M + i) for i, v in enumerate(variables)) for i, con in enumerate(constraints): for var in identify_variables(con.expr, include_fixed=include_fixed): if var in var_node_map: @@ -115,19 +115,21 @@ def get_structural_incidence_matrix(variables, constraints, include_fixed=True): Entries are 1.0. """ - _check_unindexed(variables+constraints) + _check_unindexed(variables + constraints) N, M = len(variables), len(constraints) var_idx_map = ComponentMap((v, i) for i, v in enumerate(variables)) rows = [] cols = [] for i, con in enumerate(constraints): - cols.extend(var_idx_map[v] for v in - identify_variables(con.expr, include_fixed=include_fixed) - if v in var_idx_map) - rows.extend([i]*(len(cols) - len(rows))) + cols.extend( + var_idx_map[v] + for v in identify_variables(con.expr, include_fixed=include_fixed) + if v in var_idx_map + ) + rows.extend([i] * (len(cols) - len(rows))) assert len(rows) == len(cols) - data = [1.0]*len(rows) - matrix = sp.sparse.coo_matrix( (data, (rows, cols)), shape=(M, N) ) + data = [1.0] * len(rows) + matrix = sp.sparse.coo_matrix((data, (rows, cols)), shape=(M, N)) return matrix @@ -155,14 +157,13 @@ class IncidenceGraphInterface(object): """ def __init__( - self, - model=None, - active=True, - include_fixed=False, - include_inequality=True, - ): - """ - """ + self, + model=None, + active=True, + include_fixed=False, + include_inequality=True, + ): + """ """ # If the user gives us a model or an NLP, we assume they want us # to cache the incidence matrix for fast analysis of submatrices # later on. @@ -188,9 +189,11 @@ def __init__( self.variables = nlp.get_pyomo_variables() self.constraints = nlp.get_pyomo_constraints() self.var_index_map = ComponentMap( - (var, idx) for idx, var in enumerate(self.variables)) + (var, idx) for idx, var in enumerate(self.variables) + ) self.con_index_map = ComponentMap( - (con, idx) for idx, con in enumerate(self.constraints)) + (con, idx) for idx, con in enumerate(self.constraints) + ) if include_inequality: self.incidence_matrix = nlp.evaluate_jacobian() else: @@ -198,8 +201,8 @@ def __init__( elif isinstance(model, Block): self.cached = IncidenceMatrixType.STRUCTURAL self.constraints = [ - con for con in - model.component_data_objects(Constraint, active=active) + con + for con in model.component_data_objects(Constraint, active=active) if include_inequality or isinstance(con.expr, EqualityExpression) ] self.variables = list( @@ -208,19 +211,20 @@ def __init__( ) ) self.var_index_map = ComponentMap( - (var, i) for i, var in enumerate(self.variables)) + (var, i) for i, var in enumerate(self.variables) + ) self.con_index_map = ComponentMap( - (con, i) for i, con in enumerate(self.constraints)) + (con, i) for i, con in enumerate(self.constraints) + ) self.incidence_matrix = get_structural_incidence_matrix( - self.variables, - self.constraints, - ) + self.variables, + self.constraints, + ) else: raise TypeError( "Unsupported type for incidence matrix. Expected " - "%s or %s but got %s." - % (PyomoNLP, Block, type(model)) - ) + "%s or %s but got %s." % (PyomoNLP, Block, type(model)) + ) self.row_block_map = None self.col_block_map = None @@ -228,36 +232,34 @@ def __init__( def _validate_input(self, variables, constraints): if variables is None: if self.cached is IncidenceMatrixType.NONE: - raise ValueError( - "Neither variables nor a model have been provided." - ) + raise ValueError("Neither variables nor a model have been provided.") else: variables = self.variables if constraints is None: if self.cached is IncidenceMatrixType.NONE: - raise ValueError( - "Neither constraints nor a model have been provided." - ) + raise ValueError("Neither constraints nor a model have been provided.") else: constraints = self.constraints - _check_unindexed(variables+constraints) + _check_unindexed(variables + constraints) return variables, constraints def _extract_submatrix(self, variables, constraints): # Assumes variables and constraints are valid if self.cached is IncidenceMatrixType.NONE: return get_structural_incidence_matrix( - variables, - constraints, - include_fixed=False, - ) + variables, + constraints, + include_fixed=False, + ) else: N, M = len(variables), len(constraints) - old_new_var_indices = dict((self.var_index_map[v], i) - for i, v in enumerate(variables)) - old_new_con_indices = dict((self.con_index_map[c], i) - for i, c in enumerate(constraints)) + old_new_var_indices = dict( + (self.var_index_map[v], i) for i, v in enumerate(variables) + ) + old_new_con_indices = dict( + (self.con_index_map[c], i) for i, c in enumerate(constraints) + ) coo = self.incidence_matrix new_row = [] new_col = [] @@ -268,9 +270,9 @@ def _extract_submatrix(self, variables, constraints): new_col.append(old_new_var_indices[c]) new_data.append(e) return sp.sparse.coo_matrix( - (new_data, (new_row, new_col)), - shape=(M, N), - ) + (new_data, (new_row, new_col)), + shape=(M, N), + ) def maximum_matching(self, variables=None, constraints=None): """ @@ -283,8 +285,7 @@ def maximum_matching(self, variables=None, constraints=None): matching = maximum_matching(matrix.tocoo()) # Matching maps row (constraint) indices to column (variable) indices - return ComponentMap((constraints[i], variables[j]) - for i, j in matching.items()) + return ComponentMap((constraints[i], variables[j]) for i, j in matching.items()) def get_connected_components(self, variables=None, constraints=None): """ @@ -317,10 +318,12 @@ def block_triangularize(self, variables=None, constraints=None): # future. self.row_block_map = row_block_map self.col_block_map = col_block_map - con_block_map = ComponentMap((constraints[i], idx) - for i, idx in row_block_map.items()) - var_block_map = ComponentMap((variables[j], idx) - for j, idx in col_block_map.items()) + con_block_map = ComponentMap( + (constraints[i], idx) for i, idx in row_block_map.items() + ) + var_block_map = ComponentMap( + (variables[j], idx) for j, idx in col_block_map.items() + ) # Switch the order of the maps here to match the method call. # Hopefully this does not get too confusing... return var_block_map, con_block_map @@ -370,11 +373,11 @@ def dulmage_mendelsohn(self, variables=None, constraints=None): row_partition, col_partition = dulmage_mendelsohn(matrix.tocoo()) con_partition = RowPartition( - *[[constraints[i] for i in subset] for subset in row_partition] - ) + *[[constraints[i] for i in subset] for subset in row_partition] + ) var_partition = ColPartition( - *[[variables[i] for i in subset] for subset in col_partition] - ) + *[[variables[i] for i in subset] for subset in col_partition] + ) # Switch the order of the maps here to match the method call. # Hopefully this does not get too confusing... return var_partition, con_partition @@ -411,9 +414,7 @@ def remove_nodes(self, nodes, constraints=None): to_exclude.update(constraints) vars_to_include = [v for v in self.variables if v not in to_exclude] cons_to_include = [c for c in self.constraints if c not in to_exclude] - incidence_matrix = self._extract_submatrix( - vars_to_include, cons_to_include - ) + incidence_matrix = self._extract_submatrix(vars_to_include, cons_to_include) # update attributes self.variables = vars_to_include self.constraints = cons_to_include diff --git a/pyomo/contrib/incidence_analysis/matching.py b/pyomo/contrib/incidence_analysis/matching.py index 25906958103..d2a66e92d53 100644 --- a/pyomo/contrib/incidence_analysis/matching.py +++ b/pyomo/contrib/incidence_analysis/matching.py @@ -11,9 +11,10 @@ from pyomo.common.dependencies import networkx as nx + def maximum_matching(matrix): """ - Returns a maximum matching of the rows and columns of the + Returns a maximum matching of the rows and columns of the matrix as a dict from row indices to column indices. """ nxb = nx.algorithms.bipartite @@ -26,11 +27,11 @@ def maximum_matching(matrix): # Check assumptions regarding from_biadjacency_matrix function: for i in range(M): # First M nodes in graph correspond to rows - assert bg.nodes[i]['bipartite'] == 0 + assert bg.nodes[i]["bipartite"] == 0 - for j in range(M, M+N): + for j in range(M, M + N): # Last N nodes in graph correspond to columns - assert bg.nodes[j]['bipartite'] == 1 + assert bg.nodes[j]["bipartite"] == 1 # If the matrix is block diagonal, the graph will be disconnected. # This is fine, but we need to separate into connected components @@ -38,7 +39,5 @@ def maximum_matching(matrix): conn_comp = [bg.subgraph(c) for c in nxc.connected_components(bg)] matchings = [nxb.maximum_matching(c) for c in conn_comp] # If n0 < M, then n1 >= M. n0 is the row index, n1-M is the column index - max_matching = { - n0: n1-M for m in matchings for n0, n1 in m.items() if n0 < M - } + max_matching = {n0: n1 - M for m in matchings for n0, n1 in m.items() if n0 < M} return max_matching diff --git a/pyomo/contrib/incidence_analysis/triangularize.py b/pyomo/contrib/incidence_analysis/triangularize.py index 7c5d8865c52..356c6e0e9aa 100644 --- a/pyomo/contrib/incidence_analysis/triangularize.py +++ b/pyomo/contrib/incidence_analysis/triangularize.py @@ -39,10 +39,10 @@ def block_triangularize(matrix, matching=None): M, N = matrix.shape if M != N: - raise ValueError("block_triangularize does not currently " - "support non-square matrices. Got matrix with shape %s." - % (matrix.shape,) - ) + raise ValueError( + "block_triangularize does not currently " + "support non-square matrices. Got matrix with shape %s." % (matrix.shape,) + ) bg = from_biadjacency_matrix(matrix) if matching is None: @@ -50,10 +50,11 @@ def block_triangularize(matrix, matching=None): len_matching = len(matching) if len_matching != M: - raise ValueError("block_triangularize only supports matrices " - "that have a perfect matching of rows and columns. " - "Cardinality of maximal matching is %s" % len_matching - ) + raise ValueError( + "block_triangularize only supports matrices " + "that have a perfect matching of rows and columns. " + "Cardinality of maximal matching is %s" % len_matching + ) # Construct directed graph of rows dg = nx.DiGraph() @@ -109,7 +110,7 @@ def get_blocks_from_maps(row_block_map, col_block_map): Arguments --------- row_block_map: dict - Dict mapping each row coordinate to the coordinate of the + Dict mapping each row coordinate to the coordinate of the block it belongs to col_block_map: dict @@ -160,8 +161,6 @@ def get_diagonal_blocks(matrix, matching=None): diagonal blocks. """ - row_block_map, col_block_map = block_triangularize( - matrix, matching=matching - ) + row_block_map, col_block_map = block_triangularize(matrix, matching=matching) block_rows, block_cols = get_blocks_from_maps(row_block_map, col_block_map) return block_rows, block_cols diff --git a/pyomo/contrib/incidence_analysis/util.py b/pyomo/contrib/incidence_analysis/util.py index 9e592dad922..0e5a7d05cf9 100644 --- a/pyomo/contrib/incidence_analysis/util.py +++ b/pyomo/contrib/incidence_analysis/util.py @@ -15,19 +15,19 @@ from pyomo.core.expr.visitor import identify_variables from pyomo.util.calc_var_value import calculate_variable_from_constraint from pyomo.util.subsystems import ( - create_subsystem_block, - TemporarySubsystemManager, - generate_subsystem_blocks, - ) + create_subsystem_block, + TemporarySubsystemManager, + generate_subsystem_blocks, +) from pyomo.contrib.incidence_analysis.interface import IncidenceGraphInterface def generate_strongly_connected_components( - constraints, - variables=None, - include_fixed=False, - ): - """ Performs a block triangularization of the incidence matrix + constraints, + variables=None, + include_fixed=False, +): + """Performs a block triangularization of the incidence matrix of the provided constraints and variables, and yields a block that contains the constraints and variables of each diagonal block (strongly connected component). @@ -55,9 +55,9 @@ def generate_strongly_connected_components( variables = [] for con in constraints: for var in identify_variables( - con.expr, - include_fixed=include_fixed, - ): + con.expr, + include_fixed=include_fixed, + ): if var not in var_set: variables.append(var) var_set.add(var) @@ -65,9 +65,9 @@ def generate_strongly_connected_components( assert len(variables) == len(constraints) igraph = IncidenceGraphInterface() var_block_map, con_block_map = igraph.block_triangularize( - variables=variables, - constraints=constraints, - ) + variables=variables, + constraints=constraints, + ) blocks = set(var_block_map.values()) n_blocks = len(blocks) var_blocks = [[] for b in range(n_blocks)] @@ -78,21 +78,21 @@ def generate_strongly_connected_components( con_blocks[b].append(con) subsets = list(zip(con_blocks, var_blocks)) for block, inputs in generate_subsystem_blocks( - subsets, - include_fixed=include_fixed, - ): + subsets, + include_fixed=include_fixed, + ): # TODO: How does len scale for reference-to-list? assert len(block.vars) == len(block.cons) yield (block, inputs) def solve_strongly_connected_components( - block, - solver=None, - solve_kwds=None, - calc_var_kwds=None, - ): - """ This function solves a square block of variables and equality + block, + solver=None, + solve_kwds=None, + calc_var_kwds=None, +): + """This function solves a square block of variables and equality constraints by solving strongly connected components individually. Strongly connected components (of the directed graph of constraints obtained from a perfect matching of variables and constraints) are @@ -139,9 +139,9 @@ def solve_strongly_connected_components( res_list = [] for scc, inputs in generate_strongly_connected_components( - constraints, - variables, - ): + constraints, + variables, + ): with TemporarySubsystemManager(to_fix=inputs): if len(scc.vars) == 1: results = calculate_variable_from_constraint( @@ -158,9 +158,8 @@ def solve_strongly_connected_components( raise RuntimeError( "An external solver is required if block has strongly\n" "connected components of size greater than one (is not " - "a DAG).\nGot an SCC with components: \n%s\n%s" - % (vars, cons) - ) + "a DAG).\nGot an SCC with components: \n%s\n%s" % (vars, cons) + ) results = solver.solve(scc, **solve_kwds) res_list.append(results) return res_list From 1b8a2b591bb7a343f79db2a631a9d9232a323ed9 Mon Sep 17 00:00:00 2001 From: robbybp Date: Mon, 30 Jan 2023 23:06:16 -0700 Subject: [PATCH 02/25] comments on planned implementation using a graph instead of a matrix --- pyomo/contrib/incidence_analysis/interface.py | 26 +++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py index 8ccd88c9c2a..5d144c63761 100644 --- a/pyomo/contrib/incidence_analysis/interface.py +++ b/pyomo/contrib/incidence_analysis/interface.py @@ -45,6 +45,7 @@ class IncidenceMatrixType(enum.Enum): NONE = 0 STRUCTURAL = 1 NUMERIC = 2 + # GRAPH = 3 def _check_unindexed(complist): @@ -172,6 +173,7 @@ def __init__( if model is None: self.cached = IncidenceMatrixType.NONE elif isinstance(model, PyomoNLP): + # TODO: If PyomoNLP is provided, cache a weighted incidence graph if not active: raise ValueError( "Cannot get the Jacobian of inactive constraints from the " @@ -216,9 +218,14 @@ def __init__( self.con_index_map = ComponentMap( (con, i) for i, con in enumerate(self.constraints) ) - self.incidence_matrix = get_structural_incidence_matrix( + #self.incidence_matrix = get_structural_incidence_matrix( + # self.variables, + # self.constraints, + #) + self.incidence_graph = get_incidence_graph( self.variables, self.constraints, + # TODO: include_fixed=include_fixed? ) else: raise TypeError( @@ -245,6 +252,7 @@ def _validate_input(self, variables, constraints): return variables, constraints def _extract_submatrix(self, variables, constraints): + # TODO: _extract_subgraph # Assumes variables and constraints are valid if self.cached is IncidenceMatrixType.NONE: return get_structural_incidence_matrix( @@ -253,13 +261,16 @@ def _extract_submatrix(self, variables, constraints): include_fixed=False, ) else: - N, M = len(variables), len(constraints) + N = len(variables) + M = len(constraints) old_new_var_indices = dict( (self.var_index_map[v], i) for i, v in enumerate(variables) ) old_new_con_indices = dict( (self.con_index_map[c], i) for i, c in enumerate(constraints) ) + # FIXME: This will fail if I don't have an incidence matrix + # cached. coo = self.incidence_matrix new_row = [] new_col = [] @@ -282,6 +293,7 @@ def maximum_matching(self, variables=None, constraints=None): variables, constraints = self._validate_input(variables, constraints) matrix = self._extract_submatrix(variables, constraints) + # TODO: This should just operate on the incidence graph matching = maximum_matching(matrix.tocoo()) # Matching maps row (constraint) indices to column (variable) indices @@ -296,6 +308,8 @@ def get_connected_components(self, variables=None, constraints=None): variables, constraints = self._validate_input(variables, constraints) matrix = self._extract_submatrix(variables, constraints) + # TODO: This could operate on a graph, but needs to know which nodes + # are variables/constraints row_blocks, col_blocks = get_independent_submatrices(matrix.tocoo()) con_blocks = [[constraints[i] for i in block] for block in row_blocks] var_blocks = [[variables[j] for j in block] for block in col_blocks] @@ -313,6 +327,9 @@ def block_triangularize(self, variables=None, constraints=None): variables, constraints = self._validate_input(variables, constraints) matrix = self._extract_submatrix(variables, constraints) + # TODO: block_triangularize does not quite make sense for incidence + # graphs, and it is unclear what the graph analog is. (Ordered SCCs + # of projection?) We may need to convert to a matrix here. row_block_map, col_block_map = block_triangularize(matrix.tocoo()) # Cache maps in case we want to get diagonal blocks quickly in the # future. @@ -328,6 +345,7 @@ def block_triangularize(self, variables=None, constraints=None): # Hopefully this does not get too confusing... return var_block_map, con_block_map + # TODO: deprecate this method and replace it with a better name. def get_diagonal_blocks(self, variables=None, constraints=None): """ Returns the diagonal blocks in a block triangularization of the @@ -344,6 +362,8 @@ def get_diagonal_blocks(self, variables=None, constraints=None): variables, constraints = self._validate_input(variables, constraints) matrix = self._extract_submatrix(variables, constraints) + # TODO: Again, this functionality does not really make sense for + # general bipartite graphs... if self.row_block_map is None or self.col_block_map is None: block_rows, block_cols = get_diagonal_blocks(matrix) else: @@ -371,6 +391,7 @@ def dulmage_mendelsohn(self, variables=None, constraints=None): variables, constraints = self._validate_input(variables, constraints) matrix = self._extract_submatrix(variables, constraints) + # TODO: Can just call Dulmage-Mendelsohn on the graph row_partition, col_partition = dulmage_mendelsohn(matrix.tocoo()) con_partition = RowPartition( *[[constraints[i] for i in subset] for subset in row_partition] @@ -414,6 +435,7 @@ def remove_nodes(self, nodes, constraints=None): to_exclude.update(constraints) vars_to_include = [v for v in self.variables if v not in to_exclude] cons_to_include = [c for c in self.constraints if c not in to_exclude] + # TODO: use _extract_subgraph here. incidence_matrix = self._extract_submatrix(vars_to_include, cons_to_include) # update attributes self.variables = vars_to_include From 78cf760f1248a978b8064f3364a4fb0c0ecaa727 Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Tue, 31 Jan 2023 15:01:24 -0700 Subject: [PATCH 03/25] update to cache an incidence graph instead of incidence matrix --- pyomo/contrib/incidence_analysis/interface.py | 113 ++++++++++++++++-- 1 file changed, 103 insertions(+), 10 deletions(-) diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py index 5d144c63761..8adecb4fcb2 100644 --- a/pyomo/contrib/incidence_analysis/interface.py +++ b/pyomo/contrib/incidence_analysis/interface.py @@ -58,8 +58,15 @@ def _check_unindexed(complist): def get_incidence_graph(variables, constraints, include_fixed=True): - """ - This function gets the incidence graph of Pyomo variables and constraints. + """Return the bipartite incidence graph of Pyomo variables and constraints. + + Each node in the returned graph is an integer. The convention is that, + for a graph with N variables and M constraints, nodes 0 through M-1 + correspond to constraints and nodes M through M+N-1 correspond to variables. + Nodes correspond to variables and constraints in the provided orders. + For consistency with NetworkX's "convention", constraint nodes are tagged + with `bipartite=0` while variable nodes are tagged with `bipartite=1`, + although these attributes are not used. Arguments: ---------- @@ -76,7 +83,8 @@ def get_incidence_graph(variables, constraints, include_fixed=True): """ _check_unindexed(variables + constraints) - N, M = len(variables), len(constraints) + N = len(variables) + M = len(constraints) graph = nx.Graph() graph.add_nodes_from(range(M), bipartite=0) graph.add_nodes_from(range(M, M + N), bipartite=1) @@ -88,6 +96,30 @@ def get_incidence_graph(variables, constraints, include_fixed=True): return graph +# TODO: Test this function +def extract_bipartite_subgraph(graph, nodes0, nodes1): + """ + """ + subgraph = nx.Graph() + sub_M = len(nodes0) + sub_N = len(nodes1) + subgraph.add_nodes_from(range(sub_M), bipartite=0) + subgraph.add_nodes_from(range(sub_M, sub_M + sub_N), bipartite=1) + + old_new_map = {} + for i, node in enumerate(nodes0 + nodes1): + if node in old_new_map: + raise RuntimeError("Node %s provided more than once.") + old_new_map[node] = i + + for (node1, node2) in graph.edges(): + if node1 in old_new_map and node2 in old_new_map: + new_node_1 = old_new_map[node1] + new_node_2 = old_new_map[node2] + subgraph.add_edge(new_node_1, new_node_2) + return subgraph + + def _generate_variables_in_constraints(constraints, include_fixed=False): known_vars = ComponentSet() for con in constraints: @@ -189,7 +221,10 @@ def __init__( nlp = model self.cached = IncidenceMatrixType.NUMERIC self.variables = nlp.get_pyomo_variables() - self.constraints = nlp.get_pyomo_constraints() + self.constraints = [ + con for con in nlp.get_pyomo_constraints() + if include_inequality or isinstance(con.expr, EqualityExpression) + ] self.var_index_map = ComponentMap( (var, idx) for idx, var in enumerate(self.variables) ) @@ -197,9 +232,18 @@ def __init__( (con, idx) for idx, con in enumerate(self.constraints) ) if include_inequality: - self.incidence_matrix = nlp.evaluate_jacobian() + incidence_matrix = nlp.evaluate_jacobian() else: - self.incidence_matrix = nlp.evaluate_jacobian_eq() + incidence_matrix = nlp.evaluate_jacobian_eq() + # + # TODO: If a PyomoNLP is provided, convert it to a weighted + # incidence graph. + # What advantage does a matrix have over a graph? It can be sent + # to linear algebra routines more directly. But this class doesn't + # do this, so I see no advantage to storing a Jacobian. + # + nxb = nx.algorithms.bipartite + self.incidence_graph = nxb.from_biadjacency_matrix(incidence_matrix) elif isinstance(model, Block): self.cached = IncidenceMatrixType.STRUCTURAL self.constraints = [ @@ -222,6 +266,11 @@ def __init__( # self.variables, # self.constraints, #) + + # + # If a model is provided, cache a graph instead of a matrix. + # This allows for faster lookups. + # self.incidence_graph = get_incidence_graph( self.variables, self.constraints, @@ -285,16 +334,58 @@ def _extract_submatrix(self, variables, constraints): shape=(M, N), ) + def _extract_subgraph(self, variables, constraints): + if self.cached is IncidenceMatrixType.NONE: + return get_incidence_graph( + # Does include_fixed matter here if I'm providing the variables? + variables, constraints, include_fixed=False + ) + else: + constraint_nodes = [self.con_index_map[con] for con in constraints] + M = len(self.con_index_map) + variable_nodes = [M + self.var_index_map[var] for var in variables] + subgraph = extract_bipartite_subgraph( + self.incidence_graph, constraint_nodes, variable_nodes + ) + return subgraph + + @property + def incidence_matrix(self): + if self.cached == IncidenceMatrixType.NONE: + return None + else: + M = len(self.constraints) + N = len(self.variables) + row = [] + col = [] + data = [] + # Here we assume that the incidence graph is bipartite with nodes + # 0 through M-1 forming one of the bipartite sets. + for i in range(M): + assert self.incidence_graph.nodes[i]["bipartite"] == 0 + for j in self.incidence_graph[i]: + assert self.incidence_graph.nodes[j]["bipartite"] == 1 + row.append(i) + col.append(j-M) + data.append(1.0) + return sp.sparse.coo_matrix( + (data, (row, col)), + shape=(M, N), + ) + def maximum_matching(self, variables=None, constraints=None): """ Returns a maximal matching between the constraints and variables, in terms of a map from constraints to variables. """ variables, constraints = self._validate_input(variables, constraints) - matrix = self._extract_submatrix(variables, constraints) + #matrix = self._extract_submatrix(variables, constraints) + graph = self._extract_subgraph(variables, constraints) # TODO: This should just operate on the incidence graph - matching = maximum_matching(matrix.tocoo()) + #matching = maximum_matching(matrix.tocoo()) + con_nodes = list(range(len(constraints))) + matching = maximum_matching(graph, top_nodes=con_nodes) # Matching maps row (constraint) indices to column (variable) indices return ComponentMap((constraints[i], variables[j]) for i, j in matching.items()) @@ -436,11 +527,13 @@ def remove_nodes(self, nodes, constraints=None): vars_to_include = [v for v in self.variables if v not in to_exclude] cons_to_include = [c for c in self.constraints if c not in to_exclude] # TODO: use _extract_subgraph here. - incidence_matrix = self._extract_submatrix(vars_to_include, cons_to_include) + #incidence_matrix = self._extract_submatrix(vars_to_include, cons_to_include) + incidence_graph = self._extract_subgraph(vars_to_include, cons_to_include) # update attributes self.variables = vars_to_include self.constraints = cons_to_include - self.incidence_matrix = incidence_matrix + #self.incidence_matrix = incidence_matrix + self.incidence_graph = incidence_graph self.var_index_map = ComponentMap( (var, i) for i, var in enumerate(self.variables) ) From 2e9f92f1aba6f67b656792ceb717bdbe3b1ddd48 Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Tue, 31 Jan 2023 15:01:48 -0700 Subject: [PATCH 04/25] update maximum_matching to accept a bipartite graph as well as matrix --- pyomo/contrib/incidence_analysis/matching.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/pyomo/contrib/incidence_analysis/matching.py b/pyomo/contrib/incidence_analysis/matching.py index d2a66e92d53..3a82db05f9f 100644 --- a/pyomo/contrib/incidence_analysis/matching.py +++ b/pyomo/contrib/incidence_analysis/matching.py @@ -12,7 +12,7 @@ from pyomo.common.dependencies import networkx as nx -def maximum_matching(matrix): +def maximum_matching(matrix_or_graph, top_nodes=None): """ Returns a maximum matching of the rows and columns of the matrix as a dict from row indices to column indices. @@ -21,8 +21,18 @@ def maximum_matching(matrix): nxc = nx.algorithms.components from_biadjacency_matrix = nxb.matrix.from_biadjacency_matrix - M, N = matrix.shape - bg = from_biadjacency_matrix(matrix) + if isinstance(matrix_or_graph, nx.Graph): + M = len(top_nodes) + N = len(matrix_or_graph.nodes) - M + bg = matrix_or_graph + else: + # Assume something SciPy-sparse compatible was provided. + if top_nodes is not None: + raise RuntimeError( + "top_nodes argument cannot be used if a matrix is provided" + ) + M, N = matrix_or_graph.shape + bg = from_biadjacency_matrix(matrix_or_graph) # Check assumptions regarding from_biadjacency_matrix function: for i in range(M): @@ -37,6 +47,8 @@ def maximum_matching(matrix): # This is fine, but we need to separate into connected components # for NetworkX to not complain. conn_comp = [bg.subgraph(c) for c in nxc.connected_components(bg)] + # NOTE: We could also provide top_nodes to maximum_matching to avoid + # ambiguity. matchings = [nxb.maximum_matching(c) for c in conn_comp] # If n0 < M, then n1 >= M. n0 is the row index, n1-M is the column index max_matching = {n0: n1 - M for m in matchings for n0, n1 in m.items() if n0 < M} From be47836f967450a68b6ba021635f6e10043cac8f Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Tue, 31 Jan 2023 16:32:20 -0700 Subject: [PATCH 05/25] add get_bipartite_incidence_graph to __init__ and start tests --- pyomo/contrib/incidence_analysis/__init__.py | 5 +- pyomo/contrib/incidence_analysis/interface.py | 10 +- .../tests/test_interface.py | 118 ++++++++++++++++++ 3 files changed, 130 insertions(+), 3 deletions(-) diff --git a/pyomo/contrib/incidence_analysis/__init__.py b/pyomo/contrib/incidence_analysis/__init__.py index 364046d51da..4dbfc66d25f 100644 --- a/pyomo/contrib/incidence_analysis/__init__.py +++ b/pyomo/contrib/incidence_analysis/__init__.py @@ -1,6 +1,9 @@ from .triangularize import block_triangularize from .matching import maximum_matching -from .interface import IncidenceGraphInterface +from .interface import ( + IncidenceGraphInterface, + get_bipartite_incidence_graph, +) from .util import ( generate_strongly_connected_components, solve_strongly_connected_components, diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py index 8adecb4fcb2..333776f9114 100644 --- a/pyomo/contrib/incidence_analysis/interface.py +++ b/pyomo/contrib/incidence_analysis/interface.py @@ -58,6 +58,12 @@ def _check_unindexed(complist): def get_incidence_graph(variables, constraints, include_fixed=True): + return get_bipartite_incidence_graph( + variables, constraints, include_fixed=include_fixed + ) + + +def get_bipartite_incidence_graph(variables, constraints, include_fixed=True): """Return the bipartite incidence graph of Pyomo variables and constraints. Each node in the returned graph is an integer. The convention is that, @@ -271,7 +277,7 @@ def __init__( # If a model is provided, cache a graph instead of a matrix. # This allows for faster lookups. # - self.incidence_graph = get_incidence_graph( + self.incidence_graph = get_bipartite_incidence_graph( self.variables, self.constraints, # TODO: include_fixed=include_fixed? @@ -336,7 +342,7 @@ def _extract_submatrix(self, variables, constraints): def _extract_subgraph(self, variables, constraints): if self.cached is IncidenceMatrixType.NONE: - return get_incidence_graph( + return get_bipartite_incidence_graph( # Does include_fixed matter here if I'm providing the variables? variables, constraints, include_fixed=False ) diff --git a/pyomo/contrib/incidence_analysis/tests/test_interface.py b/pyomo/contrib/incidence_analysis/tests/test_interface.py index c4f4a1e0e8a..85f36e54c1b 100644 --- a/pyomo/contrib/incidence_analysis/tests/test_interface.py +++ b/pyomo/contrib/incidence_analysis/tests/test_interface.py @@ -18,6 +18,7 @@ get_structural_incidence_matrix, get_numeric_incidence_matrix, get_incidence_graph, + get_bipartite_incidence_graph, ) from pyomo.contrib.incidence_analysis.matching import maximum_matching from pyomo.contrib.incidence_analysis.triangularize import block_triangularize @@ -32,6 +33,7 @@ if scipy_available: from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP if networkx_available: + import networkx as nx from networkx.algorithms.bipartite.matrix import from_biadjacency_matrix from pyomo.contrib.pynumero.asl import AmplInterface @@ -1238,5 +1240,121 @@ def test_include_inequality_nlp(self): self.assertEqual(igraph.incidence_matrix.shape, (12, 8)) +@unittest.skipUnless(networkx_available, "networkx is not available.") +@unittest.skipUnless(scipy_available, "scipy is not available.") +class TestGetIncidenceGraph(unittest.TestCase): + + def make_test_model(self): + m = pyo.ConcreteModel() + m.I = pyo.Set(initialize=[1, 2, 3, 4]) + m.v = pyo.Var(m.I, bounds=(0, None)) + m.eq1 = pyo.Constraint(expr=m.v[1]**2 + m.v[2]**2 == 1.0) + m.eq2 = pyo.Constraint(expr=m.v[1] + 2.0 == m.v[3]) + m.ineq1 = pyo.Constraint(expr=m.v[2] - m.v[3]**0.5 + m.v[4]**2 <= 1.0) + m.ineq2 = pyo.Constraint(expr=m.v[2]*m.v[4] >= 1.0) + m.ineq3 = pyo.Constraint(expr=m.v[1] >= m.v[4]**4) + m.obj = pyo.Objective(expr=-m.v[1] - m.v[2] + m.v[3]**2 + m.v[4]**2) + return m + + def test_bipartite_incidence_graph(self): + m = self.make_test_model() + constraints = [m.eq1, m.eq2, m.ineq1, m.ineq2, m.ineq3] + variables = list(m.v.values()) + graph = get_bipartite_incidence_graph(variables, constraints) + + # Nodes: + # 0: m.eq1 + # 1: m.eq2 + # 2: m.ineq1 + # 3: m.ineq2 + # 4: m.ineq3 + # 5: m.v[1] + # 6: m.v[2] + # 7: m.v[3] + # 8: m.v[4] + + # Assert some basic structure + self.assertEqual(len(graph.nodes), 9) + self.assertEqual(len(graph.edges), 11) + self.assertTrue(nx.algorithms.bipartite.is_bipartite(graph)) + + # Assert that the "adjacency list" is what we expect + self.assertEqual(set(graph[0]), {5, 6}) + self.assertEqual(set(graph[1]), {5, 7}) + self.assertEqual(set(graph[2]), {6, 7, 8}) + self.assertEqual(set(graph[3]), {6, 8}) + self.assertEqual(set(graph[4]), {5, 8}) + self.assertEqual(set(graph[5]), {0, 1, 4}) + self.assertEqual(set(graph[6]), {0, 2, 3}) + self.assertEqual(set(graph[7]), {1, 2}) + self.assertEqual(set(graph[8]), {2, 3, 4}) + + def test_unused_var(self): + m = self.make_test_model() + constraints = [m.eq1, m.eq2] + variables = list(m.v.values()) + graph = get_bipartite_incidence_graph(variables, constraints) + + # Nodes: + # 0: m.eq1 + # 1: m.eq2 + # 2: m.v[1] + # 3: m.v[2] + # 4: m.v[3] + # 5: m.v[4] + + self.assertEqual(len(graph.nodes), 6) + self.assertEqual(len(graph.edges), 4) + self.assertTrue(nx.algorithms.bipartite.is_bipartite(graph)) + + # Assert that the "adjacency list" is what we expect + self.assertEqual(set(graph[0]), {2, 3}) + self.assertEqual(set(graph[1]), {2, 4}) + self.assertEqual(set(graph[2]), {0, 1}) + self.assertEqual(set(graph[3]), {0}) + self.assertEqual(set(graph[4]), {1}) + self.assertEqual(set(graph[5]), set()) + + def test_fixed_vars(self): + m = self.make_test_model() + constraints = [m.eq1, m.eq2, m.ineq1, m.ineq2, m.ineq3] + variables = list(m.v.values()) + m.v[1].fix() + m.v[4].fix() + + # Slightly odd situation where we provide fixed variables, but + # then tell the graph to not include them. Nodes will be created + # for these vars, but they will not have any edges. + graph = get_bipartite_incidence_graph( + variables, constraints, include_fixed=False + ) + + # Nodes: + # 0: m.eq1 + # 1: m.eq2 + # 2: m.ineq1 + # 3: m.ineq2 + # 4: m.ineq3 + # 5: m.v[1] + # 6: m.v[2] + # 7: m.v[3] + # 8: m.v[4] + + self.assertEqual(len(graph.nodes), 9) + self.assertEqual(len(graph.edges), 5) + self.assertTrue(nx.algorithms.bipartite.is_bipartite(graph)) + + # Assert that the "adjacency list" is what we expect + self.assertEqual(set(graph[0]), {6}) + self.assertEqual(set(graph[1]), {7}) + self.assertEqual(set(graph[2]), {6, 7}) + self.assertEqual(set(graph[3]), {6}) + self.assertEqual(set(graph[4]), set()) + self.assertEqual(set(graph[5]), set()) + self.assertEqual(set(graph[6]), {0, 2, 3}) + self.assertEqual(set(graph[7]), {1, 2}) + self.assertEqual(set(graph[8]), set()) + + if __name__ == "__main__": unittest.main() From a7ab473e58a42bbca0c89061936929c4252d2289 Mon Sep 17 00:00:00 2001 From: robbybp Date: Tue, 31 Jan 2023 18:13:39 -0700 Subject: [PATCH 06/25] add tests and docstring for extract_bipartite_subgraph --- pyomo/contrib/incidence_analysis/interface.py | 36 +++++++++++++- .../tests/test_interface.py | 47 +++++++++++++++++++ 2 files changed, 81 insertions(+), 2 deletions(-) diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py index 333776f9114..c93d122ab87 100644 --- a/pyomo/contrib/incidence_analysis/interface.py +++ b/pyomo/contrib/incidence_analysis/interface.py @@ -102,9 +102,33 @@ def get_bipartite_incidence_graph(variables, constraints, include_fixed=True): return graph -# TODO: Test this function def extract_bipartite_subgraph(graph, nodes0, nodes1): - """ + """Return the bipartite subgraph of a graph. + + Two lists of nodes to project onto must be provided. These will correspond + to the "bipartite sets" in the subgraph. If the two sets provided have + M and N nodes, the subgraph will have nodes 0 through M+N, with the first + M corresponding to the first set provided and the last M corresponding + to the second set. + + Parameters + ---------- + graph: NetworkX Graph + The graph from which a subgraph is extracted + nodes0: list + A list of nodes in the original graph that will form the first + bipartite set of the projected graph (and have ``bipartite=0`` + nodes1: list + A list of nodes in the original graph that will form the second + bipartite set of the projected graph (and have ``bipartite=1`` + + Returns + ------- + subgraph: NetworkX Graph + Graph containing integer nodes corresponding to positions in the + provided lists, with edges where corresponding nodes are adjacent + in the original graph. + """ subgraph = nx.Graph() sub_M = len(nodes0) @@ -122,6 +146,14 @@ def extract_bipartite_subgraph(graph, nodes0, nodes1): if node1 in old_new_map and node2 in old_new_map: new_node_1 = old_new_map[node1] new_node_2 = old_new_map[node2] + if ( + subgraph.nodes[new_node_1]["bipartite"] + == subgraph.nodes[new_node_2]["bipartite"] + ): + raise RuntimeError( + "Subgraph is not bipartite. Found an edge between nodes" + " %s and %s (in the original graph)." % (node1, node2) + ) subgraph.add_edge(new_node_1, new_node_2) return subgraph diff --git a/pyomo/contrib/incidence_analysis/tests/test_interface.py b/pyomo/contrib/incidence_analysis/tests/test_interface.py index 85f36e54c1b..9d7d36aee73 100644 --- a/pyomo/contrib/incidence_analysis/tests/test_interface.py +++ b/pyomo/contrib/incidence_analysis/tests/test_interface.py @@ -19,6 +19,7 @@ get_numeric_incidence_matrix, get_incidence_graph, get_bipartite_incidence_graph, + extract_bipartite_subgraph, ) from pyomo.contrib.incidence_analysis.matching import maximum_matching from pyomo.contrib.incidence_analysis.triangularize import block_triangularize @@ -1355,6 +1356,52 @@ def test_fixed_vars(self): self.assertEqual(set(graph[7]), {1, 2}) self.assertEqual(set(graph[8]), set()) + def test_extract_subgraph(self): + m = self.make_test_model() + constraints = [m.eq1, m.eq2, m.ineq1, m.ineq2, m.ineq3] + variables = list(m.v.values()) + graph = get_bipartite_incidence_graph(variables, constraints) + + sg_cons = [0, 2] + sg_vars = [i + len(constraints) for i in [2, 0, 3]] + + subgraph = extract_bipartite_subgraph(graph, sg_cons, sg_vars) + + # Subgraph nodes: + # 0: m.eq1 + # 1: m.ineq1 + # 2: m.v[3] + # 3: m.v[1] + # 4: m.v[4] + + self.assertEqual(len(subgraph.nodes), 5) + self.assertEqual(len(subgraph.edges), 3) + self.assertTrue(nx.algorithms.bipartite.is_bipartite(subgraph)) + + self.assertEqual(set(subgraph[0]), {3}) + self.assertEqual(set(subgraph[1]), {2, 4}) + self.assertEqual(set(subgraph[2]), {1}) + self.assertEqual(set(subgraph[3]), {0}) + self.assertEqual(set(subgraph[4]), {1}) + + def test_extract_exceptions(self): + m = self.make_test_model() + constraints = [m.eq1, m.eq2, m.ineq1, m.ineq2, m.ineq3] + variables = list(m.v.values()) + graph = get_bipartite_incidence_graph(variables, constraints) + + sg_cons = [0, 2, 5] + sg_vars = [i + len(constraints) for i in [2, 3]] + msg = "Subgraph is not bipartite" + with self.assertRaisesRegex(RuntimeError, msg): + subgraph = extract_bipartite_subgraph(graph, sg_cons, sg_vars) + + sg_cons = [0, 2, 5] + sg_vars = [i + len(constraints) for i in [2, 0, 3]] + msg = "provided more than once" + with self.assertRaisesRegex(RuntimeError, msg): + subgraph = extract_bipartite_subgraph(graph, sg_cons, sg_vars) + if __name__ == "__main__": unittest.main() From 8a62a96568b0394714463782b5c8aa46a45c977b Mon Sep 17 00:00:00 2001 From: robbybp Date: Tue, 31 Jan 2023 19:02:27 -0700 Subject: [PATCH 07/25] get_adjacent_to method and tests --- pyomo/contrib/incidence_analysis/interface.py | 41 ++++++++++++++++- .../tests/test_interface.py | 46 +++++++++++++++++++ 2 files changed, 86 insertions(+), 1 deletion(-) diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py index c93d122ab87..3cfe6e50047 100644 --- a/pyomo/contrib/incidence_analysis/interface.py +++ b/pyomo/contrib/incidence_analysis/interface.py @@ -339,7 +339,6 @@ def _validate_input(self, variables, constraints): return variables, constraints def _extract_submatrix(self, variables, constraints): - # TODO: _extract_subgraph # Assumes variables and constraints are valid if self.cached is IncidenceMatrixType.NONE: return get_structural_incidence_matrix( @@ -411,6 +410,46 @@ def incidence_matrix(self): shape=(M, N), ) + def get_adjacent_to(self, component): + """Return a list of components adjacent to the provided component + in the cached bipartite incidence graph of variables and constraints + + Parameters + ---------- + component: ComponentData + The variable or constraint data object whose adjacent components + are returned + + Returns + ------- + list of ComponentData + List of constraint or variable data objects adjacent to the + provided component + + """ + if self.cached == IncidenceMatrixType.NONE: + raise RuntimeError( + "Cannot get components adjacent to %s if an incidence graph" + " is not cached." % component + ) + _check_unindexed([component]) + M = len(self.constraints) + N = len(self.variables) + if component in self.var_index_map: + vnode = M + self.var_index_map[component] + adj = self.incidence_graph[vnode] + adj_comps = [self.constraints[i] for i in adj] + elif component in self.con_index_map: + cnode = self.con_index_map[component] + adj = self.incidence_graph[cnode] + adj_comps = [self.variables[j-M] for j in adj] + else: + raise RuntimeError( + "Cannot find component %s in the cached incidence graph." + % component + ) + return adj_comps + def maximum_matching(self, variables=None, constraints=None): """ Returns a maximal matching between the constraints and variables, diff --git a/pyomo/contrib/incidence_analysis/tests/test_interface.py b/pyomo/contrib/incidence_analysis/tests/test_interface.py index 9d7d36aee73..08bfffcf493 100644 --- a/pyomo/contrib/incidence_analysis/tests/test_interface.py +++ b/pyomo/contrib/incidence_analysis/tests/test_interface.py @@ -1403,5 +1403,51 @@ def test_extract_exceptions(self): subgraph = extract_bipartite_subgraph(graph, sg_cons, sg_vars) +@unittest.skipUnless(networkx_available, "networkx is not available.") +@unittest.skipUnless(scipy_available, "scipy is not available.") +class TestGetAdjacent(unittest.TestCase): + + def test_get_adjacent_to_var(self): + m = make_degenerate_solid_phase_model() + igraph = IncidenceGraphInterface(m) + adj_cons = igraph.get_adjacent_to(m.rho) + self.assertEqual( + ComponentSet(adj_cons), + ComponentSet([ + m.holdup_eqn[1], + m.holdup_eqn[2], + m.holdup_eqn[3], + m.density_eqn, + ]), + ) + + def test_get_adjacent_to_con(self): + m = make_degenerate_solid_phase_model() + igraph = IncidenceGraphInterface(m) + adj_vars = igraph.get_adjacent_to(m.density_eqn) + self.assertEqual( + ComponentSet(adj_vars), + ComponentSet([ + m.x[1], + m.x[2], + m.x[3], + m.rho, + ]), + ) + + def test_get_adjacent_exceptions(self): + m = make_degenerate_solid_phase_model() + igraph = IncidenceGraphInterface() + msg = "Cannot get components adjacent to" + with self.assertRaisesRegex(RuntimeError, msg): + adj_vars = igraph.get_adjacent_to(m.density_eqn) + + m.x[1].fix() + igraph = IncidenceGraphInterface(m, include_fixed=False) + msg = "Cannot find component" + with self.assertRaisesRegex(RuntimeError, msg): + adj_cons = igraph.get_adjacent_to(m.x[1]) + + if __name__ == "__main__": unittest.main() From f9c9f04c673a81a78f9a8d1f130b9d701a521b8e Mon Sep 17 00:00:00 2001 From: robbybp Date: Tue, 31 Jan 2023 20:02:00 -0700 Subject: [PATCH 08/25] update methods to use graph instead of matrix --- pyomo/contrib/incidence_analysis/interface.py | 44 +++++++++---------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py index 3cfe6e50047..bef138b9cb9 100644 --- a/pyomo/contrib/incidence_analysis/interface.py +++ b/pyomo/contrib/incidence_analysis/interface.py @@ -456,15 +456,10 @@ def maximum_matching(self, variables=None, constraints=None): in terms of a map from constraints to variables. """ variables, constraints = self._validate_input(variables, constraints) - #matrix = self._extract_submatrix(variables, constraints) graph = self._extract_subgraph(variables, constraints) - - # TODO: This should just operate on the incidence graph - #matching = maximum_matching(matrix.tocoo()) con_nodes = list(range(len(constraints))) matching = maximum_matching(graph, top_nodes=con_nodes) - # Matching maps row (constraint) indices to column (variable) indices - + # Matching maps constraint indices to variable indices return ComponentMap((constraints[i], variables[j]) for i, j in matching.items()) def get_connected_components(self, variables=None, constraints=None): @@ -474,15 +469,21 @@ def get_connected_components(self, variables=None, constraints=None): and constraints. """ variables, constraints = self._validate_input(variables, constraints) - matrix = self._extract_submatrix(variables, constraints) + graph = self._extract_subgraph(variables, constraints) + nxc = nx.algorithms.components + M = len(constraints) + N = len(variables) + connected_components = list(nxc.connected_components(graph)) + + con_blocks = [ + sorted([i for i in comp if i < M]) for comp in connected_components + ] + con_blocks = [[constraints[i] for i in block] for block in con_blocks] + var_blocks = [ + sorted([j for j in comp if j >= M]) for comp in connected_components + ] + var_blocks = [[variables[i-M] for i in block] for block in var_blocks] - # TODO: This could operate on a graph, but needs to know which nodes - # are variables/constraints - row_blocks, col_blocks = get_independent_submatrices(matrix.tocoo()) - con_blocks = [[constraints[i] for i in block] for block in row_blocks] - var_blocks = [[variables[j] for j in block] for block in col_blocks] - # Switch the order of the partitions here to match the method call. - # Hopefully this does not get too confusing... return var_blocks, con_blocks def block_triangularize(self, variables=None, constraints=None): @@ -557,15 +558,17 @@ def dulmage_mendelsohn(self, variables=None, constraints=None): """ variables, constraints = self._validate_input(variables, constraints) - matrix = self._extract_submatrix(variables, constraints) - - # TODO: Can just call Dulmage-Mendelsohn on the graph - row_partition, col_partition = dulmage_mendelsohn(matrix.tocoo()) + graph = self._extract_subgraph(variables, constraints) + M = len(constraints) + top_nodes = list(range(M)) + row_partition, col_partition = dulmage_mendelsohn( + graph, top_nodes=top_nodes + ) con_partition = RowPartition( *[[constraints[i] for i in subset] for subset in row_partition] ) var_partition = ColPartition( - *[[variables[i] for i in subset] for subset in col_partition] + *[[variables[i-M] for i in subset] for subset in col_partition] ) # Switch the order of the maps here to match the method call. # Hopefully this does not get too confusing... @@ -603,13 +606,10 @@ def remove_nodes(self, nodes, constraints=None): to_exclude.update(constraints) vars_to_include = [v for v in self.variables if v not in to_exclude] cons_to_include = [c for c in self.constraints if c not in to_exclude] - # TODO: use _extract_subgraph here. - #incidence_matrix = self._extract_submatrix(vars_to_include, cons_to_include) incidence_graph = self._extract_subgraph(vars_to_include, cons_to_include) # update attributes self.variables = vars_to_include self.constraints = cons_to_include - #self.incidence_matrix = incidence_matrix self.incidence_graph = incidence_graph self.var_index_map = ComponentMap( (var, i) for i, var in enumerate(self.variables) From 73033eb426c901a47fb41b8d12d04a02cd012f7d Mon Sep 17 00:00:00 2001 From: robbybp Date: Tue, 31 Jan 2023 21:06:12 -0700 Subject: [PATCH 09/25] update maximum_matching to accept a bipartite graph with fewer assumptions --- pyomo/contrib/incidence_analysis/interface.py | 8 +- pyomo/contrib/incidence_analysis/matching.py | 75 +++++++++++++------ 2 files changed, 60 insertions(+), 23 deletions(-) diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py index bef138b9cb9..5dc16d507f9 100644 --- a/pyomo/contrib/incidence_analysis/interface.py +++ b/pyomo/contrib/incidence_analysis/interface.py @@ -459,8 +459,12 @@ def maximum_matching(self, variables=None, constraints=None): graph = self._extract_subgraph(variables, constraints) con_nodes = list(range(len(constraints))) matching = maximum_matching(graph, top_nodes=con_nodes) - # Matching maps constraint indices to variable indices - return ComponentMap((constraints[i], variables[j]) for i, j in matching.items()) + # Matching maps constraint nodes to variable nodes. Here we need to + # know the convention according to which the graph was constructed. + M = len(constraints) + return ComponentMap( + (constraints[i], variables[j-M]) for i, j in matching.items() + ) def get_connected_components(self, variables=None, constraints=None): """ diff --git a/pyomo/contrib/incidence_analysis/matching.py b/pyomo/contrib/incidence_analysis/matching.py index 3a82db05f9f..e9ba4e3fb00 100644 --- a/pyomo/contrib/incidence_analysis/matching.py +++ b/pyomo/contrib/incidence_analysis/matching.py @@ -13,43 +13,76 @@ def maximum_matching(matrix_or_graph, top_nodes=None): - """ - Returns a maximum matching of the rows and columns of the - matrix as a dict from row indices to column indices. + """Return a maximum cardinality matching of the provided matrix or + bipartite graph + + If a matrix is provided, the matching is returned as a map from row + indices to column indices. If a bipartite graph is provided, a list of + "top nodes" must be provided as well. These correspond to one of the + "bipartite sets". The matching is then returned as a map from "top nodes" + to the other set of nodes. + + Parameters + ---------- + matrix_or_graph: SciPy sparse matrix or NetworkX Graph + The matrix or graph whose maximum matching will be computed + top_nodes: list + Integer nodes representing a bipartite set in a graph. Must be provided + if and only if a NetworkX Graph is provided. + + Returns + ------- + max_matching: dict + Dict mapping from integer nodes in the first bipartite set (row + indices) to nodes in the second (column indices). + """ nxb = nx.algorithms.bipartite nxc = nx.algorithms.components from_biadjacency_matrix = nxb.matrix.from_biadjacency_matrix if isinstance(matrix_or_graph, nx.Graph): + graph_provided = True + if top_nodes is None: + raise RuntimeError( + "top_nodes argument must be set if a graph is provided." + ) M = len(top_nodes) N = len(matrix_or_graph.nodes) - M bg = matrix_or_graph + if not nxb.is_bipartite(bg): + raise RuntimeError("Provided graph is not bipartite.") else: + graph_provided = False # Assume something SciPy-sparse compatible was provided. if top_nodes is not None: raise RuntimeError( "top_nodes argument cannot be used if a matrix is provided" ) M, N = matrix_or_graph.shape + top_nodes = list(range(M)) bg = from_biadjacency_matrix(matrix_or_graph) - # Check assumptions regarding from_biadjacency_matrix function: - for i in range(M): - # First M nodes in graph correspond to rows - assert bg.nodes[i]["bipartite"] == 0 - - for j in range(M, M + N): - # Last N nodes in graph correspond to columns - assert bg.nodes[j]["bipartite"] == 1 - - # If the matrix is block diagonal, the graph will be disconnected. - # This is fine, but we need to separate into connected components - # for NetworkX to not complain. - conn_comp = [bg.subgraph(c) for c in nxc.connected_components(bg)] - # NOTE: We could also provide top_nodes to maximum_matching to avoid - # ambiguity. - matchings = [nxb.maximum_matching(c) for c in conn_comp] - # If n0 < M, then n1 >= M. n0 is the row index, n1-M is the column index - max_matching = {n0: n1 - M for m in matchings for n0, n1 in m.items() if n0 < M} + # Check assumptions regarding from_biadjacency_matrix function: + for i in range(M): + # First M nodes in graph correspond to rows + assert bg.nodes[i]["bipartite"] == 0 + + for j in range(M, M + N): + # Last N nodes in graph correspond to columns + assert bg.nodes[j]["bipartite"] == 1 + + matching = nxb.maximum_matching(bg, top_nodes=top_nodes) + if graph_provided: + top_node_set = set(top_nodes) + # If a graph was provided, we return a mapping from "top nodes" + # to their matched "bottom nodes" + max_matching = { + n0: n1 for n0, n1 in matching.items() if n0 in top_node_set + } + else: + # If a matrix was provided, we return a mapping from row indices + # to column indices. + max_matching = {n0: n1-M for n0, n1 in matching.items() if n0 < M} + return max_matching From 47933451959c3aa99e265f7d81a2042d9d6bce77 Mon Sep 17 00:00:00 2001 From: robbybp Date: Tue, 31 Jan 2023 22:18:28 -0700 Subject: [PATCH 10/25] allow block_triangularize to accept a bipartite graph --- .../incidence_analysis/triangularize.py | 53 +++++++++++++++---- 1 file changed, 42 insertions(+), 11 deletions(-) diff --git a/pyomo/contrib/incidence_analysis/triangularize.py b/pyomo/contrib/incidence_analysis/triangularize.py index 356c6e0e9aa..f1a2be737fb 100644 --- a/pyomo/contrib/incidence_analysis/triangularize.py +++ b/pyomo/contrib/incidence_analysis/triangularize.py @@ -13,7 +13,10 @@ from pyomo.common.dependencies import networkx as nx -def block_triangularize(matrix, matching=None): +# TODO: get_scc_of_projection(bipartite_graph, top_nodes, matching) + + +def block_triangularize(matrix, matching=None, top_nodes=None): """ Computes the necessary information to permute a matrix to block-lower triangular form, i.e. a partition of rows and columns into an ordered @@ -37,16 +40,37 @@ def block_triangularize(matrix, matching=None): nxd = nx.algorithms.dag from_biadjacency_matrix = nxb.matrix.from_biadjacency_matrix - M, N = matrix.shape + if isinstance(matrix, nx.Graph): + graph_provided = True + if top_nodes is None: + raise RuntimeError( + "top_nodes argument must be set if a graph is provided." + ) + bg = matrix + M = len(top_nodes) + N = len(bg.nodes) - M + if not nxb.is_bipartite(bg): + raise RuntimeError("Provided graph is not bipartite.") + else: + graph_provided = False + M, N = matrix.shape + bg = from_biadjacency_matrix(matrix) + if M != N: raise ValueError( "block_triangularize does not currently " - "support non-square matrices. Got matrix with shape %s." % (matrix.shape,) + "support non-square matrices. Got matrix with shape %s." + % ((M, N),) ) - bg = from_biadjacency_matrix(matrix) if matching is None: - matching = maximum_matching(matrix) + matching = maximum_matching(matrix, top_nodes=top_nodes) + + if not graph_provided: + # If we provided a matrix, the matching maps row indices to column + # indices. Update the matching to map to nodes in the graph derived + # from the matrix. + matching = {r: c+M for r, c in matching.items()} len_matching = len(matching) if len_matching != M: @@ -58,10 +82,12 @@ def block_triangularize(matrix, matching=None): # Construct directed graph of rows dg = nx.DiGraph() - dg.add_nodes_from(range(M)) + if graph_provided: + dg.add_nodes_from(top_nodes) + else: + dg.add_nodes_from(range(M)) for n in dg.nodes: - col_idx = matching[n] - col_node = col_idx + M + col_node = matching[n] # For all rows that share this column for neighbor in bg[col_node]: if neighbor != n: @@ -93,11 +119,16 @@ def block_triangularize(matrix, matching=None): # ^ This maps row indices to the blocks they belong to. # Invert the matching to map row indices to column indices - col_row_map = {c: r for r, c in matching.items()} + if graph_provided: + col_row_map = {c: r for r, c in matching.items()} + col_block_map = { + c: row_block_map[col_row_map[c]] for c in matching.values() + } + else: + col_row_map = {c-M: r for r, c in matching.items()} + col_block_map = {c: row_block_map[col_row_map[c]] for c in range(N)} assert len(col_row_map) == M - col_block_map = {c: row_block_map[col_row_map[c]] for c in range(N)} - return row_block_map, col_block_map From e6d9c82eb98c01162898cb49574e7abb52603efa Mon Sep 17 00:00:00 2001 From: robbybp Date: Tue, 31 Jan 2023 22:20:57 -0700 Subject: [PATCH 11/25] update block_triangularize to use the cached graph instead of a matrix, and remove some comments --- pyomo/contrib/incidence_analysis/interface.py | 57 +++++++------------ 1 file changed, 22 insertions(+), 35 deletions(-) diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py index 5dc16d507f9..ab9b71b7b53 100644 --- a/pyomo/contrib/incidence_analysis/interface.py +++ b/pyomo/contrib/incidence_analysis/interface.py @@ -236,14 +236,13 @@ def __init__( ): """ """ # If the user gives us a model or an NLP, we assume they want us - # to cache the incidence matrix for fast analysis of submatrices - # later on. + # to cache the incidence graph for fast analysis later on. # WARNING: This cache will become invalid if the user alters their - # model. + # model. if model is None: self.cached = IncidenceMatrixType.NONE + self.incidence_graph = None elif isinstance(model, PyomoNLP): - # TODO: If PyomoNLP is provided, cache a weighted incidence graph if not active: raise ValueError( "Cannot get the Jacobian of inactive constraints from the " @@ -273,20 +272,14 @@ def __init__( incidence_matrix = nlp.evaluate_jacobian() else: incidence_matrix = nlp.evaluate_jacobian_eq() - # - # TODO: If a PyomoNLP is provided, convert it to a weighted - # incidence graph. - # What advantage does a matrix have over a graph? It can be sent - # to linear algebra routines more directly. But this class doesn't - # do this, so I see no advantage to storing a Jacobian. - # nxb = nx.algorithms.bipartite self.incidence_graph = nxb.from_biadjacency_matrix(incidence_matrix) elif isinstance(model, Block): self.cached = IncidenceMatrixType.STRUCTURAL self.constraints = [ - con - for con in model.component_data_objects(Constraint, active=active) + con for con in model.component_data_objects( + Constraint, active=active + ) if include_inequality or isinstance(con.expr, EqualityExpression) ] self.variables = list( @@ -300,15 +293,6 @@ def __init__( self.con_index_map = ComponentMap( (con, i) for i, con in enumerate(self.constraints) ) - #self.incidence_matrix = get_structural_incidence_matrix( - # self.variables, - # self.constraints, - #) - - # - # If a model is provided, cache a graph instead of a matrix. - # This allows for faster lookups. - # self.incidence_graph = get_bipartite_incidence_graph( self.variables, self.constraints, @@ -325,12 +309,12 @@ def __init__( def _validate_input(self, variables, constraints): if variables is None: - if self.cached is IncidenceMatrixType.NONE: + if self.incidence_graph is None: raise ValueError("Neither variables nor a model have been provided.") else: variables = self.variables if constraints is None: - if self.cached is IncidenceMatrixType.NONE: + if self.incidence_graph is None: raise ValueError("Neither constraints nor a model have been provided.") else: constraints = self.constraints @@ -340,7 +324,7 @@ def _validate_input(self, variables, constraints): def _extract_submatrix(self, variables, constraints): # Assumes variables and constraints are valid - if self.cached is IncidenceMatrixType.NONE: + if self.incidence_graph is None: return get_structural_incidence_matrix( variables, constraints, @@ -372,7 +356,7 @@ def _extract_submatrix(self, variables, constraints): ) def _extract_subgraph(self, variables, constraints): - if self.cached is IncidenceMatrixType.NONE: + if self.incidence_graph is None: return get_bipartite_incidence_graph( # Does include_fixed matter here if I'm providing the variables? variables, constraints, include_fixed=False @@ -388,7 +372,7 @@ def _extract_subgraph(self, variables, constraints): @property def incidence_matrix(self): - if self.cached == IncidenceMatrixType.NONE: + if self.incidence_graph is None: return None else: M = len(self.constraints) @@ -427,7 +411,7 @@ def get_adjacent_to(self, component): provided component """ - if self.cached == IncidenceMatrixType.NONE: + if self.incidence_graph is None: raise RuntimeError( "Cannot get components adjacent to %s if an incidence graph" " is not cached." % component @@ -498,14 +482,17 @@ def block_triangularize(self, variables=None, constraints=None): of the incidence matrix. """ variables, constraints = self._validate_input(variables, constraints) - matrix = self._extract_submatrix(variables, constraints) + graph = self._extract_subgraph(variables, constraints) - # TODO: block_triangularize does not quite make sense for incidence - # graphs, and it is unclear what the graph analog is. (Ordered SCCs - # of projection?) We may need to convert to a matrix here. - row_block_map, col_block_map = block_triangularize(matrix.tocoo()) + M = len(constraints) + top_nodes = list(range(M)) + cnode_block_map, vnode_block_map = block_triangularize( + graph, top_nodes=top_nodes + ) # Cache maps in case we want to get diagonal blocks quickly in the # future. + row_block_map = cnode_block_map + col_block_map = {j-M: idx for j, idx in vnode_block_map.items()} self.row_block_map = row_block_map self.col_block_map = col_block_map con_block_map = ComponentMap( @@ -518,7 +505,7 @@ def block_triangularize(self, variables=None, constraints=None): # Hopefully this does not get too confusing... return var_block_map, con_block_map - # TODO: deprecate this method and replace it with a better name. + # TODO: Update this method with a new name and deprecate old name. def get_diagonal_blocks(self, variables=None, constraints=None): """ Returns the diagonal blocks in a block triangularization of the @@ -601,7 +588,7 @@ def remove_nodes(self, nodes, constraints=None): """ if constraints is None: constraints = [] - if self.cached is IncidenceMatrixType.NONE: + if self.incidence_graph is None: raise RuntimeError( "Attempting to remove variables and constraints from cached " "incidence matrix,\nbut no incidence matrix has been cached." From 847d5b844c5e58e39f84dc137a6036de9c1b2e0c Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Wed, 1 Feb 2023 12:12:21 -0700 Subject: [PATCH 12/25] refactor block_triangularize to use an intermediate get_scc_of_projection function --- pyomo/contrib/incidence_analysis/interface.py | 10 +- .../incidence_analysis/triangularize.py | 177 +++++++++--------- 2 files changed, 98 insertions(+), 89 deletions(-) diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py index ab9b71b7b53..757b7955811 100644 --- a/pyomo/contrib/incidence_analysis/interface.py +++ b/pyomo/contrib/incidence_analysis/interface.py @@ -26,6 +26,7 @@ get_independent_submatrices, ) from pyomo.contrib.incidence_analysis.triangularize import ( + get_scc_of_projection, block_triangularize, get_diagonal_blocks, get_blocks_from_maps, @@ -485,10 +486,13 @@ def block_triangularize(self, variables=None, constraints=None): graph = self._extract_subgraph(variables, constraints) M = len(constraints) - top_nodes = list(range(M)) - cnode_block_map, vnode_block_map = block_triangularize( - graph, top_nodes=top_nodes + con_nodes = list(range(M)) + cnode_block_map, vnode_block_map = get_scc_of_projection( + graph, con_nodes ) + #cnode_block_map, vnode_block_map = block_triangularize( + # graph, top_nodes=con_nodes + #) # Cache maps in case we want to get diagonal blocks quickly in the # future. row_block_map = cnode_block_map diff --git a/pyomo/contrib/incidence_analysis/triangularize.py b/pyomo/contrib/incidence_analysis/triangularize.py index f1a2be737fb..5f909d36364 100644 --- a/pyomo/contrib/incidence_analysis/triangularize.py +++ b/pyomo/contrib/incidence_analysis/triangularize.py @@ -10,126 +10,131 @@ # ___________________________________________________________________________ from pyomo.contrib.incidence_analysis.matching import maximum_matching +from pyomo.contrib.incidence_analysis.common.dulmage_mendelsohn import ( + # TODO: The fact that we import this function here suggests it should be + # promoted. + _get_projected_digraph, +) from pyomo.common.dependencies import networkx as nx -# TODO: get_scc_of_projection(bipartite_graph, top_nodes, matching) +def get_scc_of_projection(graph, top_nodes, matching=None): + """Return the topologically ordered strongly connected components of a + bipartite graph, projected with respect to a perfect matching + The provided undirected bipartite graph is projected into a directed graph + on the set of "top nodes" by treating "matched edges" as in-edges and + "unmatched edges" as out-edges. Then the strongly connected components of + the directed graph are computed. These strongly connected components are + unique, regardless of the choice of perfect matching. The strongly connected + components form a directed acyclic graph, and are returned in a topological + order. The order is unique, as ambiguities are resolved "lexicographically". -def block_triangularize(matrix, matching=None, top_nodes=None): - """ - Computes the necessary information to permute a matrix to block-lower - triangular form, i.e. a partition of rows and columns into an ordered - set of diagonal blocks in such a permutation. - - Arguments - --------- - matrix: A SciPy sparse matrix - matching: A perfect matching of rows and columns, in the form of a dict - mapping row indices to column indices + Parameters + ---------- + graph: NetworkX Graph + A bipartite graph + top_nodes: list + One of the bipartite sets in the graph + matching: dict + Maps each node in ``top_nodes`` to its matched node Returns ------- - Two dicts. The first maps each row index to the index of its block in a - block-lower triangular permutation of the matrix. The second maps each - column index to the index of its block in a block-lower triangular - permutation of the matrix. + """ nxb = nx.algorithms.bipartite nxc = nx.algorithms.components nxd = nx.algorithms.dag - from_biadjacency_matrix = nxb.matrix.from_biadjacency_matrix - - if isinstance(matrix, nx.Graph): - graph_provided = True - if top_nodes is None: - raise RuntimeError( - "top_nodes argument must be set if a graph is provided." - ) - bg = matrix - M = len(top_nodes) - N = len(bg.nodes) - M - if not nxb.is_bipartite(bg): - raise RuntimeError("Provided graph is not bipartite.") - else: - graph_provided = False - M, N = matrix.shape - bg = from_biadjacency_matrix(matrix) - + if not nxb.is_bipartite(graph): + raise RuntimeError("Provided graph is not bipartite.") + M = len(top_nodes) + N = len(graph.nodes) - M if M != N: raise ValueError( - "block_triangularize does not currently " - "support non-square matrices. Got matrix with shape %s." - % ((M, N),) + "get_scc_of_projection does not support bipartite graphs with" + " bipartite sets of different cardinalities. Got sizes %s and" + " %s." % (M, N) ) - if matching is None: - matching = maximum_matching(matrix, top_nodes=top_nodes) - - if not graph_provided: - # If we provided a matrix, the matching maps row indices to column - # indices. Update the matching to map to nodes in the graph derived - # from the matrix. - matching = {r: c+M for r, c in matching.items()} - - len_matching = len(matching) - if len_matching != M: + # This matching maps top nodes to "other nodes" *and* other nodes + # back to top nodes. + matching = nxb.maximum_matching(graph, top_nodes=top_nodes) + if len(matching) != 2*M: raise ValueError( - "block_triangularize only supports matrices " - "that have a perfect matching of rows and columns. " - "Cardinality of maximal matching is %s" % len_matching + "get_scc_of_projection does not support bipartite graphs without" + " a perfect matching. Got a graph with %s nodes per bipartite set" + " and a matching of cardinality %s." % (M, (len(matching)/2)) ) - - # Construct directed graph of rows - dg = nx.DiGraph() - if graph_provided: - dg.add_nodes_from(top_nodes) - else: - dg.add_nodes_from(range(M)) - for n in dg.nodes: - col_node = matching[n] - # For all rows that share this column - for neighbor in bg[col_node]: - if neighbor != n: - # Add an edge towards this column's matched row - dg.add_edge(neighbor, n) - - # Partition the rows into strongly connected components (diagonal blocks) + dg = _get_projected_digraph(graph, matching, top_nodes) scc_list = list(nxc.strongly_connected_components(dg)) + n_scc = len(scc_list) node_scc_map = {n: idx for idx, scc in enumerate(scc_list) for n in scc} # Now we need to put the SCCs in the right order. We do this by performing # a topological sort on the DAG of SCCs. dag = nx.DiGraph() - for i, c in enumerate(scc_list): - dag.add_node(i) + dag.add_nodes_from(range(n_scc)) for n in dg.nodes: source_scc = node_scc_map[n] for neighbor in dg[n]: target_scc = node_scc_map[neighbor] if target_scc != source_scc: - dag.add_edge(target_scc, source_scc) + dag.add_edge(source_scc, target_scc) + + #dag.add_edge(target_scc, source_scc) # Reverse direction of edge. This corresponds to creating # a block lower triangular matrix. scc_order = list(nxd.lexicographical_topological_sort(dag)) - scc_block_map = {c: i for i, c in enumerate(scc_order)} - row_block_map = {n: scc_block_map[c] for n, c in node_scc_map.items()} - # ^ This maps row indices to the blocks they belong to. - - # Invert the matching to map row indices to column indices - if graph_provided: - col_row_map = {c: r for r, c in matching.items()} - col_block_map = { - c: row_block_map[col_row_map[c]] for c in matching.values() - } - else: - col_row_map = {c-M: r for r, c in matching.items()} - col_block_map = {c: row_block_map[col_row_map[c]] for c in range(N)} - assert len(col_row_map) == M - - return row_block_map, col_block_map + # Reverse for compatibility + scc_order.reverse() + + scc_idx_map = {c: i for i, c in enumerate(scc_order)} + top_idx_map = {n: scc_idx_map[c] for n, c in node_scc_map.items()} + bot_idx_map = {matching[n]: scc_idx_map[c] for n, c in node_scc_map.items()} + return top_idx_map, bot_idx_map + + +def block_triangularize(matrix, matching=None, top_nodes=None): + """ + Computes the necessary information to permute a matrix to block-lower + triangular form, i.e. a partition of rows and columns into an ordered + set of diagonal blocks in such a permutation. + + Arguments + --------- + matrix: A SciPy sparse matrix + matching: A perfect matching of rows and columns, in the form of a dict + mapping row indices to column indices + + Returns + ------- + Two dicts. The first maps each row index to the index of its block in a + block-lower triangular permutation of the matrix. The second maps each + column index to the index of its block in a block-lower triangular + permutation of the matrix. + """ + nxb = nx.algorithms.bipartite + nxc = nx.algorithms.components + nxd = nx.algorithms.dag + from_biadjacency_matrix = nxb.matrix.from_biadjacency_matrix + M, N = matrix.shape + if M != N: + raise ValueError( + "block_triangularize does not currently " + "support non-square matrices. Got matrix with shape %s." + % ((M, N),) + ) + graph = from_biadjacency_matrix(matrix) + row_nodes = list(range(M)) + rnode_idx_map, cnode_idx_map = get_scc_of_projection( + graph, row_nodes, matching=matching + ) + row_idx_map = rnode_idx_map + col_idx_map = {n-M: idx for n, idx in cnode_idx_map.items()} + return row_idx_map, col_idx_map def get_blocks_from_maps(row_block_map, col_block_map): From bf1d8a27dd2353c8e14cd3d5c65fb600e95c1e4e Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Wed, 1 Feb 2023 13:00:26 -0700 Subject: [PATCH 13/25] update get_scc_of_projection to return a list of lists of tuples --- pyomo/contrib/incidence_analysis/interface.py | 13 ++++------ .../incidence_analysis/triangularize.py | 25 ++++++++++--------- 2 files changed, 18 insertions(+), 20 deletions(-) diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py index 757b7955811..1626eb06491 100644 --- a/pyomo/contrib/incidence_analysis/interface.py +++ b/pyomo/contrib/incidence_analysis/interface.py @@ -487,16 +487,13 @@ def block_triangularize(self, variables=None, constraints=None): M = len(constraints) con_nodes = list(range(M)) - cnode_block_map, vnode_block_map = get_scc_of_projection( - graph, con_nodes - ) - #cnode_block_map, vnode_block_map = block_triangularize( - # graph, top_nodes=con_nodes - #) + sccs = get_scc_of_projection(graph, con_nodes) + row_idx_map = {r: idx for idx, scc in enumerate(sccs) for r, _ in scc} + col_idx_map = {c-M: idx for idx, scc in enumerate(sccs) for _, c in scc} # Cache maps in case we want to get diagonal blocks quickly in the # future. - row_block_map = cnode_block_map - col_block_map = {j-M: idx for j, idx in vnode_block_map.items()} + row_block_map = row_idx_map + col_block_map = {j-M: idx for j, idx in col_idx_map.items()} self.row_block_map = row_block_map self.col_block_map = col_block_map con_block_map = ComponentMap( diff --git a/pyomo/contrib/incidence_analysis/triangularize.py b/pyomo/contrib/incidence_analysis/triangularize.py index 5f909d36364..2473f23b62a 100644 --- a/pyomo/contrib/incidence_analysis/triangularize.py +++ b/pyomo/contrib/incidence_analysis/triangularize.py @@ -91,10 +91,14 @@ def get_scc_of_projection(graph, top_nodes, matching=None): # Reverse for compatibility scc_order.reverse() - scc_idx_map = {c: i for i, c in enumerate(scc_order)} - top_idx_map = {n: scc_idx_map[c] for n, c in node_scc_map.items()} - bot_idx_map = {matching[n]: scc_idx_map[c] for n, c in node_scc_map.items()} - return top_idx_map, bot_idx_map + # The "natural" return type, here, is a list of lists. Each inner list + # is an SCC, and contains tuples of nodes. The "top node", and its matched + # "bottom node". + ordered_node_subsets = [ + sorted([(i, matching[i]) for i in scc_list[scc_idx]]) + for scc_idx in scc_order + ] + return ordered_node_subsets def block_triangularize(matrix, matching=None, top_nodes=None): @@ -123,17 +127,14 @@ def block_triangularize(matrix, matching=None, top_nodes=None): M, N = matrix.shape if M != N: raise ValueError( - "block_triangularize does not currently " - "support non-square matrices. Got matrix with shape %s." - % ((M, N),) + "block_triangularize does not currently support non-square" + " matrices. Got matrix with shape %s." % ((M, N),) ) graph = from_biadjacency_matrix(matrix) row_nodes = list(range(M)) - rnode_idx_map, cnode_idx_map = get_scc_of_projection( - graph, row_nodes, matching=matching - ) - row_idx_map = rnode_idx_map - col_idx_map = {n-M: idx for n, idx in cnode_idx_map.items()} + sccs = get_scc_of_projection(graph, row_nodes, matching=matching) + row_idx_map = {r: idx for idx, scc in enumerate(sccs) for r, _ in scc} + col_idx_map = {c-M: idx for idx, scc in enumerate(sccs) for _, c in scc} return row_idx_map, col_idx_map From 1bdea8f0d94a5552a4eb01790efc2b376dd681e8 Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Wed, 1 Feb 2023 13:04:06 -0700 Subject: [PATCH 14/25] update docstring --- pyomo/contrib/incidence_analysis/triangularize.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pyomo/contrib/incidence_analysis/triangularize.py b/pyomo/contrib/incidence_analysis/triangularize.py index 2473f23b62a..1d44b7b9c52 100644 --- a/pyomo/contrib/incidence_analysis/triangularize.py +++ b/pyomo/contrib/incidence_analysis/triangularize.py @@ -41,6 +41,10 @@ def get_scc_of_projection(graph, top_nodes, matching=None): Returns ------- + list of lists + The outer list is a list of strongly connected components. Each + strongly connected component is a list of tuples of matched nodes. + The first node is a "top node", and the second is an "other node". """ nxb = nx.algorithms.bipartite @@ -82,10 +86,6 @@ def get_scc_of_projection(graph, top_nodes, matching=None): if target_scc != source_scc: dag.add_edge(source_scc, target_scc) - #dag.add_edge(target_scc, source_scc) - # Reverse direction of edge. This corresponds to creating - # a block lower triangular matrix. - scc_order = list(nxd.lexicographical_topological_sort(dag)) # Reverse for compatibility From 584bcbbab1b5dc5cf27db44c4235188fb761451d Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Wed, 1 Feb 2023 13:40:08 -0700 Subject: [PATCH 15/25] tests for get_scc_of_projection --- .../tests/test_triangularize.py | 114 ++++++++++++++++-- .../incidence_analysis/triangularize.py | 4 +- 2 files changed, 108 insertions(+), 10 deletions(-) diff --git a/pyomo/contrib/incidence_analysis/tests/test_triangularize.py b/pyomo/contrib/incidence_analysis/tests/test_triangularize.py index f7542fcb6dc..9165f26f462 100644 --- a/pyomo/contrib/incidence_analysis/tests/test_triangularize.py +++ b/pyomo/contrib/incidence_analysis/tests/test_triangularize.py @@ -12,20 +12,118 @@ import random from pyomo.contrib.incidence_analysis.matching import maximum_matching from pyomo.contrib.incidence_analysis.triangularize import ( - block_triangularize, - get_diagonal_blocks, - ) + get_scc_of_projection, + block_triangularize, + get_diagonal_blocks, +) from pyomo.common.dependencies import ( - scipy, - scipy_available, - networkx_available, - ) + scipy, + scipy_available, + networkx as nx, + networkx_available, +) if scipy_available: sps = scipy.sparse +if networkx_available: + nxb = nx.algorithms.bipartite import pyomo.common.unittest as unittest +@unittest.skipUnless(networkx_available, "networkx is not available") +@unittest.skipUnless(scipy_available, "scipy is not available") +class TestGetSCCOfProjection(unittest.TestCase): + + def test_graph_decomposable_tridiagonal_shuffled(self): + """ + This is the same graph as in test_decomposable_tridiagonal_shuffled + below, but now we convert the matrix into a bipartite graph and + use get_scc_of_projection. + + The matrix decomposes into 2x2 blocks: + |x x | + |x x | + | x x x | + | x x | + | x x| + """ + N = 11 + row = [] + col = [] + data = [] + + # Diagonal + row.extend(range(N)) + col.extend(range(N)) + data.extend(1 for _ in range(N)) + + # Below diagonal + row.extend(range(1, N)) + col.extend(range(N-1)) + data.extend(1 for _ in range(N-1)) + + # Above diagonal + row.extend(i for i in range(N-1) if not i%2) + col.extend(i+1 for i in range(N-1) if not i%2) + data.extend(1 for i in range(N-1) if not i%2) + + # Same results hold after applying a random permutation. + row_perm = list(range(N)) + col_perm = list(range(N)) + random.shuffle(row_perm) + random.shuffle(col_perm) + + row = [row_perm[i] for i in row] + col = [col_perm[j] for j in col] + + matrix = sps.coo_matrix((data, (row, col)), shape=(N, N)) + graph = nxb.matrix.from_biadjacency_matrix(matrix) + row_nodes = list(range(N)) + sccs = get_scc_of_projection(graph, row_nodes) + + self.assertEqual(len(sccs), (N+1)//2) + + for i in range((N+1)//2): + # Note that these rows and cols are in the permuted space + rows = set(r for r, _ in sccs[i]) + cols = set(c-N for _, c in sccs[i]) + + pred_rows = {row_perm[2*i]} + pred_cols = {col_perm[2*i]} + + if 2*i+1 < N: + pred_rows.add(row_perm[2*i+1]) + pred_cols.add(col_perm[2*i+1]) + + self.assertEqual(pred_rows, rows) + self.assertEqual(pred_cols, cols) + + def test_scc_exceptions(self): + graph = nx.Graph() + graph.add_nodes_from(range(3)) + graph.add_edges_from([(0, 1), (0, 2), (1, 2)]) + top_nodes = [0] + msg = "graph is not bipartite" + with self.assertRaisesRegex(RuntimeError, msg): + sccs = get_scc_of_projection(graph, top_nodes=top_nodes) + + graph = nx.Graph() + graph.add_nodes_from(range(3)) + graph.add_edges_from([(0, 1), (0, 2)]) + top_nodes[0] + msg = "bipartite sets of different cardinalities" + with self.assertRaisesRegex(RuntimeError, msg): + sccs = get_scc_of_projection(graph, top_nodes=top_nodes) + + graph = nx.Graph() + graph.add_nodes_from(range(4)) + graph.add_edges_from([(0, 1), (0, 2)]) + top_nodes = [0, 3] + msg = "without a perfect matching" + with self.assertRaisesRegex(RuntimeError, msg): + sccs = get_scc_of_projection(graph, top_nodes=top_nodes) + + @unittest.skipUnless(networkx_available, "networkx is not available") @unittest.skipUnless(scipy_available, "scipy is not available") class TestTriangularize(unittest.TestCase): @@ -37,7 +135,7 @@ def test_low_rank_exception(self): matrix = sps.coo_matrix((data, (row, col)), shape=(N, N)) - with self.assertRaises(ValueError) as exc: + with self.assertRaises(RuntimeError) as exc: row_block_map, col_block_map = block_triangularize(matrix) self.assertIn('perfect matching', str(exc.exception)) diff --git a/pyomo/contrib/incidence_analysis/triangularize.py b/pyomo/contrib/incidence_analysis/triangularize.py index 1d44b7b9c52..90b5c1cdbda 100644 --- a/pyomo/contrib/incidence_analysis/triangularize.py +++ b/pyomo/contrib/incidence_analysis/triangularize.py @@ -55,7 +55,7 @@ def get_scc_of_projection(graph, top_nodes, matching=None): M = len(top_nodes) N = len(graph.nodes) - M if M != N: - raise ValueError( + raise RuntimeError( "get_scc_of_projection does not support bipartite graphs with" " bipartite sets of different cardinalities. Got sizes %s and" " %s." % (M, N) @@ -65,7 +65,7 @@ def get_scc_of_projection(graph, top_nodes, matching=None): # back to top nodes. matching = nxb.maximum_matching(graph, top_nodes=top_nodes) if len(matching) != 2*M: - raise ValueError( + raise RuntimeError( "get_scc_of_projection does not support bipartite graphs without" " a perfect matching. Got a graph with %s nodes per bipartite set" " and a matching of cardinality %s." % (M, (len(matching)/2)) From 6d1106a15978fdba8e83e9cca0375ccd50595378 Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Wed, 1 Feb 2023 13:43:50 -0700 Subject: [PATCH 16/25] re-add self.cached=NUMERIC, even though it doesnt do anything anymore --- pyomo/contrib/incidence_analysis/interface.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py index 1626eb06491..63ead9cc52c 100644 --- a/pyomo/contrib/incidence_analysis/interface.py +++ b/pyomo/contrib/incidence_analysis/interface.py @@ -46,7 +46,6 @@ class IncidenceMatrixType(enum.Enum): NONE = 0 STRUCTURAL = 1 NUMERIC = 2 - # GRAPH = 3 def _check_unindexed(complist): @@ -244,6 +243,7 @@ def __init__( self.cached = IncidenceMatrixType.NONE self.incidence_graph = None elif isinstance(model, PyomoNLP): + self.cached = IncidenceMatrixType.NUMERIC if not active: raise ValueError( "Cannot get the Jacobian of inactive constraints from the " From cd829b5ac5434900259922c0757fb9e7377cfb50 Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Wed, 1 Feb 2023 14:14:48 -0700 Subject: [PATCH 17/25] remove extra self.cached=NUMERIC; turns out I never removed it --- pyomo/contrib/incidence_analysis/interface.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py index 63ead9cc52c..081e0c956ef 100644 --- a/pyomo/contrib/incidence_analysis/interface.py +++ b/pyomo/contrib/incidence_analysis/interface.py @@ -243,7 +243,6 @@ def __init__( self.cached = IncidenceMatrixType.NONE self.incidence_graph = None elif isinstance(model, PyomoNLP): - self.cached = IncidenceMatrixType.NUMERIC if not active: raise ValueError( "Cannot get the Jacobian of inactive constraints from the " From f3feb5d72c4acd26dfa0bb1ba0dcd421c301741d Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Wed, 1 Feb 2023 14:28:45 -0700 Subject: [PATCH 18/25] use len(self.constraints) as this is more clear than len(con_index_map) --- pyomo/contrib/incidence_analysis/interface.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py index 081e0c956ef..46d5c0073f0 100644 --- a/pyomo/contrib/incidence_analysis/interface.py +++ b/pyomo/contrib/incidence_analysis/interface.py @@ -363,7 +363,10 @@ def _extract_subgraph(self, variables, constraints): ) else: constraint_nodes = [self.con_index_map[con] for con in constraints] - M = len(self.con_index_map) + + # Note that this is the number of constraints in the original graph, + # not the subgraph. + M = len(self.constraints) variable_nodes = [M + self.var_index_map[var] for var in variables] subgraph = extract_bipartite_subgraph( self.incidence_graph, constraint_nodes, variable_nodes From 84f43de7e39524721cc374e439f2d8b672b9eb68 Mon Sep 17 00:00:00 2001 From: Robert Parker Date: Wed, 1 Feb 2023 14:30:41 -0700 Subject: [PATCH 19/25] test to exercise _extract_subgraph when a matrix is cached --- .../tests/test_interface.py | 26 +++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/pyomo/contrib/incidence_analysis/tests/test_interface.py b/pyomo/contrib/incidence_analysis/tests/test_interface.py index 08bfffcf493..4668e6344a7 100644 --- a/pyomo/contrib/incidence_analysis/tests/test_interface.py +++ b/pyomo/contrib/incidence_analysis/tests/test_interface.py @@ -1449,5 +1449,31 @@ def test_get_adjacent_exceptions(self): adj_cons = igraph.get_adjacent_to(m.x[1]) +@unittest.skipUnless(networkx_available, "networkx is not available.") +@unittest.skipUnless(scipy_available, "scipy is not available.") +class TestInterface(unittest.TestCase): + + def test_subgraph_with_fewer_var_or_con(self): + m = pyo.ConcreteModel() + m.I = pyo.Set(initialize=[1, 2]) + m.v = pyo.Var(m.I) + m.eq1 = pyo.Constraint(expr=m.v[1] + m.v[2] == 1) + m.ineq1 = pyo.Constraint(expr=m.v[1] - m.v[2] <= 2) + + # Defensively set include_inequality=True, which is the current + # default, in case this default changes. + igraph = IncidenceGraphInterface(m, include_inequality=True) + + variables = list(m.v.values()) + constraints = [m.ineq1] + matching = igraph.maximum_matching(variables, constraints) + self.assertEqual(len(matching), 1) + + variables = [m.v[2]] + constraints = [m.eq1, m.ineq1] + matching = igraph.maximum_matching(variables, constraints) + self.assertEqual(len(matching), 1) + + if __name__ == "__main__": unittest.main() From f968efc87dab98f37c44347b494e33244ec7ff30 Mon Sep 17 00:00:00 2001 From: robbybp Date: Thu, 2 Feb 2023 16:53:29 -0700 Subject: [PATCH 20/25] fix typos and remove extraneous include_fixed arguments --- pyomo/contrib/incidence_analysis/interface.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py index 46d5c0073f0..a9efe164bec 100644 --- a/pyomo/contrib/incidence_analysis/interface.py +++ b/pyomo/contrib/incidence_analysis/interface.py @@ -107,8 +107,8 @@ def extract_bipartite_subgraph(graph, nodes0, nodes1): Two lists of nodes to project onto must be provided. These will correspond to the "bipartite sets" in the subgraph. If the two sets provided have - M and N nodes, the subgraph will have nodes 0 through M+N, with the first - M corresponding to the first set provided and the last M corresponding + M and N nodes, the subgraph will have nodes 0 through M+N-1, with the first + M corresponding to the first set provided and the last N corresponding to the second set. Parameters @@ -117,10 +117,10 @@ def extract_bipartite_subgraph(graph, nodes0, nodes1): The graph from which a subgraph is extracted nodes0: list A list of nodes in the original graph that will form the first - bipartite set of the projected graph (and have ``bipartite=0`` + bipartite set of the projected graph (and have ``bipartite=0``) nodes1: list A list of nodes in the original graph that will form the second - bipartite set of the projected graph (and have ``bipartite=1`` + bipartite set of the projected graph (and have ``bipartite=1``) Returns ------- @@ -296,7 +296,8 @@ def __init__( self.incidence_graph = get_bipartite_incidence_graph( self.variables, self.constraints, - # TODO: include_fixed=include_fixed? + # Note that include_fixed is not necessary here. We have + # already checked this condition above. ) else: raise TypeError( @@ -358,8 +359,9 @@ def _extract_submatrix(self, variables, constraints): def _extract_subgraph(self, variables, constraints): if self.incidence_graph is None: return get_bipartite_incidence_graph( - # Does include_fixed matter here if I'm providing the variables? - variables, constraints, include_fixed=False + variables, constraints + # Note that, as variables are explicitly specified, there + # is no need for an include_fixed argument. ) else: constraint_nodes = [self.con_index_map[con] for con in constraints] From 0d0d08be8cc652116a8b06f495579dae5921aa42 Mon Sep 17 00:00:00 2001 From: robbybp Date: Thu, 2 Feb 2023 17:07:26 -0700 Subject: [PATCH 21/25] remove include_fixed=False in _extract_submatrix --- pyomo/contrib/incidence_analysis/interface.py | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py index a9efe164bec..2d8fd7a13b0 100644 --- a/pyomo/contrib/incidence_analysis/interface.py +++ b/pyomo/contrib/incidence_analysis/interface.py @@ -326,11 +326,9 @@ def _validate_input(self, variables, constraints): def _extract_submatrix(self, variables, constraints): # Assumes variables and constraints are valid if self.incidence_graph is None: - return get_structural_incidence_matrix( - variables, - constraints, - include_fixed=False, - ) + # Note that, as variables are explicitly specified, there + # is no need for an include_fixed argument. + return get_structural_incidence_matrix(variables, constraints) else: N = len(variables) M = len(constraints) @@ -358,11 +356,9 @@ def _extract_submatrix(self, variables, constraints): def _extract_subgraph(self, variables, constraints): if self.incidence_graph is None: - return get_bipartite_incidence_graph( - variables, constraints - # Note that, as variables are explicitly specified, there - # is no need for an include_fixed argument. - ) + # Note that, as variables are explicitly specified, there + # is no need for an include_fixed argument. + return get_bipartite_incidence_graph(variables, constraints) else: constraint_nodes = [self.con_index_map[con] for con in constraints] From 361b1ebd8323fc15a161737ef0d15d44b048884f Mon Sep 17 00:00:00 2001 From: robbybp Date: Thu, 2 Feb 2023 17:23:01 -0700 Subject: [PATCH 22/25] update get_diagonal_blocks to use get_scc_of_projection and remove unused _extract_submatrix --- pyomo/contrib/incidence_analysis/interface.py | 49 +++---------------- 1 file changed, 6 insertions(+), 43 deletions(-) diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py index 2d8fd7a13b0..cc6b6ab9428 100644 --- a/pyomo/contrib/incidence_analysis/interface.py +++ b/pyomo/contrib/incidence_analysis/interface.py @@ -323,37 +323,6 @@ def _validate_input(self, variables, constraints): _check_unindexed(variables + constraints) return variables, constraints - def _extract_submatrix(self, variables, constraints): - # Assumes variables and constraints are valid - if self.incidence_graph is None: - # Note that, as variables are explicitly specified, there - # is no need for an include_fixed argument. - return get_structural_incidence_matrix(variables, constraints) - else: - N = len(variables) - M = len(constraints) - old_new_var_indices = dict( - (self.var_index_map[v], i) for i, v in enumerate(variables) - ) - old_new_con_indices = dict( - (self.con_index_map[c], i) for i, c in enumerate(constraints) - ) - # FIXME: This will fail if I don't have an incidence matrix - # cached. - coo = self.incidence_matrix - new_row = [] - new_col = [] - new_data = [] - for r, c, e in zip(coo.row, coo.col, coo.data): - if r in old_new_con_indices and c in old_new_var_indices: - new_row.append(old_new_con_indices[r]) - new_col.append(old_new_var_indices[c]) - new_data.append(e) - return sp.sparse.coo_matrix( - (new_data, (new_row, new_col)), - shape=(M, N), - ) - def _extract_subgraph(self, variables, constraints): if self.incidence_graph is None: # Note that, as variables are explicitly specified, there @@ -521,18 +490,12 @@ def get_diagonal_blocks(self, variables=None, constraints=None): """ variables, constraints = self._validate_input(variables, constraints) - matrix = self._extract_submatrix(variables, constraints) - - # TODO: Again, this functionality does not really make sense for - # general bipartite graphs... - if self.row_block_map is None or self.col_block_map is None: - block_rows, block_cols = get_diagonal_blocks(matrix) - else: - block_rows, block_cols = get_blocks_from_maps( - self.row_block_map, self.col_block_map - ) - block_cons = [[constraints[i] for i in block] for block in block_rows] - block_vars = [[variables[i] for i in block] for block in block_cols] + graph = self._extract_subgraph(variables, constraints) + M = len(constraints) + con_nodes = list(range(M)) + sccs = get_scc_of_projection(graph, con_nodes) + block_cons = [[constraints[i] for i, _ in scc] for scc in sccs] + block_vars = [[variables[j-M] for _, j in scc] for scc in sccs] return block_vars, block_cons def dulmage_mendelsohn(self, variables=None, constraints=None): From 62b78b44be3f10f2280528e506659d020dbb72ec Mon Sep 17 00:00:00 2001 From: robbybp Date: Thu, 2 Feb 2023 17:27:59 -0700 Subject: [PATCH 23/25] comment on unused row_block_map attributes --- pyomo/contrib/incidence_analysis/interface.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pyomo/contrib/incidence_analysis/interface.py b/pyomo/contrib/incidence_analysis/interface.py index cc6b6ab9428..6ddb1c546d0 100644 --- a/pyomo/contrib/incidence_analysis/interface.py +++ b/pyomo/contrib/incidence_analysis/interface.py @@ -305,6 +305,9 @@ def __init__( "%s or %s but got %s." % (PyomoNLP, Block, type(model)) ) + # Note that these attributes are no longer used in get_diagonal_blocks, + # and should be removed. However, they are tested, so we will leave + # them in until a breaking change. self.row_block_map = None self.col_block_map = None From 0448c32a5fb3f750fcfa7a5dd6c236e24cd27d98 Mon Sep 17 00:00:00 2001 From: robbybp Date: Thu, 2 Feb 2023 18:00:24 -0700 Subject: [PATCH 24/25] remove unused top_nodes argument from block_triangularize --- pyomo/contrib/incidence_analysis/triangularize.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/contrib/incidence_analysis/triangularize.py b/pyomo/contrib/incidence_analysis/triangularize.py index 90b5c1cdbda..a4096ae0a04 100644 --- a/pyomo/contrib/incidence_analysis/triangularize.py +++ b/pyomo/contrib/incidence_analysis/triangularize.py @@ -101,7 +101,7 @@ def get_scc_of_projection(graph, top_nodes, matching=None): return ordered_node_subsets -def block_triangularize(matrix, matching=None, top_nodes=None): +def block_triangularize(matrix, matching=None): """ Computes the necessary information to permute a matrix to block-lower triangular form, i.e. a partition of rows and columns into an ordered From 5661e84e1cfecb041332558d2437422e7340ded8 Mon Sep 17 00:00:00 2001 From: robbybp Date: Thu, 2 Feb 2023 18:41:43 -0700 Subject: [PATCH 25/25] update docstring for correctness, then change location of edge reversal to match docstring --- .../contrib/incidence_analysis/triangularize.py | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/pyomo/contrib/incidence_analysis/triangularize.py b/pyomo/contrib/incidence_analysis/triangularize.py index a4096ae0a04..8007cd8cb9b 100644 --- a/pyomo/contrib/incidence_analysis/triangularize.py +++ b/pyomo/contrib/incidence_analysis/triangularize.py @@ -23,13 +23,17 @@ def get_scc_of_projection(graph, top_nodes, matching=None): bipartite graph, projected with respect to a perfect matching The provided undirected bipartite graph is projected into a directed graph - on the set of "top nodes" by treating "matched edges" as in-edges and - "unmatched edges" as out-edges. Then the strongly connected components of + on the set of "top nodes" by treating "matched edges" as out-edges and + "unmatched edges" as in-edges. Then the strongly connected components of the directed graph are computed. These strongly connected components are unique, regardless of the choice of perfect matching. The strongly connected components form a directed acyclic graph, and are returned in a topological order. The order is unique, as ambiguities are resolved "lexicographically". + The "direction" of the projection (where matched edges are out-edges) + leads to a block *lower* triangular permutation when the top nodes + correspond to *rows* in the bipartite graph of a matrix. + Parameters ---------- graph: NetworkX Graph @@ -70,7 +74,11 @@ def get_scc_of_projection(graph, top_nodes, matching=None): " a perfect matching. Got a graph with %s nodes per bipartite set" " and a matching of cardinality %s." % (M, (len(matching)/2)) ) - dg = _get_projected_digraph(graph, matching, top_nodes) + + # _get_projected_digraph treats matched edges as "in-edges", so we + # reverse the direction of edges here. + dg = _get_projected_digraph(graph, matching, top_nodes).reverse() + scc_list = list(nxc.strongly_connected_components(dg)) n_scc = len(scc_list) node_scc_map = {n: idx for idx, scc in enumerate(scc_list) for n in scc} @@ -88,9 +96,6 @@ def get_scc_of_projection(graph, top_nodes, matching=None): scc_order = list(nxd.lexicographical_topological_sort(dag)) - # Reverse for compatibility - scc_order.reverse() - # The "natural" return type, here, is a list of lists. Each inner list # is an SCC, and contains tuples of nodes. The "top node", and its matched # "bottom node".