From c2b42b5c99726340466bb9488c50076f494c7eba Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Wed, 13 Mar 2024 13:43:43 -0600 Subject: [PATCH 01/65] Initializing dual transform --- pyomo/core/plugins/transform/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pyomo/core/plugins/transform/__init__.py b/pyomo/core/plugins/transform/__init__.py index 21e762047ca..7bf3e3229c2 100644 --- a/pyomo/core/plugins/transform/__init__.py +++ b/pyomo/core/plugins/transform/__init__.py @@ -24,3 +24,4 @@ import pyomo.core.plugins.transform.add_slack_vars import pyomo.core.plugins.transform.scaling import pyomo.core.plugins.transform.logical_to_linear +import pyomo.core.plugins.transform.lp_dual From 8fccf370042040cfbd0927278cf27895c60ced64 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Wed, 13 Mar 2024 13:45:26 -0600 Subject: [PATCH 02/65] Whoops, adding LP dual transform and tests --- pyomo/core/plugins/transform/lp_dual.py | 65 +++++++++++++++++++++++++ pyomo/core/tests/unit/test_lp_dual.py | 65 +++++++++++++++++++++++++ 2 files changed, 130 insertions(+) create mode 100644 pyomo/core/plugins/transform/lp_dual.py create mode 100644 pyomo/core/tests/unit/test_lp_dual.py diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py new file mode 100644 index 00000000000..88182e30d7f --- /dev/null +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -0,0 +1,65 @@ +# ___________________________________________________________________________ +# +# 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. +# ___________________________________________________________________________ + +from pyomo.common.autoslots import AutoSlots +from pyomo.common.collections import ComponentMap +from pyomo.common.errors import MouseTrap +from pyomo.common.dependencies import scipy +from pyomo.core import ( + ConcreteModel, Var, Constraint, Objective, TransformationFactory, + NonPositiveReals, maximize +) +from pyomo.opt import WriterFactory + +@TransformationFactory.register( + 'core.lp_dual', 'Generate the linear programming dual of the given model') +class LinearProgrammingDual(object): + def apply_to(self, model, **options): + raise MouseTrap( + "The 'core.lp_dual' transformation does not currently implement " + "apply_to since it is a bit ambiguous what it means to take a dual " + "in place. Please use 'create_using' and do what you wish with the " + "returned model." + ) + + def create_using(self, model, ostream=None, **options): + """Take linear programming dual of a model + + Returns + ------- + ConcreteModel containing linear programming dual + + Parameters + ---------- + model: ConcreteModel + The concrete Pyomo model to take the dual of + + ostream: None + This is provided for API compatibility with other writers + and is ignored here. + + """ + std_form = WriterFactory('compile_standard_form').write(model, + nonnegative_vars=True) + dual = ConcreteModel(name="%s dual" % model.name) + A_transpose = scipy.sparse.csc_matrix.transpose(std_form.A) + rows = range(A_transpose.shape[0]) + cols = range(A_transpose.shape[1]) + dual.x = Var(cols, domain=NonPositiveReals) + dual.constraints = Constraint(rows) + for i in rows: + dual.constraints[i] = sum(A_transpose[i, j]*dual.x[j] for j in cols) <= \ + std_form.c[0, i] + + dual.obj = Objective(expr=sum(std_form.rhs[j]*dual.x[j] for j in cols), + sense=maximize) + + return dual diff --git a/pyomo/core/tests/unit/test_lp_dual.py b/pyomo/core/tests/unit/test_lp_dual.py new file mode 100644 index 00000000000..487ae01d877 --- /dev/null +++ b/pyomo/core/tests/unit/test_lp_dual.py @@ -0,0 +1,65 @@ +# ___________________________________________________________________________ +# +# 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. +# ___________________________________________________________________________ + +from pyomo.common.dependencies import scipy_available +import pyomo.common.unittest as unittest +from pyomo.environ import ( + ConcreteModel, + Constraint, + NonNegativeReals, + NonPositiveReals, + Objective, + Reals, + Suffix, + TerminationCondition, + TransformationFactory, + value, + Var, +) +from pyomo.opt import SolverFactory, WriterFactory + +from pytest import set_trace + +@unittest.skipUnless(scipy_available, "Scipy not available") +class TestLPDual(unittest.TestCase): + @unittest.skipUnless(SolverFactory('gurobi').available(exception_flag=False) and + SolverFactory('gurobi').license_is_valid(), + "Gurobi is not available") + def test_lp_dual_solve(self): + m = ConcreteModel() + m.x = Var(domain=NonNegativeReals) + m.y = Var(domain=NonPositiveReals) + m.z = Var(domain=Reals) + + m.obj = Objective(expr=m.x + 2*m.y - 3*m.z) + m.c1 = Constraint(expr=-4*m.x - 2*m.y - m.z <= -5) + m.c2 = Constraint(expr=m.x + m.y <= 3) + m.c3 = Constraint(expr=- m.y - m.z <= -4.2) + m.c4 = Constraint(expr=m.z <= 42) + m.dual = Suffix(direction=Suffix.IMPORT) + + lp_dual = TransformationFactory('core.lp_dual') + dual = lp_dual.create_using(m) + dual.dual = Suffix(direction=Suffix.IMPORT) + + opt = SolverFactory('gurobi') + results = opt.solve(m) + self.assertEqual(results.solver.termination_condition, + TerminationCondition.optimal) + results = opt.solve(dual) + self.assertEqual(results.solver.termination_condition, + TerminationCondition.optimal) + + self.assertAlmostEqual(value(m.obj), value(dual.obj)) + for idx, cons in enumerate([m.c1, m.c2, m.c3, m.c4]): + self.assertAlmostEqual(value(dual.x[idx]), value(m.dual[cons])) + # for idx, (mult, v) in enumerate([(1, m.x), (-1, m.y), (1, m.z)]): + # self.assertAlmostEqual(mult*value(v), value(dual.dual[dual_cons])) From ab58d805b1aa212e196046f9d585525be280f497 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Wed, 22 May 2024 12:41:38 -0600 Subject: [PATCH 03/65] First draft of 'mixed' LP dual as well as adding option for parameterized dual --- pyomo/core/plugins/transform/lp_dual.py | 126 +++++++++++++++++++++--- pyomo/core/tests/unit/test_lp_dual.py | 17 ++++ 2 files changed, 130 insertions(+), 13 deletions(-) diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index 88182e30d7f..c38c7c89f99 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -11,17 +11,64 @@ from pyomo.common.autoslots import AutoSlots from pyomo.common.collections import ComponentMap +from pyomo.common.config import ConfigDict, ConfigValue from pyomo.common.errors import MouseTrap from pyomo.common.dependencies import scipy from pyomo.core import ( - ConcreteModel, Var, Constraint, Objective, TransformationFactory, - NonPositiveReals, maximize + ConcreteModel, + Var, + Constraint, + Objective, + TransformationFactory, + NonNegativeReals, + NonPositiveReals, + maximize, + minimize, + Reals, ) from pyomo.opt import WriterFactory + +# ESJ: TODO: copied from FME basically, should centralize. +def var_list(x): + if x.ctype is Var: + if not x.is_indexed(): + return ComponentSet([x]) + ans = ComponentSet() + for j in x.index_set(): + ans.add(x[j]) + return ans + elif hasattr(x, '__iter__'): + ans = ComponentSet() + for i in x: + ans.update(vars_to_eliminate_list(i)) + return ans + else: + raise ValueError("Expected Var or list of Vars.\n\tReceived %s" % type(x)) + + @TransformationFactory.register( - 'core.lp_dual', 'Generate the linear programming dual of the given model') + 'core.lp_dual', 'Generate the linear programming dual of the given model' +) class LinearProgrammingDual(object): + CONFIG = ConfigDict("core.lp_dual") + CONFIG.declare( + 'parameterize_wrt', + ConfigValue( + default=None, + domain=var_list, + description="Vars to treat as data for the purposes of taking the dual", + doc=""" + Optional list of Vars to be treated as data while taking the LP dual. + + For example, if this is the dual of the inner problem in a multilevel + optimization problem, then the outer problem's Vars would be specified + in this list since they are not variables from the perspective of the + inner problem. + """, + ), + ) + def apply_to(self, model, **options): raise MouseTrap( "The 'core.lp_dual' transformation does not currently implement " @@ -30,7 +77,7 @@ def apply_to(self, model, **options): "returned model." ) - def create_using(self, model, ostream=None, **options): + def create_using(self, model, ostream=None, **kwds): """Take linear programming dual of a model Returns @@ -47,19 +94,72 @@ def create_using(self, model, ostream=None, **options): and is ignored here. """ - std_form = WriterFactory('compile_standard_form').write(model, - nonnegative_vars=True) + config = self.CONFIG(kwds.pop('options', {})) + config.set_value(kwds) + + if config.parameterize_wrt is None: + return self._take_dual(model) + + return self._take_parameterized_dual(model, config.parameterize_wrt) + + def _take_dual(self, model): + std_form = WriterFactory('compile_standard_form').write( + model, mixed_form=True, set_sense=None + ) + if len(std_form.objectives) != 1: + raise ValueError( + "Model '%s' has n o objective or multiple active objectives. Cannot " + "take dual with more than one objective!" % model.name + ) + primal_sense = std_form.objectives[0].sense + dual = ConcreteModel(name="%s dual" % model.name) A_transpose = scipy.sparse.csc_matrix.transpose(std_form.A) rows = range(A_transpose.shape[0]) cols = range(A_transpose.shape[1]) - dual.x = Var(cols, domain=NonPositiveReals) + dual.x = Var(cols, domain=NonNegativeReals) + for j, (primal_cons, ineq) in enumerate(std_form.rows): + if primal_sense is minimize and ineq == 1: + dual.x[j].domain = NonPositiveReals + elif primal_sense is maximzie and ineq == -1: + dual.x[j].domain = NonPositiveReals + from pytest import set_trace + set_trace() + dual.constraints = Constraint(rows) - for i in rows: - dual.constraints[i] = sum(A_transpose[i, j]*dual.x[j] for j in cols) <= \ - std_form.c[0, i] - - dual.obj = Objective(expr=sum(std_form.rhs[j]*dual.x[j] for j in cols), - sense=maximize) + for i, primal in enumerate(std_form.columns): + if primal_sense is minimize: + if primal.domain is NonNegativeReals: + dual.constraints[i] = ( + sum(A_transpose[i, j] * dual.x[j] for j in cols) + >= std_form.c[0, i] + ) + elif primal.domain is NonPositiveReals: + dual.constraints[i] = ( + sum(A_transpose[i, j] * dual.x[j] for j in cols) + <= std_form.c[0, i] + ) + else: + if primal.domain is NonNegativeReals: + dual.constraints[i] = ( + sum(A_transpose[i, j] * dual.x[j] for j in cols) + <= std_form.c[0, i] + ) + elif primal.domain is NonPositiveReals: + dual.constraints[i] = ( + sum(A_transpose[i, j] * dual.x[j] for j in cols) + >= std_form.c[0, i] + ) + if primal.domain is Reals: + dual.constraints[i] = ( + sum(A_transpose[i, j] * dual.x[j] for j in cols) == std_form.c[0, i] + ) + + dual.obj = Objective( + expr=sum(std_form.rhs[j] * dual.x[j] for j in cols), sense=-primal_sense + ) return dual + + def _take_parameterized_dual(self, model, wrt): + pass diff --git a/pyomo/core/tests/unit/test_lp_dual.py b/pyomo/core/tests/unit/test_lp_dual.py index 487ae01d877..c3bfe74afe4 100644 --- a/pyomo/core/tests/unit/test_lp_dual.py +++ b/pyomo/core/tests/unit/test_lp_dual.py @@ -63,3 +63,20 @@ def test_lp_dual_solve(self): self.assertAlmostEqual(value(dual.x[idx]), value(m.dual[cons])) # for idx, (mult, v) in enumerate([(1, m.x), (-1, m.y), (1, m.z)]): # self.assertAlmostEqual(mult*value(v), value(dual.dual[dual_cons])) + + + def test_lp_dual(self): + m = ConcreteModel() + m.x = Var(domain=NonNegativeReals) + m.y = Var(domain=NonPositiveReals) + m.z = Var(domain=Reals) + + m.obj = Objective(expr=m.x + 2*m.y - 3*m.z) + m.c1 = Constraint(expr=-4*m.x - 2*m.y - m.z <= -5) + m.c2 = Constraint(expr=m.x + m.y >= 3) + m.c3 = Constraint(expr=- m.y - m.z == -4.2) + m.c4 = Constraint(expr=m.z <= 42) + m.dual = Suffix(direction=Suffix.IMPORT) + + lp_dual = TransformationFactory('core.lp_dual') + dual = lp_dual.create_using(m) From 16fdff487b344b7a3ec5494ffc445a58a65c15b4 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Thu, 23 May 2024 13:18:42 -0600 Subject: [PATCH 04/65] Beginning of mixed form LP dal, no parameterized, with tests that don't work --- pyomo/core/plugins/transform/lp_dual.py | 63 +++++++++++++++++++++++-- pyomo/core/tests/unit/test_lp_dual.py | 47 ++++++++++++++++++ 2 files changed, 107 insertions(+), 3 deletions(-) diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index c38c7c89f99..59ca6c89f71 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -16,6 +16,7 @@ from pyomo.common.dependencies import scipy from pyomo.core import ( ConcreteModel, + Block, Var, Constraint, Objective, @@ -47,6 +48,18 @@ def var_list(x): raise ValueError("Expected Var or list of Vars.\n\tReceived %s" % type(x)) +class _LPDualData(AutoSlots.Mixin): + __slots__ = ('primal_var', 'dual_var', 'primal_constraint', 'dual_constraint') + def __init__(self): + self.primal_var = {} + self.dual_var = {} + self.primal_constraint = ComponentMap() + self.dual_constraint = ComponentMap() + + +Block.register_private_data_initializer(_LPDualData) + + @TransformationFactory.register( 'core.lp_dual', 'Generate the linear programming dual of the given model' ) @@ -118,13 +131,17 @@ def _take_dual(self, model): rows = range(A_transpose.shape[0]) cols = range(A_transpose.shape[1]) dual.x = Var(cols, domain=NonNegativeReals) + trans_info = model.private_data() for j, (primal_cons, ineq) in enumerate(std_form.rows): if primal_sense is minimize and ineq == 1: dual.x[j].domain = NonPositiveReals - elif primal_sense is maximzie and ineq == -1: + elif primal_sense is maximize and ineq == -1: dual.x[j].domain = NonPositiveReals - from pytest import set_trace - set_trace() + if ineq == 0: + # equality + dual.x[j].domain = Reals + trans_info.primal_constraint[dual.x[j]] = primal_cons + trans_info.dual_var[primal_cons] = dual.x[j] dual.constraints = Constraint(rows) for i, primal in enumerate(std_form.columns): @@ -154,6 +171,8 @@ def _take_dual(self, model): dual.constraints[i] = ( sum(A_transpose[i, j] * dual.x[j] for j in cols) == std_form.c[0, i] ) + trans_info.dual_constraint[primal] = dual.constraints[i] + trans_info.primal_var[dual.constraints[i]] = primal dual.obj = Objective( expr=sum(std_form.rhs[j] * dual.x[j] for j in cols), sense=-primal_sense @@ -163,3 +182,41 @@ def _take_dual(self, model): def _take_parameterized_dual(self, model, wrt): pass + + def get_primal_constraint(self, model, dual_var): + primal_constraint = model.private_data().primal_constraint + if dual_var in primal_constraint: + return primal_constraint[dual_var] + else: + raise ValueError( + "It does not appear that Var '%s' is a dual variable on model '%s'" + % (dual_var.name, model.name) + ) + + def get_dual_constraint(self, model, primal_var): + dual_constraint = model.private_data().dual_constraint + if primal_var in dual_constraint: + return dual_constraint[primal_var] + else: + raise ValueError( + "It does not appear that Var '%s' is a primal variable from model '%s'" + % (primal_var.name, model.name) + ) + + def get_primal_var(self, model, dual_constraint): + primal_var = model.private_data().primal_var + if dual_constraint in primal_var: + return primal_var[dual_constraint] + else: + raise ValueError( + "It does not appear that Constraint '%s' is a dual constraint on " + "model '%s'" % (dual_constraint.name, model.name)) + + def get_dual_var(self, model, primal_constraint): + dual_var = model.private_data().dual_var + if primal_constraint in dual_var: + return dual_var[primal_constraint] + else: + raise ValueError( + "It does not appear that Constraint '%s' is a primal constraint from " + "model '%s'" % (primal_constraint.name, model.name)) diff --git a/pyomo/core/tests/unit/test_lp_dual.py b/pyomo/core/tests/unit/test_lp_dual.py index c3bfe74afe4..1bfdcf4c6cc 100644 --- a/pyomo/core/tests/unit/test_lp_dual.py +++ b/pyomo/core/tests/unit/test_lp_dual.py @@ -80,3 +80,50 @@ def test_lp_dual(self): lp_dual = TransformationFactory('core.lp_dual') dual = lp_dual.create_using(m) + + alpha = lp_dual.get_dual_var(m.c1) + beta = lp_dual.get_dual_var(m.c2) + lamb = lp_dual.get_dual_var(m.c3) + xi = lp_dual.get_dual_var(m.c4) + + self.assertIs(lp_dual.get_primal_constraint[alpha], m.c1) + self.assertIs(lp_dual.get_primal_constraint[beta], m.c2) + self.assertIs(lp_dual.get_primal_constraint[lamb], m.c3) + self.assertIs(lp_dual.get_primal_constraint[xi], m.c4) + + dx = lp_dual.get_dual_constraint[m.x] + dy = lp_dual.get_dual_constraint[m.y] + dz = lp_dual.get_dual_constraint[m.z] + + self.assertIs(lp_dual.get_primal_var[dx], m.x) + self.assertIs(lp_dual.get_primal_var[dy], m.y) + self.assertIs(lp_dual.get_primal_var[dz], m.z) + + self.assertIsInstance(alpha, Var) + self.assertIs(alpha.domain, NonPositiveReals) + self.assertEqual(alpha.ub, 0) + self.assertIsNone(alpha.lb) + self.assertIsInstance(beta, Var) + self.assertIs(alpha.domain, NonNegativeReals) + self.assertEqual(alpha.lb, 0) + self.assertIsNone(alpha.ub) + self.assertIsInstance(lamb, Var) + self.assertIs(lamb.domain, Reals) + self.assertIsNone(lamb.ub) + self.assertIsNone(lamb.lb) + self.assertIsInstance(xi, Var) + self.assertIs(xi.domain, NonPositiveReals) + self.assertEqual(xi.ub, 0) + self.assertIsNone(xi.lb) + + self.assertIsInstance(dx, Constraint) + self.assertIsInstance(dy, Constraint) + self.assertIsInstance(dz, Constraint) + + assertExpressionsEqual( + self, + dx.expr, + -4 * alpha + beta <= 1 + ) + + # TODO: map objectives, and test them From a1cb793b2fa200e9a0352ca245564f140b17ca6f Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Wed, 12 Jun 2024 11:54:15 -0600 Subject: [PATCH 05/65] Fixing a bug with the direction of constraints being transposed, also with where the mappings are stored --- pyomo/core/plugins/transform/lp_dual.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index 59ca6c89f71..ce2f3e4f971 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -131,7 +131,7 @@ def _take_dual(self, model): rows = range(A_transpose.shape[0]) cols = range(A_transpose.shape[1]) dual.x = Var(cols, domain=NonNegativeReals) - trans_info = model.private_data() + trans_info = dual.private_data() for j, (primal_cons, ineq) in enumerate(std_form.rows): if primal_sense is minimize and ineq == 1: dual.x[j].domain = NonPositiveReals @@ -149,23 +149,23 @@ def _take_dual(self, model): if primal.domain is NonNegativeReals: dual.constraints[i] = ( sum(A_transpose[i, j] * dual.x[j] for j in cols) - >= std_form.c[0, i] + <= std_form.c[0, i] ) elif primal.domain is NonPositiveReals: dual.constraints[i] = ( sum(A_transpose[i, j] * dual.x[j] for j in cols) - <= std_form.c[0, i] + >= std_form.c[0, i] ) else: if primal.domain is NonNegativeReals: dual.constraints[i] = ( sum(A_transpose[i, j] * dual.x[j] for j in cols) - <= std_form.c[0, i] + >= std_form.c[0, i] ) elif primal.domain is NonPositiveReals: dual.constraints[i] = ( sum(A_transpose[i, j] * dual.x[j] for j in cols) - >= std_form.c[0, i] + <= std_form.c[0, i] ) if primal.domain is Reals: dual.constraints[i] = ( From eb44531fd81de0f74d6b73097b707679c6558e40 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Wed, 12 Jun 2024 12:23:31 -0600 Subject: [PATCH 06/65] Making Constraint expressions prettier and testing them --- pyomo/core/plugins/transform/lp_dual.py | 37 +++++++------ pyomo/core/tests/unit/test_lp_dual.py | 72 ++++++++++++++----------- 2 files changed, 60 insertions(+), 49 deletions(-) diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index ce2f3e4f971..3afead1736c 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -28,6 +28,7 @@ Reals, ) from pyomo.opt import WriterFactory +from pyomo.repn.standard_repn import isclose_const # ESJ: TODO: copied from FME basically, should centralize. @@ -145,32 +146,30 @@ def _take_dual(self, model): dual.constraints = Constraint(rows) for i, primal in enumerate(std_form.columns): + lhs = 0 + for j in cols: + coef = A_transpose[i, j] + if not coef: + continue + elif isclose_const(coef, 1.0): + lhs += dual.x[j] + elif isclose_const(coef, -1.0): + lhs -= dual.x[j] + else: + lhs += float(A_transpose[i, j]) * dual.x[j] + if primal_sense is minimize: if primal.domain is NonNegativeReals: - dual.constraints[i] = ( - sum(A_transpose[i, j] * dual.x[j] for j in cols) - <= std_form.c[0, i] - ) + dual.constraints[i] = lhs <= float(std_form.c[0, i]) elif primal.domain is NonPositiveReals: - dual.constraints[i] = ( - sum(A_transpose[i, j] * dual.x[j] for j in cols) - >= std_form.c[0, i] - ) + dual.constraints[i] = lhs >= float(std_form.c[0, i]) else: if primal.domain is NonNegativeReals: - dual.constraints[i] = ( - sum(A_transpose[i, j] * dual.x[j] for j in cols) - >= std_form.c[0, i] - ) + dual.constraints[i] = lhs >= float(std_form.c[0, i]) elif primal.domain is NonPositiveReals: - dual.constraints[i] = ( - sum(A_transpose[i, j] * dual.x[j] for j in cols) - <= std_form.c[0, i] - ) + dual.constraints[i] = lhs <= float(std_form.c[0, i]) if primal.domain is Reals: - dual.constraints[i] = ( - sum(A_transpose[i, j] * dual.x[j] for j in cols) == std_form.c[0, i] - ) + dual.constraints[i] = lhs == float(std_form.c[0, i]) trans_info.dual_constraint[primal] = dual.constraints[i] trans_info.primal_var[dual.constraints[i]] = primal diff --git a/pyomo/core/tests/unit/test_lp_dual.py b/pyomo/core/tests/unit/test_lp_dual.py index 1bfdcf4c6cc..0b9b53d65f4 100644 --- a/pyomo/core/tests/unit/test_lp_dual.py +++ b/pyomo/core/tests/unit/test_lp_dual.py @@ -14,6 +14,7 @@ from pyomo.environ import ( ConcreteModel, Constraint, + maximize, NonNegativeReals, NonPositiveReals, Objective, @@ -24,9 +25,9 @@ value, Var, ) +from pyomo.core.expr.compare import assertExpressionsEqual from pyomo.opt import SolverFactory, WriterFactory -from pytest import set_trace @unittest.skipUnless(scipy_available, "Scipy not available") class TestLPDual(unittest.TestCase): @@ -81,49 +82,60 @@ def test_lp_dual(self): lp_dual = TransformationFactory('core.lp_dual') dual = lp_dual.create_using(m) - alpha = lp_dual.get_dual_var(m.c1) - beta = lp_dual.get_dual_var(m.c2) - lamb = lp_dual.get_dual_var(m.c3) - xi = lp_dual.get_dual_var(m.c4) + alpha = lp_dual.get_dual_var(dual, m.c1) + beta = lp_dual.get_dual_var(dual, m.c2) + lamb = lp_dual.get_dual_var(dual, m.c3) + xi = lp_dual.get_dual_var(dual, m.c4) - self.assertIs(lp_dual.get_primal_constraint[alpha], m.c1) - self.assertIs(lp_dual.get_primal_constraint[beta], m.c2) - self.assertIs(lp_dual.get_primal_constraint[lamb], m.c3) - self.assertIs(lp_dual.get_primal_constraint[xi], m.c4) + self.assertIs(lp_dual.get_primal_constraint(dual, alpha), m.c1) + self.assertIs(lp_dual.get_primal_constraint(dual, beta), m.c2) + self.assertIs(lp_dual.get_primal_constraint(dual, lamb), m.c3) + self.assertIs(lp_dual.get_primal_constraint(dual, xi), m.c4) - dx = lp_dual.get_dual_constraint[m.x] - dy = lp_dual.get_dual_constraint[m.y] - dz = lp_dual.get_dual_constraint[m.z] + dx = lp_dual.get_dual_constraint(dual, m.x) + dy = lp_dual.get_dual_constraint(dual, m.y) + dz = lp_dual.get_dual_constraint(dual, m.z) - self.assertIs(lp_dual.get_primal_var[dx], m.x) - self.assertIs(lp_dual.get_primal_var[dy], m.y) - self.assertIs(lp_dual.get_primal_var[dz], m.z) + self.assertIs(lp_dual.get_primal_var(dual, dx), m.x) + self.assertIs(lp_dual.get_primal_var(dual, dy), m.y) + self.assertIs(lp_dual.get_primal_var(dual, dz), m.z) - self.assertIsInstance(alpha, Var) - self.assertIs(alpha.domain, NonPositiveReals) + self.assertIs(alpha.ctype, Var) + self.assertEqual(alpha.domain, NonPositiveReals) self.assertEqual(alpha.ub, 0) self.assertIsNone(alpha.lb) - self.assertIsInstance(beta, Var) - self.assertIs(alpha.domain, NonNegativeReals) - self.assertEqual(alpha.lb, 0) - self.assertIsNone(alpha.ub) - self.assertIsInstance(lamb, Var) - self.assertIs(lamb.domain, Reals) + self.assertIs(beta.ctype, Var) + self.assertEqual(beta.domain, NonNegativeReals) + self.assertEqual(beta.lb, 0) + self.assertIsNone(beta.ub) + self.assertIs(lamb.ctype, Var) + self.assertEqual(lamb.domain, Reals) self.assertIsNone(lamb.ub) self.assertIsNone(lamb.lb) - self.assertIsInstance(xi, Var) - self.assertIs(xi.domain, NonPositiveReals) + self.assertIs(xi.ctype, Var) + self.assertEqual(xi.domain, NonPositiveReals) self.assertEqual(xi.ub, 0) self.assertIsNone(xi.lb) - self.assertIsInstance(dx, Constraint) - self.assertIsInstance(dy, Constraint) - self.assertIsInstance(dz, Constraint) + self.assertIs(dx.ctype, Constraint) + self.assertIs(dy.ctype, Constraint) + self.assertIs(dz.ctype, Constraint) assertExpressionsEqual( self, dx.expr, - -4 * alpha + beta <= 1 + -4.0 * alpha + beta <= 1.0 + ) + assertExpressionsEqual( + self, dy.expr, -2.0 * alpha + beta - lamb >= 2.0 + ) + assertExpressionsEqual( + self, dz.expr, - alpha - 1.0 * lamb + xi == -3.0 ) - # TODO: map objectives, and test them + dual_obj = dual.obj + self.assertIsInstance(dual_obj, Objective) + self.assertEqual(dual_obj.sense, maximize) + assertExpressionsEqual( + self, dual_obj.expr, -5 * alpha + 3 * beta - 4.2 * lamb + 42 * xi + ) From 8420e058fd1d70bf423cc53441ee247dfd240453 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Wed, 12 Jun 2024 12:26:22 -0600 Subject: [PATCH 07/65] black --- pyomo/core/plugins/transform/lp_dual.py | 9 +++-- pyomo/core/tests/unit/test_lp_dual.py | 47 +++++++++++-------------- 2 files changed, 27 insertions(+), 29 deletions(-) diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index 3afead1736c..a8e559cffee 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -51,6 +51,7 @@ def var_list(x): class _LPDualData(AutoSlots.Mixin): __slots__ = ('primal_var', 'dual_var', 'primal_constraint', 'dual_constraint') + def __init__(self): self.primal_var = {} self.dual_var = {} @@ -143,7 +144,7 @@ def _take_dual(self, model): dual.x[j].domain = Reals trans_info.primal_constraint[dual.x[j]] = primal_cons trans_info.dual_var[primal_cons] = dual.x[j] - + dual.constraints = Constraint(rows) for i, primal in enumerate(std_form.columns): lhs = 0 @@ -209,7 +210,8 @@ def get_primal_var(self, model, dual_constraint): else: raise ValueError( "It does not appear that Constraint '%s' is a dual constraint on " - "model '%s'" % (dual_constraint.name, model.name)) + "model '%s'" % (dual_constraint.name, model.name) + ) def get_dual_var(self, model, primal_constraint): dual_var = model.private_data().dual_var @@ -218,4 +220,5 @@ def get_dual_var(self, model, primal_constraint): else: raise ValueError( "It does not appear that Constraint '%s' is a primal constraint from " - "model '%s'" % (primal_constraint.name, model.name)) + "model '%s'" % (primal_constraint.name, model.name) + ) diff --git a/pyomo/core/tests/unit/test_lp_dual.py b/pyomo/core/tests/unit/test_lp_dual.py index 0b9b53d65f4..b94d34282aa 100644 --- a/pyomo/core/tests/unit/test_lp_dual.py +++ b/pyomo/core/tests/unit/test_lp_dual.py @@ -12,7 +12,7 @@ from pyomo.common.dependencies import scipy_available import pyomo.common.unittest as unittest from pyomo.environ import ( - ConcreteModel, + ConcreteModel, Constraint, maximize, NonNegativeReals, @@ -31,19 +31,21 @@ @unittest.skipUnless(scipy_available, "Scipy not available") class TestLPDual(unittest.TestCase): - @unittest.skipUnless(SolverFactory('gurobi').available(exception_flag=False) and - SolverFactory('gurobi').license_is_valid(), - "Gurobi is not available") + @unittest.skipUnless( + SolverFactory('gurobi').available(exception_flag=False) + and SolverFactory('gurobi').license_is_valid(), + "Gurobi is not available", + ) def test_lp_dual_solve(self): m = ConcreteModel() m.x = Var(domain=NonNegativeReals) m.y = Var(domain=NonPositiveReals) m.z = Var(domain=Reals) - m.obj = Objective(expr=m.x + 2*m.y - 3*m.z) - m.c1 = Constraint(expr=-4*m.x - 2*m.y - m.z <= -5) + m.obj = Objective(expr=m.x + 2 * m.y - 3 * m.z) + m.c1 = Constraint(expr=-4 * m.x - 2 * m.y - m.z <= -5) m.c2 = Constraint(expr=m.x + m.y <= 3) - m.c3 = Constraint(expr=- m.y - m.z <= -4.2) + m.c3 = Constraint(expr=-m.y - m.z <= -4.2) m.c4 = Constraint(expr=m.z <= 42) m.dual = Suffix(direction=Suffix.IMPORT) @@ -53,11 +55,13 @@ def test_lp_dual_solve(self): opt = SolverFactory('gurobi') results = opt.solve(m) - self.assertEqual(results.solver.termination_condition, - TerminationCondition.optimal) + self.assertEqual( + results.solver.termination_condition, TerminationCondition.optimal + ) results = opt.solve(dual) - self.assertEqual(results.solver.termination_condition, - TerminationCondition.optimal) + self.assertEqual( + results.solver.termination_condition, TerminationCondition.optimal + ) self.assertAlmostEqual(value(m.obj), value(dual.obj)) for idx, cons in enumerate([m.c1, m.c2, m.c3, m.c4]): @@ -65,17 +69,16 @@ def test_lp_dual_solve(self): # for idx, (mult, v) in enumerate([(1, m.x), (-1, m.y), (1, m.z)]): # self.assertAlmostEqual(mult*value(v), value(dual.dual[dual_cons])) - def test_lp_dual(self): m = ConcreteModel() m.x = Var(domain=NonNegativeReals) m.y = Var(domain=NonPositiveReals) m.z = Var(domain=Reals) - m.obj = Objective(expr=m.x + 2*m.y - 3*m.z) - m.c1 = Constraint(expr=-4*m.x - 2*m.y - m.z <= -5) + m.obj = Objective(expr=m.x + 2 * m.y - 3 * m.z) + m.c1 = Constraint(expr=-4 * m.x - 2 * m.y - m.z <= -5) m.c2 = Constraint(expr=m.x + m.y >= 3) - m.c3 = Constraint(expr=- m.y - m.z == -4.2) + m.c3 = Constraint(expr=-m.y - m.z == -4.2) m.c4 = Constraint(expr=m.z <= 42) m.dual = Suffix(direction=Suffix.IMPORT) @@ -121,17 +124,9 @@ def test_lp_dual(self): self.assertIs(dy.ctype, Constraint) self.assertIs(dz.ctype, Constraint) - assertExpressionsEqual( - self, - dx.expr, - -4.0 * alpha + beta <= 1.0 - ) - assertExpressionsEqual( - self, dy.expr, -2.0 * alpha + beta - lamb >= 2.0 - ) - assertExpressionsEqual( - self, dz.expr, - alpha - 1.0 * lamb + xi == -3.0 - ) + assertExpressionsEqual(self, dx.expr, -4.0 * alpha + beta <= 1.0) + assertExpressionsEqual(self, dy.expr, -2.0 * alpha + beta - lamb >= 2.0) + assertExpressionsEqual(self, dz.expr, -alpha - 1.0 * lamb + xi == -3.0) dual_obj = dual.obj self.assertIsInstance(dual_obj, Objective) From d702d57b82b7103c5c9692fc6d5ed9ac43937767 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Wed, 12 Jun 2024 12:50:09 -0600 Subject: [PATCH 08/65] Adding a test that we can recover the primal from the dual with the transformation --- pyomo/core/tests/unit/test_lp_dual.py | 53 +++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/pyomo/core/tests/unit/test_lp_dual.py b/pyomo/core/tests/unit/test_lp_dual.py index b94d34282aa..eff9487f3a1 100644 --- a/pyomo/core/tests/unit/test_lp_dual.py +++ b/pyomo/core/tests/unit/test_lp_dual.py @@ -15,6 +15,7 @@ ConcreteModel, Constraint, maximize, + minimize, NonNegativeReals, NonPositiveReals, Objective, @@ -134,3 +135,55 @@ def test_lp_dual(self): assertExpressionsEqual( self, dual_obj.expr, -5 * alpha + 3 * beta - 4.2 * lamb + 42 * xi ) + + ## + # now go the other way and recover the primal + ## + + primal = lp_dual.create_using(dual) + + x = lp_dual.get_dual_var(primal, dx) + y = lp_dual.get_dual_var(primal, dy) + z = lp_dual.get_dual_var(primal, dz) + + self.assertIs(x.ctype, Var) + self.assertEqual(x.domain, NonNegativeReals) + self.assertEqual(x.lb, 0) + self.assertIsNone(x.ub) + self.assertIs(y.ctype, Var) + self.assertEqual(y.domain, NonPositiveReals) + self.assertIsNone(y.lb) + self.assertEqual(y.ub, 0) + self.assertIs(z.ctype, Var) + self.assertEqual(z.domain, Reals) + self.assertIsNone(z.lb) + self.assertIsNone(z.ub) + + self.assertIs(lp_dual.get_primal_constraint(primal, x), dx) + self.assertIs(lp_dual.get_primal_constraint(primal, y), dy) + self.assertIs(lp_dual.get_primal_constraint(primal, z), dz) + + dalpha = lp_dual.get_dual_constraint(primal, alpha) + dbeta = lp_dual.get_dual_constraint(primal, beta) + dlambda = lp_dual.get_dual_constraint(primal, lamb) + dxi = lp_dual.get_dual_constraint(primal, xi) + + self.assertIs(lp_dual.get_primal_var(primal, dalpha), alpha) + self.assertIs(lp_dual.get_primal_var(primal, dbeta), beta) + self.assertIs(lp_dual.get_primal_var(primal, dlambda), lamb) + self.assertIs(lp_dual.get_primal_var(primal, dxi), xi) + + self.assertIs(dalpha.ctype, Constraint) + self.assertIs(dbeta.ctype, Constraint) + self.assertIs(dlambda.ctype, Constraint) + self.assertIs(dxi.ctype, Constraint) + + assertExpressionsEqual(self, dalpha.expr, -4.0 * x - 2.0 * y - z <= -5.0) + assertExpressionsEqual(self, dbeta.expr, x + y >= 3.0) + assertExpressionsEqual(self, dlambda.expr, -y - z == -4.2) + assertExpressionsEqual(self, dxi.expr, z <= 42.0) + + primal_obj = primal.obj + self.assertIsInstance(primal_obj, Objective) + self.assertEqual(primal_obj.sense, minimize) + assertExpressionsEqual(self, primal_obj.expr, x + 2.0 * y - 3.0 * z) From ad0ef8e548f0e97d499aa5438a5392651b7ef297 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Wed, 12 Jun 2024 12:56:03 -0600 Subject: [PATCH 09/65] Testing that dual values reported by solver and what we get from solving transformed model match --- pyomo/core/tests/unit/test_lp_dual.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/pyomo/core/tests/unit/test_lp_dual.py b/pyomo/core/tests/unit/test_lp_dual.py index eff9487f3a1..945b1e3a600 100644 --- a/pyomo/core/tests/unit/test_lp_dual.py +++ b/pyomo/core/tests/unit/test_lp_dual.py @@ -65,10 +65,14 @@ def test_lp_dual_solve(self): ) self.assertAlmostEqual(value(m.obj), value(dual.obj)) - for idx, cons in enumerate([m.c1, m.c2, m.c3, m.c4]): - self.assertAlmostEqual(value(dual.x[idx]), value(m.dual[cons])) - # for idx, (mult, v) in enumerate([(1, m.x), (-1, m.y), (1, m.z)]): - # self.assertAlmostEqual(mult*value(v), value(dual.dual[dual_cons])) + for cons in [m.c1, m.c2, m.c3, m.c4]: + self.assertAlmostEqual( + value(lp_dual.get_dual_var(dual, cons)), value(m.dual[cons]) + ) + for v in [m.x, m.y, m.z]: + self.assertAlmostEqual( + value(v), value(dual.dual[lp_dual.get_dual_constraint(dual, v)]) + ) def test_lp_dual(self): m = ConcreteModel() From 06c87c998e0ed65f387ee7c7ba32ce1164b81e4b Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Fri, 21 Jun 2024 11:54:14 -0600 Subject: [PATCH 10/65] Not bothering to transpose primal constraint matrix because we don't really have to --- pyomo/core/plugins/transform/lp_dual.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index a8e559cffee..50a073a2be5 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -129,9 +129,11 @@ def _take_dual(self, model): primal_sense = std_form.objectives[0].sense dual = ConcreteModel(name="%s dual" % model.name) - A_transpose = scipy.sparse.csc_matrix.transpose(std_form.A) - rows = range(A_transpose.shape[0]) - cols = range(A_transpose.shape[1]) + # This is a csc matrix, so we'll skipping transposing and just work off + # of the folumns + A = std_form.A + rows = range(A.shape[1]) + cols = range(A.shape[0]) dual.x = Var(cols, domain=NonNegativeReals) trans_info = dual.private_data() for j, (primal_cons, ineq) in enumerate(std_form.rows): @@ -149,7 +151,7 @@ def _take_dual(self, model): for i, primal in enumerate(std_form.columns): lhs = 0 for j in cols: - coef = A_transpose[i, j] + coef = A[j, i] if not coef: continue elif isclose_const(coef, 1.0): @@ -157,7 +159,7 @@ def _take_dual(self, model): elif isclose_const(coef, -1.0): lhs -= dual.x[j] else: - lhs += float(A_transpose[i, j]) * dual.x[j] + lhs += float(A[j, i]) * dual.x[j] if primal_sense is minimize: if primal.domain is NonNegativeReals: From e237311f028a9706085ebad1940526d56cd40e20 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Fri, 21 Jun 2024 12:28:06 -0600 Subject: [PATCH 11/65] Centralizing the var list domain validator in pyomo/util --- .../fme/fourier_motzkin_elimination.py | 20 ++--------- pyomo/core/plugins/transform/lp_dual.py | 21 ++---------- pyomo/util/var_list_domain.py | 33 +++++++++++++++++++ 3 files changed, 37 insertions(+), 37 deletions(-) create mode 100644 pyomo/util/var_list_domain.py diff --git a/pyomo/contrib/fme/fourier_motzkin_elimination.py b/pyomo/contrib/fme/fourier_motzkin_elimination.py index 4636450c58e..cec16e656ac 100644 --- a/pyomo/contrib/fme/fourier_motzkin_elimination.py +++ b/pyomo/contrib/fme/fourier_motzkin_elimination.py @@ -30,6 +30,7 @@ from pyomo.repn.standard_repn import generate_standard_repn from pyomo.common.collections import ComponentMap, ComponentSet from pyomo.opt import TerminationCondition +from pyomo.util.var_list_domain import var_component_set import logging @@ -57,23 +58,6 @@ def _check_var_bounds_filter(constraint): return True -def vars_to_eliminate_list(x): - if isinstance(x, (Var, VarData)): - if not x.is_indexed(): - return ComponentSet([x]) - ans = ComponentSet() - for j in x.index_set(): - ans.add(x[j]) - return ans - elif hasattr(x, '__iter__'): - ans = ComponentSet() - for i in x: - ans.update(vars_to_eliminate_list(i)) - return ans - else: - raise ValueError("Expected Var or list of Vars.\n\tReceived %s" % type(x)) - - def gcd(a, b): while b != 0: a, b = b, a % b @@ -111,7 +95,7 @@ class Fourier_Motzkin_Elimination_Transformation(Transformation): 'vars_to_eliminate', ConfigValue( default=None, - domain=vars_to_eliminate_list, + domain=var_component_set, description="Continuous variable or list of continuous variables to " "project out of the model", doc=""" diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index 50a073a2be5..6654ef6961c 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -29,24 +29,7 @@ ) from pyomo.opt import WriterFactory from pyomo.repn.standard_repn import isclose_const - - -# ESJ: TODO: copied from FME basically, should centralize. -def var_list(x): - if x.ctype is Var: - if not x.is_indexed(): - return ComponentSet([x]) - ans = ComponentSet() - for j in x.index_set(): - ans.add(x[j]) - return ans - elif hasattr(x, '__iter__'): - ans = ComponentSet() - for i in x: - ans.update(vars_to_eliminate_list(i)) - return ans - else: - raise ValueError("Expected Var or list of Vars.\n\tReceived %s" % type(x)) +from pyomo.util.var_list_domain import var_component_set class _LPDualData(AutoSlots.Mixin): @@ -71,7 +54,7 @@ class LinearProgrammingDual(object): 'parameterize_wrt', ConfigValue( default=None, - domain=var_list, + domain=var_component_set, description="Vars to treat as data for the purposes of taking the dual", doc=""" Optional list of Vars to be treated as data while taking the LP dual. diff --git a/pyomo/util/var_list_domain.py b/pyomo/util/var_list_domain.py new file mode 100644 index 00000000000..f890b308d76 --- /dev/null +++ b/pyomo/util/var_list_domain.py @@ -0,0 +1,33 @@ +# ___________________________________________________________________________ +# +# 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. +# ___________________________________________________________________________ + +from pyomo.common.collections import ComponentSet +from pyomo.core import Var + +def var_component_set(x): + """ + For domain validation in ConfigDicts: Takes singletone or iterable argument 'x' + of Vars and converts it to a ComponentSet of Vars. + """ + if hasattr(x, 'ctype') and x.ctype is Var: + if not x.is_indexed(): + return ComponentSet([x]) + ans = ComponentSet() + for j in x.index_set(): + ans.add(x[j]) + return ans + elif hasattr(x, '__iter__'): + ans = ComponentSet() + for i in x: + ans.update(var_component_set(i)) + return ans + else: + raise ValueError("Expected Var or iterable of Vars.\n\tReceived %s" % type(x)) From d3fdccb30501484f94bb779d1af33fcc56c43f4e Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Fri, 21 Jun 2024 16:37:59 -0600 Subject: [PATCH 12/65] Modularizing standard form a little so that I can override the parts that assume that the data is numeric. Creating my own csr and csc classes that are almost complete --- pyomo/core/plugins/transform/lp_dual.py | 19 +++++++++--------- pyomo/core/tests/unit/test_lp_dual.py | 22 ++++++++++++++++++++- pyomo/repn/plugins/__init__.py | 1 + pyomo/repn/plugins/standard_form.py | 26 ++++++++++++++++--------- 4 files changed, 48 insertions(+), 20 deletions(-) diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index 6654ef6961c..c105cd671d6 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -96,14 +96,16 @@ def create_using(self, model, ostream=None, **kwds): config.set_value(kwds) if config.parameterize_wrt is None: - return self._take_dual(model) - - return self._take_parameterized_dual(model, config.parameterize_wrt) + std_form = WriterFactory('compile_standard_form').write( + model, mixed_form=True, set_sense=None + ) + else: + std_form = WriterFactory('parameterized_standard_form_compiler').write( + model, wrt=config.parameterize_wrt, mixed_form=True, set_sense=None + ) + return self._take_dual(model, std_form) - def _take_dual(self, model): - std_form = WriterFactory('compile_standard_form').write( - model, mixed_form=True, set_sense=None - ) + def _take_dual(self, model, std_form): if len(std_form.objectives) != 1: raise ValueError( "Model '%s' has n o objective or multiple active objectives. Cannot " @@ -165,9 +167,6 @@ def _take_dual(self, model): return dual - def _take_parameterized_dual(self, model, wrt): - pass - def get_primal_constraint(self, model, dual_var): primal_constraint = model.private_data().primal_constraint if dual_var in primal_constraint: diff --git a/pyomo/core/tests/unit/test_lp_dual.py b/pyomo/core/tests/unit/test_lp_dual.py index 945b1e3a600..733d11275e2 100644 --- a/pyomo/core/tests/unit/test_lp_dual.py +++ b/pyomo/core/tests/unit/test_lp_dual.py @@ -12,6 +12,7 @@ from pyomo.common.dependencies import scipy_available import pyomo.common.unittest as unittest from pyomo.environ import ( + Binary, ConcreteModel, Constraint, maximize, @@ -37,6 +38,7 @@ class TestLPDual(unittest.TestCase): and SolverFactory('gurobi').license_is_valid(), "Gurobi is not available", ) + def test_lp_dual_solve(self): m = ConcreteModel() m.x = Var(domain=NonNegativeReals) @@ -85,7 +87,6 @@ def test_lp_dual(self): m.c2 = Constraint(expr=m.x + m.y >= 3) m.c3 = Constraint(expr=-m.y - m.z == -4.2) m.c4 = Constraint(expr=m.z <= 42) - m.dual = Suffix(direction=Suffix.IMPORT) lp_dual = TransformationFactory('core.lp_dual') dual = lp_dual.create_using(m) @@ -191,3 +192,22 @@ def test_lp_dual(self): self.assertIsInstance(primal_obj, Objective) self.assertEqual(primal_obj.sense, minimize) assertExpressionsEqual(self, primal_obj.expr, x + 2.0 * y - 3.0 * z) + + def test_parameterized_linear_dual(self): + m = ConcreteModel() + + m.outer1 = Var(domain=Binary) + m.outer = Var([2, 3], domain=Binary) + + m.x = Var(domain=NonNegativeReals) + m.y = Var(domain=NonPositiveReals) + m.z = Var(domain=Reals) + + m.obj = Objective(expr=m.x + 2 * m.y - 3 * m.outer[3] * m.z) + m.c1 = Constraint(expr=-4 * m.x - 2 * m.y - m.z <= -5*m.outer1) + m.c2 = Constraint(expr=m.x + m.outer[2]*m.y >= 3) + m.c3 = Constraint(expr=-m.y - m.z == -4.2) + m.c4 = Constraint(expr=m.z <= 42) + + lp_dual = TransformationFactory('core.lp_dual') + dual = lp_dual.create_using(m, parameterize_wrt=[m.outer1, m.outer]) diff --git a/pyomo/repn/plugins/__init__.py b/pyomo/repn/plugins/__init__.py index d3804c55106..aa9490d8e11 100644 --- a/pyomo/repn/plugins/__init__.py +++ b/pyomo/repn/plugins/__init__.py @@ -19,6 +19,7 @@ def load(): import pyomo.repn.plugins.lp_writer import pyomo.repn.plugins.nl_writer import pyomo.repn.plugins.standard_form + import pyomo.repn.plugins.parameterized_standard_form from pyomo.opt import WriterFactory diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index e684829e2f4..bafedb618b6 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -248,6 +248,16 @@ class _LinearStandardFormCompiler_impl(object): def __init__(self, config): self.config = config + def _get_visitor(self, var_map, var_order, sorter): + return LinearRepnVisitor({}, var_map, var_order, sorter) + + def _get_data_list(self, linear_repn): + return np.fromiter(linear_repn.values(), float, len(linear_repn)) + + def _compile_matrix(self, data, index, index_ptr, nrows, ncols): + return scipy.sparse.csr_array( (data, index, index_ptr), [nrows, ncols] + ).tocsc() + def write(self, model): timing_logger = logging.getLogger('pyomo.common.timing.writer') timer = TicTocTimer(logger=timing_logger) @@ -293,7 +303,7 @@ def write(self, model): 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, sorter) + visitor = self._get_visitor(var_map, var_order, sorter) timer.toc('Initialized column order', level=logging.DEBUG) @@ -344,7 +354,7 @@ def write(self, model): "cannot be compiled to standard (linear) form." ) N = len(repn.linear) - obj_data.append(np.fromiter(repn.linear.values(), float, N)) + obj_data.append(self._get_data_list(repn.linear)) obj_offset.append(repn.constant) if set_sense is not None and set_sense != obj.sense: obj_data[-1] *= -1 @@ -405,7 +415,7 @@ def write(self, model): if mixed_form: N = len(repn.linear) - _data = np.fromiter(repn.linear.values(), float, N) + _data = self._get_data_list(repn.linear) _index = np.fromiter(map(var_order.__getitem__, repn.linear), float, N) if ub == lb: rows.append(RowEntry(con, 0)) @@ -478,15 +488,13 @@ def write(self, model): if obj_data: obj_data = np.concatenate(obj_data) obj_index = np.concatenate(obj_index) - c = scipy.sparse.csr_array( - (obj_data, obj_index, obj_index_ptr), [len(obj_index_ptr) - 1, len(columns)] - ).tocsc() + c = self._compile_matrix(obj_data, obj_index, obj_index_ptr, + len(obj_index_ptr) - 1, len(columns)) if rows: con_data = np.concatenate(con_data) con_index = np.concatenate(con_index) - A = scipy.sparse.csr_array( - (con_data, con_index, con_index_ptr), [len(rows), len(columns)] - ).tocsc() + A = self._compile_matrix(con_data, con_index, con_index_ptr, len(rows), + len(columns)) # Some variables in the var_map may not actually appear in the # objective or constraints (e.g., added from col_order, or From 73e07f01bdf06e62093b36490402e427a6167b34 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 24 Jun 2024 10:09:15 -0600 Subject: [PATCH 13/65] Draft of scipy-style CSR and CSC matrices that accomodate pyomo expressions in the data --- pyomo/repn/plugins/standard_form.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index bafedb618b6..d510cc53a0a 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -503,11 +503,14 @@ def write(self, model): # at the index pointer list (an O(num_var) operation). c_ip = c.indptr A_ip = A.indptr + print(c_ip) + print(A_ip) 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 + print(active_var_mask) augmented_mask = np.concatenate((active_var_mask, [True])) reduced_A_indptr = A.indptr[augmented_mask] nCol = len(reduced_A_indptr) - 1 From 27e09f5396357ee67959780163e923eeb3737f41 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 24 Jun 2024 10:09:36 -0600 Subject: [PATCH 14/65] Rewriting LP dual code to actually use a CSC matrix sanely --- pyomo/core/plugins/transform/lp_dual.py | 67 +++++++++++++------------ 1 file changed, 35 insertions(+), 32 deletions(-) diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index c105cd671d6..6ce2292a663 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -108,7 +108,7 @@ def create_using(self, model, ostream=None, **kwds): def _take_dual(self, model, std_form): if len(std_form.objectives) != 1: raise ValueError( - "Model '%s' has n o objective or multiple active objectives. Cannot " + "Model '%s' has no objective or multiple active objectives. Cannot " "take dual with more than one objective!" % model.name ) primal_sense = std_form.objectives[0].sense @@ -117,53 +117,56 @@ def _take_dual(self, model, std_form): # This is a csc matrix, so we'll skipping transposing and just work off # of the folumns A = std_form.A - rows = range(A.shape[1]) - cols = range(A.shape[0]) - dual.x = Var(cols, domain=NonNegativeReals) + dual_rows = range(A.shape[1]) + dual_cols = range(A.shape[0]) + dual.x = Var(dual_cols, domain=NonNegativeReals) trans_info = dual.private_data() for j, (primal_cons, ineq) in enumerate(std_form.rows): - if primal_sense is minimize and ineq == 1: + # maximize is -1 and minimize is +1 and ineq is +1 for <= and -1 for + # >=, so we need to change domain to NonPositiveReals if the product + # of these is +1. + if primal_sense * ineq == 1: dual.x[j].domain = NonPositiveReals - elif primal_sense is maximize and ineq == -1: - dual.x[j].domain = NonPositiveReals - if ineq == 0: + elif ineq == 0: # equality dual.x[j].domain = Reals trans_info.primal_constraint[dual.x[j]] = primal_cons trans_info.dual_var[primal_cons] = dual.x[j] - dual.constraints = Constraint(rows) + dual.constraints = Constraint(dual_rows) for i, primal in enumerate(std_form.columns): lhs = 0 - for j in cols: - coef = A[j, i] - if not coef: - continue - elif isclose_const(coef, 1.0): - lhs += dual.x[j] - elif isclose_const(coef, -1.0): - lhs -= dual.x[j] - else: - lhs += float(A[j, i]) * dual.x[j] - - if primal_sense is minimize: + for j in range(A.indptr[i], A.indptr[i+1]): + coef = A.data[j] + primal_row = A.indices[j] + lhs += coef * dual.x[primal_row] + + # coef = A[j, i] + # if not coef: + # continue + # elif isclose_const(coef, 1.0): + # lhs += dual.x[j] + # elif isclose_const(coef, -1.0): + # lhs -= dual.x[j] + # else: + # lhs += float(A[j, i]) * dual.x[j] + if primal.domain is Reals: + dual.constraints[i] = lhs == std_form.c[0, i] + elif primal_sense is minimize: if primal.domain is NonNegativeReals: - dual.constraints[i] = lhs <= float(std_form.c[0, i]) - elif primal.domain is NonPositiveReals: - dual.constraints[i] = lhs >= float(std_form.c[0, i]) + dual.constraints[i] = lhs <= std_form.c[0, i] + else: # primal.domain is NonPositiveReals + dual.constraints[i] = lhs >= std_form.c[0, i] else: if primal.domain is NonNegativeReals: - dual.constraints[i] = lhs >= float(std_form.c[0, i]) - elif primal.domain is NonPositiveReals: - dual.constraints[i] = lhs <= float(std_form.c[0, i]) - if primal.domain is Reals: - dual.constraints[i] = lhs == float(std_form.c[0, i]) + dual.constraints[i] = lhs >= std_form.c[0, i] + else: # primal.domain is NonPositiveReals + dual.constraints[i] = lhs <= std_form.c[0, i] trans_info.dual_constraint[primal] = dual.constraints[i] trans_info.primal_var[dual.constraints[i]] = primal - dual.obj = Objective( - expr=sum(std_form.rhs[j] * dual.x[j] for j in cols), sense=-primal_sense - ) + dual.obj = Objective( expr=sum(std_form.rhs[j] * dual.x[j] for j in + dual_cols), sense=-primal_sense ) return dual From fd276d0b12fd3f4a947cd2e40441fc9034bfbde0 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 24 Jun 2024 10:10:01 -0600 Subject: [PATCH 15/65] whoops, missing file from two commits ago --- .../plugins/parameterized_standard_form.py | 175 ++++++++++++++++++ 1 file changed, 175 insertions(+) create mode 100644 pyomo/repn/plugins/parameterized_standard_form.py diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py new file mode 100644 index 00000000000..16ffa12800c --- /dev/null +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -0,0 +1,175 @@ +# ___________________________________________________________________________ +# +# 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. +# ___________________________________________________________________________ + +from pyomo.common.config import ( + ConfigValue, + document_kwargs_from_configdict +) +from pyomo.common.dependencies import numpy as np +from pyomo.common.gc_manager import PauseGC + +from pyomo.opt import WriterFactory +from pyomo.repn.parameterized_linear import ParameterizedLinearRepnVisitor +from pyomo.repn.plugins.standard_form import ( + LinearStandardFormInfo, + LinearStandardFormCompiler, + _LinearStandardFormCompiler_impl +) +from pyomo.util.var_list_domain import var_component_set + + +@WriterFactory.register( + 'parameterized_standard_form_compiler', + 'Compile an LP to standard form (`min cTx s.t. Ax <= b`) treating some ' + 'variables as data (e.g., variables decided by the outer problem in a ' + 'bilevel optimization problem).' +) +class ParameterizedLinearStandardFormCompiler(LinearStandardFormCompiler): + CONFIG = LinearStandardFormCompiler.CONFIG() + CONFIG.declare( + 'wrt', + ConfigValue( + default=None, + domain=var_component_set, + description="Vars to treat as data for the purposes of compiling" + "the standard form", + doc=""" + Optional list of Vars to be treated as data while compiling the + standard form. + + For example, if this is the standard form of an inner problem in a + multilevel optimization problem, then the outer problem's Vars would + be specified in this list since they are not variables from the + perspective of the inner problem. + """, + ), + ) + + @document_kwargs_from_configdict(CONFIG) + def write(self, model, ostream=None, **options): + """Convert a model to standard form (`min cTx s.t. Ax <= b`) treating the + Vars specified in 'wrt' as data + + 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 _ParameterizedLinearStandardFormCompiler_impl(config).write(model) + + +class _ParameterizedLinearStandardFormCompiler_impl(_LinearStandardFormCompiler_impl): + def _get_visitor(self, var_map, var_order, sorter): + wrt = self.config.wrt + if wrt is None: + wrt = [] + return ParameterizedLinearRepnVisitor({}, var_map, var_order, sorter, wrt=wrt) + + def _get_data_list(self, linear_repn): + # override this to not attempt conversion to float since that will fail + # on the Pyomo expressions + return [v for v in linear_repn.values()] + + def _compile_matrix(self, data, index, index_ptr, nrows, ncols): + return _CSRMatrix(data, index, index_ptr, nrows, ncols).tocsc() + + +class _CSRMatrix(object): + def __init__(self, data, col_index, row_index_ptr, nrows, ncols): + # We store the indices and index pointers as numpy arrays for the sake + # of duck typing, but not the data because that can contain Pyomo + # expressions, so we just use a list. + self.data = data + self.indices = np.array(col_index) + self.indptr = np.array(row_index_ptr) + self.shape = (nrows, ncols) + + def tocsc(self): + # Implements the same algorithm as scipy's csr_tocsc function + csr_data = self.data + col_index = self.indices + row_index_ptr = self.indptr + nrows = self.shape[0] + + num_nonzeros = len(csr_data) + csc_data = [None for x in csr_data] + row_index = np.empty(num_nonzeros)#[None for i in range(num_nonzeros)] + # tally the nonzeros in each column + col_index_ptr = np.zeros(num_nonzeros)#[0 for i in range(num_nonzeros)] + for i in col_index: + col_index_ptr[int(i)] += 1 + + # cumulative sum the tally to get the column index pointer + cum_sum = 0 + for i, tally in enumerate(col_index_ptr): + col_index_ptr[i] = cum_sum + cum_sum += tally + # we leave off the last entry (the number of nonzeros) because we are + # going to do the cute scipy thing and it will get 'added' at the end + # when we shift right (by which I mean it will conveniently already be + # there) + + # Now we are actually going to mess up what we just did while we + # construct the row index, but it's beautiful because "mess up" just + # means we increment everything by one, so we can fix it at the end. + for row in range(nrows): + for j in range(row_index_ptr[row], row_index_ptr[row + 1]): + col = int(col_index[j]) + dest = int(col_index_ptr[col]) + row_index[dest] = row + # Note that the data changes order because now we are looking + # for nonzeros through the columns rather than through the rows. + csc_data[dest] = csr_data[j] + + col_index_ptr[col] += 1 + + # fix the column index pointer by inserting 0 at the beginning--the rest + # is right because we know each entry got incremented by 1 in the loop + # above. + col_index_ptr = np.insert(col_index_ptr, 0, 0) + + return _CSCMatrix(csc_data, row_index, col_index_ptr, *self.shape) + + +class _CSCMatrix(object): + def __init__(self, data, row_index, col_index_ptr, nrows, ncols): + # We store the indices and index pointers as numpy arrays for the sake + # of duck typing, but not the data because that can contain Pyomo + # expressions, so we just use a list. + self.data = data + self.indices = np.array(row_index) + self.indptr = np.array(col_index_ptr) + self.shape = (nrows, ncols) + + +if __name__ == '__main__': + A = _CSRMatrix([5, 8, 3, 6], [0, 1, 2, 1], [0, 1, 2, 3, 4], 4, 4) + + thing = A.tocsc() + + print(thing.data) + print(thing.indices) + print(thing.indptr) From ff3808003907452f855080834e89efaa4c3883a0 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Tue, 25 Jun 2024 09:42:23 -0600 Subject: [PATCH 16/65] Implementing todense, making the data structures in CSC and CSR numpy arrays --- .../plugins/parameterized_standard_form.py | 76 +++++++++++++------ pyomo/repn/plugins/standard_form.py | 28 +++---- 2 files changed, 65 insertions(+), 39 deletions(-) diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index 16ffa12800c..86f1147f226 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -91,22 +91,26 @@ def _get_visitor(self, var_map, var_order, sorter): def _get_data_list(self, linear_repn): # override this to not attempt conversion to float since that will fail # on the Pyomo expressions - return [v for v in linear_repn.values()] + return np.array([v for v in linear_repn.values()]) - def _compile_matrix(self, data, index, index_ptr, nrows, ncols): + def _csc_matrix_from_csr(self, data, index, index_ptr, nrows, ncols): return _CSRMatrix(data, index, index_ptr, nrows, ncols).tocsc() + def _csc_matrix(self, data, index, index_ptr, nrows, ncols): + return _CSCMatrix(data, index, index_ptr, nrows, ncols) -class _CSRMatrix(object): - def __init__(self, data, col_index, row_index_ptr, nrows, ncols): - # We store the indices and index pointers as numpy arrays for the sake - # of duck typing, but not the data because that can contain Pyomo - # expressions, so we just use a list. - self.data = data - self.indices = np.array(col_index) - self.indptr = np.array(row_index_ptr) +class _SparseMatrixBase(object): + def __init__(self, data, indices, indptr, nrows, ncols): + self.data = np.array(data) + self.indices = np.array(indices, dtype=int) + self.indptr = np.array(indptr, dtype=int) self.shape = (nrows, ncols) + def __eq__(self, other): + return self.todense() == other + + +class _CSRMatrix(_SparseMatrixBase): def tocsc(self): # Implements the same algorithm as scipy's csr_tocsc function csr_data = self.data @@ -115,10 +119,10 @@ def tocsc(self): nrows = self.shape[0] num_nonzeros = len(csr_data) - csc_data = [None for x in csr_data] - row_index = np.empty(num_nonzeros)#[None for i in range(num_nonzeros)] + csc_data = np.empty(csr_data.shape[0], dtype=object) + row_index = np.empty(num_nonzeros, dtype=int) # tally the nonzeros in each column - col_index_ptr = np.zeros(num_nonzeros)#[0 for i in range(num_nonzeros)] + col_index_ptr = np.zeros(self.shape[1], dtype=int) for i in col_index: col_index_ptr[int(i)] += 1 @@ -137,8 +141,8 @@ def tocsc(self): # means we increment everything by one, so we can fix it at the end. for row in range(nrows): for j in range(row_index_ptr[row], row_index_ptr[row + 1]): - col = int(col_index[j]) - dest = int(col_index_ptr[col]) + col = col_index[j] + dest = col_index_ptr[col] row_index[dest] = row # Note that the data changes order because now we are looking # for nonzeros through the columns rather than through the rows. @@ -151,18 +155,37 @@ def tocsc(self): # above. col_index_ptr = np.insert(col_index_ptr, 0, 0) - return _CSCMatrix(csc_data, row_index, col_index_ptr, *self.shape) + return _CSCMatrix(csc_data, row_index, col_index_ptr, *self.shape) + def todense(self): + nrows = self.shape[0] + col_index = self.indices + row_index_ptr = self.indptr + data = self.data -class _CSCMatrix(object): - def __init__(self, data, row_index, col_index_ptr, nrows, ncols): - # We store the indices and index pointers as numpy arrays for the sake - # of duck typing, but not the data because that can contain Pyomo - # expressions, so we just use a list. - self.data = data - self.indices = np.array(row_index) - self.indptr = np.array(col_index_ptr) - self.shape = (nrows, ncols) + dense = np.zeros(self.shape, dtype=object) + + for row in range(nrows): + for j in range(row_index_ptr[row], row_index_ptr[j + 1]): + dense[row, col_index[j]] = data[j] + + return dense + + +class _CSCMatrix(_SparseMatrixBase): + def todense(self): + ncols = self.shape[1] + row_index = self.indices + col_index_ptr = self.indptr + data = self.data + + dense = np.zeros(self.shape, dtype=object) + + for col in range(ncols): + for j in range(col_index_ptr[col], col_index_ptr[col + 1]): + dense[row_index[j], col] = data[j] + + return dense if __name__ == '__main__': @@ -173,3 +196,6 @@ def __init__(self, data, row_index, col_index_ptr, nrows, ncols): print(thing.data) print(thing.indices) print(thing.indptr) + + dense = thing.todense() + print(dense) diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index d510cc53a0a..792a295e71c 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -254,10 +254,13 @@ def _get_visitor(self, var_map, var_order, sorter): def _get_data_list(self, linear_repn): return np.fromiter(linear_repn.values(), float, len(linear_repn)) - def _compile_matrix(self, data, index, index_ptr, nrows, ncols): + def _csc_matrix_from_csr(self, data, index, index_ptr, nrows, ncols): return scipy.sparse.csr_array( (data, index, index_ptr), [nrows, ncols] ).tocsc() + def _csc_matrix(self, data, index, index_ptr, nrows, ncols): + return scipy.sparse.csc_array((data, index, index_ptr), [nrows, ncols]) + def write(self, model): timing_logger = logging.getLogger('pyomo.common.timing.writer') timer = TicTocTimer(logger=timing_logger) @@ -463,7 +466,7 @@ 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) + _data = self._get_data_list(repn.linear) _index = np.fromiter(map(var_order.__getitem__, repn.linear), float, N) if ub is not None: rows.append(RowEntry(con, 1)) @@ -488,13 +491,13 @@ def write(self, model): if obj_data: obj_data = np.concatenate(obj_data) obj_index = np.concatenate(obj_index) - c = self._compile_matrix(obj_data, obj_index, obj_index_ptr, - len(obj_index_ptr) - 1, len(columns)) + c = self._csc_matrix_from_csr(obj_data, obj_index, obj_index_ptr, + len(obj_index_ptr) - 1, len(columns)) if rows: con_data = np.concatenate(con_data) con_index = np.concatenate(con_index) - A = self._compile_matrix(con_data, con_index, con_index_ptr, len(rows), - len(columns)) + A = self._csc_matrix_from_csr(con_data, con_index, con_index_ptr, len(rows), + len(columns)) # Some variables in the var_map may not actually appear in the # objective or constraints (e.g., added from col_order, or @@ -503,28 +506,25 @@ def write(self, model): # at the index pointer list (an O(num_var) operation). c_ip = c.indptr A_ip = A.indptr - print(c_ip) - print(A_ip) 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 - print(active_var_mask) 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): columns = [v for k, v in zip(active_var_mask, columns) if k] - c = scipy.sparse.csc_array( - (c.data, c.indices, c.indptr[augmented_mask]), [c.shape[0], nCol] - ) + c = self._csc_matrix(c.data, c.indices, c.indptr[augmented_mask], + c.shape[0], nCol) # active_var_idx[-1] = len(columns) - A = scipy.sparse.csc_array( - (A.data, A.indices, reduced_A_indptr), [A.shape[0], nCol] + A = self._csc_matrix( + A.data, A.indices, reduced_A_indptr, A.shape[0], nCol ) if self.config.nonnegative_vars: + # ESJ TODO: Need to rewrite this for parameterized too c, A, columns, eliminated_vars = _csc_to_nonnegative_vars(c, A, columns) else: eliminated_vars = [] From 8d582c907acb36e3c91332b91bd5994e84d452ca Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Tue, 25 Jun 2024 09:42:42 -0600 Subject: [PATCH 17/65] Start on unit tests for parameterized standard form --- .../tests/test_parameterized_standard_form.py | 129 ++++++++++++++++++ 1 file changed, 129 insertions(+) create mode 100644 pyomo/repn/tests/test_parameterized_standard_form.py diff --git a/pyomo/repn/tests/test_parameterized_standard_form.py b/pyomo/repn/tests/test_parameterized_standard_form.py new file mode 100644 index 00000000000..ec107503acd --- /dev/null +++ b/pyomo/repn/tests/test_parameterized_standard_form.py @@ -0,0 +1,129 @@ +# ___________________________________________________________________________ +# +# 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. +# ___________________________________________________________________________ +# + +from pyomo.common.dependencies import numpy as np, scipy_available, numpy_available +import pyomo.common.unittest as unittest + +from pyomo.environ import ( + ConcreteModel, + Constraint, + Var, +) +from pyomo.core.expr import ( + MonomialTermExpression, + NegationExpression, + ProductExpression +) +from pyomo.core.expr.compare import assertExpressionsEqual + +from pyomo.repn.plugins.parameterized_standard_form import ( + ParameterizedLinearStandardFormCompiler, + _CSRMatrix, + _CSCMatrix, +) + +@unittest.skipUnless( + numpy_available & scipy_available, "CSC and CSR representations require " + "scipy and numpy" +) +class TestSparseMatrixRepresentations(unittest.TestCase): + def test_csr_to_csc_only_data(self): + A = _CSRMatrix([5, 8, 3, 6], [0, 1, 2, 1], [0, 1, 2, 3, 4], 4, 4) + thing = A.tocsc() + + self.assertTrue(np.all(thing.data == np.array([5, 8, 6, 3]))) + self.assertTrue(np.all(thing.indices == np.array([0, 1, 3, 2]))) + self.assertTrue(np.all(thing.indptr == np.array([0, 1, 3, 4, 4]))) + + def test_csr_to_csc_pyomo_exprs(self): + m = ConcreteModel() + m.x = Var() + m.y = Var() + + A = _CSRMatrix([5, 8 * m.x, 3 * m.x * m.y ** 2, 6], [0, 1, 2, 1], + [0, 1, 2, 3, 4], 4, 4) + thing = A.tocsc() + + self.assertEqual(thing.data[0], 5) + assertExpressionsEqual( + self, thing.data[1], 8 * m.x) + self.assertEqual(thing.data[2], 6) + assertExpressionsEqual( + self, thing.data[3], 3 * m.x * m.y ** 2) + self.assertEqual(thing.data.shape, (4,)) + + self.assertTrue(np.all(thing.indices == np.array([0, 1, 3, 2]))) + self.assertTrue(np.all(thing.indptr == np.array([0, 1, 3, 4, 4]))) + + def test_csr_to_csc_empty_matrix(self): + A = _CSRMatrix([], [], [0], 0, 4) + thing = A.tocsc() + + self.assertEqual(thing.data.size, 0) + self.assertEqual(thing.indices.size, 0) + self.assertEqual(thing.shape, (0, 4)) + self.assertTrue(np.all(thing.indptr == np.zeros(5))) + + +def assertExpressionArraysEqual(self, A, B): + self.assertEqual(A.shape, B.shape) + for i in range(A.shape[0]): + for j in range(A.shape[1]): + assertExpressionsEqual(self, A[i, j], B[i, j]) + + +def assertExpressionListsEqual(self, A, B): + self.assertEqual(len(A), len(B)) + for i, a in enumerate(A): + assertExpressionsEqual(self, a, B[i]) + + +@unittest.skipUnless( + numpy_available & scipy_available, "Parameterized standard form requires " + "scipy and numpy" +) +class TestParameterizedStandardFormCompiler(unittest.TestCase): + def test_linear_model(self): + m = ConcreteModel() + m.x = Var() + m.y = Var([1, 2, 3]) + m.c = Constraint(expr=m.x + 2 * m.y[1] >= 3) + m.d = Constraint(expr=m.y[1] + 4 * m.y[3] <= 5) + + repn = ParameterizedLinearStandardFormCompiler().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]))) + self.assertEqual(repn.rows, [(m.c, -1), (m.d, 1)]) + self.assertEqual(repn.columns, [m.x, m.y[1], m.y[3]]) + + def test_parameterized_linear_model(self): + m = ConcreteModel() + m.x = Var() + m.y = Var([1, 2, 3]) + m.data = Var([1, 2]) + m.more_data = Var() + m.c = Constraint(expr=m.x + 2 * m.data[1] * m.data[2] * m.y[1] >= 3) + m.d = Constraint(expr=m.y[1] + 4 * m.y[3] <= 5 * m.more_data) + + repn = ParameterizedLinearStandardFormCompiler().write(m, wrt=[m.data, + m.more_data]) + + self.assertTrue(np.all(repn.c == np.array([0, 0, 0]))) + assertExpressionArraysEqual(self, repn.A.todense(), + np.array([[-1, NegationExpression((ProductExpression([MonomialTermExpression([2, m.data[1]]), m.data[2]]),)), 0], + [0, 1, 4]])) + assertExpressionListsEqual(self, repn.rhs, + [-3, 5 * m.more_data]) + self.assertEqual(repn.rows, [(m.c, -1), (m.d, 1)]) + self.assertEqual(repn.columns, [m.x, m.y[1], m.y[3]]) From 4c32954831d3575a7c0969d9750bafed5d2344cd Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Tue, 25 Jun 2024 09:44:57 -0600 Subject: [PATCH 18/65] Blacking --- .../plugins/parameterized_standard_form.py | 16 +++-- pyomo/repn/plugins/standard_form.py | 22 +++---- .../tests/test_parameterized_standard_form.py | 61 +++++++++++-------- 3 files changed, 55 insertions(+), 44 deletions(-) diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index 86f1147f226..bd4faeb2d14 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -9,10 +9,7 @@ # This software is distributed under the 3-clause BSD License. # ___________________________________________________________________________ -from pyomo.common.config import ( - ConfigValue, - document_kwargs_from_configdict -) +from pyomo.common.config import ConfigValue, document_kwargs_from_configdict from pyomo.common.dependencies import numpy as np from pyomo.common.gc_manager import PauseGC @@ -21,7 +18,7 @@ from pyomo.repn.plugins.standard_form import ( LinearStandardFormInfo, LinearStandardFormCompiler, - _LinearStandardFormCompiler_impl + _LinearStandardFormCompiler_impl, ) from pyomo.util.var_list_domain import var_component_set @@ -30,7 +27,7 @@ 'parameterized_standard_form_compiler', 'Compile an LP to standard form (`min cTx s.t. Ax <= b`) treating some ' 'variables as data (e.g., variables decided by the outer problem in a ' - 'bilevel optimization problem).' + 'bilevel optimization problem).', ) class ParameterizedLinearStandardFormCompiler(LinearStandardFormCompiler): CONFIG = LinearStandardFormCompiler.CONFIG() @@ -99,6 +96,7 @@ def _csc_matrix_from_csr(self, data, index, index_ptr, nrows, ncols): def _csc_matrix(self, data, index, index_ptr, nrows, ncols): return _CSCMatrix(data, index, index_ptr, nrows, ncols) + class _SparseMatrixBase(object): def __init__(self, data, indices, indptr, nrows, ncols): self.data = np.array(data) @@ -149,13 +147,13 @@ def tocsc(self): csc_data[dest] = csr_data[j] col_index_ptr[col] += 1 - + # fix the column index pointer by inserting 0 at the beginning--the rest # is right because we know each entry got incremented by 1 in the loop # above. col_index_ptr = np.insert(col_index_ptr, 0, 0) - return _CSCMatrix(csc_data, row_index, col_index_ptr, *self.shape) + return _CSCMatrix(csc_data, row_index, col_index_ptr, *self.shape) def todense(self): nrows = self.shape[0] @@ -184,7 +182,7 @@ def todense(self): for col in range(ncols): for j in range(col_index_ptr[col], col_index_ptr[col + 1]): dense[row_index[j], col] = data[j] - + return dense diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index 792a295e71c..78c76f54d60 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -255,8 +255,7 @@ def _get_data_list(self, linear_repn): return np.fromiter(linear_repn.values(), float, len(linear_repn)) def _csc_matrix_from_csr(self, data, index, index_ptr, nrows, ncols): - return scipy.sparse.csr_array( (data, index, index_ptr), [nrows, ncols] - ).tocsc() + return scipy.sparse.csr_array((data, index, index_ptr), [nrows, ncols]).tocsc() def _csc_matrix(self, data, index, index_ptr, nrows, ncols): return scipy.sparse.csc_array((data, index, index_ptr), [nrows, ncols]) @@ -491,13 +490,15 @@ def write(self, model): if obj_data: obj_data = np.concatenate(obj_data) obj_index = np.concatenate(obj_index) - c = self._csc_matrix_from_csr(obj_data, obj_index, obj_index_ptr, - len(obj_index_ptr) - 1, len(columns)) + c = self._csc_matrix_from_csr( + obj_data, obj_index, obj_index_ptr, len(obj_index_ptr) - 1, len(columns) + ) if rows: con_data = np.concatenate(con_data) con_index = np.concatenate(con_index) - A = self._csc_matrix_from_csr(con_data, con_index, con_index_ptr, len(rows), - len(columns)) + A = self._csc_matrix_from_csr( + con_data, con_index, con_index_ptr, len(rows), len(columns) + ) # Some variables in the var_map may not actually appear in the # objective or constraints (e.g., added from col_order, or @@ -516,12 +517,11 @@ def write(self, model): nCol = len(reduced_A_indptr) - 1 if nCol != len(columns): columns = [v for k, v in zip(active_var_mask, columns) if k] - c = self._csc_matrix(c.data, c.indices, c.indptr[augmented_mask], - c.shape[0], nCol) - # active_var_idx[-1] = len(columns) - A = self._csc_matrix( - A.data, A.indices, reduced_A_indptr, A.shape[0], nCol + c = self._csc_matrix( + c.data, c.indices, c.indptr[augmented_mask], c.shape[0], nCol ) + # active_var_idx[-1] = len(columns) + A = self._csc_matrix(A.data, A.indices, reduced_A_indptr, A.shape[0], nCol) if self.config.nonnegative_vars: # ESJ TODO: Need to rewrite this for parameterized too diff --git a/pyomo/repn/tests/test_parameterized_standard_form.py b/pyomo/repn/tests/test_parameterized_standard_form.py index ec107503acd..c3153562c11 100644 --- a/pyomo/repn/tests/test_parameterized_standard_form.py +++ b/pyomo/repn/tests/test_parameterized_standard_form.py @@ -13,15 +13,11 @@ from pyomo.common.dependencies import numpy as np, scipy_available, numpy_available import pyomo.common.unittest as unittest -from pyomo.environ import ( - ConcreteModel, - Constraint, - Var, -) +from pyomo.environ import ConcreteModel, Constraint, Var from pyomo.core.expr import ( MonomialTermExpression, NegationExpression, - ProductExpression + ProductExpression, ) from pyomo.core.expr.compare import assertExpressionsEqual @@ -31,9 +27,10 @@ _CSCMatrix, ) + @unittest.skipUnless( - numpy_available & scipy_available, "CSC and CSR representations require " - "scipy and numpy" + numpy_available & scipy_available, + "CSC and CSR representations require scipy and numpy", ) class TestSparseMatrixRepresentations(unittest.TestCase): def test_csr_to_csc_only_data(self): @@ -49,16 +46,15 @@ def test_csr_to_csc_pyomo_exprs(self): m.x = Var() m.y = Var() - A = _CSRMatrix([5, 8 * m.x, 3 * m.x * m.y ** 2, 6], [0, 1, 2, 1], - [0, 1, 2, 3, 4], 4, 4) + A = _CSRMatrix( + [5, 8 * m.x, 3 * m.x * m.y**2, 6], [0, 1, 2, 1], [0, 1, 2, 3, 4], 4, 4 + ) thing = A.tocsc() self.assertEqual(thing.data[0], 5) - assertExpressionsEqual( - self, thing.data[1], 8 * m.x) + assertExpressionsEqual(self, thing.data[1], 8 * m.x) self.assertEqual(thing.data[2], 6) - assertExpressionsEqual( - self, thing.data[3], 3 * m.x * m.y ** 2) + assertExpressionsEqual(self, thing.data[3], 3 * m.x * m.y**2) self.assertEqual(thing.data.shape, (4,)) self.assertTrue(np.all(thing.indices == np.array([0, 1, 3, 2]))) @@ -72,7 +68,7 @@ def test_csr_to_csc_empty_matrix(self): self.assertEqual(thing.indices.size, 0) self.assertEqual(thing.shape, (0, 4)) self.assertTrue(np.all(thing.indptr == np.zeros(5))) - + def assertExpressionArraysEqual(self, A, B): self.assertEqual(A.shape, B.shape) @@ -88,8 +84,8 @@ def assertExpressionListsEqual(self, A, B): @unittest.skipUnless( - numpy_available & scipy_available, "Parameterized standard form requires " - "scipy and numpy" + numpy_available & scipy_available, + "Parameterized standard form requires scipy and numpy", ) class TestParameterizedStandardFormCompiler(unittest.TestCase): def test_linear_model(self): @@ -116,14 +112,31 @@ def test_parameterized_linear_model(self): m.c = Constraint(expr=m.x + 2 * m.data[1] * m.data[2] * m.y[1] >= 3) m.d = Constraint(expr=m.y[1] + 4 * m.y[3] <= 5 * m.more_data) - repn = ParameterizedLinearStandardFormCompiler().write(m, wrt=[m.data, - m.more_data]) + repn = ParameterizedLinearStandardFormCompiler().write( + m, wrt=[m.data, m.more_data] + ) self.assertTrue(np.all(repn.c == np.array([0, 0, 0]))) - assertExpressionArraysEqual(self, repn.A.todense(), - np.array([[-1, NegationExpression((ProductExpression([MonomialTermExpression([2, m.data[1]]), m.data[2]]),)), 0], - [0, 1, 4]])) - assertExpressionListsEqual(self, repn.rhs, - [-3, 5 * m.more_data]) + assertExpressionArraysEqual( + self, + repn.A.todense(), + np.array( + [ + [ + -1, + NegationExpression( + ( + ProductExpression( + [MonomialTermExpression([2, m.data[1]]), m.data[2]] + ), + ) + ), + 0, + ], + [0, 1, 4], + ] + ), + ) + assertExpressionListsEqual(self, repn.rhs, [-3, 5 * m.more_data]) self.assertEqual(repn.rows, [(m.c, -1), (m.d, 1)]) self.assertEqual(repn.columns, [m.x, m.y[1], m.y[3]]) From 8f876846e078056e0307b4b3a29e972758ca3bba Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Tue, 25 Jun 2024 10:14:27 -0600 Subject: [PATCH 19/65] More unit tests (and whitespace) --- .../tests/test_parameterized_standard_form.py | 69 +++++++++++++++++++ 1 file changed, 69 insertions(+) diff --git a/pyomo/repn/tests/test_parameterized_standard_form.py b/pyomo/repn/tests/test_parameterized_standard_form.py index c3153562c11..fada879e180 100644 --- a/pyomo/repn/tests/test_parameterized_standard_form.py +++ b/pyomo/repn/tests/test_parameterized_standard_form.py @@ -140,3 +140,72 @@ def test_parameterized_linear_model(self): assertExpressionListsEqual(self, repn.rhs, [-3, 5 * m.more_data]) self.assertEqual(repn.rows, [(m.c, -1), (m.d, 1)]) self.assertEqual(repn.columns, [m.x, m.y[1], m.y[3]]) + + def test_parameterized_almost_dense_linear_model(self): + m = ConcreteModel() + m.x = Var() + m.y = Var([1, 2, 3]) + m.data = Var([1, 2]) + m.more_data = Var() + m.c = Constraint( + expr=m.x + 2 * m.y[1] + 4 * m.y[3] + m.more_data >= 10 * m.data[1] ** 2 + ) + m.d = Constraint(expr=5 * m.x + 6 * m.y[1] + 8 * m.data[2] * m.y[3] <= 20) + + repn = ParameterizedLinearStandardFormCompiler().write( + m, wrt=[m.data, m.more_data] + ) + + self.assertTrue(np.all(repn.c == np.array([0, 0, 0]))) + # m.c gets interpretted as a <= Constraint, and you can't really blame + # pyomo for that because it's not parameterized yet. So that's why this + # differs from the test in test_standard_form.py + assertExpressionArraysEqual( + self, repn.A.todense(), np.array([[-1, -2, -4], [5, 6, 8 * m.data[2]]]) + ) + assertExpressionListsEqual( + self, repn.rhs, [-(10 * m.data[1] ** 2 - m.more_data), 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_parameterized_linear_model_row_col_order(self): + m = ConcreteModel() + m.x = Var() + m.y = Var([1, 2, 3]) + m.data = Var([1, 2]) + m.more_data = Var() + m.c = Constraint(expr=m.x + 2 * m.data[1] * m.data[2] * m.y[1] >= 3) + m.d = Constraint(expr=m.y[1] + 4 * m.y[3] <= 5 * m.more_data) + + repn = ParameterizedLinearStandardFormCompiler().write( + m, + wrt=[m.data, m.more_data], + 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]))) + assertExpressionArraysEqual( + self, + repn.A.todense(), + np.array( + [ + [4, 0, 1], + [ + 0, + -1, + NegationExpression( + ( + ProductExpression( + [MonomialTermExpression([2, m.data[1]]), m.data[2]] + ), + ) + ), + ], + ] + ), + ) + assertExpressionListsEqual(self, repn.rhs, np.array([5 * m.more_data, -3])) + self.assertEqual(repn.rows, [(m.d, 1), (m.c, -1)]) + self.assertEqual(repn.columns, [m.y[3], m.x, m.y[1]]) From ced3faa46654e4bf7a9f00369726c5a14975e2ca Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Tue, 25 Jun 2024 10:23:51 -0600 Subject: [PATCH 20/65] Adding a test for todense and removing my hacky debugging in parameterized_standard_form itself. --- pyomo/repn/plugins/parameterized_standard_form.py | 15 +-------------- .../tests/test_parameterized_standard_form.py | 8 ++++++++ 2 files changed, 9 insertions(+), 14 deletions(-) diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index bd4faeb2d14..6dd187a76ef 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -164,7 +164,7 @@ def todense(self): dense = np.zeros(self.shape, dtype=object) for row in range(nrows): - for j in range(row_index_ptr[row], row_index_ptr[j + 1]): + for j in range(row_index_ptr[row], row_index_ptr[row + 1]): dense[row, col_index[j]] = data[j] return dense @@ -184,16 +184,3 @@ def todense(self): dense[row_index[j], col] = data[j] return dense - - -if __name__ == '__main__': - A = _CSRMatrix([5, 8, 3, 6], [0, 1, 2, 1], [0, 1, 2, 3, 4], 4, 4) - - thing = A.tocsc() - - print(thing.data) - print(thing.indices) - print(thing.indptr) - - dense = thing.todense() - print(dense) diff --git a/pyomo/repn/tests/test_parameterized_standard_form.py b/pyomo/repn/tests/test_parameterized_standard_form.py index fada879e180..d9b345a4626 100644 --- a/pyomo/repn/tests/test_parameterized_standard_form.py +++ b/pyomo/repn/tests/test_parameterized_standard_form.py @@ -69,6 +69,13 @@ def test_csr_to_csc_empty_matrix(self): self.assertEqual(thing.shape, (0, 4)) self.assertTrue(np.all(thing.indptr == np.zeros(5))) + def test_todense(self): + A = _CSRMatrix([5, 8, 3, 6], [0, 1, 2, 1], [0, 1, 2, 3, 4], 4, 4) + dense = np.array([[5, 0, 0, 0], [0, 8, 0, 0], [0, 0, 3, 0], [0, 6, 0, 0]]) + + self.assertTrue(np.all(A.todense() == dense)) + self.assertTrue(np.all(A.tocsc().todense() == dense)) + def assertExpressionArraysEqual(self, A, B): self.assertEqual(A.shape, B.shape) @@ -209,3 +216,4 @@ def test_parameterized_linear_model_row_col_order(self): assertExpressionListsEqual(self, repn.rhs, np.array([5 * m.more_data, -3])) self.assertEqual(repn.rows, [(m.d, 1), (m.c, -1)]) self.assertEqual(repn.columns, [m.y[3], m.x, m.y[1]]) + From 0c9567a44e2dc9a2f55a858b105d4afc06c978aa Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Tue, 25 Jun 2024 11:33:38 -0600 Subject: [PATCH 21/65] Making _csc_to_nonnegative_vars a method on the class so that it can use the correct converstion to a CSC matrix depending on the circumstances, testing the conversion to nonnegative vars --- pyomo/repn/plugins/standard_form.py | 119 +++++++++--------- .../tests/test_parameterized_standard_form.py | 66 +++++++++- 2 files changed, 124 insertions(+), 61 deletions(-) diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index 78c76f54d60..ebdcd298f4a 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -524,8 +524,8 @@ def write(self, model): A = self._csc_matrix(A.data, A.indices, reduced_A_indptr, A.shape[0], nCol) if self.config.nonnegative_vars: - # ESJ TODO: Need to rewrite this for parameterized too - c, A, columns, eliminated_vars = _csc_to_nonnegative_vars(c, A, columns) + c, A, columns, eliminated_vars = self._csc_to_nonnegative_vars(c, A, + columns) else: eliminated_vars = [] @@ -535,43 +535,56 @@ def write(self, model): 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 + def _csc_to_nonnegative_vars(self, 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'_pos_{i}', domain=v.domain, bounds=(0, ub)) + 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) @@ -579,28 +592,14 @@ def _csc_to_nonnegative_vars(c, A, columns): 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 + + nCol = len(new_columns) + c = self._csc_matrix( + np.concatenate(new_c_data), np.concatenate(new_c_indices), new_c_indptr, + c.shape[0], nCol, + ) + A = self._csc_matrix( + 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_parameterized_standard_form.py b/pyomo/repn/tests/test_parameterized_standard_form.py index d9b345a4626..f7d11774be3 100644 --- a/pyomo/repn/tests/test_parameterized_standard_form.py +++ b/pyomo/repn/tests/test_parameterized_standard_form.py @@ -13,7 +13,14 @@ from pyomo.common.dependencies import numpy as np, scipy_available, numpy_available import pyomo.common.unittest as unittest -from pyomo.environ import ConcreteModel, Constraint, Var +from pyomo.environ import ( + ConcreteModel, + Constraint, + inequality, + Objective, + maximize, + Var, +) from pyomo.core.expr import ( MonomialTermExpression, NegationExpression, @@ -76,17 +83,26 @@ def test_todense(self): self.assertTrue(np.all(A.todense() == dense)) self.assertTrue(np.all(A.tocsc().todense() == dense)) + A = _CSRMatrix([5, 6, 7, 2, 1, 1.5], [0, 1, 1, 2, 3, 1], [0, 2, 4, 5, 6], 4, 4) + dense = np.array([[5, 6, 0, 0], [0, 7, 2, 0], [0, 0, 0, 1], [0, 1.5, 0, 0]]) + self.assertTrue(np.all(A.todense() == dense)) + self.assertTrue(np.all(A.tocsc().todense() == dense)) + def assertExpressionArraysEqual(self, A, B): self.assertEqual(A.shape, B.shape) for i in range(A.shape[0]): for j in range(A.shape[1]): + print(A[i, j]) + print(B[i, j]) assertExpressionsEqual(self, A[i, j], B[i, j]) def assertExpressionListsEqual(self, A, B): self.assertEqual(len(A), len(B)) for i, a in enumerate(A): + print(a) + print(B[i]) assertExpressionsEqual(self, a, B[i]) @@ -217,3 +233,51 @@ def test_parameterized_linear_model_row_col_order(self): self.assertEqual(repn.rows, [(m.d, 1), (m.c, -1)]) self.assertEqual(repn.columns, [m.y[3], m.x, m.y[1]]) + def test_nonnegative_vars(self): + m = ConcreteModel() + m.x = Var() + m.y = Var([0, 1, 3], bounds=lambda m, i: (-1 * (i % 2) * 5, 10 - 12 * (i // 2))) + m.data = Var([1, 2]) + m.more_data = Var() + m.c = Constraint(expr=m.data[1] ** 2 * m.x + 2 * m.y[1] >= 3 * m.more_data) + m.d = Constraint(expr=m.y[1] + 4 * m.y[3] <= 5 + m.data[2]) + m.e = Constraint(expr=inequality(-2, m.y[0] + 1 + 6 * m.y[1], 7)) + m.f = Constraint(expr=m.x + (m.data[2] + m.data[1] ** 3) * m.y[0] + 2 == 10) + m.o = Objective([1, 3], rule=lambda m, i: m.x + i * 5 * m.more_data * m.y[i]) + m.o[1].sense = maximize + + col_order = [m.x, m.y[0], m.y[1], m.y[3]] + repn = ParameterizedLinearStandardFormCompiler().write( + m, wrt=[m.data, m.more_data], nonnegative_vars=True, column_order=col_order + ) + + ref = np.array( + [ + [ + NegationExpression((ProductExpression((-1, m.data[1] ** 2)),)), + ProductExpression((-1, m.data[1] ** 2)), + 0, + 2, + -2, + 0, + ], + [0, 0, 0, -1, 1, -4], + [0, 0, 1, -6, 6, 0], + [0, 0, -1, 6, -6, 0], + [-1, 1, m.data[2] + m.data[1] ** 3, 0, 0, 0], + [1, -1, -(m.data[2] + m.data[1] ** 3), 0, 0, 0], + ] + ) + assertExpressionArraysEqual(self, repn.A.todense(), ref) + assertExpressionListsEqual( + self, + repn.b, + [ + -3 * m.more_data, + NegationExpression((ProductExpression((-1, 5 + m.data[2])),)), + 6, + 3, + 8, + -8, + ], + ) From d0d9abf3bcd16afbb8665403c503c356795db515 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Fri, 28 Jun 2024 08:13:21 -0600 Subject: [PATCH 22/65] Testing all the config variations for parameterized standard form --- .../tests/test_parameterized_standard_form.py | 163 +++++++++++++++++- 1 file changed, 157 insertions(+), 6 deletions(-) diff --git a/pyomo/repn/tests/test_parameterized_standard_form.py b/pyomo/repn/tests/test_parameterized_standard_form.py index f7d11774be3..e351cda970a 100644 --- a/pyomo/repn/tests/test_parameterized_standard_form.py +++ b/pyomo/repn/tests/test_parameterized_standard_form.py @@ -93,16 +93,12 @@ def assertExpressionArraysEqual(self, A, B): self.assertEqual(A.shape, B.shape) for i in range(A.shape[0]): for j in range(A.shape[1]): - print(A[i, j]) - print(B[i, j]) assertExpressionsEqual(self, A[i, j], B[i, j]) def assertExpressionListsEqual(self, A, B): self.assertEqual(len(A), len(B)) for i, a in enumerate(A): - print(a) - print(B[i]) assertExpressionsEqual(self, a, B[i]) @@ -233,24 +229,46 @@ def test_parameterized_linear_model_row_col_order(self): self.assertEqual(repn.rows, [(m.d, 1), (m.c, -1)]) self.assertEqual(repn.columns, [m.y[3], m.x, m.y[1]]) - def test_nonnegative_vars(self): + def make_model(self, do_not_flip_c=False): m = ConcreteModel() m.x = Var() m.y = Var([0, 1, 3], bounds=lambda m, i: (-1 * (i % 2) * 5, 10 - 12 * (i // 2))) m.data = Var([1, 2]) m.more_data = Var() - m.c = Constraint(expr=m.data[1] ** 2 * m.x + 2 * m.y[1] >= 3 * m.more_data) + if do_not_flip_c: + # [ESJ: 06/24]: I should have done this sooner, but if you write c + # this way, it gets interpretted as a >= constraint, which matches + # the standard_form tests and makes life much easier. Unforuntately + # I wrote a lot of tests before I thought of this, so I'm leaving in + # both variations for the moment. + m.c = Constraint(expr=m.data[1] ** 2 * m.x + 2 * m.y[1] - 3 * m.more_data + >= 0) + else: + m.c = Constraint(expr=m.data[1] ** 2 * m.x + 2 * m.y[1] >= 3 * m.more_data) m.d = Constraint(expr=m.y[1] + 4 * m.y[3] <= 5 + m.data[2]) m.e = Constraint(expr=inequality(-2, m.y[0] + 1 + 6 * m.y[1], 7)) m.f = Constraint(expr=m.x + (m.data[2] + m.data[1] ** 3) * m.y[0] + 2 == 10) m.o = Objective([1, 3], rule=lambda m, i: m.x + i * 5 * m.more_data * m.y[i]) m.o[1].sense = maximize + return m + + def test_nonnegative_vars(self): + m = self.make_model() col_order = [m.x, m.y[0], m.y[1], m.y[3]] repn = ParameterizedLinearStandardFormCompiler().write( m, wrt=[m.data, m.more_data], nonnegative_vars=True, column_order=col_order ) + # m.c comes back opposite how it does in test_standard_form, but that's + # not unexpected. + 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( [ [ @@ -281,3 +299,136 @@ def test_nonnegative_vars(self): -8, ], ) + + c_ref = np.array( + [ + [1, -1, 0, 5 * m.more_data, -5 * m.more_data, 0], + [-1, 1, 0, 0, 0, -15 * m.more_data] + ] + ) + assertExpressionArraysEqual(self, repn.c.todense(), c_ref) + + def test_slack_form(self): + m = self.make_model() + col_order = [m.x, m.y[0], m.y[1], m.y[3]] + repn = ParameterizedLinearStandardFormCompiler().write( + m, wrt=[m.data, m.more_data], 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'], + ) + # m.c is flipped again, so the bounds on _slack_0 are flipped + self.assertEqual( + list(v.bounds for v in repn.x), + [(None, None), (0, 10), (-5, 10), (-5, -2), (0, None), (0, None), (-9, 0)], + ) + ref = np.array( + [ + [ProductExpression((-1, m.data[1] ** 2)), 0, - 2, 0, 1, 0, 0], + [0, 0, 1, 4, 0, 1, 0], + [0, 1, 6, 0, 0, 0, 1], + [1, m.data[2] + m.data[1] ** 3, 0, 0, 0, 0, 0], + ] + ) + assertExpressionArraysEqual(self, repn.A.todense(), ref) + assertExpressionListsEqual( + self, repn.b, + np.array([-3 * m.more_data, + NegationExpression((ProductExpression((-1, 5 + m.data[2])),)), -3, 8])) + c_ref = np.array( + [ + [-1, 0, -5 * m.more_data, 0, 0, 0, 0], + [1, 0, 0, 15 * m.more_data, 0, 0, 0] + ] + ) + assertExpressionArraysEqual(self, repn.c.todense(), c_ref) + + def test_mixed_form(self): + m = self.make_model() + col_order = [m.x, m.y[0], m.y[1], m.y[3]] + repn = ParameterizedLinearStandardFormCompiler().write( + m, wrt=[m.data, m.more_data], mixed_form=True, column_order=col_order + ) + + # m.c gets is opposite again + 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( + [[ProductExpression((-1, m.data[1] ** 2)), 0, - 2, 0], + [0, 0, 1, 4], + [0, 1, 6, 0], + [0, 1, 6, 0], + [1, m.data[2] + m.data[1] ** 3, 0, 0]] + ) + assertExpressionArraysEqual(self, repn.A.todense(), ref) + assertExpressionListsEqual( + self, repn.b, + np.array([- 3 * m.more_data, NegationExpression((ProductExpression((-1, 5 + m.data[2])),)), 6, -3, 8]) + ) + ref_c = np.array([ + [-1, 0, -5 * m.more_data, 0], + [1, 0, 0, 15 * m.more_data] + ]) + assertExpressionArraysEqual(self, repn.c.todense(), ref_c) + + def test_slack_form_nonnegative_vars(self): + m = self.make_model(do_not_flip_c=True) + col_order = [m.x, m.y[0], m.y[1], m.y[3]] + repn = ParameterizedLinearStandardFormCompiler().write( + m, wrt=[m.data, m.more_data], 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( + [ + [-m.data[1] ** 2, m.data[1] ** 2, 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, m.data[2] + m.data[1] ** 3, 0, 0, 0, 0, 0, 0], + ] + ) + assertExpressionArraysEqual(self, repn.A.todense(), ref) + assertExpressionListsEqual(self, repn.b, np.array([3 * m.more_data, NegationExpression((ProductExpression((-1, 5 + m.data[2])),)), -3, 8])) + c_ref = np.array([ + [1, -1, 0, 5 * m.more_data, -5 * m.more_data, 0, 0, 0, 0], + [-1, 1, 0, 0, 0, -15 * m.more_data, 0, 0, 0] + ]) + assertExpressionArraysEqual(self, repn.c.todense(), c_ref) From 39f6e400febc4232aa9b9a414b9ff462857408da Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Fri, 28 Jun 2024 08:14:20 -0600 Subject: [PATCH 23/65] Black --- pyomo/repn/plugins/standard_form.py | 19 ++-- .../tests/test_parameterized_standard_form.py | 94 ++++++++++++------- 2 files changed, 75 insertions(+), 38 deletions(-) diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index ebdcd298f4a..8d97451684e 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -524,8 +524,9 @@ def write(self, model): A = self._csc_matrix(A.data, A.indices, reduced_A_indptr, A.shape[0], nCol) if self.config.nonnegative_vars: - c, A, columns, eliminated_vars = self._csc_to_nonnegative_vars(c, A, - columns) + c, A, columns, eliminated_vars = self._csc_to_nonnegative_vars( + c, A, columns + ) else: eliminated_vars = [] @@ -595,11 +596,17 @@ def _csc_to_nonnegative_vars(self, c, A, columns): nCol = len(new_columns) c = self._csc_matrix( - np.concatenate(new_c_data), np.concatenate(new_c_indices), new_c_indptr, - c.shape[0], nCol, + np.concatenate(new_c_data), + np.concatenate(new_c_indices), + new_c_indptr, + c.shape[0], + nCol, ) A = self._csc_matrix( - np.concatenate(new_A_data), np.concatenate(new_A_indices), new_A_indptr, - A.shape[0], nCol, + 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_parameterized_standard_form.py b/pyomo/repn/tests/test_parameterized_standard_form.py index e351cda970a..aa3c4b317bc 100644 --- a/pyomo/repn/tests/test_parameterized_standard_form.py +++ b/pyomo/repn/tests/test_parameterized_standard_form.py @@ -241,8 +241,9 @@ def make_model(self, do_not_flip_c=False): # the standard_form tests and makes life much easier. Unforuntately # I wrote a lot of tests before I thought of this, so I'm leaving in # both variations for the moment. - m.c = Constraint(expr=m.data[1] ** 2 * m.x + 2 * m.y[1] - 3 * m.more_data - >= 0) + m.c = Constraint( + expr=m.data[1] ** 2 * m.x + 2 * m.y[1] - 3 * m.more_data >= 0 + ) else: m.c = Constraint(expr=m.data[1] ** 2 * m.x + 2 * m.y[1] >= 3 * m.more_data) m.d = Constraint(expr=m.y[1] + 4 * m.y[3] <= 5 + m.data[2]) @@ -303,7 +304,7 @@ def test_nonnegative_vars(self): c_ref = np.array( [ [1, -1, 0, 5 * m.more_data, -5 * m.more_data, 0], - [-1, 1, 0, 0, 0, -15 * m.more_data] + [-1, 1, 0, 0, 0, -15 * m.more_data], ] ) assertExpressionArraysEqual(self, repn.c.todense(), c_ref) @@ -312,10 +313,9 @@ def test_slack_form(self): m = self.make_model() col_order = [m.x, m.y[0], m.y[1], m.y[3]] repn = ParameterizedLinearStandardFormCompiler().write( - m, wrt=[m.data, m.more_data], slack_form=True, - column_order=col_order + m, wrt=[m.data, m.more_data], 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)), @@ -328,7 +328,7 @@ def test_slack_form(self): ) ref = np.array( [ - [ProductExpression((-1, m.data[1] ** 2)), 0, - 2, 0, 1, 0, 0], + [ProductExpression((-1, m.data[1] ** 2)), 0, -2, 0, 1, 0, 0], [0, 0, 1, 4, 0, 1, 0], [0, 1, 6, 0, 0, 0, 1], [1, m.data[2] + m.data[1] ** 3, 0, 0, 0, 0, 0], @@ -336,13 +336,21 @@ def test_slack_form(self): ) assertExpressionArraysEqual(self, repn.A.todense(), ref) assertExpressionListsEqual( - self, repn.b, - np.array([-3 * m.more_data, - NegationExpression((ProductExpression((-1, 5 + m.data[2])),)), -3, 8])) + self, + repn.b, + np.array( + [ + -3 * m.more_data, + NegationExpression((ProductExpression((-1, 5 + m.data[2])),)), + -3, + 8, + ] + ), + ) c_ref = np.array( [ [-1, 0, -5 * m.more_data, 0, 0, 0, 0], - [1, 0, 0, 15 * m.more_data, 0, 0, 0] + [1, 0, 0, 15 * m.more_data, 0, 0, 0], ] ) assertExpressionArraysEqual(self, repn.c.todense(), c_ref) @@ -355,37 +363,46 @@ def test_mixed_form(self): ) # m.c gets is opposite again - self.assertEqual( - repn.rows, [(m.c, 1), (m.d, 1), (m.e, 1), (m.e, -1), (m.f, 0)] - ) + 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( - [[ProductExpression((-1, m.data[1] ** 2)), 0, - 2, 0], - [0, 0, 1, 4], - [0, 1, 6, 0], - [0, 1, 6, 0], - [1, m.data[2] + m.data[1] ** 3, 0, 0]] + [ + [ProductExpression((-1, m.data[1] ** 2)), 0, -2, 0], + [0, 0, 1, 4], + [0, 1, 6, 0], + [0, 1, 6, 0], + [1, m.data[2] + m.data[1] ** 3, 0, 0], + ] ) assertExpressionArraysEqual(self, repn.A.todense(), ref) assertExpressionListsEqual( - self, repn.b, - np.array([- 3 * m.more_data, NegationExpression((ProductExpression((-1, 5 + m.data[2])),)), 6, -3, 8]) + self, + repn.b, + np.array( + [ + -3 * m.more_data, + NegationExpression((ProductExpression((-1, 5 + m.data[2])),)), + 6, + -3, + 8, + ] + ), ) - ref_c = np.array([ - [-1, 0, -5 * m.more_data, 0], - [1, 0, 0, 15 * m.more_data] - ]) + ref_c = np.array([[-1, 0, -5 * m.more_data, 0], [1, 0, 0, 15 * m.more_data]]) assertExpressionArraysEqual(self, repn.c.todense(), ref_c) def test_slack_form_nonnegative_vars(self): m = self.make_model(do_not_flip_c=True) col_order = [m.x, m.y[0], m.y[1], m.y[3]] repn = ParameterizedLinearStandardFormCompiler().write( - m, wrt=[m.data, m.more_data], slack_form=True, nonnegative_vars=True, - column_order=col_order + m, + wrt=[m.data, m.more_data], + 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)]) @@ -426,9 +443,22 @@ def test_slack_form_nonnegative_vars(self): ] ) assertExpressionArraysEqual(self, repn.A.todense(), ref) - assertExpressionListsEqual(self, repn.b, np.array([3 * m.more_data, NegationExpression((ProductExpression((-1, 5 + m.data[2])),)), -3, 8])) - c_ref = np.array([ - [1, -1, 0, 5 * m.more_data, -5 * m.more_data, 0, 0, 0, 0], - [-1, 1, 0, 0, 0, -15 * m.more_data, 0, 0, 0] - ]) + assertExpressionListsEqual( + self, + repn.b, + np.array( + [ + 3 * m.more_data, + NegationExpression((ProductExpression((-1, 5 + m.data[2])),)), + -3, + 8, + ] + ), + ) + c_ref = np.array( + [ + [1, -1, 0, 5 * m.more_data, -5 * m.more_data, 0, 0, 0, 0], + [-1, 1, 0, 0, 0, -15 * m.more_data, 0, 0, 0], + ] + ) assertExpressionArraysEqual(self, repn.c.todense(), c_ref) From 7173461704c32fab7f82ee3491123f45d416e752 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Tue, 2 Jul 2024 10:59:29 -0600 Subject: [PATCH 24/65] Renaming _get_data_list to _to_vector --- pyomo/repn/plugins/parameterized_standard_form.py | 8 +++++++- pyomo/repn/plugins/standard_form.py | 8 ++++---- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index 6dd187a76ef..e0c11a95b94 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -85,7 +85,7 @@ def _get_visitor(self, var_map, var_order, sorter): wrt = [] return ParameterizedLinearRepnVisitor({}, var_map, var_order, sorter, wrt=wrt) - def _get_data_list(self, linear_repn): + def _to_vector(self, linear_repn): # override this to not attempt conversion to float since that will fail # on the Pyomo expressions return np.array([v for v in linear_repn.values()]) @@ -169,6 +169,12 @@ def todense(self): return dense + def sum_duplicates(self): + pass + + def eliminate_zeros(self): + pass + class _CSCMatrix(_SparseMatrixBase): def todense(self): diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index 8d97451684e..25ff1c41fd9 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -251,7 +251,7 @@ def __init__(self, config): def _get_visitor(self, var_map, var_order, sorter): return LinearRepnVisitor({}, var_map, var_order, sorter) - def _get_data_list(self, linear_repn): + def _to_vector(self, linear_repn): return np.fromiter(linear_repn.values(), float, len(linear_repn)) def _csc_matrix_from_csr(self, data, index, index_ptr, nrows, ncols): @@ -356,7 +356,7 @@ def write(self, model): "cannot be compiled to standard (linear) form." ) N = len(repn.linear) - obj_data.append(self._get_data_list(repn.linear)) + obj_data.append(self._to_vector(repn.linear)) obj_offset.append(repn.constant) if set_sense is not None and set_sense != obj.sense: obj_data[-1] *= -1 @@ -417,7 +417,7 @@ def write(self, model): if mixed_form: N = len(repn.linear) - _data = self._get_data_list(repn.linear) + _data = self._to_vector(repn.linear) _index = np.fromiter(map(var_order.__getitem__, repn.linear), float, N) if ub == lb: rows.append(RowEntry(con, 0)) @@ -465,7 +465,7 @@ def write(self, model): con_index_ptr.append(con_index_ptr[-1] + len(_index)) else: N = len(repn.linear) - _data = self._get_data_list(repn.linear) + _data = self._to_vector(repn.linear) _index = np.fromiter(map(var_order.__getitem__, repn.linear), float, N) if ub is not None: rows.append(RowEntry(con, 1)) From 89defeb013446379862ccb2086ce1b926c71cb50 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Tue, 2 Jul 2024 15:36:39 -0600 Subject: [PATCH 25/65] Generalizing standard form to handle parameterized standard form again --- pyomo/repn/parameterized_linear.py | 2 +- .../plugins/parameterized_standard_form.py | 15 ++++---- pyomo/repn/plugins/standard_form.py | 36 ++++++++++++------- .../tests/test_parameterized_standard_form.py | 7 ++-- 4 files changed, 38 insertions(+), 22 deletions(-) diff --git a/pyomo/repn/parameterized_linear.py b/pyomo/repn/parameterized_linear.py index ebd6803e8c3..f86dfd52ba7 100644 --- a/pyomo/repn/parameterized_linear.py +++ b/pyomo/repn/parameterized_linear.py @@ -205,7 +205,7 @@ def _before_var(visitor, child): # We aren't treating this Var as a Var for the purposes of this walker return False, (_PSEUDO_CONSTANT, child) # This is a normal situation - ParameterizedLinearBeforeChildDispatcher._record_var(visitor, child) + visitor.before_child_dispatcher.record_var(visitor, child) ans = visitor.Result() ans.linear[_id] = 1 return False, (ExprType.LINEAR, ans) diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index 2131208508d..39e25954680 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -96,6 +96,9 @@ def _csc_matrix_from_csr(self, data, index, index_ptr, nrows, ncols): def _csc_matrix(self, data, index, index_ptr, nrows, ncols): return _CSCMatrix(data, index, index_ptr, nrows, ncols) + def _csr_matrix(self, data, index, index_ptr, nrows, ncols): + return _CSRMatrix(data, index, index_ptr, nrows, ncols) + class _SparseMatrixBase(object): def __init__(self, data, indices, indptr, nrows, ncols): @@ -169,12 +172,6 @@ def todense(self): return dense - def sum_duplicates(self): - pass - - def eliminate_zeros(self): - pass - class _CSCMatrix(_SparseMatrixBase): def todense(self): @@ -190,3 +187,9 @@ def todense(self): dense[row_index[j], col] = data[j] return dense + + def sum_duplicates(self): + pass + + def eliminate_zeros(self): + pass diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index 7790a3f96eb..b2c69d6c4e2 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -262,6 +262,9 @@ def _csc_matrix_from_csr(self, data, index, index_ptr, nrows, ncols): def _csc_matrix(self, data, index, index_ptr, nrows, ncols): return scipy.sparse.csc_array((data, index, index_ptr), [nrows, ncols]) + def _csr_matrix(self, data, index, index_ptr, nrows, ncols): + return scipy.sparse.csr_array((data, index, index_ptr), [nrows, ncols]) + def write(self, model): timing_logger = logging.getLogger('pyomo.common.timing.writer') timer = TicTocTimer(logger=timing_logger) @@ -572,25 +575,32 @@ def write(self, model): ) # obj_index = list(itertools.chain(*obj_index)) obj_index_ptr = np.array(obj_index_ptr, dtype=np.int32) - c = scipy.sparse.csr_array( - (obj_data, obj_index, obj_index_ptr), [len(obj_index_ptr) - 1, nCol] - ) + c = self._csr_matrix(obj_data, obj_index, obj_index_ptr, + len(obj_index_ptr) - 1, nCol) + # c = scipy.sparse.csr_array( + # (obj_data, obj_index, obj_index_ptr), [len(obj_index_ptr) - 1, nCol] + # ) + c = c.tocsc() c.sum_duplicates() c.eliminate_zeros() - c = c.tocsc() if rows: - con_data = np.fromiter( - itertools.chain.from_iterable(con_data), np.float64, con_nnz - ) + con_data = self._to_vector(itertools.chain.from_iterable(con_data), + con_nnz, np.float64) + # con_data = np.fromiter( + # itertools.chain.from_iterable(con_data), np.float64, con_nnz + # ) # con_data = list(itertools.chain(*con_data)) - con_index = np.fromiter( - itertools.chain.from_iterable(con_index), np.int32, con_nnz - ) + # con_index = np.fromiter( + # itertools.chain.from_iterable(con_index), np.int32, con_nnz + # ) + con_index = self._to_vector( + itertools.chain.from_iterable(con_index), con_nnz, np.int32) # con_index = list(itertools.chain(*con_index)) con_index_ptr = np.array(con_index_ptr, dtype=np.int32) - A = scipy.sparse.csr_array( - (con_data, con_index, con_index_ptr), [len(rows), nCol] - ) + A = self._csr_matrix(con_data, con_index, con_index_ptr, len(rows), nCol) + # A = scipy.sparse.csr_array( + # (con_data, con_index, con_index_ptr), [len(rows), nCol] + # ) A = A.tocsc() A.sum_duplicates() A.eliminate_zeros() diff --git a/pyomo/repn/tests/test_parameterized_standard_form.py b/pyomo/repn/tests/test_parameterized_standard_form.py index aa3c4b317bc..cde395ba7a4 100644 --- a/pyomo/repn/tests/test_parameterized_standard_form.py +++ b/pyomo/repn/tests/test_parameterized_standard_form.py @@ -26,7 +26,10 @@ NegationExpression, ProductExpression, ) -from pyomo.core.expr.compare import assertExpressionsEqual +from pyomo.core.expr.compare import ( + assertExpressionsEqual, + assertExpressionsStructurallyEqual +) from pyomo.repn.plugins.parameterized_standard_form import ( ParameterizedLinearStandardFormCompiler, @@ -93,7 +96,7 @@ def assertExpressionArraysEqual(self, A, B): self.assertEqual(A.shape, B.shape) for i in range(A.shape[0]): for j in range(A.shape[1]): - assertExpressionsEqual(self, A[i, j], B[i, j]) + assertExpressionsStructurallyEqual(self, A[i, j], B[i, j]) def assertExpressionListsEqual(self, A, B): From 18c848c920adef0915cf3d240227b0f023589836 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Tue, 2 Jul 2024 15:37:37 -0600 Subject: [PATCH 26/65] black --- pyomo/repn/plugins/standard_form.py | 93 ++++++++++--------- .../tests/test_parameterized_standard_form.py | 2 +- 2 files changed, 51 insertions(+), 44 deletions(-) diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index b2c69d6c4e2..986d5309d33 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -361,11 +361,11 @@ def write(self, model): offset, linear_index, linear_data, _, _ = ( template_visitor.expand_expression(obj, obj.template_expr()) ) -# <<<<<<< HEAD -# N = len(repn.linear) -# obj_data.append(self._to_vector(repn.linear)) -# obj_offset.append(repn.constant) -# ======= + # <<<<<<< HEAD + # N = len(repn.linear) + # obj_data.append(self._to_vector(repn.linear)) + # obj_offset.append(repn.constant) + # ======= N = len(linear_index) obj_index.append(map(var_order.__getitem__, linear_index)) obj_data.append(linear_data) @@ -386,10 +386,12 @@ def write(self, model): ) obj_nnz += N -#>>>>>>> templatized-writer + # >>>>>>> templatized-writer if set_sense is not None and set_sense != obj.sense: # ESJ TODO: This is gonna need _to_vector I think - obj_data[-1] = - self._to_vector(obj_data[-1], N, float) #-np.fromiter(obj_data[-1], float, N) + obj_data[-1] = -self._to_vector( + obj_data[-1], N, float + ) # -np.fromiter(obj_data[-1], float, N) obj_offset[-1] *= -1 obj_index_ptr.append(obj_index_ptr[-1] + N) if with_debug_timing: @@ -458,15 +460,15 @@ def write(self, model): ) if mixed_form: -# <<<<<<< HEAD -# N = len(repn.linear) -# _data = self._to_vector(repn.linear) -# _index = np.fromiter(map(var_order.__getitem__, repn.linear), float, N) -# if ub == lb: -# ======= + # <<<<<<< HEAD + # N = len(repn.linear) + # _data = self._to_vector(repn.linear) + # _index = np.fromiter(map(var_order.__getitem__, repn.linear), float, N) + # if ub == lb: + # ======= if lb == ub: con_nnz += N -#>>>>>>> templatized-writer + # >>>>>>> templatized-writer rows.append(RowEntry(con, 0)) rhs.append(ub - offset) con_data.append(linear_data) @@ -517,12 +519,12 @@ def write(self, model): con_index.append(linear_index) con_index_ptr.append(con_nnz) else: -# <<<<<<< HEAD -# N = len(repn.linear) -# _data = self._to_vector(repn.linear) -# _index = np.fromiter(map(var_order.__getitem__, repn.linear), float, N) -# ======= -# >>>>>>> templatized-writer + # <<<<<<< HEAD + # N = len(repn.linear) + # _data = self._to_vector(repn.linear) + # _index = np.fromiter(map(var_order.__getitem__, repn.linear), float, N) + # ======= + # >>>>>>> templatized-writer if ub is not None: if lb is not None: linear_index = list(linear_index) @@ -550,22 +552,23 @@ def write(self, model): nCol = len(columns) # Convert the compiled data to scipy sparse matrices if obj_data: -# <<<<<<< HEAD -# obj_data = np.concatenate(obj_data) -# obj_index = np.concatenate(obj_index) -# c = self._csc_matrix_from_csr( -# obj_data, obj_index, obj_index_ptr, len(obj_index_ptr) - 1, len(columns) -# ) -# if rows: -# con_data = np.concatenate(con_data) -# con_index = np.concatenate(con_index) -# A = self._csc_matrix_from_csr( -# con_data, con_index, con_index_ptr, len(rows), len(columns) -# ) -# ======= + # <<<<<<< HEAD + # obj_data = np.concatenate(obj_data) + # obj_index = np.concatenate(obj_index) + # c = self._csc_matrix_from_csr( + # obj_data, obj_index, obj_index_ptr, len(obj_index_ptr) - 1, len(columns) + # ) + # if rows: + # con_data = np.concatenate(con_data) + # con_index = np.concatenate(con_index) + # A = self._csc_matrix_from_csr( + # con_data, con_index, con_index_ptr, len(rows), len(columns) + # ) + # ======= # ESJ TODO: This is gonna need _to_vector... - obj_data = self._to_vector(itertools.chain.from_iterable(obj_data), - obj_nnz, np.float64) + obj_data = self._to_vector( + itertools.chain.from_iterable(obj_data), obj_nnz, np.float64 + ) # obj_data = np.fromiter( # itertools.chain.from_iterable(obj_data), np.float64, obj_nnz # ) @@ -575,8 +578,9 @@ def write(self, model): ) # obj_index = list(itertools.chain(*obj_index)) obj_index_ptr = np.array(obj_index_ptr, dtype=np.int32) - c = self._csr_matrix(obj_data, obj_index, obj_index_ptr, - len(obj_index_ptr) - 1, nCol) + c = self._csr_matrix( + obj_data, obj_index, obj_index_ptr, len(obj_index_ptr) - 1, nCol + ) # c = scipy.sparse.csr_array( # (obj_data, obj_index, obj_index_ptr), [len(obj_index_ptr) - 1, nCol] # ) @@ -584,8 +588,9 @@ def write(self, model): c.sum_duplicates() c.eliminate_zeros() if rows: - con_data = self._to_vector(itertools.chain.from_iterable(con_data), - con_nnz, np.float64) + con_data = self._to_vector( + itertools.chain.from_iterable(con_data), con_nnz, np.float64 + ) # con_data = np.fromiter( # itertools.chain.from_iterable(con_data), np.float64, con_nnz # ) @@ -594,7 +599,8 @@ def write(self, model): # itertools.chain.from_iterable(con_index), np.int32, con_nnz # ) con_index = self._to_vector( - itertools.chain.from_iterable(con_index), con_nnz, np.int32) + itertools.chain.from_iterable(con_index), con_nnz, np.int32 + ) # con_index = list(itertools.chain(*con_index)) con_index_ptr = np.array(con_index_ptr, dtype=np.int32) A = self._csr_matrix(con_data, con_index, con_index_ptr, len(rows), nCol) @@ -607,7 +613,7 @@ def write(self, model): if with_debug_timing: timer.toc('Formed matrices', level=logging.DEBUG) -#>>>>>>> templatized-writer + # >>>>>>> templatized-writer # Some variables in the var_map may not actually appear in the # objective or constraints (e.g., added from col_order, or @@ -630,8 +636,9 @@ def write(self, model): c.data, c.indices, c.indptr[augmented_mask], c.shape[0], len(columns) ) # active_var_idx[-1] = len(columns) - A = self._csc_matrix(A.data, A.indices, reduced_A_indptr, - A.shape[0], len(columns)) + A = self._csc_matrix( + A.data, A.indices, reduced_A_indptr, A.shape[0], len(columns) + ) if with_debug_timing: timer.toc('Eliminated %s unused columns', nCol, level=logging.DEBUG) diff --git a/pyomo/repn/tests/test_parameterized_standard_form.py b/pyomo/repn/tests/test_parameterized_standard_form.py index cde395ba7a4..5ae2889ce83 100644 --- a/pyomo/repn/tests/test_parameterized_standard_form.py +++ b/pyomo/repn/tests/test_parameterized_standard_form.py @@ -28,7 +28,7 @@ ) from pyomo.core.expr.compare import ( assertExpressionsEqual, - assertExpressionsStructurallyEqual + assertExpressionsStructurallyEqual, ) from pyomo.repn.plugins.parameterized_standard_form import ( From 6548f293ff73137fb0728654590d2dbe38254369 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Tue, 2 Jul 2024 15:39:58 -0600 Subject: [PATCH 27/65] NFC: cleaning up comments --- pyomo/repn/plugins/standard_form.py | 59 +---------------------------- 1 file changed, 1 insertion(+), 58 deletions(-) diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index 986d5309d33..21c958f2cae 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -361,11 +361,6 @@ def write(self, model): offset, linear_index, linear_data, _, _ = ( template_visitor.expand_expression(obj, obj.template_expr()) ) - # <<<<<<< HEAD - # N = len(repn.linear) - # obj_data.append(self._to_vector(repn.linear)) - # obj_offset.append(repn.constant) - # ======= N = len(linear_index) obj_index.append(map(var_order.__getitem__, linear_index)) obj_data.append(linear_data) @@ -374,8 +369,6 @@ def write(self, model): repn = visitor.walk_expression(obj.expr) N = len(repn.linear) obj_index.append(map(var_order.__getitem__, repn.linear)) - # ESJ TODO: This miiiight need a _to_vector call but I think it - # doesn't anymore? obj_data.append(repn.linear.values()) obj_offset.append(repn.constant) @@ -386,12 +379,8 @@ def write(self, model): ) obj_nnz += N - # >>>>>>> templatized-writer if set_sense is not None and set_sense != obj.sense: - # ESJ TODO: This is gonna need _to_vector I think - obj_data[-1] = -self._to_vector( - obj_data[-1], N, float - ) # -np.fromiter(obj_data[-1], float, N) + obj_data[-1] = -self._to_vector(obj_data[-1], N, float) obj_offset[-1] *= -1 obj_index_ptr.append(obj_index_ptr[-1] + N) if with_debug_timing: @@ -460,12 +449,6 @@ def write(self, model): ) if mixed_form: - # <<<<<<< HEAD - # N = len(repn.linear) - # _data = self._to_vector(repn.linear) - # _index = np.fromiter(map(var_order.__getitem__, repn.linear), float, N) - # if ub == lb: - # ======= if lb == ub: con_nnz += N # >>>>>>> templatized-writer @@ -519,12 +502,6 @@ def write(self, model): con_index.append(linear_index) con_index_ptr.append(con_nnz) else: - # <<<<<<< HEAD - # N = len(repn.linear) - # _data = self._to_vector(repn.linear) - # _index = np.fromiter(map(var_order.__getitem__, repn.linear), float, N) - # ======= - # >>>>>>> templatized-writer if ub is not None: if lb is not None: linear_index = list(linear_index) @@ -552,38 +529,16 @@ def write(self, model): nCol = len(columns) # Convert the compiled data to scipy sparse matrices if obj_data: - # <<<<<<< HEAD - # obj_data = np.concatenate(obj_data) - # obj_index = np.concatenate(obj_index) - # c = self._csc_matrix_from_csr( - # obj_data, obj_index, obj_index_ptr, len(obj_index_ptr) - 1, len(columns) - # ) - # if rows: - # con_data = np.concatenate(con_data) - # con_index = np.concatenate(con_index) - # A = self._csc_matrix_from_csr( - # con_data, con_index, con_index_ptr, len(rows), len(columns) - # ) - # ======= - # ESJ TODO: This is gonna need _to_vector... obj_data = self._to_vector( itertools.chain.from_iterable(obj_data), obj_nnz, np.float64 ) - # obj_data = np.fromiter( - # itertools.chain.from_iterable(obj_data), np.float64, obj_nnz - # ) - # obj_data = list(itertools.chain(*obj_data)) obj_index = np.fromiter( itertools.chain.from_iterable(obj_index), np.int32, obj_nnz ) - # obj_index = list(itertools.chain(*obj_index)) obj_index_ptr = np.array(obj_index_ptr, dtype=np.int32) c = self._csr_matrix( obj_data, obj_index, obj_index_ptr, len(obj_index_ptr) - 1, nCol ) - # c = scipy.sparse.csr_array( - # (obj_data, obj_index, obj_index_ptr), [len(obj_index_ptr) - 1, nCol] - # ) c = c.tocsc() c.sum_duplicates() c.eliminate_zeros() @@ -591,29 +546,17 @@ def write(self, model): con_data = self._to_vector( itertools.chain.from_iterable(con_data), con_nnz, np.float64 ) - # con_data = np.fromiter( - # itertools.chain.from_iterable(con_data), np.float64, con_nnz - # ) - # con_data = list(itertools.chain(*con_data)) - # con_index = np.fromiter( - # itertools.chain.from_iterable(con_index), np.int32, con_nnz - # ) con_index = self._to_vector( itertools.chain.from_iterable(con_index), con_nnz, np.int32 ) - # con_index = list(itertools.chain(*con_index)) con_index_ptr = np.array(con_index_ptr, dtype=np.int32) A = self._csr_matrix(con_data, con_index, con_index_ptr, len(rows), nCol) - # A = scipy.sparse.csr_array( - # (con_data, con_index, con_index_ptr), [len(rows), nCol] - # ) A = A.tocsc() A.sum_duplicates() A.eliminate_zeros() if with_debug_timing: timer.toc('Formed matrices', level=logging.DEBUG) - # >>>>>>> templatized-writer # Some variables in the var_map may not actually appear in the # objective or constraints (e.g., added from col_order, or From 470de80249a0df2bd7cfa52d64220f9fd5686439 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 8 Jul 2024 13:04:12 -0600 Subject: [PATCH 28/65] Adding a copy of scipy's sum duplicates code, but I think we can cheat a little if we assume ordered, so this isn't done or tested --- .../plugins/parameterized_standard_form.py | 20 +++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index 39e25954680..d8be759f2ad 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -189,7 +189,23 @@ def todense(self): return dense def sum_duplicates(self): - pass + nnz = 0 + ncols = self.shape[1] + col_end = 0 + for i in range(ncols): + jj = col_end + col_end = self.row_index[i + 1] + while jj < col_end: + j = self.index[jj] + x = self.data[jj] + jj += 1 + while jj < col_end and self.index[jj] == j: + x += self.data[jj] + jj += 1 + self.row_index[nnz] = j + self.data[nnz] = x + nnz += 1 + self.indptr[i + 1] = nnz def eliminate_zeros(self): - pass + raise NotImplementedError From 3edc9643b017cc92156b654268a9166e8feeb568 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 8 Jul 2024 15:35:27 -0600 Subject: [PATCH 29/65] Removing an unused method --- pyomo/repn/plugins/parameterized_standard_form.py | 3 --- pyomo/repn/plugins/standard_form.py | 4 ---- 2 files changed, 7 deletions(-) diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index d8be759f2ad..0efef7c402e 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -90,9 +90,6 @@ def _to_vector(self, data, N, vector_type): # on the Pyomo expressions return np.array([v for v in data]) - def _csc_matrix_from_csr(self, data, index, index_ptr, nrows, ncols): - return _CSRMatrix(data, index, index_ptr, nrows, ncols).tocsc() - def _csc_matrix(self, data, index, index_ptr, nrows, ncols): return _CSCMatrix(data, index, index_ptr, nrows, ncols) diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index 79d6f731baf..ff4403523cf 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -256,9 +256,6 @@ def _get_visitor(self, var_map, var_order, sorter): def _to_vector(self, data, N, vector_type): return np.fromiter(data, vector_type, N) - def _csc_matrix_from_csr(self, data, index, index_ptr, nrows, ncols): - return scipy.sparse.csr_array((data, index, index_ptr), [nrows, ncols]).tocsc() - def _csc_matrix(self, data, index, index_ptr, nrows, ncols): return scipy.sparse.csc_array((data, index, index_ptr), [nrows, ncols]) @@ -451,7 +448,6 @@ def write(self, model): if mixed_form: if lb == ub: con_nnz += N - # >>>>>>> templatized-writer rows.append(RowEntry(con, 0)) rhs.append(ub - offset) con_data.append(linear_data) From d86e00b6079eed3912149e57f1af2071a2c590fe Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Wed, 10 Jul 2024 21:44:33 -0600 Subject: [PATCH 30/65] Adjusting parameterized standard form to match the changes in standard form with attributes rather than methods, implementing eliminate_zeros --- .../plugins/parameterized_standard_form.py | 88 ++++++++++++------- .../tests/test_parameterized_standard_form.py | 10 +-- 2 files changed, 61 insertions(+), 37 deletions(-) diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index 0efef7c402e..94beb76c9f1 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -12,6 +12,7 @@ from pyomo.common.config import ConfigValue, document_kwargs_from_configdict from pyomo.common.dependencies import numpy as np from pyomo.common.gc_manager import PauseGC +from pyomo.common.numeric_types import native_numeric_types from pyomo.opt import WriterFactory from pyomo.repn.parameterized_linear import ParameterizedLinearRepnVisitor @@ -78,27 +79,11 @@ def write(self, model, ostream=None, **options): return _ParameterizedLinearStandardFormCompiler_impl(config).write(model) -class _ParameterizedLinearStandardFormCompiler_impl(_LinearStandardFormCompiler_impl): - def _get_visitor(self, var_map, var_order, sorter): - wrt = self.config.wrt - if wrt is None: - wrt = [] - return ParameterizedLinearRepnVisitor({}, var_map, var_order, sorter, wrt=wrt) - - def _to_vector(self, data, N, vector_type): - # override this to not attempt conversion to float since that will fail - # on the Pyomo expressions - return np.array([v for v in data]) - - def _csc_matrix(self, data, index, index_ptr, nrows, ncols): - return _CSCMatrix(data, index, index_ptr, nrows, ncols) - - def _csr_matrix(self, data, index, index_ptr, nrows, ncols): - return _CSRMatrix(data, index, index_ptr, nrows, ncols) - - class _SparseMatrixBase(object): - def __init__(self, data, indices, indptr, nrows, ncols): + def __init__(self, matrix_data, shape): + (data, indices, indptr) = matrix_data + (nrows, ncols) = shape + self.data = np.array(data) self.indices = np.array(indices, dtype=int) self.indptr = np.array(indptr, dtype=int) @@ -153,7 +138,7 @@ def tocsc(self): # above. col_index_ptr = np.insert(col_index_ptr, 0, 0) - return _CSCMatrix(csc_data, row_index, col_index_ptr, *self.shape) + return _CSCMatrix((csc_data, row_index, col_index_ptr), self.shape) def todense(self): nrows = self.shape[0] @@ -186,23 +171,62 @@ def todense(self): return dense def sum_duplicates(self): - nnz = 0 ncols = self.shape[1] + row_index = self.indices + col_index_ptr = self.indptr + data = self.data + + num_non_zeros = 0 col_end = 0 for i in range(ncols): jj = col_end - col_end = self.row_index[i + 1] + col_end = col_index_ptr[i + 1] while jj < col_end: - j = self.index[jj] - x = self.data[jj] + j = row_index[jj] + x = data[jj] jj += 1 - while jj < col_end and self.index[jj] == j: - x += self.data[jj] + while jj < col_end and row_index[jj] == j: + x += data[jj] jj += 1 - self.row_index[nnz] = j - self.data[nnz] = x - nnz += 1 - self.indptr[i + 1] = nnz + row_index[num_non_zeros] = j + data[num_non_zeros] = x + num_non_zeros += 1 + col_index_ptr[i + 1] = num_non_zeros def eliminate_zeros(self): - raise NotImplementedError + ncols = self.shape[1] + row_index = self.indices + col_index_ptr = self.indptr + data = self.data + + num_non_zeros = 0 + col_end = 0 + for i in range(ncols): + jj = col_end + col_end = col_index_ptr[i + 1] + while jj < col_end: + j = row_index[jj] + x = data[jj] + if x.__class__ not in native_numeric_types or x != 0: + row_index[num_non_zeros] = j + data[num_non_zeros] = x + num_non_zeros += 1 + jj += 1 + col_index_ptr[i + 1] = num_non_zeros + + +class _ParameterizedLinearStandardFormCompiler_impl(_LinearStandardFormCompiler_impl): + _csc_matrix = _CSCMatrix + _csr_matrix = _CSRMatrix + + def _get_visitor(self, subexpression_cache, var_map, var_order, sorter): + wrt = self.config.wrt + if wrt is None: + wrt = [] + return ParameterizedLinearRepnVisitor(subexpression_cache, var_map, + var_order, sorter, wrt=wrt) + + def _to_vector(self, data, N, vector_type): + # override this to not attempt conversion to float since that will fail + # on the Pyomo expressions + return np.array([v for v in data]) diff --git a/pyomo/repn/tests/test_parameterized_standard_form.py b/pyomo/repn/tests/test_parameterized_standard_form.py index 5ae2889ce83..6d549a14f65 100644 --- a/pyomo/repn/tests/test_parameterized_standard_form.py +++ b/pyomo/repn/tests/test_parameterized_standard_form.py @@ -44,7 +44,7 @@ ) class TestSparseMatrixRepresentations(unittest.TestCase): def test_csr_to_csc_only_data(self): - A = _CSRMatrix([5, 8, 3, 6], [0, 1, 2, 1], [0, 1, 2, 3, 4], 4, 4) + A = _CSRMatrix(([5, 8, 3, 6], [0, 1, 2, 1], [0, 1, 2, 3, 4]), [4, 4]) thing = A.tocsc() self.assertTrue(np.all(thing.data == np.array([5, 8, 6, 3]))) @@ -57,7 +57,7 @@ def test_csr_to_csc_pyomo_exprs(self): m.y = Var() A = _CSRMatrix( - [5, 8 * m.x, 3 * m.x * m.y**2, 6], [0, 1, 2, 1], [0, 1, 2, 3, 4], 4, 4 + ([5, 8 * m.x, 3 * m.x * m.y**2, 6], [0, 1, 2, 1], [0, 1, 2, 3, 4]), [4, 4] ) thing = A.tocsc() @@ -71,7 +71,7 @@ def test_csr_to_csc_pyomo_exprs(self): self.assertTrue(np.all(thing.indptr == np.array([0, 1, 3, 4, 4]))) def test_csr_to_csc_empty_matrix(self): - A = _CSRMatrix([], [], [0], 0, 4) + A = _CSRMatrix(([], [], [0]), [0, 4]) thing = A.tocsc() self.assertEqual(thing.data.size, 0) @@ -80,13 +80,13 @@ def test_csr_to_csc_empty_matrix(self): self.assertTrue(np.all(thing.indptr == np.zeros(5))) def test_todense(self): - A = _CSRMatrix([5, 8, 3, 6], [0, 1, 2, 1], [0, 1, 2, 3, 4], 4, 4) + A = _CSRMatrix(([5, 8, 3, 6], [0, 1, 2, 1], [0, 1, 2, 3, 4]), [4, 4]) dense = np.array([[5, 0, 0, 0], [0, 8, 0, 0], [0, 0, 3, 0], [0, 6, 0, 0]]) self.assertTrue(np.all(A.todense() == dense)) self.assertTrue(np.all(A.tocsc().todense() == dense)) - A = _CSRMatrix([5, 6, 7, 2, 1, 1.5], [0, 1, 1, 2, 3, 1], [0, 2, 4, 5, 6], 4, 4) + A = _CSRMatrix(([5, 6, 7, 2, 1, 1.5], [0, 1, 1, 2, 3, 1], [0, 2, 4, 5, 6]), [4, 4]) dense = np.array([[5, 6, 0, 0], [0, 7, 2, 0], [0, 0, 0, 1], [0, 1.5, 0, 0]]) self.assertTrue(np.all(A.todense() == dense)) self.assertTrue(np.all(A.tocsc().todense() == dense)) From c5698978b1ad54e44ce79eb2d4919e0a76afb233 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Wed, 10 Jul 2024 21:52:09 -0600 Subject: [PATCH 31/65] black --- pyomo/repn/linear_template.py | 1 - pyomo/repn/plugins/parameterized_standard_form.py | 7 ++++--- pyomo/repn/plugins/standard_form.py | 4 ++-- pyomo/repn/tests/test_parameterized_standard_form.py | 4 +++- 4 files changed, 9 insertions(+), 7 deletions(-) diff --git a/pyomo/repn/linear_template.py b/pyomo/repn/linear_template.py index 7a00a1d5cf3..995ce370e55 100644 --- a/pyomo/repn/linear_template.py +++ b/pyomo/repn/linear_template.py @@ -266,7 +266,6 @@ def _before_named_expression(self, visitor, child): raise NotImplementedError() - def _handle_getitem(visitor, node, comp, *args): expr = comp[1][tuple(arg[1] for arg in args)] if comp[0] is _CONSTANT: diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index 94beb76c9f1..fb266284741 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -213,7 +213,7 @@ def eliminate_zeros(self): num_non_zeros += 1 jj += 1 col_index_ptr[i + 1] = num_non_zeros - + class _ParameterizedLinearStandardFormCompiler_impl(_LinearStandardFormCompiler_impl): _csc_matrix = _CSCMatrix @@ -223,8 +223,9 @@ def _get_visitor(self, subexpression_cache, var_map, var_order, sorter): wrt = self.config.wrt if wrt is None: wrt = [] - return ParameterizedLinearRepnVisitor(subexpression_cache, var_map, - var_order, sorter, wrt=wrt) + return ParameterizedLinearRepnVisitor( + subexpression_cache, var_map, var_order, sorter, wrt=wrt + ) def _to_vector(self, data, N, vector_type): # override this to not attempt conversion to float since that will fail diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index ef2f8ba9332..cd73ca96695 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -349,7 +349,7 @@ def write(self, model): obj_index_ptr = [0] for obj in objectives: if hasattr(obj, 'template_expr'): - offset, linear_index, linear_data, _, _ = ( + (offset, linear_index, linear_data, _, _) = ( template_visitor.expand_expression(obj, obj.template_expr()) ) N = len(linear_index) @@ -398,7 +398,7 @@ def write(self, model): last_parent = con._component if hasattr(con, 'template_expr'): - offset, linear_index, linear_data, lb, ub = ( + (offset, linear_index, linear_data, lb, ub) = ( template_visitor.expand_expression(con, con.template_expr()) ) N = len(linear_data) diff --git a/pyomo/repn/tests/test_parameterized_standard_form.py b/pyomo/repn/tests/test_parameterized_standard_form.py index 6d549a14f65..56ed16342b2 100644 --- a/pyomo/repn/tests/test_parameterized_standard_form.py +++ b/pyomo/repn/tests/test_parameterized_standard_form.py @@ -86,7 +86,9 @@ def test_todense(self): self.assertTrue(np.all(A.todense() == dense)) self.assertTrue(np.all(A.tocsc().todense() == dense)) - A = _CSRMatrix(([5, 6, 7, 2, 1, 1.5], [0, 1, 1, 2, 3, 1], [0, 2, 4, 5, 6]), [4, 4]) + A = _CSRMatrix( + ([5, 6, 7, 2, 1, 1.5], [0, 1, 1, 2, 3, 1], [0, 2, 4, 5, 6]), [4, 4] + ) dense = np.array([[5, 6, 0, 0], [0, 7, 2, 0], [0, 0, 0, 1], [0, 1.5, 0, 0]]) self.assertTrue(np.all(A.todense() == dense)) self.assertTrue(np.all(A.tocsc().todense() == dense)) From 70d1fbe7da6bfd43abad50c617188ee8f0ad2ea4 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Thu, 11 Jul 2024 09:30:12 -0600 Subject: [PATCH 32/65] Converting the c matrix dense for now to get things working --- pyomo/core/plugins/transform/lp_dual.py | 20 ++++++-------------- pyomo/core/tests/unit/test_lp_dual.py | 19 +++++++++++-------- 2 files changed, 17 insertions(+), 22 deletions(-) diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index 6ce2292a663..7c49b762ad5 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -117,6 +117,7 @@ def _take_dual(self, model, std_form): # This is a csc matrix, so we'll skipping transposing and just work off # of the folumns A = std_form.A + c = std_form.c.todense() dual_rows = range(A.shape[1]) dual_cols = range(A.shape[0]) dual.x = Var(dual_cols, domain=NonNegativeReals) @@ -141,27 +142,18 @@ def _take_dual(self, model, std_form): primal_row = A.indices[j] lhs += coef * dual.x[primal_row] - # coef = A[j, i] - # if not coef: - # continue - # elif isclose_const(coef, 1.0): - # lhs += dual.x[j] - # elif isclose_const(coef, -1.0): - # lhs -= dual.x[j] - # else: - # lhs += float(A[j, i]) * dual.x[j] if primal.domain is Reals: - dual.constraints[i] = lhs == std_form.c[0, i] + dual.constraints[i] = lhs == c[0, i] elif primal_sense is minimize: if primal.domain is NonNegativeReals: - dual.constraints[i] = lhs <= std_form.c[0, i] + dual.constraints[i] = lhs <= c[0, i] else: # primal.domain is NonPositiveReals - dual.constraints[i] = lhs >= std_form.c[0, i] + dual.constraints[i] = lhs >= c[0, i] else: if primal.domain is NonNegativeReals: - dual.constraints[i] = lhs >= std_form.c[0, i] + dual.constraints[i] = lhs >= c[0, i] else: # primal.domain is NonPositiveReals - dual.constraints[i] = lhs <= std_form.c[0, i] + dual.constraints[i] = lhs <= c[0, i] trans_info.dual_constraint[primal] = dual.constraints[i] trans_info.primal_var[dual.constraints[i]] = primal diff --git a/pyomo/core/tests/unit/test_lp_dual.py b/pyomo/core/tests/unit/test_lp_dual.py index 733d11275e2..76b17a8677f 100644 --- a/pyomo/core/tests/unit/test_lp_dual.py +++ b/pyomo/core/tests/unit/test_lp_dual.py @@ -27,7 +27,10 @@ value, Var, ) -from pyomo.core.expr.compare import assertExpressionsEqual +from pyomo.core.expr.compare import ( + assertExpressionsEqual, + assertExpressionsStructurallyEqual, +) from pyomo.opt import SolverFactory, WriterFactory @@ -130,9 +133,9 @@ def test_lp_dual(self): self.assertIs(dy.ctype, Constraint) self.assertIs(dz.ctype, Constraint) - assertExpressionsEqual(self, dx.expr, -4.0 * alpha + beta <= 1.0) - assertExpressionsEqual(self, dy.expr, -2.0 * alpha + beta - lamb >= 2.0) - assertExpressionsEqual(self, dz.expr, -alpha - 1.0 * lamb + xi == -3.0) + assertExpressionsStructurallyEqual(self, dx.expr, -4.0 * alpha + beta <= 1.0) + assertExpressionsStructurallyEqual(self, dy.expr, -2.0 * alpha + beta - lamb >= 2.0) + assertExpressionsStructurallyEqual(self, dz.expr, -alpha - 1.0 * lamb + xi == -3.0) dual_obj = dual.obj self.assertIsInstance(dual_obj, Objective) @@ -183,10 +186,10 @@ def test_lp_dual(self): self.assertIs(dlambda.ctype, Constraint) self.assertIs(dxi.ctype, Constraint) - assertExpressionsEqual(self, dalpha.expr, -4.0 * x - 2.0 * y - z <= -5.0) - assertExpressionsEqual(self, dbeta.expr, x + y >= 3.0) - assertExpressionsEqual(self, dlambda.expr, -y - z == -4.2) - assertExpressionsEqual(self, dxi.expr, z <= 42.0) + assertExpressionsStructurallyEqual(self, dalpha.expr, -4.0 * x - 2.0 * y - z <= -5.0) + assertExpressionsStructurallyEqual(self, dbeta.expr, x + y >= 3.0) + assertExpressionsStructurallyEqual(self, dlambda.expr, -y - z == -4.2) + assertExpressionsStructurallyEqual(self, dxi.expr, z <= 42.0) primal_obj = primal.obj self.assertIsInstance(primal_obj, Objective) From d7d3b7470faf3a78a02b774772567f3dae549252 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Thu, 11 Jul 2024 09:42:55 -0600 Subject: [PATCH 33/65] Adding a test of a parameterized dual --- pyomo/core/tests/unit/test_lp_dual.py | 71 ++++++++++++++++++++++++--- 1 file changed, 65 insertions(+), 6 deletions(-) diff --git a/pyomo/core/tests/unit/test_lp_dual.py b/pyomo/core/tests/unit/test_lp_dual.py index 76b17a8677f..b46d3902c8c 100644 --- a/pyomo/core/tests/unit/test_lp_dual.py +++ b/pyomo/core/tests/unit/test_lp_dual.py @@ -41,7 +41,6 @@ class TestLPDual(unittest.TestCase): and SolverFactory('gurobi').license_is_valid(), "Gurobi is not available", ) - def test_lp_dual_solve(self): m = ConcreteModel() m.x = Var(domain=NonNegativeReals) @@ -134,8 +133,12 @@ def test_lp_dual(self): self.assertIs(dz.ctype, Constraint) assertExpressionsStructurallyEqual(self, dx.expr, -4.0 * alpha + beta <= 1.0) - assertExpressionsStructurallyEqual(self, dy.expr, -2.0 * alpha + beta - lamb >= 2.0) - assertExpressionsStructurallyEqual(self, dz.expr, -alpha - 1.0 * lamb + xi == -3.0) + assertExpressionsStructurallyEqual( + self, dy.expr, -2.0 * alpha + beta - lamb >= 2.0 + ) + assertExpressionsStructurallyEqual( + self, dz.expr, -alpha - 1.0 * lamb + xi == -3.0 + ) dual_obj = dual.obj self.assertIsInstance(dual_obj, Objective) @@ -186,7 +189,9 @@ def test_lp_dual(self): self.assertIs(dlambda.ctype, Constraint) self.assertIs(dxi.ctype, Constraint) - assertExpressionsStructurallyEqual(self, dalpha.expr, -4.0 * x - 2.0 * y - z <= -5.0) + assertExpressionsStructurallyEqual( + self, dalpha.expr, -4.0 * x - 2.0 * y - z <= -5.0 + ) assertExpressionsStructurallyEqual(self, dbeta.expr, x + y >= 3.0) assertExpressionsStructurallyEqual(self, dlambda.expr, -y - z == -4.2) assertExpressionsStructurallyEqual(self, dxi.expr, z <= 42.0) @@ -207,10 +212,64 @@ def test_parameterized_linear_dual(self): m.z = Var(domain=Reals) m.obj = Objective(expr=m.x + 2 * m.y - 3 * m.outer[3] * m.z) - m.c1 = Constraint(expr=-4 * m.x - 2 * m.y - m.z <= -5*m.outer1) - m.c2 = Constraint(expr=m.x + m.outer[2]*m.y >= 3) + m.c1 = Constraint(expr=-4 * m.x - 2 * m.y - m.z <= -5 * m.outer1) + m.c2 = Constraint(expr=m.x + m.outer[2] * m.y >= 3) m.c3 = Constraint(expr=-m.y - m.z == -4.2) m.c4 = Constraint(expr=m.z <= 42) lp_dual = TransformationFactory('core.lp_dual') dual = lp_dual.create_using(m, parameterize_wrt=[m.outer1, m.outer]) + + alpha = lp_dual.get_dual_var(dual, m.c1) + beta = lp_dual.get_dual_var(dual, m.c2) + lamb = lp_dual.get_dual_var(dual, m.c3) + mu = lp_dual.get_dual_var(dual, m.c4) + + self.assertIs(lp_dual.get_primal_constraint(dual, alpha), m.c1) + self.assertIs(lp_dual.get_primal_constraint(dual, beta), m.c2) + self.assertIs(lp_dual.get_primal_constraint(dual, lamb), m.c3) + self.assertIs(lp_dual.get_primal_constraint(dual, mu), m.c4) + + dx = lp_dual.get_dual_constraint(dual, m.x) + dy = lp_dual.get_dual_constraint(dual, m.y) + dz = lp_dual.get_dual_constraint(dual, m.z) + + self.assertIs(lp_dual.get_primal_var(dual, dx), m.x) + self.assertIs(lp_dual.get_primal_var(dual, dy), m.y) + self.assertIs(lp_dual.get_primal_var(dual, dz), m.z) + + self.assertIs(alpha.ctype, Var) + self.assertEqual(alpha.domain, NonPositiveReals) + self.assertEqual(alpha.ub, 0) + self.assertIsNone(alpha.lb) + self.assertIs(beta.ctype, Var) + self.assertEqual(beta.domain, NonNegativeReals) + self.assertEqual(beta.lb, 0) + self.assertIsNone(beta.ub) + self.assertIs(lamb.ctype, Var) + self.assertEqual(lamb.domain, Reals) + self.assertIsNone(lamb.ub) + self.assertIsNone(lamb.lb) + self.assertIs(mu.ctype, Var) + self.assertEqual(mu.domain, NonPositiveReals) + self.assertEqual(mu.ub, 0) + self.assertIsNone(mu.lb) + + self.assertIs(dx.ctype, Constraint) + self.assertIs(dy.ctype, Constraint) + self.assertIs(dz.ctype, Constraint) + + assertExpressionsStructurallyEqual(self, dx.expr, -4.0 * alpha + beta <= 1.0) + assertExpressionsStructurallyEqual( + self, dy.expr, -2.0 * alpha + m.outer[2] * beta - lamb >= 2.0 + ) + assertExpressionsStructurallyEqual( + self, dz.expr, -alpha - 1.0 * lamb + mu == -3.0 * m.outer[3] + ) + + dual_obj = dual.obj + self.assertIsInstance(dual_obj, Objective) + self.assertEqual(dual_obj.sense, maximize) + assertExpressionsEqual( + self, dual_obj.expr, -5 * m.outer1 * alpha + 3 * beta - 4.2 * lamb + 42 * mu + ) From 82326af3acbe4fc9a92e374a9b5f91ac085cfecd Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Thu, 11 Jul 2024 09:46:33 -0600 Subject: [PATCH 34/65] black --- pyomo/core/plugins/transform/lp_dual.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index 7c49b762ad5..8777415ad71 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -137,7 +137,7 @@ def _take_dual(self, model, std_form): dual.constraints = Constraint(dual_rows) for i, primal in enumerate(std_form.columns): lhs = 0 - for j in range(A.indptr[i], A.indptr[i+1]): + for j in range(A.indptr[i], A.indptr[i + 1]): coef = A.data[j] primal_row = A.indices[j] lhs += coef * dual.x[primal_row] @@ -152,13 +152,15 @@ def _take_dual(self, model, std_form): else: if primal.domain is NonNegativeReals: dual.constraints[i] = lhs >= c[0, i] - else: # primal.domain is NonPositiveReals + else: # primal.domain is NonPositiveReals dual.constraints[i] = lhs <= c[0, i] trans_info.dual_constraint[primal] = dual.constraints[i] trans_info.primal_var[dual.constraints[i]] = primal - dual.obj = Objective( expr=sum(std_form.rhs[j] * dual.x[j] for j in - dual_cols), sense=-primal_sense ) + dual.obj = Objective( + expr=sum(std_form.rhs[j] * dual.x[j] for j in dual_cols), + sense=-primal_sense, + ) return dual From 6f5790a50fd23c8cbe893d4480d0d3f989249f2f Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 4 Nov 2024 15:12:08 -0700 Subject: [PATCH 35/65] Converting to var_recorder argument for parameterized standard form --- pyomo/repn/linear_template.py | 1 + pyomo/repn/plugins/parameterized_standard_form.py | 7 ++++--- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/pyomo/repn/linear_template.py b/pyomo/repn/linear_template.py index e6f70e7b59f..55cca7d6c81 100644 --- a/pyomo/repn/linear_template.py +++ b/pyomo/repn/linear_template.py @@ -241,6 +241,7 @@ def _before_component(self, visitor, child): def _before_named_expression(self, visitor, child): raise MouseTrap("We do not yet support Expression components") + def _handle_getitem(visitor, node, comp, *args): expr = comp[1][tuple(arg[1] for arg in args)] if comp[0] is _CONSTANT: diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index fb266284741..36bbe038a00 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -76,7 +76,8 @@ def write(self, model, ostream=None, **options): # representation generates (and disposes of) a large number of # small objects. with PauseGC(): - return _ParameterizedLinearStandardFormCompiler_impl(config).write(model) + return _ParameterizedLinearStandardFormCompiler_impl(config).write( + model) class _SparseMatrixBase(object): @@ -219,12 +220,12 @@ class _ParameterizedLinearStandardFormCompiler_impl(_LinearStandardFormCompiler_ _csc_matrix = _CSCMatrix _csr_matrix = _CSRMatrix - def _get_visitor(self, subexpression_cache, var_map, var_order, sorter): + def _get_visitor(self, subexpression_cache, var_recorder): wrt = self.config.wrt if wrt is None: wrt = [] return ParameterizedLinearRepnVisitor( - subexpression_cache, var_map, var_order, sorter, wrt=wrt + subexpression_cache, wrt=wrt, var_recorder=var_recorder ) def _to_vector(self, data, N, vector_type): From 139d307797427c9c859f78c0101489b2b457a497 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 4 Nov 2024 15:18:29 -0700 Subject: [PATCH 36/65] black --- pyomo/repn/plugins/parameterized_standard_form.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index 36bbe038a00..4367f80e739 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -76,8 +76,7 @@ def write(self, model, ostream=None, **options): # representation generates (and disposes of) a large number of # small objects. with PauseGC(): - return _ParameterizedLinearStandardFormCompiler_impl(config).write( - model) + return _ParameterizedLinearStandardFormCompiler_impl(config).write(model) class _SparseMatrixBase(object): From 2f40eb93f942c614c9f3ac4f47a5e3461424883f Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 4 Nov 2024 15:20:42 -0700 Subject: [PATCH 37/65] More black --- pyomo/util/var_list_domain.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pyomo/util/var_list_domain.py b/pyomo/util/var_list_domain.py index f890b308d76..15c0f775b57 100644 --- a/pyomo/util/var_list_domain.py +++ b/pyomo/util/var_list_domain.py @@ -12,6 +12,7 @@ from pyomo.common.collections import ComponentSet from pyomo.core import Var + def var_component_set(x): """ For domain validation in ConfigDicts: Takes singletone or iterable argument 'x' From 5a25e2a4f34436201d8d7dbf13363b65196a18d0 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 4 Nov 2024 15:32:39 -0700 Subject: [PATCH 38/65] NFC: fixing comment typo --- pyomo/core/plugins/transform/lp_dual.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index 8777415ad71..bd94cc96852 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -114,8 +114,8 @@ def _take_dual(self, model, std_form): primal_sense = std_form.objectives[0].sense dual = ConcreteModel(name="%s dual" % model.name) - # This is a csc matrix, so we'll skipping transposing and just work off - # of the folumns + # This is a csc matrix, so we'll skip transposing and just work off + # of the columns A = std_form.A c = std_form.c.todense() dual_rows = range(A.shape[1]) From c2c5c0be2798d745014b033153bf3667fc17919d Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 4 Nov 2024 16:07:52 -0700 Subject: [PATCH 39/65] Adding docstrings for mapping API and clearning up some prepositional mix-ups --- pyomo/core/plugins/transform/lp_dual.py | 60 ++++++++++++++++++++++++- 1 file changed, 58 insertions(+), 2 deletions(-) diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index bd94cc96852..a6c41d41053 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -165,6 +165,20 @@ def _take_dual(self, model, std_form): return dual def get_primal_constraint(self, model, dual_var): + """Return the primal constraint corresponding to 'dual_var' + + Returns + ------- + Constraint + + Parameters + ---------- + model: ConcreteModel + A dual model returned from the 'core.lp_dual' transformation + dual_var: Var + A dual variable on 'model' + + """ primal_constraint = model.private_data().primal_constraint if dual_var in primal_constraint: return primal_constraint[dual_var] @@ -175,16 +189,44 @@ def get_primal_constraint(self, model, dual_var): ) def get_dual_constraint(self, model, primal_var): + """Return the dual constraint corresponding to 'primal_var' + + Returns + ------- + Constraint + + Parameters + ---------- + model: ConcreteModel + A primal model passed as an argument to the 'core.lp_dual' transformation + primal_var: Var + A primal variable on 'model' + + """ dual_constraint = model.private_data().dual_constraint if primal_var in dual_constraint: return dual_constraint[primal_var] else: raise ValueError( - "It does not appear that Var '%s' is a primal variable from model '%s'" + "It does not appear that Var '%s' is a primal variable on model '%s'" % (primal_var.name, model.name) ) def get_primal_var(self, model, dual_constraint): + """Return the primal variable corresponding to 'dual_constraint' + + Returns + ------- + Var + + Parameters + ---------- + model: ConcreteModel + A dual model returned from the 'core.lp_dual' transformation + dual_constraint: Constraint + A constraint on 'model' + + """ primal_var = model.private_data().primal_var if dual_constraint in primal_var: return primal_var[dual_constraint] @@ -195,11 +237,25 @@ def get_primal_var(self, model, dual_constraint): ) def get_dual_var(self, model, primal_constraint): + """Return the dual variable corresponding to 'primal_constraint' + + Returns + ------- + Var + + Parameters + ---------- + model: ConcreteModel + A primal model passed as an argument to the 'core.lp_dual' transformation + primal_constraint: Constraint + A constraint on 'model' + + """ dual_var = model.private_data().dual_var if primal_constraint in dual_var: return dual_var[primal_constraint] else: raise ValueError( - "It does not appear that Constraint '%s' is a primal constraint from " + "It does not appear that Constraint '%s' is a primal constraint on " "model '%s'" % (primal_constraint.name, model.name) ) From 7238c310055b29102c168940223eac037a0ce39f Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 4 Nov 2024 16:08:22 -0700 Subject: [PATCH 40/65] Adding tests for paraemeterized dual from solving primal and dual, adding tests for API error messages --- pyomo/core/tests/unit/test_lp_dual.py | 165 ++++++++++++++++++++++---- 1 file changed, 145 insertions(+), 20 deletions(-) diff --git a/pyomo/core/tests/unit/test_lp_dual.py b/pyomo/core/tests/unit/test_lp_dual.py index b46d3902c8c..34ebf5c741b 100644 --- a/pyomo/core/tests/unit/test_lp_dual.py +++ b/pyomo/core/tests/unit/test_lp_dual.py @@ -36,26 +36,10 @@ @unittest.skipUnless(scipy_available, "Scipy not available") class TestLPDual(unittest.TestCase): - @unittest.skipUnless( - SolverFactory('gurobi').available(exception_flag=False) - and SolverFactory('gurobi').license_is_valid(), - "Gurobi is not available", - ) - def test_lp_dual_solve(self): - m = ConcreteModel() - m.x = Var(domain=NonNegativeReals) - m.y = Var(domain=NonPositiveReals) - m.z = Var(domain=Reals) + def check_primal_dual_solns(self, m, dual): + lp_dual = TransformationFactory('core.lp_dual') - m.obj = Objective(expr=m.x + 2 * m.y - 3 * m.z) - m.c1 = Constraint(expr=-4 * m.x - 2 * m.y - m.z <= -5) - m.c2 = Constraint(expr=m.x + m.y <= 3) - m.c3 = Constraint(expr=-m.y - m.z <= -4.2) - m.c4 = Constraint(expr=m.z <= 42) m.dual = Suffix(direction=Suffix.IMPORT) - - lp_dual = TransformationFactory('core.lp_dual') - dual = lp_dual.create_using(m) dual.dual = Suffix(direction=Suffix.IMPORT) opt = SolverFactory('gurobi') @@ -78,6 +62,28 @@ def test_lp_dual_solve(self): value(v), value(dual.dual[lp_dual.get_dual_constraint(dual, v)]) ) + @unittest.skipUnless( + SolverFactory('gurobi').available(exception_flag=False) + and SolverFactory('gurobi').license_is_valid(), + "Gurobi is not available", + ) + def test_lp_dual_solve(self): + m = ConcreteModel() + m.x = Var(domain=NonNegativeReals) + m.y = Var(domain=NonPositiveReals) + m.z = Var(domain=Reals) + + m.obj = Objective(expr=m.x + 2 * m.y - 3 * m.z) + m.c1 = Constraint(expr=-4 * m.x - 2 * m.y - m.z <= -5) + m.c2 = Constraint(expr=m.x + m.y <= 3) + m.c3 = Constraint(expr=-m.y - m.z <= -4.2) + m.c4 = Constraint(expr=m.z <= 42) + + lp_dual = TransformationFactory('core.lp_dual') + dual = lp_dual.create_using(m) + + self.check_primal_dual_solns(m, dual) + def test_lp_dual(self): m = ConcreteModel() m.x = Var(domain=NonNegativeReals) @@ -201,8 +207,8 @@ def test_lp_dual(self): self.assertEqual(primal_obj.sense, minimize) assertExpressionsEqual(self, primal_obj.expr, x + 2.0 * y - 3.0 * z) - def test_parameterized_linear_dual(self): - m = ConcreteModel() + def get_bilevel_model(self): + m = ConcreteModel(name='primal') m.outer1 = Var(domain=Binary) m.outer = Var([2, 3], domain=Binary) @@ -217,6 +223,11 @@ def test_parameterized_linear_dual(self): m.c3 = Constraint(expr=-m.y - m.z == -4.2) m.c4 = Constraint(expr=m.z <= 42) + return m + + def test_parameterized_linear_dual(self): + m = self.get_bilevel_model() + lp_dual = TransformationFactory('core.lp_dual') dual = lp_dual.create_using(m, parameterize_wrt=[m.outer1, m.outer]) @@ -273,3 +284,117 @@ def test_parameterized_linear_dual(self): assertExpressionsEqual( self, dual_obj.expr, -5 * m.outer1 * alpha + 3 * beta - 4.2 * lamb + 42 * mu ) + + @unittest.skipUnless( + SolverFactory('gurobi').available(exception_flag=False) + and SolverFactory('gurobi').license_is_valid(), + "Gurobi is not available", + ) + def test_solve_parameterized_lp_dual(self): + m = self.get_bilevel_model() + + lp_dual = TransformationFactory('core.lp_dual') + dual = lp_dual.create_using(m, parameterize_wrt=[m.outer1, m.outer]) + + # We just check half of the possible permutations since we're calling a + # solver twice for all of these: + m.outer1.fix(1) + m.outer[2].fix(1) + m.outer[3].fix(1) + + self.check_primal_dual_solns(m, dual) + + m.outer1.fix(0) + m.outer[2].fix(1) + m.outer[3].fix(0) + + self.check_primal_dual_solns(m, dual) + + m.outer1.fix(0) + m.outer[2].fix(0) + m.outer[3].fix(0) + + self.check_primal_dual_solns(m, dual) + + m.outer1.fix(1) + m.outer[2].fix(1) + m.outer[3].fix(0) + + self.check_primal_dual_solns(m, dual) + + def test_multiple_obj_error(self): + m = self.get_bilevel_model() + m.obj.deactivate() + + lp_dual = TransformationFactory('core.lp_dual') + + with self.assertRaisesRegex( + ValueError, + "Model 'primal' has no objective or multiple active objectives. " + "Cannot take dual with more than one objective!" + ): + dual = lp_dual.create_using(m, parameterize_wrt=[m.outer1, m.outer]) + + m.obj.activate() + m.obj2 = Objective(expr=m.outer1 + m.outer[3]) + + with self.assertRaisesRegex( + ValueError, + "Model 'primal' has no objective or multiple active objectives. " + "Cannot take dual with more than one objective!" + ): + dual = lp_dual.create_using(m, parameterize_wrt=[m.outer1, m.outer]) + + def test_primal_constraint_map_error(self): + m = self.get_bilevel_model() + + lp_dual = TransformationFactory('core.lp_dual') + dual = lp_dual.create_using(m, parameterize_wrt=[m.outer1, m.outer]) + + with self.assertRaisesRegex( + ValueError, + "It does not appear that Var 'x' is a dual variable on model " + "'primal dual'" + ): + thing = lp_dual.get_primal_constraint(dual, m.x) + + def test_dual_constraint_map_error(self): + m = self.get_bilevel_model() + + lp_dual = TransformationFactory('core.lp_dual') + dual = lp_dual.create_using(m, parameterize_wrt=[m.outer1, m.outer]) + + with self.assertRaisesRegex( + ValueError, + "It does not appear that Var 'outer1' is a primal variable on model " + "'primal'" + ): + thing = lp_dual.get_dual_constraint(m, m.outer1) + + def test_primal_var_map_error(self): + m = self.get_bilevel_model() + + lp_dual = TransformationFactory('core.lp_dual') + dual = lp_dual.create_using(m, parameterize_wrt=[m.outer1, m.outer]) + + with self.assertRaisesRegex( + ValueError, + "It does not appear that Constraint 'c1' is a dual constraint " + "on model 'primal dual'" + ): + thing = lp_dual.get_primal_var(dual, m.c1) + + def test_dual_var_map_error(self): + m = self.get_bilevel_model() + + lp_dual = TransformationFactory('core.lp_dual') + dual = lp_dual.create_using(m, parameterize_wrt=[m.outer1, m.outer]) + + m.c_new = Constraint(expr=m.x + m.y <= 35) + + with self.assertRaisesRegex( + ValueError, + "It does not appear that Constraint 'c_new' is a primal constraint " + "on model 'primal'" + ): + thing = lp_dual.get_dual_var(m, m.c_new) From d914e52992d796f521a56adbc5157905f007a3c2 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 4 Nov 2024 16:09:21 -0700 Subject: [PATCH 41/65] Black --- pyomo/core/plugins/transform/lp_dual.py | 8 ++--- pyomo/core/tests/unit/test_lp_dual.py | 48 ++++++++++++------------- 2 files changed, 28 insertions(+), 28 deletions(-) diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index a6c41d41053..1a2d20efd69 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -177,7 +177,7 @@ def get_primal_constraint(self, model, dual_var): A dual model returned from the 'core.lp_dual' transformation dual_var: Var A dual variable on 'model' - + """ primal_constraint = model.private_data().primal_constraint if dual_var in primal_constraint: @@ -201,7 +201,7 @@ def get_dual_constraint(self, model, primal_var): A primal model passed as an argument to the 'core.lp_dual' transformation primal_var: Var A primal variable on 'model' - + """ dual_constraint = model.private_data().dual_constraint if primal_var in dual_constraint: @@ -225,7 +225,7 @@ def get_primal_var(self, model, dual_constraint): A dual model returned from the 'core.lp_dual' transformation dual_constraint: Constraint A constraint on 'model' - + """ primal_var = model.private_data().primal_var if dual_constraint in primal_var: @@ -249,7 +249,7 @@ def get_dual_var(self, model, primal_constraint): A primal model passed as an argument to the 'core.lp_dual' transformation primal_constraint: Constraint A constraint on 'model' - + """ dual_var = model.private_data().dual_var if primal_constraint in dual_var: diff --git a/pyomo/core/tests/unit/test_lp_dual.py b/pyomo/core/tests/unit/test_lp_dual.py index 34ebf5c741b..33e3e9287cd 100644 --- a/pyomo/core/tests/unit/test_lp_dual.py +++ b/pyomo/core/tests/unit/test_lp_dual.py @@ -327,11 +327,11 @@ def test_multiple_obj_error(self): m.obj.deactivate() lp_dual = TransformationFactory('core.lp_dual') - + with self.assertRaisesRegex( - ValueError, - "Model 'primal' has no objective or multiple active objectives. " - "Cannot take dual with more than one objective!" + ValueError, + "Model 'primal' has no objective or multiple active objectives. " + "Cannot take dual with more than one objective!", ): dual = lp_dual.create_using(m, parameterize_wrt=[m.outer1, m.outer]) @@ -339,9 +339,9 @@ def test_multiple_obj_error(self): m.obj2 = Objective(expr=m.outer1 + m.outer[3]) with self.assertRaisesRegex( - ValueError, - "Model 'primal' has no objective or multiple active objectives. " - "Cannot take dual with more than one objective!" + ValueError, + "Model 'primal' has no objective or multiple active objectives. " + "Cannot take dual with more than one objective!", ): dual = lp_dual.create_using(m, parameterize_wrt=[m.outer1, m.outer]) @@ -350,11 +350,11 @@ def test_primal_constraint_map_error(self): lp_dual = TransformationFactory('core.lp_dual') dual = lp_dual.create_using(m, parameterize_wrt=[m.outer1, m.outer]) - + with self.assertRaisesRegex( - ValueError, - "It does not appear that Var 'x' is a dual variable on model " - "'primal dual'" + ValueError, + "It does not appear that Var 'x' is a dual variable on model " + "'primal dual'", ): thing = lp_dual.get_primal_constraint(dual, m.x) @@ -363,24 +363,24 @@ def test_dual_constraint_map_error(self): lp_dual = TransformationFactory('core.lp_dual') dual = lp_dual.create_using(m, parameterize_wrt=[m.outer1, m.outer]) - + with self.assertRaisesRegex( - ValueError, - "It does not appear that Var 'outer1' is a primal variable on model " - "'primal'" + ValueError, + "It does not appear that Var 'outer1' is a primal variable on model " + "'primal'", ): thing = lp_dual.get_dual_constraint(m, m.outer1) - + def test_primal_var_map_error(self): m = self.get_bilevel_model() lp_dual = TransformationFactory('core.lp_dual') dual = lp_dual.create_using(m, parameterize_wrt=[m.outer1, m.outer]) - + with self.assertRaisesRegex( - ValueError, - "It does not appear that Constraint 'c1' is a dual constraint " - "on model 'primal dual'" + ValueError, + "It does not appear that Constraint 'c1' is a dual constraint " + "on model 'primal dual'", ): thing = lp_dual.get_primal_var(dual, m.c1) @@ -391,10 +391,10 @@ def test_dual_var_map_error(self): dual = lp_dual.create_using(m, parameterize_wrt=[m.outer1, m.outer]) m.c_new = Constraint(expr=m.x + m.y <= 35) - + with self.assertRaisesRegex( - ValueError, - "It does not appear that Constraint 'c_new' is a primal constraint " - "on model 'primal'" + ValueError, + "It does not appear that Constraint 'c_new' is a primal constraint " + "on model 'primal'", ): thing = lp_dual.get_dual_var(m, m.c_new) From d92d4df0a3bb89f58d9430e01b27fe0b54978e72 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 4 Nov 2024 16:13:10 -0700 Subject: [PATCH 42/65] NFC: fixing two typos --- pyomo/repn/tests/test_parameterized_standard_form.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyomo/repn/tests/test_parameterized_standard_form.py b/pyomo/repn/tests/test_parameterized_standard_form.py index 56ed16342b2..0b92071bc26 100644 --- a/pyomo/repn/tests/test_parameterized_standard_form.py +++ b/pyomo/repn/tests/test_parameterized_standard_form.py @@ -181,7 +181,7 @@ def test_parameterized_almost_dense_linear_model(self): ) self.assertTrue(np.all(repn.c == np.array([0, 0, 0]))) - # m.c gets interpretted as a <= Constraint, and you can't really blame + # m.c gets interpreted as a <= Constraint, and you can't really blame # pyomo for that because it's not parameterized yet. So that's why this # differs from the test in test_standard_form.py assertExpressionArraysEqual( @@ -242,7 +242,7 @@ def make_model(self, do_not_flip_c=False): m.more_data = Var() if do_not_flip_c: # [ESJ: 06/24]: I should have done this sooner, but if you write c - # this way, it gets interpretted as a >= constraint, which matches + # this way, it gets interpreted as a >= constraint, which matches # the standard_form tests and makes life much easier. Unforuntately # I wrote a lot of tests before I thought of this, so I'm leaving in # both variations for the moment. From 662fbfb9992a26750f6d7746f0db3e6d9410aee4 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Tue, 5 Nov 2024 13:05:53 -0700 Subject: [PATCH 43/65] Reverting unnecessary formatting changes in standard form --- pyomo/repn/plugins/standard_form.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index ff336bbb5c0..ac85d9a0656 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -378,7 +378,7 @@ def write(self, model): obj_index_ptr = [0] for obj in objectives: if hasattr(obj, 'template_expr'): - (offset, linear_index, linear_data, _, _) = ( + offset, linear_index, linear_data, _, _ = ( template_visitor.expand_expression(obj, obj.template_expr()) ) N = len(linear_index) @@ -427,7 +427,7 @@ def write(self, model): last_parent = con._component if hasattr(con, 'template_expr'): - (offset, linear_index, linear_data, lb, ub) = ( + offset, linear_index, linear_data, lb, ub = ( template_visitor.expand_expression(con, con.template_expr()) ) N = len(linear_data) From e4f31938087bae036c9a0e2f1aa8dfe5fe5da752 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Wed, 6 Nov 2024 08:42:58 -0700 Subject: [PATCH 44/65] Implementing Miranda's suggestion for the API mapping generalization --- pyomo/core/plugins/transform/lp_dual.py | 72 +++++++++++++++---------- 1 file changed, 44 insertions(+), 28 deletions(-) diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index 1a2d20efd69..8689faa4158 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -164,6 +164,41 @@ def _take_dual(self, model, std_form): return dual + def _get_corresponding_component(self, model, model_type, component, + component_type, mapping): + """Return the corresponding component based on the provided mapping. + + Parameters + ---------- + model: ConcreteModel + The model containing the mappings. + component: Var or Constraint + The primal/dual component for which we want to find the corresponding + dual/primal component. + component_type: str + A string indicating whether the component is 'dual' or 'primal'. + mapping: dict + The mapping to look up the corresponding component. + + Returns + ------- + Var or Constraint + The corresponding component. + + Raises + ------ + ValueError + If the component is not found in the mapping. + """ + if component in mapping: + return mapping[component] + else: + raise ValueError( + "It does not appear that %s '%s' is a %s %s on model '%s'" + % (component_type, component.name, model_type, 'variable' if component_type == 'Var' else 'constraint', + model.name) + ) + def get_primal_constraint(self, model, dual_var): """Return the primal constraint corresponding to 'dual_var' @@ -180,13 +215,8 @@ def get_primal_constraint(self, model, dual_var): """ primal_constraint = model.private_data().primal_constraint - if dual_var in primal_constraint: - return primal_constraint[dual_var] - else: - raise ValueError( - "It does not appear that Var '%s' is a dual variable on model '%s'" - % (dual_var.name, model.name) - ) + return self._get_corresponding_component(model, 'dual', dual_var, 'Var', + primal_constraint) def get_dual_constraint(self, model, primal_var): """Return the dual constraint corresponding to 'primal_var' @@ -204,13 +234,8 @@ def get_dual_constraint(self, model, primal_var): """ dual_constraint = model.private_data().dual_constraint - if primal_var in dual_constraint: - return dual_constraint[primal_var] - else: - raise ValueError( - "It does not appear that Var '%s' is a primal variable on model '%s'" - % (primal_var.name, model.name) - ) + return self._get_corresponding_component(model, 'primal', primal_var, 'Var', + dual_constraint) def get_primal_var(self, model, dual_constraint): """Return the primal variable corresponding to 'dual_constraint' @@ -228,13 +253,8 @@ def get_primal_var(self, model, dual_constraint): """ primal_var = model.private_data().primal_var - if dual_constraint in primal_var: - return primal_var[dual_constraint] - else: - raise ValueError( - "It does not appear that Constraint '%s' is a dual constraint on " - "model '%s'" % (dual_constraint.name, model.name) - ) + return self._get_corresponding_component(model, 'dual', dual_constraint, + 'Constraint', primal_var) def get_dual_var(self, model, primal_constraint): """Return the dual variable corresponding to 'primal_constraint' @@ -252,10 +272,6 @@ def get_dual_var(self, model, primal_constraint): """ dual_var = model.private_data().dual_var - if primal_constraint in dual_var: - return dual_var[primal_constraint] - else: - raise ValueError( - "It does not appear that Constraint '%s' is a primal constraint on " - "model '%s'" % (primal_constraint.name, model.name) - ) + return self._get_corresponding_component(model, 'primal', + primal_constraint, + 'Constraint', dual_var) From 7539e81b5f6e1032eb576bdfba10fd22b37c4cc7 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Wed, 6 Nov 2024 08:43:47 -0700 Subject: [PATCH 45/65] black --- pyomo/core/plugins/transform/lp_dual.py | 35 ++++++++++++++++--------- 1 file changed, 22 insertions(+), 13 deletions(-) diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index 8689faa4158..e9e5fbc0cff 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -164,8 +164,9 @@ def _take_dual(self, model, std_form): return dual - def _get_corresponding_component(self, model, model_type, component, - component_type, mapping): + def _get_corresponding_component( + self, model, model_type, component, component_type, mapping + ): """Return the corresponding component based on the provided mapping. Parameters @@ -195,8 +196,13 @@ def _get_corresponding_component(self, model, model_type, component, else: raise ValueError( "It does not appear that %s '%s' is a %s %s on model '%s'" - % (component_type, component.name, model_type, 'variable' if component_type == 'Var' else 'constraint', - model.name) + % ( + component_type, + component.name, + model_type, + 'variable' if component_type == 'Var' else 'constraint', + model.name, + ) ) def get_primal_constraint(self, model, dual_var): @@ -215,8 +221,9 @@ def get_primal_constraint(self, model, dual_var): """ primal_constraint = model.private_data().primal_constraint - return self._get_corresponding_component(model, 'dual', dual_var, 'Var', - primal_constraint) + return self._get_corresponding_component( + model, 'dual', dual_var, 'Var', primal_constraint + ) def get_dual_constraint(self, model, primal_var): """Return the dual constraint corresponding to 'primal_var' @@ -234,8 +241,9 @@ def get_dual_constraint(self, model, primal_var): """ dual_constraint = model.private_data().dual_constraint - return self._get_corresponding_component(model, 'primal', primal_var, 'Var', - dual_constraint) + return self._get_corresponding_component( + model, 'primal', primal_var, 'Var', dual_constraint + ) def get_primal_var(self, model, dual_constraint): """Return the primal variable corresponding to 'dual_constraint' @@ -253,8 +261,9 @@ def get_primal_var(self, model, dual_constraint): """ primal_var = model.private_data().primal_var - return self._get_corresponding_component(model, 'dual', dual_constraint, - 'Constraint', primal_var) + return self._get_corresponding_component( + model, 'dual', dual_constraint, 'Constraint', primal_var + ) def get_dual_var(self, model, primal_constraint): """Return the dual variable corresponding to 'primal_constraint' @@ -272,6 +281,6 @@ def get_dual_var(self, model, primal_constraint): """ dual_var = model.private_data().dual_var - return self._get_corresponding_component(model, 'primal', - primal_constraint, - 'Constraint', dual_var) + return self._get_corresponding_component( + model, 'primal', primal_constraint, 'Constraint', dual_var + ) From cfb97017d3e162631999a634addb5f736265fe53 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Wed, 6 Nov 2024 11:20:21 -0700 Subject: [PATCH 46/65] Citing the scipy algorithms I'm reimplementing in docstrings --- .../repn/plugins/parameterized_standard_form.py | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index 4367f80e739..6c1cd42d7f3 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -95,7 +95,9 @@ def __eq__(self, other): class _CSRMatrix(_SparseMatrixBase): def tocsc(self): - # Implements the same algorithm as scipy's csr_tocsc function + """Implements the same algorithm as scipy's csr_tocsc function from + sparsetools. + """ csr_data = self.data col_index = self.indices row_index_ptr = self.indptr @@ -141,6 +143,9 @@ def tocsc(self): return _CSCMatrix((csc_data, row_index, col_index_ptr), self.shape) def todense(self): + """Implements the algorithm from scipy's csr_todense function + in sparsetools. + """ nrows = self.shape[0] col_index = self.indices row_index_ptr = self.indptr @@ -157,6 +162,9 @@ def todense(self): class _CSCMatrix(_SparseMatrixBase): def todense(self): + """Implements the algorithm from scipy's csr_todense function + in sparsetools. + """ ncols = self.shape[1] row_index = self.indices col_index_ptr = self.indptr @@ -171,6 +179,9 @@ def todense(self): return dense def sum_duplicates(self): + """Implements the algorithm from scipy's csr_sum_duplicates function + in sparsetools. + """ ncols = self.shape[1] row_index = self.indices col_index_ptr = self.indptr @@ -194,6 +205,9 @@ def sum_duplicates(self): col_index_ptr[i + 1] = num_non_zeros def eliminate_zeros(self): + """Implements the algorithm from scipy's csr_eliminate_zeros function + in sparsetools. + """ ncols = self.shape[1] row_index = self.indices col_index_ptr = self.indptr From 317802c8f0a36415070112170f9ecf75a7df24af Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Wed, 6 Nov 2024 13:55:26 -0700 Subject: [PATCH 47/65] Taking John's suggestions for comments about scipy beauty --- .../plugins/parameterized_standard_form.py | 28 +++++++++++++------ 1 file changed, 19 insertions(+), 9 deletions(-) diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index 6c1cd42d7f3..62fa358cee1 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -116,14 +116,20 @@ def tocsc(self): for i, tally in enumerate(col_index_ptr): col_index_ptr[i] = cum_sum cum_sum += tally - # we leave off the last entry (the number of nonzeros) because we are - # going to do the cute scipy thing and it will get 'added' at the end - # when we shift right (by which I mean it will conveniently already be - # there) + # We have now initialized the col_index_ptr to the *starting* position + # of each column in the data vector. Note that col_index_ptr is only + # num_cols long: we have ignored the last entry in the standard CSC + # col_index_ptr (the total number of nonzeros). This will get resolved + # below when we shift this vector by one position. # Now we are actually going to mess up what we just did while we - # construct the row index, but it's beautiful because "mess up" just - # means we increment everything by one, so we can fix it at the end. + # construct the row index: We can imagine that col_index_ptr holds the + # position of the *next* nonzero in each column, so each time we move a + # data element into a column we will increment that col_index_ptr by + # one. This is beautiful because by "messing up" the col_index_pointer, + # we are just transforming the vector of *starting* indices for each + # column to a vector of *ending* indices (actually 1 past the last + # index) of each column. Thank you, scipy. for row in range(nrows): for j in range(row_index_ptr[row], row_index_ptr[row + 1]): col = col_index[j] @@ -135,9 +141,13 @@ def tocsc(self): col_index_ptr[col] += 1 - # fix the column index pointer by inserting 0 at the beginning--the rest - # is right because we know each entry got incremented by 1 in the loop - # above. + # Fix the column index pointer by inserting 0 at the beginning. The + # col_index_ptr currently holds pointers to 1 past the last element of + # each column, which is really the starting index for the next + # column. Inserting the 0 (the starting index for the firsst column) + # shifts everything by one column, "converting" the vector to the + # starting indices of each column, and extending the vector length to + # num_cols + 1 (as is expected by the CSC matrix). col_index_ptr = np.insert(col_index_ptr, 0, 0) return _CSCMatrix((csc_data, row_index, col_index_ptr), self.shape) From 7a296c8183cedeb04dc9fee702526d1d870643d6 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Fri, 8 Nov 2024 11:33:24 -0700 Subject: [PATCH 48/65] Renaming the parameterized standard form compiler --- pyomo/core/plugins/transform/lp_dual.py | 2 +- pyomo/repn/plugins/parameterized_standard_form.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index e9e5fbc0cff..2dc2bff7213 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -100,7 +100,7 @@ def create_using(self, model, ostream=None, **kwds): model, mixed_form=True, set_sense=None ) else: - std_form = WriterFactory('parameterized_standard_form_compiler').write( + std_form = WriterFactory('compile_parameterized_standard_form').write( model, wrt=config.parameterize_wrt, mixed_form=True, set_sense=None ) return self._take_dual(model, std_form) diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index 62fa358cee1..ee844a0875d 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -25,7 +25,7 @@ @WriterFactory.register( - 'parameterized_standard_form_compiler', + 'compile_parameterized_standard_form', 'Compile an LP to standard form (`min cTx s.t. Ax <= b`) treating some ' 'variables as data (e.g., variables decided by the outer problem in a ' 'bilevel optimization problem).', From 669e6d58591f5b4c49aae128225a244321a8499d Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Sat, 9 Nov 2024 21:25:32 -0700 Subject: [PATCH 49/65] Generalizing the VarData list domain, fixing the error for multiple active objectives --- .../fme/fourier_motzkin_elimination.py | 4 +- pyomo/core/plugins/transform/lp_dual.py | 8 +-- pyomo/core/tests/unit/test_lp_dual.py | 4 +- .../plugins/parameterized_standard_form.py | 9 +++- pyomo/util/config_domains.py | 53 +++++++++++++++++++ pyomo/util/var_list_domain.py | 34 ------------ 6 files changed, 68 insertions(+), 44 deletions(-) create mode 100644 pyomo/util/config_domains.py delete mode 100644 pyomo/util/var_list_domain.py diff --git a/pyomo/contrib/fme/fourier_motzkin_elimination.py b/pyomo/contrib/fme/fourier_motzkin_elimination.py index cec16e656ac..0bf5abf7221 100644 --- a/pyomo/contrib/fme/fourier_motzkin_elimination.py +++ b/pyomo/contrib/fme/fourier_motzkin_elimination.py @@ -30,7 +30,7 @@ from pyomo.repn.standard_repn import generate_standard_repn from pyomo.common.collections import ComponentMap, ComponentSet from pyomo.opt import TerminationCondition -from pyomo.util.var_list_domain import var_component_set +from pyomo.util.config_domains import ComponentDataList import logging @@ -95,7 +95,7 @@ class Fourier_Motzkin_Elimination_Transformation(Transformation): 'vars_to_eliminate', ConfigValue( default=None, - domain=var_component_set, + domain=ComponentDataList(Var), description="Continuous variable or list of continuous variables to " "project out of the model", doc=""" diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index 2dc2bff7213..81c23e1d874 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -29,7 +29,7 @@ ) from pyomo.opt import WriterFactory from pyomo.repn.standard_repn import isclose_const -from pyomo.util.var_list_domain import var_component_set +from pyomo.util.config_domains import ComponentDataList class _LPDualData(AutoSlots.Mixin): @@ -54,7 +54,7 @@ class LinearProgrammingDual(object): 'parameterize_wrt', ConfigValue( default=None, - domain=var_component_set, + domain=ComponentDataList(Var), description="Vars to treat as data for the purposes of taking the dual", doc=""" Optional list of Vars to be treated as data while taking the LP dual. @@ -108,8 +108,8 @@ def create_using(self, model, ostream=None, **kwds): def _take_dual(self, model, std_form): if len(std_form.objectives) != 1: raise ValueError( - "Model '%s' has no objective or multiple active objectives. Cannot " - "take dual with more than one objective!" % model.name + "Model '%s' has no objective or multiple active objectives. Can " + "only take dual with exactly one active objective!" % model.name ) primal_sense = std_form.objectives[0].sense diff --git a/pyomo/core/tests/unit/test_lp_dual.py b/pyomo/core/tests/unit/test_lp_dual.py index 33e3e9287cd..5feacfe5c61 100644 --- a/pyomo/core/tests/unit/test_lp_dual.py +++ b/pyomo/core/tests/unit/test_lp_dual.py @@ -331,7 +331,7 @@ def test_multiple_obj_error(self): with self.assertRaisesRegex( ValueError, "Model 'primal' has no objective or multiple active objectives. " - "Cannot take dual with more than one objective!", + "Can only take dual with exactly one active objective!", ): dual = lp_dual.create_using(m, parameterize_wrt=[m.outer1, m.outer]) @@ -341,7 +341,7 @@ def test_multiple_obj_error(self): with self.assertRaisesRegex( ValueError, "Model 'primal' has no objective or multiple active objectives. " - "Cannot take dual with more than one objective!", + "Can only take dual with exactly one active objective!", ): dual = lp_dual.create_using(m, parameterize_wrt=[m.outer1, m.outer]) diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index ee844a0875d..ed924a241be 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -13,6 +13,7 @@ from pyomo.common.dependencies import numpy as np from pyomo.common.gc_manager import PauseGC from pyomo.common.numeric_types import native_numeric_types +from pyomo.core import Var from pyomo.opt import WriterFactory from pyomo.repn.parameterized_linear import ParameterizedLinearRepnVisitor @@ -21,7 +22,7 @@ LinearStandardFormCompiler, _LinearStandardFormCompiler_impl, ) -from pyomo.util.var_list_domain import var_component_set +from pyomo.util.config_domains import ComponentDataList @WriterFactory.register( @@ -36,7 +37,7 @@ class ParameterizedLinearStandardFormCompiler(LinearStandardFormCompiler): 'wrt', ConfigValue( default=None, - domain=var_component_set, + domain=ComponentDataList(Var), description="Vars to treat as data for the purposes of compiling" "the standard form", doc=""" @@ -191,6 +192,10 @@ def todense(self): def sum_duplicates(self): """Implements the algorithm from scipy's csr_sum_duplicates function in sparsetools. + + Note that this only removes duplicates that are adjacent, so it will remove + all duplicates if the incoming CSC matrix has sorted indices. (In particular + this will be true if it was just convered from CSR). """ ncols = self.shape[1] row_index = self.indices diff --git a/pyomo/util/config_domains.py b/pyomo/util/config_domains.py new file mode 100644 index 00000000000..5353698d338 --- /dev/null +++ b/pyomo/util/config_domains.py @@ -0,0 +1,53 @@ +# ___________________________________________________________________________ +# +# 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. +# ___________________________________________________________________________ + +from pyomo.common.collections import ComponentSet + + +class ComponentDataList(): + """ComponentDataList(ctype) + Domain validation class that accepts singleton or iterable arguments and + compiles them into a ComponentSet, verifying that they are all ComponentDatas + of type 'ctype.' + + Parameters + ---------- + ctype: The component type of the list + + Raises + ------ + ValueError if all of the arguments are not of type 'ctype' + """ + def __init__(self, ctype): + self._ctype = ctype + + def __call__(self, x): + if hasattr(x, 'ctype') and x.ctype is self._ctype: + if not x.is_indexed(): + return ComponentSet([x]) + ans = ComponentSet() + for j in x.index_set(): + ans.add(x[j]) + return ans + elif hasattr(x, '__iter__'): + ans = ComponentSet() + for i in x: + ans.update(self(i)) + return ans + else: + _ctype_name = str(self._ctype) + raise ValueError( + f"Expected {_ctype_name} or iterable of " + f"{_ctype_name}s.\n\tReceived {type(x)}") + + def domain_name(self): + _ctype_name = str(self._ctype) + return f'ComponentDataList({_ctype_name})' diff --git a/pyomo/util/var_list_domain.py b/pyomo/util/var_list_domain.py deleted file mode 100644 index 15c0f775b57..00000000000 --- a/pyomo/util/var_list_domain.py +++ /dev/null @@ -1,34 +0,0 @@ -# ___________________________________________________________________________ -# -# 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. -# ___________________________________________________________________________ - -from pyomo.common.collections import ComponentSet -from pyomo.core import Var - - -def var_component_set(x): - """ - For domain validation in ConfigDicts: Takes singletone or iterable argument 'x' - of Vars and converts it to a ComponentSet of Vars. - """ - if hasattr(x, 'ctype') and x.ctype is Var: - if not x.is_indexed(): - return ComponentSet([x]) - ans = ComponentSet() - for j in x.index_set(): - ans.add(x[j]) - return ans - elif hasattr(x, '__iter__'): - ans = ComponentSet() - for i in x: - ans.update(var_component_set(i)) - return ans - else: - raise ValueError("Expected Var or iterable of Vars.\n\tReceived %s" % type(x)) From 4b8ec54904881d2a22a4a14d6fa855ce614b5df6 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Sat, 9 Nov 2024 21:32:28 -0700 Subject: [PATCH 50/65] Emma learns numpy--making c vector a vector... --- pyomo/core/plugins/transform/lp_dual.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index 81c23e1d874..8829ac31325 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -117,7 +117,7 @@ def _take_dual(self, model, std_form): # This is a csc matrix, so we'll skip transposing and just work off # of the columns A = std_form.A - c = std_form.c.todense() + c = std_form.c.todense().ravel() dual_rows = range(A.shape[1]) dual_cols = range(A.shape[0]) dual.x = Var(dual_cols, domain=NonNegativeReals) @@ -143,17 +143,17 @@ def _take_dual(self, model, std_form): lhs += coef * dual.x[primal_row] if primal.domain is Reals: - dual.constraints[i] = lhs == c[0, i] + dual.constraints[i] = lhs == c[i] elif primal_sense is minimize: if primal.domain is NonNegativeReals: - dual.constraints[i] = lhs <= c[0, i] + dual.constraints[i] = lhs <= c[i] else: # primal.domain is NonPositiveReals - dual.constraints[i] = lhs >= c[0, i] + dual.constraints[i] = lhs >= c[i] else: if primal.domain is NonNegativeReals: - dual.constraints[i] = lhs >= c[0, i] + dual.constraints[i] = lhs >= c[i] else: # primal.domain is NonPositiveReals - dual.constraints[i] = lhs <= c[0, i] + dual.constraints[i] = lhs <= c[i] trans_info.dual_constraint[primal] = dual.constraints[i] trans_info.primal_var[dual.constraints[i]] = primal From 9e0705716643d570369bfacc1011f2a476e72814 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Sat, 9 Nov 2024 21:34:06 -0700 Subject: [PATCH 51/65] Black --- pyomo/util/config_domains.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pyomo/util/config_domains.py b/pyomo/util/config_domains.py index 5353698d338..4a7dac0934b 100644 --- a/pyomo/util/config_domains.py +++ b/pyomo/util/config_domains.py @@ -12,7 +12,7 @@ from pyomo.common.collections import ComponentSet -class ComponentDataList(): +class ComponentDataList: """ComponentDataList(ctype) Domain validation class that accepts singleton or iterable arguments and compiles them into a ComponentSet, verifying that they are all ComponentDatas @@ -26,6 +26,7 @@ class ComponentDataList(): ------ ValueError if all of the arguments are not of type 'ctype' """ + def __init__(self, ctype): self._ctype = ctype @@ -46,7 +47,8 @@ def __call__(self, x): _ctype_name = str(self._ctype) raise ValueError( f"Expected {_ctype_name} or iterable of " - f"{_ctype_name}s.\n\tReceived {type(x)}") + f"{_ctype_name}s.\n\tReceived {type(x)}" + ) def domain_name(self): _ctype_name = str(self._ctype) From 1c489824001b3361afcdf81d8b1df4c6112fb55b Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Sat, 9 Nov 2024 21:40:02 -0700 Subject: [PATCH 52/65] Reverting the generalization of the dual-to-primal mapping API --- pyomo/core/plugins/transform/lp_dual.py | 81 +++++++++---------------- 1 file changed, 28 insertions(+), 53 deletions(-) diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index 8829ac31325..5e916d66435 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -164,47 +164,6 @@ def _take_dual(self, model, std_form): return dual - def _get_corresponding_component( - self, model, model_type, component, component_type, mapping - ): - """Return the corresponding component based on the provided mapping. - - Parameters - ---------- - model: ConcreteModel - The model containing the mappings. - component: Var or Constraint - The primal/dual component for which we want to find the corresponding - dual/primal component. - component_type: str - A string indicating whether the component is 'dual' or 'primal'. - mapping: dict - The mapping to look up the corresponding component. - - Returns - ------- - Var or Constraint - The corresponding component. - - Raises - ------ - ValueError - If the component is not found in the mapping. - """ - if component in mapping: - return mapping[component] - else: - raise ValueError( - "It does not appear that %s '%s' is a %s %s on model '%s'" - % ( - component_type, - component.name, - model_type, - 'variable' if component_type == 'Var' else 'constraint', - model.name, - ) - ) - def get_primal_constraint(self, model, dual_var): """Return the primal constraint corresponding to 'dual_var' @@ -221,9 +180,13 @@ def get_primal_constraint(self, model, dual_var): """ primal_constraint = model.private_data().primal_constraint - return self._get_corresponding_component( - model, 'dual', dual_var, 'Var', primal_constraint - ) + if dual_var in primal_constraint: + return primal_constraint[dual_var] + else: + raise ValueError( + "It does not appear that Var '%s' is a dual variable on model '%s'" + % (dual_var.name, model.name) + ) def get_dual_constraint(self, model, primal_var): """Return the dual constraint corresponding to 'primal_var' @@ -241,9 +204,13 @@ def get_dual_constraint(self, model, primal_var): """ dual_constraint = model.private_data().dual_constraint - return self._get_corresponding_component( - model, 'primal', primal_var, 'Var', dual_constraint - ) + if primal_var in dual_constraint: + return dual_constraint[primal_var] + else: + raise ValueError( + "It does not appear that Var '%s' is a primal variable on model '%s'" + % (primal_var.name, model.name) + ) def get_primal_var(self, model, dual_constraint): """Return the primal variable corresponding to 'dual_constraint' @@ -261,9 +228,13 @@ def get_primal_var(self, model, dual_constraint): """ primal_var = model.private_data().primal_var - return self._get_corresponding_component( - model, 'dual', dual_constraint, 'Constraint', primal_var - ) + if dual_constraint in primal_var: + return primal_var[dual_constraint] + else: + raise ValueError( + "It does not appear that Constraint '%s' is a dual constraint on " + "model '%s'" % (dual_constraint.name, model.name) + ) def get_dual_var(self, model, primal_constraint): """Return the dual variable corresponding to 'primal_constraint' @@ -281,6 +252,10 @@ def get_dual_var(self, model, primal_constraint): """ dual_var = model.private_data().dual_var - return self._get_corresponding_component( - model, 'primal', primal_constraint, 'Constraint', dual_var - ) + if primal_constraint in dual_var: + return dual_var[primal_constraint] + else: + raise ValueError( + "It does not appear that Constraint '%s' is a primal constraint on " + "model '%s'" % (primal_constraint.name, model.name) + ) From f3e4acd43b75549f65a2d70b695798d3715e2bf9 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Sat, 9 Nov 2024 21:42:25 -0700 Subject: [PATCH 53/65] NFC: fixing comment typo --- pyomo/repn/plugins/parameterized_standard_form.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index ed924a241be..04ce6be601f 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -195,7 +195,7 @@ def sum_duplicates(self): Note that this only removes duplicates that are adjacent, so it will remove all duplicates if the incoming CSC matrix has sorted indices. (In particular - this will be true if it was just convered from CSR). + this will be true if it was just converted from CSR). """ ncols = self.shape[1] row_index = self.indices From 25749b02bc0b371d1845f4cd522f1cafdec6d51a Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Sat, 9 Nov 2024 21:45:24 -0700 Subject: [PATCH 54/65] Switching apply_to error to NotImplementedError --- pyomo/core/plugins/transform/lp_dual.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index 5e916d66435..a9c9f61d561 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -12,7 +12,6 @@ from pyomo.common.autoslots import AutoSlots from pyomo.common.collections import ComponentMap from pyomo.common.config import ConfigDict, ConfigValue -from pyomo.common.errors import MouseTrap from pyomo.common.dependencies import scipy from pyomo.core import ( ConcreteModel, @@ -68,9 +67,9 @@ class LinearProgrammingDual(object): ) def apply_to(self, model, **options): - raise MouseTrap( - "The 'core.lp_dual' transformation does not currently implement " - "apply_to since it is a bit ambiguous what it means to take a dual " + raise NotImplementedError( + "The 'core.lp_dual' transformation does not implement " + "apply_to since it is ambiguous what it means to take a dual " "in place. Please use 'create_using' and do what you wish with the " "returned model." ) From d9c2a7375b8d053ada0d8c6415807583b4eb8815 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 11 Nov 2024 13:07:18 -0700 Subject: [PATCH 55/65] Adding a test for sum_duplicates --- pyomo/repn/plugins/parameterized_standard_form.py | 5 +++++ .../tests/test_parameterized_standard_form.py | 15 +++++++++++++++ 2 files changed, 20 insertions(+) diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index 04ce6be601f..1235914e63c 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -218,6 +218,11 @@ def sum_duplicates(self): data[num_non_zeros] = x num_non_zeros += 1 col_index_ptr[i + 1] = num_non_zeros + + # [ESJ 11/11/24]: I'm not 100% sure how scipy handles this, but we need + # to remove the "extra" entries from the data and row_index arrays. + self.data = data[:num_non_zeros] + self.row_index = row_index[:num_non_zeros] def eliminate_zeros(self): """Implements the algorithm from scipy's csr_eliminate_zeros function diff --git a/pyomo/repn/tests/test_parameterized_standard_form.py b/pyomo/repn/tests/test_parameterized_standard_form.py index 0b92071bc26..05c7d796bc4 100644 --- a/pyomo/repn/tests/test_parameterized_standard_form.py +++ b/pyomo/repn/tests/test_parameterized_standard_form.py @@ -93,6 +93,21 @@ def test_todense(self): self.assertTrue(np.all(A.todense() == dense)) self.assertTrue(np.all(A.tocsc().todense() == dense)) + def test_sum_duplicates(self): + A = _CSCMatrix(([4, 5], [1, 1], [0, 0, 2, 2]), [3, 3]) + self.assertTrue(np.all(A.data == [4, 5])) + self.assertTrue(np.all(A.indptr == [0, 0, 2, 2])) + self.assertTrue(np.all(A.indices == [1, 1])) + + A.sum_duplicates() + + self.assertTrue(np.all(A.data == [9,])) + self.assertTrue(np.all(A.indptr == [0, 0, 1, 1])) + self.assertTrue(np.all(A.indices == [1,])) + + dense = np.array([[0, 0, 0], [0, 9, 0], [0, 0, 0]]) + self.assertTrue(np.all(A.todense() == dense)) + def assertExpressionArraysEqual(self, A, B): self.assertEqual(A.shape, B.shape) From fc5d5dd3868c3b703f315c0abd8801f9a2a8d9f1 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 11 Nov 2024 13:08:14 -0700 Subject: [PATCH 56/65] black --- pyomo/repn/plugins/parameterized_standard_form.py | 2 +- pyomo/repn/tests/test_parameterized_standard_form.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index 1235914e63c..174e86a6ffb 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -218,7 +218,7 @@ def sum_duplicates(self): data[num_non_zeros] = x num_non_zeros += 1 col_index_ptr[i + 1] = num_non_zeros - + # [ESJ 11/11/24]: I'm not 100% sure how scipy handles this, but we need # to remove the "extra" entries from the data and row_index arrays. self.data = data[:num_non_zeros] diff --git a/pyomo/repn/tests/test_parameterized_standard_form.py b/pyomo/repn/tests/test_parameterized_standard_form.py index 05c7d796bc4..dbcfaf8310b 100644 --- a/pyomo/repn/tests/test_parameterized_standard_form.py +++ b/pyomo/repn/tests/test_parameterized_standard_form.py @@ -100,10 +100,10 @@ def test_sum_duplicates(self): self.assertTrue(np.all(A.indices == [1, 1])) A.sum_duplicates() - - self.assertTrue(np.all(A.data == [9,])) + + self.assertTrue(np.all(A.data == [9])) self.assertTrue(np.all(A.indptr == [0, 0, 1, 1])) - self.assertTrue(np.all(A.indices == [1,])) + self.assertTrue(np.all(A.indices == [1])) dense = np.array([[0, 0, 0], [0, 9, 0], [0, 0, 0]]) self.assertTrue(np.all(A.todense() == dense)) From 2a5e1aaf8e83c2a63b9b1d5e0ae3179e4ca44552 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 11 Nov 2024 13:31:40 -0700 Subject: [PATCH 57/65] Adding some input validation and a test for it --- .../plugins/parameterized_standard_form.py | 20 +++++++++++++++++++ .../tests/test_parameterized_standard_form.py | 17 ++++++++++++++++ 2 files changed, 37 insertions(+) diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index 174e86a6ffb..cb26572d2d6 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -95,6 +95,16 @@ def __eq__(self, other): class _CSRMatrix(_SparseMatrixBase): + def __init__(self, matrix_data, shape): + super().__init__(matrix_data, shape) + if len(self.indptr) != self.shape[0] + 1: + raise ValueError( + "Shape specifies the number of rows as %s but the index " + "pointer has length %s. The index pointer must have length " + "nrows + 1: Check the 'shape' and 'matrix_data' arguments." + % (self.shape[0], len(self.indptr)) + ) + def tocsc(self): """Implements the same algorithm as scipy's csr_tocsc function from sparsetools. @@ -172,6 +182,16 @@ def todense(self): class _CSCMatrix(_SparseMatrixBase): + def __init__(self, matrix_data, shape): + super().__init__(matrix_data, shape) + if len(self.indptr) != self.shape[1] + 1: + raise ValueError( + "Shape specifies the number of columns as %s but the index " + "pointer has length %s. The index pointer must have length " + "ncols + 1: Check the 'shape' and 'matrix_data' arguments." + % (self.shape[1], len(self.indptr)) + ) + def todense(self): """Implements the algorithm from scipy's csr_todense function in sparsetools. diff --git a/pyomo/repn/tests/test_parameterized_standard_form.py b/pyomo/repn/tests/test_parameterized_standard_form.py index dbcfaf8310b..b23ccfed955 100644 --- a/pyomo/repn/tests/test_parameterized_standard_form.py +++ b/pyomo/repn/tests/test_parameterized_standard_form.py @@ -108,6 +108,23 @@ def test_sum_duplicates(self): dense = np.array([[0, 0, 0], [0, 9, 0], [0, 0, 0]]) self.assertTrue(np.all(A.todense() == dense)) + def test_invalid_sparse_matrix_input(self): + with self.assertRaisesRegex( + ValueError, + r"Shape specifies the number of rows as 3 but the index " + r"pointer has length 2. The index pointer must have length " + r"nrows \+ 1: Check the 'shape' and 'matrix_data' arguments." + ): + A = _CSRMatrix(([4, 5], [1, 1], [1, 1]), shape=(3, 3)) + + with self.assertRaisesRegex( + ValueError, + r"Shape specifies the number of columns as 3 but the index " + r"pointer has length 2. The index pointer must have length " + r"ncols \+ 1: Check the 'shape' and 'matrix_data' arguments." + ): + A = _CSCMatrix(([4, 5], [1, 1], [1, 1]), shape=(3, 3)) + def assertExpressionArraysEqual(self, A, B): self.assertEqual(A.shape, B.shape) From 42ab07522626d4f2babe5d9cd26c890cf10ecda3 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 11 Nov 2024 13:32:05 -0700 Subject: [PATCH 58/65] black --- pyomo/repn/tests/test_parameterized_standard_form.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyomo/repn/tests/test_parameterized_standard_form.py b/pyomo/repn/tests/test_parameterized_standard_form.py index b23ccfed955..4a88ea4f883 100644 --- a/pyomo/repn/tests/test_parameterized_standard_form.py +++ b/pyomo/repn/tests/test_parameterized_standard_form.py @@ -113,7 +113,7 @@ def test_invalid_sparse_matrix_input(self): ValueError, r"Shape specifies the number of rows as 3 but the index " r"pointer has length 2. The index pointer must have length " - r"nrows \+ 1: Check the 'shape' and 'matrix_data' arguments." + r"nrows \+ 1: Check the 'shape' and 'matrix_data' arguments.", ): A = _CSRMatrix(([4, 5], [1, 1], [1, 1]), shape=(3, 3)) @@ -121,7 +121,7 @@ def test_invalid_sparse_matrix_input(self): ValueError, r"Shape specifies the number of columns as 3 but the index " r"pointer has length 2. The index pointer must have length " - r"ncols \+ 1: Check the 'shape' and 'matrix_data' arguments." + r"ncols \+ 1: Check the 'shape' and 'matrix_data' arguments.", ): A = _CSCMatrix(([4, 5], [1, 1], [1, 1]), shape=(3, 3)) From b694bbe8b5b6f8f9d0321217a17f1643d10f8767 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 11 Nov 2024 16:37:00 -0700 Subject: [PATCH 59/65] Changing ComponentDataList to ComponentDataSet and taking most of John's suggestions --- .../fme/fourier_motzkin_elimination.py | 4 +- pyomo/core/plugins/transform/lp_dual.py | 4 +- .../plugins/parameterized_standard_form.py | 4 +- pyomo/util/config_domains.py | 57 +++++++++++-------- 4 files changed, 40 insertions(+), 29 deletions(-) diff --git a/pyomo/contrib/fme/fourier_motzkin_elimination.py b/pyomo/contrib/fme/fourier_motzkin_elimination.py index 0bf5abf7221..021650e8f9a 100644 --- a/pyomo/contrib/fme/fourier_motzkin_elimination.py +++ b/pyomo/contrib/fme/fourier_motzkin_elimination.py @@ -30,7 +30,7 @@ from pyomo.repn.standard_repn import generate_standard_repn from pyomo.common.collections import ComponentMap, ComponentSet from pyomo.opt import TerminationCondition -from pyomo.util.config_domains import ComponentDataList +from pyomo.util.config_domains import ComponentDataSet import logging @@ -95,7 +95,7 @@ class Fourier_Motzkin_Elimination_Transformation(Transformation): 'vars_to_eliminate', ConfigValue( default=None, - domain=ComponentDataList(Var), + domain=ComponentDataSet(Var), description="Continuous variable or list of continuous variables to " "project out of the model", doc=""" diff --git a/pyomo/core/plugins/transform/lp_dual.py b/pyomo/core/plugins/transform/lp_dual.py index a9c9f61d561..82f27a879ea 100644 --- a/pyomo/core/plugins/transform/lp_dual.py +++ b/pyomo/core/plugins/transform/lp_dual.py @@ -28,7 +28,7 @@ ) from pyomo.opt import WriterFactory from pyomo.repn.standard_repn import isclose_const -from pyomo.util.config_domains import ComponentDataList +from pyomo.util.config_domains import ComponentDataSet class _LPDualData(AutoSlots.Mixin): @@ -53,7 +53,7 @@ class LinearProgrammingDual(object): 'parameterize_wrt', ConfigValue( default=None, - domain=ComponentDataList(Var), + domain=ComponentDataSet(Var), description="Vars to treat as data for the purposes of taking the dual", doc=""" Optional list of Vars to be treated as data while taking the LP dual. diff --git a/pyomo/repn/plugins/parameterized_standard_form.py b/pyomo/repn/plugins/parameterized_standard_form.py index cb26572d2d6..541b859caab 100644 --- a/pyomo/repn/plugins/parameterized_standard_form.py +++ b/pyomo/repn/plugins/parameterized_standard_form.py @@ -22,7 +22,7 @@ LinearStandardFormCompiler, _LinearStandardFormCompiler_impl, ) -from pyomo.util.config_domains import ComponentDataList +from pyomo.util.config_domains import ComponentDataSet @WriterFactory.register( @@ -37,7 +37,7 @@ class ParameterizedLinearStandardFormCompiler(LinearStandardFormCompiler): 'wrt', ConfigValue( default=None, - domain=ComponentDataList(Var), + domain=ComponentDataSet(Var), description="Vars to treat as data for the purposes of compiling" "the standard form", doc=""" diff --git a/pyomo/util/config_domains.py b/pyomo/util/config_domains.py index 4a7dac0934b..c7bf407c4f6 100644 --- a/pyomo/util/config_domains.py +++ b/pyomo/util/config_domains.py @@ -10,46 +10,57 @@ # ___________________________________________________________________________ from pyomo.common.collections import ComponentSet +from typing import Sequence -class ComponentDataList: - """ComponentDataList(ctype) +class ComponentDataSet: + """ComponentDataSet(ctype) Domain validation class that accepts singleton or iterable arguments and compiles them into a ComponentSet, verifying that they are all ComponentDatas of type 'ctype.' Parameters ---------- - ctype: The component type of the list + ctype: The component type(s) of the list Raises ------ - ValueError if all of the arguments are not of type 'ctype' + ValueError if all of the arguments are not of a type in 'ctype' """ - def __init__(self, ctype): - self._ctype = ctype + if isinstance(ctype, Sequence): + self._ctypes = set(ctype) + else: + self._ctypes = set([ctype]) def __call__(self, x): - if hasattr(x, 'ctype') and x.ctype is self._ctype: - if not x.is_indexed(): - return ComponentSet([x]) - ans = ComponentSet() - for j in x.index_set(): - ans.add(x[j]) - return ans - elif hasattr(x, '__iter__'): - ans = ComponentSet() - for i in x: - ans.update(self(i)) - return ans + return ComponentSet(self._process(x)) + + def _process(self, x): + _names = ', '.join(ct.__name__ for ct in self._ctypes) + if hasattr(x, 'ctype'): + if x.ctype not in self._ctypes: + raise ValueError( + f"Expected component or iterable of one " + f"of the following ctypes: " + f"{_names}.\n\tReceived {type(x)}" + ) + if x.is_indexed(): + yield from x.values() + else: + yield x + elif isinstance(x, (Sequence, ComponentSet)): + for y in x: + yield from self._process(y) else: - _ctype_name = str(self._ctype) raise ValueError( - f"Expected {_ctype_name} or iterable of " - f"{_ctype_name}s.\n\tReceived {type(x)}" + f"Expected component or iterable of one " + f"of the following ctypes: " + f"{_names}.\n\tReceived {type(x)}" ) def domain_name(self): - _ctype_name = str(self._ctype) - return f'ComponentDataList({_ctype_name})' + _names = ', '.join(ct.__name__ for ct in self._ctypes) + if len(self._ctypes) > 1: + _names = '[' + _names + ']' + return f"ComponentDataSet({_names})" From ac49dda69ef701f7bb46bf9d738159b41d4f6adb Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Mon, 11 Nov 2024 16:37:39 -0700 Subject: [PATCH 60/65] black --- pyomo/util/config_domains.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pyomo/util/config_domains.py b/pyomo/util/config_domains.py index c7bf407c4f6..b5f45ffe6cc 100644 --- a/pyomo/util/config_domains.py +++ b/pyomo/util/config_domains.py @@ -27,6 +27,7 @@ class ComponentDataSet: ------ ValueError if all of the arguments are not of a type in 'ctype' """ + def __init__(self, ctype): if isinstance(ctype, Sequence): self._ctypes = set(ctype) From d0b813fe0525078591849c012cde22361909eb1a Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Tue, 12 Nov 2024 09:29:06 -0700 Subject: [PATCH 61/65] NFC: Updating docstring for ComponentDataSet --- pyomo/util/config_domains.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyomo/util/config_domains.py b/pyomo/util/config_domains.py index b5f45ffe6cc..e5131748aba 100644 --- a/pyomo/util/config_domains.py +++ b/pyomo/util/config_domains.py @@ -21,7 +21,8 @@ class ComponentDataSet: Parameters ---------- - ctype: The component type(s) of the list + ctype: Either a single component type an iterable of component types to + admit in the set Raises ------ From 729631d7c90f8d460b2750fafec143d1b085e0f3 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Tue, 12 Nov 2024 13:42:42 -0700 Subject: [PATCH 62/65] Adding tests for the ComponentDataSet config domain --- pyomo/util/config_domains.py | 8 +- pyomo/util/tests/test_config_domains.py | 107 ++++++++++++++++++++++++ 2 files changed, 113 insertions(+), 2 deletions(-) create mode 100644 pyomo/util/tests/test_config_domains.py diff --git a/pyomo/util/config_domains.py b/pyomo/util/config_domains.py index e5131748aba..73d0c6aba9c 100644 --- a/pyomo/util/config_domains.py +++ b/pyomo/util/config_domains.py @@ -39,7 +39,9 @@ def __call__(self, x): return ComponentSet(self._process(x)) def _process(self, x): - _names = ', '.join(ct.__name__ for ct in self._ctypes) + # Ordering for determinism + _ctypes = sorted([ct.__name__ for ct in self._ctypes]) + _names = ', '.join(_ctypes) if hasattr(x, 'ctype'): if x.ctype not in self._ctypes: raise ValueError( @@ -62,7 +64,9 @@ def _process(self, x): ) def domain_name(self): - _names = ', '.join(ct.__name__ for ct in self._ctypes) + # Ordering for determinism + _ctypes = sorted([ct.__name__ for ct in self._ctypes]) + _names = ', '.join(_ctypes) if len(self._ctypes) > 1: _names = '[' + _names + ']' return f"ComponentDataSet({_names})" diff --git a/pyomo/util/tests/test_config_domains.py b/pyomo/util/tests/test_config_domains.py new file mode 100644 index 00000000000..8660a4802d5 --- /dev/null +++ b/pyomo/util/tests/test_config_domains.py @@ -0,0 +1,107 @@ +# ___________________________________________________________________________ +# +# 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. +# ___________________________________________________________________________ + +from pyomo.common.collections import ComponentSet +from pyomo.common.config import ConfigDict, ConfigValue +import pyomo.common.unittest as unittest +from pyomo.core import ( + Block, + ConcreteModel, + Constraint, + Objective, + Var, +) +from pyomo.util.config_domains import ComponentDataSet + +def ComponentSetConfig(): + CONFIG = ConfigDict() + CONFIG.declare( + 'var_set', + ConfigValue( + default=None, + domain=ComponentDataSet(Var), + doc="VarDataSet" + ) + ) + CONFIG.declare( + 'var_and_constraint_set', + ConfigValue( + default=None, + domain=ComponentDataSet(ctype=(Var, Constraint)), + doc="VarAndConstraintSet" + ) + ) + return CONFIG + + +def a_model(): + m = ConcreteModel() + m.x = Var() + m.y = Var([1, 2]) + m.c = Constraint(expr=m.x + m.y[1] <= 3) + m.c2 = Constraint(expr=m.y[2] >= 7) + m.obj = Objective(expr=m.x + m.y[1] + m.y[2]) + m.b = Block() + + return m + + +class TestComponentDataSetDomain(unittest.TestCase): + def test_var_set(self): + m = a_model() + config = ComponentSetConfig() + self.assertIsNone(config.var_set) + config.var_set = ComponentSet([m.x, m.y]) + self.assertIsInstance(config.var_set, ComponentSet) + self.assertEqual(len(config.var_set), 3) + for v in [m.x, m.y[1], m.y[2]]: + self.assertIn(v, config.var_set) + + with self.assertRaisesRegex( + ValueError, + ".*Expected component or iterable of one " + "of the following ctypes: Var.\n\t" + "Received " + ): + config.var_set = ComponentSet([m.y, m.c]) + + def test_var_and_constraint_set(self): + m = a_model() + config = ComponentSetConfig() + self.assertIsNone(config.var_and_constraint_set) + config.var_and_constraint_set = ComponentSet([m.x, m.c]) + self.assertIsInstance(config.var_and_constraint_set, ComponentSet) + self.assertEqual(len(config.var_and_constraint_set), 2) + for v in [m.x, m.c]: + self.assertIn(v, config.var_and_constraint_set) + + with self.assertRaisesRegex( + ValueError, + ".*Expected component or iterable of one " + "of the following ctypes: Constraint, Var.\n\t" + "Received " + ): + config.var_and_constraint_set = ComponentSet([m.y, m.c, m.b]) + + with self.assertRaisesRegex( + ValueError, + ".*Expected component or iterable of one " + "of the following ctypes: Constraint, Var.\n\t" + "Received " + ): + config.var_and_constraint_set = ComponentSet([3, m.y, m.c]) + + def test_domain_name(self): + config = ComponentSetConfig() + self.assertEqual(config.get("var_set").domain_name(), + "ComponentDataSet(Var)") + self.assertEqual(config.get("var_and_constraint_set").domain_name(), + "ComponentDataSet([Constraint, Var])") From 51c57d265d32abd3fba707308b67b872d78565ea Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Tue, 12 Nov 2024 13:43:57 -0700 Subject: [PATCH 63/65] Black --- pyomo/util/tests/test_config_domains.py | 34 ++++++++++--------------- 1 file changed, 13 insertions(+), 21 deletions(-) diff --git a/pyomo/util/tests/test_config_domains.py b/pyomo/util/tests/test_config_domains.py index 8660a4802d5..9955950fcf5 100644 --- a/pyomo/util/tests/test_config_domains.py +++ b/pyomo/util/tests/test_config_domains.py @@ -12,32 +12,23 @@ from pyomo.common.collections import ComponentSet from pyomo.common.config import ConfigDict, ConfigValue import pyomo.common.unittest as unittest -from pyomo.core import ( - Block, - ConcreteModel, - Constraint, - Objective, - Var, -) +from pyomo.core import Block, ConcreteModel, Constraint, Objective, Var from pyomo.util.config_domains import ComponentDataSet + def ComponentSetConfig(): CONFIG = ConfigDict() CONFIG.declare( 'var_set', - ConfigValue( - default=None, - domain=ComponentDataSet(Var), - doc="VarDataSet" - ) + ConfigValue(default=None, domain=ComponentDataSet(Var), doc="VarDataSet"), ) CONFIG.declare( 'var_and_constraint_set', ConfigValue( default=None, domain=ComponentDataSet(ctype=(Var, Constraint)), - doc="VarAndConstraintSet" - ) + doc="VarAndConstraintSet", + ), ) return CONFIG @@ -69,7 +60,7 @@ def test_var_set(self): ValueError, ".*Expected component or iterable of one " "of the following ctypes: Var.\n\t" - "Received " + "Received ", ): config.var_set = ComponentSet([m.y, m.c]) @@ -87,7 +78,7 @@ def test_var_and_constraint_set(self): ValueError, ".*Expected component or iterable of one " "of the following ctypes: Constraint, Var.\n\t" - "Received " + "Received ", ): config.var_and_constraint_set = ComponentSet([m.y, m.c, m.b]) @@ -95,13 +86,14 @@ def test_var_and_constraint_set(self): ValueError, ".*Expected component or iterable of one " "of the following ctypes: Constraint, Var.\n\t" - "Received " + "Received ", ): config.var_and_constraint_set = ComponentSet([3, m.y, m.c]) def test_domain_name(self): config = ComponentSetConfig() - self.assertEqual(config.get("var_set").domain_name(), - "ComponentDataSet(Var)") - self.assertEqual(config.get("var_and_constraint_set").domain_name(), - "ComponentDataSet([Constraint, Var])") + self.assertEqual(config.get("var_set").domain_name(), "ComponentDataSet(Var)") + self.assertEqual( + config.get("var_and_constraint_set").domain_name(), + "ComponentDataSet([Constraint, Var])", + ) From 5ec2dc78dbfe6abdf58ab56176dbe7b7dd15f193 Mon Sep 17 00:00:00 2001 From: Emma Johnson Date: Tue, 12 Nov 2024 13:48:36 -0700 Subject: [PATCH 64/65] NFC: Fixing a typo in the ComponentDataSet docstring --- pyomo/util/config_domains.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pyomo/util/config_domains.py b/pyomo/util/config_domains.py index 73d0c6aba9c..7daf97b794e 100644 --- a/pyomo/util/config_domains.py +++ b/pyomo/util/config_domains.py @@ -21,8 +21,7 @@ class ComponentDataSet: Parameters ---------- - ctype: Either a single component type an iterable of component types to - admit in the set + ctype: Either a single component type or an iterable of component types Raises ------ From 5e8913f8a214acabfcac584370fde78d0d16b1a7 Mon Sep 17 00:00:00 2001 From: John Siirola Date: Wed, 13 Nov 2024 01:17:59 -0700 Subject: [PATCH 65/65] Update ComponentDataSet for efficiency, more permissible inputs --- pyomo/util/config_domains.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/pyomo/util/config_domains.py b/pyomo/util/config_domains.py index 7daf97b794e..86a38bea37b 100644 --- a/pyomo/util/config_domains.py +++ b/pyomo/util/config_domains.py @@ -38,11 +38,10 @@ def __call__(self, x): return ComponentSet(self._process(x)) def _process(self, x): - # Ordering for determinism - _ctypes = sorted([ct.__name__ for ct in self._ctypes]) - _names = ', '.join(_ctypes) if hasattr(x, 'ctype'): if x.ctype not in self._ctypes: + # Ordering for determinism + _names = ', '.join(sorted([ct.__name__ for ct in self._ctypes])) raise ValueError( f"Expected component or iterable of one " f"of the following ctypes: " @@ -52,10 +51,12 @@ def _process(self, x): yield from x.values() else: yield x - elif isinstance(x, (Sequence, ComponentSet)): + elif hasattr(x, '__iter__'): for y in x: yield from self._process(y) else: + # Ordering for determinism + _names = ', '.join(sorted([ct.__name__ for ct in self._ctypes])) raise ValueError( f"Expected component or iterable of one " f"of the following ctypes: "