Repository: qutech/qupulse (language: python)

Issue:
Move drivers into separate projects
The goal is to make driver tinkering independent of each other and indepentend of the frontend.

Currently the driver code of each device lives in two places. A file `_program` contains the programming logic and a file in `hardware/awgs` contains the device communication.

The `_program` package and its contents are currently not considered a stable interface. We need to declare (parts of) `Loop` as stable for a somewhat independent driver implementation to make sense.

 - [ ] Create `qupulse.program.loop` (Simon)
 - [ ] Create the project `qupulse-hdawg-idx` for the `playWaveIndexed` based driver and move `seqc.py` and `hdawg.py` there
 - [ ] Create the project `qupulse-hdawg-cmd` for the command table based driver

Open questions:
 - Where should the projects be located? Gitlab? Here?
 - Should we avoid code duplication between the two projects? Here the plan is to duplicate everything to avoid conflicts and make workarounds for different firmware versions easier.

@Nomos11 @maxbeer99 What are your thoughts on this plan?


Discussion / hints:
Regarding loop: I know too little about that to judge which parts can be deemed "stable" or not.
Regarding the questions:
The first version lives in Gitlab already and could possibly stay there for internal consumption until some issues have been discussed (?) (referring e.g. to reocurring upload issues out of our control and possible optimizations)
I would not have a strong aversion to code duplication; but workarounds have previously been made in separate branches for individual experimental setups anyway...

Fix patch:
diff --git a/changes.d/779.feature b/changes.d/779.feature
new file mode 100644
index 00000000..edae2b25
--- /dev/null
+++ b/changes.d/779.feature
@@ -0,0 +1,1 @@
+Promote parts of the private subpackage `qupulse._program` to the public subpackage `qupulse.program`, i.e. `loop`, `volatile`, `transformation` and `waveforms`. This allows external packages/drivers to rely on stability of the `Loop` class.
diff --git a/doc/source/concepts/program.rst b/doc/source/concepts/program.rst
index 0b34df0a..3d4172bc 100644
--- a/doc/source/concepts/program.rst
+++ b/doc/source/concepts/program.rst
@@ -4,14 +4,19 @@ Instantiated Pulse: Program
 ---------------------------
 
 In qupulse an instantiated pulse template is called a program as it is something that an arbitrary waveform generator
-(AWG) can execute/playback. It is created by the `create_program` method of the pulse template which returns a hardware
-independent representation which is currently of type ``Loop``. Opposed to the `PulseTemplate` interfaces the interface of the program is currently not covered by the qupulse backward compatibility and stability guarantee.
-This is reflected by the fact that it lives in the private module ``qupulse._program._loop``.
+(AWG) can execute/playback.
+It is created by the ``create_program`` method of the pulse template which returns a hardware
+independent representation which is of type ``Loop``.
+This ``Loop`` object is the root node of a tree ``Loop``s of arbitrary depth.
+Each node consists of a repetition count and either a waveform or a sequence of nodes which are repeated that many times.
+Iterations like the ```ForLoopPT`` cannot be represented natively but are unrolled into a sequence of items.
+The repetition count is currently the only property of a program that can be defined as volatile. This means that the AWG driver tries to upload the program in a way, where the repetition count can quickly be changed. This is implemented via the ```VolatileRepetitionCount`` class.
 
-There is no description of the details of the program object here to avoid outdated documentation.
+There is no description of the details of the program object here to avoid duplicated and outdated documentation.
 The documentation is in the docstrings of the source code.
 The program can be thought of as compact representation of a mapping :math:`\{t | 0 \le t \le t_{\texttt{duration}}} \rightarrow \mathbb{R}^n` from the time while the program lasts :math:´t´ to an n-dimensional voltage space :math:´\mathbb{R}^n´.
 The dimensions are named by the channel names.
 
-The subpackage ``_qupulse._program`` also contains hardware specific translations of the programs for example a
-transpiler to Zurich Instruments sequencing C in ``_qupulse._program.seqc``.
+The ``Loop`` class and its constituents ``Waveform`` and ``VolatileRepetitionCount`` are defined in the ``qupulse.program`` subpackage and it's submodules.
+The private subpackage ``qupulse._program`` contains AWG driver internals that can change with any release, for example a
+transpiler to Zurich Instruments sequencing C in ``qupulse._program.seqc``.
diff --git a/qupulse/_program/_loop.py b/qupulse/_program/_loop.py
index 20ec1204..b5ad53de 100644
--- a/qupulse/_program/_loop.py
+++ b/qupulse/_program/_loop.py
@@ -1,724 +1,9 @@
-from typing import Union, Dict, Iterable, Tuple, cast, List, Optional, Generator, Mapping
-from collections import defaultdict
-from enum import Enum
-import warnings
-import bisect
-import reprlib
+"""Backwards compatibility link to qupulse.program.loop"""
 
-import numpy as np
-import sympy.ntheory
+from qupulse.program.loop import *
 
-from qupulse._program.waveforms import Waveform, ConstantWaveform
-from qupulse._program.volatile import VolatileRepetitionCount, VolatileProperty
+import qupulse.program.loop
 
-from qupulse.utils import is_integer
-from qupulse.utils.types import TimeType, MeasurementWindow
-from qupulse.utils.tree import Node, is_tree_circular
-from qupulse.utils.numeric import smallest_factor_ge
+__all__ = qupulse.program.loop.__all__
 
-from qupulse._program.waveforms import SequenceWaveform, RepetitionWaveform
-
-__all__ = ['Loop', 'make_compatible', 'MakeCompatibleWarning']
-
-
-class Loop(Node):
-    MAX_REPR_SIZE = 2000
-    __slots__ = ('_waveform', '_measurements', '_repetition_definition', '_cached_body_duration')
-
-    """This class represents a initialized (sub-)program as a tree. Each Loop of a valid program has a repetition count
-    and either a waveform or a sequence of loops as children.
-    
-    A Loop can have associated measurements which are also repeated.
-    """
-    def __init__(self,
-                 parent: Union['Loop', None] = None,
-                 children: Iterable['Loop'] = (),
-                 waveform: Optional[Waveform] = None,
-                 measurements: Optional[List[MeasurementWindow]] = None,
-                 repetition_count: Union[int, VolatileRepetitionCount] = 1):
-        """Initialize a new loop
-
-        Args:
-            parent: Forwarded to Node.__init__
-            children: Forwarded to Node.__init__
-            waveform: "Payload"
-            measurements: Associated measurements
-            repetition_count: The children / waveform are repeated this often
-        """
-        super().__init__(parent=parent, children=children)
-
-        self._waveform = waveform
-        self._measurements = measurements
-        self._repetition_definition = repetition_count
-        self._cached_body_duration = None
-        assert isinstance(repetition_count, VolatileRepetitionCount) or is_integer(repetition_count)
-        assert isinstance(waveform, (type(None), Waveform))
-
-    def __eq__(self, other: 'Loop') -> bool:
-        if type(self) == type(other):
-            return (self._repetition_definition == other._repetition_definition and
-                    self.waveform == other.waveform and
-                    (self._measurements or None) == (other._measurements or None) and
-                    len(self) == len(other) and
-                    all(self_child == other_child for self_child, other_child in zip(self, other)))
-        else:
-            return NotImplemented
-
-    def append_child(self, loop: Optional['Loop'] = None, **kwargs) -> None:
-        """Append a child to this loop. Either an existing Loop object or a newly created from kwargs
-
-        Args:
-            loop: loop to append
-            **kwargs: Child is constructed with these kwargs
-
-        Raises:
-            ValueError: if called with loop and kwargs
-        """
-        if loop is not None:
-            if kwargs:
-                raise ValueError("Cannot pass a Loop object and Loop constructor arguments at the same time in "
-                                 "append_child")
-            arg = (loop,)
-        else:
-            arg = (kwargs,)
-        super().__setitem__(slice(len(self), len(self)), arg)
-        self._invalidate_duration(body_duration_increment=self[-1].duration)
-
-    def _invalidate_duration(self, body_duration_increment=None):
-        if self._cached_body_duration is not None:
-            if body_duration_increment is not None:
-                self._cached_body_duration += body_duration_increment
-            else:
-                self._cached_body_duration = None
-        if self.parent:
-            if body_duration_increment is not None:
-                self.parent._invalidate_duration(body_duration_increment=body_duration_increment*self.repetition_count)
-            else:
-                self.parent._invalidate_duration()
-
-    def add_measurements(self, measurements: Iterable[MeasurementWindow]):
-        """Add measurements offset by the current body duration i.e. to the END of the current loop
-
-        Args:
-            measurements: Measurements to add
-        """
-        body_duration = float(self.body_duration)
-        if body_duration == 0:
-            measurements = measurements
-        else:
-            measurements = ((mw_name, begin+body_duration, length) for mw_name, begin, length in measurements)
-
-        if self._measurements is None:
-            self._measurements = list(measurements)
-        else:
-            self._measurements.extend(measurements)
-
-    @property
-    def waveform(self) -> Waveform:
-        return self._waveform
-
-    @waveform.setter
-    def waveform(self, val) -> None:
-        self._waveform = val
-        self._invalidate_duration()
-
-    @property
-    def body_duration(self) -> TimeType:
-        if self._cached_body_duration is None:
-            if self.is_leaf():
-                if self.waveform:
-                    self._cached_body_duration = self.waveform.duration
-                else:
-                    self._cached_body_duration = TimeType.from_fraction(0, 1)
-            else:
-                self._cached_body_duration = sum(child.duration for child in self)
-        return self._cached_body_duration
-
-    @property
-    def duration(self) -> TimeType:
-        return self.body_duration * TimeType.from_fraction(self.repetition_count, 1)
-
-    @property
-    def volatile_repetition(self) -> Optional[VolatileProperty]:
-        return getattr(self._repetition_definition, 'volatile_property', None)
-
-    @property
-    def repetition_definition(self) -> Union[int, VolatileRepetitionCount]:
-        return self._repetition_definition
-
-    @repetition_definition.setter
-    def repetition_definition(self, new_definition: Union[int, VolatileRepetitionCount]):
-        self._repetition_definition = new_definition
-
-    @property
-    def repetition_count(self) -> int:
-        return int(self._repetition_definition)
-
-    @repetition_count.setter
-    def repetition_count(self, val: int) -> None:
-        assert isinstance(val, (int, float))
-        new_repetition = int(val)
-        if abs(new_repetition - val) > 1e-10:
-            raise ValueError('Repetition count was not an integer')
-        self._repetition_definition = new_repetition
-
-    def unroll(self) -> None:
-        if self.is_leaf():
-            raise RuntimeError('Leaves cannot be unrolled')
-        if self.volatile_repetition:
-            warnings.warn("Unrolling a Loop with volatile repetition count", VolatileModificationWarning)
-
-        i = self.parent_index
-        self.parent[i:i+1] = (child.copy_tree_structure(new_parent=self.parent)
-                              for _ in range(self.repetition_count)
-                              for child in self)
-        self.parent.assert_tree_integrity()
-
-    def __setitem__(self, idx, value):
-        super().__setitem__(idx, value)
-        self._invalidate_duration()
-
-    def unroll_children(self) -> None:
-        if self.volatile_repetition:
-            warnings.warn("Unrolling a Loop with volatile repetition count", VolatileModificationWarning)
-        old_children = self.children
-        self[:] = (child.copy_tree_structure()
-                   for _ in range(self.repetition_count)
-                   for child in old_children)
-        self.repetition_count = 1
-        self.assert_tree_integrity()
-
-    def encapsulate(self) -> None:
-        """Add a nesting level by moving self to its children."""
-        self[:] = [Loop(children=self,
-                        repetition_count=self._repetition_definition,
-                        waveform=self._waveform,
-                        measurements=self._measurements)]
-        self.repetition_count = 1
-        self._waveform = None
-        self._measurements = None
-        self.assert_tree_integrity()
-
-    def _get_str(self, first_prefix, other_prefixes) -> Generator[str, None, None]:
-        if self.is_leaf():
-            yield '%sEXEC %r %d times' % (first_prefix, self._waveform, self.repetition_count)
-        else:
-            yield '%sLOOP %d times:' % (first_prefix, self.repetition_count)
-
-            for elem in self:
-                yield from cast(Loop, elem)._get_str(other_prefixes + '  ->', other_prefixes + '    ')
-
-    @reprlib.recursive_repr()
-    def __str__(self) -> str:
-        str_len = 0
-        repr_list = []
-        for sub_repr in self._get_str('', ''):
-            str_len += len(sub_repr)
-
-            if self.MAX_REPR_SIZE and str_len > self.MAX_REPR_SIZE:
-                repr_list.append('...')
-                break
-            else:
-                repr_list.append(sub_repr)
-        return '\n'.join(repr_list)
-
-    @reprlib.recursive_repr()
-    def __repr__(self):
-        children = self.children
-        waveform = self.waveform
-        measurements = self._measurements
-        repetition_count = self.repetition_count
-
-        max_repr_size = self.MAX_REPR_SIZE
-
-        reprs = {}
-        if children:
-            if len(children) < max_repr_size // len('Loop(...)'):
-                reprs['children'] = repr(list(children))
-            else:
-                reprs['children'] = '[...]'
-
-        if waveform:
-            waveform_repr = repr(waveform)
-            if len(waveform_repr) >= max_repr_size:
-                waveform_repr = '...'
-            reprs['waveform'] = waveform_repr
-
-        if measurements:
-            meas_repr = repr(measurements)
-            if len(meas_repr) >= max_repr_size:
-                meas_repr = '[...]'
-            reprs['measurements'] = meas_repr
-
-        if repetition_count != 1:
-            reprs['repetition_count'] = repr(repetition_count)
-
-        kwargs = ', '.join(f'{attr}={val}' for attr, val in reprs.items())
-
-        type_name = type(self).__name__
-        max_kwargs = max(self.MAX_REPR_SIZE - len(type_name) - 2, 0)
-
-        if len(kwargs) > max_kwargs:
-            kwargs = '...'
-
-        return f'{type(self).__name__}({kwargs})'
-
-    def copy_tree_structure(self, new_parent: Union['Loop', bool]=False) -> 'Loop':
-        return type(self)(parent=self.parent if new_parent is False else new_parent,
-                          waveform=self._waveform,
-                          repetition_count=self._repetition_definition,
-                          measurements=None if self._measurements is None else list(self._measurements),
-                          children=(child.copy_tree_structure() for child in self))
-
-    def _get_measurement_windows(self, drop: bool) -> Mapping[str, np.ndarray]:
-        """Private implementation of get_measurement_windows with a slightly different data format for easier tiling.
-
-        Args:
-            drop: Drops the measurements from the Loop i.e. the Loop will no longer have measurements attached after
-            collecting them
-
-        Returns:
-             A dictionary (measurement_name -> array) with begin == array[:, 0] and length == array[:, 1]
-        """
-        temp_meas_windows = defaultdict(list)
-        if self._measurements:
-            for (mw_name, begin, length) in self._measurements:
-                temp_meas_windows[mw_name].append((begin, length))
-
-            for mw_name, begin_length_list in temp_meas_windows.items():
-                temp_meas_windows[mw_name] = [np.asarray(begin_length_list, dtype=float)]
-
-            if drop:
-                self._measurements = None
-
-        # calculate duration together with meas windows in the same iteration
-        if self.is_leaf():
-            body_duration = float(self.body_duration)
-        else:
-            offset = TimeType(0)
-            for child in self:
-                for mw_name, begins_length_array in child._get_measurement_windows(drop).items():
-                    begins_length_array[:, 0] += float(offset)
-                    temp_meas_windows[mw_name].append(begins_length_array)
-                offset += child.duration
-
-            body_duration = float(offset)
-
-        # formatting like this for easier debugging
-        result = {}
-
-        # repeat and add repetition based offset
-        for mw_name, begin_length_list in temp_meas_windows.items():
-            result[mw_name] = _repeat_loop_measurements(begin_length_list, self.repetition_count, body_duration)
-
-        return result
-
-    def get_measurement_windows(self, drop=False) -> Dict[str, Tuple[np.ndarray, np.ndarray]]:
-        """Iterates over all children and collect the begin and length arrays of each measurement window.
-
-        Args:
-            drop: Drops the measurements from the Loop i.e. the Loop will no longer have measurements attached after
-            collecting them.
-
-        Returns:
-            A dictionary (measurement_name -> (begin, length)) with begin and length being :class:`numpy.ndarray`
-        """
-        return {mw_name: (begin_length_list[:, 0], begin_length_list[:, 1])
-                for mw_name, begin_length_list in self._get_measurement_windows(drop=drop).items()}
-
-    def split_one_child(self, child_index=None) -> None:
-        """Take the last child that has a repetition count larger one, decrease it's repetition count and insert a copy
-        with repetition cout one after it"""
-        if child_index is not None:
-            if self[child_index].repetition_count < 2:
-                raise ValueError('Cannot split child {} as the repetition count is not larger 1')
-
-        else:
-            # we cannot reverse enumerate
-            n_child = len(self) - 1
-            for reverse_idx, child in enumerate(reversed(self)):
-                if child.repetition_count > 1:
-                    forward_idx = n_child - reverse_idx
-                    if not child.volatile_repetition:
-                        child_index = forward_idx
-                        break
-                    elif child_index is None:
-                        child_index = forward_idx
-            else:
-                if child_index is None:
-                    raise RuntimeError('There is no child with repetition count > 1')
-
-        if self[child_index].volatile_repetition:
-            warnings.warn("Splitting a child with volatile repetition count", VolatileModificationWarning)
-
-        new_child = self[child_index].copy_tree_structure()
-        new_child.repetition_count = 1
-
-        self[child_index].repetition_count -= 1
-
-        self[child_index+1:child_index+1] = (new_child,)
-        self.assert_tree_integrity()
-
-    def flatten_and_balance(self, depth: int) -> None:
-        """Modifies the program so all tree branches have the same depth.
-
-        Args:
-            depth: Target depth of the program
-        """
-        i = 0
-        while i < len(self):
-            # only used by type checker
-            sub_program = cast(Loop, self[i])
-
-            if sub_program.depth() < depth - 1:
-                # increase nesting because the subprogram is not deep enough
-                sub_program.encapsulate()
-
-            elif not sub_program.is_balanced():
-                # balance the sub program. We revisit it in the next iteration (no change of i )
-                # because it might modify self. While writing this comment I am not sure this is true. 14.01.2020 Simon
-                sub_program.flatten_and_balance(depth - 1)
-
-            elif sub_program.depth() == depth - 1:
-                # subprogram is balanced with the correct depth
-                i += 1
-
-            elif sub_program._has_single_child_that_can_be_merged():
-                # subprogram is balanced but to deep and has no measurements -> we can "lift" the sub-sub-program
-                # TODO: There was a len(sub_sub_program) == 1 check here that I cannot explain
-                sub_program._merge_single_child()
-
-            elif not sub_program.is_leaf():
-                # subprogram is balanced but too deep
-                sub_program.unroll()
-
-            else:
-                # we land in this case if the function gets called with depth == 0 and the current subprogram is a leaf
-                i += 1
-
-    def _has_single_child_that_can_be_merged(self) -> bool:
-        if len(self) == 1:
-            child = cast(Loop, self[0])
-            return not self._measurements or (child.repetition_count == 1 and not child.volatile_repetition)
-        else:
-            return False
-
-    def _merge_single_child(self):
-        """Lift the single child to current level. Requires _has_single_child_that_can_be_merged to be true"""
-        assert len(self) == 1, "bug: _merge_single_child called on loop with len != 1"
-        child = cast(Loop, self[0])
-
-        # if the child has a fixed repetition count of 1 the measurements can be merged
-        mergable_measurements = child.repetition_count == 1 and not child.volatile_repetition
-
-        assert not self._measurements or mergable_measurements, "bug: _merge_single_child called on loop with measurements"
-        assert not self._waveform, "bug: _merge_single_child called on loop with children and waveform"
-
-        measurements = child._measurements
-        if self._measurements:
-            if measurements:
-                measurements.extend(self._measurements)
-            else:
-                measurements = self._measurements
-
-        if not self.volatile_repetition and not child.volatile_repetition:
-            # simple integer multiplication
-            repetition_definition = self.repetition_count * child.repetition_count
-        elif not self.volatile_repetition:
-            repetition_definition = child._repetition_definition * self.repetition_count
-        elif not child.volatile_repetition:
-            repetition_definition = self._repetition_definition * child.repetition_count
-        else:
-            # create a new expression that depends on both
-            expression = 'parent_repetition_count * child_repetition_count'
-            repetition_definition = VolatileRepetitionCount.operation(
-                expression=expression,
-                parent_repetition_count=self._repetition_definition,
-                child_repetition_count=child._repetition_definition)
-
-        self[:] = iter(child)
-        self._waveform = child._waveform
-        self._repetition_definition = repetition_definition
-        self._measurements = measurements
-        self._invalidate_duration()
-        return True
-
-    def cleanup(self, actions=('remove_empty_loops', 'merge_single_child')):
-        """Apply the specified actions to cleanup the Loop.
-
-        remove_empty_loops: Remove loops with no children and no waveform (a DroppedMeasurementWarning is issued)
-        merge_single_child: see `_try_merge_single_child` documentation
-
-        Warnings:
-            DroppedMeasurementWarning: Likely a bug in qupulse. TODO: investigate whether there are usecases
-        """
-        if 'remove_empty_loops' in actions:
-            new_children = []
-            for child in self:
-                child = cast(Loop, child)
-                if child.is_leaf():
-                    if child.waveform is None:
-                        if child._measurements:
-                            warnings.warn("Dropping measurement since there is no waveform attached",
-                                          category=DroppedMeasurementWarning)
-                    else:
-                        new_children.append(child)
-
-                else:
-                    child.cleanup(actions)
-                    if child.waveform or not child.is_leaf():
-                        new_children.append(child)
-
-                    elif child._measurements:
-                        warnings.warn("Dropping measurement since there is no waveform in children",
-                                      category=DroppedMeasurementWarning)
-
-            if len(self) != len(new_children):
-                self[:] = new_children
-
-        else:
-            # only do the recursive call
-            for child in self:
-                child.cleanup(actions)
-
-        if 'merge_single_child' in actions and self._has_single_child_that_can_be_merged():
-            self._merge_single_child()
-    
-    def get_duration_structure(self) -> Tuple[int, Union[TimeType, tuple]]:
-        if self.is_leaf():
-            return self.repetition_count, self.waveform.duration
-        else:
-            return self.repetition_count, tuple(child.get_duration_structure() for child in self)
-
-    def reverse_inplace(self):
-        if self.is_leaf():
-            self._waveform = self._waveform.reversed()
-        else:
-            self._reverse_children()
-            for child in self:
-                child.reverse_inplace()
-        if self._measurements:
-            duration = self.duration
-            self._measurements = [
-                (name, duration - (begin + length), length)
-                for name, begin, length in self._measurements
-            ]
-
-
-class ChannelSplit(Exception):
-    def __init__(self, channel_sets):
-        self.channel_sets = channel_sets
-
-
-def to_waveform(program: Loop) -> Waveform:
-    if program.is_leaf():
-        if program.repetition_count == 1:
-            return program.waveform
-        else:
-            return RepetitionWaveform.from_repetition_count(program.waveform, program.repetition_count)
-    else:
-        if len(program) == 1:
-            sequenced_waveform = to_waveform(cast(Loop, program[0]))
-        else:
-            sequenced_waveform = SequenceWaveform.from_sequence(
-                [to_waveform(cast(Loop, sub_program))
-                 for sub_program in program])
-        if program.repetition_count > 1:
-            return RepetitionWaveform.from_repetition_count(sequenced_waveform, program.repetition_count)
-        else:
-            return sequenced_waveform
-
-
-class _CompatibilityLevel(Enum):
-    compatible = 0
-    action_required = 1
-    incompatible_too_short = 2
-    incompatible_fraction = 3
-    incompatible_quantum = 4
-
-    def is_incompatible(self) -> bool:
-        return self in (self.incompatible_fraction, self.incompatible_quantum, self.incompatible_too_short)
-
-
-def _is_compatible(program: Loop, min_len: int, quantum: int, sample_rate: TimeType) -> _CompatibilityLevel:
-    """ check whether program loop is compatible with awg requirements
-        possible reasons for incompatibility:
-            program shorter than minimum length
-            program duration not an integer
-            program duration not a multiple of quantum """
-    program_duration_in_samples = program.duration * sample_rate
-
-    if program_duration_in_samples.denominator != 1:
-        return _CompatibilityLevel.incompatible_fraction
-
-    if program_duration_in_samples < min_len:
-        return _CompatibilityLevel.incompatible_too_short
-
-    if program_duration_in_samples % quantum > 0:
-        return _CompatibilityLevel.incompatible_quantum
-
-    if program.is_leaf():
-        waveform_duration_in_samples = program.body_duration * sample_rate
-        if waveform_duration_in_samples < min_len or (waveform_duration_in_samples / quantum).denominator != 1:
-            if program.volatile_repetition:
-                warnings.warn("_is_compatible requires an action which drops volatility.",
-                              category=VolatileModificationWarning)
-            return _CompatibilityLevel.action_required
-        else:
-            return _CompatibilityLevel.compatible
-    else:
-        if all(_is_compatible(cast(Loop, sub_program), min_len, quantum, sample_rate) == _CompatibilityLevel.compatible
-               for sub_program in program):
-            return _CompatibilityLevel.compatible
-        else:
-            if program.volatile_repetition:
-                warnings.warn("_is_compatible requires an action which drops volatility.",
-                              category=VolatileModificationWarning)
-            return _CompatibilityLevel.action_required
-
-
-def _make_compatible(program: Loop, min_len: int, quantum: int, sample_rate: TimeType) -> None:
-    if program.is_leaf():
-        program.waveform = to_waveform(program.copy_tree_structure())
-        program.repetition_count = 1
-    else:
-        comp_levels = [_is_compatible(cast(Loop, sub_program), min_len, quantum, sample_rate)
-                       for sub_program in program]
-
-        if any(comp_level.is_incompatible() for comp_level in comp_levels):
-            single_run = program.duration * sample_rate / program.repetition_count
-            if (single_run / quantum).denominator == 1 and single_run >= min_len:
-                # it is enough to concatenate all children
-                new_repetition_definition = program.repetition_definition
-                program.repetition_count = 1
-            else:
-                # we need to concatenate all children and unroll
-                new_repetition_definition = 1
-
-            program.waveform = to_waveform(program.copy_tree_structure())
-            program.repetition_definition = new_repetition_definition
-            program[:] = []
-            return
-        else:
-            for sub_program, comp_level in zip(program, comp_levels):
-                if comp_level == _CompatibilityLevel.action_required:
-                    _make_compatible(sub_program, min_len, quantum, sample_rate)
-
-
-def make_compatible(program: Loop, minimal_waveform_length: int, waveform_quantum: int, sample_rate: TimeType):
-    """ check program for compatibility to AWG requirements, make it compatible if necessary and  possible"""
-    comp_level = _is_compatible(program,
-                                min_len=minimal_waveform_length,
-                                quantum=waveform_quantum,
-                                sample_rate=sample_rate)
-    if comp_level == _CompatibilityLevel.incompatible_fraction:
-        raise ValueError('The program duration in samples {} is not an integer'.format(program.duration * sample_rate))
-    if comp_level == _CompatibilityLevel.incompatible_too_short:
-        raise ValueError('The program is too short to be a valid waveform. \n'
-                         ' program duration in samples: {} \n'
-                         ' minimal length: {}'.format(program.duration * sample_rate, minimal_waveform_length))
-    if comp_level == _CompatibilityLevel.incompatible_quantum:
-        raise ValueError('The program duration in samples {} '
-                         'is not a multiple of quantum {}'.format(program.duration * sample_rate, waveform_quantum))
-
-    elif comp_level == _CompatibilityLevel.action_required:
-        warnings.warn("qupulse will now concatenate waveforms to make the pulse/program compatible with the chosen AWG."
-                      " This might take some time. If you need this pulse more often it makes sense to write it in a "
-                      "way which is more AWG friendly.", MakeCompatibleWarning)
-
-        _make_compatible(program,
-                         min_len=minimal_waveform_length,
-                         quantum=waveform_quantum,
-                         sample_rate=sample_rate)
-
-    else:
-        assert comp_level == _CompatibilityLevel.compatible
-
-
-def roll_constant_waveforms(program: Loop, minimal_waveform_quanta: int, waveform_quantum: int, sample_rate: TimeType):
-    """This function finds waveforms in program that can be replaced with repetitions of shorter waveforms and replaces
-    them. Complexity O(N_waveforms). Drops measurements because they are not correctly handled here for performance
-    reasons.
-
-    This is possible if:
-     - The waveform is constant on all channels
-     - waveform.duration * sample_rate / waveform_quantum has a factor that is bigger than minimal_waveform_quanta
-
-    Args:
-        program:
-        minimal_waveform_quanta:
-        waveform_quantum:
-        sample_rate:
-
-    Warnings:
-        DroppedMeasurementWarning: This warning is raised if a measurement is dropped.
-    """
-    if program._measurements:
-        warnings.warn("Dropping measurements. Remove measurements before calling roll_constant_waveforms by calling"
-                      " get_measurement_windows(drop=True).", category=DroppedMeasurementWarning)
-        program._measurements = None
-
-    waveform = program.waveform
-
-    if waveform is None:
-        for child in program:
-            roll_constant_waveforms(child, minimal_waveform_quanta, waveform_quantum, sample_rate)
-    else:
-        waveform_quanta = (waveform.duration * sample_rate) // waveform_quantum
-
-        # example
-        # waveform_quanta = 15
-        # minimal_waveform_quanta = 2
-        # => repetition_count = 5, new_waveform_quanta = 3
-        if waveform_quanta < minimal_waveform_quanta * 2:
-            # there is no way to roll this waveform because it is too short
-            return
-
-        const_values = waveform.constant_value_dict()
-        if const_values is None:
-            # The waveform is not constant
-            return
-
-        new_waveform_quanta = smallest_factor_ge(waveform_quanta, min_factor=minimal_waveform_quanta)
-        if new_waveform_quanta == waveform_quanta:
-            # the waveform duration in samples has no suitable factor
-            # TODO: Option to insert multiple Loop objects
-            return
-
-        additional_repetition_count = waveform_quanta // new_waveform_quanta
-
-        new_waveform = ConstantWaveform.from_mapping(
-            duration=waveform_quantum * new_waveform_quanta / sample_rate,
-            constant_values=const_values)
-
-        # use the private properties to avoid invalidating the duration cache of the parent loop
-        program._repetition_definition = program.repetition_definition * additional_repetition_count
-        program._waveform = new_waveform
-
-
-def _repeat_loop_measurements(begin_length_list: List[np.ndarray],
-                              repetition_count: int,
-                              body_duration: float
-                              ) -> np.ndarray:
-    temp_begin_length_array = np.concatenate(begin_length_list)
-
-    begin_length_array = np.tile(temp_begin_length_array, (repetition_count, 1))
-
-    shaped_begin_length_array = np.reshape(begin_length_array, (repetition_count, -1, 2))
-
-    shaped_begin_length_array[:, :, 0] += (np.arange(repetition_count) * body_duration)[:, np.newaxis]
-
-    return begin_length_array
-
-
-class MakeCompatibleWarning(ResourceWarning):
-    pass
-
-
-class VolatileModificationWarning(RuntimeWarning):
-    """This warning is emitted if the colatile part of a program gets modified. This might imply that the volatile
-    parameter cannot be change anymore."""
-
-
-class DroppedMeasurementWarning(RuntimeWarning):
-    """This warning is emitted if a measurement was dropped because there was no waveform attached."""
+del qupulse
diff --git a/qupulse/_program/seqc.py b/qupulse/_program/seqc.py
index 2924268b..989ffdc5 100644
--- a/qupulse/_program/seqc.py
+++ b/qupulse/_program/seqc.py
@@ -33,9 +33,9 @@
 
 from qupulse.utils.types import ChannelID, TimeType
 from qupulse.utils import replace_multiple, grouper
-from qupulse._program.waveforms import Waveform
-from qupulse._program._loop import Loop
-from qupulse._program.volatile import VolatileRepetitionCount, VolatileProperty
+from qupulse.program.waveforms import Waveform
+from qupulse.program.loop import Loop
+from qupulse.program.volatile import VolatileRepetitionCount, VolatileProperty
 from qupulse.hardware.awgs.base import ProgramEntry
 from qupulse.hardware.util import zhinst_voltage_to_uint16
 
diff --git a/qupulse/_program/tabor.py b/qupulse/_program/tabor.py
index 71ffe439..b92c76f1 100644
--- a/qupulse/_program/tabor.py
+++ b/qupulse/_program/tabor.py
@@ -12,9 +12,9 @@
 from qupulse.utils.types import ChannelID, TimeType
 from qupulse.hardware.awgs.base import ProgramEntry
 from qupulse.hardware.util import get_sample_times, voltage_to_uint16, find_positions
-from qupulse._program.waveforms import Waveform
-from qupulse._program._loop import Loop
-from qupulse._program.volatile import VolatileRepetitionCount, VolatileProperty
+from qupulse.program.waveforms import Waveform
+from qupulse.program.loop import Loop
+from qupulse.program.volatile import VolatileRepetitionCount, VolatileProperty
 
 assert(sys.byteorder == 'little')
 
diff --git a/qupulse/_program/transformation.py b/qupulse/_program/transformation.py
index 66a1641f..c6a0d0de 100644
--- a/qupulse/_program/transformation.py
+++ b/qupulse/_program/transformation.py
@@ -1,400 +1,7 @@
-from typing import Any, Mapping, Set, Tuple, Sequence, AbstractSet, Union, TYPE_CHECKING, Hashable
-from abc import abstractmethod
-from numbers import Real
+from qupulse.program.transformation import *
 
-import numpy as np
+import qupulse.program.transformation
 
-from qupulse import ChannelID
-from qupulse.comparable import Comparable
-from qupulse.utils.types import SingletonABCMeta, frozendict
-from qupulse.expressions import ExpressionScalar
+__all__ = qupulse.program.transformation.__all__
 
-
-_TrafoValue = Union[Real, ExpressionScalar]
-
-
-class Transformation(Comparable):
-    _identity_singleton = None
-    """Transforms numeric time-voltage values for multiple channels to other time-voltage values. The number and names
-     of input and output channels might differ."""
-
-    @abstractmethod
-    def __call__(self, time: Union[np.ndarray, float],
-                 data: Mapping[ChannelID, Union[np.ndarray, float]]) -> Mapping[ChannelID, Union[np.ndarray, float]]:
-        """Apply transformation to data
-        Args:
-            time:
-            data:
-
-        Returns:
-            transformed: A DataFrame that has been transformed with index == output_channels
-        """
-
-    @abstractmethod
-    def get_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        """Return the channel identifiers"""
-
-    @abstractmethod
-    def get_input_channels(self, output_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        """Channels that are required for getting data for the requested output channel"""
-
-    def chain(self, next_transformation: 'Transformation') -> 'Transformation':
-        if next_transformation is IdentityTransformation():
-            return self
-        else:
-            return chain_transformations(self, next_transformation)
-
-    def is_constant_invariant(self):
-        """Signals if the transformation always maps constants to constants."""
-        return False
-
-    def get_constant_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        return frozenset()
-
-
-class IdentityTransformation(Transformation, metaclass=SingletonABCMeta):
-    def __call__(self, time: Union[np.ndarray, float],
-                 data: Mapping[ChannelID, Union[np.ndarray, float]]) -> Mapping[ChannelID, Union[np.ndarray, float]]:
-        return data
-
-    def get_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        return input_channels
-
-    @property
-    def compare_key(self) -> None:
-        return None
-
-    def get_input_channels(self, output_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        return output_channels
-
-    def chain(self, next_transformation: Transformation) -> Transformation:
-        return next_transformation
-
-    def __repr__(self):
-        return 'IdentityTransformation()'
-
-    def is_constant_invariant(self):
-        """Signals if the transformation always maps constants to constants."""
-        return True
-
-    def get_constant_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        return input_channels
-
-
-class ChainedTransformation(Transformation):
-    def __init__(self, *transformations: Transformation):
-        self._transformations = transformations
-
-    @property
-    def transformations(self) -> Tuple[Transformation, ...]:
-        return self._transformations
-
-    def get_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        for transformation in self._transformations:
-            input_channels = transformation.get_output_channels(input_channels)
-        return input_channels
-
-    def get_input_channels(self, output_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        for transformation in reversed(self._transformations):
-            output_channels = transformation.get_input_channels(output_channels)
-        return output_channels
-
-    def __call__(self, time: Union[np.ndarray, float],
-                 data: Mapping[ChannelID, Union[np.ndarray, float]]) -> Mapping[ChannelID, Union[np.ndarray, float]]:
-        for transformation in self._transformations:
-            data = transformation(time, data)
-        return data
-
-    @property
-    def compare_key(self) -> Tuple[Transformation, ...]:
-        return self._transformations
-
-    def chain(self, next_transformation) -> Transformation:
-        return chain_transformations(*self.transformations, next_transformation)
-
-    def __repr__(self):
-        return f'{type(self).__name__}{self._transformations!r}'
-
-    def is_constant_invariant(self):
-        """Signals if the transformation always maps constants to constants."""
-        return all(trafo.is_constant_invariant() for trafo in self._transformations)
-
-    def get_constant_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        for trafo in self._transformations:
-            input_channels = trafo.get_constant_output_channels(input_channels)
-        return input_channels
-
-
-class LinearTransformation(Transformation):
-    def __init__(self,
-                 transformation_matrix: np.ndarray,
-                 input_channels: Sequence[ChannelID],
-                 output_channels: Sequence[ChannelID]):
-        """
-
-        Args:
-            transformation_matrix: Matrix describing the transformation with shape (output_channels, input_channels)
-            input_channels: Channel ids of the columns
-            output_channels: Channel ids of the rows
-        """
-        transformation_matrix = np.asarray(transformation_matrix)
-
-        if transformation_matrix.shape != (len(output_channels), len(input_channels)):
-            raise ValueError('Shape of transformation matrix does not match to the given channels')
-
-        output_sorter = np.argsort(output_channels)
-        transformation_matrix = transformation_matrix[output_sorter, :]
-
-        input_sorter = np.argsort(input_channels)
-        transformation_matrix = transformation_matrix[:, input_sorter]
-
-        self._matrix = transformation_matrix
-        self._input_channels = tuple(sorted(input_channels))
-        self._output_channels = tuple(sorted(output_channels))
-        self._input_channels_set = frozenset(self._input_channels)
-        self._output_channels_set = frozenset(self._output_channels)
-
-    def __call__(self, time: Union[np.ndarray, float],
-                 data: Mapping[ChannelID, Union[np.ndarray, float]]) -> Mapping[ChannelID, Union[np.ndarray, float]]:
-        data_out = {forwarded_channel: data[forwarded_channel]
-                    for forwarded_channel in set(data.keys()).difference(self._input_channels)}
-
-        if len(data_out) == len(data):
-            # only forwarded data
-            return data_out
-
-        try:
-            data_in = np.stack([data[in_channel] for in_channel in self._input_channels])
-        except KeyError as error:
-            raise KeyError('Invalid input channels', set(data.keys()), set(self._input_channels)) from error
-
-        transformed_data = self._matrix @ data_in
-
-        for idx, out_channel in enumerate(self._output_channels):
-            data_out[out_channel] = transformed_data[idx, ...]
-
-        return data_out
-
-    def get_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        if not input_channels >= self._input_channels_set:
-            # input_channels is not a superset of the required input channels
-            raise KeyError('Invalid input channels', input_channels, self._input_channels_set)
-
-        return (input_channels - self._input_channels_set) | self._output_channels_set
-
-    def get_input_channels(self, output_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        forwarded = output_channels - self._output_channels_set
-        if not forwarded.isdisjoint(self._input_channels):
-            raise KeyError('Is input channel', forwarded & self._input_channels_set)
-        elif output_channels.isdisjoint(self._output_channels):
-            return output_channels
-        else:
-            return forwarded | self._input_channels_set
-
-    @property
-    def compare_key(self) -> Tuple[Tuple[ChannelID], Tuple[ChannelID], bytes]:
-        return self._input_channels, self._output_channels, self._matrix.tobytes()
-
-    def __repr__(self):
-        return ('LinearTransformation('
-                'transformation_matrix={transformation_matrix},'
-                'input_channels={input_channels},'
-                'output_channels={output_channels})').format(transformation_matrix=self._matrix.tolist(),
-                                                             input_channels=self._input_channels,
-                                                             output_channels=self._output_channels)
-
-    def is_constant_invariant(self):
-        """Signals if the transformation always maps constants to constants."""
-        return True
-
-    def get_constant_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        return input_channels
-
-
-class OffsetTransformation(Transformation):
-    def __init__(self, offsets: Mapping[ChannelID, _TrafoValue]):
-        """Adds an offset to each channel specified in offsets.
-
-        Channels not in offsets are forewarded
-
-        Args:
-            offsets: Channel -> offset mapping
-        """
-        self._offsets = frozendict(offsets)
-        assert _are_valid_transformation_expressions(self._offsets), f"Not valid transformation expressions: {self._offsets}"
-
-    def __call__(self, time: Union[np.ndarray, float],
-                 data: Mapping[ChannelID, Union[np.ndarray, float]]) -> Mapping[ChannelID, Union[np.ndarray, float]]:
-        offsets = _instantiate_expression_dict(time, self._offsets)
-        return {channel: channel_values + offsets[channel] if channel in offsets else channel_values
-                for channel, channel_values in data.items()}
-
-    def get_input_channels(self, output_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        return output_channels
-
-    def get_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        return input_channels
-
-    @property
-    def compare_key(self) -> Hashable:
-        return self._offsets
-
-    def __repr__(self):
-        return f'{type(self).__name__}({dict(self._offsets)!r})'
-
-    def is_constant_invariant(self):
-        """Signals if the transformation always maps constants to constants."""
-        return not _has_time_dependent_values(self._offsets)
-
-    def get_constant_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        return _get_constant_output_channels(self._offsets, input_channels)
-
-
-class ScalingTransformation(Transformation):
-    def __init__(self, factors: Mapping[ChannelID, _TrafoValue]):
-        self._factors = frozendict(factors)
-        assert _are_valid_transformation_expressions(self._factors), f"Not valid transformation expressions: {self._factors}"
-
-    def __call__(self, time: Union[np.ndarray, float],
-                 data: Mapping[ChannelID, Union[np.ndarray, float]]) -> Mapping[ChannelID, Union[np.ndarray, float]]:
-        factors = _instantiate_expression_dict(time, self._factors)
-        return {channel: channel_values * factors[channel] if channel in factors else channel_values
-                for channel, channel_values in data.items()}
-
-    def get_input_channels(self, output_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        return output_channels
-
-    def get_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        return input_channels
-
-    @property
-    def compare_key(self) -> Hashable:
-        return self._factors
-
-    def __repr__(self):
-        return f'{type(self).__name__}({dict(self._factors)!r})'
-
-    def is_constant_invariant(self):
-        """Signals if the transformation always maps constants to constants."""
-        return not _has_time_dependent_values(self._factors)
-
-    def get_constant_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        return _get_constant_output_channels(self._factors, input_channels)
-
-
-try:
-    if TYPE_CHECKING:
-        import pandas
-        PandasDataFrameType = pandas.DataFrame
-    else:
-        PandasDataFrameType = Any
-
-    def linear_transformation_from_pandas(transformation: PandasDataFrameType) -> LinearTransformation:
-        """ Creates a LinearTransformation object out of a pandas data frame.
-
-        Args:
-            transformation (pandas.DataFrame): The pandas.DataFrame object out of which a LinearTransformation will be formed.
-
-        Returns:
-            the created LinearTransformation instance
-        """
-        return LinearTransformation(transformation.values, transformation.columns, transformation.index)
-
-    LinearTransformation.from_pandas = linear_transformation_from_pandas
-except ImportError:
-    pass
-
-
-class ParallelChannelTransformation(Transformation):
-    def __init__(self, channels: Mapping[ChannelID, _TrafoValue]):
-        """Set channel values to given values regardless their former existence. The values can be time dependent
-        expressions.
-
-        Args:
-            channels: Channels present in this map are set to the given value.
-        """
-        self._channels: Mapping[ChannelID, _TrafoValue] = frozendict(channels.items())
-        assert _are_valid_transformation_expressions(self._channels), f"Not valid transformation expressions: {self._channels}"
-
-    def __call__(self, time: Union[np.ndarray, float],
-                 data: Mapping[ChannelID, Union[np.ndarray, float]]) -> Mapping[ChannelID, Union[np.ndarray, float]]:
-        overwritten = self._instantiated_values(time)
-        return {**data, **overwritten}
-
-    def _instantiated_values(self, time):
-        scope = {'t': time}
-        return {channel: value.evaluate_in_scope(scope) if hasattr(value, 'evaluate_in_scope') else np.full_like(time, fill_value=value, dtype=float)
-                for channel, value in self._channels.items()}
-
-    @property
-    def compare_key(self) -> Hashable:
-        return self._channels
-
-    def get_input_channels(self, output_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        return output_channels - self._channels.keys()
-
-    def get_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        return input_channels | self._channels.keys()
-
-    def __repr__(self):
-        return f'{type(self).__name__}({dict(self._channels)!r})'
-
-    def is_constant_invariant(self):
-        """Signals if the transformation always maps constants to constants."""
-        return not _has_time_dependent_values(self._channels)
-
-    def get_constant_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-        output_channels = set(input_channels)
-        for ch, value in self._channels.items():
-            if hasattr(value, 'variables'):
-                output_channels.discard(ch)
-            else:
-                output_channels.add(ch)
-
-        return output_channels
-
-
-def chain_transformations(*transformations: Transformation) -> Transformation:
-    parsed_transformations = []
-    for transformation in transformations:
-        if transformation is IdentityTransformation() or transformation is None:
-            pass
-        elif isinstance(transformation, ChainedTransformation):
-            parsed_transformations.extend(transformation.transformations)
-        else:
-            parsed_transformations.append(transformation)
-    if len(parsed_transformations) == 0:
-        return IdentityTransformation()
-    elif len(parsed_transformations) == 1:
-        return parsed_transformations[0]
-    else:
-        return ChainedTransformation(*parsed_transformations)
-
-
-def _instantiate_expression_dict(time, expressions: Mapping[str, _TrafoValue]) -> Mapping[str, Union[Real, np.ndarray]]:
-    scope = {'t': time}
-    modified_expressions = {}
-    for name, value in expressions.items():
-        if hasattr(value, 'evaluate_in_scope'):
-            modified_expressions[name] = value.evaluate_in_scope(scope)
-    if modified_expressions:
-        return {**expressions, **modified_expressions}
-    else:
-        return expressions
-
-
-def _has_time_dependent_values(expressions: Mapping[ChannelID, _TrafoValue]) -> bool:
-    return any(hasattr(value, 'variables')
-               for value in expressions.values())
-
-
-def _get_constant_output_channels(expressions: Mapping[ChannelID, _TrafoValue],
-                                  constant_input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
-    return {ch
-            for ch in constant_input_channels
-            if not hasattr(expressions.get(ch, None), 'variables')}
-
-def _are_valid_transformation_expressions(expressions: Mapping[ChannelID, _TrafoValue]) -> bool:
-    return all(expr.variables == ('t',)
-               for expr in expressions.values()
-               if hasattr(expr, 'variables'))
+del qupulse
diff --git a/qupulse/_program/volatile.py b/qupulse/_program/volatile.py
index ab76e8ec..ddfe2aa1 100644
--- a/qupulse/_program/volatile.py
+++ b/qupulse/_program/volatile.py
@@ -1,70 +1,7 @@
-from typing import NamedTuple, Mapping
-import warnings
-import numbers
+from qupulse.program.volatile import *
 
+import qupulse.program.volatile
 
-from qupulse.parameter_scope import Scope, MappedScope, JointScope
-from qupulse.expressions import Expression, ExpressionScalar
-from qupulse.utils.types import FrozenDict, FrozenMapping
-from qupulse.utils import is_integer
+__all__ = qupulse.program.volatile.__all__
 
-
-VolatileProperty = NamedTuple('VolatileProperty', [('expression', Expression),
-                                                   ('dependencies', FrozenMapping[str, Expression])])
-VolatileProperty.__doc__ = """Hashable representation of a volatile program property. It does not contain the concrete
-value. Using the dependencies attribute to calculate the value might yield unexpected results."""
-
-
-class VolatileValue:
-    """Not hashable"""
-
-    def __init__(self, expression: ExpressionScalar, scope: Scope):
-        self._expression = expression
-        self._scope = scope
-
-    @property
-    def volatile_property(self) -> VolatileProperty:
-        dependencies = self._scope.get_volatile_parameters()
-        dependencies = FrozenDict({parameter_name: dependencies[parameter_name]
-                                   for parameter_name in self._expression.variables
-                                   if parameter_name in dependencies})
-        return VolatileProperty(expression=self._expression, dependencies=dependencies)
-
-    @classmethod
-    def operation(cls, expression, **operands):
-        expression = Expression(expression)
-        assert set(expression.variables) == operands.keys()
-
-        scope = JointScope(FrozenDict(
-            {operand_name: MappedScope(operand._scope, FrozenDict({operand_name: operand._expression}))
-             for operand_name, operand in operands.items()}
-        ))
-        return cls(expression, scope)
-    
-    def __sub__(self, other: int):
-        return type(self)(self._expression - other, self._scope)
-
-    def __mul__(self, other: int):
-        return type(self)(self._expression * other, self._scope)
-
-
-class VolatileRepetitionCount(VolatileValue):
-    def __int__(self):
-        value = self._expression.evaluate_in_scope(self._scope)
-        if not is_integer(value):
-            warnings.warn("Repetition count is no integer. Rounding might lead to unexpected results.")
-        value = int(round(value))
-        if value < 0:
-            warnings.warn("Repetition count is negative. Clamping lead to unexpected results.")
-            value = 0
-        return value
-
-    def update_volatile_dependencies(self, new_constants: Mapping[str, numbers.Number]) -> int:
-        self._scope = self._scope.change_constants(new_constants)
-        return int(self)
-
-    def __eq__(self, other):
-        if type(self) == type(other):
-            return self._scope is other._scope and self._expression == other._expression
-        else:
-            return NotImplemented
+del qupulse
diff --git a/qupulse/_program/waveforms.py b/qupulse/_program/waveforms.py
index a173f3bf..5038c398 100644
--- a/qupulse/_program/waveforms.py
+++ b/qupulse/_program/waveforms.py
@@ -1,1233 +1,9 @@
-"""This module contains all waveform classes
+"""Backwards compatibility link to qupulse.program.waveforms"""
 
-Classes:
-    - Waveform: An instantiated pulse which can be sampled to a raw voltage value array.
-"""
+from qupulse.program.waveforms import *
 
-import itertools
-import operator
-import warnings
-from abc import ABCMeta, abstractmethod
-from numbers import Real
-from typing import (
-    AbstractSet, Any, FrozenSet, Iterable, Mapping, NamedTuple, Sequence, Set,
-    Tuple, Union, cast, Optional, List, Hashable)
-from weakref import WeakValueDictionary, ref
+import qupulse.program.waveforms
 
-import numpy as np
+__all__ = qupulse.program.waveforms.__all__
 
-from qupulse import ChannelID
-from qupulse._program.transformation import Transformation
-from qupulse.utils import checked_int_cast, isclose
-from qupulse.utils.types import TimeType, time_from_float
-from qupulse.utils.performance import is_monotonic
-from qupulse.comparable import Comparable
-from qupulse.expressions import ExpressionScalar
-from qupulse.pulses.interpolation import InterpolationStrategy
-from qupulse.utils import checked_int_cast, isclose
-from qupulse.utils.types import TimeType, time_from_float, FrozenDict
-from qupulse._program.transformation import Transformation
-from qupulse.utils import pairwise
-
-class ConstantFunctionPulseTemplateWarning(UserWarning):
-    """  This warning indicates a constant waveform is constructed from a FunctionPulseTemplate """
-    pass
-
-__all__ = ["Waveform", "TableWaveform", "TableWaveformEntry", "FunctionWaveform", "SequenceWaveform",
-           "MultiChannelWaveform", "RepetitionWaveform", "TransformingWaveform", "ArithmeticWaveform",
-           "ConstantFunctionPulseTemplateWarning"]
-
-PULSE_TO_WAVEFORM_ERROR = None  # error margin in pulse template to waveform conversion
-
-#  these are private because there probably will be changes here
-_ALLOCATION_FUNCTION = np.full_like  # pre_allocated = ALLOCATION_FUNCTION(sample_times, **ALLOCATION_FUNCTION_KWARGS)
-_ALLOCATION_FUNCTION_KWARGS = dict(fill_value=np.nan, dtype=float)
-
-
-def _to_time_type(duration: Real) -> TimeType:
-    if isinstance(duration, TimeType):
-        return duration
-    else:
-        return time_from_float(float(duration), absolute_error=PULSE_TO_WAVEFORM_ERROR)
-
-
-class Waveform(Comparable, metaclass=ABCMeta):
-    """Represents an instantiated PulseTemplate which can be sampled to retrieve arrays of voltage
-    values for the hardware."""
-
-    __sampled_cache = WeakValueDictionary()
-
-    __slots__ = ('_duration',)
-
-    def __init__(self, duration: TimeType):
-        self._duration = duration
-
-    @property
-    def duration(self) -> TimeType:
-        """The duration of the waveform in time units."""
-        return self._duration
-
-    @abstractmethod
-    def unsafe_sample(self,
-                      channel: ChannelID,
-                      sample_times: np.ndarray,
-                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
-        """Sample the waveform at given sample times.
-
-        The unsafe means that there are no sanity checks performed. The provided sample times are assumed to be
-        monotonously increasing and lie in the range of [0, waveform.duration]
-
-        Args:
-            sample_times: Times at which this Waveform will be sampled.
-            output_array: Has to be either None or an array of the same size and type as sample_times. If
-                not None, the sampled values will be written here and this array will be returned
-        Result:
-            The sampled values of this Waveform at the provided sample times. Has the same number of
-            elements as sample_times.
-        """
-
-    def get_sampled(self,
-                    channel: ChannelID,
-                    sample_times: np.ndarray,
-                    output_array: Union[np.ndarray, None] = None) -> np.ndarray:
-        """A wrapper to the unsafe_sample method which caches the result. This method enforces the constrains
-        unsafe_sample expects and caches the result to save memory.
-
-        Args:
-            sample_times: Times at which this Waveform will be sampled.
-            output_array: Has to be either None or an array of the same size and type as sample_times. If an array is
-                given, the sampled values will be written into the given array and it will be returned. Otherwise, a new
-                array will be created and cached to save memory.
-
-        Result:
-            The sampled values of this Waveform at the provided sample times. Is `output_array` if provided
-        """
-        if len(sample_times) == 0:
-            if output_array is None:
-                return np.zeros_like(sample_times, dtype=float)
-            elif len(output_array) == len(sample_times):
-                return output_array
-            else:
-                raise ValueError('Output array length and sample time length are different')
-
-        if not is_monotonic(sample_times):
-            raise ValueError('The sample times are not monotonously increasing')
-        if sample_times[0] < 0 or sample_times[-1] > float(self.duration):
-            raise ValueError(f'The sample times [{sample_times[0]}, ..., {sample_times[-1]}] are not in the range'
-                             f' [0, duration={float(self.duration)}]')
-        if channel not in self.defined_channels:
-            raise KeyError('Channel not defined in this waveform: {}'.format(channel))
-
-        constant_value = self.constant_value(channel)
-        if constant_value is None:
-            if output_array is None:
-                # cache the result to save memory
-                result = self.unsafe_sample(channel, sample_times)
-                result.flags.writeable = False
-                key = hash(bytes(result))
-                if key not in self.__sampled_cache:
-                    self.__sampled_cache[key] = result
-                return self.__sampled_cache[key]
-            else:
-                if len(output_array) != len(sample_times):
-                    raise ValueError('Output array length and sample time length are different')
-                # use the user provided memory
-                return self.unsafe_sample(channel=channel,
-                                          sample_times=sample_times,
-                                          output_array=output_array)
-        else:
-            if output_array is None:
-                output_array = np.full_like(sample_times, fill_value=constant_value, dtype=float)
-            else:
-                output_array[:] = constant_value
-            return output_array
-
-    @property
-    @abstractmethod
-    def defined_channels(self) -> AbstractSet[ChannelID]:
-        """The channels this waveform should played on. Use
-            :func:`~qupulse.pulses.instructions.get_measurement_windows` to get a waveform for a subset of these."""
-
-    @abstractmethod
-    def unsafe_get_subset_for_channels(self, channels: AbstractSet[ChannelID]) -> 'Waveform':
-        """Unsafe version of :func:`~qupulse.pulses.instructions.get_measurement_windows`."""
-
-    def get_subset_for_channels(self, channels: AbstractSet[ChannelID]) -> 'Waveform':
-        """Get a waveform that only describes the channels contained in `channels`.
-
-        Args:
-            channels: A channel set the return value should confine to.
-
-        Raises:
-            KeyError: If `channels` is not a subset of the waveform's defined channels.
-
-        Returns:
-            A waveform with waveform.defined_channels == `channels`
-        """
-        if not channels <= self.defined_channels:
-            raise KeyError('Channels not defined on waveform: {}'.format(channels))
-        if channels == self.defined_channels:
-            return self
-        return self.unsafe_get_subset_for_channels(channels=channels)
-
-    def is_constant(self) -> bool:
-        """Convenience function to check if all channels are constant. The result is equal to
-        `all(waveform.constant_value(ch) is not None for ch in waveform.defined_channels)` but might be more performant.
-
-        Returns:
-            True if all channels have constant values.
-        """
-        return self.constant_value_dict() is not None
-
-    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
-        result = {ch: self.constant_value(ch) for ch in self.defined_channels}
-        if None in result.values():
-            return None
-        else:
-            return result
-
-    def constant_value(self, channel: ChannelID) -> Optional[float]:
-        """Checks if the requested channel has a constant value and returns it if so.
-
-        Guarantee that this assertion passes for every t in waveform duration:
-        >>> assert waveform.constant_value(channel) is None or waveform.constant_value(t) = waveform.get_sampled(channel, t)
-
-        Args:
-            channel: The channel to check
-
-        Returns:
-            None if there is no guarantee that the channel is constant. The value otherwise.
-        """
-        return None
-
-    def __neg__(self):
-        return FunctorWaveform.from_functor(self, {ch: np.negative for ch in self.defined_channels})
-
-    def __pos__(self):
-        return self
-
-    def _sort_key_for_channels(self) -> Sequence[Tuple[str, int]]:
-        """Makes reproducible sorting by defined channels possible"""
-        return sorted((ch, 0) if isinstance(ch, str) else ('', ch) for ch in self.defined_channels)
-
-    def reversed(self) -> 'Waveform':
-        """Returns a reversed version of this waveform."""
-        # We don't check for constness here because const waveforms are supposed to override this method
-        return ReversedWaveform(self)
-
-
-class TableWaveformEntry(NamedTuple('TableWaveformEntry', [('t', Real),
-                                                           ('v', float),
-                                                           ('interp', InterpolationStrategy)])):
-    def __init__(self, t: float, v: float, interp: InterpolationStrategy):
-        if not callable(interp):
-            raise TypeError('{} is neither callable nor of type InterpolationStrategy'.format(interp))
-
-    def __repr__(self):
-        return f'{type(self).__name__}(t={self.t!r}, v={self.v!r}, interp={self.interp!r})'
-
-
-class TableWaveform(Waveform):
-    EntryInInit = Union[TableWaveformEntry, Tuple[float, float, InterpolationStrategy]]
-
-    """Waveform obtained from instantiating a TablePulseTemplate."""
-
-    __slots__ = ('_table', '_channel_id')
-
-    def __init__(self,
-                 channel: ChannelID,
-                 waveform_table: Tuple[TableWaveformEntry, ...]) -> None:
-        """Create a new TableWaveform instance.
-
-        Args:
-            waveform_table: A tuple of instantiated and validated table entries
-        """
-        if not isinstance(waveform_table, tuple):
-            warnings.warn("Please use a tuple of TableWaveformEntry to construct TableWaveform directly",
-                          category=DeprecationWarning)
-            waveform_table = self._validate_input(waveform_table)
-
-        super().__init__(duration=_to_time_type(waveform_table[-1].t))
-
-        self._table = waveform_table
-        self._channel_id = channel
-
-    @staticmethod
-    def _validate_input(input_waveform_table: Sequence[EntryInInit]) -> Union[Tuple[Real, Real],
-                                                                              List[TableWaveformEntry]]:
-        """ Checks that:
-         - the time is increasing,
-         - there are at least two entries
-
-        Optimizations:
-          - removes subsequent entries with same time or voltage values.
-          - checks if the complete waveform is constant. Returns a (duration, value) tuple if this is the case
-
-        Raises:
-            ValueError:
-              - there are less than two entries
-              - the entries are not ordered in time
-              - Any time is negative
-              - The total length is zero
-
-        Returns:
-            A list of de-duplicated table entries
-            OR
-            A (duration, value) tuple if the waveform is constant
-        """
-        # we use an iterator here to avoid duplicate work and be maximally efficient for short tables
-        # We never use StopIteration to abort iteration. It always signifies an error.
-        input_iter = iter(input_waveform_table)
-        try:
-            first_t, first_v, first_interp = next(input_iter)
-        except StopIteration:
-            raise ValueError("Waveform table mut not be empty")
-
-        if first_t != 0.0:
-            raise ValueError('First time entry is not zero.')
-
-        previous_t = 0.0
-        previous_v = first_v
-        output_waveform_table = [TableWaveformEntry(0.0, first_v, first_interp)]
-
-        try:
-            t, v, interp = next(input_iter)
-        except StopIteration:
-            raise ValueError("Waveform table has less than two entries.")
-        if t < 0:
-            raise ValueError('Negative time values are not allowed.')
-
-        # constant_v is None <=> the waveform is constant until up to the current entry
-        constant_v = interp.constant_value((previous_t, previous_v), (t, v))
-
-        for next_t, next_v, next_interp in input_iter:
-            if next_t < t:
-                if next_t < 0:
-                    raise ValueError('Negative time values are not allowed.')
-                else:
-                    raise ValueError('Times are not increasing.')
-
-            if constant_v is not None and interp.constant_value((t, v), (next_t, next_v)) != constant_v:
-                constant_v = None
-
-            if (previous_t != t or t != next_t) and (previous_v != v or v != next_v):
-                # the time and the value differ both either from the next or the previous
-                # otherwise we skip the entry
-                previous_t = t
-                previous_v = v
-                output_waveform_table.append(TableWaveformEntry(t, v, interp))
-
-            t, v, interp = next_t, next_v, next_interp
-
-        # Until now, we only checked that the time does not decrease. We require an increase because duration == 0
-        # waveforms are ill-formed. t is now the time of the last entry.
-        if t == 0:
-            raise ValueError('Last time entry is zero.')
-
-        if constant_v is not None:
-            # the waveform is constant
-            return t, constant_v
-        else:
-            # we must still add the last entry to the table
-            output_waveform_table.append(TableWaveformEntry(t, v, interp))
-            return output_waveform_table
-
-    def is_constant(self) -> bool:
-        # only correct if `from_table` is used
-        return False
-
-    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
-        # only correct if `from_table` is used
-        return None
-
-    @classmethod
-    def from_table(cls, channel: ChannelID, table: Sequence[EntryInInit]) -> Union['TableWaveform', 'ConstantWaveform']:
-        table = cls._validate_input(table)
-        if isinstance(table, tuple):
-            duration, amplitude = table
-            return ConstantWaveform(duration=duration, amplitude=amplitude, channel=channel)
-        else:
-            return TableWaveform(channel, tuple(table))
-
-    @property
-    def compare_key(self) -> Any:
-        return self._channel_id, self._table
-
-    def unsafe_sample(self,
-                      channel: ChannelID,
-                      sample_times: np.ndarray,
-                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
-        if output_array is None:
-            output_array = _ALLOCATION_FUNCTION(sample_times, **_ALLOCATION_FUNCTION_KWARGS)
-
-        if PULSE_TO_WAVEFORM_ERROR:
-            # we need to replace the last entry's t with self.duration
-            *entries, last = self._table
-            entries.append(TableWaveformEntry(float(self.duration), last.v, last.interp))
-        else:
-            entries = self._table
-
-        for entry1, entry2 in pairwise(entries):
-            indices = slice(np.searchsorted(sample_times, entry1.t, 'left'),
-                            np.searchsorted(sample_times, entry2.t, 'right'))
-            output_array[indices] = \
-                entry2.interp((float(entry1.t), entry1.v),
-                              (float(entry2.t), entry2.v),
-                              sample_times[indices])
-        return output_array
-
-    @property
-    def defined_channels(self) -> AbstractSet[ChannelID]:
-        return {self._channel_id}
-
-    def unsafe_get_subset_for_channels(self, channels: AbstractSet[ChannelID]) -> 'Waveform':
-        return self
-
-    def __repr__(self):
-        return f'{type(self).__name__}(channel={self._channel_id!r}, waveform_table={self._table!r})'
-
-
-class ConstantWaveform(Waveform):
-
-    # TODO: remove
-    _is_constant_waveform = True
-
-    __slots__ = ('_amplitude', '_channel')
-
-    def __init__(self, duration: Real, amplitude: Any, channel: ChannelID):
-        """ Create a qupulse waveform corresponding to a ConstantPulseTemplate """
-        super().__init__(duration=_to_time_type(duration))
-        self._amplitude = amplitude
-        self._channel = channel
-
-    @classmethod
-    def from_mapping(cls, duration: Real, constant_values: Mapping[ChannelID, float]) -> Union['ConstantWaveform',
-                                                                                               'MultiChannelWaveform']:
-        """Construct a ConstantWaveform or a MultiChannelWaveform of ConstantWaveforms with given duration and values"""
-        assert constant_values
-        duration = _to_time_type(duration)
-        if len(constant_values) == 1:
-            (channel, amplitude), = constant_values.items()
-            return cls(duration, amplitude=amplitude, channel=channel)
-        else:
-            return MultiChannelWaveform([cls(duration, amplitude=amplitude, channel=channel)
-                                         for channel, amplitude in constant_values.items()])
-
-    def is_constant(self) -> bool:
-        return True
-
-    def constant_value(self, channel: ChannelID) -> Optional[float]:
-        assert channel == self._channel
-        return self._amplitude
-
-    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
-        return {self._channel: self._amplitude}
-
-    @property
-    def defined_channels(self) -> AbstractSet[ChannelID]:
-        """The channels this waveform should played on. Use
-            :func:`~qupulse.pulses.instructions.get_measurement_windows` to get a waveform for a subset of these."""
-
-        return {self._channel}
-
-    @property
-    def compare_key(self) -> Tuple[Any, ...]:
-        return self._duration, self._amplitude, self._channel
-
-    def unsafe_sample(self,
-                      channel: ChannelID,
-                      sample_times: np.ndarray,
-                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
-        if output_array is None:
-            return np.full_like(sample_times, fill_value=self._amplitude, dtype=float)
-        else:
-            output_array[:] = self._amplitude
-            return output_array
-
-    def unsafe_get_subset_for_channels(self, channels: Set[ChannelID]) -> Waveform:
-        """Unsafe version of :func:`~qupulse.pulses.instructions.get_measurement_windows`."""
-        return self
-
-    def __repr__(self):
-        return f"{type(self).__name__}(duration={self.duration!r}, "\
-               f"amplitude={self._amplitude!r}, channel={self._channel!r})"
-
-    def reversed(self) -> 'Waveform':
-        return self
-
-
-class FunctionWaveform(Waveform):
-    """Waveform obtained from instantiating a FunctionPulseTemplate."""
-
-    __slots__ = ('_expression', '_channel_id')
-
-    def __init__(self, expression: ExpressionScalar,
-                 duration: float,
-                 channel: ChannelID) -> None:
-        """Creates a new FunctionWaveform instance.
-
-        Args:
-            expression: The function represented by this FunctionWaveform
-                as a mathematical expression where 't' denotes the time variable. It must not have other variables
-            duration: The duration of the waveform
-            measurement_windows: A list of measurement windows
-            channel: The channel this waveform is played on
-        """
-
-        if set(expression.variables) - set('t'):
-            raise ValueError('FunctionWaveforms may not depend on anything but "t"')
-        elif not expression.variables:
-            warnings.warn("Constant FunctionWaveform is not recommended as the constant propagation will be suboptimal",
-                          category=ConstantFunctionPulseTemplateWarning)
-        super().__init__(duration=_to_time_type(duration))
-        self._expression = expression
-        self._channel_id = channel
-
-    @classmethod
-    def from_expression(cls, expression: ExpressionScalar, duration: float, channel: ChannelID) -> Union['FunctionWaveform', ConstantWaveform]:
-        if expression.variables:
-            return cls(expression, duration, channel)
-        else:
-            return ConstantWaveform(amplitude=expression.evaluate_numeric(), duration=duration, channel=channel)
-
-    def is_constant(self) -> bool:
-        # only correct if `from_expression` is used
-        return False
-
-    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
-        # only correct if `from_expression` is used
-        return None
-
-    @property
-    def defined_channels(self) -> AbstractSet[ChannelID]:
-        return {self._channel_id}
-
-    @property
-    def compare_key(self) -> Any:
-        return self._channel_id, self._expression, self._duration
-
-    @property
-    def duration(self) -> TimeType:
-        return self._duration
-
-    def unsafe_sample(self,
-                      channel: ChannelID,
-                      sample_times: np.ndarray,
-                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
-        evaluated = self._expression.evaluate_numeric(t=sample_times)
-        if output_array is None:
-            if self._expression.variables:
-                return evaluated.astype(float)
-            else:
-                return np.full_like(sample_times, fill_value=float(evaluated), dtype=float)
-        else:
-            output_array[:] = evaluated
-            return output_array
-
-    def unsafe_get_subset_for_channels(self, channels: AbstractSet[ChannelID]) -> Waveform:
-        return self
-
-    def __repr__(self):
-        return f"{type(self).__name__}(duration={self.duration!r}, "\
-               f"expression={self._expression!r}, channel={self._channel_id!r})"
-
-
-class SequenceWaveform(Waveform):
-    """This class allows putting multiple PulseTemplate together in one waveform on the hardware."""
-
-    __slots__ = ('_sequenced_waveforms', )
-
-    def __init__(self, sub_waveforms: Iterable[Waveform]):
-        """Use Waveform.from_sequence for optimal construction
-
-        :param subwaveforms: All waveforms must have the same defined channels
-        """
-        if not sub_waveforms:
-            raise ValueError(
-                "SequenceWaveform cannot be constructed without channel waveforms."
-            )
-
-        # do not fail on iterators although we do not allow them as an argument
-        sequenced_waveforms = tuple(sub_waveforms)
-
-        super().__init__(duration=sum(waveform.duration for waveform in sequenced_waveforms))
-        self._sequenced_waveforms = sequenced_waveforms
-
-        defined_channels = self._sequenced_waveforms[0].defined_channels
-        if not all(waveform.defined_channels == defined_channels
-                   for waveform in itertools.islice(self._sequenced_waveforms, 1, None)):
-            for waveform in self._sequenced_waveforms[1:]:
-                 if not waveform.defined_channels == self.defined_channels:
-                     print(f"SequenceWaveform: defined channels {self.defined_channels} do not match {waveform.defined_channels} ")
-            raise ValueError(
-                "SequenceWaveform cannot be constructed from waveforms of different"
-                "defined channels."
-            )
-
-    @classmethod
-    def from_sequence(cls, waveforms: Sequence['Waveform']) -> 'Waveform':
-        """Returns a waveform the represents the given sequence of waveforms. Applies some optimizations."""
-        assert waveforms, "Sequence must not be empty"
-        if len(waveforms) == 1:
-            return waveforms[0]
-
-        flattened = []
-        constant_values = waveforms[0].constant_value_dict()
-        for wf in waveforms:
-            if constant_values and constant_values != wf.constant_value_dict():
-                constant_values = None
-            if isinstance(wf, cls):
-                flattened.extend(wf.sequenced_waveforms)
-            else:
-                flattened.append(wf)
-        if constant_values is None:
-            return cls(sub_waveforms=flattened)
-        else:
-            duration = sum(wf.duration for wf in flattened)
-            return ConstantWaveform.from_mapping(duration, constant_values)
-
-    def is_constant(self) -> bool:
-        # only correct if from_sequence is used for construction
-        return False
-
-    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
-        # only correct if from_sequence is used for construction
-        return None
-
-    def constant_value(self, channel: ChannelID) -> Optional[float]:
-        v = None
-        for wf in self._sequenced_waveforms:
-            wf_cv = wf.constant_value(channel)
-            if wf_cv is None:
-                return None
-            elif wf_cv == v:
-                continue
-            elif v is None:
-                v = wf_cv
-            else:
-                assert v != wf_cv
-                return None
-        return v
-
-    @property
-    def defined_channels(self) -> AbstractSet[ChannelID]:
-        return self._sequenced_waveforms[0].defined_channels
-
-    def unsafe_sample(self,
-                      channel: ChannelID,
-                      sample_times: np.ndarray,
-                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
-        if output_array is None:
-            output_array = _ALLOCATION_FUNCTION(sample_times, **_ALLOCATION_FUNCTION_KWARGS)
-        time = 0
-        for subwaveform in self._sequenced_waveforms:
-            # before you change anything here, make sure to understand the difference between basic and advanced
-            # indexing in numpy and their copy/reference behaviour
-            end = time + subwaveform.duration
-
-            indices = slice(*np.searchsorted(sample_times, (float(time), float(end)), 'left'))
-            subwaveform.unsafe_sample(channel=channel,
-                                      sample_times=sample_times[indices]-np.float64(time),
-                                      output_array=output_array[indices])
-            time = end
-        return output_array
-
-    @property
-    def compare_key(self) -> Tuple[Waveform]:
-        return self._sequenced_waveforms
-
-    @property
-    def duration(self) -> TimeType:
-        return self._duration
-
-    def unsafe_get_subset_for_channels(self, channels: AbstractSet[ChannelID]) -> 'Waveform':
-        return SequenceWaveform.from_sequence([
-            sub_waveform.unsafe_get_subset_for_channels(channels & sub_waveform.defined_channels)
-            for sub_waveform in self._sequenced_waveforms if sub_waveform.defined_channels & channels])
-
-    @property
-    def sequenced_waveforms(self) -> Sequence[Waveform]:
-        return self._sequenced_waveforms
-
-    def __repr__(self):
-        return f"{type(self).__name__}({self._sequenced_waveforms})"
-
-
-class MultiChannelWaveform(Waveform):
-    """A MultiChannelWaveform is a Waveform object that allows combining arbitrary Waveform objects
-    to into a single waveform defined for several channels.
-
-    The number of channels used by the MultiChannelWaveform object is the sum of the channels used
-    by the Waveform objects it consists of.
-
-    MultiChannelWaveform allows an arbitrary mapping of channels defined by the Waveforms it
-    consists of and the channels it defines. For example, if the MultiChannelWaveform consists
-    of a two Waveform objects A and B which define two channels each, then the channels of the
-    MultiChannelWaveform may be 0: A.1, 1: B.0, 2: B.1, 3: A.0 where A.0 means channel 0 of Waveform
-    object A.
-
-    The following constraints must hold:
-     - The durations of all Waveform objects must be equal.
-     - The channel mapping must be sane, i.e., no channel of the MultiChannelWaveform must be
-        assigned more than one channel of any Waveform object it consists of
-    """
-
-    __slots__ = ('_sub_waveforms', '_defined_channels')
-
-    def __init__(self, sub_waveforms: List[Waveform]) -> None:
-        """Create a new MultiChannelWaveform instance.
-        Use `MultiChannelWaveform.from_parallel` for optimal construction.
-
-        Requires a list of subwaveforms in the form (Waveform, List(int)) where the list defines
-        the channel mapping, i.e., a value y at index x in the list means that channel x of the
-        subwaveform will be mapped to channel y of this MultiChannelWaveform object.
-
-        Args:
-            sub_waveforms: The list of sub waveforms of this
-                MultiChannelWaveform. List might get sorted!
-        Raises:
-            ValueError, if a channel mapping is out of bounds of the channels defined by this
-                MultiChannelWaveform
-            ValueError, if several subwaveform channels are assigned to a single channel of this
-                MultiChannelWaveform
-            ValueError, if subwaveforms have inconsistent durations
-        """
-
-        if not sub_waveforms:
-            raise ValueError(
-                "MultiChannelWaveform cannot be constructed without channel waveforms."
-            )
-
-        # sort the waveforms with their defined channels to make compare key reproducible
-        if not isinstance(sub_waveforms, list):
-            sub_waveforms = list(sub_waveforms)
-        sub_waveforms.sort(key=lambda wf: wf._sort_key_for_channels())
-
-        super().__init__(duration=sub_waveforms[0].duration)
-        self._sub_waveforms = tuple(sub_waveforms)
-
-        defined_channels = set()
-        for waveform in self._sub_waveforms:
-            if waveform.defined_channels & defined_channels:
-                raise ValueError('Channel may not be defined in multiple waveforms',
-                                 waveform.defined_channels & defined_channels)
-            defined_channels |= waveform.defined_channels
-        self._defined_channels = frozenset(defined_channels)
-
-        if not all(isclose(waveform.duration, self.duration) for waveform in self._sub_waveforms[1:]):
-            # meaningful error message:
-            durations = {}
-
-            for waveform in self._sub_waveforms:
-                for duration, channels in durations.items():
-                    if isclose(waveform.duration, duration):
-                        channels.update(waveform.defined_channels)
-                        break
-                else:
-                    durations[waveform.duration] = set(waveform.defined_channels)
-
-            raise ValueError(
-                "MultiChannelWaveform cannot be constructed from channel waveforms of different durations.",
-                durations
-            )
-
-    @staticmethod
-    def from_parallel(waveforms: Sequence[Waveform]) -> Waveform:
-        assert waveforms, "ARgument must not be empty"
-        if len(waveforms) == 1:
-            return waveforms[0]
-
-        # we do not look at constant values here because there is no benefit. We would need to construct a new
-        # MultiChannelWaveform anyways
-
-        # avoid unnecessary multi channel nesting
-        flattened = []
-        for waveform in waveforms:
-            if isinstance(waveform, MultiChannelWaveform):
-                flattened.extend(waveform._sub_waveforms)
-            else:
-                flattened.append(waveform)
-
-        return MultiChannelWaveform(flattened)
-
-    def is_constant(self) -> bool:
-        return all(wf.is_constant() for wf in self._sub_waveforms)
-
-    def constant_value(self, channel: ChannelID) -> Optional[float]:
-        return self[channel].constant_value(channel)
-
-    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
-        d = {}
-        for wf in self._sub_waveforms:
-            wf_d = wf.constant_value_dict()
-            if wf_d is None:
-                return None
-            else:
-                d.update(wf_d)
-        return d
-
-    @property
-    def duration(self) -> TimeType:
-        return self._sub_waveforms[0].duration
-
-    def __getitem__(self, key: ChannelID) -> Waveform:
-        for waveform in self._sub_waveforms:
-            if key in waveform.defined_channels:
-                return waveform
-        raise KeyError('Unknown channel ID: {}'.format(key), key)
-
-    @property
-    def defined_channels(self) -> AbstractSet[ChannelID]:
-        return self._defined_channels
-
-    @property
-    def compare_key(self) -> Any:
-        # sort with channels
-        return self._sub_waveforms
-
-    def unsafe_sample(self,
-                      channel: ChannelID,
-                      sample_times: np.ndarray,
-                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
-        return self[channel].unsafe_sample(channel, sample_times, output_array)
-
-    def unsafe_get_subset_for_channels(self, channels: AbstractSet[ChannelID]) -> 'Waveform':
-        relevant_sub_waveforms = [swf for swf in self._sub_waveforms if swf.defined_channels & channels]
-        if len(relevant_sub_waveforms) == 1:
-            return relevant_sub_waveforms[0].get_subset_for_channels(channels)
-        elif len(relevant_sub_waveforms) > 1:
-            return MultiChannelWaveform.from_parallel(
-                [sub_waveform.get_subset_for_channels(channels & sub_waveform.defined_channels)
-                 for sub_waveform in relevant_sub_waveforms])
-        else:
-            raise KeyError('Unknown channels: {}'.format(channels))
-
-    def __repr__(self):
-        return f"{type(self).__name__}({self._sub_waveforms!r})"
-
-
-class RepetitionWaveform(Waveform):
-    """This class allows putting multiple PulseTemplate together in one waveform on the hardware."""
-
-    __slots__ = ('_body', '_repetition_count')
-
-    def __init__(self, body: Waveform, repetition_count: int):
-        repetition_count = checked_int_cast(repetition_count)
-        if repetition_count < 1 or not isinstance(repetition_count, int):
-            raise ValueError('Repetition count must be an integer >0')
-
-        super().__init__(duration=body.duration * repetition_count)
-        self._body = body
-        self._repetition_count = repetition_count
-
-    @classmethod
-    def from_repetition_count(cls, body: Waveform, repetition_count: int) -> Waveform:
-        constant_values = body.constant_value_dict()
-        if constant_values is None:
-            return RepetitionWaveform(body, repetition_count)
-        else:
-            return ConstantWaveform.from_mapping(body.duration * repetition_count, constant_values)
-
-    @property
-    def defined_channels(self) -> AbstractSet[ChannelID]:
-        return self._body.defined_channels
-
-    def unsafe_sample(self,
-                      channel: ChannelID,
-                      sample_times: np.ndarray,
-                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
-        if output_array is None:
-            output_array = _ALLOCATION_FUNCTION(sample_times, **_ALLOCATION_FUNCTION_KWARGS)
-        body_duration = self._body.duration
-        time = 0
-        for _ in range(self._repetition_count):
-            end = time + body_duration
-            indices = slice(*np.searchsorted(sample_times, (float(time), float(end)), 'left'))
-            self._body.unsafe_sample(channel=channel,
-                                     sample_times=sample_times[indices] - float(time),
-                                     output_array=output_array[indices])
-            time = end
-        return output_array
-
-    @property
-    def compare_key(self) -> Tuple[Any, int]:
-        return self._body.compare_key, self._repetition_count
-
-    def unsafe_get_subset_for_channels(self, channels: AbstractSet[ChannelID]) -> Waveform:
-        return RepetitionWaveform.from_repetition_count(
-            body=self._body.unsafe_get_subset_for_channels(channels),
-            repetition_count=self._repetition_count)
-
-    def is_constant(self) -> bool:
-        return self._body.is_constant()
-
-    def constant_value(self, channel: ChannelID) -> Optional[float]:
-        return self._body.constant_value(channel)
-
-    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
-        return self._body.constant_value_dict()
-
-    def __repr__(self):
-        return f"{type(self).__name__}(body={self._body!r}, repetition_count={self._repetition_count!r})"
-
-
-class TransformingWaveform(Waveform):
-    __slots__ = ('_inner_waveform', '_transformation', '_cached_data', '_cached_times')
-
-    def __init__(self, inner_waveform: Waveform, transformation: Transformation):
-        """"""
-        super(TransformingWaveform, self).__init__(duration=inner_waveform.duration)
-        self._inner_waveform = inner_waveform
-        self._transformation = transformation
-
-        # cache data of inner channels based identified and invalidated by the sample times
-        self._cached_data = None
-        self._cached_times = lambda: None
-
-    @classmethod
-    def from_transformation(cls, inner_waveform: Waveform, transformation: Transformation) -> Waveform:
-        constant_values = inner_waveform.constant_value_dict()
-
-        if constant_values is None or not transformation.is_constant_invariant():
-            return cls(inner_waveform, transformation)
-
-        transformed_constant_values = {key: float(value) for key, value in transformation(0., constant_values).items()}
-        return ConstantWaveform.from_mapping(inner_waveform.duration, transformed_constant_values)
-
-    def is_constant(self) -> bool:
-        # only true if `from_transformation` was used
-        return False
-
-    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
-        # only true if `from_transformation` was used
-        return None
-
-    def constant_value(self, channel: ChannelID) -> Optional[float]:
-        if not self._transformation.is_constant_invariant():
-            return None
-        in_channels = self._transformation.get_input_channels({channel})
-        in_values = {ch: self._inner_waveform.constant_value(ch) for ch in in_channels}
-        if any(val is None for val in in_values.values()):
-            return None
-        else:
-            return self._transformation(0., in_values)[channel]
-
-    @property
-    def inner_waveform(self) -> Waveform:
-        return self._inner_waveform
-
-    @property
-    def transformation(self) -> Transformation:
-        return self._transformation
-
-    @property
-    def defined_channels(self) -> AbstractSet[ChannelID]:
-        return self.transformation.get_output_channels(self.inner_waveform.defined_channels)
-
-    @property
-    def compare_key(self) -> Tuple[Waveform, Transformation]:
-        return self.inner_waveform, self.transformation
-
-    def unsafe_get_subset_for_channels(self, channels: Set[ChannelID]) -> 'SubsetWaveform':
-        return SubsetWaveform(self, channel_subset=channels)
-
-    def unsafe_sample(self,
-                      channel: ChannelID,
-                      sample_times: np.ndarray,
-                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
-        if self._cached_times() is not sample_times:
-            self._cached_data = dict()
-            self._cached_times = ref(sample_times)
-
-        if channel not in self._cached_data:
-
-            inner_channels = self.transformation.get_input_channels({channel})
-
-            inner_data = {inner_channel: self.inner_waveform.unsafe_sample(inner_channel, sample_times)
-                          for inner_channel in inner_channels}
-
-            outer_data = self.transformation(sample_times, inner_data)
-            self._cached_data.update(outer_data)
-
-        if output_array is None:
-            output_array = self._cached_data[channel]
-        else:
-            output_array[:] = self._cached_data[channel]
-
-        return output_array
-
-
-class SubsetWaveform(Waveform):
-    __slots__ = ('_inner_waveform', '_channel_subset')
-
-    def __init__(self, inner_waveform: Waveform, channel_subset: Set[ChannelID]):
-        super().__init__(duration=inner_waveform.duration)
-        self._inner_waveform = inner_waveform
-        self._channel_subset = frozenset(channel_subset)
-
-    @property
-    def inner_waveform(self) -> Waveform:
-        return self._inner_waveform
-
-    @property
-    def defined_channels(self) -> FrozenSet[ChannelID]:
-        return self._channel_subset
-
-    @property
-    def compare_key(self) -> Tuple[frozenset, Waveform]:
-        return self.defined_channels, self.inner_waveform
-
-    def unsafe_get_subset_for_channels(self, channels: Set[ChannelID]) -> Waveform:
-        return self.inner_waveform.get_subset_for_channels(channels)
-
-    def unsafe_sample(self,
-                      channel: ChannelID,
-                      sample_times: np.ndarray,
-                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
-        return self.inner_waveform.unsafe_sample(channel, sample_times, output_array)
-
-    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
-        d = self._inner_waveform.constant_value_dict()
-        if d is not None:
-            return {ch: d[ch] for ch in self._channel_subset}
-
-    def constant_value(self, channel: ChannelID) -> Optional[float]:
-        if channel not in self._channel_subset:
-            raise KeyError(channel)
-        return self._inner_waveform.constant_value(channel)
-
-
-class ArithmeticWaveform(Waveform):
-    """Channels only present in one waveform have the operations neutral element on the other."""
-
-    numpy_operator_map = {'+': np.add,
-                          '-': np.subtract}
-    operator_map = {'+': operator.add,
-                    '-': operator.sub}
-
-    rhs_only_map = {'+': operator.pos,
-                    '-': operator.neg}
-    numpy_rhs_only_map = {'+': np.positive,
-                          '-': np.negative}
-
-    __slots__ = ('_lhs', '_rhs', '_arithmetic_operator')
-
-    def __init__(self,
-                 lhs: Waveform,
-                 arithmetic_operator: str,
-                 rhs: Waveform):
-        super().__init__(duration=lhs.duration)
-        self._lhs = lhs
-        self._rhs = rhs
-        self._arithmetic_operator = arithmetic_operator
-
-        assert np.isclose(float(self._lhs.duration), float(self._rhs.duration))
-        assert arithmetic_operator in self.operator_map
-
-    @classmethod
-    def from_operator(cls, lhs: Waveform, arithmetic_operator: str, rhs: Waveform):
-        # one could optimize rhs_cv to being only created if lhs_cv is not None but this makes the code harder to read
-        lhs_cv = lhs.constant_value_dict()
-        rhs_cv = rhs.constant_value_dict()
-        if lhs_cv is None or rhs_cv is None:
-            return cls(lhs, arithmetic_operator, rhs)
-
-        else:
-            constant_values = dict(lhs_cv)
-            op = cls.operator_map[arithmetic_operator]
-            rhs_op = cls.rhs_only_map[arithmetic_operator]
-
-            for ch, rhs_val in rhs_cv.items():
-                if ch in constant_values:
-                    constant_values[ch] = op(constant_values[ch], rhs_val)
-                else:
-                    constant_values[ch] = rhs_op(rhs_val)
-
-            duration = lhs.duration
-            assert isclose(duration, rhs.duration)
-
-            return ConstantWaveform.from_mapping(duration, constant_values)
-
-    def constant_value(self, channel: ChannelID) -> Optional[float]:
-        if channel not in self._rhs.defined_channels:
-            return self._lhs.constant_value(channel)
-        rhs = self._rhs.constant_value(channel)
-        if rhs is None:
-            return None
-
-        if channel in self._lhs.defined_channels:
-            lhs = self._lhs.constant_value(channel)
-            if lhs is None:
-                return None
-
-            return self.operator_map[self._arithmetic_operator](lhs, rhs)
-        else:
-            return self.rhs_only_map[self._arithmetic_operator](rhs)
-
-    def is_constant(self) -> bool:
-        # only correct if from_operator is used
-        return False
-
-    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
-        # only correct if from_operator is used
-        return None
-
-    @property
-    def lhs(self) -> Waveform:
-        return self._lhs
-
-    @property
-    def rhs(self) -> Waveform:
-        return self._rhs
-
-    @property
-    def arithmetic_operator(self) -> str:
-        return self._arithmetic_operator
-
-    @property
-    def duration(self) -> TimeType:
-        return self._lhs.duration
-
-    @property
-    def defined_channels(self) -> AbstractSet[ChannelID]:
-        return self._lhs.defined_channels | self._rhs.defined_channels
-
-    def unsafe_sample(self,
-                      channel: ChannelID,
-                      sample_times: np.ndarray,
-                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
-        if channel in self._lhs.defined_channels:
-            lhs = self._lhs.unsafe_sample(channel=channel, sample_times=sample_times, output_array=output_array)
-        else:
-            lhs = None
-
-        if channel in self._rhs.defined_channels:
-            rhs = self._rhs.unsafe_sample(channel=channel, sample_times=sample_times,
-                                          output_array=None if lhs is not None else output_array)
-        else:
-            rhs = None
-
-        if rhs is not None and lhs is not None:
-            arithmetic_operator = self.numpy_operator_map[self._arithmetic_operator]
-            if output_array is None:
-                output_array = lhs
-            return arithmetic_operator(lhs, rhs, out=output_array)
-
-        else:
-            if lhs is None:
-                assert rhs is not None, "channel %r not in defined channels (internal bug)" % channel
-                return self.numpy_rhs_only_map[self._arithmetic_operator](rhs, out=output_array)
-            else:
-                return lhs
-
-    def unsafe_get_subset_for_channels(self, channels: Set[ChannelID]) -> Waveform:
-        # TODO: optimization possible
-        return SubsetWaveform(self, channels)
-
-    @property
-    def compare_key(self) -> Tuple[str, Waveform, Waveform]:
-        return self._arithmetic_operator, self._lhs, self._rhs
-
-
-class FunctorWaveform(Waveform):
-    # TODO: Use Protocol to enforce that it accepts second argument has the keyword out
-    Functor = callable
-
-    __slots__ = ('_inner_waveform', '_functor')
-
-    """Apply a channel wise functor that works inplace to all results. The functor must accept two arguments"""
-    def __init__(self, inner_waveform: Waveform, functor: Mapping[ChannelID, Functor]):
-        super(FunctorWaveform, self).__init__(duration=inner_waveform.duration)
-        self._inner_waveform = inner_waveform
-        self._functor = dict(functor.items())
-
-        assert set(functor.keys()) == inner_waveform.defined_channels, ("There is no default identity mapping (yet)."
-                                                                        "File an issue on github if you need it.")
-
-    @classmethod
-    def from_functor(cls, inner_waveform: Waveform, functor: Mapping[ChannelID, Functor]):
-        constant_values = inner_waveform.constant_value_dict()
-        if constant_values is None:
-            return FunctorWaveform(inner_waveform, functor)
-
-        funced_constant_values = {ch: functor[ch](val) for ch, val in constant_values.items()}
-        return ConstantWaveform.from_mapping(inner_waveform.duration, funced_constant_values)
-
-    def is_constant(self) -> bool:
-        # only correct if `from_functor` was used
-        return False
-
-    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
-        # only correct if `from_functor` was used
-        return None
-
-    def constant_value(self, channel: ChannelID) -> Optional[float]:
-        inner = self._inner_waveform.constant_value(channel)
-        if inner is None:
-            return None
-        else:
-            return self._functor[channel](inner)
-
-    @property
-    def defined_channels(self) -> AbstractSet[ChannelID]:
-        return self._inner_waveform.defined_channels
-
-    def unsafe_sample(self,
-                      channel: ChannelID,
-                      sample_times: np.ndarray,
-                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
-        inner_output = self._inner_waveform.unsafe_sample(channel, sample_times, output_array)
-        return self._functor[channel](inner_output, out=inner_output)
-
-    def unsafe_get_subset_for_channels(self, channels: Set[ChannelID]) -> Waveform:
-        return FunctorWaveform.from_functor(
-            self._inner_waveform.unsafe_get_subset_for_channels(channels),
-            {ch: self._functor[ch] for ch in channels})
-
-    @property
-    def compare_key(self) -> Tuple[Waveform, FrozenSet]:
-        return self._inner_waveform, frozenset(self._functor.items())
-
-
-class ReversedWaveform(Waveform):
-    """Reverses the inner waveform in time."""
-
-    __slots__ = ('_inner',)
-
-    def __init__(self, inner: Waveform):
-        super().__init__(duration=inner.duration)
-        self._inner = inner
-
-    @classmethod
-    def from_to_reverse(cls, inner: Waveform) -> Waveform:
-        if inner.constant_value_dict():
-            return inner
-        else:
-            return cls(inner)
-
-    def unsafe_sample(self, channel: ChannelID, sample_times: np.ndarray,
-                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
-        inner_sample_times = (float(self.duration) - sample_times)[::-1]
-        if output_array is None:
-            return self._inner.unsafe_sample(channel, inner_sample_times, None)[::-1]
-        else:
-            inner_output_array = output_array[::-1]
-            inner_output_array = self._inner.unsafe_sample(channel, inner_sample_times, output_array=inner_output_array)
-            if inner_output_array.base not in (output_array, output_array.base):
-                # TODO: is there a guarantee by numpy we never end up here?
-                output_array[:] = inner_output_array[::-1]
-            return output_array
-
-    @property
-    def defined_channels(self) -> AbstractSet[ChannelID]:
-        return self._inner.defined_channels
-
-    def unsafe_get_subset_for_channels(self, channels: AbstractSet[ChannelID]) -> 'Waveform':
-        return ReversedWaveform.from_to_reverse(self._inner.unsafe_get_subset_for_channels(channels))
-
-    @property
-    def compare_key(self) -> Hashable:
-        return self._inner.compare_key
-
-    def reversed(self) -> 'Waveform':
-        return self._inner
+del qupulse
diff --git a/qupulse/hardware/awgs/base.py b/qupulse/hardware/awgs/base.py
index 49878803..10823017 100644
--- a/qupulse/hardware/awgs/base.py
+++ b/qupulse/hardware/awgs/base.py
@@ -14,8 +14,8 @@
 
 from qupulse.hardware.util import get_sample_times, not_none_indices
 from qupulse.utils.types import ChannelID
-from qupulse._program._loop import Loop
-from qupulse._program.waveforms import Waveform
+from qupulse.program.loop import Loop
+from qupulse.program.waveforms import Waveform
 from qupulse.comparable import Comparable
 from qupulse.utils.types import TimeType
 
diff --git a/qupulse/hardware/awgs/tabor.py b/qupulse/hardware/awgs/tabor.py
index d964a566..90e03866 100644
--- a/qupulse/hardware/awgs/tabor.py
+++ b/qupulse/hardware/awgs/tabor.py
@@ -11,7 +11,7 @@
 import numpy as np
 
 from qupulse.utils.types import ChannelID
-from qupulse._program._loop import Loop, make_compatible
+from qupulse.program.loop import Loop, make_compatible
 from qupulse.hardware.util import voltage_to_uint16, traced
 from qupulse.hardware.awgs.base import AWG, AWGAmplitudeOffsetHandling
 from qupulse._program.tabor import TaborSegment, TaborException, TaborProgram, PlottableProgram, TaborSequencing,\
diff --git a/qupulse/hardware/awgs/tektronix.py b/qupulse/hardware/awgs/tektronix.py
index 5865388e..eaaa7a36 100644
--- a/qupulse/hardware/awgs/tektronix.py
+++ b/qupulse/hardware/awgs/tektronix.py
@@ -11,8 +11,8 @@
 
 from qupulse.hardware.awgs.base import AWG, AWGAmplitudeOffsetHandling, ProgramOverwriteException
 from qupulse import ChannelID
-from qupulse._program._loop import Loop, make_compatible
-from qupulse._program.waveforms import Waveform as QuPulseWaveform
+from qupulse.program.loop import Loop, make_compatible
+from qupulse.program.waveforms import Waveform as QuPulseWaveform
 from qupulse.utils.types import TimeType
 from qupulse.hardware.util import voltage_to_uint16, get_sample_times, traced
 from qupulse.utils import pairwise
diff --git a/qupulse/hardware/awgs/zihdawg.py b/qupulse/hardware/awgs/zihdawg.py
index 6d6972a9..c6ec80f2 100644
--- a/qupulse/hardware/awgs/zihdawg.py
+++ b/qupulse/hardware/awgs/zihdawg.py
@@ -30,7 +30,7 @@
 import time
 
 from qupulse.utils.types import ChannelID, TimeType, time_from_float
-from qupulse._program._loop import Loop, make_compatible
+from qupulse.program.loop import Loop, make_compatible
 from qupulse._program.seqc import HDAWGProgramManager, UserRegister, WaveformFileSystem
 from qupulse.hardware.awgs.base import AWG, ChannelNotFoundException, AWGAmplitudeOffsetHandling
 from qupulse.hardware.util import traced
diff --git a/qupulse/hardware/feature_awg/channel_tuple_wrapper.py b/qupulse/hardware/feature_awg/channel_tuple_wrapper.py
index 10212153..a840ce9c 100644
--- a/qupulse/hardware/feature_awg/channel_tuple_wrapper.py
+++ b/qupulse/hardware/feature_awg/channel_tuple_wrapper.py
@@ -1,7 +1,7 @@
 from typing import Tuple, Optional, Callable, Set
 
 from qupulse import ChannelID
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 from qupulse.hardware.feature_awg.base import AWGChannelTuple
 from qupulse.hardware.feature_awg.features import ProgramManagement, VolatileParameters
 from qupulse.hardware.awgs.base import AWG
diff --git a/qupulse/hardware/feature_awg/features.py b/qupulse/hardware/feature_awg/features.py
index 1b82189c..29439f70 100644
--- a/qupulse/hardware/feature_awg/features.py
+++ b/qupulse/hardware/feature_awg/features.py
@@ -3,7 +3,7 @@
 from numbers import Real
 from enum import Enum
 
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 from qupulse.hardware.feature_awg.base import AWGDeviceFeature, AWGChannelFeature, AWGChannelTupleFeature,\
     AWGChannelTuple
 from qupulse.utils.types import ChannelID
diff --git a/qupulse/hardware/feature_awg/tabor.py b/qupulse/hardware/feature_awg/tabor.py
index 4a08ab8b..55493971 100644
--- a/qupulse/hardware/feature_awg/tabor.py
+++ b/qupulse/hardware/feature_awg/tabor.py
@@ -11,7 +11,7 @@
 import numpy as np
 
 from qupulse import ChannelID
-from qupulse._program._loop import Loop, make_compatible
+from qupulse.program.loop import Loop, make_compatible
 
 from qupulse.hardware.feature_awg.channel_tuple_wrapper import ChannelTupleAdapter
 from qupulse.hardware.feature_awg.features import ChannelSynchronization, AmplitudeOffsetHandling, VoltageRange, \
diff --git a/qupulse/hardware/setup.py b/qupulse/hardware/setup.py
index 64681574..d034e3b2 100644
--- a/qupulse/hardware/setup.py
+++ b/qupulse/hardware/setup.py
@@ -5,7 +5,7 @@
 
 from qupulse.hardware.awgs.base import AWG
 from qupulse.hardware.dacs import DAC
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 
 from qupulse.utils.types import ChannelID
 
diff --git a/qupulse/hardware/util.py b/qupulse/hardware/util.py
index f2c60dba..f7d17357 100644
--- a/qupulse/hardware/util.py
+++ b/qupulse/hardware/util.py
@@ -10,7 +10,7 @@ def traced(obj):
         """Noop traced that is used if autologging package is not available"""
         return obj
 
-from qupulse._program.waveforms import Waveform
+from qupulse.program.waveforms import Waveform
 from qupulse.utils.types import TimeType
 from qupulse.utils import pairwise
 
diff --git a/qupulse/plotting.py b/qupulse/plotting.py
index 06f4d022..9810fb6a 100644
--- a/qupulse/plotting.py
+++ b/qupulse/plotting.py
@@ -23,11 +23,11 @@
     # was deprecated in matplotlib 3.7, but we keep it around to allow this code to work with older versions
     get_cmap = plt.get_cmap
 
-from qupulse._program import waveforms
+from qupulse.program import waveforms
 from qupulse.utils.types import ChannelID, MeasurementWindow, has_type_interface
 from qupulse.pulses.pulse_template import PulseTemplate
-from qupulse._program.waveforms import Waveform
-from qupulse._program._loop import Loop, to_waveform
+from qupulse.program.waveforms import Waveform
+from qupulse.program.loop import Loop, to_waveform
 
 
 __all__ = ["render", "plot", "PlottingNotPossibleException"]
diff --git a/qupulse/program/__init__.py b/qupulse/program/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/qupulse/program/loop.py b/qupulse/program/loop.py
new file mode 100644
index 00000000..a1182676
--- /dev/null
+++ b/qupulse/program/loop.py
@@ -0,0 +1,725 @@
+import reprlib
+import warnings
+from collections import defaultdict
+from enum import Enum
+from typing import Union, Dict, Iterable, Tuple, cast, List, Optional, Generator, Mapping
+
+import numpy as np
+
+from qupulse.program.volatile import VolatileRepetitionCount, VolatileProperty
+from qupulse.program.waveforms import SequenceWaveform, RepetitionWaveform
+from qupulse.program.waveforms import Waveform, ConstantWaveform
+from qupulse.utils import is_integer
+from qupulse.utils.numeric import smallest_factor_ge
+from qupulse.utils.tree import Node
+from qupulse.utils.types import TimeType, MeasurementWindow
+
+__all__ = ['Loop', 'make_compatible', 'MakeCompatibleWarning', 'to_waveform']
+
+
+DurationStructure = Tuple[int, Union[TimeType, 'DurationStructure']]
+
+
+class Loop(Node):
+    MAX_REPR_SIZE = 2000
+    __slots__ = ('_waveform', '_measurements', '_repetition_definition', '_cached_body_duration')
+
+    """This class represents a initialized (sub-)program as a tree. Each Loop of a valid program has a repetition count
+    and either a waveform or a sequence of loops as children.
+    
+    A Loop can have associated measurements which are also repeated.
+    """
+    def __init__(self,
+                 parent: Union['Loop', None] = None,
+                 children: Iterable['Loop'] = (),
+                 waveform: Optional[Waveform] = None,
+                 measurements: Optional[List[MeasurementWindow]] = None,
+                 repetition_count: Union[int, VolatileRepetitionCount] = 1):
+        """Initialize a new loop
+
+        Args:
+            parent: Forwarded to Node.__init__
+            children: Forwarded to Node.__init__
+            waveform: "Payload"
+            measurements: Associated measurements
+            repetition_count: The children / waveform are repeated this often
+        """
+        super().__init__(parent=parent, children=children)
+
+        self._waveform = waveform
+        self._measurements = measurements
+        self._repetition_definition = repetition_count
+        self._cached_body_duration = None
+        assert isinstance(repetition_count, VolatileRepetitionCount) or is_integer(repetition_count)
+        assert isinstance(waveform, (type(None), Waveform))
+
+    def __eq__(self, other: 'Loop') -> bool:
+        if type(self) == type(other):
+            return (self._repetition_definition == other._repetition_definition and
+                    self.waveform == other.waveform and
+                    (self._measurements or None) == (other._measurements or None) and
+                    len(self) == len(other) and
+                    all(self_child == other_child for self_child, other_child in zip(self, other)))
+        else:
+            return NotImplemented
+
+    def append_child(self, loop: Optional['Loop'] = None, **kwargs) -> None:
+        """Append a child to this loop. Either an existing Loop object or a newly created from kwargs
+
+        Args:
+            loop: loop to append
+            **kwargs: Child is constructed with these kwargs
+
+        Raises:
+            ValueError: if called with loop and kwargs
+        """
+        if loop is not None:
+            if kwargs:
+                raise ValueError("Cannot pass a Loop object and Loop constructor arguments at the same time in "
+                                 "append_child")
+            arg = (loop,)
+        else:
+            arg = (kwargs,)
+        super().__setitem__(slice(len(self), len(self)), arg)
+        self._invalidate_duration(body_duration_increment=self[-1].duration)
+
+    def _invalidate_duration(self, body_duration_increment=None):
+        if self._cached_body_duration is not None:
+            if body_duration_increment is not None:
+                self._cached_body_duration += body_duration_increment
+            else:
+                self._cached_body_duration = None
+        if self.parent:
+            if body_duration_increment is not None:
+                self.parent._invalidate_duration(body_duration_increment=body_duration_increment*self.repetition_count)
+            else:
+                self.parent._invalidate_duration()
+
+    def add_measurements(self, measurements: Iterable[MeasurementWindow]):
+        """Add measurements offset by the current body duration i.e. to the END of the current loop
+
+        Args:
+            measurements: Measurements to add
+        """
+        body_duration = float(self.body_duration)
+        if body_duration == 0:
+            measurements = measurements
+        else:
+            measurements = ((mw_name, begin+body_duration, length) for mw_name, begin, length in measurements)
+
+        if self._measurements is None:
+            self._measurements = list(measurements)
+        else:
+            self._measurements.extend(measurements)
+
+    @property
+    def waveform(self) -> Waveform:
+        return self._waveform
+
+    @waveform.setter
+    def waveform(self, val) -> None:
+        self._waveform = val
+        self._invalidate_duration()
+
+    @property
+    def body_duration(self) -> TimeType:
+        if self._cached_body_duration is None:
+            if self.is_leaf():
+                if self.waveform:
+                    self._cached_body_duration = self.waveform.duration
+                else:
+                    self._cached_body_duration = TimeType.from_fraction(0, 1)
+            else:
+                self._cached_body_duration = sum(child.duration for child in self)
+        return self._cached_body_duration
+
+    @property
+    def duration(self) -> TimeType:
+        return self.body_duration * TimeType.from_fraction(self.repetition_count, 1)
+
+    @property
+    def volatile_repetition(self) -> Optional[VolatileProperty]:
+        return getattr(self._repetition_definition, 'volatile_property', None)
+
+    @property
+    def repetition_definition(self) -> Union[int, VolatileRepetitionCount]:
+        return self._repetition_definition
+
+    @repetition_definition.setter
+    def repetition_definition(self, new_definition: Union[int, VolatileRepetitionCount]):
+        self._repetition_definition = new_definition
+
+    @property
+    def repetition_count(self) -> int:
+        return int(self._repetition_definition)
+
+    @repetition_count.setter
+    def repetition_count(self, val: int) -> None:
+        assert isinstance(val, (int, float))
+        new_repetition = int(val)
+        if abs(new_repetition - val) > 1e-10:
+            raise ValueError('Repetition count was not an integer')
+        self._repetition_definition = new_repetition
+
+    def unroll(self) -> None:
+        if self.is_leaf():
+            raise RuntimeError('Leaves cannot be unrolled')
+        if self.volatile_repetition:
+            warnings.warn("Unrolling a Loop with volatile repetition count", VolatileModificationWarning)
+
+        i = self.parent_index
+        self.parent[i:i+1] = (child.copy_tree_structure(new_parent=self.parent)
+                              for _ in range(self.repetition_count)
+                              for child in self)
+        self.parent.assert_tree_integrity()
+
+    def __setitem__(self, idx, value):
+        super().__setitem__(idx, value)
+        self._invalidate_duration()
+
+    def unroll_children(self) -> None:
+        if self.volatile_repetition:
+            warnings.warn("Unrolling a Loop with volatile repetition count", VolatileModificationWarning)
+        old_children = self.children
+        self[:] = (child.copy_tree_structure()
+                   for _ in range(self.repetition_count)
+                   for child in old_children)
+        self.repetition_count = 1
+        self.assert_tree_integrity()
+
+    def encapsulate(self) -> None:
+        """Add a nesting level by moving self to its children."""
+        self[:] = [Loop(children=self,
+                        repetition_count=self._repetition_definition,
+                        waveform=self._waveform,
+                        measurements=self._measurements)]
+        self.repetition_count = 1
+        self._waveform = None
+        self._measurements = None
+        self.assert_tree_integrity()
+
+    def _get_str(self, first_prefix, other_prefixes) -> Generator[str, None, None]:
+        if self.is_leaf():
+            yield '%sEXEC %r %d times' % (first_prefix, self._waveform, self.repetition_count)
+        else:
+            yield '%sLOOP %d times:' % (first_prefix, self.repetition_count)
+
+            for elem in self:
+                yield from cast(Loop, elem)._get_str(other_prefixes + '  ->', other_prefixes + '    ')
+
+    @reprlib.recursive_repr()
+    def __str__(self) -> str:
+        str_len = 0
+        repr_list = []
+        for sub_repr in self._get_str('', ''):
+            str_len += len(sub_repr)
+
+            if self.MAX_REPR_SIZE and str_len > self.MAX_REPR_SIZE:
+                repr_list.append('...')
+                break
+            else:
+                repr_list.append(sub_repr)
+        return '\n'.join(repr_list)
+
+    @reprlib.recursive_repr()
+    def __repr__(self):
+        children = self.children
+        waveform = self.waveform
+        measurements = self._measurements
+        repetition_count = self.repetition_count
+
+        max_repr_size = self.MAX_REPR_SIZE
+
+        reprs = {}
+        if children:
+            if len(children) < max_repr_size // len('Loop(...)'):
+                reprs['children'] = repr(list(children))
+            else:
+                reprs['children'] = '[...]'
+
+        if waveform:
+            waveform_repr = repr(waveform)
+            if len(waveform_repr) >= max_repr_size:
+                waveform_repr = '...'
+            reprs['waveform'] = waveform_repr
+
+        if measurements:
+            meas_repr = repr(measurements)
+            if len(meas_repr) >= max_repr_size:
+                meas_repr = '[...]'
+            reprs['measurements'] = meas_repr
+
+        if repetition_count != 1:
+            reprs['repetition_count'] = repr(repetition_count)
+
+        kwargs = ', '.join(f'{attr}={val}' for attr, val in reprs.items())
+
+        type_name = type(self).__name__
+        max_kwargs = max(self.MAX_REPR_SIZE - len(type_name) - 2, 0)
+
+        if len(kwargs) > max_kwargs:
+            kwargs = '...'
+
+        return f'{type(self).__name__}({kwargs})'
+
+    def copy_tree_structure(self, new_parent: Union['Loop', bool]=False) -> 'Loop':
+        return type(self)(parent=self.parent if new_parent is False else new_parent,
+                          waveform=self._waveform,
+                          repetition_count=self._repetition_definition,
+                          measurements=None if self._measurements is None else list(self._measurements),
+                          children=(child.copy_tree_structure() for child in self))
+
+    def _get_measurement_windows(self, drop: bool) -> Mapping[str, np.ndarray]:
+        """Private implementation of get_measurement_windows with a slightly different data format for easier tiling.
+
+        Args:
+            drop: Drops the measurements from the Loop i.e. the Loop will no longer have measurements attached after
+            collecting them
+
+        Returns:
+             A dictionary (measurement_name -> array) with begin == array[:, 0] and length == array[:, 1]
+        """
+        temp_meas_windows = defaultdict(list)
+        if self._measurements:
+            for (mw_name, begin, length) in self._measurements:
+                temp_meas_windows[mw_name].append((begin, length))
+
+            for mw_name, begin_length_list in temp_meas_windows.items():
+                temp_meas_windows[mw_name] = [np.asarray(begin_length_list, dtype=float)]
+
+            if drop:
+                self._measurements = None
+
+        # calculate duration together with meas windows in the same iteration
+        if self.is_leaf():
+            body_duration = float(self.body_duration)
+        else:
+            offset = TimeType(0)
+            for child in self:
+                for mw_name, begins_length_array in child._get_measurement_windows(drop).items():
+                    begins_length_array[:, 0] += float(offset)
+                    temp_meas_windows[mw_name].append(begins_length_array)
+                offset += child.duration
+
+            body_duration = float(offset)
+
+        # formatting like this for easier debugging
+        result = {}
+
+        # repeat and add repetition based offset
+        for mw_name, begin_length_list in temp_meas_windows.items():
+            result[mw_name] = _repeat_loop_measurements(begin_length_list, self.repetition_count, body_duration)
+
+        return result
+
+    def get_measurement_windows(self, drop=False) -> Dict[str, Tuple[np.ndarray, np.ndarray]]:
+        """Iterates over all children and collect the begin and length arrays of each measurement window.
+
+        Args:
+            drop: Drops the measurements from the Loop i.e. the Loop will no longer have measurements attached after
+            collecting them.
+
+        Returns:
+            A dictionary (measurement_name -> (begin, length)) with begin and length being :class:`numpy.ndarray`
+        """
+        return {mw_name: (begin_length_list[:, 0], begin_length_list[:, 1])
+                for mw_name, begin_length_list in self._get_measurement_windows(drop=drop).items()}
+
+    def split_one_child(self, child_index=None) -> None:
+        """Take the last child that has a repetition count larger one, decrease it's repetition count and insert a copy
+        with repetition cout one after it"""
+        if child_index is not None:
+            if self[child_index].repetition_count < 2:
+                raise ValueError('Cannot split child {} as the repetition count is not larger 1')
+
+        else:
+            # we cannot reverse enumerate
+            n_child = len(self) - 1
+            for reverse_idx, child in enumerate(reversed(self)):
+                if child.repetition_count > 1:
+                    forward_idx = n_child - reverse_idx
+                    if not child.volatile_repetition:
+                        child_index = forward_idx
+                        break
+                    elif child_index is None:
+                        child_index = forward_idx
+            else:
+                if child_index is None:
+                    raise RuntimeError('There is no child with repetition count > 1')
+
+        if self[child_index].volatile_repetition:
+            warnings.warn("Splitting a child with volatile repetition count", VolatileModificationWarning)
+
+        new_child = self[child_index].copy_tree_structure()
+        new_child.repetition_count = 1
+
+        self[child_index].repetition_count -= 1
+
+        self[child_index+1:child_index+1] = (new_child,)
+        self.assert_tree_integrity()
+
+    def flatten_and_balance(self, depth: int) -> None:
+        """Modifies the program so all tree branches have the same depth.
+
+        Args:
+            depth: Target depth of the program
+        """
+        i = 0
+        while i < len(self):
+            # only used by type checker
+            sub_program = cast(Loop, self[i])
+
+            if sub_program.depth() < depth - 1:
+                # increase nesting because the subprogram is not deep enough
+                sub_program.encapsulate()
+
+            elif not sub_program.is_balanced():
+                # balance the sub program. We revisit it in the next iteration (no change of i )
+                # because it might modify self. While writing this comment I am not sure this is true. 14.01.2020 Simon
+                sub_program.flatten_and_balance(depth - 1)
+
+            elif sub_program.depth() == depth - 1:
+                # subprogram is balanced with the correct depth
+                i += 1
+
+            elif sub_program._has_single_child_that_can_be_merged():
+                # subprogram is balanced but to deep and has no measurements -> we can "lift" the sub-sub-program
+                # TODO: There was a len(sub_sub_program) == 1 check here that I cannot explain
+                sub_program._merge_single_child()
+
+            elif not sub_program.is_leaf():
+                # subprogram is balanced but too deep
+                sub_program.unroll()
+
+            else:
+                # we land in this case if the function gets called with depth == 0 and the current subprogram is a leaf
+                i += 1
+
+    def _has_single_child_that_can_be_merged(self) -> bool:
+        """Check if self has only once child which can cheaply be merged into self by multiplying the repetition counts.
+        """
+        if len(self) == 1:
+            child = cast(Loop, self[0])
+            return not self._measurements or (child.repetition_count == 1 and not child.volatile_repetition)
+        else:
+            return False
+
+    def _merge_single_child(self):
+        """Lift the single child to current level. Requires _has_single_child_that_can_be_merged to be true"""
+        assert len(self) == 1, "bug: _merge_single_child called on loop with len != 1"
+        child = cast(Loop, self[0])
+
+        # if the child has a fixed repetition count of 1 the measurements can be merged
+        mergable_measurements = child.repetition_count == 1 and not child.volatile_repetition
+
+        assert not self._measurements or mergable_measurements, "bug: _merge_single_child called on loop with measurements"
+        assert not self._waveform, "bug: _merge_single_child called on loop with children and waveform"
+
+        measurements = child._measurements
+        if self._measurements:
+            if measurements:
+                measurements.extend(self._measurements)
+            else:
+                measurements = self._measurements
+
+        if not self.volatile_repetition and not child.volatile_repetition:
+            # simple integer multiplication
+            repetition_definition = self.repetition_count * child.repetition_count
+        elif not self.volatile_repetition:
+            repetition_definition = child._repetition_definition * self.repetition_count
+        elif not child.volatile_repetition:
+            repetition_definition = self._repetition_definition * child.repetition_count
+        else:
+            # create a new expression that depends on both
+            expression = 'parent_repetition_count * child_repetition_count'
+            repetition_definition = VolatileRepetitionCount.operation(
+                expression=expression,
+                parent_repetition_count=self._repetition_definition,
+                child_repetition_count=child._repetition_definition)
+
+        self[:] = iter(child)
+        self._waveform = child._waveform
+        self._repetition_definition = repetition_definition
+        self._measurements = measurements
+        self._invalidate_duration()
+        return True
+
+    def cleanup(self, actions=('remove_empty_loops', 'merge_single_child')):
+        """Apply the specified actions to cleanup the Loop.
+
+        remove_empty_loops: Remove loops with no children and no waveform (a DroppedMeasurementWarning is issued)
+        merge_single_child: see `_try_merge_single_child` documentation
+
+        Warnings:
+            DroppedMeasurementWarning: Likely a bug in qupulse. TODO: investigate whether there are usecases
+        """
+        if 'remove_empty_loops' in actions:
+            new_children = []
+            for child in self:
+                child = cast(Loop, child)
+                if child.is_leaf():
+                    if child.waveform is None:
+                        if child._measurements:
+                            warnings.warn("Dropping measurement since there is no waveform attached",
+                                          category=DroppedMeasurementWarning)
+                    else:
+                        new_children.append(child)
+
+                else:
+                    child.cleanup(actions)
+                    if child.waveform or not child.is_leaf():
+                        new_children.append(child)
+
+                    elif child._measurements:
+                        warnings.warn("Dropping measurement since there is no waveform in children",
+                                      category=DroppedMeasurementWarning)
+
+            if len(self) != len(new_children):
+                self[:] = new_children
+
+        else:
+            # only do the recursive call
+            for child in self:
+                child.cleanup(actions)
+
+        if 'merge_single_child' in actions and self._has_single_child_that_can_be_merged():
+            self._merge_single_child()
+    
+    def get_duration_structure(self) -> DurationStructure:
+        """Returns a tuple that fingerprints the structure of waveform durations and repetitions of self.
+
+        One possible use case is to identify repeated duration structures and reuse the same control flow with
+        differing data.
+        """
+        if self.is_leaf():
+            return self.repetition_count, self.waveform.duration
+        else:
+            return self.repetition_count, tuple(child.get_duration_structure() for child in self)
+
+    def reverse_inplace(self):
+        if self.is_leaf():
+            self._waveform = self._waveform.reversed()
+        else:
+            self._reverse_children()
+            for child in self:
+                child.reverse_inplace()
+        if self._measurements:
+            duration = self.duration
+            self._measurements = [
+                (name, duration - (begin + length), length)
+                for name, begin, length in self._measurements
+            ]
+
+
+def to_waveform(program: Loop) -> Waveform:
+    if program.is_leaf():
+        if program.repetition_count == 1:
+            return program.waveform
+        else:
+            return RepetitionWaveform.from_repetition_count(program.waveform, program.repetition_count)
+    else:
+        if len(program) == 1:
+            sequenced_waveform = to_waveform(cast(Loop, program[0]))
+        else:
+            sequenced_waveform = SequenceWaveform.from_sequence(
+                [to_waveform(cast(Loop, sub_program))
+                 for sub_program in program])
+        if program.repetition_count > 1:
+            return RepetitionWaveform.from_repetition_count(sequenced_waveform, program.repetition_count)
+        else:
+            return sequenced_waveform
+
+
+class _CompatibilityLevel(Enum):
+    compatible = 0
+    action_required = 1
+    incompatible_too_short = 2
+    incompatible_fraction = 3
+    incompatible_quantum = 4
+
+    def is_incompatible(self) -> bool:
+        return self in (self.incompatible_fraction, self.incompatible_quantum, self.incompatible_too_short)
+
+
+def _is_compatible(program: Loop, min_len: int, quantum: int, sample_rate: TimeType) -> _CompatibilityLevel:
+    """ check whether program loop is compatible with awg requirements
+        possible reasons for incompatibility:
+            program shorter than minimum length
+            program duration not an integer
+            program duration not a multiple of quantum """
+    program_duration_in_samples = program.duration * sample_rate
+
+    if program_duration_in_samples.denominator != 1:
+        return _CompatibilityLevel.incompatible_fraction
+
+    if program_duration_in_samples < min_len:
+        return _CompatibilityLevel.incompatible_too_short
+
+    if program_duration_in_samples % quantum > 0:
+        return _CompatibilityLevel.incompatible_quantum
+
+    if program.is_leaf():
+        waveform_duration_in_samples = program.body_duration * sample_rate
+        if waveform_duration_in_samples < min_len or (waveform_duration_in_samples / quantum).denominator != 1:
+            if program.volatile_repetition:
+                warnings.warn("_is_compatible requires an action which drops volatility.",
+                              category=VolatileModificationWarning)
+            return _CompatibilityLevel.action_required
+        else:
+            return _CompatibilityLevel.compatible
+    else:
+        if all(_is_compatible(cast(Loop, sub_program), min_len, quantum, sample_rate) == _CompatibilityLevel.compatible
+               for sub_program in program):
+            return _CompatibilityLevel.compatible
+        else:
+            if program.volatile_repetition:
+                warnings.warn("_is_compatible requires an action which drops volatility.",
+                              category=VolatileModificationWarning)
+            return _CompatibilityLevel.action_required
+
+
+def _make_compatible(program: Loop, min_len: int, quantum: int, sample_rate: TimeType) -> None:
+    if program.is_leaf():
+        program.waveform = to_waveform(program.copy_tree_structure())
+        program.repetition_count = 1
+    else:
+        comp_levels = [_is_compatible(cast(Loop, sub_program), min_len, quantum, sample_rate)
+                       for sub_program in program]
+
+        if any(comp_level.is_incompatible() for comp_level in comp_levels):
+            single_run = program.duration * sample_rate / program.repetition_count
+            if (single_run / quantum).denominator == 1 and single_run >= min_len:
+                # it is enough to concatenate all children
+                new_repetition_definition = program.repetition_definition
+                program.repetition_count = 1
+            else:
+                # we need to concatenate all children and unroll
+                new_repetition_definition = 1
+
+            program.waveform = to_waveform(program.copy_tree_structure())
+            program.repetition_definition = new_repetition_definition
+            program[:] = []
+            return
+        else:
+            for sub_program, comp_level in zip(program, comp_levels):
+                if comp_level == _CompatibilityLevel.action_required:
+                    _make_compatible(sub_program, min_len, quantum, sample_rate)
+
+
+def make_compatible(program: Loop, minimal_waveform_length: int, waveform_quantum: int, sample_rate: TimeType):
+    """ check program for compatibility to AWG requirements, make it compatible if necessary and  possible"""
+    comp_level = _is_compatible(program,
+                                min_len=minimal_waveform_length,
+                                quantum=waveform_quantum,
+                                sample_rate=sample_rate)
+    if comp_level == _CompatibilityLevel.incompatible_fraction:
+        raise ValueError('The program duration in samples {} is not an integer'.format(program.duration * sample_rate))
+    if comp_level == _CompatibilityLevel.incompatible_too_short:
+        raise ValueError('The program is too short to be a valid waveform. \n'
+                         ' program duration in samples: {} \n'
+                         ' minimal length: {}'.format(program.duration * sample_rate, minimal_waveform_length))
+    if comp_level == _CompatibilityLevel.incompatible_quantum:
+        raise ValueError('The program duration in samples {} '
+                         'is not a multiple of quantum {}'.format(program.duration * sample_rate, waveform_quantum))
+
+    elif comp_level == _CompatibilityLevel.action_required:
+        warnings.warn("qupulse will now concatenate waveforms to make the pulse/program compatible with the chosen AWG."
+                      " This might take some time. If you need this pulse more often it makes sense to write it in a "
+                      "way which is more AWG friendly.", MakeCompatibleWarning)
+
+        _make_compatible(program,
+                         min_len=minimal_waveform_length,
+                         quantum=waveform_quantum,
+                         sample_rate=sample_rate)
+
+    else:
+        assert comp_level == _CompatibilityLevel.compatible
+
+
+def roll_constant_waveforms(program: Loop, minimal_waveform_quanta: int, waveform_quantum: int, sample_rate: TimeType):
+    """This function finds waveforms in program that can be replaced with repetitions of shorter waveforms and replaces
+    them. Complexity O(N_waveforms). Drops measurements because they are not correctly handled here for performance
+    reasons.
+
+    This is possible if:
+     - The waveform is constant on all channels
+     - waveform.duration * sample_rate / waveform_quantum has a factor that is bigger than minimal_waveform_quanta
+
+    Args:
+        program:
+        minimal_waveform_quanta:
+        waveform_quantum:
+        sample_rate:
+
+    Warnings:
+        DroppedMeasurementWarning: This warning is raised if a measurement is dropped.
+    """
+    if program._measurements:
+        warnings.warn("Dropping measurements. Remove measurements before calling roll_constant_waveforms by calling"
+                      " get_measurement_windows(drop=True).", category=DroppedMeasurementWarning)
+        program._measurements = None
+
+    waveform = program.waveform
+
+    if waveform is None:
+        for child in program:
+            roll_constant_waveforms(child, minimal_waveform_quanta, waveform_quantum, sample_rate)
+    else:
+        waveform_quanta = (waveform.duration * sample_rate) // waveform_quantum
+
+        # example
+        # waveform_quanta = 15
+        # minimal_waveform_quanta = 2
+        # => repetition_count = 5, new_waveform_quanta = 3
+        if waveform_quanta < minimal_waveform_quanta * 2:
+            # there is no way to roll this waveform because it is too short
+            return
+
+        const_values = waveform.constant_value_dict()
+        if const_values is None:
+            # The waveform is not constant
+            return
+
+        new_waveform_quanta = smallest_factor_ge(waveform_quanta, min_factor=minimal_waveform_quanta)
+        if new_waveform_quanta == waveform_quanta:
+            # the waveform duration in samples has no suitable factor
+            # TODO: Option to insert multiple Loop objects
+            return
+
+        additional_repetition_count = waveform_quanta // new_waveform_quanta
+
+        new_waveform = ConstantWaveform.from_mapping(
+            duration=waveform_quantum * new_waveform_quanta / sample_rate,
+            constant_values=const_values)
+
+        # use the private properties to avoid invalidating the duration cache of the parent loop
+        program._repetition_definition = program.repetition_definition * additional_repetition_count
+        program._waveform = new_waveform
+
+
+def _repeat_loop_measurements(begin_length_list: List[np.ndarray],
+                              repetition_count: int,
+                              body_duration: float
+                              ) -> np.ndarray:
+    temp_begin_length_array = np.concatenate(begin_length_list)
+
+    begin_length_array = np.tile(temp_begin_length_array, (repetition_count, 1))
+
+    shaped_begin_length_array = np.reshape(begin_length_array, (repetition_count, -1, 2))
+
+    shaped_begin_length_array[:, :, 0] += (np.arange(repetition_count) * body_duration)[:, np.newaxis]
+
+    return begin_length_array
+
+
+class MakeCompatibleWarning(ResourceWarning):
+    pass
+
+
+class VolatileModificationWarning(RuntimeWarning):
+    """This warning is emitted if the colatile part of a program gets modified. This might imply that the volatile
+    parameter cannot be change anymore."""
+
+
+class DroppedMeasurementWarning(RuntimeWarning):
+    """This warning is emitted if a measurement was dropped because there was no waveform attached."""
diff --git a/qupulse/program/transformation.py b/qupulse/program/transformation.py
new file mode 100644
index 00000000..784e8e19
--- /dev/null
+++ b/qupulse/program/transformation.py
@@ -0,0 +1,406 @@
+from typing import Any, Mapping, Set, Tuple, Sequence, AbstractSet, Union, TYPE_CHECKING, Hashable
+from abc import abstractmethod
+from numbers import Real
+
+import numpy as np
+
+from qupulse import ChannelID
+from qupulse.comparable import Comparable
+from qupulse.utils.types import SingletonABCMeta, frozendict
+from qupulse.expressions import ExpressionScalar
+
+
+_TrafoValue = Union[Real, ExpressionScalar]
+
+
+__all__ = ['Transformation', 'IdentityTransformation', 'LinearTransformation', 'ScalingTransformation',
+           'OffsetTransformation', 'ParallelChannelTransformation', 'ChainedTransformation',
+           'chain_transformations']
+
+
+class Transformation(Comparable):
+    _identity_singleton = None
+    """Transforms numeric time-voltage values for multiple channels to other time-voltage values. The number and names
+     of input and output channels might differ."""
+
+    @abstractmethod
+    def __call__(self, time: Union[np.ndarray, float],
+                 data: Mapping[ChannelID, Union[np.ndarray, float]]) -> Mapping[ChannelID, Union[np.ndarray, float]]:
+        """Apply transformation to data
+        Args:
+            time:
+            data:
+
+        Returns:
+            transformed: A DataFrame that has been transformed with index == output_channels
+        """
+
+    @abstractmethod
+    def get_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        """Return the channel identifiers"""
+
+    @abstractmethod
+    def get_input_channels(self, output_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        """Channels that are required for getting data for the requested output channel"""
+
+    def chain(self, next_transformation: 'Transformation') -> 'Transformation':
+        if next_transformation is IdentityTransformation():
+            return self
+        else:
+            return chain_transformations(self, next_transformation)
+
+    def is_constant_invariant(self):
+        """Signals if the transformation always maps constants to constants."""
+        return False
+
+    def get_constant_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        return frozenset()
+
+
+class IdentityTransformation(Transformation, metaclass=SingletonABCMeta):
+    def __call__(self, time: Union[np.ndarray, float],
+                 data: Mapping[ChannelID, Union[np.ndarray, float]]) -> Mapping[ChannelID, Union[np.ndarray, float]]:
+        return data
+
+    def get_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        return input_channels
+
+    @property
+    def compare_key(self) -> None:
+        return None
+
+    def get_input_channels(self, output_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        return output_channels
+
+    def chain(self, next_transformation: Transformation) -> Transformation:
+        return next_transformation
+
+    def __repr__(self):
+        return 'IdentityTransformation()'
+
+    def is_constant_invariant(self):
+        """Signals if the transformation always maps constants to constants."""
+        return True
+
+    def get_constant_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        return input_channels
+
+
+class ChainedTransformation(Transformation):
+    def __init__(self, *transformations: Transformation):
+        self._transformations = transformations
+
+    @property
+    def transformations(self) -> Tuple[Transformation, ...]:
+        return self._transformations
+
+    def get_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        for transformation in self._transformations:
+            input_channels = transformation.get_output_channels(input_channels)
+        return input_channels
+
+    def get_input_channels(self, output_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        for transformation in reversed(self._transformations):
+            output_channels = transformation.get_input_channels(output_channels)
+        return output_channels
+
+    def __call__(self, time: Union[np.ndarray, float],
+                 data: Mapping[ChannelID, Union[np.ndarray, float]]) -> Mapping[ChannelID, Union[np.ndarray, float]]:
+        for transformation in self._transformations:
+            data = transformation(time, data)
+        return data
+
+    @property
+    def compare_key(self) -> Tuple[Transformation, ...]:
+        return self._transformations
+
+    def chain(self, next_transformation) -> Transformation:
+        return chain_transformations(*self.transformations, next_transformation)
+
+    def __repr__(self):
+        return f'{type(self).__name__}{self._transformations!r}'
+
+    def is_constant_invariant(self):
+        """Signals if the transformation always maps constants to constants."""
+        return all(trafo.is_constant_invariant() for trafo in self._transformations)
+
+    def get_constant_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        for trafo in self._transformations:
+            input_channels = trafo.get_constant_output_channels(input_channels)
+        return input_channels
+
+
+class LinearTransformation(Transformation):
+    def __init__(self,
+                 transformation_matrix: np.ndarray,
+                 input_channels: Sequence[ChannelID],
+                 output_channels: Sequence[ChannelID]):
+        """
+
+        Args:
+            transformation_matrix: Matrix describing the transformation with shape (output_channels, input_channels)
+            input_channels: Channel ids of the columns
+            output_channels: Channel ids of the rows
+        """
+        transformation_matrix = np.asarray(transformation_matrix)
+
+        if transformation_matrix.shape != (len(output_channels), len(input_channels)):
+            raise ValueError('Shape of transformation matrix does not match to the given channels')
+
+        output_sorter = np.argsort(output_channels)
+        transformation_matrix = transformation_matrix[output_sorter, :]
+
+        input_sorter = np.argsort(input_channels)
+        transformation_matrix = transformation_matrix[:, input_sorter]
+
+        self._matrix = transformation_matrix
+        self._input_channels = tuple(sorted(input_channels))
+        self._output_channels = tuple(sorted(output_channels))
+        self._input_channels_set = frozenset(self._input_channels)
+        self._output_channels_set = frozenset(self._output_channels)
+
+    def __call__(self, time: Union[np.ndarray, float],
+                 data: Mapping[ChannelID, Union[np.ndarray, float]]) -> Mapping[ChannelID, Union[np.ndarray, float]]:
+        data_out = {forwarded_channel: data[forwarded_channel]
+                    for forwarded_channel in set(data.keys()).difference(self._input_channels)}
+
+        if len(data_out) == len(data):
+            # only forwarded data
+            return data_out
+
+        try:
+            data_in = np.stack([data[in_channel] for in_channel in self._input_channels])
+        except KeyError as error:
+            raise KeyError('Invalid input channels', set(data.keys()), set(self._input_channels)) from error
+
+        transformed_data = self._matrix @ data_in
+
+        for idx, out_channel in enumerate(self._output_channels):
+            data_out[out_channel] = transformed_data[idx, ...]
+
+        return data_out
+
+    def get_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        if not input_channels >= self._input_channels_set:
+            # input_channels is not a superset of the required input channels
+            raise KeyError('Invalid input channels', input_channels, self._input_channels_set)
+
+        return (input_channels - self._input_channels_set) | self._output_channels_set
+
+    def get_input_channels(self, output_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        forwarded = output_channels - self._output_channels_set
+        if not forwarded.isdisjoint(self._input_channels):
+            raise KeyError('Is input channel', forwarded & self._input_channels_set)
+        elif output_channels.isdisjoint(self._output_channels):
+            return output_channels
+        else:
+            return forwarded | self._input_channels_set
+
+    @property
+    def compare_key(self) -> Tuple[Tuple[ChannelID], Tuple[ChannelID], bytes]:
+        return self._input_channels, self._output_channels, self._matrix.tobytes()
+
+    def __repr__(self):
+        return ('LinearTransformation('
+                'transformation_matrix={transformation_matrix},'
+                'input_channels={input_channels},'
+                'output_channels={output_channels})').format(transformation_matrix=self._matrix.tolist(),
+                                                             input_channels=self._input_channels,
+                                                             output_channels=self._output_channels)
+
+    def is_constant_invariant(self):
+        """Signals if the transformation always maps constants to constants."""
+        return True
+
+    def get_constant_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        return input_channels
+
+
+class OffsetTransformation(Transformation):
+    def __init__(self, offsets: Mapping[ChannelID, _TrafoValue]):
+        """Adds an offset to each channel specified in offsets.
+
+        Channels not in offsets are forewarded
+
+        Args:
+            offsets: Channel -> offset mapping
+        """
+        self._offsets = frozendict(offsets)
+        assert _are_valid_transformation_expressions(self._offsets), f"Not valid transformation expressions: {self._offsets}"
+
+    def __call__(self, time: Union[np.ndarray, float],
+                 data: Mapping[ChannelID, Union[np.ndarray, float]]) -> Mapping[ChannelID, Union[np.ndarray, float]]:
+        offsets = _instantiate_expression_dict(time, self._offsets)
+        return {channel: channel_values + offsets[channel] if channel in offsets else channel_values
+                for channel, channel_values in data.items()}
+
+    def get_input_channels(self, output_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        return output_channels
+
+    def get_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        return input_channels
+
+    @property
+    def compare_key(self) -> Hashable:
+        return self._offsets
+
+    def __repr__(self):
+        return f'{type(self).__name__}({dict(self._offsets)!r})'
+
+    def is_constant_invariant(self):
+        """Signals if the transformation always maps constants to constants."""
+        return not _has_time_dependent_values(self._offsets)
+
+    def get_constant_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        return _get_constant_output_channels(self._offsets, input_channels)
+
+
+class ScalingTransformation(Transformation):
+    def __init__(self, factors: Mapping[ChannelID, _TrafoValue]):
+        self._factors = frozendict(factors)
+        assert _are_valid_transformation_expressions(self._factors), f"Not valid transformation expressions: {self._factors}"
+
+    def __call__(self, time: Union[np.ndarray, float],
+                 data: Mapping[ChannelID, Union[np.ndarray, float]]) -> Mapping[ChannelID, Union[np.ndarray, float]]:
+        factors = _instantiate_expression_dict(time, self._factors)
+        return {channel: channel_values * factors[channel] if channel in factors else channel_values
+                for channel, channel_values in data.items()}
+
+    def get_input_channels(self, output_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        return output_channels
+
+    def get_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        return input_channels
+
+    @property
+    def compare_key(self) -> Hashable:
+        return self._factors
+
+    def __repr__(self):
+        return f'{type(self).__name__}({dict(self._factors)!r})'
+
+    def is_constant_invariant(self):
+        """Signals if the transformation always maps constants to constants."""
+        return not _has_time_dependent_values(self._factors)
+
+    def get_constant_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        return _get_constant_output_channels(self._factors, input_channels)
+
+
+try:
+    if TYPE_CHECKING:
+        import pandas
+        PandasDataFrameType = pandas.DataFrame
+    else:
+        PandasDataFrameType = Any
+
+    def linear_transformation_from_pandas(transformation: PandasDataFrameType) -> LinearTransformation:
+        """ Creates a LinearTransformation object out of a pandas data frame.
+
+        Args:
+            transformation (pandas.DataFrame): The pandas.DataFrame object out of which a LinearTransformation will be formed.
+
+        Returns:
+            the created LinearTransformation instance
+        """
+        return LinearTransformation(transformation.values, transformation.columns, transformation.index)
+
+    LinearTransformation.from_pandas = linear_transformation_from_pandas
+except ImportError:
+    pass
+
+
+class ParallelChannelTransformation(Transformation):
+    def __init__(self, channels: Mapping[ChannelID, _TrafoValue]):
+        """Set channel values to given values regardless their former existence. The values can be time dependent
+        expressions.
+
+        Args:
+            channels: Channels present in this map are set to the given value.
+        """
+        self._channels: Mapping[ChannelID, _TrafoValue] = frozendict(channels.items())
+        assert _are_valid_transformation_expressions(self._channels), f"Not valid transformation expressions: {self._channels}"
+
+    def __call__(self, time: Union[np.ndarray, float],
+                 data: Mapping[ChannelID, Union[np.ndarray, float]]) -> Mapping[ChannelID, Union[np.ndarray, float]]:
+        overwritten = self._instantiated_values(time)
+        return {**data, **overwritten}
+
+    def _instantiated_values(self, time):
+        scope = {'t': time}
+        return {channel: value.evaluate_in_scope(scope) if hasattr(value, 'evaluate_in_scope') else np.full_like(time, fill_value=value, dtype=float)
+                for channel, value in self._channels.items()}
+
+    @property
+    def compare_key(self) -> Hashable:
+        return self._channels
+
+    def get_input_channels(self, output_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        return output_channels - self._channels.keys()
+
+    def get_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        return input_channels | self._channels.keys()
+
+    def __repr__(self):
+        return f'{type(self).__name__}({dict(self._channels)!r})'
+
+    def is_constant_invariant(self):
+        """Signals if the transformation always maps constants to constants."""
+        return not _has_time_dependent_values(self._channels)
+
+    def get_constant_output_channels(self, input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+        output_channels = set(input_channels)
+        for ch, value in self._channels.items():
+            if hasattr(value, 'variables'):
+                output_channels.discard(ch)
+            else:
+                output_channels.add(ch)
+
+        return output_channels
+
+
+def chain_transformations(*transformations: Transformation) -> Transformation:
+    parsed_transformations = []
+    for transformation in transformations:
+        if transformation is IdentityTransformation() or transformation is None:
+            pass
+        elif isinstance(transformation, ChainedTransformation):
+            parsed_transformations.extend(transformation.transformations)
+        else:
+            parsed_transformations.append(transformation)
+    if len(parsed_transformations) == 0:
+        return IdentityTransformation()
+    elif len(parsed_transformations) == 1:
+        return parsed_transformations[0]
+    else:
+        return ChainedTransformation(*parsed_transformations)
+
+
+def _instantiate_expression_dict(time, expressions: Mapping[str, _TrafoValue]) -> Mapping[str, Union[Real, np.ndarray]]:
+    scope = {'t': time}
+    modified_expressions = {}
+    for name, value in expressions.items():
+        if hasattr(value, 'evaluate_in_scope'):
+            modified_expressions[name] = value.evaluate_in_scope(scope)
+    if modified_expressions:
+        return {**expressions, **modified_expressions}
+    else:
+        return expressions
+
+
+def _has_time_dependent_values(expressions: Mapping[ChannelID, _TrafoValue]) -> bool:
+    return any(hasattr(value, 'variables')
+               for value in expressions.values())
+
+
+def _get_constant_output_channels(expressions: Mapping[ChannelID, _TrafoValue],
+                                  constant_input_channels: AbstractSet[ChannelID]) -> AbstractSet[ChannelID]:
+    return {ch
+            for ch in constant_input_channels
+            if not hasattr(expressions.get(ch, None), 'variables')}
+
+
+def _are_valid_transformation_expressions(expressions: Mapping[ChannelID, _TrafoValue]) -> bool:
+    return all(expr.variables == ('t',)
+               for expr in expressions.values()
+               if hasattr(expr, 'variables'))
diff --git a/qupulse/program/volatile.py b/qupulse/program/volatile.py
new file mode 100644
index 00000000..afd4692f
--- /dev/null
+++ b/qupulse/program/volatile.py
@@ -0,0 +1,73 @@
+from typing import NamedTuple, Mapping
+import warnings
+import numbers
+
+
+from qupulse.parameter_scope import Scope, MappedScope, JointScope
+from qupulse.expressions import Expression, ExpressionScalar
+from qupulse.utils.types import FrozenDict, FrozenMapping
+from qupulse.utils import is_integer
+
+
+__all__ = ['VolatileProperty', 'VolatileValue', 'VolatileRepetitionCount']
+
+
+VolatileProperty = NamedTuple('VolatileProperty', [('expression', Expression),
+                                                   ('dependencies', FrozenMapping[str, Expression])])
+VolatileProperty.__doc__ = """Hashable representation of a volatile program property. It does not contain the concrete
+value. Using the dependencies attribute to calculate the value might yield unexpected results."""
+
+
+class VolatileValue:
+    """Not hashable"""
+
+    def __init__(self, expression: ExpressionScalar, scope: Scope):
+        self._expression = expression
+        self._scope = scope
+
+    @property
+    def volatile_property(self) -> VolatileProperty:
+        dependencies = self._scope.get_volatile_parameters()
+        dependencies = FrozenDict({parameter_name: dependencies[parameter_name]
+                                   for parameter_name in self._expression.variables
+                                   if parameter_name in dependencies})
+        return VolatileProperty(expression=self._expression, dependencies=dependencies)
+
+    @classmethod
+    def operation(cls, expression, **operands):
+        expression = Expression(expression)
+        assert set(expression.variables) == operands.keys()
+
+        scope = JointScope(FrozenDict(
+            {operand_name: MappedScope(operand._scope, FrozenDict({operand_name: operand._expression}))
+             for operand_name, operand in operands.items()}
+        ))
+        return cls(expression, scope)
+    
+    def __sub__(self, other: int):
+        return type(self)(self._expression - other, self._scope)
+
+    def __mul__(self, other: int):
+        return type(self)(self._expression * other, self._scope)
+
+
+class VolatileRepetitionCount(VolatileValue):
+    def __int__(self):
+        value = self._expression.evaluate_in_scope(self._scope)
+        if not is_integer(value):
+            warnings.warn("Repetition count is no integer. Rounding might lead to unexpected results.")
+        value = int(round(value))
+        if value < 0:
+            warnings.warn("Repetition count is negative. Clamping lead to unexpected results.")
+            value = 0
+        return value
+
+    def update_volatile_dependencies(self, new_constants: Mapping[str, numbers.Number]) -> int:
+        self._scope = self._scope.change_constants(new_constants)
+        return int(self)
+
+    def __eq__(self, other):
+        if type(self) == type(other):
+            return self._scope is other._scope and self._expression == other._expression
+        else:
+            return NotImplemented
diff --git a/qupulse/program/waveforms.py b/qupulse/program/waveforms.py
new file mode 100644
index 00000000..9cdcbf4f
--- /dev/null
+++ b/qupulse/program/waveforms.py
@@ -0,0 +1,1234 @@
+"""This module contains all waveform classes
+
+Classes:
+    - Waveform: An instantiated pulse which can be sampled to a raw voltage value array.
+"""
+
+import itertools
+import operator
+import warnings
+from abc import ABCMeta, abstractmethod
+from numbers import Real
+from typing import (
+    AbstractSet, Any, FrozenSet, Iterable, Mapping, NamedTuple, Sequence, Set,
+    Tuple, Union, cast, Optional, List, Hashable)
+from weakref import WeakValueDictionary, ref
+
+import numpy as np
+
+from qupulse import ChannelID
+from qupulse.program.transformation import Transformation
+from qupulse.utils import checked_int_cast, isclose
+from qupulse.utils.types import TimeType, time_from_float
+from qupulse.utils.performance import is_monotonic
+from qupulse.comparable import Comparable
+from qupulse.expressions import ExpressionScalar
+from qupulse.pulses.interpolation import InterpolationStrategy
+from qupulse.utils import checked_int_cast, isclose
+from qupulse.utils.types import TimeType, time_from_float, FrozenDict
+from qupulse.program.transformation import Transformation
+from qupulse.utils import pairwise
+
+class ConstantFunctionPulseTemplateWarning(UserWarning):
+    """  This warning indicates a constant waveform is constructed from a FunctionPulseTemplate """
+    pass
+
+
+__all__ = ["Waveform", "TableWaveform", "TableWaveformEntry", "FunctionWaveform", "SequenceWaveform",
+           "MultiChannelWaveform", "RepetitionWaveform", "TransformingWaveform", "ArithmeticWaveform",
+           "ConstantFunctionPulseTemplateWarning", "ConstantWaveform"]
+
+PULSE_TO_WAVEFORM_ERROR = None  # error margin in pulse template to waveform conversion
+
+#  these are private because there probably will be changes here
+_ALLOCATION_FUNCTION = np.full_like  # pre_allocated = ALLOCATION_FUNCTION(sample_times, **ALLOCATION_FUNCTION_KWARGS)
+_ALLOCATION_FUNCTION_KWARGS = dict(fill_value=np.nan, dtype=float)
+
+
+def _to_time_type(duration: Real) -> TimeType:
+    if isinstance(duration, TimeType):
+        return duration
+    else:
+        return time_from_float(float(duration), absolute_error=PULSE_TO_WAVEFORM_ERROR)
+
+
+class Waveform(Comparable, metaclass=ABCMeta):
+    """Represents an instantiated PulseTemplate which can be sampled to retrieve arrays of voltage
+    values for the hardware."""
+
+    __sampled_cache = WeakValueDictionary()
+
+    __slots__ = ('_duration',)
+
+    def __init__(self, duration: TimeType):
+        self._duration = duration
+
+    @property
+    def duration(self) -> TimeType:
+        """The duration of the waveform in time units."""
+        return self._duration
+
+    @abstractmethod
+    def unsafe_sample(self,
+                      channel: ChannelID,
+                      sample_times: np.ndarray,
+                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
+        """Sample the waveform at given sample times.
+
+        The unsafe means that there are no sanity checks performed. The provided sample times are assumed to be
+        monotonously increasing and lie in the range of [0, waveform.duration]
+
+        Args:
+            sample_times: Times at which this Waveform will be sampled.
+            output_array: Has to be either None or an array of the same size and type as sample_times. If
+                not None, the sampled values will be written here and this array will be returned
+        Result:
+            The sampled values of this Waveform at the provided sample times. Has the same number of
+            elements as sample_times.
+        """
+
+    def get_sampled(self,
+                    channel: ChannelID,
+                    sample_times: np.ndarray,
+                    output_array: Union[np.ndarray, None] = None) -> np.ndarray:
+        """A wrapper to the unsafe_sample method which caches the result. This method enforces the constrains
+        unsafe_sample expects and caches the result to save memory.
+
+        Args:
+            sample_times: Times at which this Waveform will be sampled.
+            output_array: Has to be either None or an array of the same size and type as sample_times. If an array is
+                given, the sampled values will be written into the given array and it will be returned. Otherwise, a new
+                array will be created and cached to save memory.
+
+        Result:
+            The sampled values of this Waveform at the provided sample times. Is `output_array` if provided
+        """
+        if len(sample_times) == 0:
+            if output_array is None:
+                return np.zeros_like(sample_times, dtype=float)
+            elif len(output_array) == len(sample_times):
+                return output_array
+            else:
+                raise ValueError('Output array length and sample time length are different')
+
+        if not is_monotonic(sample_times):
+            raise ValueError('The sample times are not monotonously increasing')
+        if sample_times[0] < 0 or sample_times[-1] > float(self.duration):
+            raise ValueError(f'The sample times [{sample_times[0]}, ..., {sample_times[-1]}] are not in the range'
+                             f' [0, duration={float(self.duration)}]')
+        if channel not in self.defined_channels:
+            raise KeyError('Channel not defined in this waveform: {}'.format(channel))
+
+        constant_value = self.constant_value(channel)
+        if constant_value is None:
+            if output_array is None:
+                # cache the result to save memory
+                result = self.unsafe_sample(channel, sample_times)
+                result.flags.writeable = False
+                key = hash(bytes(result))
+                if key not in self.__sampled_cache:
+                    self.__sampled_cache[key] = result
+                return self.__sampled_cache[key]
+            else:
+                if len(output_array) != len(sample_times):
+                    raise ValueError('Output array length and sample time length are different')
+                # use the user provided memory
+                return self.unsafe_sample(channel=channel,
+                                          sample_times=sample_times,
+                                          output_array=output_array)
+        else:
+            if output_array is None:
+                output_array = np.full_like(sample_times, fill_value=constant_value, dtype=float)
+            else:
+                output_array[:] = constant_value
+            return output_array
+
+    @property
+    @abstractmethod
+    def defined_channels(self) -> AbstractSet[ChannelID]:
+        """The channels this waveform should played on. Use
+            :func:`~qupulse.pulses.instructions.get_measurement_windows` to get a waveform for a subset of these."""
+
+    @abstractmethod
+    def unsafe_get_subset_for_channels(self, channels: AbstractSet[ChannelID]) -> 'Waveform':
+        """Unsafe version of :func:`~qupulse.pulses.instructions.get_measurement_windows`."""
+
+    def get_subset_for_channels(self, channels: AbstractSet[ChannelID]) -> 'Waveform':
+        """Get a waveform that only describes the channels contained in `channels`.
+
+        Args:
+            channels: A channel set the return value should confine to.
+
+        Raises:
+            KeyError: If `channels` is not a subset of the waveform's defined channels.
+
+        Returns:
+            A waveform with waveform.defined_channels == `channels`
+        """
+        if not channels <= self.defined_channels:
+            raise KeyError('Channels not defined on waveform: {}'.format(channels))
+        if channels == self.defined_channels:
+            return self
+        return self.unsafe_get_subset_for_channels(channels=channels)
+
+    def is_constant(self) -> bool:
+        """Convenience function to check if all channels are constant. The result is equal to
+        `all(waveform.constant_value(ch) is not None for ch in waveform.defined_channels)` but might be more performant.
+
+        Returns:
+            True if all channels have constant values.
+        """
+        return self.constant_value_dict() is not None
+
+    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
+        result = {ch: self.constant_value(ch) for ch in self.defined_channels}
+        if None in result.values():
+            return None
+        else:
+            return result
+
+    def constant_value(self, channel: ChannelID) -> Optional[float]:
+        """Checks if the requested channel has a constant value and returns it if so.
+
+        Guarantee that this assertion passes for every t in waveform duration:
+        >>> assert waveform.constant_value(channel) is None or waveform.constant_value(t) = waveform.get_sampled(channel, t)
+
+        Args:
+            channel: The channel to check
+
+        Returns:
+            None if there is no guarantee that the channel is constant. The value otherwise.
+        """
+        return None
+
+    def __neg__(self):
+        return FunctorWaveform.from_functor(self, {ch: np.negative for ch in self.defined_channels})
+
+    def __pos__(self):
+        return self
+
+    def _sort_key_for_channels(self) -> Sequence[Tuple[str, int]]:
+        """Makes reproducible sorting by defined channels possible"""
+        return sorted((ch, 0) if isinstance(ch, str) else ('', ch) for ch in self.defined_channels)
+
+    def reversed(self) -> 'Waveform':
+        """Returns a reversed version of this waveform."""
+        # We don't check for constness here because const waveforms are supposed to override this method
+        return ReversedWaveform(self)
+
+
+class TableWaveformEntry(NamedTuple('TableWaveformEntry', [('t', Real),
+                                                           ('v', float),
+                                                           ('interp', InterpolationStrategy)])):
+    def __init__(self, t: float, v: float, interp: InterpolationStrategy):
+        if not callable(interp):
+            raise TypeError('{} is neither callable nor of type InterpolationStrategy'.format(interp))
+
+    def __repr__(self):
+        return f'{type(self).__name__}(t={self.t!r}, v={self.v!r}, interp={self.interp!r})'
+
+
+class TableWaveform(Waveform):
+    EntryInInit = Union[TableWaveformEntry, Tuple[float, float, InterpolationStrategy]]
+
+    """Waveform obtained from instantiating a TablePulseTemplate."""
+
+    __slots__ = ('_table', '_channel_id')
+
+    def __init__(self,
+                 channel: ChannelID,
+                 waveform_table: Tuple[TableWaveformEntry, ...]) -> None:
+        """Create a new TableWaveform instance.
+
+        Args:
+            waveform_table: A tuple of instantiated and validated table entries
+        """
+        if not isinstance(waveform_table, tuple):
+            warnings.warn("Please use a tuple of TableWaveformEntry to construct TableWaveform directly",
+                          category=DeprecationWarning)
+            waveform_table = self._validate_input(waveform_table)
+
+        super().__init__(duration=_to_time_type(waveform_table[-1].t))
+
+        self._table = waveform_table
+        self._channel_id = channel
+
+    @staticmethod
+    def _validate_input(input_waveform_table: Sequence[EntryInInit]) -> Union[Tuple[Real, Real],
+                                                                              List[TableWaveformEntry]]:
+        """ Checks that:
+         - the time is increasing,
+         - there are at least two entries
+
+        Optimizations:
+          - removes subsequent entries with same time or voltage values.
+          - checks if the complete waveform is constant. Returns a (duration, value) tuple if this is the case
+
+        Raises:
+            ValueError:
+              - there are less than two entries
+              - the entries are not ordered in time
+              - Any time is negative
+              - The total length is zero
+
+        Returns:
+            A list of de-duplicated table entries
+            OR
+            A (duration, value) tuple if the waveform is constant
+        """
+        # we use an iterator here to avoid duplicate work and be maximally efficient for short tables
+        # We never use StopIteration to abort iteration. It always signifies an error.
+        input_iter = iter(input_waveform_table)
+        try:
+            first_t, first_v, first_interp = next(input_iter)
+        except StopIteration:
+            raise ValueError("Waveform table mut not be empty")
+
+        if first_t != 0.0:
+            raise ValueError('First time entry is not zero.')
+
+        previous_t = 0.0
+        previous_v = first_v
+        output_waveform_table = [TableWaveformEntry(0.0, first_v, first_interp)]
+
+        try:
+            t, v, interp = next(input_iter)
+        except StopIteration:
+            raise ValueError("Waveform table has less than two entries.")
+        if t < 0:
+            raise ValueError('Negative time values are not allowed.')
+
+        # constant_v is None <=> the waveform is constant until up to the current entry
+        constant_v = interp.constant_value((previous_t, previous_v), (t, v))
+
+        for next_t, next_v, next_interp in input_iter:
+            if next_t < t:
+                if next_t < 0:
+                    raise ValueError('Negative time values are not allowed.')
+                else:
+                    raise ValueError('Times are not increasing.')
+
+            if constant_v is not None and interp.constant_value((t, v), (next_t, next_v)) != constant_v:
+                constant_v = None
+
+            if (previous_t != t or t != next_t) and (previous_v != v or v != next_v):
+                # the time and the value differ both either from the next or the previous
+                # otherwise we skip the entry
+                previous_t = t
+                previous_v = v
+                output_waveform_table.append(TableWaveformEntry(t, v, interp))
+
+            t, v, interp = next_t, next_v, next_interp
+
+        # Until now, we only checked that the time does not decrease. We require an increase because duration == 0
+        # waveforms are ill-formed. t is now the time of the last entry.
+        if t == 0:
+            raise ValueError('Last time entry is zero.')
+
+        if constant_v is not None:
+            # the waveform is constant
+            return t, constant_v
+        else:
+            # we must still add the last entry to the table
+            output_waveform_table.append(TableWaveformEntry(t, v, interp))
+            return output_waveform_table
+
+    def is_constant(self) -> bool:
+        # only correct if `from_table` is used
+        return False
+
+    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
+        # only correct if `from_table` is used
+        return None
+
+    @classmethod
+    def from_table(cls, channel: ChannelID, table: Sequence[EntryInInit]) -> Union['TableWaveform', 'ConstantWaveform']:
+        table = cls._validate_input(table)
+        if isinstance(table, tuple):
+            duration, amplitude = table
+            return ConstantWaveform(duration=duration, amplitude=amplitude, channel=channel)
+        else:
+            return TableWaveform(channel, tuple(table))
+
+    @property
+    def compare_key(self) -> Any:
+        return self._channel_id, self._table
+
+    def unsafe_sample(self,
+                      channel: ChannelID,
+                      sample_times: np.ndarray,
+                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
+        if output_array is None:
+            output_array = _ALLOCATION_FUNCTION(sample_times, **_ALLOCATION_FUNCTION_KWARGS)
+
+        if PULSE_TO_WAVEFORM_ERROR:
+            # we need to replace the last entry's t with self.duration
+            *entries, last = self._table
+            entries.append(TableWaveformEntry(float(self.duration), last.v, last.interp))
+        else:
+            entries = self._table
+
+        for entry1, entry2 in pairwise(entries):
+            indices = slice(np.searchsorted(sample_times, entry1.t, 'left'),
+                            np.searchsorted(sample_times, entry2.t, 'right'))
+            output_array[indices] = \
+                entry2.interp((float(entry1.t), entry1.v),
+                              (float(entry2.t), entry2.v),
+                              sample_times[indices])
+        return output_array
+
+    @property
+    def defined_channels(self) -> AbstractSet[ChannelID]:
+        return {self._channel_id}
+
+    def unsafe_get_subset_for_channels(self, channels: AbstractSet[ChannelID]) -> 'Waveform':
+        return self
+
+    def __repr__(self):
+        return f'{type(self).__name__}(channel={self._channel_id!r}, waveform_table={self._table!r})'
+
+
+class ConstantWaveform(Waveform):
+
+    # TODO: remove
+    _is_constant_waveform = True
+
+    __slots__ = ('_amplitude', '_channel')
+
+    def __init__(self, duration: Real, amplitude: Any, channel: ChannelID):
+        """ Create a qupulse waveform corresponding to a ConstantPulseTemplate """
+        super().__init__(duration=_to_time_type(duration))
+        self._amplitude = amplitude
+        self._channel = channel
+
+    @classmethod
+    def from_mapping(cls, duration: Real, constant_values: Mapping[ChannelID, float]) -> Union['ConstantWaveform',
+                                                                                               'MultiChannelWaveform']:
+        """Construct a ConstantWaveform or a MultiChannelWaveform of ConstantWaveforms with given duration and values"""
+        assert constant_values
+        duration = _to_time_type(duration)
+        if len(constant_values) == 1:
+            (channel, amplitude), = constant_values.items()
+            return cls(duration, amplitude=amplitude, channel=channel)
+        else:
+            return MultiChannelWaveform([cls(duration, amplitude=amplitude, channel=channel)
+                                         for channel, amplitude in constant_values.items()])
+
+    def is_constant(self) -> bool:
+        return True
+
+    def constant_value(self, channel: ChannelID) -> Optional[float]:
+        assert channel == self._channel
+        return self._amplitude
+
+    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
+        return {self._channel: self._amplitude}
+
+    @property
+    def defined_channels(self) -> AbstractSet[ChannelID]:
+        """The channels this waveform should played on. Use
+            :func:`~qupulse.pulses.instructions.get_measurement_windows` to get a waveform for a subset of these."""
+
+        return {self._channel}
+
+    @property
+    def compare_key(self) -> Tuple[Any, ...]:
+        return self._duration, self._amplitude, self._channel
+
+    def unsafe_sample(self,
+                      channel: ChannelID,
+                      sample_times: np.ndarray,
+                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
+        if output_array is None:
+            return np.full_like(sample_times, fill_value=self._amplitude, dtype=float)
+        else:
+            output_array[:] = self._amplitude
+            return output_array
+
+    def unsafe_get_subset_for_channels(self, channels: Set[ChannelID]) -> Waveform:
+        """Unsafe version of :func:`~qupulse.pulses.instructions.get_measurement_windows`."""
+        return self
+
+    def __repr__(self):
+        return f"{type(self).__name__}(duration={self.duration!r}, "\
+               f"amplitude={self._amplitude!r}, channel={self._channel!r})"
+
+    def reversed(self) -> 'Waveform':
+        return self
+
+
+class FunctionWaveform(Waveform):
+    """Waveform obtained from instantiating a FunctionPulseTemplate."""
+
+    __slots__ = ('_expression', '_channel_id')
+
+    def __init__(self, expression: ExpressionScalar,
+                 duration: float,
+                 channel: ChannelID) -> None:
+        """Creates a new FunctionWaveform instance.
+
+        Args:
+            expression: The function represented by this FunctionWaveform
+                as a mathematical expression where 't' denotes the time variable. It must not have other variables
+            duration: The duration of the waveform
+            measurement_windows: A list of measurement windows
+            channel: The channel this waveform is played on
+        """
+
+        if set(expression.variables) - set('t'):
+            raise ValueError('FunctionWaveforms may not depend on anything but "t"')
+        elif not expression.variables:
+            warnings.warn("Constant FunctionWaveform is not recommended as the constant propagation will be suboptimal",
+                          category=ConstantFunctionPulseTemplateWarning)
+        super().__init__(duration=_to_time_type(duration))
+        self._expression = expression
+        self._channel_id = channel
+
+    @classmethod
+    def from_expression(cls, expression: ExpressionScalar, duration: float, channel: ChannelID) -> Union['FunctionWaveform', ConstantWaveform]:
+        if expression.variables:
+            return cls(expression, duration, channel)
+        else:
+            return ConstantWaveform(amplitude=expression.evaluate_numeric(), duration=duration, channel=channel)
+
+    def is_constant(self) -> bool:
+        # only correct if `from_expression` is used
+        return False
+
+    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
+        # only correct if `from_expression` is used
+        return None
+
+    @property
+    def defined_channels(self) -> AbstractSet[ChannelID]:
+        return {self._channel_id}
+
+    @property
+    def compare_key(self) -> Any:
+        return self._channel_id, self._expression, self._duration
+
+    @property
+    def duration(self) -> TimeType:
+        return self._duration
+
+    def unsafe_sample(self,
+                      channel: ChannelID,
+                      sample_times: np.ndarray,
+                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
+        evaluated = self._expression.evaluate_numeric(t=sample_times)
+        if output_array is None:
+            if self._expression.variables:
+                return evaluated.astype(float)
+            else:
+                return np.full_like(sample_times, fill_value=float(evaluated), dtype=float)
+        else:
+            output_array[:] = evaluated
+            return output_array
+
+    def unsafe_get_subset_for_channels(self, channels: AbstractSet[ChannelID]) -> Waveform:
+        return self
+
+    def __repr__(self):
+        return f"{type(self).__name__}(duration={self.duration!r}, "\
+               f"expression={self._expression!r}, channel={self._channel_id!r})"
+
+
+class SequenceWaveform(Waveform):
+    """This class allows putting multiple PulseTemplate together in one waveform on the hardware."""
+
+    __slots__ = ('_sequenced_waveforms', )
+
+    def __init__(self, sub_waveforms: Iterable[Waveform]):
+        """Use Waveform.from_sequence for optimal construction
+
+        :param subwaveforms: All waveforms must have the same defined channels
+        """
+        if not sub_waveforms:
+            raise ValueError(
+                "SequenceWaveform cannot be constructed without channel waveforms."
+            )
+
+        # do not fail on iterators although we do not allow them as an argument
+        sequenced_waveforms = tuple(sub_waveforms)
+
+        super().__init__(duration=sum(waveform.duration for waveform in sequenced_waveforms))
+        self._sequenced_waveforms = sequenced_waveforms
+
+        defined_channels = self._sequenced_waveforms[0].defined_channels
+        if not all(waveform.defined_channels == defined_channels
+                   for waveform in itertools.islice(self._sequenced_waveforms, 1, None)):
+            for waveform in self._sequenced_waveforms[1:]:
+                 if not waveform.defined_channels == self.defined_channels:
+                     print(f"SequenceWaveform: defined channels {self.defined_channels} do not match {waveform.defined_channels} ")
+            raise ValueError(
+                "SequenceWaveform cannot be constructed from waveforms of different"
+                "defined channels."
+            )
+
+    @classmethod
+    def from_sequence(cls, waveforms: Sequence['Waveform']) -> 'Waveform':
+        """Returns a waveform the represents the given sequence of waveforms. Applies some optimizations."""
+        assert waveforms, "Sequence must not be empty"
+        if len(waveforms) == 1:
+            return waveforms[0]
+
+        flattened = []
+        constant_values = waveforms[0].constant_value_dict()
+        for wf in waveforms:
+            if constant_values and constant_values != wf.constant_value_dict():
+                constant_values = None
+            if isinstance(wf, cls):
+                flattened.extend(wf.sequenced_waveforms)
+            else:
+                flattened.append(wf)
+        if constant_values is None:
+            return cls(sub_waveforms=flattened)
+        else:
+            duration = sum(wf.duration for wf in flattened)
+            return ConstantWaveform.from_mapping(duration, constant_values)
+
+    def is_constant(self) -> bool:
+        # only correct if from_sequence is used for construction
+        return False
+
+    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
+        # only correct if from_sequence is used for construction
+        return None
+
+    def constant_value(self, channel: ChannelID) -> Optional[float]:
+        v = None
+        for wf in self._sequenced_waveforms:
+            wf_cv = wf.constant_value(channel)
+            if wf_cv is None:
+                return None
+            elif wf_cv == v:
+                continue
+            elif v is None:
+                v = wf_cv
+            else:
+                assert v != wf_cv
+                return None
+        return v
+
+    @property
+    def defined_channels(self) -> AbstractSet[ChannelID]:
+        return self._sequenced_waveforms[0].defined_channels
+
+    def unsafe_sample(self,
+                      channel: ChannelID,
+                      sample_times: np.ndarray,
+                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
+        if output_array is None:
+            output_array = _ALLOCATION_FUNCTION(sample_times, **_ALLOCATION_FUNCTION_KWARGS)
+        time = 0
+        for subwaveform in self._sequenced_waveforms:
+            # before you change anything here, make sure to understand the difference between basic and advanced
+            # indexing in numpy and their copy/reference behaviour
+            end = time + subwaveform.duration
+
+            indices = slice(*np.searchsorted(sample_times, (float(time), float(end)), 'left'))
+            subwaveform.unsafe_sample(channel=channel,
+                                      sample_times=sample_times[indices]-np.float64(time),
+                                      output_array=output_array[indices])
+            time = end
+        return output_array
+
+    @property
+    def compare_key(self) -> Tuple[Waveform]:
+        return self._sequenced_waveforms
+
+    @property
+    def duration(self) -> TimeType:
+        return self._duration
+
+    def unsafe_get_subset_for_channels(self, channels: AbstractSet[ChannelID]) -> 'Waveform':
+        return SequenceWaveform.from_sequence([
+            sub_waveform.unsafe_get_subset_for_channels(channels & sub_waveform.defined_channels)
+            for sub_waveform in self._sequenced_waveforms if sub_waveform.defined_channels & channels])
+
+    @property
+    def sequenced_waveforms(self) -> Sequence[Waveform]:
+        return self._sequenced_waveforms
+
+    def __repr__(self):
+        return f"{type(self).__name__}({self._sequenced_waveforms})"
+
+
+class MultiChannelWaveform(Waveform):
+    """A MultiChannelWaveform is a Waveform object that allows combining arbitrary Waveform objects
+    to into a single waveform defined for several channels.
+
+    The number of channels used by the MultiChannelWaveform object is the sum of the channels used
+    by the Waveform objects it consists of.
+
+    MultiChannelWaveform allows an arbitrary mapping of channels defined by the Waveforms it
+    consists of and the channels it defines. For example, if the MultiChannelWaveform consists
+    of a two Waveform objects A and B which define two channels each, then the channels of the
+    MultiChannelWaveform may be 0: A.1, 1: B.0, 2: B.1, 3: A.0 where A.0 means channel 0 of Waveform
+    object A.
+
+    The following constraints must hold:
+     - The durations of all Waveform objects must be equal.
+     - The channel mapping must be sane, i.e., no channel of the MultiChannelWaveform must be
+        assigned more than one channel of any Waveform object it consists of
+    """
+
+    __slots__ = ('_sub_waveforms', '_defined_channels')
+
+    def __init__(self, sub_waveforms: List[Waveform]) -> None:
+        """Create a new MultiChannelWaveform instance.
+        Use `MultiChannelWaveform.from_parallel` for optimal construction.
+
+        Requires a list of subwaveforms in the form (Waveform, List(int)) where the list defines
+        the channel mapping, i.e., a value y at index x in the list means that channel x of the
+        subwaveform will be mapped to channel y of this MultiChannelWaveform object.
+
+        Args:
+            sub_waveforms: The list of sub waveforms of this
+                MultiChannelWaveform. List might get sorted!
+        Raises:
+            ValueError, if a channel mapping is out of bounds of the channels defined by this
+                MultiChannelWaveform
+            ValueError, if several subwaveform channels are assigned to a single channel of this
+                MultiChannelWaveform
+            ValueError, if subwaveforms have inconsistent durations
+        """
+
+        if not sub_waveforms:
+            raise ValueError(
+                "MultiChannelWaveform cannot be constructed without channel waveforms."
+            )
+
+        # sort the waveforms with their defined channels to make compare key reproducible
+        if not isinstance(sub_waveforms, list):
+            sub_waveforms = list(sub_waveforms)
+        sub_waveforms.sort(key=lambda wf: wf._sort_key_for_channels())
+
+        super().__init__(duration=sub_waveforms[0].duration)
+        self._sub_waveforms = tuple(sub_waveforms)
+
+        defined_channels = set()
+        for waveform in self._sub_waveforms:
+            if waveform.defined_channels & defined_channels:
+                raise ValueError('Channel may not be defined in multiple waveforms',
+                                 waveform.defined_channels & defined_channels)
+            defined_channels |= waveform.defined_channels
+        self._defined_channels = frozenset(defined_channels)
+
+        if not all(isclose(waveform.duration, self.duration) for waveform in self._sub_waveforms[1:]):
+            # meaningful error message:
+            durations = {}
+
+            for waveform in self._sub_waveforms:
+                for duration, channels in durations.items():
+                    if isclose(waveform.duration, duration):
+                        channels.update(waveform.defined_channels)
+                        break
+                else:
+                    durations[waveform.duration] = set(waveform.defined_channels)
+
+            raise ValueError(
+                "MultiChannelWaveform cannot be constructed from channel waveforms of different durations.",
+                durations
+            )
+
+    @staticmethod
+    def from_parallel(waveforms: Sequence[Waveform]) -> Waveform:
+        assert waveforms, "ARgument must not be empty"
+        if len(waveforms) == 1:
+            return waveforms[0]
+
+        # we do not look at constant values here because there is no benefit. We would need to construct a new
+        # MultiChannelWaveform anyways
+
+        # avoid unnecessary multi channel nesting
+        flattened = []
+        for waveform in waveforms:
+            if isinstance(waveform, MultiChannelWaveform):
+                flattened.extend(waveform._sub_waveforms)
+            else:
+                flattened.append(waveform)
+
+        return MultiChannelWaveform(flattened)
+
+    def is_constant(self) -> bool:
+        return all(wf.is_constant() for wf in self._sub_waveforms)
+
+    def constant_value(self, channel: ChannelID) -> Optional[float]:
+        return self[channel].constant_value(channel)
+
+    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
+        d = {}
+        for wf in self._sub_waveforms:
+            wf_d = wf.constant_value_dict()
+            if wf_d is None:
+                return None
+            else:
+                d.update(wf_d)
+        return d
+
+    @property
+    def duration(self) -> TimeType:
+        return self._sub_waveforms[0].duration
+
+    def __getitem__(self, key: ChannelID) -> Waveform:
+        for waveform in self._sub_waveforms:
+            if key in waveform.defined_channels:
+                return waveform
+        raise KeyError('Unknown channel ID: {}'.format(key), key)
+
+    @property
+    def defined_channels(self) -> AbstractSet[ChannelID]:
+        return self._defined_channels
+
+    @property
+    def compare_key(self) -> Any:
+        # sort with channels
+        return self._sub_waveforms
+
+    def unsafe_sample(self,
+                      channel: ChannelID,
+                      sample_times: np.ndarray,
+                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
+        return self[channel].unsafe_sample(channel, sample_times, output_array)
+
+    def unsafe_get_subset_for_channels(self, channels: AbstractSet[ChannelID]) -> 'Waveform':
+        relevant_sub_waveforms = [swf for swf in self._sub_waveforms if swf.defined_channels & channels]
+        if len(relevant_sub_waveforms) == 1:
+            return relevant_sub_waveforms[0].get_subset_for_channels(channels)
+        elif len(relevant_sub_waveforms) > 1:
+            return MultiChannelWaveform.from_parallel(
+                [sub_waveform.get_subset_for_channels(channels & sub_waveform.defined_channels)
+                 for sub_waveform in relevant_sub_waveforms])
+        else:
+            raise KeyError('Unknown channels: {}'.format(channels))
+
+    def __repr__(self):
+        return f"{type(self).__name__}({self._sub_waveforms!r})"
+
+
+class RepetitionWaveform(Waveform):
+    """This class allows putting multiple PulseTemplate together in one waveform on the hardware."""
+
+    __slots__ = ('_body', '_repetition_count')
+
+    def __init__(self, body: Waveform, repetition_count: int):
+        repetition_count = checked_int_cast(repetition_count)
+        if repetition_count < 1 or not isinstance(repetition_count, int):
+            raise ValueError('Repetition count must be an integer >0')
+
+        super().__init__(duration=body.duration * repetition_count)
+        self._body = body
+        self._repetition_count = repetition_count
+
+    @classmethod
+    def from_repetition_count(cls, body: Waveform, repetition_count: int) -> Waveform:
+        constant_values = body.constant_value_dict()
+        if constant_values is None:
+            return RepetitionWaveform(body, repetition_count)
+        else:
+            return ConstantWaveform.from_mapping(body.duration * repetition_count, constant_values)
+
+    @property
+    def defined_channels(self) -> AbstractSet[ChannelID]:
+        return self._body.defined_channels
+
+    def unsafe_sample(self,
+                      channel: ChannelID,
+                      sample_times: np.ndarray,
+                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
+        if output_array is None:
+            output_array = _ALLOCATION_FUNCTION(sample_times, **_ALLOCATION_FUNCTION_KWARGS)
+        body_duration = self._body.duration
+        time = 0
+        for _ in range(self._repetition_count):
+            end = time + body_duration
+            indices = slice(*np.searchsorted(sample_times, (float(time), float(end)), 'left'))
+            self._body.unsafe_sample(channel=channel,
+                                     sample_times=sample_times[indices] - float(time),
+                                     output_array=output_array[indices])
+            time = end
+        return output_array
+
+    @property
+    def compare_key(self) -> Tuple[Any, int]:
+        return self._body.compare_key, self._repetition_count
+
+    def unsafe_get_subset_for_channels(self, channels: AbstractSet[ChannelID]) -> Waveform:
+        return RepetitionWaveform.from_repetition_count(
+            body=self._body.unsafe_get_subset_for_channels(channels),
+            repetition_count=self._repetition_count)
+
+    def is_constant(self) -> bool:
+        return self._body.is_constant()
+
+    def constant_value(self, channel: ChannelID) -> Optional[float]:
+        return self._body.constant_value(channel)
+
+    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
+        return self._body.constant_value_dict()
+
+    def __repr__(self):
+        return f"{type(self).__name__}(body={self._body!r}, repetition_count={self._repetition_count!r})"
+
+
+class TransformingWaveform(Waveform):
+    __slots__ = ('_inner_waveform', '_transformation', '_cached_data', '_cached_times')
+
+    def __init__(self, inner_waveform: Waveform, transformation: Transformation):
+        """"""
+        super(TransformingWaveform, self).__init__(duration=inner_waveform.duration)
+        self._inner_waveform = inner_waveform
+        self._transformation = transformation
+
+        # cache data of inner channels based identified and invalidated by the sample times
+        self._cached_data = None
+        self._cached_times = lambda: None
+
+    @classmethod
+    def from_transformation(cls, inner_waveform: Waveform, transformation: Transformation) -> Waveform:
+        constant_values = inner_waveform.constant_value_dict()
+
+        if constant_values is None or not transformation.is_constant_invariant():
+            return cls(inner_waveform, transformation)
+
+        transformed_constant_values = {key: float(value) for key, value in transformation(0., constant_values).items()}
+        return ConstantWaveform.from_mapping(inner_waveform.duration, transformed_constant_values)
+
+    def is_constant(self) -> bool:
+        # only true if `from_transformation` was used
+        return False
+
+    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
+        # only true if `from_transformation` was used
+        return None
+
+    def constant_value(self, channel: ChannelID) -> Optional[float]:
+        if not self._transformation.is_constant_invariant():
+            return None
+        in_channels = self._transformation.get_input_channels({channel})
+        in_values = {ch: self._inner_waveform.constant_value(ch) for ch in in_channels}
+        if any(val is None for val in in_values.values()):
+            return None
+        else:
+            return self._transformation(0., in_values)[channel]
+
+    @property
+    def inner_waveform(self) -> Waveform:
+        return self._inner_waveform
+
+    @property
+    def transformation(self) -> Transformation:
+        return self._transformation
+
+    @property
+    def defined_channels(self) -> AbstractSet[ChannelID]:
+        return self.transformation.get_output_channels(self.inner_waveform.defined_channels)
+
+    @property
+    def compare_key(self) -> Tuple[Waveform, Transformation]:
+        return self.inner_waveform, self.transformation
+
+    def unsafe_get_subset_for_channels(self, channels: Set[ChannelID]) -> 'SubsetWaveform':
+        return SubsetWaveform(self, channel_subset=channels)
+
+    def unsafe_sample(self,
+                      channel: ChannelID,
+                      sample_times: np.ndarray,
+                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
+        if self._cached_times() is not sample_times:
+            self._cached_data = dict()
+            self._cached_times = ref(sample_times)
+
+        if channel not in self._cached_data:
+
+            inner_channels = self.transformation.get_input_channels({channel})
+
+            inner_data = {inner_channel: self.inner_waveform.unsafe_sample(inner_channel, sample_times)
+                          for inner_channel in inner_channels}
+
+            outer_data = self.transformation(sample_times, inner_data)
+            self._cached_data.update(outer_data)
+
+        if output_array is None:
+            output_array = self._cached_data[channel]
+        else:
+            output_array[:] = self._cached_data[channel]
+
+        return output_array
+
+
+class SubsetWaveform(Waveform):
+    __slots__ = ('_inner_waveform', '_channel_subset')
+
+    def __init__(self, inner_waveform: Waveform, channel_subset: Set[ChannelID]):
+        super().__init__(duration=inner_waveform.duration)
+        self._inner_waveform = inner_waveform
+        self._channel_subset = frozenset(channel_subset)
+
+    @property
+    def inner_waveform(self) -> Waveform:
+        return self._inner_waveform
+
+    @property
+    def defined_channels(self) -> FrozenSet[ChannelID]:
+        return self._channel_subset
+
+    @property
+    def compare_key(self) -> Tuple[frozenset, Waveform]:
+        return self.defined_channels, self.inner_waveform
+
+    def unsafe_get_subset_for_channels(self, channels: Set[ChannelID]) -> Waveform:
+        return self.inner_waveform.get_subset_for_channels(channels)
+
+    def unsafe_sample(self,
+                      channel: ChannelID,
+                      sample_times: np.ndarray,
+                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
+        return self.inner_waveform.unsafe_sample(channel, sample_times, output_array)
+
+    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
+        d = self._inner_waveform.constant_value_dict()
+        if d is not None:
+            return {ch: d[ch] for ch in self._channel_subset}
+
+    def constant_value(self, channel: ChannelID) -> Optional[float]:
+        if channel not in self._channel_subset:
+            raise KeyError(channel)
+        return self._inner_waveform.constant_value(channel)
+
+
+class ArithmeticWaveform(Waveform):
+    """Channels only present in one waveform have the operations neutral element on the other."""
+
+    numpy_operator_map = {'+': np.add,
+                          '-': np.subtract}
+    operator_map = {'+': operator.add,
+                    '-': operator.sub}
+
+    rhs_only_map = {'+': operator.pos,
+                    '-': operator.neg}
+    numpy_rhs_only_map = {'+': np.positive,
+                          '-': np.negative}
+
+    __slots__ = ('_lhs', '_rhs', '_arithmetic_operator')
+
+    def __init__(self,
+                 lhs: Waveform,
+                 arithmetic_operator: str,
+                 rhs: Waveform):
+        super().__init__(duration=lhs.duration)
+        self._lhs = lhs
+        self._rhs = rhs
+        self._arithmetic_operator = arithmetic_operator
+
+        assert np.isclose(float(self._lhs.duration), float(self._rhs.duration))
+        assert arithmetic_operator in self.operator_map
+
+    @classmethod
+    def from_operator(cls, lhs: Waveform, arithmetic_operator: str, rhs: Waveform):
+        # one could optimize rhs_cv to being only created if lhs_cv is not None but this makes the code harder to read
+        lhs_cv = lhs.constant_value_dict()
+        rhs_cv = rhs.constant_value_dict()
+        if lhs_cv is None or rhs_cv is None:
+            return cls(lhs, arithmetic_operator, rhs)
+
+        else:
+            constant_values = dict(lhs_cv)
+            op = cls.operator_map[arithmetic_operator]
+            rhs_op = cls.rhs_only_map[arithmetic_operator]
+
+            for ch, rhs_val in rhs_cv.items():
+                if ch in constant_values:
+                    constant_values[ch] = op(constant_values[ch], rhs_val)
+                else:
+                    constant_values[ch] = rhs_op(rhs_val)
+
+            duration = lhs.duration
+            assert isclose(duration, rhs.duration)
+
+            return ConstantWaveform.from_mapping(duration, constant_values)
+
+    def constant_value(self, channel: ChannelID) -> Optional[float]:
+        if channel not in self._rhs.defined_channels:
+            return self._lhs.constant_value(channel)
+        rhs = self._rhs.constant_value(channel)
+        if rhs is None:
+            return None
+
+        if channel in self._lhs.defined_channels:
+            lhs = self._lhs.constant_value(channel)
+            if lhs is None:
+                return None
+
+            return self.operator_map[self._arithmetic_operator](lhs, rhs)
+        else:
+            return self.rhs_only_map[self._arithmetic_operator](rhs)
+
+    def is_constant(self) -> bool:
+        # only correct if from_operator is used
+        return False
+
+    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
+        # only correct if from_operator is used
+        return None
+
+    @property
+    def lhs(self) -> Waveform:
+        return self._lhs
+
+    @property
+    def rhs(self) -> Waveform:
+        return self._rhs
+
+    @property
+    def arithmetic_operator(self) -> str:
+        return self._arithmetic_operator
+
+    @property
+    def duration(self) -> TimeType:
+        return self._lhs.duration
+
+    @property
+    def defined_channels(self) -> AbstractSet[ChannelID]:
+        return self._lhs.defined_channels | self._rhs.defined_channels
+
+    def unsafe_sample(self,
+                      channel: ChannelID,
+                      sample_times: np.ndarray,
+                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
+        if channel in self._lhs.defined_channels:
+            lhs = self._lhs.unsafe_sample(channel=channel, sample_times=sample_times, output_array=output_array)
+        else:
+            lhs = None
+
+        if channel in self._rhs.defined_channels:
+            rhs = self._rhs.unsafe_sample(channel=channel, sample_times=sample_times,
+                                          output_array=None if lhs is not None else output_array)
+        else:
+            rhs = None
+
+        if rhs is not None and lhs is not None:
+            arithmetic_operator = self.numpy_operator_map[self._arithmetic_operator]
+            if output_array is None:
+                output_array = lhs
+            return arithmetic_operator(lhs, rhs, out=output_array)
+
+        else:
+            if lhs is None:
+                assert rhs is not None, "channel %r not in defined channels (internal bug)" % channel
+                return self.numpy_rhs_only_map[self._arithmetic_operator](rhs, out=output_array)
+            else:
+                return lhs
+
+    def unsafe_get_subset_for_channels(self, channels: Set[ChannelID]) -> Waveform:
+        # TODO: optimization possible
+        return SubsetWaveform(self, channels)
+
+    @property
+    def compare_key(self) -> Tuple[str, Waveform, Waveform]:
+        return self._arithmetic_operator, self._lhs, self._rhs
+
+
+class FunctorWaveform(Waveform):
+    # TODO: Use Protocol to enforce that it accepts second argument has the keyword out
+    Functor = callable
+
+    __slots__ = ('_inner_waveform', '_functor')
+
+    """Apply a channel wise functor that works inplace to all results. The functor must accept two arguments"""
+    def __init__(self, inner_waveform: Waveform, functor: Mapping[ChannelID, Functor]):
+        super(FunctorWaveform, self).__init__(duration=inner_waveform.duration)
+        self._inner_waveform = inner_waveform
+        self._functor = dict(functor.items())
+
+        assert set(functor.keys()) == inner_waveform.defined_channels, ("There is no default identity mapping (yet)."
+                                                                        "File an issue on github if you need it.")
+
+    @classmethod
+    def from_functor(cls, inner_waveform: Waveform, functor: Mapping[ChannelID, Functor]):
+        constant_values = inner_waveform.constant_value_dict()
+        if constant_values is None:
+            return FunctorWaveform(inner_waveform, functor)
+
+        funced_constant_values = {ch: functor[ch](val) for ch, val in constant_values.items()}
+        return ConstantWaveform.from_mapping(inner_waveform.duration, funced_constant_values)
+
+    def is_constant(self) -> bool:
+        # only correct if `from_functor` was used
+        return False
+
+    def constant_value_dict(self) -> Optional[Mapping[ChannelID, float]]:
+        # only correct if `from_functor` was used
+        return None
+
+    def constant_value(self, channel: ChannelID) -> Optional[float]:
+        inner = self._inner_waveform.constant_value(channel)
+        if inner is None:
+            return None
+        else:
+            return self._functor[channel](inner)
+
+    @property
+    def defined_channels(self) -> AbstractSet[ChannelID]:
+        return self._inner_waveform.defined_channels
+
+    def unsafe_sample(self,
+                      channel: ChannelID,
+                      sample_times: np.ndarray,
+                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
+        inner_output = self._inner_waveform.unsafe_sample(channel, sample_times, output_array)
+        return self._functor[channel](inner_output, out=inner_output)
+
+    def unsafe_get_subset_for_channels(self, channels: Set[ChannelID]) -> Waveform:
+        return FunctorWaveform.from_functor(
+            self._inner_waveform.unsafe_get_subset_for_channels(channels),
+            {ch: self._functor[ch] for ch in channels})
+
+    @property
+    def compare_key(self) -> Tuple[Waveform, FrozenSet]:
+        return self._inner_waveform, frozenset(self._functor.items())
+
+
+class ReversedWaveform(Waveform):
+    """Reverses the inner waveform in time."""
+
+    __slots__ = ('_inner',)
+
+    def __init__(self, inner: Waveform):
+        super().__init__(duration=inner.duration)
+        self._inner = inner
+
+    @classmethod
+    def from_to_reverse(cls, inner: Waveform) -> Waveform:
+        if inner.constant_value_dict():
+            return inner
+        else:
+            return cls(inner)
+
+    def unsafe_sample(self, channel: ChannelID, sample_times: np.ndarray,
+                      output_array: Union[np.ndarray, None] = None) -> np.ndarray:
+        inner_sample_times = (float(self.duration) - sample_times)[::-1]
+        if output_array is None:
+            return self._inner.unsafe_sample(channel, inner_sample_times, None)[::-1]
+        else:
+            inner_output_array = output_array[::-1]
+            inner_output_array = self._inner.unsafe_sample(channel, inner_sample_times, output_array=inner_output_array)
+            if inner_output_array.base not in (output_array, output_array.base):
+                # TODO: is there a guarantee by numpy we never end up here?
+                output_array[:] = inner_output_array[::-1]
+            return output_array
+
+    @property
+    def defined_channels(self) -> AbstractSet[ChannelID]:
+        return self._inner.defined_channels
+
+    def unsafe_get_subset_for_channels(self, channels: AbstractSet[ChannelID]) -> 'Waveform':
+        return ReversedWaveform.from_to_reverse(self._inner.unsafe_get_subset_for_channels(channels))
+
+    @property
+    def compare_key(self) -> Hashable:
+        return self._inner.compare_key
+
+    def reversed(self) -> 'Waveform':
+        return self._inner
diff --git a/qupulse/pulses/arithmetic_pulse_template.py b/qupulse/pulses/arithmetic_pulse_template.py
index dc9066f8..46b61df4 100644
--- a/qupulse/pulses/arithmetic_pulse_template.py
+++ b/qupulse/pulses/arithmetic_pulse_template.py
@@ -14,8 +14,8 @@
 from qupulse.utils.types import ChannelID
 from qupulse.pulses.measurement import MeasurementWindow
 from qupulse.pulses.pulse_template import AtomicPulseTemplate, PulseTemplate
-from qupulse._program.waveforms import Waveform, ArithmeticWaveform, TransformingWaveform
-from qupulse._program.transformation import Transformation, ScalingTransformation, OffsetTransformation,\
+from qupulse.program.waveforms import Waveform, ArithmeticWaveform, TransformingWaveform
+from qupulse.program.transformation import Transformation, ScalingTransformation, OffsetTransformation,\
     IdentityTransformation
 
 
diff --git a/qupulse/pulses/constant_pulse_template.py b/qupulse/pulses/constant_pulse_template.py
index 24ae1684..192ee82b 100644
--- a/qupulse/pulses/constant_pulse_template.py
+++ b/qupulse/pulses/constant_pulse_template.py
@@ -10,7 +10,7 @@
 import numbers
 from typing import Any, Dict, List, Optional, Union, Mapping, AbstractSet
 
-from qupulse._program.waveforms import ConstantWaveform
+from qupulse.program.waveforms import ConstantWaveform
 from qupulse.utils.types import TimeType, ChannelID
 from qupulse.utils import cached_property
 from qupulse.expressions import ExpressionScalar, ExpressionLike
diff --git a/qupulse/pulses/function_pulse_template.py b/qupulse/pulses/function_pulse_template.py
index fa02feaf..24d98fbe 100644
--- a/qupulse/pulses/function_pulse_template.py
+++ b/qupulse/pulses/function_pulse_template.py
@@ -17,7 +17,7 @@
 from qupulse.utils.types import ChannelID, TimeType, time_from_float
 from qupulse.pulses.parameters import ParameterConstrainer, ParameterConstraint
 from qupulse.pulses.pulse_template import AtomicPulseTemplate, MeasurementDeclaration
-from qupulse._program.waveforms import FunctionWaveform
+from qupulse.program.waveforms import FunctionWaveform
 
 
 __all__ = ["FunctionPulseTemplate"]
diff --git a/qupulse/pulses/loop_pulse_template.py b/qupulse/pulses/loop_pulse_template.py
index 101705eb..72e8a602 100644
--- a/qupulse/pulses/loop_pulse_template.py
+++ b/qupulse/pulses/loop_pulse_template.py
@@ -13,13 +13,13 @@
 from qupulse.parameter_scope import Scope, MappedScope, DictScope
 from qupulse.utils.types import FrozenDict, FrozenMapping
 
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 
 from qupulse.expressions import ExpressionScalar, ExpressionVariableMissingException, Expression
 from qupulse.utils import checked_int_cast, cached_property
 from qupulse.pulses.parameters import InvalidParameterNameException, ParameterConstrainer, ParameterNotProvidedException
 from qupulse.pulses.pulse_template import PulseTemplate, ChannelID, AtomicPulseTemplate
-from qupulse._program.waveforms import SequenceWaveform as ForLoopWaveform
+from qupulse.program.waveforms import SequenceWaveform as ForLoopWaveform
 from qupulse.pulses.measurement import MeasurementDefiner, MeasurementDeclaration
 from qupulse.pulses.range import ParametrizedRange, RangeScope
 
diff --git a/qupulse/pulses/mapping_pulse_template.py b/qupulse/pulses/mapping_pulse_template.py
index 9134da25..c39aa846 100644
--- a/qupulse/pulses/mapping_pulse_template.py
+++ b/qupulse/pulses/mapping_pulse_template.py
@@ -8,8 +8,8 @@
 from qupulse.parameter_scope import Scope, MappedScope
 from qupulse.pulses.pulse_template import PulseTemplate, MappingTuple
 from qupulse.pulses.parameters import ParameterNotProvidedException, ParameterConstrainer
-from qupulse._program.waveforms import Waveform
-from qupulse._program._loop import Loop
+from qupulse.program.waveforms import Waveform
+from qupulse.program.loop import Loop
 from qupulse.serialization import Serializer, PulseRegistryType
 
 __all__ = [
diff --git a/qupulse/pulses/multi_channel_pulse_template.py b/qupulse/pulses/multi_channel_pulse_template.py
index beb12b68..8b19d152 100644
--- a/qupulse/pulses/multi_channel_pulse_template.py
+++ b/qupulse/pulses/multi_channel_pulse_template.py
@@ -17,8 +17,8 @@
 from qupulse.utils import isclose
 from qupulse.utils.sympy import almost_equal, Sympifyable
 from qupulse.utils.types import ChannelID, TimeType
-from qupulse._program.waveforms import MultiChannelWaveform, Waveform, TransformingWaveform
-from qupulse._program.transformation import ParallelChannelTransformation, Transformation, chain_transformations
+from qupulse.program.waveforms import MultiChannelWaveform, Waveform, TransformingWaveform
+from qupulse.program.transformation import ParallelChannelTransformation, Transformation, chain_transformations
 from qupulse.pulses.pulse_template import PulseTemplate, AtomicPulseTemplate
 from qupulse.pulses.mapping_pulse_template import MappingPulseTemplate, MappingTuple
 from qupulse.pulses.parameters import ParameterConstrainer
diff --git a/qupulse/pulses/point_pulse_template.py b/qupulse/pulses/point_pulse_template.py
index b5c5cfa0..5972e474 100644
--- a/qupulse/pulses/point_pulse_template.py
+++ b/qupulse/pulses/point_pulse_template.py
@@ -9,7 +9,7 @@
 from qupulse.utils.sympy import IndexedBroadcast
 from qupulse.utils.types import ChannelID
 from qupulse.expressions import Expression, ExpressionScalar
-from qupulse._program.waveforms import TableWaveform, TableWaveformEntry
+from qupulse.program.waveforms import TableWaveform, TableWaveformEntry
 from qupulse.pulses.parameters import ParameterConstraint, ParameterConstrainer
 from qupulse.pulses.pulse_template import AtomicPulseTemplate, MeasurementDeclaration
 from qupulse.pulses.table_pulse_template import TableEntry, EntryInInit
diff --git a/qupulse/pulses/pulse_template.py b/qupulse/pulses/pulse_template.py
index c8b598e8..57f756e2 100644
--- a/qupulse/pulses/pulse_template.py
+++ b/qupulse/pulses/pulse_template.py
@@ -19,10 +19,10 @@
 from qupulse.utils import forced_hash
 from qupulse.serialization import Serializable
 from qupulse.expressions import ExpressionScalar, Expression, ExpressionLike
-from qupulse._program._loop import Loop, to_waveform
-from qupulse._program.transformation import Transformation, IdentityTransformation, ChainedTransformation, chain_transformations
+from qupulse.program.loop import Loop, to_waveform
+from qupulse.program.transformation import Transformation, IdentityTransformation, ChainedTransformation, chain_transformations
 
-from qupulse._program.waveforms import Waveform, TransformingWaveform
+from qupulse.program.waveforms import Waveform, TransformingWaveform
 from qupulse.pulses.measurement import MeasurementDefiner, MeasurementDeclaration
 from qupulse.parameter_scope import Scope, DictScope
 
diff --git a/qupulse/pulses/repetition_pulse_template.py b/qupulse/pulses/repetition_pulse_template.py
index 809a91ce..674e4cdd 100644
--- a/qupulse/pulses/repetition_pulse_template.py
+++ b/qupulse/pulses/repetition_pulse_template.py
@@ -8,7 +8,8 @@
 import numpy as np
 
 from qupulse.serialization import Serializer, PulseRegistryType
-from qupulse._program._loop import Loop, VolatileRepetitionCount
+from qupulse.program.loop import Loop
+from qupulse.program.volatile import VolatileRepetitionCount
 from qupulse.parameter_scope import Scope
 
 from qupulse.utils.types import ChannelID
diff --git a/qupulse/pulses/sequence_pulse_template.py b/qupulse/pulses/sequence_pulse_template.py
index 1ef1bed0..7a8eafee 100644
--- a/qupulse/pulses/sequence_pulse_template.py
+++ b/qupulse/pulses/sequence_pulse_template.py
@@ -8,14 +8,14 @@
 import warnings
 
 from qupulse.serialization import Serializer, PulseRegistryType
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 from qupulse.parameter_scope import Scope
 from qupulse.utils import cached_property
 from qupulse.utils.types import MeasurementWindow, ChannelID, TimeType
 from qupulse.pulses.pulse_template import PulseTemplate, AtomicPulseTemplate
 from qupulse.pulses.parameters import ConstraintLike, ParameterConstrainer
 from qupulse.pulses.mapping_pulse_template import MappingPulseTemplate, MappingTuple
-from qupulse._program.waveforms import SequenceWaveform
+from qupulse.program.waveforms import SequenceWaveform
 from qupulse.pulses.measurement import MeasurementDeclaration, MeasurementDefiner
 from qupulse.expressions import Expression, ExpressionScalar
 
diff --git a/qupulse/pulses/table_pulse_template.py b/qupulse/pulses/table_pulse_template.py
index 873bbdf9..973536e8 100644
--- a/qupulse/pulses/table_pulse_template.py
+++ b/qupulse/pulses/table_pulse_template.py
@@ -23,7 +23,7 @@
 from qupulse.pulses.pulse_template import AtomicPulseTemplate, MeasurementDeclaration
 from qupulse.pulses.interpolation import InterpolationStrategy, LinearInterpolationStrategy, \
     HoldInterpolationStrategy, JumpInterpolationStrategy
-from qupulse._program.waveforms import TableWaveform, TableWaveformEntry
+from qupulse.program.waveforms import TableWaveform, TableWaveformEntry
 from qupulse.expressions import ExpressionScalar, Expression
 from qupulse.pulses.multi_channel_pulse_template import MultiChannelWaveform
 
diff --git a/qupulse/pulses/time_reversal_pulse_template.py b/qupulse/pulses/time_reversal_pulse_template.py
index a4758d1a..35a1884b 100644
--- a/qupulse/pulses/time_reversal_pulse_template.py
+++ b/qupulse/pulses/time_reversal_pulse_template.py
@@ -1,8 +1,8 @@
 from typing import Optional, Set, Dict, Union
 
 from qupulse import ChannelID
-from qupulse._program._loop import Loop
-from qupulse._program.waveforms import Waveform
+from qupulse.program.loop import Loop
+from qupulse.program.waveforms import Waveform
 from qupulse.serialization import PulseRegistryType
 from qupulse.expressions import ExpressionScalar
 


Test patch:
diff --git a/tests/_program/loop_tests.py b/tests/_program/loop_tests.py
index a3ff48aa..55d4e465 100644
--- a/tests/_program/loop_tests.py
+++ b/tests/_program/loop_tests.py
@@ -11,12 +11,11 @@
 from qupulse.parameter_scope import DictScope
 
 from qupulse.utils.types import TimeType, time_from_float
-from qupulse._program.volatile import VolatileRepetitionCount
-from qupulse._program._loop import Loop, _make_compatible, _is_compatible, _CompatibilityLevel,\
+from qupulse.program.volatile import VolatileRepetitionCount
+from qupulse.program.loop import Loop, _make_compatible, _is_compatible, _CompatibilityLevel,\
     RepetitionWaveform, SequenceWaveform, make_compatible, MakeCompatibleWarning, DroppedMeasurementWarning,\
     VolatileModificationWarning, roll_constant_waveforms
-from qupulse._program._loop import Loop, _make_compatible, _is_compatible, _CompatibilityLevel,\
-    RepetitionWaveform, SequenceWaveform, make_compatible, MakeCompatibleWarning, ConstantWaveform
+from qupulse.program.waveforms import *
 from tests.pulses.sequencing_dummies import DummyWaveform
 from qupulse.pulses.multi_channel_pulse_template import MultiChannelWaveform
 
@@ -488,27 +487,27 @@ def test_make_compatible(self):
                           sample_rate=TimeType.from_float(1.))
         priv_kwargs = dict(min_len=5, quantum=10, sample_rate=TimeType.from_float(1.))
 
-        with mock.patch('qupulse._program._loop._is_compatible',
+        with mock.patch('qupulse.program.loop._is_compatible',
                         return_value=_CompatibilityLevel.incompatible_too_short) as mocked:
             with self.assertRaisesRegex(ValueError, 'too short'):
                 make_compatible(program, **pub_kwargs)
             mocked.assert_called_once_with(program, **priv_kwargs)
 
-        with mock.patch('qupulse._program._loop._is_compatible',
+        with mock.patch('qupulse.program.loop._is_compatible',
                         return_value=_CompatibilityLevel.incompatible_fraction) as mocked:
             with self.assertRaisesRegex(ValueError, 'not an integer'):
                 make_compatible(program, **pub_kwargs)
             mocked.assert_called_once_with(program, **priv_kwargs)
 
-        with mock.patch('qupulse._program._loop._is_compatible',
+        with mock.patch('qupulse.program.loop._is_compatible',
                         return_value=_CompatibilityLevel.incompatible_quantum) as mocked:
             with self.assertRaisesRegex(ValueError, 'not a multiple of quantum'):
                 make_compatible(program, **pub_kwargs)
             mocked.assert_called_once_with(program, **priv_kwargs)
 
-        with mock.patch('qupulse._program._loop._is_compatible',
+        with mock.patch('qupulse.program.loop._is_compatible',
                         return_value=_CompatibilityLevel.action_required) as is_compat:
-            with mock.patch('qupulse._program._loop._make_compatible') as make_compat:
+            with mock.patch('qupulse.program.loop._make_compatible') as make_compat:
                 make_compatible(program, **pub_kwargs)
 
                 is_compat.assert_called_once_with(program, **priv_kwargs)
diff --git a/tests/_program/seqc_tests.py b/tests/_program/seqc_tests.py
index 729dacb1..75644fdf 100644
--- a/tests/_program/seqc_tests.py
+++ b/tests/_program/seqc_tests.py
@@ -13,7 +13,7 @@
 from qupulse.expressions import ExpressionScalar
 from qupulse.parameter_scope import DictScope
 
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 from qupulse._program.waveforms import ConstantWaveform
 from qupulse._program.seqc import BinaryWaveform, loop_to_seqc, WaveformPlayback, Repeat, SteppingRepeat, Scope,\
     to_node_clusters, find_sharable_waveforms, mark_sharable_waveforms, UserRegisterManager, HDAWGProgramManager,\
diff --git a/tests/_program/tabor_tests.py b/tests/_program/tabor_tests.py
index 9021d99a..919147ec 100644
--- a/tests/_program/tabor_tests.py
+++ b/tests/_program/tabor_tests.py
@@ -18,7 +18,7 @@
 
 from qupulse._program.tabor import TaborException, TaborProgram, find_place_for_segments_in_memory,\
     TaborSegment, TaborSequencing, PlottableProgram, TableDescription, make_combined_wave, TableEntry
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 from qupulse._program.volatile import VolatileRepetitionCount
 from qupulse.hardware.util import voltage_to_uint16
 from qupulse.utils.types import TimeType
diff --git a/tests/_program/transformation_tests.py b/tests/_program/transformation_tests.py
index 3f366482..f9422f1f 100644
--- a/tests/_program/transformation_tests.py
+++ b/tests/_program/transformation_tests.py
@@ -5,7 +5,7 @@
 
 
 from qupulse.expressions import ExpressionScalar
-from qupulse._program.transformation import LinearTransformation, Transformation, IdentityTransformation,\
+from qupulse.program.transformation import LinearTransformation, Transformation, IdentityTransformation,\
     ChainedTransformation, ParallelChannelTransformation, chain_transformations, OffsetTransformation,\
     ScalingTransformation
 
@@ -43,7 +43,7 @@ def test_chain(self):
 
         self.assertIs(trafo.chain(IdentityTransformation()), trafo)
 
-        with mock.patch('qupulse._program.transformation.chain_transformations',
+        with mock.patch('qupulse.program.transformation.chain_transformations',
                         return_value='asd') as chain_transformations:
             self.assertEqual(trafo.chain(trafo), 'asd')
             chain_transformations.assert_called_once_with(trafo, trafo)
@@ -277,7 +277,7 @@ def test_chain(self):
         trafo = TransformationStub()
         chained = ChainedTransformation(*trafos)
 
-        with mock.patch('qupulse._program.transformation.chain_transformations',
+        with mock.patch('qupulse.program.transformation.chain_transformations',
                         return_value='asd') as chain_transformations:
             self.assertEqual(chained.chain(trafo), 'asd')
             chain_transformations.assert_called_once_with(*trafos, trafo)
diff --git a/tests/_program/waveforms_tests.py b/tests/_program/waveforms_tests.py
index d84b2c76..c62fceb3 100644
--- a/tests/_program/waveforms_tests.py
+++ b/tests/_program/waveforms_tests.py
@@ -7,10 +7,10 @@
 from qupulse.utils.types import TimeType
 from qupulse.pulses.interpolation import HoldInterpolationStrategy, LinearInterpolationStrategy,\
     JumpInterpolationStrategy
-from qupulse._program.waveforms import MultiChannelWaveform, RepetitionWaveform, SequenceWaveform,\
+from qupulse.program.waveforms import MultiChannelWaveform, RepetitionWaveform, SequenceWaveform,\
     TableWaveformEntry, TableWaveform, TransformingWaveform, SubsetWaveform, ArithmeticWaveform, ConstantWaveform,\
     Waveform, FunctorWaveform, FunctionWaveform, ReversedWaveform
-from qupulse._program.transformation import LinearTransformation
+from qupulse.program.transformation import LinearTransformation
 from qupulse.expressions import ExpressionScalar, Expression
 
 from tests.pulses.sequencing_dummies import DummyWaveform, DummyInterpolationStrategy
diff --git a/tests/hardware/base_tests.py b/tests/hardware/base_tests.py
index a3f02961..db658573 100644
--- a/tests/hardware/base_tests.py
+++ b/tests/hardware/base_tests.py
@@ -5,7 +5,7 @@
 import numpy as np
 
 from qupulse.utils.types import TimeType
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 from qupulse.hardware.awgs.base import ProgramEntry
 
 from tests.pulses.sequencing_dummies import DummyWaveform
diff --git a/tests/hardware/feature_awg/awg_new_driver_base_tests.py b/tests/hardware/feature_awg/awg_new_driver_base_tests.py
index 12070e1c..1845ab1e 100644
--- a/tests/hardware/feature_awg/awg_new_driver_base_tests.py
+++ b/tests/hardware/feature_awg/awg_new_driver_base_tests.py
@@ -3,7 +3,7 @@
 import warnings
 
 from qupulse import ChannelID
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 from qupulse.hardware.feature_awg import channel_tuple_wrapper
 from qupulse.hardware.feature_awg.base import AWGDevice, AWGChannel, AWGChannelTuple, AWGMarkerChannel
 from qupulse.hardware.feature_awg.features import ChannelSynchronization, ProgramManagement, VoltageRange, \
diff --git a/tests/hardware/setup_tests.py b/tests/hardware/setup_tests.py
index d36467e2..fd3ba3c1 100644
--- a/tests/hardware/setup_tests.py
+++ b/tests/hardware/setup_tests.py
@@ -4,7 +4,7 @@
 import numpy as np
 
 from qupulse.hardware.setup import HardwareSetup, PlaybackChannel, MarkerChannel, MeasurementMask
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 
 from tests.pulses.sequencing_dummies import DummyWaveform
 
diff --git a/tests/hardware/tabor_dummy_based_tests.py b/tests/hardware/tabor_dummy_based_tests.py
index 7a49e0b6..98588b5e 100644
--- a/tests/hardware/tabor_dummy_based_tests.py
+++ b/tests/hardware/tabor_dummy_based_tests.py
@@ -10,7 +10,7 @@
 
 from qupulse.hardware.awgs.base import AWGAmplitudeOffsetHandling
 from qupulse.utils.types import TimeType
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 from qupulse._program.tabor import TableDescription, TimeType, TableEntry, TaborSegment, TaborProgram,\
     make_combined_wave, TaborSequencing
 from qupulse._program.waveforms import ConstantWaveform
diff --git a/tests/hardware/tabor_tests.py b/tests/hardware/tabor_tests.py
index 6226b2c6..e9bf6e5b 100644
--- a/tests/hardware/tabor_tests.py
+++ b/tests/hardware/tabor_tests.py
@@ -9,7 +9,7 @@
 
 from qupulse.hardware.awgs.tabor import TaborException, TaborProgram, \
     TaborSegment, TaborSequencing, with_configuration_guard, PlottableProgram
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 from qupulse.hardware.util import voltage_to_uint16
 
 from tests.pulses.sequencing_dummies import DummyWaveform
diff --git a/tests/hardware/tektronix_tests.py b/tests/hardware/tektronix_tests.py
index 773876e2..aed67d79 100644
--- a/tests/hardware/tektronix_tests.py
+++ b/tests/hardware/tektronix_tests.py
@@ -13,7 +13,7 @@
 
 from qupulse.hardware.awgs.tektronix import TektronixAWG, TektronixProgram, parse_program, _make_binary_waveform,\
     voltage_to_uint16, WaveformEntry, WaveformStorage
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 from qupulse.utils.types import TimeType
 from tests.pulses.sequencing_dummies import DummyWaveform
 from qupulse._program.waveforms import MultiChannelWaveform
diff --git a/tests/hardware/zihdawg_tests.py b/tests/hardware/zihdawg_tests.py
index 9fe675ea..2cb6bbaf 100644
--- a/tests/hardware/zihdawg_tests.py
+++ b/tests/hardware/zihdawg_tests.py
@@ -26,7 +26,7 @@
         raise unittest.SkipTest("zhinst not present") from err
 
 from qupulse.utils.types import TimeType
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 from tests.pulses.sequencing_dummies import DummyWaveform
 from qupulse.hardware.awgs.zihdawg import HDAWGChannelGroup, HDAWGRepresentation, HDAWGValueError, UserRegister,\
     ELFManager, HDAWGChannelGrouping, SingleDeviceChannelGroup
diff --git a/tests/pulses/bug_tests.py b/tests/pulses/bug_tests.py
index 6be11fcf..44fd4967 100644
--- a/tests/pulses/bug_tests.py
+++ b/tests/pulses/bug_tests.py
@@ -13,7 +13,7 @@
 
 from qupulse.plotting import plot
 
-from qupulse._program._loop import to_waveform
+from qupulse.program.loop import to_waveform
 from qupulse.utils import isclose
 
 class BugTests(unittest.TestCase):
@@ -77,7 +77,7 @@ def test_issue_584_uninitialized_table_sample(self):
         """issue 584"""
         d = 598.3333333333334 - 480
         tpt = TablePulseTemplate(entries={'P': [(0, 1.0, 'hold'), (d, 1.0, 'hold')]})
-        with mock.patch('qupulse._program.waveforms.PULSE_TO_WAVEFORM_ERROR', 1e-6):
+        with mock.patch('qupulse.program.waveforms.PULSE_TO_WAVEFORM_ERROR', 1e-6):
             wf = to_waveform(tpt.create_program())
             self.assertTrue(isclose(d, wf.duration, abs_tol=1e-6))
 
diff --git a/tests/pulses/constant_pulse_template_tests.py b/tests/pulses/constant_pulse_template_tests.py
index f0e37486..197aa101 100644
--- a/tests/pulses/constant_pulse_template_tests.py
+++ b/tests/pulses/constant_pulse_template_tests.py
@@ -1,14 +1,14 @@
 import unittest
 
 import qupulse.plotting
-import qupulse._program.waveforms
+import qupulse.program.waveforms
 import qupulse.utils.sympy
 from qupulse.pulses import TablePT, FunctionPT, AtomicMultiChannelPT, MappingPT
 from qupulse.pulses.multi_channel_pulse_template import AtomicMultiChannelPulseTemplate
 from qupulse.plotting import plot
 from qupulse.pulses.sequence_pulse_template import SequencePulseTemplate
-from qupulse._program._loop import make_compatible
-from qupulse._program.waveforms import ConstantWaveform
+from qupulse.program.loop import make_compatible
+from qupulse.program.waveforms import ConstantWaveform
 
 from qupulse.serialization import DictBackend, PulseStorage
 from qupulse.pulses.constant_pulse_template import ConstantPulseTemplate, ExpressionScalar, TimeType
@@ -44,10 +44,10 @@ def test_zero_duration(self):
         self.assertEqual(pulse.duration, 12)
 
     def test_regression_duration_conversion(self):
-        old_value = qupulse._program.waveforms.PULSE_TO_WAVEFORM_ERROR
+        old_value = qupulse.program.waveforms.PULSE_TO_WAVEFORM_ERROR
 
         try:
-            qupulse._program.waveforms.PULSE_TO_WAVEFORM_ERROR = 1e-6
+            qupulse.program.waveforms.PULSE_TO_WAVEFORM_ERROR = 1e-6
             for duration_in_samples in [64, 936320, 24615392]:
                 p = ConstantPulseTemplate(duration_in_samples / 2.4, {'a': 0})
                 number_of_samples = p.create_program().duration * 2.4
@@ -57,19 +57,19 @@ def test_regression_duration_conversion(self):
                 p2 = ConstantPulseTemplate((duration_in_samples + 1) / 2.4, {'a': 0})
                 self.assertNotEqual(p.create_program().duration, p2.create_program().duration)
         finally:
-            qupulse._program.waveforms.PULSE_TO_WAVEFORM_ERROR = old_value
+            qupulse.program.waveforms.PULSE_TO_WAVEFORM_ERROR = old_value
 
     def test_regression_duration_conversion_functionpt(self):
-        old_value = qupulse._program.waveforms.PULSE_TO_WAVEFORM_ERROR
+        old_value = qupulse.program.waveforms.PULSE_TO_WAVEFORM_ERROR
 
         try:
-            qupulse._program.waveforms.PULSE_TO_WAVEFORM_ERROR = 1e-6
+            qupulse.program.waveforms.PULSE_TO_WAVEFORM_ERROR = 1e-6
             for duration_in_samples in [64, 2000, 936320]:
                 p = FunctionPT('1', duration_expression=duration_in_samples / 2.4, channel='a')
                 number_of_samples = p.create_program().duration * 2.4
                 self.assertEqual(number_of_samples.denominator, 1)
         finally:
-            qupulse._program.waveforms.PULSE_TO_WAVEFORM_ERROR = old_value
+            qupulse.program.waveforms.PULSE_TO_WAVEFORM_ERROR = old_value
 
     def test_regression_template_combination(self):
         old_value = qupulse.utils.sympy.SYMPY_DURATION_ERROR_MARGIN
diff --git a/tests/pulses/loop_pulse_template_tests.py b/tests/pulses/loop_pulse_template_tests.py
index 0df6a2a0..e5380607 100644
--- a/tests/pulses/loop_pulse_template_tests.py
+++ b/tests/pulses/loop_pulse_template_tests.py
@@ -10,7 +10,7 @@
 from qupulse.pulses.parameters import InvalidParameterNameException, ParameterConstraintViolation,\
     ParameterNotProvidedException, ParameterConstraint
 
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 
 from tests.pulses.sequencing_dummies import DummyPulseTemplate, MeasurementWindowTestCase, DummyWaveform
 from tests.serialization_dummies import DummySerializer
diff --git a/tests/pulses/mapping_pulse_template_tests.py b/tests/pulses/mapping_pulse_template_tests.py
index 366360d9..b3da1174 100644
--- a/tests/pulses/mapping_pulse_template_tests.py
+++ b/tests/pulses/mapping_pulse_template_tests.py
@@ -8,7 +8,7 @@
     AmbiguousMappingException, MappingCollisionException
 from qupulse.pulses.parameters import ParameterConstraintViolation, ParameterConstraint, ParameterNotProvidedException
 from qupulse.expressions import Expression
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 
 from tests.pulses.sequencing_dummies import DummyPulseTemplate, MeasurementWindowTestCase, DummyWaveform
 from tests.serialization_tests import SerializableTests
diff --git a/tests/pulses/plotting_tests.py b/tests/pulses/plotting_tests.py
index 0e96a2c9..fdceac21 100644
--- a/tests/pulses/plotting_tests.py
+++ b/tests/pulses/plotting_tests.py
@@ -9,7 +9,7 @@
 from qupulse.plotting import PlottingNotPossibleException, render, plot
 from qupulse.pulses.table_pulse_template import TablePulseTemplate
 from qupulse.pulses.sequence_pulse_template import SequencePulseTemplate
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 
 from tests.pulses.sequencing_dummies import DummyWaveform, DummyPulseTemplate
 
diff --git a/tests/pulses/pulse_template_tests.py b/tests/pulses/pulse_template_tests.py
index 64e5451d..f8c0e241 100644
--- a/tests/pulses/pulse_template_tests.py
+++ b/tests/pulses/pulse_template_tests.py
@@ -12,7 +12,7 @@
     TimeReversalPT, AtomicMultiChannelPT
 from qupulse.pulses.pulse_template import AtomicPulseTemplate, PulseTemplate, UnknownVolatileParameter
 from qupulse.pulses.multi_channel_pulse_template import MultiChannelWaveform
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 
 from qupulse._program.transformation import Transformation
 from qupulse._program.waveforms import TransformingWaveform
diff --git a/tests/pulses/repetition_pulse_template_tests.py b/tests/pulses/repetition_pulse_template_tests.py
index b59cef82..5f83fc9d 100644
--- a/tests/pulses/repetition_pulse_template_tests.py
+++ b/tests/pulses/repetition_pulse_template_tests.py
@@ -5,7 +5,7 @@
 from qupulse.parameter_scope import Scope, DictScope
 from qupulse.utils.types import FrozenDict
 
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 from qupulse.expressions import Expression, ExpressionScalar
 from qupulse.pulses import ConstantPT
 from qupulse.pulses.repetition_pulse_template import RepetitionPulseTemplate,ParameterNotIntegerException
diff --git a/tests/pulses/sequence_pulse_template_tests.py b/tests/pulses/sequence_pulse_template_tests.py
index 4985fb8d..3a625eb3 100644
--- a/tests/pulses/sequence_pulse_template_tests.py
+++ b/tests/pulses/sequence_pulse_template_tests.py
@@ -7,7 +7,7 @@
 from qupulse.pulses.sequence_pulse_template import SequencePulseTemplate, SequenceWaveform
 from qupulse.pulses.mapping_pulse_template import MappingPulseTemplate
 from qupulse.pulses.parameters import ParameterConstraint, ParameterConstraintViolation, ParameterNotProvidedException
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 
 from tests.pulses.sequencing_dummies import DummyPulseTemplate, DummyWaveform, MeasurementWindowTestCase
 from tests.serialization_dummies import DummySerializer
diff --git a/tests/pulses/sequencing_dummies.py b/tests/pulses/sequencing_dummies.py
index 889aea57..18ca487c 100644
--- a/tests/pulses/sequencing_dummies.py
+++ b/tests/pulses/sequencing_dummies.py
@@ -8,7 +8,7 @@
 
 """LOCAL IMPORTS"""
 from qupulse.parameter_scope import Scope
-from qupulse._program._loop import Loop
+from qupulse.program.loop import Loop
 from qupulse.utils.types import MeasurementWindow, ChannelID, TimeType, time_from_float
 from qupulse.serialization import Serializer
 from qupulse._program.waveforms import Waveform
diff --git a/tests/qctoolkit_alias_tests.py b/tests/qctoolkit_alias_tests.py
index 5c7a5d83..1ae44677 100644
--- a/tests/qctoolkit_alias_tests.py
+++ b/tests/qctoolkit_alias_tests.py
@@ -10,7 +10,7 @@ def test_alias(self):
         self.assertIs(qctoolkit.pulses.TablePT, qupulse.pulses.TablePT)
 
     def test_class_identity(self):
-        from qupulse._program._loop import Loop as Loop_qu
-        from qctoolkit._program._loop import Loop as Loop_qc
+        from qupulse.program.loop import Loop as Loop_qu
+        from qctoolkit.program.loop import Loop as Loop_qc
 
         self.assertIs(Loop_qc, Loop_qu)