|
| 1 | +# ___________________________________________________________________________ |
| 2 | +# |
| 3 | +# Pyomo: Python Optimization Modeling Objects |
| 4 | +# Copyright (c) 2008-2024 |
| 5 | +# National Technology and Engineering Solutions of Sandia, LLC |
| 6 | +# Under the terms of Contract DE-NA0003525 with National Technology and |
| 7 | +# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain |
| 8 | +# rights in this software. |
| 9 | +# This software is distributed under the 3-clause BSD License. |
| 10 | +# ___________________________________________________________________________ |
| 11 | +"""Module for visualizing results of incidence graph or matrix analysis |
| 12 | +
|
| 13 | +""" |
| 14 | +from pyomo.contrib.incidence_analysis.config import IncidenceOrder |
| 15 | +from pyomo.contrib.incidence_analysis.interface import ( |
| 16 | + IncidenceGraphInterface, |
| 17 | + get_structural_incidence_matrix, |
| 18 | +) |
| 19 | +from pyomo.common.dependencies import matplotlib |
| 20 | + |
| 21 | + |
| 22 | +def _partition_variables_and_constraints( |
| 23 | + model, order=IncidenceOrder.dulmage_mendelsohn_upper, **kwds |
| 24 | +): |
| 25 | + """Partition variables and constraints in an incidence graph""" |
| 26 | + igraph = IncidenceGraphInterface(model, **kwds) |
| 27 | + vdmp, cdmp = igraph.dulmage_mendelsohn() |
| 28 | + |
| 29 | + ucv = vdmp.unmatched + vdmp.underconstrained |
| 30 | + ucc = cdmp.underconstrained |
| 31 | + |
| 32 | + ocv = vdmp.overconstrained |
| 33 | + occ = cdmp.overconstrained + cdmp.unmatched |
| 34 | + |
| 35 | + ucvblocks, uccblocks = igraph.get_connected_components( |
| 36 | + variables=ucv, constraints=ucc |
| 37 | + ) |
| 38 | + ocvblocks, occblocks = igraph.get_connected_components( |
| 39 | + variables=ocv, constraints=occ |
| 40 | + ) |
| 41 | + wcvblocks, wccblocks = igraph.block_triangularize( |
| 42 | + variables=vdmp.square, constraints=cdmp.square |
| 43 | + ) |
| 44 | + # By default, we block-*lower* triangularize. By default, however, we want |
| 45 | + # the Dulmage-Mendelsohn decomposition to be block-*upper* triangular. |
| 46 | + wcvblocks.reverse() |
| 47 | + wccblocks.reverse() |
| 48 | + vpartition = [ucvblocks, wcvblocks, ocvblocks] |
| 49 | + cpartition = [uccblocks, wccblocks, occblocks] |
| 50 | + |
| 51 | + if order == IncidenceOrder.dulmage_mendelsohn_lower: |
| 52 | + # If a block-lower triangular matrix was requested, we need to reverse |
| 53 | + # both the inner and outer partitions |
| 54 | + vpartition.reverse() |
| 55 | + cpartition.reverse() |
| 56 | + for vb in vpartition: |
| 57 | + vb.reverse() |
| 58 | + for cb in cpartition: |
| 59 | + cb.reverse() |
| 60 | + |
| 61 | + return vpartition, cpartition |
| 62 | + |
| 63 | + |
| 64 | +def _get_rectangle_around_coords(ij1, ij2, linewidth=2, linestyle="-"): |
| 65 | + i1, j1 = ij1 |
| 66 | + i2, j2 = ij2 |
| 67 | + buffer = 0.5 |
| 68 | + ll_corner = (min(i1, i2) - buffer, min(j1, j2) - buffer) |
| 69 | + width = abs(i1 - i2) + 2 * buffer |
| 70 | + height = abs(j1 - j2) + 2 * buffer |
| 71 | + rect = matplotlib.patches.Rectangle( |
| 72 | + ll_corner, |
| 73 | + width, |
| 74 | + height, |
| 75 | + clip_on=False, |
| 76 | + fill=False, |
| 77 | + edgecolor="orange", |
| 78 | + linewidth=linewidth, |
| 79 | + linestyle=linestyle, |
| 80 | + ) |
| 81 | + return rect |
| 82 | + |
| 83 | + |
| 84 | +def spy_dulmage_mendelsohn( |
| 85 | + model, |
| 86 | + *, |
| 87 | + incidence_kwds=None, |
| 88 | + order=IncidenceOrder.dulmage_mendelsohn_upper, |
| 89 | + highlight_coarse=True, |
| 90 | + highlight_fine=True, |
| 91 | + skip_wellconstrained=False, |
| 92 | + ax=None, |
| 93 | + linewidth=2, |
| 94 | + spy_kwds=None, |
| 95 | +): |
| 96 | + """Plot sparsity structure in Dulmage-Mendelsohn order on Matplotlib axes |
| 97 | +
|
| 98 | + This is a wrapper around the Matplotlib ``Axes.spy`` method for plotting |
| 99 | + an incidence matrix in Dulmage-Mendelsohn order, with coarse and/or fine |
| 100 | + partitions highlighted. The coarse partition refers to the under-constrained, |
| 101 | + over-constrained, and well-constrained subsystems, while the fine partition |
| 102 | + refers to block diagonal or block triangular partitions of the former |
| 103 | + subsystems. |
| 104 | +
|
| 105 | + Parameters |
| 106 | + ---------- |
| 107 | +
|
| 108 | + model: ``ConcreteModel`` |
| 109 | + Input model to plot sparsity structure of |
| 110 | +
|
| 111 | + incidence_kwds: dict, optional |
| 112 | + Config options for ``IncidenceGraphInterface`` |
| 113 | +
|
| 114 | + order: ``IncidenceOrder``, optional |
| 115 | + Order in which to plot sparsity structure. Default is |
| 116 | + ``IncidenceOrder.dulmage_mendelsohn_upper`` for a block-upper triangular |
| 117 | + matrix. Set to ``IncidenceOrder.dulmage_mendelsohn_lower`` for a |
| 118 | + block-lower triangular matrix. |
| 119 | +
|
| 120 | + highlight_coarse: bool, optional |
| 121 | + Whether to draw a rectangle around the coarse partition. Default True |
| 122 | +
|
| 123 | + highlight_fine: bool, optional |
| 124 | + Whether to draw a rectangle around the fine partition. Default True |
| 125 | +
|
| 126 | + skip_wellconstrained: bool, optional |
| 127 | + Whether to skip highlighting the well-constrained subsystem of the |
| 128 | + coarse partition. Default False |
| 129 | +
|
| 130 | + ax: ``matplotlib.pyplot.Axes``, optional |
| 131 | + Axes object on which to plot. If not provided, new figure |
| 132 | + and axes are created. |
| 133 | +
|
| 134 | + linewidth: int, optional |
| 135 | + Line width of for rectangle used to highlight. Default 2 |
| 136 | +
|
| 137 | + spy_kwds: dict, optional |
| 138 | + Keyword arguments for ``Axes.spy`` |
| 139 | +
|
| 140 | + Returns |
| 141 | + ------- |
| 142 | +
|
| 143 | + fig: ``matplotlib.pyplot.Figure`` or ``None`` |
| 144 | + Figure on which the sparsity structure is plotted. ``None`` if axes |
| 145 | + are provided |
| 146 | +
|
| 147 | + ax: ``matplotlib.pyplot.Axes`` |
| 148 | + Axes on which the sparsity structure is plotted |
| 149 | +
|
| 150 | + """ |
| 151 | + plt = matplotlib.pyplot |
| 152 | + if incidence_kwds is None: |
| 153 | + incidence_kwds = {} |
| 154 | + if spy_kwds is None: |
| 155 | + spy_kwds = {} |
| 156 | + |
| 157 | + vpart, cpart = _partition_variables_and_constraints(model, order=order) |
| 158 | + vpart_fine = sum(vpart, []) |
| 159 | + cpart_fine = sum(cpart, []) |
| 160 | + vorder = sum(vpart_fine, []) |
| 161 | + corder = sum(cpart_fine, []) |
| 162 | + |
| 163 | + imat = get_structural_incidence_matrix(vorder, corder) |
| 164 | + nvar = len(vorder) |
| 165 | + ncon = len(corder) |
| 166 | + |
| 167 | + if ax is None: |
| 168 | + fig, ax = plt.subplots() |
| 169 | + else: |
| 170 | + fig = None |
| 171 | + |
| 172 | + markersize = spy_kwds.pop("markersize", None) |
| 173 | + if markersize is None: |
| 174 | + # At 10000 vars/cons, we want markersize=0.2 |
| 175 | + # At 20 vars/cons, we want markersize=10 |
| 176 | + # We assume we want a linear relationship between 1/nvar |
| 177 | + # and the markersize. |
| 178 | + markersize = (10.0 - 0.2) / (1 / 20 - 1 / 10000) * ( |
| 179 | + 1 / max(nvar, ncon) - 1 / 10000 |
| 180 | + ) + 0.2 |
| 181 | + |
| 182 | + ax.spy(imat, markersize=markersize, **spy_kwds) |
| 183 | + ax.tick_params(length=0) |
| 184 | + if highlight_coarse: |
| 185 | + start = (0, 0) |
| 186 | + for i, (vblocks, cblocks) in enumerate(zip(vpart, cpart)): |
| 187 | + # Get the total number of variables/constraints in this part |
| 188 | + # of the coarse partition |
| 189 | + nv = sum(len(vb) for vb in vblocks) |
| 190 | + nc = sum(len(cb) for cb in cblocks) |
| 191 | + stop = (start[0] + nv - 1, start[1] + nc - 1) |
| 192 | + if not (i == 1 and skip_wellconstrained) and nv > 0 and nc > 0: |
| 193 | + # Regardless of whether we are plotting in upper or lower |
| 194 | + # triangular order, the well-constrained subsystem is at |
| 195 | + # position 1 |
| 196 | + # |
| 197 | + # The get-rectangle function doesn't look good if we give it |
| 198 | + # an "empty region" to box. |
| 199 | + ax.add_patch( |
| 200 | + _get_rectangle_around_coords(start, stop, linewidth=linewidth) |
| 201 | + ) |
| 202 | + start = (stop[0] + 1, stop[1] + 1) |
| 203 | + |
| 204 | + if highlight_fine: |
| 205 | + # Use dashed lines to distinguish inner from outer partitions |
| 206 | + # if we are highlighting both |
| 207 | + linestyle = "--" if highlight_coarse else "-" |
| 208 | + start = (0, 0) |
| 209 | + for vb, cb in zip(vpart_fine, cpart_fine): |
| 210 | + stop = (start[0] + len(vb) - 1, start[1] + len(cb) - 1) |
| 211 | + # Note that the subset's we're boxing here can't be empty. |
| 212 | + ax.add_patch( |
| 213 | + _get_rectangle_around_coords( |
| 214 | + start, stop, linestyle=linestyle, linewidth=linewidth |
| 215 | + ) |
| 216 | + ) |
| 217 | + start = (stop[0] + 1, stop[1] + 1) |
| 218 | + |
| 219 | + return fig, ax |
0 commit comments