Repository: TeamGraphix/graphix

Issue:
Remove `vop` from the `M` commands
**Is your feature request related to a problem? Please describe.**
(to be worked on after merging #155): `vop` info is redundant as it can be converted to measurement planes and angle changes.
Modify them instead of adding`vop`, in e.g. `measure_pauli` `traspile(opt=True)`, etc. remove `vop` from measurement command definition.

we may deprecate `opt=True` transpile so we may want to do this after that.

Discussion:
shinich1: @thierry-martinez thanks! this looks really nice approach. examples rendered well. we might want to add `Clifford` to API docs page, `docs/source/data.rst` (I can do this after this is merged).
Can @EarlMilktea also take a quick look?
shinich1: @thierry-martinez I think we can merge #182 before this, so that ruff will be added to checks (resolving that may auto-resolve some of above comments by @EarlMilktea).
EarlMilktea: @thierry-martinez Could you resolve the remaining three threads? Once ready, I will approve merge.
thierry-martinez: Oops, sorry for having missed them (they were folded in "3 hidden more conversations")! All of them should be fixed now. Thanks!

Patch:
diff --git a/graphix/command.py b/graphix/command.py
index e7bdf42..fcd2fbd 100644
--- a/graphix/command.py
+++ b/graphix/command.py
@@ -5,8 +5,11 @@ from __future__ import annotations
 import abc
 import enum
 
+import numpy as np
 from pydantic import BaseModel
 
+import graphix.clifford
+from graphix.clifford import Clifford
 from graphix.pauli import Plane
 
 Node = int
@@ -51,7 +54,29 @@ class M(Command):
     angle: float = 0.0
     s_domain: list[Node] = []
     t_domain: list[Node] = []
-    vop: int = 0
+
+    def clifford(self, clifford: Clifford) -> M:
+        s_domain = self.s_domain
+        t_domain = self.t_domain
+        for gate in clifford.hsz:
+            if gate == graphix.clifford.I:
+                pass
+            elif gate == graphix.clifford.H:
+                t_domain, s_domain = s_domain, t_domain
+            elif gate == graphix.clifford.S:
+                t_domain = s_domain + t_domain
+            elif gate == graphix.clifford.Z:
+                pass
+            else:
+                raise RuntimeError(f"{gate} should be either I, H, S or Z.")
+        update = graphix.pauli.MeasureUpdate.compute(self.plane, False, False, clifford)
+        return M(
+            node=self.node,
+            plane=update.new_plane,
+            angle=self.angle * update.coeff + update.add_term / np.pi,
+            s_domain=s_domain,
+            t_domain=t_domain,
+        )
 
 
 class E(Command):
diff --git a/graphix/pattern.py b/graphix/pattern.py
index 8f5b785..5c3dbe0 100644
--- a/graphix/pattern.py
+++ b/graphix/pattern.py
@@ -5,9 +5,11 @@ ref: V. Danos, E. Kashefi and P. Panangaden. J. ACM 54.2 8 (2007)
 from __future__ import annotations
 
 from copy import deepcopy
+from dataclasses import dataclass
 
 import networkx as nx
 import numpy as np
+import typing_extensions
 
 import graphix.clifford
 import graphix.pauli
@@ -16,6 +18,7 @@ from graphix.clifford import CLIFFORD_CONJ, CLIFFORD_MEASURE, CLIFFORD_TO_QASM3
 from graphix.device_interface import PatternRunner
 from graphix.gflow import find_flow, find_gflow, get_layers
 from graphix.graphsim.graphstate import GraphState
+from graphix.pauli import Plane
 from graphix.simulator import PatternSimulator
 from graphix.visualization import GraphVisualizer
 
@@ -246,7 +249,7 @@ class Pattern:
                 count += 1
                 print(
                     f"M, node = {cmd.node}, plane = {cmd.plane}, angle(pi) = {cmd.angle}, "
-                    + f"s_domain = {cmd.s_domain}, t_domain = {cmd.t_domain}, Clifford index = {cmd.vop}"
+                    + f"s_domain = {cmd.s_domain}, t_domain = {cmd.t_domain}"
                 )
             elif cmd.kind == command.CommandKind.X and (command.CommandKind.X in target):
                 count += 1
@@ -288,7 +291,7 @@ class Pattern:
         def fresh_node():
             return {
                 "seq": [],
-                "Mprop": [None, None, [], [], 0],
+                "Mprop": [None, None, [], []],
                 "Xsignal": [],
                 "Xsignals": [],
                 "Zsignal": [],
@@ -306,7 +309,7 @@ class Pattern:
                 node_prop[cmd.nodes[1]]["seq"].append(cmd.nodes[0])
                 node_prop[cmd.nodes[0]]["seq"].append(cmd.nodes[1])
             elif kind == command.CommandKind.M:
-                node_prop[cmd.node]["Mprop"] = [cmd.plane, cmd.angle, cmd.s_domain, cmd.t_domain, cmd.vop]
+                node_prop[cmd.node]["Mprop"] = [cmd.plane, cmd.angle, cmd.s_domain, cmd.t_domain]
                 node_prop[cmd.node]["seq"].append(-1)
                 morder.append(cmd.node)
             elif kind == command.CommandKind.X:
@@ -383,7 +386,7 @@ class Pattern:
         except StopIteration:
             return True
 
-    def shift_signals(self, method="local"):
+    def shift_signals(self, method="local") -> dict[int, list[int]]:
         """Performs signal shifting procedure
         Extract the t-dependence of the measurement into 'S' commands
         and commute them to the end of the command sequence where it can be removed.
@@ -400,13 +403,19 @@ class Pattern:
             'global' shift_signals is executed on a conventional Pattern sequence.
             'local' shift_signals is done on a LocalPattern class which is faster but results in equivalent pattern.
             defaults to 'local'
+
+        Returns
+        -------
+        swapped_dict : dict[int, list[int]]
+            for each node, the signal that have been shifted if the outcome is
+            swapped by the shift.
         """
         if method == "local":
             localpattern = self.get_local_pattern()
-            localpattern.shift_signals()
+            swapped_dict = localpattern.shift_signals()
             self.__seq = localpattern.get_pattern().__seq
         elif method == "global":
-            self.extract_signals()
+            swapped_dict = self.extract_signals()
             target = self._find_op_to_be_moved(command.CommandKind.S, rev=True)
             while target is not None:
                 if target == len(self.__seq) - 1:
@@ -428,6 +437,7 @@ class Pattern:
                 target += 1
         else:
             raise ValueError("Invalid method")
+        return swapped_dict
 
     def _find_op_to_be_moved(self, op: command.CommandKind, rev=False, skipnum=0):
         """Internal method for pattern modification.
@@ -498,11 +508,7 @@ class Pattern:
         X = self.__seq[target]
         M = self.__seq[target + 1]
         if X.node == M.node:
-            vop = M.vop
-            if M.plane == graphix.pauli.Plane.YZ or vop == 6:
-                M.t_domain.extend(X.domain)
-            elif M.plane == graphix.pauli.Plane.XY:
-                M.s_domain.extend(X.domain)
+            M.s_domain.extend(X.domain)
             self.__seq.pop(target)  # del X
             return True
         else:
@@ -523,11 +529,7 @@ class Pattern:
         Z = self.__seq[target]
         M = self.__seq[target + 1]
         if Z.node == M.node:
-            vop = M.vop
-            if M.plane == graphix.pauli.Plane.YZ or vop == 6:
-                M.s_domain.extend(Z.domain)
-            elif M.plane == graphix.pauli.Plane.XY:
-                M.t_domain.extend(Z.domain)
+            M.t_domain.extend(Z.domain)
             self.__seq.pop(target)  # del Z
             return True
         else:
@@ -710,22 +712,25 @@ class Pattern:
             self._commute_with_preceding(target)
             target -= 1
 
-    def extract_signals(self):
+    def extract_signals(self) -> dict[int, list[int]]:
         """Extracts 't' domain of measurement commands, turn them into
         signal 'S' commands and add to the command sequence.
         This is used for shift_signals() method.
         """
+        signal_dict = {}
         pos = 0
         while pos < len(self.__seq):
             if self.__seq[pos].kind == command.CommandKind.M:
                 cmd: command.M = self.__seq[pos]
-                if cmd.plane == graphix.pauli.Plane.XY:
-                    node = cmd.node
-                    if cmd.t_domain:
-                        self.__seq.insert(pos + 1, command.S(node=node, domain=cmd.t_domain))
-                        cmd.t_domain = []
-                        pos += 1
+                extracted_signal = extract_signal(cmd.plane, cmd.s_domain, cmd.t_domain)
+                if extracted_signal.signal:
+                    self.__seq.insert(pos + 1, command.S(node=cmd.node, domain=extracted_signal.signal))
+                    cmd.s_domain = extracted_signal.s_domain
+                    cmd.t_domain = extracted_signal.t_domain
+                    pos += 1
+                signal_dict[cmd.node] = extracted_signal.signal
             pos += 1
+        return signal_dict
 
     def _get_dependency(self):
         """Get dependency (byproduct correction & dependent measurement)
@@ -987,10 +992,7 @@ class Pattern:
         meas_plane = dict()
         for cmd in self.__seq:
             if cmd.kind == command.CommandKind.M:
-                mplane = cmd.plane
-                cliff = graphix.clifford.get(cmd.vop)
-                new_axes = [cliff.measure(graphix.pauli.Pauli.from_axis(axis)).axis for axis in mplane.axes]
-                meas_plane[cmd.node] = graphix.pauli.Plane.from_axes(*new_axes)
+                meas_plane[cmd.node] = cmd.plane
         return meas_plane
 
     def get_angles(self):
@@ -1077,14 +1079,8 @@ class Pattern:
         vops = dict()
         for cmd in self.__seq:
             if cmd.kind == command.CommandKind.M:
-                if cmd.vop == 0:
-                    if include_identity:
-                        vops[cmd.node] = cmd.vop
-                else:
-                    if conj:
-                        vops[cmd.node] = CLIFFORD_CONJ[cmd.vop]
-                    else:
-                        vops[cmd.node] = cmd.vop
+                if include_identity:
+                    vops[cmd.node] = cmd.vop
             elif cmd.kind == command.CommandKind.C:
                 if cmd.cliff_index == 0:
                     if include_identity:
@@ -1377,7 +1373,7 @@ class Pattern:
 
         if flow_from_pattern:
             vis.visualize_from_pattern(
-                pattern=deepcopy(self),
+                pattern=self.copy(),
                 show_pauli_measurement=show_pauli_measurement,
                 show_local_clifford=show_local_clifford,
                 show_measurement_planes=show_measurement_planes,
@@ -1422,6 +1418,16 @@ class Pattern:
                 for line in cmd_to_qasm3(cmd):
                     file.write(line)
 
+    def copy(self) -> Pattern:
+        result = self.__new__(self.__class__)
+        result.__seq = [cmd.model_copy() for cmd in self.__seq]
+        result.__input_nodes = self.__input_nodes.copy()
+        result.__output_nodes = self.__output_nodes.copy()
+        result.__Nnode = self.__Nnode
+        result._pauli_preprocessed = self._pauli_preprocessed
+        result.results = self.results.copy()
+        return result
+
 
 class CommandNode:
     """A node decorated with a distributed command sequence.
@@ -1439,7 +1445,7 @@ class CommandNode:
         Z: -3
         C: -4
     Mprop : list
-        attributes for a measurement command. consists of [meas_plane, angle, s_domain, t_domain, vop]
+        attributes for a measurement command. consists of [meas_plane, angle, s_domain, t_domain]
     result : int
         measurement result of the node
     Xsignal : list
@@ -1448,8 +1454,6 @@ class CommandNode:
         signal domain. Xsignals may contains lists. For standardization, this variable is used.
     Zsignal : list
         signal domain
-    vop : int
-        value for clifford index
     input : bool
         whether the node is an input or not
     output : bool
@@ -1496,7 +1500,6 @@ class CommandNode:
         self.Xsignal = Xsignal
         self.Xsignals = Xsignals
         self.Zsignal = Zsignal  # appeared at most e + 1
-        self.vop = Mprop[4] if len(Mprop) == 5 else 0
         self.input = is_input
         self.output = is_output
 
@@ -1549,10 +1552,7 @@ class CommandNode:
             self.Xsignal = combined_Xsignal
             self.Xsignals = [combined_Xsignal]
         else:
-            if self.Mprop[0] == graphix.pauli.Plane.YZ or self.vop == 6:
-                self.Mprop[3] = xor_combination_list(combined_Xsignal, self.Mprop[3])
-            elif self.Mprop[0] == graphix.pauli.Plane.XY:
-                self.Mprop[2] = xor_combination_list(combined_Xsignal, self.Mprop[2])
+            self.Mprop[2] = xor_combination_list(combined_Xsignal, self.Mprop[2])
             self.Xsignal = []
             self.Xsignals = []
         return EXcommutated_nodes
@@ -1566,10 +1566,7 @@ class CommandNode:
         if self.output and z_in_seq:
             self.seq.append(-3)
         else:
-            if self.Mprop[0] == graphix.pauli.Plane.YZ or self.vop == 6:
-                self.Mprop[2] = xor_combination_list(self.Zsignal, self.Mprop[2])
-            elif self.Mprop[0] == graphix.pauli.Plane.XY:
-                self.Mprop[3] = xor_combination_list(self.Zsignal, self.Mprop[3])
+            self.Mprop[3] = xor_combination_list(self.Zsignal, self.Mprop[3])
             self.Zsignal = []
 
     def _add_Z(self, pair, signal):
@@ -1619,7 +1616,6 @@ class CommandNode:
                 angle=self.Mprop[1],
                 s_domain=self.Mprop[2],
                 t_domain=self.Mprop[3],
-                vop=self.Mprop[4],
             )
         elif cmd == -2:
             if self.seq.count(-2) > 1:
@@ -1753,12 +1749,19 @@ class LocalPattern:
             for dependent_node in dependent_node_dicts["Z"]:
                 self.signal_destination[dependent_node]["Z"] |= {index}
 
-    def shift_signals(self):
+    def shift_signals(self) -> dict[int, list[int]]:
         """Shift signals to the back based on signal destinations."""
         self.collect_signal_destination()
+        signal_dict = {}
         for node_index in self.morder + self.output_nodes:
-            signal = self.nodes[node_index].Mprop[3]
-            self.nodes[node_index].Mprop[3] = []
+            node = self.nodes[node_index]
+            if node.Mprop[0] is None:
+                continue
+            extracted_signal = extract_signal(node.Mprop[0], node.Mprop[2], node.Mprop[3])
+            signal_dict[node_index] = extracted_signal.signal
+            signal = extracted_signal.signal
+            self.nodes[node_index].Mprop[2] = extracted_signal.s_domain
+            self.nodes[node_index].Mprop[3] = extracted_signal.t_domain
             for signal_label, destinated_nodes in self.signal_destination[node_index].items():
                 for destinated_node in destinated_nodes:
                     node = self.nodes[destinated_node]
@@ -1772,6 +1775,7 @@ class LocalPattern:
                         node.Zsignal += signal
                     else:
                         raise ValueError(f"Invalid signal label: {signal_label}")
+        return signal_dict
 
     def get_graph(self):
         """Get a graph from a local pattern
@@ -1956,10 +1960,7 @@ def measure_pauli(pattern, leave_input, copy=False, use_rustworkx=False):
     for cmd in pattern:
         if cmd.kind == command.CommandKind.M:
             if cmd.node in graph_state.nodes:
-                cmd_new = deepcopy(cmd)
-                new_clifford_ = vops[cmd.node]
-                cmd_new.vop = new_clifford_
-                new_seq.append(cmd_new)
+                new_seq.append(cmd.clifford(graphix.clifford.get(vops[cmd.node])))
     for index in pattern.output_nodes:
         new_clifford_ = vops[index]
         if new_clifford_ != 0:
@@ -2192,3 +2193,29 @@ def assert_permutation(original, user):
             node_set.remove(node)
         else:
             raise ValueError(f"{node} appears twice")
+
+
+@dataclass
+class ExtractedSignal:
+    """
+    Return data structure for `extract_signal`.
+    """
+
+    s_domain: list[int]
+    "New `s_domain` for the measure command."
+
+    t_domain: list[int]
+    "New `t_domain` for the measure command."
+
+    signal: list[int]
+    "Domain for the shift command."
+
+
+def extract_signal(plane: Plane, s_domain: list[int], t_domain: list[int]) -> ExtractedSignal:
+    if plane == Plane.XY:
+        return ExtractedSignal(s_domain=s_domain, t_domain=[], signal=t_domain)
+    if plane == Plane.XZ:
+        return ExtractedSignal(s_domain=[], t_domain=s_domain + t_domain, signal=s_domain)
+    if plane == Plane.YZ:
+        return ExtractedSignal(s_domain=[], t_domain=t_domain, signal=s_domain)
+    typing_extensions.assert_never(plane)
diff --git a/graphix/sim/base_backend.py b/graphix/sim/base_backend.py
index 57f00d5..0b9838c 100644
--- a/graphix/sim/base_backend.py
+++ b/graphix/sim/base_backend.py
@@ -22,19 +22,19 @@ def perform_measure(
     if pr_calc:
         op_mat = op_mat_from_result(vec, False)
         prob_0 = state.expectation_single(op_mat, qubit)
-        result = rng.rand() > prob_0
+        result = rng.random() > prob_0
         if result:
             op_mat = op_mat_from_result(vec, True)
     else:
         # choose the measurement result randomly
-        result = np.random.choice([0, 1])
+        result = rng.choice([0, 1])
         op_mat = op_mat_from_result(vec, result)
     state.evolve_single(op_mat, qubit)
     return result
 
 
 class Backend:
-    def __init__(self, pr_calc: bool = True):
+    def __init__(self, pr_calc: bool = True, rng: np.random.Generator | None = None):
         """
         Parameters
         ----------
@@ -43,19 +43,22 @@ class Backend:
                 if False, measurements yield results 0/1 with 50% probabilities each.
         """
         # whether to compute the probability
+        if rng is None:
+            self.__rng = np.random.default_rng()
+        else:
+            self.__rng = rng
         self.pr_calc = pr_calc
 
     def _perform_measure(self, cmd: graphix.command.M):
         s_signal = np.sum([self.results[j] for j in cmd.s_domain])
         t_signal = np.sum([self.results[j] for j in cmd.t_domain])
         angle = cmd.angle * np.pi
-        vop = cmd.vop
         measure_update = graphix.pauli.MeasureUpdate.compute(
-            cmd.plane, s_signal % 2 == 1, t_signal % 2 == 1, graphix.clifford.TABLE[vop]
+            cmd.plane, s_signal % 2 == 1, t_signal % 2 == 1, graphix.clifford.I
         )
         angle = angle * measure_update.coeff + measure_update.add_term
         loc = self.node_index.index(cmd.node)
-        result = perform_measure(loc, measure_update.new_plane, angle, self.state, np.random, self.pr_calc)
+        result = perform_measure(loc, measure_update.new_plane, angle, self.state, self.__rng, self.pr_calc)
         self.results[cmd.node] = result
         self.node_index.remove(cmd.node)
         return loc
diff --git a/graphix/sim/statevec.py b/graphix/sim/statevec.py
index 3b2652e..f45bef9 100644
--- a/graphix/sim/statevec.py
+++ b/graphix/sim/statevec.py
@@ -27,6 +27,7 @@ class StatevectorBackend(graphix.sim.base_backend.Backend):
         input_state: Data = graphix.states.BasicStates.PLUS,
         max_qubit_num=20,
         pr_calc=True,
+        rng=None,
     ):
         """
         Parameters
@@ -57,7 +58,7 @@ class StatevectorBackend(graphix.sim.base_backend.Backend):
         self.max_qubit_num = max_qubit_num
         if pattern.max_space() > max_qubit_num:
             raise ValueError("Pattern.max_space is larger than max_qubit_num. Increase max_qubit_num and try again.")
-        super().__init__(pr_calc)
+        super().__init__(pr_calc, rng)
 
         # initialize input qubits to desired init_state
         self.add_nodes(pattern.input_nodes, input_state)
diff --git a/graphix/sim/tensornet.py b/graphix/sim/tensornet.py
index e246dee..802eb81 100644
--- a/graphix/sim/tensornet.py
+++ b/graphix/sim/tensornet.py
@@ -5,8 +5,10 @@ from copy import deepcopy
 
 import numpy as np
 import quimb.tensor as qtn
+import typing_extensions
 from quimb.tensor import Tensor, TensorNetwork
 
+import graphix.clifford
 from graphix import command
 from graphix.clifford import CLIFFORD, CLIFFORD_CONJ, CLIFFORD_MUL
 from graphix.ops import Ops
@@ -151,16 +153,32 @@ class TensorNetworkBackend:
             self.results[cmd.node] = result
             buffer = 2**0.5
 
-        # extract signals for adaptive angle
-        s_signal = np.sum(self.results[j] for j in cmd.s_domain)
-        t_signal = np.sum(self.results[j] for j in cmd.t_domain)
+        plane = Plane.XY
         angle = cmd.angle * np.pi
-        vop = cmd.vop
+        s_domain = cmd.s_domain
+        t_domain = cmd.t_domain
+        if cmd.plane == Plane.XY:
+            vop = 0
+        elif cmd.plane == Plane.YZ:
+            vop = graphix.clifford.H.index
+            angle = -angle
+            s_domain, t_domain = t_domain, s_domain
+        elif cmd.plane == Plane.XZ:
+            vop = (graphix.clifford.H @ graphix.clifford.S).index
+            angle = -angle
+            t_domain = t_domain + s_domain
+            s_domain, t_domain = t_domain, s_domain
+        else:
+            typing_extensions.assert_never(cmd.plane)
+
+        # extract signals for adaptive angle
+        s_signal = sum(self.results[j] for j in s_domain)
+        t_signal = sum(self.results[j] for j in t_domain)
         if int(s_signal % 2) == 1:
             vop = CLIFFORD_MUL[1, vop]
         if int(t_signal % 2) == 1:
             vop = CLIFFORD_MUL[3, vop]
-        proj_vec = proj_basis(angle, vop=vop, plane=cmd.plane, choice=result)
+        proj_vec = proj_basis(angle, vop=vop, plane=plane, choice=result)
 
         # buffer is necessary for maintaing the norm invariant
         proj_vec = proj_vec * buffer
diff --git a/graphix/transpiler.py b/graphix/transpiler.py
index 50c699b..58edd94 100644
--- a/graphix/transpiler.py
+++ b/graphix/transpiler.py
@@ -16,9 +16,12 @@ import graphix.pauli
 import graphix.sim.base_backend
 import graphix.sim.statevec
 from graphix import command, instruction
+from graphix.clifford import H
 from graphix.command import CommandKind, E, M, N, X, Z
 from graphix.ops import Ops
 from graphix.pattern import Pattern
+from graphix.pauli import Plane
+from graphix.sim.statevec import Statevec
 
 
 @dataclasses.dataclass
@@ -39,11 +42,11 @@ class SimulateResult:
     """
     The result of a simulation.
 
-    statevec : :class:`graphix.sim.statevec.Statevec` object
+    statevec : :class:`Statevec` object
     classical_measures : tuple[int,...], classical measures
     """
 
-    statevec: graphix.sim.statevec.Statevec
+    statevec: Statevec
     classical_measures: tuple[int, ...]
 
 
@@ -247,7 +250,7 @@ class Circuit:
         assert qubit in self.active_qubits
         self.instruction.append(instruction.I(target=qubit))
 
-    def m(self, qubit: int, plane: graphix.pauli.Plane, angle: float):
+    def m(self, qubit: int, plane: Plane, angle: float):
         """measure a quantum qubit
 
         The measured qubit cannot be used afterwards.
@@ -256,7 +259,7 @@ class Circuit:
         ---------
         qubit : int
             target qubit
-        plane : graphix.pauli.Plane
+        plane : Plane
         angle : float
         """
         assert qubit in self.active_qubits
@@ -816,7 +819,7 @@ class Circuit:
         if correction_instr.target == rx_instr.target:
             if correction_instr.kind == instruction.InstructionKind.ZC:
                 # add to the s-domain
-                self._M[rx_instr.meas_index].s_domain.extend(correction_instr.domain)
+                extend_domain(self._M[rx_instr.meas_index], correction_instr.domain)
                 self._commute_with_following(target)
             else:
                 self._commute_with_following(target)
@@ -833,7 +836,7 @@ class Circuit:
         assert ry_instr.kind == instruction.InstructionKind.RY
         if correction_instr.target == ry_instr.target:
             # add to the s-domain
-            self._M[ry_instr.meas_index].s_domain.extend(correction_instr.domain)
+            extend_domain(self._M[ry_instr.meas_index], correction_instr.domain)
             self._commute_with_following(target)
         else:
             self._commute_with_following(target)
@@ -849,7 +852,7 @@ class Circuit:
         if correction_instr.target == rz_instr.target:
             if correction_instr.kind == instruction.InstructionKind.XC:
                 # add to the s-domain
-                self._M[rz_instr.meas_index].s_domain.extend(correction_instr.domain)
+                extend_domain(self._M[rz_instr.meas_index], correction_instr.domain)
                 self._commute_with_following(target)
             else:
                 self._commute_with_following(target)
@@ -869,7 +872,7 @@ class Circuit:
             cond2 = correction_instr.target == rzz_instr.target
             if cond or cond2:
                 # add to the s-domain
-                self._M[rzz_instr.meas_index].s_domain.extend(correction_instr.domain)
+                extend_domain(self._M[rzz_instr.meas_index], correction_instr.domain)
         self._commute_with_following(target)
 
     def _commute_with_following(self, target: int):
@@ -988,14 +991,14 @@ class Circuit:
         return control_node, ancilla[1], seq
 
     @classmethod
-    def _m_command(self, input_node: int, plane: graphix.pauli.Plane, angle: float):
+    def _m_command(self, input_node: int, plane: Plane, angle: float):
         """MBQC commands for measuring qubit
 
         Parameters
         ---------
         input_node : int
             target node on graph
-        plane : graphix.pauli.Plane
+        plane : Plane
             plane of the measure
         angle : float
             angle of the measure (unit: pi radian)
@@ -1266,7 +1269,7 @@ class Circuit:
         """
         seq = [N(node=ancilla)]
         seq.append(E(nodes=(input_node, ancilla)))
-        seq.append(M(node=ancilla, angle=-angle / np.pi, vop=6))
+        seq.append(M(node=ancilla, angle=-angle / np.pi).clifford(H))
         seq.append(Z(node=input_node, domain=[ancilla]))
         return input_node, seq
 
@@ -1297,7 +1300,7 @@ class Circuit:
         seq = [N(node=ancilla)]
         seq.append(E(nodes=(control_node, ancilla)))
         seq.append(E(nodes=(target_node, ancilla)))
-        seq.append(M(node=ancilla, angle=-angle / np.pi, vop=6))
+        seq.append(M(node=ancilla, angle=-angle / np.pi).clifford(H))
         seq.append(Z(node=control_node, domain=[ancilla]))
         seq.append(Z(node=target_node, domain=[ancilla]))
         return control_node, target_node, seq
@@ -1537,13 +1540,13 @@ class Circuit:
         seq.append(E(nodes=(ancilla[8], ancilla[10])))
         seq.append(M(node=target_node))
         seq.append(M(node=control_node1))
-        seq.append(M(node=ancilla[0], angle=-1.75, s_domain=[target_node], vop=6))
+        seq.append(M(node=ancilla[0], angle=-1.75, s_domain=[target_node]).clifford(H))
         seq.append(M(node=ancilla[8], s_domain=[control_node1]))
-        seq.append(M(node=ancilla[2], angle=-0.25, s_domain=[target_node, ancilla[8]], vop=6))
+        seq.append(M(node=ancilla[2], angle=-0.25, s_domain=[target_node, ancilla[8]]).clifford(H))
         seq.append(M(node=control_node2, angle=-0.25))
-        seq.append(M(node=ancilla[3], angle=-1.75, s_domain=[ancilla[8], target_node], vop=6))
-        seq.append(M(node=ancilla[4], angle=-1.75, s_domain=[ancilla[8]], vop=6))
-        seq.append(M(node=ancilla[1], angle=-0.25, s_domain=[ancilla[8]], vop=6))
+        seq.append(M(node=ancilla[3], angle=-1.75, s_domain=[ancilla[8], target_node]).clifford(H))
+        seq.append(M(node=ancilla[4], angle=-1.75, s_domain=[ancilla[8]]).clifford(H))
+        seq.append(M(node=ancilla[1], angle=-0.25, s_domain=[ancilla[8]]).clifford(H))
         seq.append(
             M(
                 node=ancilla[5],
@@ -1616,9 +1619,9 @@ class Circuit:
         """
 
         if input_state is None:
-            state = graphix.sim.statevec.Statevec(nqubit=self.width)
+            state = Statevec(nqubit=self.width)
         else:
-            state = graphix.sim.statevec.Statevec(nqubit=self.width, data=input_state)
+            state = Statevec(nqubit=self.width, data=input_state)
 
         classical_measures = []
 
@@ -1660,3 +1663,10 @@ class Circuit:
                 raise ValueError(f"Unknown instruction: {instr}")
 
         return SimulateResult(state, classical_measures)
+
+
+def extend_domain(measure: M, domain: list[int]) -> None:
+    if measure.plane == Plane.XY:
+        measure.s_domain.extend(domain)
+    else:
+        measure.t_domain.extend(domain)


Test patch:
diff --git a/tests/test_pattern.py b/tests/test_pattern.py
index 14580cf..9205558 100644
--- a/tests/test_pattern.py
+++ b/tests/test_pattern.py
@@ -1,18 +1,23 @@
 from __future__ import annotations
 
+import collections.abc
+import itertools
 import sys
-from typing import TYPE_CHECKING, Literal
+import typing
+from typing import TYPE_CHECKING
 
 import numpy as np
 import pytest
 
+import graphix.clifford
 import graphix.ops
 import graphix.pauli
 import graphix.sim.base_backend
 import graphix.states
 import tests.random_circuit as rc
-from graphix.command import M, N
+from graphix.command import E, M, N
 from graphix.pattern import CommandNode, Pattern
+from graphix.pauli import Plane
 from graphix.sim.density_matrix import DensityMatrix
 from graphix.sim.statevec import Statevec
 from graphix.simulator import PatternSimulator
@@ -33,6 +38,17 @@ def compare_backend_result_with_statevec(backend: str, backend_state, statevec:
         raise NotImplementedError(backend)
 
 
+Outcome = typing.Literal[0, 1]
+
+
+class IterGenerator:
+    def __init__(self, it: collections.abc.Iterable[Outcome]) -> None:
+        self.__it = iter(it)
+
+    def choice(self, _outcomes: list[Outcome]) -> Outcome:
+        return next(self.__it)
+
+
 class TestPattern:
     # this fails without behaviour modification
     def test_manual_generation(self) -> None:
@@ -81,7 +97,7 @@ class TestPattern:
 
     @pytest.mark.filterwarnings("ignore:Simulating using densitymatrix backend with no noise.")
     @pytest.mark.parametrize("backend", ["statevector", "densitymatrix", "tensornetwork"])
-    def test_empty_output_nodes(self, backend: Literal["statevector", "densitymatrix", "tensornetwork"]) -> None:
+    def test_empty_output_nodes(self, backend: typing.Literal["statevector", "densitymatrix", "tensornetwork"]) -> None:
         pattern = Pattern(input_nodes=[0])
         pattern.add(M(node=0, angle=0.5))
 
@@ -311,7 +327,7 @@ class TestPattern:
         vop_list = [0, 5, 6]  # [identity, S gate, H gate]
         pattern = Pattern(input_nodes=list(range(len(preset_meas_plane))))
         for i in range(len(preset_meas_plane)):
-            pattern.add(M(node=i, plane=preset_meas_plane[i], vop=vop_list[i % 3]))
+            pattern.add(M(node=i, plane=preset_meas_plane[i]).clifford(graphix.clifford.get(vop_list[i % 3])))
         ref_meas_plane = {
             0: graphix.pauli.Plane.XY,
             1: graphix.pauli.Plane.XY,
@@ -326,6 +342,29 @@ class TestPattern:
         meas_plane = pattern.get_meas_plane()
         assert meas_plane == ref_meas_plane
 
+    @pytest.mark.parametrize("plane", Plane)
+    @pytest.mark.parametrize("method", ["local", "global"])
+    def test_shift_signals_plane(self, plane: Plane, method: str) -> None:
+        pattern = Pattern(input_nodes=[0])
+        for i in (1, 2, 3):
+            pattern.add(N(node=i))
+            pattern.add(E(nodes=[0, i]))
+        pattern.add(M(node=0, angle=0.5))
+        pattern.add(M(node=1, angle=0.5))
+        pattern.add(M(node=2, angle=0.5, plane=plane, s_domain=[0], t_domain=[1]))
+        pattern_ref = pattern.copy()
+        pattern.standardize(method="global")
+        signal_dict = pattern.shift_signals(method=method)
+        # Test for every possible outcome of each measure
+        for outcomes_ref in itertools.product(*([[0, 1]] * 3)):
+            state_ref = pattern_ref.simulate_pattern(pr_calc=False, rng=IterGenerator(iter(outcomes_ref)))
+            outcomes_p = [
+                1 - outcome if sum(outcomes_ref[i] for i in swapped) % 2 == 1 else outcome
+                for outcome, swapped in zip(outcomes_ref, signal_dict.values())
+            ]
+            state_p = pattern.simulate_pattern(pr_calc=False, rng=IterGenerator(iter(outcomes_p)))
+            assert np.abs(np.dot(state_p.flatten().conjugate(), state_ref.flatten())) == pytest.approx(1)
+
 
 def cp(circuit: Circuit, theta: float, control: int, target: int) -> None:
     """Controlled rotation gate, decomposed"""