Merge pull request #448 from qutech/issues/447_render_error_float_pre…

…cision

Fix #447

ReleaseNotes.txt

@@ -7,6 +7,8 @@
 - Pulse Templates:
     - `MappingPulseTemplate`:
         - Raise a ValueError if more than one inner channel is mapped to the same outer channel
+    - Plotting:
+        - Make `plotting.render` behaviour and return value consistent between calls with `InstructionBlock` and `Loop`. Render now always returns 3 arguments.
 
 ## 0.3 ##
 

qupulse/pulses/plotting.py

@@ -14,10 +14,10 @@
 import operator
 import itertools
 
-from qupulse.utils.types import ChannelID, MeasurementWindow, has_type_interface
+from qupulse.utils.types import ChannelID, MeasurementWindow, has_type_interface, TimeType
 from qupulse.pulses.pulse_template import PulseTemplate
 from qupulse.pulses.parameters import Parameter
-from qupulse._program.waveforms import Waveform
+from qupulse._program.waveforms import Waveform, SequenceWaveform
 from qupulse._program.instructions import EXECInstruction, STOPInstruction, AbstractInstructionBlock, \
     REPJInstruction, MEASInstruction, GOTOInstruction, InstructionPointer
 from qupulse._program._loop import Loop, to_waveform
@@ -50,7 +50,7 @@ def iter_waveforms(instruction_block: AbstractInstructionBlock,
 
 
 def iter_instruction_block(instruction_block: AbstractInstructionBlock,
-                           extract_measurements: bool) -> Tuple[list, list, Real]:
+                           extract_measurements: bool) -> Tuple[list, list, TimeType]:
     """Iterates over the instructions contained in an InstructionBlock (thus simulating execution).
 
     In effect, this function simulates the execution of the control flow represented by the passed InstructionBlock
@@ -72,7 +72,7 @@ def iter_instruction_block(instruction_block: AbstractInstructionBlock,
     block_stack = [(enumerate(instruction_block), None)]
     waveforms = []
     measurements = []
-    time = 0
+    time = TimeType(0)
 
     while block_stack:
         block, expected_return = block_stack.pop()
@@ -107,11 +107,11 @@ def iter_instruction_block(instruction_block: AbstractInstructionBlock,
 
 
 def render(program: Union[AbstractInstructionBlock, Loop],
-           sample_rate: Real=10.0,
-           render_measurements: bool=False,
-           time_slice: Tuple[Real, Real]=None) -> Union[Tuple[np.ndarray, Dict[ChannelID, np.ndarray]],
-                                                        Tuple[np.ndarray, Dict[ChannelID, np.ndarray],
-                                                              List[MeasurementWindow]]]:
+           sample_rate: Real = 10.0,
+           render_measurements: bool = False,
+           time_slice: Tuple[Real, Real] = None,
+           plot_channels: Optional[Set[ChannelID]] = None) -> Tuple[np.ndarray, Dict[ChannelID, np.ndarray],
+                                                                    List[MeasurementWindow]]:
     """'Renders' a pulse program.
 
         Samples all contained waveforms into an array according to the control flow of the program.
@@ -121,7 +121,8 @@ def render(program: Union[AbstractInstructionBlock, Loop],
                 old-fashioned InstructionBlock.
             sample_rate: The sample rate in GHz.
             render_measurements: If True, the third return value is a list of measurement windows.
-            time_slice: The time slice to be plotted. If None, the entire pulse will be shown.
+            time_slice: The time slice to be rendered. If None, the entire pulse will be shown.
+            plot_channels: Only channels in this set are rendered. If None, all will.
 
         Returns:
             A tuple (times, values, measurements). times is a numpy.ndarray of dimensions sample_count where
@@ -131,92 +132,84 @@ def render(program: Union[AbstractInstructionBlock, Loop],
             (name, start_time, duration).
         """
     if has_type_interface(program, Loop):
-        return _render_loop(program, sample_rate=sample_rate,
-                            render_measurements=render_measurements, time_slice=time_slice)
+        waveform, measurements = _render_loop(program, render_measurements=render_measurements)
     elif has_type_interface(program, AbstractInstructionBlock):
         warnings.warn("InstructionBlock API is deprecated", DeprecationWarning)
         if time_slice is not None:
             raise ValueError("Keyword argument time_slice is not supported when rendering instruction blocks")
-        return _render_instruction_block(program, sample_rate=sample_rate, render_measurements=render_measurements)
+        waveform, measurements = _render_instruction_block(program, render_measurements=render_measurements)
     else:
         raise ValueError('Cannot render an object of type %r' % type(program), program)
 
+    if waveform is None:
+        return np.array([]), dict(), measurements
 
-def _render_instruction_block(sequence: AbstractInstructionBlock,
-                              sample_rate: Real=10.0,
-                              render_measurements=False) -> Union[Tuple[np.ndarray, Dict[ChannelID, np.ndarray]],
-                                                                  Tuple[np.ndarray, Dict[ChannelID, np.ndarray],
-                                                                        List[MeasurementWindow]]]:
-    """The specific implementation of render for InstructionBlock arguments."""
+    if plot_channels is None:
+        channels = waveform.defined_channels
+    else:
+        channels = waveform.defined_channels & plot_channels
 
-    waveforms, measurements, total_time = iter_instruction_block(sequence, render_measurements)
-    if not waveforms:
-        return np.empty(0), dict()
+    if time_slice is None:
+        start_time, end_time = 0, waveform.duration
+
+    elif time_slice[1] < time_slice[0] or time_slice[0] < 0 or time_slice[1] < 0:
+        raise ValueError("time_slice is not valid.")
+
+    else:
+        start_time, end_time, *_ = time_slice
 
-    channels = waveforms[0].defined_channels
+        # filter measurement windows
+        measurements = [(name, begin, length)
+                        for name, begin, length in measurements
+                        if begin < end_time and begin + length > start_time]
 
-    # add one sample to see the end of the waveform
-    sample_count = total_time * sample_rate + 1
+    sample_count = (end_time - start_time) * sample_rate + 1
+    if sample_count < 2:
+        raise PlottingNotPossibleException(pulse=None,
+                                           description='cannot render sequence with less than 2 data points')
     if not float(sample_count).is_integer():
-        warnings.warn('Sample count not whole number. Casted to integer.')
-    times = np.linspace(0, float(total_time), num=int(sample_count), dtype=float)
-    # move the last sample inside the waveform
+        warnings.warn("Sample count is not an integer. Will be rounded (this changes the sample rate).")
+
+    times = np.linspace(float(start_time), float(end_time), num=int(sample_count), dtype=float)
     times[-1] = np.nextafter(times[-1], times[-2])
 
-    voltages = dict((ch, np.empty(len(times))) for ch in channels)
-    offset = 0
-    for waveform in waveforms:
-        wf_end = offset + waveform.duration
-        indices = slice(*np.searchsorted(times, (offset, wf_end)))
-        sample_times = times[indices] - float(offset)
-        for channel in channels:
-            output_array = voltages[channel][indices]
-            waveform.get_sampled(channel=channel,
-                                 sample_times=sample_times,
-                                 output_array=output_array)
-        assert(output_array.shape == sample_times.shape)
-        offset = wf_end
+    voltages = {ch: np.empty_like(times)
+                for ch in channels}
+    for ch, ch_voltage in voltages.items():
+        waveform.get_sampled(channel=ch, sample_times=times, output_array=ch_voltage)
+
     return times, voltages, measurements
 
 
-def _render_loop(loop: Loop,
-                 sample_rate: Real,
-                 render_measurements: bool,
-                 time_slice: Tuple[Real, Real] = None) -> Union[Tuple[np.ndarray, Dict[ChannelID, np.ndarray]],
-                                                                Tuple[np.ndarray, Dict[ChannelID, np.ndarray],
-                                                                List[MeasurementWindow]]]:
-    """The specific implementation of render for Loop arguments."""
-    waveform = to_waveform(loop)
-    channels = waveform.defined_channels
+def _render_instruction_block(sequence: AbstractInstructionBlock,
+                              render_measurements=False) -> Tuple[Optional[Waveform], List[MeasurementWindow]]:
+    """Transform program into single waveform and measurement windows. The specific implementation of render for
+    InstructionBlock arguments."""
 
-    if time_slice is None:
-        time_slice = (0, waveform.duration)
-    elif time_slice[1] < time_slice[0] or time_slice[0] < 0 or time_slice[1] < 0:
-        raise ValueError("time_slice is not valid.")
+    waveforms, measurements, total_time = iter_instruction_block(sequence, render_measurements)
+    if not waveforms:
+        return None, measurements
+    sequence_waveform = SequenceWaveform(waveforms)
 
-    sample_count = (time_slice[1] - time_slice[0]) * sample_rate + 1
-    if sample_count<2:
-        raise PlottingNotPossibleException(pulse = None, description = 'cannot render sequence with less than 2 data points')
-    times = np.linspace(float(time_slice[0]), float(time_slice[1]), num=int(sample_count), dtype=float)
-    times[-1] = np.nextafter(times[-1], times[-2])
+    return sequence_waveform, measurements
 
-    voltages = {}
 
-    for ch in channels:
-        voltages[ch] = waveform.get_sampled(ch, times)
+def _render_loop(loop: Loop,
+                 render_measurements: bool,) -> Tuple[Waveform, List[MeasurementWindow]]:
+    """Transform program into single waveform and measurement windows.
+    The specific implementation of render for Loop arguments."""
+    waveform = to_waveform(loop)
 
     if render_measurements:
         measurement_dict = loop.get_measurement_windows()
         measurement_list = []
         for name, (begins, lengths) in measurement_dict.items():
-            measurement_list.extend(m
-                                    for m in zip(itertools.repeat(name), begins, lengths)
-                                    if m[1]+m[2] > time_slice[0] and m[1] < time_slice[1])
+            measurement_list.extend(zip(itertools.repeat(name), begins, lengths))
         measurements = sorted(measurement_list, key=operator.itemgetter(1))
     else:
         measurements = []
 
-    return times, voltages, measurements
+    return waveform, measurements
 
 
 def plot(pulse: PulseTemplate,
@@ -271,7 +264,7 @@ def plot(pulse: PulseTemplate,
     if program is not None:
         times, voltages, measurements = render(program,
                                                sample_rate,
-                                               render_measurements=plot_measurements,
+                                               render_measurements=bool(plot_measurements),
                                                time_slice=time_slice)
     else:
         times, voltages, measurements = np.array([]), dict(), []

tests/pulses/plotting_tests.py

@@ -25,8 +25,9 @@ def test_render_unsupported_instructions(self) -> None:
             render(block)
 
     def test_render_no_waveforms(self) -> None:
-        time, channel_data = render(InstructionBlock())
+        time, channel_data, measurements = render(InstructionBlock())
         self.assertEqual(channel_data, dict())
+        self.assertEqual(measurements, [])
         numpy.testing.assert_equal(time, numpy.empty(0))
 
     def test_iter_waveforms(self) -> None:
@@ -300,6 +301,29 @@ def test_plot_empty_pulse(self) -> None:
         with self.assertWarnsRegex(UserWarning, "empty", msg="plot() did not issue a warning for an empty pulse"):
             plot(pt, dict(), show=False)
 
+    def test_bug_447(self):
+        """Code from https://github.com/qutech/qupulse/issues/447"""
+        TablePT = TablePulseTemplate
+        SequencePT = SequencePulseTemplate
+        Sequencing = Sequencer
+
+        period = 8.192004194306148e-05
+        repetitions = 80
+        sampling_rate = 1e7
+        sec_to_ns = 1e9
+
+        table_pt = TablePT({'test': [(0, 0), (period * sec_to_ns, 0, 'linear')]})
+
+        sequencer = Sequencing()
+        template = SequencePT(*((table_pt,) * repetitions))
+        channels = template.defined_channels
+        sequencer.push(template, dict(), channel_mapping={ch: ch for ch in channels},
+                       window_mapping={w: w for w in template.measurement_names})
+        instructions = sequencer.build()
+
+        with self.assertWarns(UserWarning):
+            (_, voltages, _) = render(instructions, sampling_rate / sec_to_ns)
+
 
 class PlottingNotPossibleExceptionTests(unittest.TestCase):