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