From bae88522907a423fcca0c47f8bc6f22b2d63cbf4 Mon Sep 17 00:00:00 2001 From: John Siirola Date: Sun, 26 Nov 2023 15:45:10 -0700 Subject: [PATCH 1/9] NFC: clean up docstring --- pyomo/repn/plugins/nl_writer.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pyomo/repn/plugins/nl_writer.py b/pyomo/repn/plugins/nl_writer.py index 8dfaf0d7f42..6282dc95ce5 100644 --- a/pyomo/repn/plugins/nl_writer.py +++ b/pyomo/repn/plugins/nl_writer.py @@ -165,9 +165,9 @@ class NLWriterInfo(object): eliminated_vars: List[Tuple[_VarData, NumericExpression]] The list of variables in the model that were eliminated by the - presolve. each entry is a 2-tuple of (:py:class:`_VarData`, - :py:class`NumericExpression`|`float`). the list is ordered in - the necessary order for correct evaluation (i.e., all variables + presolve. Each entry is a 2-tuple of (:py:class:`_VarData`, + :py:class`NumericExpression`|`float`). The list is in the + necessary order for correct evaluation (i.e., all variables appearing in the expression must either have been sent to the solver, or appear *earlier* in this list. From 8262bd5430654ea4f8bef3454f1a0106fab70371 Mon Sep 17 00:00:00 2001 From: John Siirola Date: Sun, 26 Nov 2023 15:47:02 -0700 Subject: [PATCH 2/9] Initial implementation and tests for generating standard linear forms --- pyomo/repn/plugins/__init__.py | 1 + pyomo/repn/plugins/standard_form.py | 540 +++++++++++++++++++++++++ pyomo/repn/tests/test_standard_form.py | 267 ++++++++++++ 3 files changed, 808 insertions(+) create mode 100644 pyomo/repn/plugins/standard_form.py create mode 100644 pyomo/repn/tests/test_standard_form.py diff --git a/pyomo/repn/plugins/__init__.py b/pyomo/repn/plugins/__init__.py index f1e8270b8c7..56b221d3129 100644 --- a/pyomo/repn/plugins/__init__.py +++ b/pyomo/repn/plugins/__init__.py @@ -18,6 +18,7 @@ def load(): import pyomo.repn.plugins.gams_writer import pyomo.repn.plugins.lp_writer import pyomo.repn.plugins.nl_writer + import pyomo.repn.plugins.standard_form from pyomo.opt import WriterFactory diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py new file mode 100644 index 00000000000..6e74faca7d1 --- /dev/null +++ b/pyomo/repn/plugins/standard_form.py @@ -0,0 +1,540 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2022 +# 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 collections +import logging +from operator import attrgetter, neg + +from pyomo.common.config import ( + ConfigBlock, + ConfigValue, + InEnum, + document_kwargs_from_configdict, +) +from pyomo.common.dependencies import scipy, numpy as np +from pyomo.common.gc_manager import PauseGC +from pyomo.common.timing import TicTocTimer + +from pyomo.core.base import ( + Block, + Objective, + Constraint, + Var, + Param, + Expression, + SOSConstraint, + SortComponents, + Suffix, + SymbolMap, + minimize, +) +from pyomo.core.base.component import ActiveComponent +from pyomo.core.base.label import LPFileLabeler, NumericLabeler +from pyomo.opt import WriterFactory +from pyomo.repn.linear import LinearRepnVisitor +from pyomo.repn.quadratic import QuadraticRepnVisitor +from pyomo.repn.util import ( + FileDeterminism, + FileDeterminism_to_SortComponents, + categorize_valid_components, + initialize_var_map_from_column_order, + int_float, + ordered_active_constraints, +) + +### FIXME: Remove the following as soon as non-active components no +### longer report active==True +from pyomo.core.base import Set, RangeSet, ExternalFunction +from pyomo.network import Port + +logger = logging.getLogger(__name__) +inf = float('inf') +neg_inf = float('-inf') + + +RowEntry = collections.namedtuple('RowEntry', ['constraint', 'bound_type']) + + +# TODO: make a proper base class +class LinearStandardFormInfo(object): + """Return type for LinearStandardFormCompiler.write() + + Attributes + ---------- + c : scipy.sparse.csr_array + + The objective coefficients. Note that this is a sparse array + and may contain multiple rows (for multiobjective problems). The + objectives may be calculated by "c @ x" + + A : scipy.sparse.csc_array + + The constraint coefficients. The constraint bodies may be + calculated by "A @ x" + + rhs : numpy.ndarray + + The constraint right-hand sides. + + rows : List[Tuple[_ConstraintData, int]] + + The list of Pyomo constraint objects corresponding to the rows + in `A`. Each element in the list is a 2-tuple of + (_ConstraintData, row_multiplier). The `row_multiplier` will be + +/- 1 (indicating if the row was multiplied by -1 (corresponding + to a constraint lower bound or +1 (upper bound). + + columns : List[_VarData] + + The list of Pyomo variable objects corresponding to columns in + the `A` and `c` matricies. + + eliminated_vars: List[Tuple[_VarData, NumericExpression]] + + The list of variables from the original model that do not appear + in the standard form (usually because they were replaced by + nonnegative variables). Each entry is a 2-tuple of + (:py:class:`_VarData`, :py:class`NumericExpression`|`float`). + The list is in the necessary order for correct evaluation (i.e., + all variables appearing in the expression must either have + appeared in the standard form, or appear *earlier* in this list. + + """ + + def __init__(self, c, A, rhs, rows, columns, eliminated_vars): + self.c = c + self.A = A + self.rhs = rhs + self.rows = rows + self.columns = columns + self.eliminated_vars = eliminated_vars + + @property + def x(self): + return self.columns + + @property + def b(self): + return self.rhs + + +@WriterFactory.register( + 'compile_standard_form', 'Compile an LP to standard form (`min cTx s.t. Ax <= b`)' +) +class LinearStandardFormCompiler(object): + CONFIG = ConfigBlock('compile_standard_form') + CONFIG.declare( + 'nonnegative_vars', + ConfigValue( + default=False, + domain=bool, + description='Convert all variables to be nonnegative variables', + ), + ) + CONFIG.declare( + 'slack_form', + ConfigValue( + default=False, + domain=bool, + description='Add slack variables and return `min cTx s.t. Ax == b`', + ), + ) + CONFIG.declare( + 'show_section_timing', + ConfigValue( + default=False, + domain=bool, + description='Print timing after writing each section of the LP file', + ), + ) + CONFIG.declare( + 'file_determinism', + ConfigValue( + default=FileDeterminism.ORDERED, + domain=InEnum(FileDeterminism), + description='How much effort to ensure result is deterministic', + doc=""" + How much effort do we want to put into ensuring the + resulting matricies are produced deterministically: + NONE (0) : None + ORDERED (10): rely on underlying component ordering (default) + SORT_INDICES (20) : sort keys of indexed components + SORT_SYMBOLS (30) : sort keys AND sort names (not declaration order) + """, + ), + ) + CONFIG.declare( + 'row_order', + ConfigValue( + default=None, + description='Preferred constraint ordering', + doc=""" + List of constraints in the order that they should appear in the + LP file. Unspecified constraints will appear at the end.""", + ), + ) + CONFIG.declare( + 'column_order', + ConfigValue( + default=None, + description='Preferred variable ordering', + doc=""" + List of variables in the order that they should appear in + the LP file. Note that this is only a suggestion, as the LP + file format is row-major and the columns are inferred from + the order in which variables appear in the objective + followed by each constraint.""", + ), + ) + + def __init__(self): + self.config = self.CONFIG() + + @document_kwargs_from_configdict(CONFIG) + def write(self, model, ostream=None, **options): + """Convert a model to standard form (`min cTx s.t. Ax <= b`) + + Returns + ------- + LinearStandardFormInfo + + Parameters + ---------- + model: ConcreteModel + The concrete Pyomo model to write out. + + ostream: None + This is provided for API compatibility with other writers + and is ignored here. + + """ + config = self.config(options) + + # Pause the GC, as the walker that generates the compiled LP + # representation generates (and disposes of) a large number of + # small objects. + with PauseGC(): + return _LinearStandardFormCompiler_impl(config).write(model) + + +class _LinearStandardFormCompiler_impl(object): + def __init__(self, config): + self.config = config + + def write(self, model): + timing_logger = logging.getLogger('pyomo.common.timing.writer') + timer = TicTocTimer(logger=timing_logger) + with_debug_timing = ( + timing_logger.isEnabledFor(logging.DEBUG) and timing_logger.hasHandlers() + ) + + sorter = FileDeterminism_to_SortComponents(self.config.file_determinism) + component_map, unknown = categorize_valid_components( + model, + active=True, + sort=sorter, + valid={ + Block, + Constraint, + Var, + Param, + Expression, + # FIXME: Non-active components should not report as Active + ExternalFunction, + Set, + RangeSet, + Port, + # TODO: Piecewise, Complementarity + }, + targets={Suffix, Objective}, + ) + if unknown: + raise ValueError( + "The model ('%s') contains the following active components " + "that the LP compiler does not know how to process:\n\t%s" + % ( + model.name, + "\n\t".join( + "%s:\n\t\t%s" % (k, "\n\t\t".join(map(attrgetter('name'), v))) + for k, v in unknown.items() + ), + ) + ) + + self.var_map = var_map = {} + initialize_var_map_from_column_order(model, self.config, var_map) + var_order = {_id: i for i, _id in enumerate(var_map)} + + visitor = LinearRepnVisitor({}, var_map, var_order) + + timer.toc('Initialized column order', level=logging.DEBUG) + + # We don't export any suffix information to the Standard Form + # + if component_map[Suffix]: + suffixesByName = {} + for block in component_map[Suffix]: + for suffix in block.component_objects( + Suffix, active=True, descend_into=False, sort=sorter + ): + if not suffix.export_enabled() or not suffix: + continue + name = suffix.local_name + if name in suffixesByName: + suffixesByName[name].append(suffix) + else: + suffixesByName[name] = [suffix] + for name, suffixes in suffixesByName.items(): + n = len(suffixes) + plural = 's' if n > 1 else '' + logger.warning( + f"EXPORT Suffix '{name}' found on {n} block{plural}:\n " + + "\n ".join(s.name for s in suffixes) + + "\nStandard Form compiler ignores export suffixes. Skipping." + ) + + # + # Process objective + # + if not component_map[Objective]: + objectives = [Objective(expr=1)] + objectives[0].construct() + else: + objectives = [] + for blk in component_map[Objective]: + objectives.extend( + blk.component_data_objects( + Objective, active=True, descend_into=False, sort=sorter + ) + ) + obj_data = [] + obj_index = [] + obj_index_ptr = [0] + for i, obj in enumerate(objectives): + repn = visitor.walk_expression(obj.expr) + if repn.nonlinear is not None: + raise ValueError( + f"Model objective ({obj.name}) contains nonlinear terms that " + "cannot be compiled to standard (linear) form." + ) + obj_data.extend(repn.linear.values()) + obj_index.extend(map(var_order.__getitem__, repn.linear)) + obj_index_ptr.append(len(obj_index)) + if with_debug_timing: + timer.toc('Objective %s', obj, level=logging.DEBUG) + + # + # Tabulate constraints + # + slack_form = self.config.slack_form + rows = [] + rhs = [] + con_data = [] + con_index = [] + con_index_ptr = [0] + last_parent = None + for con in ordered_active_constraints(model, self.config): + if with_debug_timing and con.parent_component() is not last_parent: + if last_parent is not None: + timer.toc('Constraint %s', last_parent, level=logging.DEBUG) + last_parent = con.parent_component() + # Note: Constraint.lb/ub guarantee a return value that is + # either a (finite) native_numeric_type, or None + lb = con.lb + ub = con.ub + + repn = visitor.walk_expression(con.body) + + if lb is None and ub is None: + # Note: you *cannot* output trivial (unbounded) + # constraints in matrix format. I suppose we could add a + # slack variable, but that seems rather silly. + continue + if repn.nonlinear is not None: + raise ValueError( + f"Model constraint ({con.name}) contains nonlinear terms that " + "cannot be compiled to standard (linear) form." + ) + + # Pull out the constant: we will move it to the bounds + offset = repn.constant + repn.constant = 0 + + if not repn.linear: + if (lb is None or lb <= offset) and (ub is None or ub >= offset): + continue + raise InfeasibleError( + f"model contains a trivially infeasible constraint, '{con.name}'" + ) + + if slack_form: + _data = list(repn.linear.values()) + _index = list(map(var_order.__getitem__, repn.linear)) + if lb == ub: # TODO: add tolerance? + rhs.append(ub - offset) + else: + # add slack variable + v = Var(name=f'_slack_{len(rhs)}', bounds=(None, None)) + v.construct() + if lb is None: + rhs.append(ub - offset) + v.lb = 0 + else: + rhs.append(lb - offset) + v.ub = 0 + if ub is not None: + v.lb = lb - ub + var_map[id(v)] = v + var_order[id(v)] = slack_col = len(var_order) + _data.append(1) + _index.append(slack_col) + rows.append(RowEntry(con, 1)) + con_data.append(np.array(_data)) + con_index.append(np.array(_index)) + con_index_ptr.append(con_index_ptr[-1] + len(_index)) + else: + if ub is not None: + rows.append(RowEntry(con, 1)) + rhs.append(ub - offset) + con_data.append(np.array(list(repn.linear.values()))) + con_index.append( + np.array(list(map(var_order.__getitem__, repn.linear))) + ) + con_index_ptr.append(con_index_ptr[-1] + len(con_index[-1])) + if lb is not None: + rows.append(RowEntry(con, -1)) + rhs.append(offset - lb) + con_data.append(np.array(list(map(neg, repn.linear.values())))) + con_index.append( + np.array(list(map(var_order.__getitem__, repn.linear))) + ) + con_index_ptr.append(con_index_ptr[-1] + len(con_index[-1])) + + if with_debug_timing: + # report the last constraint + timer.toc('Constraint %s', last_parent, level=logging.DEBUG) + + # Get the variable list + columns = list(var_map.values()) + # Convert the compiled data to scipy sparse matricies + c = scipy.sparse.csr_array( + (obj_data, obj_index, obj_index_ptr), [len(obj_index_ptr) - 1, len(columns)] + ).tocsc() + A = scipy.sparse.csr_array( + (np.concatenate(con_data), np.concatenate(con_index), con_index_ptr), + [len(rows), len(columns)], + ).tocsc() + + # Some variables in the var_map may not actually have been + # written out to the LP file (e.g., added from col_order, or + # multiplied by 0 in the expressions). The easiest way to check + # for empty columns is to convert from CSR to CSC and then look + # 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) + 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]) + c = scipy.sparse.csc_array( + (c.data, c.indices, c.indptr.take(active_var_idx)), [c.shape[0], nCol] + ) + active_var_idx[-1] = A.indptr[-1] + A = scipy.sparse.csc_array( + (A.data, A.indices, A.indptr.take(active_var_idx)), [A.shape[0], nCol] + ) + + if self.config.nonnegative_vars: + c, A, columns, eliminated_vars = _csc_to_nonnegative_vars(c, A, columns) + else: + eliminated_vars = [] + + info = LinearStandardFormInfo(c, A, rhs, rows, columns, eliminated_vars) + timer.toc("Generated linear standard form representation", delta=False) + return info + + +def _csc_to_nonnegative_vars(c, A, columns): + eliminated_vars = [] + new_columns = [] + new_c_data = [] + new_c_indices = [] + new_c_indptr = [0] + new_A_data = [] + new_A_indices = [] + new_A_indptr = [0] + for i, v in enumerate(columns): + lb, ub = v.bounds + if lb is None or lb < 0: + name = v.name + new_columns.append( + Var( + name=f'_neg_{i}', + domain=v.domain, + bounds=(0, None if lb is None else -lb), + ) + ) + new_columns[-1].construct() + s, e = A.indptr[i : i + 2] + new_A_data.append(-A.data[s:e]) + new_A_indices.append(A.indices[s:e]) + new_A_indptr.append(new_A_indptr[-1] + e - s) + s, e = c.indptr[i : i + 2] + new_c_data.append(-c.data[s:e]) + new_c_indices.append(c.indices[s:e]) + new_c_indptr.append(new_c_indptr[-1] + e - s) + if ub is None or ub > 0: + # Crosses 0; split into 2 vars + new_columns.append( + Var(name=f'_pos_{i}', domain=v.domain, bounds=(0, ub)) + ) + new_columns[-1].construct() + s, e = A.indptr[i : i + 2] + new_A_data.append(A.data[s:e]) + new_A_indices.append(A.indices[s:e]) + new_A_indptr.append(new_A_indptr[-1] + e - s) + s, e = c.indptr[i : i + 2] + new_c_data.append(c.data[s:e]) + new_c_indices.append(c.indices[s:e]) + new_c_indptr.append(new_c_indptr[-1] + e - s) + eliminated_vars.append((v, new_columns[-1] - new_columns[-2])) + else: + new_columns[-1].lb = -ub + eliminated_vars.append((v, -new_columns[-1])) + else: # lb >= 0 + new_columns.append(v) + s, e = A.indptr[i : i + 2] + new_A_data.append(A.data[s:e]) + new_A_indices.append(A.indices[s:e]) + new_A_indptr.append(new_A_indptr[-1] + e - s) + s, e = c.indptr[i : i + 2] + new_c_data.append(c.data[s:e]) + new_c_indices.append(c.indices[s:e]) + new_c_indptr.append(new_c_indptr[-1] + e - s) + + nCol = len(new_columns) + c = scipy.sparse.csc_array( + (np.concatenate(new_c_data), np.concatenate(new_c_indices), new_c_indptr), + [c.shape[0], nCol], + ) + A = scipy.sparse.csc_array( + (np.concatenate(new_A_data), np.concatenate(new_A_indices), new_A_indptr), + [A.shape[0], nCol], + ) + return c, A, new_columns, eliminated_vars diff --git a/pyomo/repn/tests/test_standard_form.py b/pyomo/repn/tests/test_standard_form.py new file mode 100644 index 00000000000..6cd95f78c5b --- /dev/null +++ b/pyomo/repn/tests/test_standard_form.py @@ -0,0 +1,267 @@ +# ___________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2022 +# 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 + +import pyomo.environ as pyo + +from pyomo.common.dependencies import numpy as np, scipy_available, numpy_available +from pyomo.common.log import LoggingIntercept +from pyomo.repn.plugins.standard_form import LinearStandardFormCompiler + +for sol in ['glpk', 'cbc', 'gurobi', 'cplex', 'xpress']: + linear_solver = pyo.SolverFactory(sol) + if linear_solver.available(): + break +else: + linear_solver = None + + +@unittest.skipUnless( + scipy_available & numpy_available, "standard_form requires scipy and numpy" +) +class TestLinearStandardFormCompiler(unittest.TestCase): + def test_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] >= 3) + m.d = pyo.Constraint(expr=m.y[1] + 4 * m.y[3] <= 5) + + 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, 0], [0, 1, 4]]))) + self.assertTrue(np.all(repn.rhs == np.array([-3, 5]))) + + def test_linear_model_row_col_order(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] >= 3) + m.d = pyo.Constraint(expr=m.y[1] + 4 * m.y[3] <= 5) + + repn = LinearStandardFormCompiler().write( + m, column_order=[m.y[3], m.y[2], m.x, m.y[1]], row_order=[m.d, m.c] + ) + + 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]))) + + def test_suffix_warning(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] >= 3) + m.d = pyo.Constraint(expr=m.y[1] + 4 * m.y[3] <= 5) + + m.dual = pyo.Suffix(direction=pyo.Suffix.EXPORT) + m.b = pyo.Block() + m.b.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT_EXPORT) + + with LoggingIntercept() as LOG: + repn = LinearStandardFormCompiler().write(m) + self.assertEqual(LOG.getvalue(), "") + + m.dual[m.c] = 5 + with LoggingIntercept() as LOG: + repn = LinearStandardFormCompiler().write(m) + self.assertEqual( + LOG.getvalue(), + "EXPORT Suffix 'dual' found on 1 block:\n" + " dual\n" + "Standard Form compiler ignores export suffixes. Skipping.\n", + ) + + m.b.dual[m.d] = 1 + with LoggingIntercept() as LOG: + repn = LinearStandardFormCompiler().write(m) + self.assertEqual( + LOG.getvalue(), + "EXPORT Suffix 'dual' found on 2 blocks:\n" + " dual\n" + " b.dual\n" + "Standard Form compiler ignores export suffixes. Skipping.\n", + ) + + def _verify_solution(self, soln, repn, eq): + # clear out any old solution + for v, val in soln: + v.value = None + for v in repn.x: + v.value = None + + x = np.array(repn.x, dtype=object) + ax = repn.A.todense() @ x + + def c_rule(m, i): + if eq: + return ax[i] == repn.b[i] + else: + return ax[i] <= repn.b[i] + + test_model = pyo.ConcreteModel() + test_model.o = pyo.Objective(expr=repn.c[[1], :].todense()[0] @ x) + test_model.c = pyo.Constraint(range(len(repn.b)), rule=c_rule) + pyo.SolverFactory('glpk').solve(test_model, tee=True) + + # Propagate any solution back to the original variables + for v, expr in repn.eliminated_vars: + v.value = pyo.value(expr) + self.assertEqual(*zip(*((v.value, val) for v, val in soln))) + + @unittest.skipIf( + linear_solver is None, 'verifying results requires a linear solver' + ) + def test_alternative_forms(self): + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.y = pyo.Var( + [0, 1, 3], bounds=lambda m, i: (-1 * (i % 2) * 5, 10 - 12 * (i // 2)) + ) + m.c = pyo.Constraint(expr=m.x + 2 * m.y[1] >= 3) + m.d = pyo.Constraint(expr=m.y[1] + 4 * m.y[3] <= 5) + m.e = pyo.Constraint(expr=pyo.inequality(-2, m.y[0] + 1 + 6 * m.y[1], 7)) + m.f = pyo.Constraint(expr=m.x + m.y[0] + 2 == 10) + m.o = pyo.Objective([1, 3], rule=lambda m, i: m.x + i * 5 * m.y[i]) + + col_order = [m.x, m.y[0], m.y[1], m.y[3]] + + m.o[1].deactivate() + pyo.SolverFactory('glpk').solve(m) + m.o[1].activate() + soln = [(v, v.value) for v in col_order] + + repn = LinearStandardFormCompiler().write(m, column_order=col_order) + + self.assertEqual( + repn.rows, [(m.c, -1), (m.d, 1), (m.e, 1), (m.e, -1), (m.f, 1), (m.f, -1)] + ) + self.assertEqual(repn.x, [m.x, m.y[0], m.y[1], m.y[3]]) + ref = np.array( + [ + [-1, 0, -2, 0], + [0, 0, 1, 4], + [0, 1, 6, 0], + [0, -1, -6, 0], + [1, 1, 0, 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, -8]))) + self.assertTrue(np.all(repn.c == np.array([[1, 0, 5, 0], [1, 0, 0, 15]]))) + self._verify_solution(soln, repn, False) + + repn = LinearStandardFormCompiler().write( + m, nonnegative_vars=True, column_order=col_order + ) + + self.assertEqual( + repn.rows, [(m.c, -1), (m.d, 1), (m.e, 1), (m.e, -1), (m.f, 1), (m.f, -1)] + ) + self.assertEqual( + list(map(str, repn.x)), + ['_neg_0', '_pos_0', 'y[0]', '_neg_2', '_pos_2', '_neg_3'], + ) + ref = np.array( + [ + [1, -1, 0, 2, -2, 0], + [0, 0, 0, -1, 1, -4], + [0, 0, 1, -6, 6, 0], + [0, 0, -1, 6, -6, 0], + [-1, 1, 1, 0, 0, 0], + [1, -1, -1, 0, 0, 0], + ] + ) + self.assertTrue(np.all(repn.A == ref)) + self.assertTrue(np.all(repn.b == np.array([-3, 5, 6, 3, 8, -8]))) + self.assertTrue( + np.all(repn.c == np.array([[-1, 1, 0, -5, 5, 0], [-1, 1, 0, 0, 0, -15]])) + ) + self._verify_solution(soln, repn, False) + + repn = LinearStandardFormCompiler().write( + m, slack_form=True, column_order=col_order + ) + + self.assertEqual(repn.rows, [(m.c, 1), (m.d, 1), (m.e, 1), (m.f, 1)]) + self.assertEqual( + list(map(str, repn.x)), + ['x', 'y[0]', 'y[1]', 'y[3]', '_slack_0', '_slack_1', '_slack_2'], + ) + self.assertEqual( + list(v.bounds for v in repn.x), + [(None, None), (0, 10), (-5, 10), (-5, -2), (None, 0), (0, None), (-9, 0)], + ) + ref = np.array( + [ + [1, 0, 2, 0, 1, 0, 0], + [0, 0, 1, 4, 0, 1, 0], + [0, 1, 6, 0, 0, 0, 1], + [1, 1, 0, 0, 0, 0, 0], + ] + ) + self.assertTrue(np.all(repn.A == ref)) + self.assertTrue(np.all(repn.b == np.array([3, 5, -3, 8]))) + self.assertTrue( + np.all(repn.c == np.array([[1, 0, 5, 0, 0, 0, 0], [1, 0, 0, 15, 0, 0, 0]])) + ) + self._verify_solution(soln, repn, True) + + repn = LinearStandardFormCompiler().write( + m, slack_form=True, nonnegative_vars=True, column_order=col_order + ) + + self.assertEqual(repn.rows, [(m.c, 1), (m.d, 1), (m.e, 1), (m.f, 1)]) + self.assertEqual( + list(map(str, repn.x)), + [ + '_neg_0', + '_pos_0', + 'y[0]', + '_neg_2', + '_pos_2', + '_neg_3', + '_neg_4', + '_slack_1', + '_neg_6', + ], + ) + self.assertEqual( + list(v.bounds for v in repn.x), + [ + (0, None), + (0, None), + (0, 10), + (0, 5), + (0, 10), + (2, 5), + (0, None), + (0, None), + (0, 9), + ], + ) + ref = np.array( + [ + [-1, 1, 0, -2, 2, 0, -1, 0, 0], + [0, 0, 0, -1, 1, -4, 0, 1, 0], + [0, 0, 1, -6, 6, 0, 0, 0, -1], + [-1, 1, 1, 0, 0, 0, 0, 0, 0], + ] + ) + self.assertTrue(np.all(repn.A == ref)) + self.assertTrue(np.all(repn.b == np.array([3, 5, -3, 8]))) + ref = np.array([[-1, 1, 0, -5, 5, 0, 0, 0, 0], [-1, 1, 0, 0, 0, -15, 0, 0, 0]]) + self.assertTrue(np.all(repn.c == ref)) + self._verify_solution(soln, repn, True) From 38df7032d1c521709cb0778417a90b07b594ad31 Mon Sep 17 00:00:00 2001 From: John Siirola Date: Sun, 26 Nov 2023 16:42:49 -0700 Subject: [PATCH 3/9] cleanup unneeded definitions/imports --- pyomo/repn/plugins/standard_form.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index 6e74faca7d1..13d5a005910 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -30,23 +30,18 @@ Var, Param, Expression, - SOSConstraint, SortComponents, Suffix, SymbolMap, minimize, ) -from pyomo.core.base.component import ActiveComponent -from pyomo.core.base.label import LPFileLabeler, NumericLabeler from pyomo.opt import WriterFactory from pyomo.repn.linear import LinearRepnVisitor -from pyomo.repn.quadratic import QuadraticRepnVisitor from pyomo.repn.util import ( FileDeterminism, FileDeterminism_to_SortComponents, categorize_valid_components, initialize_var_map_from_column_order, - int_float, ordered_active_constraints, ) @@ -56,9 +51,6 @@ from pyomo.network import Port logger = logging.getLogger(__name__) -inf = float('inf') -neg_inf = float('-inf') - RowEntry = collections.namedtuple('RowEntry', ['constraint', 'bound_type']) From 3eda0b50bce7563f86fcc7c6760b7b7c1ba8c168 Mon Sep 17 00:00:00 2001 From: John Siirola Date: Sun, 26 Nov 2023 16:43:32 -0700 Subject: [PATCH 4/9] Switch to np.fromiter for performance --- pyomo/repn/plugins/standard_form.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index 13d5a005910..9f40199ff84 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -393,22 +393,25 @@ def write(self, model): con_index.append(np.array(_index)) con_index_ptr.append(con_index_ptr[-1] + len(_index)) else: + N = len(repn.linear) if ub is not None: rows.append(RowEntry(con, 1)) rhs.append(ub - offset) - con_data.append(np.array(list(repn.linear.values()))) + con_data.append(np.fromiter(repn.linear.values(), float, N)) con_index.append( - np.array(list(map(var_order.__getitem__, repn.linear))) + np.fromiter(map(var_order.__getitem__, repn.linear), float, N) ) - con_index_ptr.append(con_index_ptr[-1] + len(con_index[-1])) + con_index_ptr.append(con_index_ptr[-1] + N) if lb is not None: rows.append(RowEntry(con, -1)) rhs.append(offset - lb) - con_data.append(np.array(list(map(neg, repn.linear.values())))) + con_data.append( + np.fromiter(map(neg, repn.linear.values()), float, N) + ) con_index.append( - np.array(list(map(var_order.__getitem__, repn.linear))) + np.fromiter(map(var_order.__getitem__, repn.linear), float, N) ) - con_index_ptr.append(con_index_ptr[-1] + len(con_index[-1])) + con_index_ptr.append(con_index_ptr[-1] + N) if with_debug_timing: # report the last constraint From 6eb570452e2176ba5214973cbcd96370c2376df3 Mon Sep 17 00:00:00 2001 From: John Siirola Date: Sun, 26 Nov 2023 16:44:08 -0700 Subject: [PATCH 5/9] Convert maximize to minimize --- pyomo/repn/plugins/standard_form.py | 16 +++++++++++----- pyomo/repn/tests/test_standard_form.py | 11 +++++++---- 2 files changed, 18 insertions(+), 9 deletions(-) diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index 9f40199ff84..8e5e5dbb942 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -33,7 +33,7 @@ SortComponents, Suffix, SymbolMap, - minimize, + maximize, ) from pyomo.opt import WriterFactory from pyomo.repn.linear import LinearRepnVisitor @@ -317,9 +317,14 @@ def write(self, model): f"Model objective ({obj.name}) contains nonlinear terms that " "cannot be compiled to standard (linear) form." ) - obj_data.extend(repn.linear.values()) - obj_index.extend(map(var_order.__getitem__, repn.linear)) - obj_index_ptr.append(len(obj_index)) + N = len(repn.linear) + obj_data.append(np.fromiter(repn.linear.values(), float, N)) + if obj.sense == maximize: + obj_data[-1] *= -1 + obj_index.append( + np.fromiter(map(var_order.__getitem__, repn.linear), float, N) + ) + obj_index_ptr.append(obj_index_ptr[-1] + N) if with_debug_timing: timer.toc('Objective %s', obj, level=logging.DEBUG) @@ -421,7 +426,8 @@ def write(self, model): columns = list(var_map.values()) # Convert the compiled data to scipy sparse matricies c = scipy.sparse.csr_array( - (obj_data, obj_index, obj_index_ptr), [len(obj_index_ptr) - 1, len(columns)] + (np.concatenate(obj_data), np.concatenate(obj_index), obj_index_ptr), + [len(obj_index_ptr) - 1, len(columns)], ).tocsc() A = scipy.sparse.csr_array( (np.concatenate(con_data), np.concatenate(con_index), con_index_ptr), diff --git a/pyomo/repn/tests/test_standard_form.py b/pyomo/repn/tests/test_standard_form.py index 6cd95f78c5b..b805d18b303 100644 --- a/pyomo/repn/tests/test_standard_form.py +++ b/pyomo/repn/tests/test_standard_form.py @@ -134,6 +134,7 @@ def test_alternative_forms(self): m.e = pyo.Constraint(expr=pyo.inequality(-2, m.y[0] + 1 + 6 * m.y[1], 7)) m.f = pyo.Constraint(expr=m.x + m.y[0] + 2 == 10) m.o = pyo.Objective([1, 3], rule=lambda m, i: m.x + i * 5 * m.y[i]) + m.o[1].sense = pyo.maximize col_order = [m.x, m.y[0], m.y[1], m.y[3]] @@ -160,7 +161,7 @@ def test_alternative_forms(self): ) self.assertTrue(np.all(repn.A == ref)) self.assertTrue(np.all(repn.b == np.array([-3, 5, 6, 3, 8, -8]))) - self.assertTrue(np.all(repn.c == np.array([[1, 0, 5, 0], [1, 0, 0, 15]]))) + self.assertTrue(np.all(repn.c == np.array([[-1, 0, -5, 0], [1, 0, 0, 15]]))) self._verify_solution(soln, repn, False) repn = LinearStandardFormCompiler().write( @@ -187,7 +188,7 @@ def test_alternative_forms(self): self.assertTrue(np.all(repn.A == ref)) self.assertTrue(np.all(repn.b == np.array([-3, 5, 6, 3, 8, -8]))) self.assertTrue( - np.all(repn.c == np.array([[-1, 1, 0, -5, 5, 0], [-1, 1, 0, 0, 0, -15]])) + np.all(repn.c == np.array([[1, -1, 0, 5, -5, 0], [-1, 1, 0, 0, 0, -15]])) ) self._verify_solution(soln, repn, False) @@ -215,7 +216,9 @@ def test_alternative_forms(self): self.assertTrue(np.all(repn.A == ref)) self.assertTrue(np.all(repn.b == np.array([3, 5, -3, 8]))) self.assertTrue( - np.all(repn.c == np.array([[1, 0, 5, 0, 0, 0, 0], [1, 0, 0, 15, 0, 0, 0]])) + np.all( + repn.c == np.array([[-1, 0, -5, 0, 0, 0, 0], [1, 0, 0, 15, 0, 0, 0]]) + ) ) self._verify_solution(soln, repn, True) @@ -262,6 +265,6 @@ def test_alternative_forms(self): ) self.assertTrue(np.all(repn.A == ref)) self.assertTrue(np.all(repn.b == np.array([3, 5, -3, 8]))) - ref = np.array([[-1, 1, 0, -5, 5, 0, 0, 0, 0], [-1, 1, 0, 0, 0, -15, 0, 0, 0]]) + ref = np.array([[1, -1, 0, 5, -5, 0, 0, 0, 0], [-1, 1, 0, 0, 0, -15, 0, 0, 0]]) self.assertTrue(np.all(repn.c == ref)) self._verify_solution(soln, repn, True) From f2c7f53a37adb361f6c2a074956a832e5d38ed74 Mon Sep 17 00:00:00 2001 From: John Siirola Date: Sun, 26 Nov 2023 16:44:32 -0700 Subject: [PATCH 6/9] Fix checks for solver availability --- pyomo/repn/tests/test_standard_form.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/pyomo/repn/tests/test_standard_form.py b/pyomo/repn/tests/test_standard_form.py index b805d18b303..d2c096efd79 100644 --- a/pyomo/repn/tests/test_standard_form.py +++ b/pyomo/repn/tests/test_standard_form.py @@ -20,12 +20,11 @@ for sol in ['glpk', 'cbc', 'gurobi', 'cplex', 'xpress']: linear_solver = pyo.SolverFactory(sol) - if linear_solver.available(): + if linear_solver.available(exception_flag=False): break else: linear_solver = None - @unittest.skipUnless( scipy_available & numpy_available, "standard_form requires scipy and numpy" ) @@ -113,7 +112,7 @@ def c_rule(m, i): test_model = pyo.ConcreteModel() test_model.o = pyo.Objective(expr=repn.c[[1], :].todense()[0] @ x) test_model.c = pyo.Constraint(range(len(repn.b)), rule=c_rule) - pyo.SolverFactory('glpk').solve(test_model, tee=True) + linear_solver.solve(test_model, tee=True) # Propagate any solution back to the original variables for v, expr in repn.eliminated_vars: @@ -139,7 +138,7 @@ def test_alternative_forms(self): col_order = [m.x, m.y[0], m.y[1], m.y[3]] m.o[1].deactivate() - pyo.SolverFactory('glpk').solve(m) + linear_solver.solve(m) m.o[1].activate() soln = [(v, v.value) for v in col_order] From 5ac9da1c9fa1dd850bbf56fc1f8ed4a890275266 Mon Sep 17 00:00:00 2001 From: John Siirola Date: Sun, 26 Nov 2023 17:12:49 -0700 Subject: [PATCH 7/9] NFC: apply black --- pyomo/repn/tests/test_standard_form.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pyomo/repn/tests/test_standard_form.py b/pyomo/repn/tests/test_standard_form.py index d2c096efd79..d186f28dab8 100644 --- a/pyomo/repn/tests/test_standard_form.py +++ b/pyomo/repn/tests/test_standard_form.py @@ -25,6 +25,7 @@ else: linear_solver = None + @unittest.skipUnless( scipy_available & numpy_available, "standard_form requires scipy and numpy" ) From ced74d4cbb0e11d29ead70e04079a4412942e20a Mon Sep 17 00:00:00 2001 From: John Siirola Date: Sun, 26 Nov 2023 18:42:35 -0700 Subject: [PATCH 8/9] simplify construction of numpy vectors --- pyomo/repn/plugins/standard_form.py | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index 8e5e5dbb942..2b2c81b2ca9 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -11,7 +11,7 @@ import collections import logging -from operator import attrgetter, neg +from operator import attrgetter from pyomo.common.config import ( ConfigBlock, @@ -399,23 +399,19 @@ def write(self, model): con_index_ptr.append(con_index_ptr[-1] + len(_index)) else: 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 is not None: rows.append(RowEntry(con, 1)) rhs.append(ub - offset) - con_data.append(np.fromiter(repn.linear.values(), float, N)) - con_index.append( - np.fromiter(map(var_order.__getitem__, repn.linear), float, N) - ) + 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(offset - lb) - con_data.append( - np.fromiter(map(neg, repn.linear.values()), float, N) - ) - con_index.append( - np.fromiter(map(var_order.__getitem__, repn.linear), float, N) - ) + con_data.append(-_data) + con_index.append(_index) con_index_ptr.append(con_index_ptr[-1] + N) if with_debug_timing: From a6662cd37d03f8776a4136c16840e42a9a89f1b0 Mon Sep 17 00:00:00 2001 From: John Siirola Date: Mon, 27 Nov 2023 16:32:04 -0700 Subject: [PATCH 9/9] NFC: fix errors / typos in comments and docstrings --- pyomo/repn/plugins/standard_form.py | 30 ++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index 2b2c81b2ca9..8440c1ab92d 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -81,13 +81,13 @@ class LinearStandardFormInfo(object): The list of Pyomo constraint objects corresponding to the rows in `A`. Each element in the list is a 2-tuple of (_ConstraintData, row_multiplier). The `row_multiplier` will be - +/- 1 (indicating if the row was multiplied by -1 (corresponding - to a constraint lower bound or +1 (upper bound). + +/- 1 indicating if the row was multiplied by -1 (corresponding + to a constraint lower bound) or +1 (upper bound). columns : List[_VarData] The list of Pyomo variable objects corresponding to columns in - the `A` and `c` matricies. + the `A` and `c` matrices. eliminated_vars: List[Tuple[_VarData, NumericExpression]] @@ -144,7 +144,7 @@ class LinearStandardFormCompiler(object): ConfigValue( default=False, domain=bool, - description='Print timing after writing each section of the LP file', + description='Print timing after each stage of the compilation process', ), ) CONFIG.declare( @@ -155,7 +155,7 @@ class LinearStandardFormCompiler(object): description='How much effort to ensure result is deterministic', doc=""" How much effort do we want to put into ensuring the - resulting matricies are produced deterministically: + resulting matrices are produced deterministically: NONE (0) : None ORDERED (10): rely on underlying component ordering (default) SORT_INDICES (20) : sort keys of indexed components @@ -169,8 +169,9 @@ class LinearStandardFormCompiler(object): default=None, description='Preferred constraint ordering', doc=""" - List of constraints in the order that they should appear in the - LP file. Unspecified constraints will appear at the end.""", + List of constraints in the order that they should appear in + the resulting `A` matrix. Unspecified constraints will + appear at the end.""", ), ) CONFIG.declare( @@ -180,10 +181,8 @@ class LinearStandardFormCompiler(object): description='Preferred variable ordering', doc=""" List of variables in the order that they should appear in - the LP file. Note that this is only a suggestion, as the LP - file format is row-major and the columns are inferred from - the order in which variables appear in the objective - followed by each constraint.""", + the compiled representation. Unspecified variables will be + appended to the end of this list.""", ), ) @@ -251,7 +250,8 @@ def write(self, model): if unknown: raise ValueError( "The model ('%s') contains the following active components " - "that the LP compiler does not know how to process:\n\t%s" + "that the Linear Standard Form compiler does not know how to " + "process:\n\t%s" % ( model.name, "\n\t".join( @@ -420,7 +420,7 @@ def write(self, model): # Get the variable list columns = list(var_map.values()) - # Convert the compiled data to scipy sparse matricies + # Convert the compiled data to scipy sparse matrices c = scipy.sparse.csr_array( (np.concatenate(obj_data), np.concatenate(obj_index), obj_index_ptr), [len(obj_index_ptr) - 1, len(columns)], @@ -430,8 +430,8 @@ def write(self, model): [len(rows), len(columns)], ).tocsc() - # Some variables in the var_map may not actually have been - # written out to the LP file (e.g., added from col_order, or + # Some variables in the var_map may not actually appear in the + # objective or constraints (e.g., added from col_order, or # multiplied by 0 in the expressions). The easiest way to check # for empty columns is to convert from CSR to CSC and then look # at the index pointer list (an O(num_var) operation).