Issues/423 ppt integral bug

ReleaseNotes.txt

@@ -11,6 +11,8 @@
         - Make duration equality check approximate (numeric tolerance)
     - Plotting:
         - Add `time_slice` keyword argument to render() and plot()
+    - `PointPulseTemplate`:
+        - Fixed bug in integral evaluation
 
 - Expressions:
     - Make ExpressionScalar hashable

qupulse/pulses/interpolation.py

@@ -79,7 +79,7 @@ def __call__(self,
 
     @property
     def integral(self) -> ExpressionScalar:
-        return ExpressionScalar('0.5 * (t1-t0) * (v0 + v1)')
+        return ExpressionScalar('(t1-t0) * (v0 + v1) / 2')
 
     @property
     def expression(self) -> ExpressionScalar:

qupulse/pulses/point_pulse_template.py

@@ -3,8 +3,10 @@
 import itertools
 import numbers
 
+import sympy
 import numpy as np
 
+from qupulse.utils.sympy import Broadcast
 from qupulse.utils.types import ChannelID
 from qupulse.expressions import Expression, ExpressionScalar
 from qupulse.pulses.conditions import Condition
@@ -148,12 +150,16 @@ def requires_stop(self,
     def integral(self) -> Dict[ChannelID, ExpressionScalar]:
         expressions = {channel: 0 for channel in self._channels}
         for first_entry, second_entry in zip(self._entries[:-1], self._entries[1:]):
-            substitutions = {'t0': ExpressionScalar(first_entry.t).sympified_expression,
-                             't1': ExpressionScalar(second_entry.t).sympified_expression}
+            substitutions = {'t0': first_entry.t.sympified_expression,
+                             't1': second_entry.t.sympified_expression}
+
+            v0 = sympy.IndexedBase(Broadcast(first_entry.v.underlying_expression, (len(self.defined_channels),)))
+            v1 = sympy.IndexedBase(Broadcast(second_entry.v.underlying_expression, (len(self.defined_channels),)))
 
             for i, channel in enumerate(self._channels):
-                substitutions['v0'] = ExpressionScalar(first_entry.v[i]).sympified_expression
-                substitutions['v1'] = ExpressionScalar(second_entry.v[i]).sympified_expression
+                substitutions['v0'] = v0[i]
+                substitutions['v1'] = v1[i]
+
                 expressions[channel] += first_entry.interp.integral.sympified_expression.subs(substitutions)
 
         expressions = {c: ExpressionScalar(expressions[c]) for c in expressions}

qupulse/utils/sympy.py

@@ -65,6 +65,25 @@ def __contains__(self, k) -> bool:
         return True
 
 
+class Broadcast(sympy.Function):
+    """Broadcast x to the specified shape using numpy.broadcast_to
+
+    Examples:
+        >>> bc = Broadcast('a', (3,))
+        >>> assert bc.subs({'a': 2}) == sympy.Array([2, 2, 2])
+        >>> assert bc.subs({'a': (1, 2, 3)}) == sympy.Array([1, 2, 3])
+    """
+
+    @classmethod
+    def eval(cls, x, shape) -> Optional[sympy.Array]:
+        if hasattr(shape, 'free_symbols') and shape.free_symbols:
+            # cannot do anything
+            return None
+
+        if hasattr(x, '__len__') or not x.free_symbols:
+            return sympy.Array(numpy.broadcast_to(x, shape))
+
+
 class Len(sympy.Function):
     nargs = 1
 
@@ -78,7 +97,8 @@ def eval(cls, arg) -> Optional[sympy.Integer]:
 
 
 sympify_namespace = {'len': Len,
-                     'Len': Len}
+                     'Len': Len,
+                     'Broadcast': Broadcast}
 
 
 def numpy_compatible_mul(*args) -> Union[sympy.Mul, sympy.Array]:
@@ -199,7 +219,7 @@ def recursive_substitution(expression: sympy.Expr,
 _numpy_environment = {**_base_environment, **numpy.__dict__}
 _sympy_environment = {**_base_environment, **sympy.__dict__}
 
-_lambdify_modules = [{'ceiling': numpy_compatible_ceiling}, 'numpy', _special_functions]
+_lambdify_modules = [{'ceiling': numpy_compatible_ceiling, 'Broadcast': numpy.broadcast_to}, 'numpy', _special_functions]
 
 
 def evaluate_compiled(expression: sympy.Expr,

tests/pulses/point_pulse_template_tests.py

@@ -76,21 +76,27 @@ def test_parameter_names(self):
 
     def test_integral(self) -> None:
         pulse = PointPulseTemplate(
-            [(1, (2, 'b'), 'linear'), (3, (0, 0), 'jump'), (4, (2, 'c'), 'hold'), (5, (8, 'd'), 'hold')],
+            [(1, (2, 'b'), 'linear'),
+             (3, (0, 0), 'jump'),
+             (4, (2, 'c'), 'hold'),
+             (5, (8, 'd'), 'hold')],
             [0, 'other_channel']
         )
-        self.assertEqual({0: ExpressionScalar(6),
-                          'other_channel': ExpressionScalar('1.0*b + 2.0*c')},
+        self.assertEqual({0: ExpressionScalar('6'),
+                          'other_channel': ExpressionScalar('b + 2*c')},
                          pulse.integral)
 
         pulse = PointPulseTemplate(
-            [(1, ('2', 'b'), 'linear'), ('t0', (0, 0), 'jump'), (4, (2, 'c'), 'hold'), ('g', (8, 'd'), 'hold')],
+            [(1, ('2', 'b'), 'linear'), ('t0', (0, 0), 'jump'), (4, (2.0, 'c'), 'hold'), ('g', (8, 'd'), 'hold')],
             ['symbolic', 1]
         )
-        self.assertEqual({'symbolic': ExpressionScalar('2.0*g - t0 - 1.0'),
-                          1: ExpressionScalar('b*(0.5*t0 - 0.5) + c*(g - 4.0) + c*(-t0 + 4.0)')},
+        self.assertEqual({'symbolic': ExpressionScalar('2.0*g - 1.0*t0 - 1.0'),
+                          1: ExpressionScalar('b*(t0 - 1) / 2 + c*(g - 4) + c*(-t0 + 4)')},
                          pulse.integral)
 
+        ppt = PointPulseTemplate([(0, 0), ('t_init', 0)], ['X', 'Y'])
+        self.assertEqual(ppt.integral, {'X': 0, 'Y': 0})
+
 
 class PointPulseTemplateSequencingTests(unittest.TestCase):
 

tests/pulses/table_pulse_template_tests.py

@@ -426,7 +426,7 @@ def test_integral(self) -> None:
                                             'symbolic': [(3, 'a', 'hold'), ('b', 4, 'linear'), ('c', Expression('d'), 'hold')]})
         self.assertEqual(pulse.integral, {0: Expression('6'),
                                           'other_channel': Expression(7),
-                                          'symbolic': Expression('(b-3)*a + 0.5 * (c-b)*(d+4)')})
+                                          'symbolic': Expression('(b-3)*a + (c-b)*(d+4) / 2')})
 
 
 class TablePulseTemplateConstraintTest(ParameterConstrainerTest):

tests/utils/sympy_tests.py

@@ -16,7 +16,7 @@
 
 from qupulse.utils.sympy import sympify as qc_sympify, substitute_with_eval, recursive_substitution, Len,\
     evaluate_lambdified, evaluate_compiled, get_most_simple_representation, get_variables, get_free_symbols,\
-    almost_equal
+    almost_equal, Broadcast
 
 
 ################################################### SUBSTITUTION #######################################################
@@ -329,3 +329,86 @@ def test_almost_equal(self):
         self.assertFalse(almost_equal(sympy.sin(a), sympy.sin(a) + 1e-14))
 
         self.assertTrue(almost_equal(sympy.sin(a), sympy.sin(a) + 1e-14, epsilon=1e-13))
+
+
+class BroadcastTests(unittest.TestCase):
+    def test_symbolic_shape(self):
+        symbolic = Broadcast(a, (b,))
+        self.assertIs(symbolic.func, Broadcast)
+        self.assertEqual(symbolic.args, (a, (b,)))
+
+        subs_b = symbolic.subs({b: 6})
+        self.assertIs(subs_b.func, Broadcast)
+        self.assertEqual(subs_b.args, (a, (6,)))
+
+        subs_a = symbolic.subs({a: 3})
+        self.assertIs(subs_a.func, Broadcast)
+        self.assertEqual(subs_a.args, (3, (b,)))
+
+        subs_both_scalar = symbolic.subs({a: 3, b: 6})
+        self.assertEqual(subs_both_scalar, sympy.Array([3, 3, 3, 3, 3, 3]))
+
+        subs_both_array = symbolic.subs({a: (1, 2, 3, 4, 5, 6), b: 6})
+        self.assertEqual(subs_both_array, sympy.Array([1, 2, 3, 4, 5, 6]))
+
+        with self.assertRaises(ValueError):
+            symbolic.subs({a: (1, 2, 3, 4, 5, 6), b: 7})
+
+    def test_scalar_broad_cast(self):
+        symbolic = Broadcast(a, (6,))
+        self.assertIs(symbolic.func, Broadcast)
+        self.assertEqual(symbolic.args, (a, (6,)))
+
+        subs_symbol = symbolic.subs({a: b})
+        self.assertIs(subs_symbol.func, Broadcast)
+        self.assertEqual(subs_symbol.args, (b, (6,)))
+
+        subs_scalar = symbolic.subs({a: 3.4})
+        self.assertEqual(subs_scalar, sympy.Array([3.4, 3.4, 3.4, 3.4, 3.4, 3.4]))
+
+        subs_symbol_vector = symbolic.subs({a: (b, 1, 2, 3, 4, 5)})
+        self.assertEqual(subs_symbol_vector, sympy.Array([b, 1, 2, 3, 4, 5]))
+
+        subs_numeric_vector = symbolic.subs({a: (0, 1, 2, 3, 4, 5)})
+        self.assertEqual(subs_numeric_vector, sympy.Array([0, 1, 2, 3, 4, 5]))
+
+        with self.assertRaises(ValueError):
+            symbolic.subs({a: (b, 4, 5)})
+
+        with self.assertRaises(ValueError):
+            symbolic.subs({a: (8, 5, 3, 5, 5, 4, 4, 5)})
+
+    def test_array_broadcast(self):
+        expected = sympy.Array([1, 2, a, b])
+
+        self.assertEqual(expected, Broadcast(list(expected), (4,)))
+        self.assertEqual(expected, Broadcast(tuple(expected), (4,)))
+        self.assertEqual(expected, Broadcast(expected, (4,)))
+
+    def test_numeric_evaluation(self):
+        symbolic = Broadcast(a, (b,))
+
+        arguments = {'a': (1, 2., 3), 'b': 3}
+        expected = np.asarray([1, 2., 3])
+        result, _ = evaluate_lambdified(symbolic, ['a', 'b'], arguments, None)
+        np.testing.assert_array_equal(expected, result)
+
+        with self.assertRaises(ValueError):
+            arguments = {'a': (1, 2., 3), 'b': 4}
+            evaluate_lambdified(symbolic, ['a', 'b'], arguments, None)
+
+        arguments = {'a': 1, 'b': 3}
+        expected = np.asarray([1, 1, 1])
+        result, _ = evaluate_lambdified(symbolic, ['a', 'b'], arguments, None)
+        np.testing.assert_array_equal(expected, result)
+
+    def test_sympification(self):
+        symbolic = Broadcast(a, (3,))
+        as_str = str(symbolic)
+
+        re_sympified = qc_sympify(as_str)
+        self.assertEqual(re_sympified, symbolic)
+
+        sympification = qc_sympify('Broadcast(a, (3,))')
+        self.assertEqual(sympification, symbolic)
+