1# ___________________________________________________________________________ 2# 3# Pyomo: Python Optimization Modeling Objects 4# Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC 5# Under the terms of Contract DE-NA0003525 with National Technology and 6# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain 7# rights in this software. 8# This software is distributed under the 3-clause BSD License. 9# ___________________________________________________________________________ 10 11import enum 12from pyomo.core.base.block import Block 13from pyomo.core.base.var import Var 14from pyomo.core.base.constraint import Constraint 15from pyomo.core.base.objective import Objective 16from pyomo.core.base.reference import Reference 17from pyomo.core.expr.visitor import identify_variables 18from pyomo.common.collections import ComponentSet, ComponentMap 19from pyomo.common.dependencies import scipy_available 20from pyomo.common.dependencies import networkx as nx 21from pyomo.contrib.incidence_analysis.matching import maximum_matching 22from pyomo.contrib.incidence_analysis.triangularize import block_triangularize 23from pyomo.contrib.incidence_analysis.dulmage_mendelsohn import ( 24 dulmage_mendelsohn, 25 ) 26if scipy_available: 27 from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP 28 import scipy as sp 29 30 31class IncidenceMatrixType(enum.Enum): 32 NONE = 0 33 STRUCTURAL = 1 34 NUMERIC = 2 35 36 37def _check_unindexed(complist): 38 for comp in complist: 39 if comp.is_indexed(): 40 raise ValueError( 41 "Variables and constraints must be unindexed " 42 "ComponentData objects. Got %s, which is indexed." 43 % comp.name 44 ) 45 46 47def get_incidence_graph(variables, constraints, include_fixed=True): 48 """ 49 This function gets the incidence graph of Pyomo variables and constraints. 50 51 Arguments: 52 ---------- 53 variables: List of Pyomo VarData objects 54 Variables that will appear in incidence graph 55 constraints: List of Pyomo ConstraintData objects 56 Constraints that will appear in incidence graph 57 include_fixed: Bool 58 Flag for whether fixed variable should be included in the incidence 59 60 Returns: 61 -------- 62 NetworkX Graph 63 64 """ 65 _check_unindexed(variables+constraints) 66 N, M = len(variables), len(constraints) 67 graph = nx.Graph() 68 graph.add_nodes_from(range(M), bipartite=0) 69 graph.add_nodes_from(range(M, M+N), bipartite=1) 70 var_node_map = ComponentMap((v, M+i) for i, v in enumerate(variables)) 71 for i, con in enumerate(constraints): 72 for var in identify_variables(con.body, include_fixed=include_fixed): 73 if var in var_node_map: 74 graph.add_edge(i, var_node_map[var]) 75 return graph 76 77 78def get_structural_incidence_matrix(variables, constraints, include_fixed=True): 79 """ 80 This function gets the incidence matrix of Pyomo constraints and variables. 81 82 Arguments 83 --------- 84 variables: List of Pyomo VarData objects 85 constraints: List of Pyomo ConstraintData objects 86 include_fixed: Bool 87 Flag for whether fixed variables should be included in the matrix 88 nonzeros 89 90 Returns 91 ------- 92 A scipy.sparse coo matrix. Rows are indices into the user-provided list of 93 constraints, columns are indices into the user-provided list of variables. 94 Entries are 1.0. 95 96 """ 97 _check_unindexed(variables+constraints) 98 N, M = len(variables), len(constraints) 99 var_idx_map = ComponentMap((v, i) for i, v in enumerate(variables)) 100 rows = [] 101 cols = [] 102 for i, con in enumerate(constraints): 103 cols.extend(var_idx_map[v] for v in 104 identify_variables(con.body, include_fixed=include_fixed) 105 if v in var_idx_map) 106 rows.extend([i]*(len(cols) - len(rows))) 107 assert len(rows) == len(cols) 108 data = [1.0]*len(rows) 109 matrix = sp.sparse.coo_matrix( (data, (rows, cols)), shape=(M, N) ) 110 return matrix 111 112 113def get_numeric_incidence_matrix(variables, constraints): 114 """ 115 This function gets the numeric incidence matrix (Jacobian) of Pyomo 116 constraints with respect to variables. 117 """ 118 # NOTE: There are several ways to get a numeric incidence matrix 119 # from a Pyomo model. This function implements a somewhat roundabout 120 # method, which is to construct a dummy Block with the necessary 121 # variables and constraints, then construct a PyNumero PyomoNLP 122 # from the block and have PyNumero evaluate the desired Jacobian 123 # via ASL. 124 comps = list(variables) + list(constraints) 125 _check_unindexed(comps) 126 M, N = len(constraints), len(variables) 127 _block = Block() 128 _block.construct() 129 _block.obj = Objective(expr=0) 130 _block.vars = Reference(variables) 131 _block.cons = Reference(constraints) 132 var_set = ComponentSet(variables) 133 other_vars = [] 134 for con in constraints: 135 for var in identify_variables(con.body, include_fixed=False): 136 # Fixed vars will be ignored by the nl file write, so 137 # there is no point to including them here. 138 # A different method of assembling this matrix, e.g. 139 # Pyomo's automatic differentiation, could support taking 140 # derivatives with respect to fixed variables. 141 if var not in var_set: 142 other_vars.append(var) 143 var_set.add(var) 144 # These variables are necessary due to the nl writer's philosophy 145 # about what constitutes a model. Note that we take derivatives with 146 # respect to them even though this is not necessary. We could fix them 147 # here to avoid doing this extra work, but that would alter the user's 148 # model, which we would rather not do. 149 _block.other_vars = Reference(other_vars) 150 _nlp = PyomoNLP(_block) 151 return _nlp.extract_submatrix_jacobian(variables, constraints) 152 153 154class IncidenceGraphInterface(object): 155 """ 156 The purpose of this class is to allow the user to easily 157 analyze graphs of variables and contraints in a Pyomo 158 model without constructing multiple PyomoNLPs. 159 """ 160 161 def __init__(self, model=None): 162 """ 163 """ 164 # If the user gives us a model or an NLP, we assume they want us 165 # to cache the incidence matrix for fast analysis of submatrices 166 # later on. 167 # WARNING: This cache will become invalid if the user alters their 168 # model. 169 if model is None: 170 self.cached = IncidenceMatrixType.NONE 171 elif isinstance(model, PyomoNLP): 172 nlp = model 173 self.cached = IncidenceMatrixType.NUMERIC 174 self.variables = nlp.get_pyomo_variables() 175 self.constraints = nlp.get_pyomo_constraints() 176 self.var_index_map = ComponentMap( 177 (var, idx) for idx, var in enumerate(self.variables)) 178 self.con_index_map = ComponentMap( 179 (con, idx) for idx, con in enumerate(self.constraints)) 180 self.incidence_matrix = nlp.evaluate_jacobian_eq() 181 elif isinstance(model, Block): 182 self.cached = IncidenceMatrixType.STRUCTURAL 183 self.variables = list(model.component_data_objects(Var)) 184 self.constraints = list(model.component_data_objects(Constraint)) 185 self.var_index_map = ComponentMap( 186 (var, i) for i, var in enumerate(self.variables)) 187 self.con_index_map = ComponentMap( 188 (con, i) for i, con in enumerate(self.constraints)) 189 self.incidence_matrix = get_structural_incidence_matrix( 190 self.variables, 191 self.constraints, 192 ) 193 else: 194 raise TypeError( 195 "Unsupported type for incidence matrix. Expected " 196 "%s or %s but got %s." 197 % (PyomoNLP, Block, type(model)) 198 ) 199 200 def _validate_input(self, variables, constraints): 201 if variables is None: 202 if self.cached is IncidenceMatrixType.NONE: 203 raise ValueError( 204 "Neither variables nor a model have been provided." 205 ) 206 else: 207 variables = self.variables 208 if constraints is None: 209 if self.cached is IncidenceMatrixType.NONE: 210 raise ValueError( 211 "Neither constraints nor a model have been provided." 212 ) 213 else: 214 constraints = self.constraints 215 216 _check_unindexed(variables+constraints) 217 return variables, constraints 218 219 def _extract_submatrix(self, variables, constraints): 220 # Assumes variables and constraints are valid 221 if self.cached is IncidenceMatrixType.NONE: 222 return get_structural_incidence_matrix( 223 variables, 224 constraints, 225 include_fixed=False, 226 ) 227 else: 228 N, M = len(variables), len(constraints) 229 old_new_var_indices = dict((self.var_index_map[v], i) 230 for i, v in enumerate(variables)) 231 old_new_con_indices = dict((self.con_index_map[c], i) 232 for i, c in enumerate(constraints)) 233 coo = self.incidence_matrix 234 new_row = [] 235 new_col = [] 236 new_data = [] 237 for r, c, e in zip(coo.row, coo.col, coo.data): 238 if r in old_new_con_indices and c in old_new_var_indices: 239 new_row.append(old_new_con_indices[r]) 240 new_col.append(old_new_var_indices[c]) 241 new_data.append(e) 242 return sp.sparse.coo_matrix( 243 (new_data, (new_row, new_col)), 244 shape=(M, N), 245 ) 246 247 def maximum_matching(self, variables=None, constraints=None): 248 """ 249 Returns a maximal matching between the constraints and variables, 250 in terms of a map from constraints to variables. 251 """ 252 variables, constraints = self._validate_input(variables, constraints) 253 matrix = self._extract_submatrix(variables, constraints) 254 255 matching = maximum_matching(matrix.tocoo()) 256 # Matching maps row (constraint) indices to column (variable) indices 257 258 return ComponentMap((constraints[i], variables[j]) 259 for i, j in matching.items()) 260 261 def block_triangularize(self, variables=None, constraints=None): 262 """ 263 Returns two ComponentMaps. A map from variables to their blocks 264 in a block triangularization of the incidence matrix, and a 265 map from constraints to their blocks in a block triangularization 266 of the incidence matrix. 267 """ 268 variables, constraints = self._validate_input(variables, constraints) 269 matrix = self._extract_submatrix(variables, constraints) 270 271 row_block_map, col_block_map = block_triangularize(matrix.tocoo()) 272 con_block_map = ComponentMap((constraints[i], idx) 273 for i, idx in row_block_map.items()) 274 var_block_map = ComponentMap((variables[j], idx) 275 for j, idx in col_block_map.items()) 276 # Switch the order of the maps here to match the method call. 277 # Hopefully this does not get too confusing... 278 return var_block_map, con_block_map 279 280 def dulmage_mendelsohn(self, variables=None, constraints=None): 281 """ 282 Returns the Dulmage-Mendelsohn partition of the incidence graph 283 of the provided variables and constraints. 284 285 Returns: 286 -------- 287 ColPartition namedtuple and RowPartition namedtuple. 288 The ColPartition is returned first to match the order of variables 289 and constraints in the method arguments. 290 These partition variables (columns) and constraints (rows) 291 into overconstrained, underconstrained, unmatched, and square. 292 293 """ 294 variables, constraints = self._validate_input(variables, constraints) 295 matrix = self._extract_submatrix(variables, constraints) 296 297 row_partition, col_partition = dulmage_mendelsohn(matrix.tocoo()) 298 con_partition = tuple( 299 [constraints[i] for i in subset] for subset in row_partition 300 ) 301 var_partition = tuple( 302 [variables[i] for i in subset] for subset in col_partition 303 ) 304 # Switch the order of the maps here to match the method call. 305 # Hopefully this does not get too confusing... 306 return var_partition, con_partition 307