diff --git a/.github/workflows/typos.toml b/.github/workflows/typos.toml index be31d21a641..7a38164898b 100644 --- a/.github/workflows/typos.toml +++ b/.github/workflows/typos.toml @@ -43,4 +43,28 @@ Sur = "Sur" # contrib package named mis and the acronym whence the name comes mis = "mis" MIS = "MIS" +# Ignore the shorthand ans for answer +ans = "ans" +# Ignore the keyword arange +arange = "arange" +# Ignore IIS +IIS = "IIS" +iis = "iis" +# Ignore PN +PN = "PN" +# Ignore hd +hd = "hd" +# Ignore opf +opf = "opf" +# Ignore FRE +FRE = "FRE" +# Ignore MCH +MCH = "MCH" +# Ignore RO +ro = "ro" +RO = "RO" +# Ignore EOF - end of file +EOF = "EOF" +# Ignore lst as shorthand for list +lst = "lst" # AS NEEDED: Add More Words Below diff --git a/README.md b/README.md index 12c3ce8ed9a..707f1a06c5a 100644 --- a/README.md +++ b/README.md @@ -71,6 +71,7 @@ version, we will remove testing for that Python version. ### Tutorials and Examples +* [Pyomo — Optimization Modeling in Python](https://link.springer.com/book/10.1007/978-3-030-68928-5) * [Pyomo Workshop Slides](https://github.com/Pyomo/pyomo-tutorials/blob/main/Pyomo-Workshop-December-2023.pdf) * [Prof. Jeffrey Kantor's Pyomo Cookbook](https://jckantor.github.io/ND-Pyomo-Cookbook/) * The [companion notebooks](https://mobook.github.io/MO-book/intro.html) diff --git a/doc/OnlineDocs/bibliography.rst b/doc/OnlineDocs/bibliography.rst index 6cbb96d3bfb..c12d3f81d8c 100644 --- a/doc/OnlineDocs/bibliography.rst +++ b/doc/OnlineDocs/bibliography.rst @@ -39,6 +39,8 @@ Bibliography John D. Siirola, Jean-Paul Watson, and David L. Woodruff. Pyomo - Optimization Modeling in Python, 3rd Edition. Vol. 67. Springer, 2021. + doi: `10.1007/978-3-030-68928-5 + `_ .. [PyomoJournal] William E. Hart, Jean-Paul Watson, David L. Woodruff. "Pyomo: modeling and solving mathematical programs in diff --git a/doc/OnlineDocs/tutorial_examples.rst b/doc/OnlineDocs/tutorial_examples.rst index 6a40949ef90..a18f9d77d42 100644 --- a/doc/OnlineDocs/tutorial_examples.rst +++ b/doc/OnlineDocs/tutorial_examples.rst @@ -3,15 +3,18 @@ Pyomo Tutorial Examples Additional Pyomo tutorials and examples can be found at the following links: -`Pyomo Workshop Slides and Exercises -`_ +* `Pyomo — Optimization Modeling in Python + `_ ([PyomoBookIII]_) -`Prof. Jeffrey Kantor's Pyomo Cookbook -`_ +* `Pyomo Workshop Slides and Exercises + `_ -The `companion notebooks `_ -for *Hands-On Mathematical Optimization with Python* +* `Prof. Jeffrey Kantor's Pyomo Cookbook + `_ -`Pyomo Gallery `_ +* The `companion notebooks `_ + for *Hands-On Mathematical Optimization with Python* + +* `Pyomo Gallery `_ diff --git a/pyomo/contrib/gdpopt/tests/test_LBB.py b/pyomo/contrib/gdpopt/tests/test_LBB.py index 273327b02a4..8a553398fa6 100644 --- a/pyomo/contrib/gdpopt/tests/test_LBB.py +++ b/pyomo/contrib/gdpopt/tests/test_LBB.py @@ -59,6 +59,7 @@ def test_infeasible_GDP(self): self.assertIsNone(m.d.disjuncts[0].indicator_var.value) self.assertIsNone(m.d.disjuncts[1].indicator_var.value) + @unittest.skipUnless(z3_available, "Z3 SAT solver is not available") def test_infeasible_GDP_check_sat(self): """Test for infeasible GDP with check_sat option True.""" m = ConcreteModel() diff --git a/pyomo/contrib/gdpopt/tests/test_gdpopt.py b/pyomo/contrib/gdpopt/tests/test_gdpopt.py index 005df56ced5..873bafabc76 100644 --- a/pyomo/contrib/gdpopt/tests/test_gdpopt.py +++ b/pyomo/contrib/gdpopt/tests/test_gdpopt.py @@ -22,7 +22,6 @@ from pyomo.common.collections import Bunch from pyomo.common.config import ConfigDict, ConfigValue from pyomo.common.fileutils import import_file, PYOMO_ROOT_DIR -from pyomo.contrib.appsi.solvers.gurobi import Gurobi from pyomo.contrib.gdpopt.create_oa_subproblems import ( add_util_block, add_disjunct_list, @@ -767,6 +766,9 @@ def test_time_limit(self): results.solver.termination_condition, TerminationCondition.maxTimeLimit ) + @unittest.skipUnless( + license_available, "No BARON license--8PP logical problem exceeds demo size" + ) def test_LOA_8PP_logical_default_init(self): """Test logic-based outer approximation with 8PP.""" exfile = import_file(join(exdir, 'eight_process', 'eight_proc_logical.py')) @@ -870,6 +872,9 @@ def test_LOA_8PP_maxBinary(self): ) ct.check_8PP_solution(self, eight_process, results) + @unittest.skipUnless( + license_available, "No BARON license--8PP logical problem exceeds demo size" + ) def test_LOA_8PP_logical_maxBinary(self): """Test logic-based OA with max_binary initialization.""" exfile = import_file(join(exdir, 'eight_process', 'eight_proc_logical.py')) @@ -1050,7 +1055,11 @@ def assert_correct_disjuncts_active( self.assertTrue(fabs(value(eight_process.profit.expr) - 68) <= 1e-2) - @unittest.skipUnless(Gurobi().available(), "APPSI Gurobi solver is not available") + @unittest.skipUnless( + SolverFactory('appsi_gurobi').available(exception_flag=False) + and SolverFactory('appsi_gurobi').license_is_valid(), + "Legacy APPSI Gurobi solver is not available", + ) def test_auto_persistent_solver(self): exfile = import_file(join(exdir, 'eight_process', 'eight_proc_model.py')) m = exfile.build_eight_process_flowsheet() @@ -1126,6 +1135,9 @@ def test_RIC_8PP_default_init(self): ) ct.check_8PP_solution(self, eight_process, results) + @unittest.skipUnless( + license_available, "No BARON license--8PP logical problem exceeds demo size" + ) def test_RIC_8PP_logical_default_init(self): """Test logic-based outer approximation with 8PP.""" exfile = import_file(join(exdir, 'eight_process', 'eight_proc_logical.py')) diff --git a/pyomo/contrib/gdpopt/util.py b/pyomo/contrib/gdpopt/util.py index 2cb70f0ea60..babe0245d57 100644 --- a/pyomo/contrib/gdpopt/util.py +++ b/pyomo/contrib/gdpopt/util.py @@ -553,6 +553,13 @@ def _add_bigm_constraint_to_transformed_model(m, constraint, block): # making a Reference to the ComponentData so that it will look like an # indexed component for now. If I redesign bigm at some point, then this # could be prettier. - bigm._transform_constraint(Reference(constraint), parent_disjunct, None, [], []) + bigm._transform_constraint( + Reference(constraint), + parent_disjunct, + None, + [], + [], + 1 - parent_disjunct.binary_indicator_var, + ) # Now get rid of it because this is a class attribute! del bigm._config diff --git a/pyomo/contrib/incidence_analysis/config.py b/pyomo/contrib/incidence_analysis/config.py index 128273b4dec..2a7734ba433 100644 --- a/pyomo/contrib/incidence_analysis/config.py +++ b/pyomo/contrib/incidence_analysis/config.py @@ -36,6 +36,13 @@ class IncidenceMethod(enum.Enum): """Use ``pyomo.repn.plugins.nl_writer.AMPLRepnVisitor``""" +class IncidenceOrder(enum.Enum): + + dulmage_mendelsohn_upper = 0 + + dulmage_mendelsohn_lower = 1 + + _include_fixed = ConfigValue( default=False, domain=bool, diff --git a/pyomo/contrib/incidence_analysis/tests/test_visualize.py b/pyomo/contrib/incidence_analysis/tests/test_visualize.py new file mode 100644 index 00000000000..7c5538b671f --- /dev/null +++ b/pyomo/contrib/incidence_analysis/tests/test_visualize.py @@ -0,0 +1,47 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2024 +# National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and +# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain +# rights in this software. +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ + +import pyomo.common.unittest as unittest +from pyomo.common.dependencies import ( + matplotlib, + matplotlib_available, + scipy_available, + networkx_available, +) +from pyomo.contrib.incidence_analysis.visualize import spy_dulmage_mendelsohn +from pyomo.contrib.incidence_analysis.tests.models_for_testing import ( + make_gas_expansion_model, + make_dynamic_model, + make_degenerate_solid_phase_model, +) + + +@unittest.skipUnless(matplotlib_available, "Matplotlib is not available") +@unittest.skipUnless(scipy_available, "SciPy is not available") +@unittest.skipUnless(networkx_available, "NetworkX is not available") +class TestSpy(unittest.TestCase): + def test_spy_dulmage_mendelsohn(self): + models = [ + make_gas_expansion_model(), + make_dynamic_model(), + make_degenerate_solid_phase_model(), + ] + for m in models: + fig, ax = spy_dulmage_mendelsohn(m) + # Note that this is a weak test. We just test that we can call the + # plot method, it doesn't raise an error, and gives us back the + # types we expect. We don't attempt to validate the resulting plot. + self.assertTrue(isinstance(fig, matplotlib.pyplot.Figure)) + self.assertTrue(isinstance(ax, matplotlib.pyplot.Axes)) + + +if __name__ == "__main__": + unittest.main() diff --git a/pyomo/contrib/incidence_analysis/visualize.py b/pyomo/contrib/incidence_analysis/visualize.py new file mode 100644 index 00000000000..af1bdbbb918 --- /dev/null +++ b/pyomo/contrib/incidence_analysis/visualize.py @@ -0,0 +1,219 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2024 +# National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and +# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain +# rights in this software. +# This software is distributed under the 3-clause BSD License. +# ___________________________________________________________________________ +"""Module for visualizing results of incidence graph or matrix analysis + +""" +from pyomo.contrib.incidence_analysis.config import IncidenceOrder +from pyomo.contrib.incidence_analysis.interface import ( + IncidenceGraphInterface, + get_structural_incidence_matrix, +) +from pyomo.common.dependencies import matplotlib + + +def _partition_variables_and_constraints( + model, order=IncidenceOrder.dulmage_mendelsohn_upper, **kwds +): + """Partition variables and constraints in an incidence graph""" + igraph = IncidenceGraphInterface(model, **kwds) + vdmp, cdmp = igraph.dulmage_mendelsohn() + + ucv = vdmp.unmatched + vdmp.underconstrained + ucc = cdmp.underconstrained + + ocv = vdmp.overconstrained + occ = cdmp.overconstrained + cdmp.unmatched + + ucvblocks, uccblocks = igraph.get_connected_components( + variables=ucv, constraints=ucc + ) + ocvblocks, occblocks = igraph.get_connected_components( + variables=ocv, constraints=occ + ) + wcvblocks, wccblocks = igraph.block_triangularize( + variables=vdmp.square, constraints=cdmp.square + ) + # By default, we block-*lower* triangularize. By default, however, we want + # the Dulmage-Mendelsohn decomposition to be block-*upper* triangular. + wcvblocks.reverse() + wccblocks.reverse() + vpartition = [ucvblocks, wcvblocks, ocvblocks] + cpartition = [uccblocks, wccblocks, occblocks] + + if order == IncidenceOrder.dulmage_mendelsohn_lower: + # If a block-lower triangular matrix was requested, we need to reverse + # both the inner and outer partitions + vpartition.reverse() + cpartition.reverse() + for vb in vpartition: + vb.reverse() + for cb in cpartition: + cb.reverse() + + return vpartition, cpartition + + +def _get_rectangle_around_coords(ij1, ij2, linewidth=2, linestyle="-"): + i1, j1 = ij1 + i2, j2 = ij2 + buffer = 0.5 + ll_corner = (min(i1, i2) - buffer, min(j1, j2) - buffer) + width = abs(i1 - i2) + 2 * buffer + height = abs(j1 - j2) + 2 * buffer + rect = matplotlib.patches.Rectangle( + ll_corner, + width, + height, + clip_on=False, + fill=False, + edgecolor="orange", + linewidth=linewidth, + linestyle=linestyle, + ) + return rect + + +def spy_dulmage_mendelsohn( + model, + *, + incidence_kwds=None, + order=IncidenceOrder.dulmage_mendelsohn_upper, + highlight_coarse=True, + highlight_fine=True, + skip_wellconstrained=False, + ax=None, + linewidth=2, + spy_kwds=None, +): + """Plot sparsity structure in Dulmage-Mendelsohn order on Matplotlib axes + + This is a wrapper around the Matplotlib ``Axes.spy`` method for plotting + an incidence matrix in Dulmage-Mendelsohn order, with coarse and/or fine + partitions highlighted. The coarse partition refers to the under-constrained, + over-constrained, and well-constrained subsystems, while the fine partition + refers to block diagonal or block triangular partitions of the former + subsystems. + + Parameters + ---------- + + model: ``ConcreteModel`` + Input model to plot sparsity structure of + + incidence_kwds: dict, optional + Config options for ``IncidenceGraphInterface`` + + order: ``IncidenceOrder``, optional + Order in which to plot sparsity structure. Default is + ``IncidenceOrder.dulmage_mendelsohn_upper`` for a block-upper triangular + matrix. Set to ``IncidenceOrder.dulmage_mendelsohn_lower`` for a + block-lower triangular matrix. + + highlight_coarse: bool, optional + Whether to draw a rectangle around the coarse partition. Default True + + highlight_fine: bool, optional + Whether to draw a rectangle around the fine partition. Default True + + skip_wellconstrained: bool, optional + Whether to skip highlighting the well-constrained subsystem of the + coarse partition. Default False + + ax: ``matplotlib.pyplot.Axes``, optional + Axes object on which to plot. If not provided, new figure + and axes are created. + + linewidth: int, optional + Line width of for rectangle used to highlight. Default 2 + + spy_kwds: dict, optional + Keyword arguments for ``Axes.spy`` + + Returns + ------- + + fig: ``matplotlib.pyplot.Figure`` or ``None`` + Figure on which the sparsity structure is plotted. ``None`` if axes + are provided + + ax: ``matplotlib.pyplot.Axes`` + Axes on which the sparsity structure is plotted + + """ + plt = matplotlib.pyplot + if incidence_kwds is None: + incidence_kwds = {} + if spy_kwds is None: + spy_kwds = {} + + vpart, cpart = _partition_variables_and_constraints(model, order=order) + vpart_fine = sum(vpart, []) + cpart_fine = sum(cpart, []) + vorder = sum(vpart_fine, []) + corder = sum(cpart_fine, []) + + imat = get_structural_incidence_matrix(vorder, corder) + nvar = len(vorder) + ncon = len(corder) + + if ax is None: + fig, ax = plt.subplots() + else: + fig = None + + markersize = spy_kwds.pop("markersize", None) + if markersize is None: + # At 10000 vars/cons, we want markersize=0.2 + # At 20 vars/cons, we want markersize=10 + # We assume we want a linear relationship between 1/nvar + # and the markersize. + markersize = (10.0 - 0.2) / (1 / 20 - 1 / 10000) * ( + 1 / max(nvar, ncon) - 1 / 10000 + ) + 0.2 + + ax.spy(imat, markersize=markersize, **spy_kwds) + ax.tick_params(length=0) + if highlight_coarse: + start = (0, 0) + for i, (vblocks, cblocks) in enumerate(zip(vpart, cpart)): + # Get the total number of variables/constraints in this part + # of the coarse partition + nv = sum(len(vb) for vb in vblocks) + nc = sum(len(cb) for cb in cblocks) + stop = (start[0] + nv - 1, start[1] + nc - 1) + if not (i == 1 and skip_wellconstrained) and nv > 0 and nc > 0: + # Regardless of whether we are plotting in upper or lower + # triangular order, the well-constrained subsystem is at + # position 1 + # + # The get-rectangle function doesn't look good if we give it + # an "empty region" to box. + ax.add_patch( + _get_rectangle_around_coords(start, stop, linewidth=linewidth) + ) + start = (stop[0] + 1, stop[1] + 1) + + if highlight_fine: + # Use dashed lines to distinguish inner from outer partitions + # if we are highlighting both + linestyle = "--" if highlight_coarse else "-" + start = (0, 0) + for vb, cb in zip(vpart_fine, cpart_fine): + stop = (start[0] + len(vb) - 1, start[1] + len(cb) - 1) + # Note that the subset's we're boxing here can't be empty. + ax.add_patch( + _get_rectangle_around_coords( + start, stop, linestyle=linestyle, linewidth=linewidth + ) + ) + start = (stop[0] + 1, stop[1] + 1) + + return fig, ax diff --git a/pyomo/contrib/latex_printer/__init__.py b/pyomo/contrib/latex_printer/__init__.py index c434b53dfe1..02eaa636a36 100644 --- a/pyomo/contrib/latex_printer/__init__.py +++ b/pyomo/contrib/latex_printer/__init__.py @@ -9,22 +9,11 @@ # This software is distributed under the 3-clause BSD License. # ___________________________________________________________________________ -# ___________________________________________________________________________ -# -# Pyomo: Python Optimization Modeling Objects -# Copyright (c) 2008-2023 -# National Technology and Engineering Solutions of Sandia, LLC -# Under the terms of Contract DE-NA0003525 with National Technology and -# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain -# rights in this software. -# This software is distributed under the 3-clause BSD License. -# ___________________________________________________________________________ - # Recommended just to build all of the appropriate things import pyomo.environ # Remove one layer of .latex_printer -# import statemnt is now: +# import statement is now: # from pyomo.contrib.latex_printer import latex_printer try: from pyomo.contrib.latex_printer.latex_printer import latex_printer diff --git a/pyomo/contrib/pynumero/algorithms/solvers/cyipopt_solver.py b/pyomo/contrib/pynumero/algorithms/solvers/cyipopt_solver.py index cdea542295b..53616298415 100644 --- a/pyomo/contrib/pynumero/algorithms/solvers/cyipopt_solver.py +++ b/pyomo/contrib/pynumero/algorithms/solvers/cyipopt_solver.py @@ -319,7 +319,13 @@ def license_is_valid(self): return True def version(self): - return tuple(int(_) for _ in cyipopt.__version__.split(".")) + def _int(x): + try: + return int(x) + except: + return x + + return tuple(_int(_) for _ in cyipopt_interface.cyipopt.__version__.split(".")) def solve(self, model, **kwds): config = self.config(kwds, preserve_implicit=True) diff --git a/pyomo/contrib/pynumero/examples/tests/test_cyipopt_examples.py b/pyomo/contrib/pynumero/examples/tests/test_cyipopt_examples.py index 408a0197382..2df43c1e797 100644 --- a/pyomo/contrib/pynumero/examples/tests/test_cyipopt_examples.py +++ b/pyomo/contrib/pynumero/examples/tests/test_cyipopt_examples.py @@ -44,11 +44,13 @@ raise unittest.SkipTest("Pynumero needs the ASL extension to run CyIpopt tests") import pyomo.contrib.pynumero.algorithms.solvers.cyipopt_solver as cyipopt_solver +from pyomo.contrib.pynumero.interfaces.cyipopt_interface import cyipopt_available -if not cyipopt_solver.cyipopt_available: +if not cyipopt_available: raise unittest.SkipTest("PyNumero needs CyIpopt installed to run CyIpopt tests") import cyipopt as cyipopt_core + example_dir = os.path.join(this_file_dir(), '..') @@ -266,6 +268,11 @@ def test_cyipopt_functor(self): s = df['ca_bal'] self.assertAlmostEqual(s.iloc[6], 0, places=3) + @unittest.skipIf( + cyipopt_solver.PyomoCyIpoptSolver().version() == (1, 4, 0), + "Terminating Ipopt through a user callback is broken in CyIpopt 1.4.0 " + "(see mechmotum/cyipopt#249)", + ) def test_cyipopt_callback_halt(self): ex = import_file( os.path.join(example_dir, 'callback', 'cyipopt_callback_halt.py') diff --git a/pyomo/core/base/component.py b/pyomo/core/base/component.py index 22c2bc4b804..9b1929daa06 100644 --- a/pyomo/core/base/component.py +++ b/pyomo/core/base/component.py @@ -368,7 +368,7 @@ def pprint(self, ostream=None, verbose=False, prefix=""): @property def name(self): - """Get the fully qualifed component name.""" + """Get the fully qualified component name.""" return self.getname(fully_qualified=True) # Adding a setter here to help users adapt to the new @@ -664,7 +664,7 @@ def getname(self, fully_qualified=False, name_buffer=None, relative_to=None): @property def name(self): - """Get the fully qualifed component name.""" + """Get the fully qualified component name.""" return self.getname(fully_qualified=True) # Allow setting a component's name if it is not owned by a parent diff --git a/pyomo/core/base/indexed_component.py b/pyomo/core/base/indexed_component.py index e1be613d666..37a62e5c4d7 100644 --- a/pyomo/core/base/indexed_component.py +++ b/pyomo/core/base/indexed_component.py @@ -731,7 +731,7 @@ def __delitem__(self, index): # this supports "del m.x[:,1]" through a simple recursive call if index.__class__ is IndexedComponent_slice: - # Assert that this slice ws just generated + # Assert that this slice was just generated assert len(index._call_stack) == 1 # Make a copy of the slicer items *before* we start # iterating over it (since we will be removing items!). diff --git a/pyomo/core/base/reference.py b/pyomo/core/base/reference.py index 2279db067a6..84cccec9749 100644 --- a/pyomo/core/base/reference.py +++ b/pyomo/core/base/reference.py @@ -579,7 +579,7 @@ def Reference(reference, ctype=NOTSET): :py:class:`IndexedComponent`. If the indices associated with wildcards in the component slice all - refer to the same :py:class:`Set` objects for all data identifed by + refer to the same :py:class:`Set` objects for all data identified by the slice, then the resulting indexed component will be indexed by the product of those sets. However, if all data do not share common set objects, or only a subset of indices in a multidimentional set diff --git a/pyomo/core/expr/base.py b/pyomo/core/expr/base.py index f506956e478..6e2066afcc5 100644 --- a/pyomo/core/expr/base.py +++ b/pyomo/core/expr/base.py @@ -360,7 +360,7 @@ def size(self): """ return visitor.sizeof_expression(self) - def _apply_operation(self, result): # pragma: no cover + def _apply_operation(self, result): """ Compute the values of this node given the values of its children. diff --git a/pyomo/core/expr/numeric_expr.py b/pyomo/core/expr/numeric_expr.py index 25d83ca20f4..9b624d2b8bd 100644 --- a/pyomo/core/expr/numeric_expr.py +++ b/pyomo/core/expr/numeric_expr.py @@ -2288,8 +2288,11 @@ def _iadd_mutablenpvsum_mutable(a, b): def _iadd_mutablenpvsum_native(a, b): if not b: return a - a._args_.append(b) - a._nargs += 1 + if a._args_ and a._args_[-1].__class__ in native_numeric_types: + a._args_[-1] += b + else: + a._args_.append(b) + a._nargs += 1 return a @@ -2301,9 +2304,7 @@ def _iadd_mutablenpvsum_npv(a, b): def _iadd_mutablenpvsum_param(a, b): if b.is_constant(): - b = b.value - if not b: - return a + return _iadd_mutablesum_native(a, b.value) a._args_.append(b) a._nargs += 1 return a @@ -2384,8 +2385,11 @@ def _iadd_mutablelinear_mutable(a, b): def _iadd_mutablelinear_native(a, b): if not b: return a - a._args_.append(b) - a._nargs += 1 + if a._args_ and a._args_[-1].__class__ in native_numeric_types: + a._args_[-1] += b + else: + a._args_.append(b) + a._nargs += 1 return a @@ -2397,9 +2401,7 @@ def _iadd_mutablelinear_npv(a, b): def _iadd_mutablelinear_param(a, b): if b.is_constant(): - b = b.value - if not b: - return a + return _iadd_mutablesum_native(a, b.value) a._args_.append(b) a._nargs += 1 return a @@ -2483,8 +2485,11 @@ def _iadd_mutablesum_mutable(a, b): def _iadd_mutablesum_native(a, b): if not b: return a - a._args_.append(b) - a._nargs += 1 + if a._args_ and a._args_[-1].__class__ in native_numeric_types: + a._args_[-1] += b + else: + a._args_.append(b) + a._nargs += 1 return a @@ -2496,9 +2501,7 @@ def _iadd_mutablesum_npv(a, b): def _iadd_mutablesum_param(a, b): if b.is_constant(): - b = b.value - if not b: - return a + return _iadd_mutablesum_native(a, b.value) a._args_.append(b) a._nargs += 1 return a diff --git a/pyomo/core/expr/template_expr.py b/pyomo/core/expr/template_expr.py index f65a1f2b9b0..d30046e9d82 100644 --- a/pyomo/core/expr/template_expr.py +++ b/pyomo/core/expr/template_expr.py @@ -19,11 +19,12 @@ from pyomo.core.expr.base import ExpressionBase, ExpressionArgs_Mixin, NPV_Mixin from pyomo.core.expr.logical_expr import BooleanExpression from pyomo.core.expr.numeric_expr import ( + ARG_TYPE, NumericExpression, - SumExpression, Numeric_NPV_Mixin, + SumExpression, + mutable_expression, register_arg_type, - ARG_TYPE, _balanced_parens, ) from pyomo.core.expr.numvalue import ( @@ -116,18 +117,10 @@ def _to_string(self, values, verbose, smap): return "%s[%s]" % (values[0], ','.join(values[1:])) def _resolve_template(self, args): - return args[0].__getitem__(tuple(args[1:])) + return args[0].__getitem__(args[1:]) def _apply_operation(self, result): - args = tuple( - ( - arg - if arg.__class__ in native_types or not arg.is_numeric_type() - else value(arg) - ) - for arg in result[1:] - ) - return result[0].__getitem__(tuple(result[1:])) + return result[0].__getitem__(result[1:]) class Numeric_GetItemExpression(GetItemExpression, NumericExpression): @@ -258,8 +251,8 @@ def nargs(self): return 2 def _apply_operation(self, result): - assert len(result) == 2 - return getattr(result[0], result[1]) + obj, attr = result + return getattr(obj, attr) def _to_string(self, values, verbose, smap): assert len(values) == 2 @@ -273,7 +266,7 @@ def _to_string(self, values, verbose, smap): return "%s.%s" % (values[0], attr) def _resolve_template(self, args): - return getattr(*tuple(args)) + return getattr(*args) class Numeric_GetAttrExpression(GetAttrExpression, NumericExpression): @@ -521,7 +514,15 @@ def _to_string(self, values, verbose, smap): return 'SUM(%s %s)' % (val, iterStr) def _resolve_template(self, args): - return SumExpression(args) + with mutable_expression() as e: + for arg in args: + e += arg + if e.nargs() > 1: + return e + elif not e.nargs(): + return 0 + else: + return e.arg(0) class IndexTemplate(NumericValue): diff --git a/pyomo/core/tests/unit/test_numeric_expr_dispatcher.py b/pyomo/core/tests/unit/test_numeric_expr_dispatcher.py index 37833d7e8a4..bb7a291e67d 100644 --- a/pyomo/core/tests/unit/test_numeric_expr_dispatcher.py +++ b/pyomo/core/tests/unit/test_numeric_expr_dispatcher.py @@ -6490,11 +6490,11 @@ def test_mutable_nvp_iadd(self): (mutable_npv, self.invalid, NotImplemented), (mutable_npv, self.asbinary, _MutableLinearExpression([10, self.bin])), (mutable_npv, self.zero, _MutableNPVSumExpression([10])), - (mutable_npv, self.one, _MutableNPVSumExpression([10, 1])), + (mutable_npv, self.one, _MutableNPVSumExpression([11])), # 4: - (mutable_npv, self.native, _MutableNPVSumExpression([10, 5])), + (mutable_npv, self.native, _MutableNPVSumExpression([15])), (mutable_npv, self.npv, _MutableNPVSumExpression([10, self.npv])), - (mutable_npv, self.param, _MutableNPVSumExpression([10, 6])), + (mutable_npv, self.param, _MutableNPVSumExpression([16])), ( mutable_npv, self.param_mut, @@ -6534,7 +6534,7 @@ def test_mutable_nvp_iadd(self): _MutableSumExpression([10] + self.mutable_l2.args), ), (mutable_npv, self.param0, _MutableNPVSumExpression([10])), - (mutable_npv, self.param1, _MutableNPVSumExpression([10, 1])), + (mutable_npv, self.param1, _MutableNPVSumExpression([11])), # 20: (mutable_npv, self.mutable_l3, _MutableNPVSumExpression([10, self.npv])), ] diff --git a/pyomo/core/tests/unit/test_numeric_expr_zerofilter.py b/pyomo/core/tests/unit/test_numeric_expr_zerofilter.py index 8e75ccc3feb..19968640a21 100644 --- a/pyomo/core/tests/unit/test_numeric_expr_zerofilter.py +++ b/pyomo/core/tests/unit/test_numeric_expr_zerofilter.py @@ -6020,11 +6020,11 @@ def test_mutable_nvp_iadd(self): (mutable_npv, self.invalid, NotImplemented), (mutable_npv, self.asbinary, _MutableLinearExpression([10, self.bin])), (mutable_npv, self.zero, _MutableNPVSumExpression([10])), - (mutable_npv, self.one, _MutableNPVSumExpression([10, 1])), + (mutable_npv, self.one, _MutableNPVSumExpression([11])), # 4: - (mutable_npv, self.native, _MutableNPVSumExpression([10, 5])), + (mutable_npv, self.native, _MutableNPVSumExpression([15])), (mutable_npv, self.npv, _MutableNPVSumExpression([10, self.npv])), - (mutable_npv, self.param, _MutableNPVSumExpression([10, 6])), + (mutable_npv, self.param, _MutableNPVSumExpression([16])), ( mutable_npv, self.param_mut, @@ -6064,7 +6064,7 @@ def test_mutable_nvp_iadd(self): _MutableSumExpression([10] + self.mutable_l2.args), ), (mutable_npv, self.param0, _MutableNPVSumExpression([10])), - (mutable_npv, self.param1, _MutableNPVSumExpression([10, 1])), + (mutable_npv, self.param1, _MutableNPVSumExpression([11])), # 20: (mutable_npv, self.mutable_l3, _MutableNPVSumExpression([10, self.npv])), ] diff --git a/pyomo/core/tests/unit/test_template_expr.py b/pyomo/core/tests/unit/test_template_expr.py index 4f255e3567a..4c872e1e11d 100644 --- a/pyomo/core/tests/unit/test_template_expr.py +++ b/pyomo/core/tests/unit/test_template_expr.py @@ -490,14 +490,14 @@ def c(m): self.assertEqual( str(resolve_template(template)), 'x[1,1,10] + ' - '(x[2,1,10] + x[2,1,20]) + ' - '(x[3,1,10] + x[3,1,20] + x[3,1,30]) + ' - '(x[1,2,10]) + ' - '(x[2,2,10] + x[2,2,20]) + ' - '(x[3,2,10] + x[3,2,20] + x[3,2,30]) + ' - '(x[1,3,10]) + ' - '(x[2,3,10] + x[2,3,20]) + ' - '(x[3,3,10] + x[3,3,20] + x[3,3,30]) <= 0', + 'x[2,1,10] + x[2,1,20] + ' + 'x[3,1,10] + x[3,1,20] + x[3,1,30] + ' + 'x[1,2,10] + ' + 'x[2,2,10] + x[2,2,20] + ' + 'x[3,2,10] + x[3,2,20] + x[3,2,30] + ' + 'x[1,3,10] + ' + 'x[2,3,10] + x[2,3,20] + ' + 'x[3,3,10] + x[3,3,20] + x[3,3,30] <= 0', ) def test_multidim_nested_sum_rule(self): @@ -566,14 +566,14 @@ def c(m): self.assertEqual( str(resolve_template(template)), 'x[1,1,10] + ' - '(x[2,1,10] + x[2,1,20]) + ' - '(x[3,1,10] + x[3,1,20] + x[3,1,30]) + ' - '(x[1,2,10]) + ' - '(x[2,2,10] + x[2,2,20]) + ' - '(x[3,2,10] + x[3,2,20] + x[3,2,30]) + ' - '(x[1,3,10]) + ' - '(x[2,3,10] + x[2,3,20]) + ' - '(x[3,3,10] + x[3,3,20] + x[3,3,30]) <= 0', + 'x[2,1,10] + x[2,1,20] + ' + 'x[3,1,10] + x[3,1,20] + x[3,1,30] + ' + 'x[1,2,10] + ' + 'x[2,2,10] + x[2,2,20] + ' + 'x[3,2,10] + x[3,2,20] + x[3,2,30] + ' + 'x[1,3,10] + ' + 'x[2,3,10] + x[2,3,20] + ' + 'x[3,3,10] + x[3,3,20] + x[3,3,30] <= 0', ) def test_multidim_nested_getattr_sum_rule(self): @@ -609,14 +609,14 @@ def c(m): self.assertEqual( str(resolve_template(template)), 'x[1,1,10] + ' - '(x[2,1,10] + x[2,1,20]) + ' - '(x[3,1,10] + x[3,1,20] + x[3,1,30]) + ' - '(x[1,2,10]) + ' - '(x[2,2,10] + x[2,2,20]) + ' - '(x[3,2,10] + x[3,2,20] + x[3,2,30]) + ' - '(x[1,3,10]) + ' - '(x[2,3,10] + x[2,3,20]) + ' - '(x[3,3,10] + x[3,3,20] + x[3,3,30]) <= 0', + 'x[2,1,10] + x[2,1,20] + ' + 'x[3,1,10] + x[3,1,20] + x[3,1,30] + ' + 'x[1,2,10] + ' + 'x[2,2,10] + x[2,2,20] + ' + 'x[3,2,10] + x[3,2,20] + x[3,2,30] + ' + 'x[1,3,10] + ' + 'x[2,3,10] + x[2,3,20] + ' + 'x[3,3,10] + x[3,3,20] + x[3,3,30] <= 0', ) def test_eval_getattr(self): diff --git a/pyomo/gdp/plugins/bigm.py b/pyomo/gdp/plugins/bigm.py index 3f450dbbd4f..d715d913db8 100644 --- a/pyomo/gdp/plugins/bigm.py +++ b/pyomo/gdp/plugins/bigm.py @@ -213,21 +213,15 @@ def _apply_to_impl(self, instance, **kwds): bigM = self._config.bigM for t in preprocessed_targets: if t.ctype is Disjunction: - self._transform_disjunctionData( - t, - t.index(), - bigM, - parent_disjunct=gdp_tree.parent(t), - root_disjunct=gdp_tree.root_disjunct(t), - ) + self._transform_disjunctionData(t, t.index(), bigM, gdp_tree) # issue warnings about anything that was in the bigM args dict that we # didn't use _warn_for_unused_bigM_args(bigM, self.used_args, logger) - def _transform_disjunctionData( - self, obj, index, bigM, parent_disjunct=None, root_disjunct=None - ): + def _transform_disjunctionData(self, obj, index, bigM, gdp_tree): + parent_disjunct = gdp_tree.parent(obj) + root_disjunct = gdp_tree.root_disjunct(obj) (transBlock, xorConstraint) = self._setup_transform_disjunctionData( obj, root_disjunct ) @@ -236,7 +230,7 @@ def _transform_disjunctionData( or_expr = 0 for disjunct in obj.disjuncts: or_expr += disjunct.binary_indicator_var - self._transform_disjunct(disjunct, bigM, transBlock) + self._transform_disjunct(disjunct, bigM, transBlock, gdp_tree) if obj.xor: xorConstraint[index] = or_expr == 1 @@ -249,7 +243,7 @@ def _transform_disjunctionData( # and deactivate for the writers obj.deactivate() - def _transform_disjunct(self, obj, bigM, transBlock): + def _transform_disjunct(self, obj, bigM, transBlock, gdp_tree): # We're not using the preprocessed list here, so this could be # inactive. We've already done the error checking in preprocessing, so # we just skip it here. @@ -261,6 +255,12 @@ def _transform_disjunct(self, obj, bigM, transBlock): relaxationBlock = self._get_disjunct_transformation_block(obj, transBlock) + indicator_expression = 0 + node = obj + while node is not None: + indicator_expression += 1 - node.binary_indicator_var + node = gdp_tree.parent_disjunct(node) + # This is crazy, but if the disjunction has been previously # relaxed, the disjunct *could* be deactivated. This is a big # deal for Hull, as it uses the component_objects / @@ -270,13 +270,21 @@ def _transform_disjunct(self, obj, bigM, transBlock): # comparing the two relaxations. # # Transform each component within this disjunct - self._transform_block_components(obj, obj, bigM, arg_list, suffix_list) + self._transform_block_components( + obj, obj, bigM, arg_list, suffix_list, indicator_expression + ) # deactivate disjunct to keep the writers happy obj._deactivate_without_fixing_indicator() def _transform_constraint( - self, obj, disjunct, bigMargs, arg_list, disjunct_suffix_list + self, + obj, + disjunct, + bigMargs, + arg_list, + disjunct_suffix_list, + indicator_expression, ): # add constraint to the transformation block, we'll transform it there. transBlock = disjunct._transformation_block() @@ -348,7 +356,13 @@ def _transform_constraint( bigm_src[c] = (lower, upper) self._add_constraint_expressions( - c, i, M, disjunct.binary_indicator_var, newConstraint, constraint_map + c, + i, + M, + disjunct.binary_indicator_var, + newConstraint, + constraint_map, + indicator_expression=indicator_expression, ) # deactivate because we relaxed diff --git a/pyomo/gdp/plugins/bigm_mixin.py b/pyomo/gdp/plugins/bigm_mixin.py index 510b36b5102..1c3fcb2c64a 100644 --- a/pyomo/gdp/plugins/bigm_mixin.py +++ b/pyomo/gdp/plugins/bigm_mixin.py @@ -232,7 +232,14 @@ def _estimate_M(self, expr, constraint): return tuple(M) def _add_constraint_expressions( - self, c, i, M, indicator_var, newConstraint, constraint_map + self, + c, + i, + M, + indicator_var, + newConstraint, + constraint_map, + indicator_expression=None, ): # Since we are both combining components from multiple blocks and using # local names, we need to make sure that the first index for @@ -244,6 +251,8 @@ def _add_constraint_expressions( # over the constraint indices, but I don't think it matters a lot.) unique = len(newConstraint) name = c.local_name + "_%s" % unique + if indicator_expression is None: + indicator_expression = 1 - indicator_var if c.lower is not None: if M[0] is None: @@ -251,7 +260,7 @@ def _add_constraint_expressions( "Cannot relax disjunctive constraint '%s' " "because M is not defined." % name ) - M_expr = M[0] * (1 - indicator_var) + M_expr = M[0] * indicator_expression newConstraint.add((name, i, 'lb'), c.lower <= c.body - M_expr) constraint_map.transformed_constraints[c].append( newConstraint[name, i, 'lb'] @@ -263,7 +272,7 @@ def _add_constraint_expressions( "Cannot relax disjunctive constraint '%s' " "because M is not defined." % name ) - M_expr = M[1] * (1 - indicator_var) + M_expr = M[1] * indicator_expression newConstraint.add((name, i, 'ub'), c.body - M_expr <= c.upper) constraint_map.transformed_constraints[c].append( newConstraint[name, i, 'ub'] diff --git a/pyomo/gdp/plugins/fix_disjuncts.py b/pyomo/gdp/plugins/fix_disjuncts.py index 44a9d91d513..172363caab7 100644 --- a/pyomo/gdp/plugins/fix_disjuncts.py +++ b/pyomo/gdp/plugins/fix_disjuncts.py @@ -52,7 +52,7 @@ class GDP_Disjunct_Fixer(Transformation): This reclassifies all disjuncts in the passed model instance as ctype Block and deactivates the constraints and disjunctions within inactive disjuncts. - In addition, it transforms relvant LogicalConstraints and BooleanVars so + In addition, it transforms relevant LogicalConstraints and BooleanVars so that the resulting model is a (MI)(N)LP (where it is only mixed-integer if the model contains integer-domain Vars or BooleanVars which were not indicator_vars of Disjuncs. diff --git a/pyomo/gdp/tests/models.py b/pyomo/gdp/tests/models.py index 0b84641899c..2995cacb450 100644 --- a/pyomo/gdp/tests/models.py +++ b/pyomo/gdp/tests/models.py @@ -840,7 +840,7 @@ def makeAnyIndexedDisjunctionOfDisjunctDatas(): build from DisjunctDatas. Identical mathematically to makeDisjunctionOfDisjunctDatas. - Used to test that the right things happen for a case where soemone + Used to test that the right things happen for a case where someone implements an algorithm which iteratively generates disjuncts and retransforms""" m = ConcreteModel() diff --git a/pyomo/gdp/tests/test_bigm.py b/pyomo/gdp/tests/test_bigm.py index c6ac49f6d36..3174a95292e 100644 --- a/pyomo/gdp/tests/test_bigm.py +++ b/pyomo/gdp/tests/test_bigm.py @@ -19,8 +19,11 @@ Set, Constraint, ComponentMap, + LogicalConstraint, + Objective, SolverFactory, Suffix, + TerminationCondition, ConcreteModel, Var, Any, @@ -1879,12 +1882,11 @@ def test_m_value_mappings(self): # many of the transformed constraints look like this, so can call this # function to test them. def check_bigM_constraint(self, cons, variable, M, indicator_var): - repn = generate_standard_repn(cons.body) - self.assertTrue(repn.is_linear()) - self.assertEqual(repn.constant, -M) - self.assertEqual(len(repn.linear_vars), 2) - ct.check_linear_coef(self, repn, variable, 1) - ct.check_linear_coef(self, repn, indicator_var, M) + assertExpressionsEqual( + self, + cons.body, + variable - float(M) * (1 - indicator_var.get_associated_binary()), + ) def check_inner_xor_constraint(self, inner_disjunction, outer_disjunct, bigm): inner_xor = inner_disjunction.algebraic_constraint @@ -1949,6 +1951,10 @@ def test_transformed_constraints(self): .binary_indicator_var, ) ), + 1, + EXPR.MonomialTermExpression( + (-1, m.disjunct[1].binary_indicator_var) + ), ] ), ) @@ -1958,37 +1964,69 @@ def test_transformed_constraints(self): ] ), ) - self.assertIsNone(cons1ub.lower) - self.assertEqual(cons1ub.upper, 0) - self.check_bigM_constraint( - cons1ub, m.z, 10, m.disjunct[1].innerdisjunct[0].indicator_var + assertExpressionsEqual( + self, + cons1ub.expr, + m.z + - 10.0 + * ( + 1 + - m.disjunct[1].innerdisjunct[0].binary_indicator_var + + 1 + - m.disjunct[1].binary_indicator_var + ) + <= 0.0, ) cons2 = bigm.get_transformed_constraints(m.disjunct[1].innerdisjunct[1].c) self.assertEqual(len(cons2), 1) cons2lb = cons2[0] - self.assertEqual(cons2lb.lower, 5) - self.assertIsNone(cons2lb.upper) - self.check_bigM_constraint( - cons2lb, m.z, -5, m.disjunct[1].innerdisjunct[1].indicator_var + assertExpressionsEqual( + self, + cons2lb.expr, + 5.0 + <= m.z + - (-5.0) + * ( + 1 + - m.disjunct[1].innerdisjunct[1].binary_indicator_var + + 1 + - m.disjunct[1].binary_indicator_var + ), ) cons3 = bigm.get_transformed_constraints(m.simpledisjunct.innerdisjunct0.c) self.assertEqual(len(cons3), 1) cons3ub = cons3[0] - self.assertEqual(cons3ub.upper, 2) - self.assertIsNone(cons3ub.lower) - self.check_bigM_constraint( - cons3ub, m.x, 7, m.simpledisjunct.innerdisjunct0.indicator_var + assertExpressionsEqual( + self, + cons3ub.expr, + m.x + - 7.0 + * ( + 1 + - m.simpledisjunct.innerdisjunct0.binary_indicator_var + + 1 + - m.simpledisjunct.binary_indicator_var + ) + <= 2.0, ) cons4 = bigm.get_transformed_constraints(m.simpledisjunct.innerdisjunct1.c) self.assertEqual(len(cons4), 1) cons4lb = cons4[0] - self.assertEqual(cons4lb.lower, 4) - self.assertIsNone(cons4lb.upper) - self.check_bigM_constraint( - cons4lb, m.x, -13, m.simpledisjunct.innerdisjunct1.indicator_var + assertExpressionsEqual( + self, + cons4lb.expr, + m.x + - (-13.0) + * ( + 1 + - m.simpledisjunct.innerdisjunct1.binary_indicator_var + + 1 + - m.simpledisjunct.binary_indicator_var + ) + >= 4.0, ) # Here we check that the xor constraint from @@ -2088,35 +2126,6 @@ def innerIndexed(d, i): m._pyomo_gdp_bigm_reformulation.relaxedDisjuncts, ) - def check_first_disjunct_constraint(self, disj1c, x, ind_var): - self.assertEqual(len(disj1c), 1) - cons = disj1c[0] - self.assertIsNone(cons.lower) - self.assertEqual(cons.upper, 1) - repn = generate_standard_repn(cons.body) - self.assertTrue(repn.is_quadratic()) - self.assertEqual(len(repn.linear_vars), 1) - self.assertEqual(len(repn.quadratic_vars), 4) - ct.check_linear_coef(self, repn, ind_var, 143) - self.assertEqual(repn.constant, -143) - for i in range(1, 5): - ct.check_squared_term_coef(self, repn, x[i], 1) - - def check_second_disjunct_constraint(self, disj2c, x, ind_var): - self.assertEqual(len(disj2c), 1) - cons = disj2c[0] - self.assertIsNone(cons.lower) - self.assertEqual(cons.upper, 1) - repn = generate_standard_repn(cons.body) - self.assertTrue(repn.is_quadratic()) - self.assertEqual(len(repn.linear_vars), 5) - self.assertEqual(len(repn.quadratic_vars), 4) - self.assertEqual(repn.constant, -63) # M = 99, so this is 36 - 99 - ct.check_linear_coef(self, repn, ind_var, 99) - for i in range(1, 5): - ct.check_squared_term_coef(self, repn, x[i], 1) - ct.check_linear_coef(self, repn, x[i], -6) - def simplify_cons(self, cons, leq): visitor = LinearRepnVisitor({}, {}, {}, None) repn = visitor.walk_expression(cons.body) @@ -2142,30 +2151,76 @@ def check_hierarchical_nested_model(self, m, bigm): # outer disjunction constraints disj1c = bigm.get_transformed_constraints(m.disj1.c) - self.check_first_disjunct_constraint(disj1c, m.x, m.disj1.binary_indicator_var) + self.assertEqual(len(disj1c), 1) + cons = disj1c[0] + assertExpressionsEqual( + self, + cons.expr, + m.x[1] ** 2 + + m.x[2] ** 2 + + m.x[3] ** 2 + + m.x[4] ** 2 + - 143.0 * (1 - m.disj1.binary_indicator_var) + <= 1.0, + ) disj2c = bigm.get_transformed_constraints(m.disjunct_block.disj2.c) - self.check_second_disjunct_constraint( - disj2c, m.x, m.disjunct_block.disj2.binary_indicator_var + self.assertEqual(len(disj2c), 1) + cons = disj2c[0] + assertExpressionsEqual( + self, + cons.expr, + (3 - m.x[1]) ** 2 + + (3 - m.x[2]) ** 2 + + (3 - m.x[3]) ** 2 + + (3 - m.x[4]) ** 2 + - 99.0 * (1 - m.disjunct_block.disj2.binary_indicator_var) + <= 1.0, ) # inner disjunction constraints innerd1c = bigm.get_transformed_constraints( m.disjunct_block.disj2.disjunction_disjuncts[0].constraint[1] ) - self.check_first_disjunct_constraint( - innerd1c, - m.x, - m.disjunct_block.disj2.disjunction_disjuncts[0].binary_indicator_var, + self.assertEqual(len(innerd1c), 1) + cons = innerd1c[0] + assertExpressionsEqual( + self, + cons.expr, + m.x[1] ** 2 + + m.x[2] ** 2 + + m.x[3] ** 2 + + m.x[4] ** 2 + - 143.0 + * ( + 1 + - m.disjunct_block.disj2.disjunction_disjuncts[0].binary_indicator_var + + 1 + - m.disjunct_block.disj2.binary_indicator_var + ) + <= 1.0, ) innerd2c = bigm.get_transformed_constraints( m.disjunct_block.disj2.disjunction_disjuncts[1].constraint[1] ) - self.check_second_disjunct_constraint( - innerd2c, - m.x, - m.disjunct_block.disj2.disjunction_disjuncts[1].binary_indicator_var, + self.assertEqual(len(innerd2c), 1) + cons = innerd2c[0] + assertExpressionsEqual( + self, + cons.expr, + (3 - m.x[1]) ** 2 + + (3 - m.x[2]) ** 2 + + (3 - m.x[3]) ** 2 + + (3 - m.x[4]) ** 2 + - 99.0 + * ( + 1 + - m.disjunct_block.disj2.disjunction_disjuncts[1].binary_indicator_var + + 1 + - m.disjunct_block.disj2.binary_indicator_var + ) + <= 1.0, ) def test_hierarchical_badly_ordered_targets(self): @@ -2193,6 +2248,46 @@ def test_decl_order_opposite_instantiation_order(self): def test_do_not_assume_nested_indicators_local(self): ct.check_do_not_assume_nested_indicators_local(self, 'gdp.bigm') + @unittest.skipUnless(gurobi_available, "Gurobi is not available") + def test_constraints_not_enforced_when_an_ancestor_indicator_is_False(self): + m = ConcreteModel() + m.x = Var(bounds=(0, 30)) + + m.left = Disjunct() + m.left.left = Disjunct() + m.left.left.c = Constraint(expr=m.x >= 10) + m.left.right = Disjunct() + m.left.right.c = Constraint(expr=m.x >= 9) + m.left.disjunction = Disjunction(expr=[m.left.left, m.left.right]) + m.right = Disjunct() + m.right.left = Disjunct() + m.right.left.c = Constraint(expr=m.x >= 11) + m.right.right = Disjunct() + m.right.right.c = Constraint(expr=m.x >= 8) + m.right.disjunction = Disjunction(expr=[m.right.left, m.right.right]) + m.disjunction = Disjunction(expr=[m.left, m.right]) + + m.equiv_left = LogicalConstraint( + expr=m.left.left.indicator_var.equivalent_to(m.right.left.indicator_var) + ) + m.equiv_right = LogicalConstraint( + expr=m.left.right.indicator_var.equivalent_to(m.right.right.indicator_var) + ) + + m.obj = Objective(expr=m.x) + + TransformationFactory('gdp.bigm').apply_to(m) + results = SolverFactory('gurobi').solve(m) + self.assertEqual( + results.solver.termination_condition, TerminationCondition.optimal + ) + self.assertTrue(value(m.right.indicator_var)) + self.assertFalse(value(m.left.indicator_var)) + self.assertTrue(value(m.right.right.indicator_var)) + self.assertFalse(value(m.right.left.indicator_var)) + self.assertTrue(value(m.left.right.indicator_var)) + self.assertAlmostEqual(value(m.x), 8) + class IndexedDisjunction(unittest.TestCase): # this tests that if the targets are a subset of the diff --git a/pyomo/gdp/util.py b/pyomo/gdp/util.py index fe11975954d..a8c6393f0b3 100644 --- a/pyomo/gdp/util.py +++ b/pyomo/gdp/util.py @@ -144,13 +144,13 @@ def parent(self, u): Arg: u : A node in the tree """ + if u in self._parent: + return self._parent[u] if u not in self._vertices: raise ValueError( "'%s' is not a vertex in the GDP tree. Cannot " "retrieve its parent." % u ) - if u in self._parent: - return self._parent[u] else: return None diff --git a/pyomo/network/foqus_graph.py b/pyomo/network/foqus_graph.py index e4cf3b92014..7c6c05256d9 100644 --- a/pyomo/network/foqus_graph.py +++ b/pyomo/network/foqus_graph.py @@ -358,9 +358,9 @@ def scc_calculation_order(self, sccNodes, ie, oe): done = False for i in range(len(sccNodes)): for j in range(len(sccNodes)): - for ine in ie[i]: - for oute in oe[j]: - if ine == oute: + for in_e in ie[i]: + for out_e in oe[j]: + if in_e == out_e: adj[j].append(i) adjR[i].append(j) done = True diff --git a/pyomo/repn/linear.py b/pyomo/repn/linear.py index d601ccbcd7c..ba08c7ef245 100644 --- a/pyomo/repn/linear.py +++ b/pyomo/repn/linear.py @@ -200,9 +200,8 @@ def _handle_negation_ANY(visitor, node, arg): _exit_node_handlers[NegationExpression] = { + None: _handle_negation_ANY, (_CONSTANT,): _handle_negation_constant, - (_LINEAR,): _handle_negation_ANY, - (_GENERAL,): _handle_negation_ANY, } # @@ -211,20 +210,18 @@ def _handle_negation_ANY(visitor, node, arg): def _handle_product_constant_constant(visitor, node, arg1, arg2): - _, arg1 = arg1 - _, arg2 = arg2 - ans = arg1 * arg2 + ans = arg1[1] * arg2[1] if ans != ans: - if not arg1 or not arg2: + if not arg1[1] or not arg2[1]: deprecation_warning( - f"Encountered {str(arg1)}*{str(arg2)} in expression tree. " + f"Encountered {str(arg1[1])}*{str(arg2[1])} in expression tree. " "Mapping the NaN result to 0 for compatibility " "with the lp_v1 writer. In the future, this NaN " "will be preserved/emitted to comply with IEEE-754.", version='6.6.0', ) - return _, 0 - return _, arg1 * arg2 + return _CONSTANT, 0 + return _CONSTANT, ans def _handle_product_constant_ANY(visitor, node, arg1, arg2): @@ -276,15 +273,12 @@ def _handle_product_nonlinear(visitor, node, arg1, arg2): _exit_node_handlers[ProductExpression] = { + None: _handle_product_nonlinear, (_CONSTANT, _CONSTANT): _handle_product_constant_constant, (_CONSTANT, _LINEAR): _handle_product_constant_ANY, (_CONSTANT, _GENERAL): _handle_product_constant_ANY, (_LINEAR, _CONSTANT): _handle_product_ANY_constant, - (_LINEAR, _LINEAR): _handle_product_nonlinear, - (_LINEAR, _GENERAL): _handle_product_nonlinear, (_GENERAL, _CONSTANT): _handle_product_ANY_constant, - (_GENERAL, _LINEAR): _handle_product_nonlinear, - (_GENERAL, _GENERAL): _handle_product_nonlinear, } _exit_node_handlers[MonomialTermExpression] = _exit_node_handlers[ProductExpression] @@ -309,15 +303,10 @@ def _handle_division_nonlinear(visitor, node, arg1, arg2): _exit_node_handlers[DivisionExpression] = { + None: _handle_division_nonlinear, (_CONSTANT, _CONSTANT): _handle_division_constant_constant, - (_CONSTANT, _LINEAR): _handle_division_nonlinear, - (_CONSTANT, _GENERAL): _handle_division_nonlinear, (_LINEAR, _CONSTANT): _handle_division_ANY_constant, - (_LINEAR, _LINEAR): _handle_division_nonlinear, - (_LINEAR, _GENERAL): _handle_division_nonlinear, (_GENERAL, _CONSTANT): _handle_division_ANY_constant, - (_GENERAL, _LINEAR): _handle_division_nonlinear, - (_GENERAL, _GENERAL): _handle_division_nonlinear, } # @@ -325,8 +314,7 @@ def _handle_division_nonlinear(visitor, node, arg1, arg2): # -def _handle_pow_constant_constant(visitor, node, *args): - arg1, arg2 = args +def _handle_pow_constant_constant(visitor, node, arg1, arg2): ans = apply_node_operation(node, (arg1[1], arg2[1])) if ans.__class__ in native_complex_types: ans = complex_number_error(ans, visitor, node) @@ -358,15 +346,10 @@ def _handle_pow_nonlinear(visitor, node, arg1, arg2): _exit_node_handlers[PowExpression] = { + None: _handle_pow_nonlinear, (_CONSTANT, _CONSTANT): _handle_pow_constant_constant, - (_CONSTANT, _LINEAR): _handle_pow_nonlinear, - (_CONSTANT, _GENERAL): _handle_pow_nonlinear, (_LINEAR, _CONSTANT): _handle_pow_ANY_constant, - (_LINEAR, _LINEAR): _handle_pow_nonlinear, - (_LINEAR, _GENERAL): _handle_pow_nonlinear, (_GENERAL, _CONSTANT): _handle_pow_ANY_constant, - (_GENERAL, _LINEAR): _handle_pow_nonlinear, - (_GENERAL, _GENERAL): _handle_pow_nonlinear, } # @@ -389,9 +372,8 @@ def _handle_unary_nonlinear(visitor, node, arg): _exit_node_handlers[UnaryFunctionExpression] = { + None: _handle_unary_nonlinear, (_CONSTANT,): _handle_unary_constant, - (_LINEAR,): _handle_unary_nonlinear, - (_GENERAL,): _handle_unary_nonlinear, } _exit_node_handlers[AbsExpression] = _exit_node_handlers[UnaryFunctionExpression] @@ -414,9 +396,8 @@ def _handle_named_ANY(visitor, node, arg1): _exit_node_handlers[Expression] = { + None: _handle_named_ANY, (_CONSTANT,): _handle_named_constant, - (_LINEAR,): _handle_named_ANY, - (_GENERAL,): _handle_named_ANY, } # @@ -449,12 +430,7 @@ def _handle_expr_if_nonlinear(visitor, node, arg1, arg2, arg3): return _GENERAL, ans -_exit_node_handlers[Expr_ifExpression] = { - (i, j, k): _handle_expr_if_nonlinear - for i in (_LINEAR, _GENERAL) - for j in (_CONSTANT, _LINEAR, _GENERAL) - for k in (_CONSTANT, _LINEAR, _GENERAL) -} +_exit_node_handlers[Expr_ifExpression] = {None: _handle_expr_if_nonlinear} for j in (_CONSTANT, _LINEAR, _GENERAL): for k in (_CONSTANT, _LINEAR, _GENERAL): _exit_node_handlers[Expr_ifExpression][_CONSTANT, j, k] = _handle_expr_if_const @@ -487,11 +463,9 @@ def _handle_equality_general(visitor, node, arg1, arg2): _exit_node_handlers[EqualityExpression] = { - (i, j): _handle_equality_general - for i in (_CONSTANT, _LINEAR, _GENERAL) - for j in (_CONSTANT, _LINEAR, _GENERAL) + None: _handle_equality_general, + (_CONSTANT, _CONSTANT): _handle_equality_const, } -_exit_node_handlers[EqualityExpression][_CONSTANT, _CONSTANT] = _handle_equality_const def _handle_inequality_const(visitor, node, arg1, arg2): @@ -517,13 +491,9 @@ def _handle_inequality_general(visitor, node, arg1, arg2): _exit_node_handlers[InequalityExpression] = { - (i, j): _handle_inequality_general - for i in (_CONSTANT, _LINEAR, _GENERAL) - for j in (_CONSTANT, _LINEAR, _GENERAL) + None: _handle_inequality_general, + (_CONSTANT, _CONSTANT): _handle_inequality_const, } -_exit_node_handlers[InequalityExpression][ - _CONSTANT, _CONSTANT -] = _handle_inequality_const def _handle_ranged_const(visitor, node, arg1, arg2, arg3): @@ -554,14 +524,9 @@ def _handle_ranged_general(visitor, node, arg1, arg2, arg3): _exit_node_handlers[RangedExpression] = { - (i, j, k): _handle_ranged_general - for i in (_CONSTANT, _LINEAR, _GENERAL) - for j in (_CONSTANT, _LINEAR, _GENERAL) - for k in (_CONSTANT, _LINEAR, _GENERAL) + None: _handle_ranged_general, + (_CONSTANT, _CONSTANT, _CONSTANT): _handle_ranged_const, } -_exit_node_handlers[RangedExpression][ - _CONSTANT, _CONSTANT, _CONSTANT -] = _handle_ranged_const class LinearBeforeChildDispatcher(BeforeChildDispatcher): @@ -754,7 +719,10 @@ def _initialize_exit_node_dispatcher(exit_handlers): exit_dispatcher = {} for cls, handlers in exit_handlers.items(): for args, fcn in handlers.items(): - exit_dispatcher[(cls, *args)] = fcn + if args is None: + exit_dispatcher[cls] = fcn + else: + exit_dispatcher[(cls, *args)] = fcn return exit_dispatcher diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index 239cd845930..ea7b6a6a9e6 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -139,6 +139,15 @@ class LinearStandardFormCompiler(object): description='Add slack variables and return `min cTx s.t. Ax == b`', ), ) + CONFIG.declare( + 'mixed_form', + ConfigValue( + default=False, + domain=bool, + description='Return A in mixed form (the comparison operator is a ' + 'mix of <=, ==, and >=)', + ), + ) CONFIG.declare( 'show_section_timing', ConfigValue( @@ -332,6 +341,9 @@ def write(self, model): # Tabulate constraints # slack_form = self.config.slack_form + mixed_form = self.config.mixed_form + if slack_form and mixed_form: + raise ValueError("cannot specify both slack_form and mixed_form") rows = [] rhs = [] con_data = [] @@ -372,7 +384,30 @@ def write(self, model): f"model contains a trivially infeasible constraint, '{con.name}'" ) - if slack_form: + if mixed_form: + N = len(repn.linear) + _data = np.fromiter(repn.linear.values(), float, N) + _index = np.fromiter(map(var_order.__getitem__, repn.linear), float, N) + if ub == lb: + rows.append(RowEntry(con, 0)) + rhs.append(ub - offset) + con_data.append(_data) + con_index.append(_index) + con_index_ptr.append(con_index_ptr[-1] + N) + else: + if ub is not None: + rows.append(RowEntry(con, 1)) + rhs.append(ub - offset) + con_data.append(_data) + con_index.append(_index) + con_index_ptr.append(con_index_ptr[-1] + N) + if lb is not None: + rows.append(RowEntry(con, -1)) + rhs.append(lb - offset) + con_data.append(_data) + con_index.append(_index) + con_index_ptr.append(con_index_ptr[-1] + N) + elif slack_form: _data = list(repn.linear.values()) _index = list(map(var_order.__getitem__, repn.linear)) if lb == ub: # TODO: add tolerance? @@ -437,24 +472,22 @@ def write(self, model): # at the index pointer list (an O(num_var) operation). c_ip = c.indptr A_ip = A.indptr - active_var_idx = list( - filter( - lambda i: A_ip[i] != A_ip[i + 1] or c_ip[i] != c_ip[i + 1], - range(len(columns)), - ) - ) - nCol = len(active_var_idx) + active_var_mask = (A_ip[1:] > A_ip[:-1]) | (c_ip[1:] > c_ip[:-1]) + + # Masks on NumPy arrays are very fast. Build the reduced A + # indptr and then check if we actually have to manipulate the + # columns + augmented_mask = np.concatenate((active_var_mask, [True])) + reduced_A_indptr = A.indptr[augmented_mask] + nCol = len(reduced_A_indptr) - 1 if nCol != len(columns): - # Note that the indptr can't just use range() because a var - # may only appear in the objectives or the constraints. - columns = list(map(columns.__getitem__, active_var_idx)) - active_var_idx.append(c.indptr[-1]) + columns = [v for k, v in zip(active_var_mask, columns) if k] c = scipy.sparse.csc_array( - (c.data, c.indices, c.indptr.take(active_var_idx)), [c.shape[0], nCol] + (c.data, c.indices, c.indptr[augmented_mask]), [c.shape[0], nCol] ) - active_var_idx[-1] = A.indptr[-1] + # active_var_idx[-1] = len(columns) A = scipy.sparse.csc_array( - (A.data, A.indices, A.indptr.take(active_var_idx)), [A.shape[0], nCol] + (A.data, A.indices, reduced_A_indptr), [A.shape[0], nCol] ) if self.config.nonnegative_vars: diff --git a/pyomo/repn/quadratic.py b/pyomo/repn/quadratic.py index 0ddfda829ed..f6e0a43623d 100644 --- a/pyomo/repn/quadratic.py +++ b/pyomo/repn/quadratic.py @@ -277,18 +277,11 @@ def _handle_product_nonlinear(visitor, node, arg1, arg2): _exit_node_handlers[ProductExpression].update( { + None: _handle_product_nonlinear, (_CONSTANT, _QUADRATIC): linear._handle_product_constant_ANY, - (_LINEAR, _QUADRATIC): _handle_product_nonlinear, - (_QUADRATIC, _QUADRATIC): _handle_product_nonlinear, - (_GENERAL, _QUADRATIC): _handle_product_nonlinear, (_QUADRATIC, _CONSTANT): linear._handle_product_ANY_constant, - (_QUADRATIC, _LINEAR): _handle_product_nonlinear, - (_QUADRATIC, _GENERAL): _handle_product_nonlinear, # Replace handler from the linear walker (_LINEAR, _LINEAR): _handle_product_linear_linear, - (_GENERAL, _GENERAL): _handle_product_nonlinear, - (_GENERAL, _LINEAR): _handle_product_nonlinear, - (_LINEAR, _GENERAL): _handle_product_nonlinear, } ) @@ -296,15 +289,7 @@ def _handle_product_nonlinear(visitor, node, arg1, arg2): # DIVISION # _exit_node_handlers[DivisionExpression].update( - { - (_CONSTANT, _QUADRATIC): linear._handle_division_nonlinear, - (_LINEAR, _QUADRATIC): linear._handle_division_nonlinear, - (_QUADRATIC, _QUADRATIC): linear._handle_division_nonlinear, - (_GENERAL, _QUADRATIC): linear._handle_division_nonlinear, - (_QUADRATIC, _CONSTANT): linear._handle_division_ANY_constant, - (_QUADRATIC, _LINEAR): linear._handle_division_nonlinear, - (_QUADRATIC, _GENERAL): linear._handle_division_nonlinear, - } + {(_QUADRATIC, _CONSTANT): linear._handle_division_ANY_constant} ) @@ -312,84 +297,42 @@ def _handle_product_nonlinear(visitor, node, arg1, arg2): # EXPONENTIATION # _exit_node_handlers[PowExpression].update( - { - (_CONSTANT, _QUADRATIC): linear._handle_pow_nonlinear, - (_LINEAR, _QUADRATIC): linear._handle_pow_nonlinear, - (_QUADRATIC, _QUADRATIC): linear._handle_pow_nonlinear, - (_GENERAL, _QUADRATIC): linear._handle_pow_nonlinear, - (_QUADRATIC, _CONSTANT): linear._handle_pow_ANY_constant, - (_QUADRATIC, _LINEAR): linear._handle_pow_nonlinear, - (_QUADRATIC, _GENERAL): linear._handle_pow_nonlinear, - } + {(_QUADRATIC, _CONSTANT): linear._handle_pow_ANY_constant} ) # # ABS and UNARY handlers # -_exit_node_handlers[AbsExpression][(_QUADRATIC,)] = linear._handle_unary_nonlinear -_exit_node_handlers[UnaryFunctionExpression][ - (_QUADRATIC,) -] = linear._handle_unary_nonlinear +# (no changes needed) # # NAMED EXPRESSION handlers # -_exit_node_handlers[Expression][(_QUADRATIC,)] = linear._handle_named_ANY +# (no changes needed) # # EXPR_IF handlers # # Note: it is easier to just recreate the entire data structure, rather # than update it -_exit_node_handlers[Expr_ifExpression] = { - (i, j, k): linear._handle_expr_if_nonlinear - for i in (_LINEAR, _QUADRATIC, _GENERAL) - for j in (_CONSTANT, _LINEAR, _QUADRATIC, _GENERAL) - for k in (_CONSTANT, _LINEAR, _QUADRATIC, _GENERAL) -} -for j in (_CONSTANT, _LINEAR, _QUADRATIC, _GENERAL): - for k in (_CONSTANT, _LINEAR, _QUADRATIC, _GENERAL): - _exit_node_handlers[Expr_ifExpression][ - _CONSTANT, j, k - ] = linear._handle_expr_if_const - -# -# RELATIONAL handlers -# -_exit_node_handlers[EqualityExpression].update( +_exit_node_handlers[Expr_ifExpression].update( { - (_CONSTANT, _QUADRATIC): linear._handle_equality_general, - (_LINEAR, _QUADRATIC): linear._handle_equality_general, - (_QUADRATIC, _QUADRATIC): linear._handle_equality_general, - (_GENERAL, _QUADRATIC): linear._handle_equality_general, - (_QUADRATIC, _CONSTANT): linear._handle_equality_general, - (_QUADRATIC, _LINEAR): linear._handle_equality_general, - (_QUADRATIC, _GENERAL): linear._handle_equality_general, + (_CONSTANT, i, _QUADRATIC): linear._handle_expr_if_const + for i in (_CONSTANT, _LINEAR, _QUADRATIC, _GENERAL) } ) -_exit_node_handlers[InequalityExpression].update( +_exit_node_handlers[Expr_ifExpression].update( { - (_CONSTANT, _QUADRATIC): linear._handle_inequality_general, - (_LINEAR, _QUADRATIC): linear._handle_inequality_general, - (_QUADRATIC, _QUADRATIC): linear._handle_inequality_general, - (_GENERAL, _QUADRATIC): linear._handle_inequality_general, - (_QUADRATIC, _CONSTANT): linear._handle_inequality_general, - (_QUADRATIC, _LINEAR): linear._handle_inequality_general, - (_QUADRATIC, _GENERAL): linear._handle_inequality_general, - } -) -_exit_node_handlers[RangedExpression].update( - { - (_CONSTANT, _QUADRATIC): linear._handle_ranged_general, - (_LINEAR, _QUADRATIC): linear._handle_ranged_general, - (_QUADRATIC, _QUADRATIC): linear._handle_ranged_general, - (_GENERAL, _QUADRATIC): linear._handle_ranged_general, - (_QUADRATIC, _CONSTANT): linear._handle_ranged_general, - (_QUADRATIC, _LINEAR): linear._handle_ranged_general, - (_QUADRATIC, _GENERAL): linear._handle_ranged_general, + (_CONSTANT, _QUADRATIC, i): linear._handle_expr_if_const + for i in (_CONSTANT, _LINEAR, _GENERAL) } ) +# +# RELATIONAL handlers +# +# (no changes needed) + class QuadraticRepnVisitor(linear.LinearRepnVisitor): Result = QuadraticRepn diff --git a/pyomo/repn/tests/test_standard_form.py b/pyomo/repn/tests/test_standard_form.py index e24195edfde..4c66ae87c41 100644 --- a/pyomo/repn/tests/test_standard_form.py +++ b/pyomo/repn/tests/test_standard_form.py @@ -42,6 +42,23 @@ def test_linear_model(self): self.assertTrue(np.all(repn.c == np.array([0, 0, 0]))) self.assertTrue(np.all(repn.A == np.array([[-1, -2, 0], [0, 1, 4]]))) self.assertTrue(np.all(repn.rhs == np.array([-3, 5]))) + self.assertEqual(repn.rows, [(m.c, -1), (m.d, 1)]) + self.assertEqual(repn.columns, [m.x, m.y[1], m.y[3]]) + + def test_almost_dense_linear_model(self): + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.y = pyo.Var([1, 2, 3]) + m.c = pyo.Constraint(expr=m.x + 2 * m.y[1] + 4 * m.y[3] >= 10) + m.d = pyo.Constraint(expr=5 * m.x + 6 * m.y[1] + 8 * m.y[3] <= 20) + + repn = LinearStandardFormCompiler().write(m) + + self.assertTrue(np.all(repn.c == np.array([0, 0, 0]))) + self.assertTrue(np.all(repn.A == np.array([[-1, -2, -4], [5, 6, 8]]))) + self.assertTrue(np.all(repn.rhs == np.array([-10, 20]))) + self.assertEqual(repn.rows, [(m.c, -1), (m.d, 1)]) + self.assertEqual(repn.columns, [m.x, m.y[1], m.y[3]]) def test_linear_model_row_col_order(self): m = pyo.ConcreteModel() @@ -57,6 +74,8 @@ def test_linear_model_row_col_order(self): self.assertTrue(np.all(repn.c == np.array([0, 0, 0]))) self.assertTrue(np.all(repn.A == np.array([[4, 0, 1], [0, -1, -2]]))) self.assertTrue(np.all(repn.rhs == np.array([5, -3]))) + self.assertEqual(repn.rows, [(m.d, 1), (m.c, -1)]) + self.assertEqual(repn.columns, [m.y[3], m.x, m.y[1]]) def test_suffix_warning(self): m = pyo.ConcreteModel() @@ -222,6 +241,28 @@ def test_alternative_forms(self): ) self._verify_solution(soln, repn, True) + repn = LinearStandardFormCompiler().write( + m, mixed_form=True, column_order=col_order + ) + + self.assertEqual( + repn.rows, [(m.c, -1), (m.d, 1), (m.e, 1), (m.e, -1), (m.f, 0)] + ) + self.assertEqual(list(map(str, repn.x)), ['x', 'y[0]', 'y[1]', 'y[3]']) + self.assertEqual( + list(v.bounds for v in repn.x), [(None, None), (0, 10), (-5, 10), (-5, -2)] + ) + ref = np.array( + [[1, 0, 2, 0], [0, 0, 1, 4], [0, 1, 6, 0], [0, 1, 6, 0], [1, 1, 0, 0]] + ) + self.assertTrue(np.all(repn.A == ref)) + self.assertTrue(np.all(repn.b == np.array([3, 5, 6, -3, 8]))) + self.assertTrue(np.all(repn.c == np.array([[-1, 0, -5, 0], [1, 0, 0, 15]]))) + # Note that the mixed_form solution is a mix of inequality and + # equality constraints, so we cannot (easily) reuse the + # _verify_solutions helper (as in the above cases): + # self._verify_solution(soln, repn, False) + repn = LinearStandardFormCompiler().write( m, slack_form=True, nonnegative_vars=True, column_order=col_order ) diff --git a/pyomo/repn/tests/test_util.py b/pyomo/repn/tests/test_util.py index b5e4cc4facf..e0fea0fb45c 100644 --- a/pyomo/repn/tests/test_util.py +++ b/pyomo/repn/tests/test_util.py @@ -718,16 +718,14 @@ class UnknownExpression(NumericExpression): DeveloperError, r".*Unexpected expression node type 'UnknownExpression'" ): end[node.__class__](None, node, *node.args) - self.assertEqual(len(end), 9) - self.assertIn(UnknownExpression, end) + self.assertEqual(len(end), 8) node = UnknownExpression((6, 7)) with self.assertRaisesRegex( DeveloperError, r".*Unexpected expression node type 'UnknownExpression'" ): end[node.__class__, 6, 7](None, node, *node.args) - self.assertEqual(len(end), 10) - self.assertIn((UnknownExpression, 6, 7), end) + self.assertEqual(len(end), 8) def test_BeforeChildDispatcher_registration(self): class BeforeChildDispatcherTester(BeforeChildDispatcher): diff --git a/pyomo/repn/util.py b/pyomo/repn/util.py index 49cca32eaf9..1b3056738bf 100644 --- a/pyomo/repn/util.py +++ b/pyomo/repn/util.py @@ -400,7 +400,15 @@ def __init__(self, *args, **kwargs): def __missing__(self, key): if type(key) is tuple: - node_class = key[0] + # Only lookup/cache argument-specific handlers for unary, + # binary and ternary operators + if len(key) <= 3: + node_class = key[0] + node_args = key[1:] + else: + node_class = key = key[0] + if node_class in self: + return self[node_class] else: node_class = key bases = node_class.__mro__ @@ -412,30 +420,31 @@ def __missing__(self, key): bases = [Expression] fcn = None for base_type in bases: - if isinstance(key, tuple): - base_key = (base_type,) + key[1:] - # Only cache handlers for unary, binary and ternary operators - cache = len(key) <= 4 - else: - base_key = base_type - cache = True - if base_key in self: - fcn = self[base_key] - elif base_type in self: + if key is not node_class: + if (base_type,) + node_args in self: + fcn = self[(base_type,) + node_args] + break + if base_type in self: fcn = self[base_type] - elif any((k[0] if type(k) is tuple else k) is base_type for k in self): - raise DeveloperError( - f"Base expression key '{base_key}' not found when inserting " - f"dispatcher for node '{node_class.__name__}' while walking " - "expression tree." - ) + break if fcn is None: - fcn = self.unexpected_expression_type - if cache: - self[key] = fcn + partial_matches = set( + k[0] for k in self if type(k) is tuple and issubclass(node_class, k[0]) + ) + for base_type in node_class.__mro__: + if node_class is not key: + key = (base_type,) + node_args + if base_type in partial_matches: + raise DeveloperError( + f"Base expression key '{key}' not found when inserting " + f"dispatcher for node '{node_class.__name__}' while walking " + "expression tree." + ) + return self.unexpected_expression_type + self[key] = fcn return fcn - def unexpected_expression_type(self, visitor, node, *arg): + def unexpected_expression_type(self, visitor, node, *args): raise DeveloperError( f"Unexpected expression node type '{type(node).__name__}' " f"found while walking expression tree in {type(visitor).__name__}." diff --git a/pyomo/solvers/plugins/solvers/GAMS.py b/pyomo/solvers/plugins/solvers/GAMS.py index be3499a2f6b..c0bab4dc23e 100644 --- a/pyomo/solvers/plugins/solvers/GAMS.py +++ b/pyomo/solvers/plugins/solvers/GAMS.py @@ -198,8 +198,8 @@ def _get_version(self): return _extract_version('') from gams import GamsWorkspace - ws = GamsWorkspace() - version = tuple(int(i) for i in ws._version.split('.')[:4]) + workspace = GamsWorkspace() + version = tuple(int(i) for i in workspace._version.split('.')[:4]) while len(version) < 4: version += (0,) return version @@ -209,8 +209,8 @@ def _run_simple_model(self, n): try: from gams import GamsWorkspace, DebugLevel - ws = GamsWorkspace(debug=DebugLevel.Off, working_directory=tmpdir) - t1 = ws.add_job_from_string(self._simple_model(n)) + workspace = GamsWorkspace(debug=DebugLevel.Off, working_directory=tmpdir) + t1 = workspace.add_job_from_string(self._simple_model(n)) t1.run() return True except: @@ -330,12 +330,12 @@ def solve(self, *args, **kwds): if tmpdir is not None and os.path.exists(tmpdir): newdir = False - ws = GamsWorkspace( + workspace = GamsWorkspace( debug=DebugLevel.KeepFiles if keepfiles else DebugLevel.Off, working_directory=tmpdir, ) - t1 = ws.add_job_from_string(output_file.getvalue()) + t1 = workspace.add_job_from_string(output_file.getvalue()) try: with OutputStream(tee=tee, logfile=logfile) as output_stream: @@ -349,7 +349,9 @@ def solve(self, *args, **kwds): # Always name working directory or delete files, # regardless of any errors. if keepfiles: - print("\nGAMS WORKING DIRECTORY: %s\n" % ws.working_directory) + print( + "\nGAMS WORKING DIRECTORY: %s\n" % workspace.working_directory + ) elif tmpdir is not None: # Garbage collect all references to t1.out_db # So that .gdx file can be deleted @@ -359,7 +361,7 @@ def solve(self, *args, **kwds): except: # Catch other errors and remove files first if keepfiles: - print("\nGAMS WORKING DIRECTORY: %s\n" % ws.working_directory) + print("\nGAMS WORKING DIRECTORY: %s\n" % workspace.working_directory) elif tmpdir is not None: # Garbage collect all references to t1.out_db # So that .gdx file can be deleted @@ -398,7 +400,9 @@ def solve(self, *args, **kwds): extract_rc = 'rc' in model_suffixes results = SolverResults() - results.problem.name = os.path.join(ws.working_directory, t1.name + '.gms') + results.problem.name = os.path.join( + workspace.working_directory, t1.name + '.gms' + ) results.problem.lower_bound = t1.out_db["OBJEST"].find_record().value results.problem.upper_bound = t1.out_db["OBJEST"].find_record().value results.problem.number_of_variables = t1.out_db["NUMVAR"].find_record().value @@ -587,7 +591,7 @@ def solve(self, *args, **kwds): results.solution.insert(soln) if keepfiles: - print("\nGAMS WORKING DIRECTORY: %s\n" % ws.working_directory) + print("\nGAMS WORKING DIRECTORY: %s\n" % workspace.working_directory) elif tmpdir is not None: # Garbage collect all references to t1.out_db # So that .gdx file can be deleted