[WIP] Add a QuadraticExpression node

pyomo/core/expr/compare.py

@@ -16,7 +16,7 @@
     NPV_ProductExpression, NPV_DivisionExpression, NPV_PowExpression,
     NPV_SumExpression, NPV_NegationExpression, NPV_UnaryFunctionExpression,
     NPV_ExternalFunctionExpression, Expr_ifExpression, AbsExpression,
-    NPV_AbsExpression, NumericValue)
+    NPV_AbsExpression, NumericValue, QuadraticExpression)
 from pyomo.core.expr.logical_expr import (
     RangedExpression, InequalityExpression, EqualityExpression
 )
@@ -32,6 +32,14 @@ def handle_linear_expression(node: LinearExpression, pn: List):
     return tuple()
 
 
+def handle_quadratic_expression(node: QuadraticExpression, pn: List):
+    pn.append((type(node), 3*node.nargs()))
+    pn.extend(node.coefs)
+    pn.extend(node.vars_1)
+    pn.extend(node.vars_2)
+    return tuple()
+
+
 def handle_expression(node: ExpressionBase, pn: List):
     pn.append((type(node), node.nargs()))
     return node.args
@@ -76,6 +84,7 @@ def handle_external_function_expression(node: ExternalFunctionExpression, pn: Li
 handler[RangedExpression] = handle_expression
 handler[InequalityExpression] = handle_expression
 handler[EqualityExpression] = handle_expression
+handler[QuadraticExpression] = handle_quadratic_expression
 
 
 class PrefixVisitor(StreamBasedExpressionVisitor):

pyomo/core/expr/current.py

@@ -66,6 +66,7 @@ class Mode(object):
                                               NPV_UnaryFunctionExpression,
                                               AbsExpression, NPV_AbsExpression,
                                               LinearExpression,
+                                              QuadraticExpression,
                                               _MutableLinearExpression,
                                               decompose_term,
                                               LinearDecompositionError,

pyomo/core/expr/numeric_expr.py

@@ -1518,6 +1518,160 @@ def _combine_expr(self, etype, _other):
         return self
 
 
+class QuadraticExpression(ExpressionBase):
+    """
+    An expression object for quadratic expressions
+
+    Args:
+        args: tuple or list
+            The children nodes
+    """
+    __slots__ = ('_coefs', '_vars_1', '_vars_2', '_args_cache_')
+
+    PRECEDENCE = 6
+
+    def __init__(self, args=None, coefs=None, vars_1=None, vars_2=None):
+        """
+        A quadratic expression of the form sum_i(c_i*x_i*y_i).
+
+        You can specify args OR (coefs, vars_1, and vars_2). If args is provided, it
+        should be a list or tuple that contains a series of
+        :py:class:`ProductExpression` objects with arguments
+        :py:class:`MonomialTermExpression` and :py:class:`_GeneralVarData`.
+        Alternatively, you can specify the lists of coefficients and variables
+        separately.
+        """
+        if args:
+            if any(arg is not None for arg in (coefs, vars_1, vars_2)):
+                raise ValueError('Cannot specify both args and any of '
+                                 '{coefs, vars_1, vars_2}')
+            self._args_ = args
+        else:
+            self._coefs = tuple(coefs)
+            self._vars_1 = tuple(vars_1)
+            self._vars_2 = tuple(vars_2)
+            self._args_cache_ = tuple()
+
+        if len(self.vars_1) != self.nargs() or len(self.vars_2) != self.nargs():
+            raise ValueError('The length of coefs, vars_1, and '
+                             'vars_2 must be the same.')
+
+    @property
+    def coefs(self):
+        return self._coefs
+
+    @property
+    def vars_1(self):
+        return self._vars_1
+
+    @property
+    def vars_2(self):
+        return self._vars_2
+
+    def nargs(self):
+        return len(self.coefs)
+
+    @property
+    def _args_(self):
+        nargs = self.nargs()
+        if len(self._args_cache_) != nargs:
+            tmp = map(MonomialTermExpression, zip(self.coefs, self.vars_1))
+            tmp = map(ProductExpression, zip(tmp, self.vars_2))
+            self._args_cache_ = tuple(tmp)
+        return self._args_cache_
+
+    @_args_.setter
+    def _args_(self, val):
+        self._args_cache_ = tuple(val)
+        coefs = list()
+        vars_1 = list()
+        vars_2 = list()
+        for term in val:
+            arg1, arg2 = term.args
+            if (type(arg1) is not MonomialTermExpression
+                    or type(arg2) in native_numeric_types
+                    or not arg2.is_variable_type()):
+                raise ValueError('When constructing a QuadraticExpression with args, '
+                                 'the args must all be ProductExpressions containing a '
+                                 'MonomialTermExpression for the first argument and a '
+                                 'variable for the second argument.')
+            c, v1 = arg1.args
+            v2 = arg2
+            coefs.append(c)
+            vars_1.append(v1)
+            vars_2.append(v2)
+        self._coefs = tuple(coefs)
+        self._vars_1 = tuple(vars_1)
+        self._vars_2 = tuple(vars_2)
+
+    def _precedence(self):
+        return QuadraticExpression.PRECEDENCE
+
+    def create_node_with_local_data(self, args, classtype=None):
+        if classtype is None:
+            if not args:
+                classtype = self.__class__
+            else:
+                correct_format = True
+                for arg in args:
+                    if (type(arg) is not ProductExpression
+                            or type(arg.args[0]) is not MonomialTermExpression
+                            or type(arg.args[1]) in native_numeric_types
+                            or not arg.args[1].is_variable_type()):
+                        correct_format = False
+                        break
+                if correct_format:
+                    classtype = self.__class__
+                else:
+                    classtype = SumExpression
+        return classtype(args)
+
+    def getname(self, *args, **kwds):
+        return 'sum'
+
+    def _compute_polynomial_degree(self, result):
+        res = 0
+        for v1, v2 in zip(self.vars_1, self.vars_2):
+            if v1.fixed:
+                if not v2.fixed:
+                    res = max(res, 1)
+            else:
+                if v2.fixed:
+                    res = max(res, 1)
+                else:
+                    res = max(res, 2)
+        return res
+
+    def is_constant(self):
+        return self.nargs() == 0
+
+    def _is_fixed(self, values=None):
+        return all(v.fixed for v in self.vars_1) and all(v.fixed for v in self.vars_2)
+
+    def is_fixed(self):
+        return self._is_fixed()
+
+    def _to_string(self, values, verbose, smap, compute_values):
+        if not values:
+            values = ['0']
+        if verbose:
+            return "%s(%s)" % (self.getname(), ', '.join(values))
+
+        for i in range(1, len(values)):
+            term = values[i]
+            if term[0] not in '+-':
+                values[i] = '+ ' + term
+            elif term[1] != ' ':
+                values[i] = term[0] + ' ' + term[1:]
+        return ' '.join(values)
+
+    def is_potentially_variable(self):
+        return self.nargs() > 0
+
+    def _apply_operation(self, result):
+        return sum(result)
+
+
 class _MutableLinearExpression(LinearExpression):
     __slots__ = ()
 

pyomo/core/tests/unit/test_numeric_expr.py

@@ -41,7 +41,7 @@
 from pyomo.core.expr.numeric_expr import (
     ExpressionBase, UnaryFunctionExpression, SumExpression, PowExpression,
     ProductExpression, NegationExpression, linear_expression,
-    MonomialTermExpression, LinearExpression, DivisionExpression,
+    MonomialTermExpression, LinearExpression, QuadraticExpression, DivisionExpression,
     NPV_NegationExpression, NPV_ProductExpression, 
     NPV_PowExpression, NPV_DivisionExpression,
     decompose_term, clone_counter, nonlinear_expression,
@@ -57,6 +57,7 @@
 from pyomo.core.expr.template_expr import IndexTemplate
 from pyomo.core.expr import expr_common
 from pyomo.core.base.var import _GeneralVarData
+from pyomo.core.expr.compare import compare_expressions
 
 from pyomo.repn import generate_standard_repn
 from pyomo.core.expr.numvalue import NumericValue
@@ -4784,6 +4785,138 @@ def test_pow_other(self):
             self.assertIs(e.__class__, PowExpression)
 
 
+class TestQuadraticExpression(unittest.TestCase):
+    def test_init_without_args(self):
+        m = ConcreteModel()
+        m.a = Set(initialize=[0,1,2])
+        m.x = Var(m.a)
+        m.y = Var(m.a)
+        coefs = [0.1, 20, 3]
+        e = QuadraticExpression(coefs=coefs, vars_1=m.x.values(), vars_2=m.y.values())
+        self.assertEqual(e._args_cache_, tuple())
+        self.assertEqual(e.coefs, tuple(coefs))
+        self.assertEqual(e.vars_1, (m.x[0], m.x[1], m.x[2]))
+        self.assertEqual(e.vars_2, (m.y[0], m.y[1], m.y[2]))
+        self.assertTrue(compare_expressions(e.args[0], 0.1*m.x[0]*m.y[0]))
+        self.assertTrue(compare_expressions(e.args[1], 20*m.x[1]*m.y[1]))
+        self.assertTrue(compare_expressions(e.args[2], 3*m.x[2]*m.y[2]))
+        self.assertEqual(len(e._args_cache_), 3)
+
+    def test_init_with_args(self):
+        m = ConcreteModel()
+        m.a = Set(initialize=[0,1,2])
+        m.x = Var(m.a)
+        m.y = Var(m.a)
+        coefs = [0.1, 20, 3]
+        e = QuadraticExpression(args=[coefs[i]*m.x[i]*m.y[i] for i in m.a])
+        self.assertEqual(len(e._args_cache_), 3)
+        self.assertEqual(e.coefs, tuple(coefs))
+        self.assertEqual(e.vars_1, (m.x[0], m.x[1], m.x[2]))
+        self.assertEqual(e.vars_2, (m.y[0], m.y[1], m.y[2]))
+        print(e.args[0])
+        self.assertTrue(compare_expressions(e.args[0], 0.1*m.x[0]*m.y[0]))
+        self.assertTrue(compare_expressions(e.args[1], 20*m.x[1]*m.y[1]))
+        self.assertTrue(compare_expressions(e.args[2], 3*m.x[2]*m.y[2]))
+
+    def test_init_errors(self):
+        m = ConcreteModel()
+        m.a = Set(initialize=[0,1,2])
+        m.x = Var(m.a)
+        m.y = Var(m.a)
+        coefs = [0.1, 20, 3]
+        with self.assertRaisesRegex(ValueError, 'Cannot specify both args and any of {coefs, vars_1, vars_2}'):
+            e = QuadraticExpression(args=[coefs[i]*m.x[i]*m.y[i] for i in m.a], coefs=coefs)
+        with self.assertRaisesRegex(ValueError, 'Cannot specify both args and any of {coefs, vars_1, vars_2}'):
+            e = QuadraticExpression(args=[coefs[i]*m.x[i]*m.y[i] for i in m.a], vars_1=m.x.values())
+        with self.assertRaisesRegex(ValueError, 'Cannot specify both args and any of {coefs, vars_1, vars_2}'):
+            e = QuadraticExpression(args=[coefs[i]*m.x[i]*m.y[i] for i in m.a], vars_2=m.x.values())
+        with self.assertRaisesRegex(ValueError, 'The length of coefs, vars_1, and vars_2 must be the same.'):
+            e = QuadraticExpression(coefs=coefs, vars_1=[m.x[0]], vars_2=m.y.values())
+        with self.assertRaisesRegex(ValueError, 'The length of coefs, vars_1, and vars_2 must be the same.'):
+            e = QuadraticExpression(coefs=coefs, vars_1=m.x.values(), vars_2=[m.y[0]])
+        with self.assertRaisesRegex(ValueError, 'The length of coefs, vars_1, and vars_2 must be the same.'):
+            e = QuadraticExpression(coefs=coefs[0:2], vars_1=m.x.values(), vars_2=m.y.values())
+        with self.assertRaisesRegex(ValueError,
+                                    'When constructing a QuadraticExpression with '
+                                    'args, the args must all be ProductExpressions '
+                                    'containing a MonomialTermExpression for the first '
+                                    'argument and a variable for the second argument.'):
+            e = QuadraticExpression(args=[coefs[i]*(m.x[i]*m.y[i]) for i in m.a])
+
+    def test_to_string(self):
+        m = ConcreteModel()
+        m.a = Set(initialize=[0,1,2,3])
+        m.x = Var(m.a)
+        m.y = Var(m.a)
+        coefs = [1, -1, 2, -2]
+        e = QuadraticExpression(coefs=[], vars_1=[], vars_2=[])
+        self.assertEqual(e.to_string(), "0")
+        e = QuadraticExpression(coefs=coefs, vars_1=m.x.values(), vars_2=m.y.values())
+        self.assertEqual(e.to_string(), "x[0]*y[0] - x[1]*y[1] + 2*x[2]*y[2] - 2*x[3]*y[3]")
+
+    def test_polynomial_degree(self):
+        m = ConcreteModel()
+        m.a = Set(initialize=[0,1,2,3])
+        m.x = Var(m.a)
+        m.y = Var(m.a)
+        coefs = [1, -1, 2, -2]
+        e = QuadraticExpression(coefs=coefs, vars_1=m.x.values(), vars_2=m.y.values())
+        self.assertEqual(e.polynomial_degree(), 2)
+        m.x.fix(1)
+        m.y.fix(2)
+        self.assertEqual(e.polynomial_degree(), 0)
+        m.x.unfix()
+        self.assertEqual(e.polynomial_degree(), 1)
+        m.x.fix()
+        m.y.unfix()
+        self.assertEqual(e.polynomial_degree(), 1)
+
+    def test_is_constant_and_fixed(self):
+        m = ConcreteModel()
+        m.a = Set(initialize=[0,1,2,3])
+        m.x = Var(m.a)
+        m.y = Var(m.a)
+        coefs = [1, -1, 2, -2]
+        e = QuadraticExpression(coefs=coefs, vars_1=m.x.values(), vars_2=m.y.values())
+        self.assertFalse(e.is_constant())
+        self.assertFalse(e.is_fixed())
+        self.assertTrue(e.is_potentially_variable())
+        m.x.fix(1)
+        m.y.fix(2)
+        self.assertFalse(e.is_constant())
+        self.assertTrue(e.is_fixed())
+        self.assertTrue(e.is_potentially_variable())
+        m.x.unfix()
+        self.assertFalse(e.is_constant())
+        self.assertFalse(e.is_fixed())
+        self.assertTrue(e.is_potentially_variable())
+        m.x.fix()
+        m.y.unfix()
+        self.assertFalse(e.is_constant())
+        self.assertFalse(e.is_fixed())
+        self.assertTrue(e.is_potentially_variable())
+
+        e = QuadraticExpression(coefs=[], vars_1=[], vars_2=[])
+        self.assertTrue(e.is_constant())
+        self.assertTrue(e.is_fixed())
+        self.assertFalse(e.is_potentially_variable())
+
+    def test_apply_operation(self):
+        m = ConcreteModel()
+        m.a = Set(initialize=[0,1,2])
+        m.x = Var(m.a)
+        m.y = Var(m.a)
+        coefs = [0.1, 20, 3]
+        e = QuadraticExpression(coefs=coefs, vars_1=m.x.values(), vars_2=m.y.values())
+        m.x[0].value = 1.2
+        m.x[1].value = 2.3
+        m.x[2].value = 3.4
+        m.y[0].value = -1.5
+        m.y[1].value = 2.6
+        m.y[2].value = 3.8
+        self.assertAlmostEqual(value(e), sum(coefs[i]*m.x[i].value*m.y[i].value for i in m.a))
+
+
 class TestNonlinearExpression(unittest.TestCase):
 
     def test_sum_other(self):