diff --git a/CHANGELOG.md b/CHANGELOG.md
index aa4b70942..a588bae2f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -15,6 +15,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   - Introduced new wrapper methods in the `Pattern` class: `graphix.pattern.Pattern.extract_partial_order_layers`, `graphix.pattern.Pattern.extract_causal_flow` and `graphix.pattern.Pattern.extract_gflow`.
   - Introduced new module `graphix.flow._partial_order` with the function :func:`compute_topological_generations`.
 
+- #392: Added `graphix.pattern.Pattern.remove_input_nodes` method which removes the input nodes from the pattern and replaces them with N commands.
+
 - #385
   - Introduced `graphix.flow.core.XZCorrections.check_well_formed` which verifies the correctness of an XZ-corrections instance and raises an exception if incorrect.
   - Added XZ-correction exceptions to module `graphix.flow.core.exceptions`.
@@ -35,6 +37,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Fixed
 
+- #392: `Pattern.remove_input_nodes` is required before the `Pattern.perform_pauli_measurements` method to ensure input nodes are removed and fixed in the |+> state.
+
 - #379: Removed unnecessary `meas_index` from API for rotation instructions `RZ`, `RY` and `RX`.
 
 - #347: Adapted existing method `graphix.opengraph.OpenGraph.isclose` to the new API introduced in #358.
diff --git a/docs/source/tutorial.rst b/docs/source/tutorial.rst
index 3d7ba2093..7364e4be6 100644
--- a/docs/source/tutorial.rst
+++ b/docs/source/tutorial.rst
@@ -175,12 +175,13 @@ It is known that quantum circuit consisting of Pauli basis states, Clifford gate
 <https://en.wikipedia.org/wiki/Gottesman%E2%80%93Knill_theorem>`_; e.g. the graph state simulator runs in :math:`\mathcal{O}(n \log n)` time).
 The Pauli measurement part of the MBQC is exactly this, and they can be preprocessed by our graph state simulator :class:`~graphix.graphsim.GraphState` - see :doc:`lc-mbqc` for more detailed description.
 
-We can call this in a line by calling :meth:`~graphix.pattern.Pattern.perform_pauli_measurements()` of :class:`~graphix.pattern.Pattern` object, which acts as the optimization routine of the measurement pattern.
+We can call this in a line by calling :meth:`~graphix.pattern.Pattern.remove_input_nodes` followed by :meth:`~graphix.pattern.Pattern.perform_pauli_measurements()` (both methods of the :class:`~graphix.pattern.Pattern` object). The first method removes the input nodes, while the second method optimizes the measurement pattern.
 We get an updated measurement pattern without Pauli measurements as follows:
 
+>>> pattern.remove_input_nodes()
 >>> pattern.perform_pauli_measurements()
 >>> pattern
-Pattern(input_nodes=[0, 1], cmds=[N(3), N(7), E((0, 3)), E((1, 3)), E((1, 7)), M(0, Plane.YZ, 0.2907266109187514), M(1, Plane.YZ, 0.01258854060311348), C(3, Clifford.I), C(7, Clifford.I), Z(3, {0, 1, 5}), Z(7, {1, 5}), X(3, {2}), X(7, {2, 4, 6})], output_nodes=[3, 7])
+Pattern(input_nodes=[], cmds=[N(0), N(1), N(3), N(7), E((0, 3)), E((1, 3)), E((1, 7)), M(0, Plane.YZ, 0.2907266109187514), M(1, Plane.YZ, 0.01258854060311348), C(3, Clifford.I), C(7, Clifford.I), Z(3, {0, 1, 5}), Z(7, {1, 5}), X(3, {2}), X(7, {2, 4, 6})], output_nodes=[3, 7])
 
 
 Notice that all measurements with angle=0 (Pauli X measurements) disappeared - this means that a part of quantum computation was *classically* (and efficiently) preprocessed such that we only need much smaller quantum resource.
@@ -210,7 +211,7 @@ We can simply call :meth:`~graphix.pattern.Pattern.minimize_space()` to reduce t
 
 >>> pattern.minimize_space()
 >>> pattern
-Pattern(input_nodes=[0, 1], cmds=[N(3), E((0, 3)), M(0, Plane.YZ, 0.11120090987081546), E((1, 3)), N(7), E((1, 7)), M(1, Plane.YZ, 0.230565199664617), C(3, Clifford.I), C(7, Clifford.I), Z(3, {0, 1, 5}), Z(7, {1, 5}), X(3, {2}), X(7, {2, 4, 6})], output_nodes=[3, 7])
+Pattern(input_nodes=[], cmds=[N(0), N(1), N(3), E((0, 3)), M(0, Plane.YZ, 0.11120090987081546), E((1, 3)), N(7), E((1, 7)), M(1, Plane.YZ, 0.230565199664617), C(3, Clifford.I), C(7, Clifford.I), Z(3, {0, 1, 5}), Z(7, {1, 5}), X(3, {2}), X(7, {2, 4, 6})], output_nodes=[3, 7])
 
 
 With the original measurement pattern, the simulation should have proceeded as follows, with maximum of four qubits on the memory.
diff --git a/examples/deutsch_jozsa.py b/examples/deutsch_jozsa.py
index 771c8ca50..4db0f055a 100644
--- a/examples/deutsch_jozsa.py
+++ b/examples/deutsch_jozsa.py
@@ -71,8 +71,9 @@
 print(pattern.to_ascii(left_to_right=True, limit=15))
 
 # %%
-# Now we preprocess all Pauli measurements
+# Now we preprocess all Pauli measurements, which requires that we move inputs to N commands
 
+pattern.remove_input_nodes()
 pattern.perform_pauli_measurements()
 print(
     pattern.to_ascii(
diff --git a/examples/mbqc_vqe.py b/examples/mbqc_vqe.py
index 269fb1ce6..ae75a7f7f 100644
--- a/examples/mbqc_vqe.py
+++ b/examples/mbqc_vqe.py
@@ -89,6 +89,7 @@ def build_mbqc_pattern(self, params: Iterable[Angle]) -> Pattern:
         pattern = circuit.transpile().pattern
         pattern.standardize()
         pattern.shift_signals()
+        pattern.remove_input_nodes()
         pattern.perform_pauli_measurements()  # Perform Pauli measurements
         return pattern
 
diff --git a/examples/qaoa.py b/examples/qaoa.py
index 51355b05c..57b2cafc9 100644
--- a/examples/qaoa.py
+++ b/examples/qaoa.py
@@ -40,6 +40,7 @@
 # %%
 # perform Pauli measurements and plot the new (minimal) graph to perform the same quantum computation
 
+pattern.remove_input_nodes()
 pattern.perform_pauli_measurements()
 pattern.draw_graph(flow_from_pattern=False)
 
diff --git a/examples/qft_with_tn.py b/examples/qft_with_tn.py
index 198772513..578bf9945 100644
--- a/examples/qft_with_tn.py
+++ b/examples/qft_with_tn.py
@@ -62,6 +62,7 @@ def qft(circuit, n):
 # %%
 # Using efficient graph state simulator `graphix.graphsim`, we can classically preprocess Pauli measurements.
 # We are currently improving the speed of this process by using rust-based graph manipulation backend.
+pattern.remove_input_nodes()
 pattern.perform_pauli_measurements()
 
 
diff --git a/examples/tn_simulation.py b/examples/tn_simulation.py
index e92ff979a..1a18d9156 100644
--- a/examples/tn_simulation.py
+++ b/examples/tn_simulation.py
@@ -70,7 +70,7 @@ def ansatz(circuit, n, gamma, beta, iterations):
 
 # %%
 # Optimizing by performing Pauli measurements in the pattern using efficient stabilizer simulator.
-
+pattern.remove_input_nodes()
 pattern.perform_pauli_measurements()
 
 # %%
@@ -181,6 +181,7 @@ def cost(params, n, ham, quantum_iter, slice_index, opt=None):
     pattern = circuit.transpile().pattern
     pattern.standardize()
     pattern.shift_signals()
+    pattern.remove_input_nodes()
     pattern.perform_pauli_measurements()
     mbqc_tn = pattern.simulate_pattern(backend="tensornetwork", graph_prep="parallel")
     exp_val = 0
diff --git a/examples/visualization.py b/examples/visualization.py
index d837c8f40..983b9657b 100644
--- a/examples/visualization.py
+++ b/examples/visualization.py
@@ -37,8 +37,8 @@
 
 # %%
 # next, show the gflow:
-
-pattern.perform_pauli_measurements(leave_input=True)
+pattern.remove_input_nodes()
+pattern.perform_pauli_measurements()
 pattern.draw_graph(flow_from_pattern=False, show_measurement_planes=True, node_distance=(1, 0.6))
 
 
diff --git a/graphix/pattern.py b/graphix/pattern.py
index bab294b7b..976f26fd9 100644
--- a/graphix/pattern.py
+++ b/graphix/pattern.py
@@ -124,6 +124,7 @@ def add(self, cmd: Command) -> None:
         cmd : :class:`graphix.command.Command`
             MBQC command.
         """
+        self._pauli_preprocessed = False
         if cmd.kind == CommandKind.N:
             self.__n_node += 1
             self.__output_nodes.append(cmd.node)
@@ -137,6 +138,7 @@ def extend(self, *cmds: Command | Iterable[Command]) -> None:
 
         :param cmds: sequences of commands
         """
+        self._pauli_preprocessed = False
         for item in cmds:
             if isinstance(item, Iterable):
                 for cmd in item:
@@ -146,6 +148,7 @@ def extend(self, *cmds: Command | Iterable[Command]) -> None:
 
     def clear(self) -> None:
         """Clear the sequence of pattern commands."""
+        self._pauli_preprocessed = False
         self.__n_node = len(self.__input_nodes)
         self.__seq = []
         self.__output_nodes = list(self.__input_nodes)
@@ -157,6 +160,7 @@ def replace(self, cmds: list[Command], input_nodes: list[int] | None = None) ->
 
         :param input_nodes: optional, list of input qubits (by default, keep the same input nodes as before)
         """
+        self._pauli_preprocessed = False
         if input_nodes is not None:
             self.__input_nodes = list(input_nodes)
         self.clear()
@@ -202,6 +206,7 @@ def compose(
         - Input (and, respectively, output) nodes in the returned pattern have the order of the pattern `self` followed by those of the pattern `other`. Merged nodes are removed.
         - If `preserve_mapping = True` and :math:`|M_1| = |I_2| = |O_2|`, then the outputs of the returned pattern are the outputs of pattern `self`, where the nth merged output is replaced by the output of pattern `other` corresponding to its nth input instead.
         """
+        self._pauli_preprocessed = False
         nodes_p1 = self.extract_nodes() | self.results.keys()  # Results contain preprocessed Pauli nodes
         nodes_p2 = other.extract_nodes() | other.results.keys()
 
@@ -1402,14 +1407,21 @@ def simulate_pattern(
         sim.run(input_state, rng=rng)
         return sim.backend.state
 
-    def perform_pauli_measurements(self, leave_input: bool = False, ignore_pauli_with_deps: bool = False) -> None:
+    def remove_input_nodes(self) -> None:
+        """Remove the input nodes from the pattern and replace them with N commands.
+
+        This removes the possibility of choosing the input state, fixing the input state to the plus state.
+        .. seealso:: :class:`graphix.command.N`
+        """
+        self.__seq[0:0] = [command.N(node=node) for node in self.input_nodes]
+        empty_nodes: list[int] = []
+        self.__input_nodes = empty_nodes
+
+    def perform_pauli_measurements(self, ignore_pauli_with_deps: bool = False) -> None:
         """Perform Pauli measurements in the pattern using efficient stabilizer simulator.
 
         Parameters
         ----------
-        leave_input : bool
-            Optional (*False* by default).
-            If *True*, measurements on input nodes are preserved as-is in the pattern.
         ignore_pauli_with_deps : bool
             Optional (*False* by default).
             If *True*, Pauli measurements with domains depending on other measures are preserved as-is in the pattern.
@@ -1418,9 +1430,9 @@ def perform_pauli_measurements(self, leave_input: bool = False, ignore_pauli_wit
         .. seealso:: :func:`measure_pauli`
 
         """
-        # if not ignore_pauli_with_deps:
-        #    self.move_pauli_measurements_to_the_front()
-        measure_pauli(self, leave_input, copy=False, ignore_pauli_with_deps=ignore_pauli_with_deps)
+        if self.input_nodes:
+            raise ValueError("Remove inputs with `self.remove_input_nodes()` before performing Pauli presimulation.")
+        measure_pauli(self, copy=False, ignore_pauli_with_deps=ignore_pauli_with_deps)
 
     def draw_graph(
         self,
@@ -1652,23 +1664,17 @@ def __str__(self) -> str:
         assert_never(self.reason)
 
 
-def measure_pauli(
-    pattern: Pattern, leave_input: bool, *, copy: bool = False, ignore_pauli_with_deps: bool = False
-) -> Pattern:
+def measure_pauli(pattern: Pattern, *, copy: bool = False, ignore_pauli_with_deps: bool = False) -> Pattern:
     """Perform Pauli measurement of a pattern by fast graph state simulator.
 
-    Uses the decorated-graph method implemented in graphix.graphsim to perform
-    the measurements in Pauli bases, and then sort remaining nodes back into
-    pattern together with Clifford commands.
+    Uses the decorated-graph method implemented in graphix.graphsim to perform the measurements in Pauli bases, and then sort remaining nodes back into
+    pattern together with Clifford commands. Users are required to ensure there are no input nodes with :func:`graphix.pattern.Pattern.remove_input_nodes` before using this function.
 
     TODO: non-XY plane measurements in original pattern
 
     Parameters
     ----------
     pattern : graphix.pattern.Pattern object
-    leave_input : bool
-        True: input nodes will not be removed
-        False: all the nodes measured in Pauli bases will be removed
     copy : bool
         True: changes will be applied to new copied object and will be returned
         False: changes will be applied to the supplied Pattern object
@@ -1684,20 +1690,20 @@ def measure_pauli(
         only returned if copy argument is True.
 
 
+    .. seealso:: :class:`graphix.pattern.Pattern.remove_input_nodes`
     .. seealso:: :class:`graphix.graphsim.GraphState`
     """
+    pat = Pattern() if copy else pattern
     standardized_pattern = optimization.StandardizedPattern.from_pattern(pattern)
     if not ignore_pauli_with_deps:
         standardized_pattern = standardized_pattern.perform_pauli_pushing()
     output_nodes = set(pattern.output_nodes)
     graph = standardized_pattern.extract_graph()
     graph_state = GraphState(nodes=graph.nodes, edges=graph.edges, vops=standardized_pattern.c_dict)
-    results: dict[int, Outcome] = {}
-    to_measure, non_pauli_meas = pauli_nodes(standardized_pattern, leave_input)
-    if not leave_input and len(list(set(pattern.input_nodes) & {i[0].node for i in to_measure})) > 0:
-        new_inputs = []
-    else:
-        new_inputs = pattern.input_nodes
+    results: dict[int, Outcome] = pat.results
+    to_measure, non_pauli_meas = pauli_nodes(standardized_pattern)
+    if not to_measure:
+        return pattern
     for cmd in to_measure:
         pattern_cmd = cmd[0]
         measurement_basis = cmd[1]
@@ -1746,7 +1752,7 @@ def measure_pauli(
     # update command sequence
     vops = graph_state.get_vops()
     new_seq: list[Command] = []
-    new_seq.extend(command.N(node=index) for index in set(graph_state.nodes) - set(new_inputs))
+    new_seq.extend(command.N(node=index) for index in set(graph_state.nodes))
     new_seq.extend(command.E(nodes=edge) for edge in graph_state.edges)
     new_seq.extend(
         cmd.clifford(Clifford(vops[cmd.node])) for cmd in standardized_pattern.m_list if cmd.node in graph_state.nodes
@@ -1758,10 +1764,7 @@ def measure_pauli(
     )
     new_seq.extend(command.Z(node=node, domain=set(domain)) for node, domain in standardized_pattern.z_dict.items())
     new_seq.extend(command.X(node=node, domain=set(domain)) for node, domain in standardized_pattern.x_dict.items())
-
-    pat = Pattern() if copy else pattern
-
-    pat.replace(new_seq, input_nodes=new_inputs)
+    pat.replace(new_seq, input_nodes=[])
     pat.reorder_output_nodes(standardized_pattern.output_nodes)
     assert pat.n_node == len(graph_state.nodes)
     pat.results = results
@@ -1769,15 +1772,12 @@ def measure_pauli(
     return pat
 
 
-def pauli_nodes(
-    pattern: optimization.StandardizedPattern, leave_input: bool
-) -> tuple[list[tuple[command.M, PauliMeasurement]], set[int]]:
+def pauli_nodes(pattern: optimization.StandardizedPattern) -> tuple[list[tuple[command.M, PauliMeasurement]], set[int]]:
     """Return the list of measurement commands that are in Pauli bases and that are not dependent on any non-Pauli measurements.
 
     Parameters
     ----------
     pattern : optimization.StandardizedPattern
-    leave_input : bool
 
     Returns
     -------
@@ -1790,7 +1790,7 @@ def pauli_nodes(
     non_pauli_node: set[int] = set()
     for cmd in pattern.m_list:
         pm = PauliMeasurement.try_from(cmd.plane, cmd.angle)  # None returned if the measurement is not in Pauli basis
-        if pm is not None and (cmd.node not in pattern.input_nodes or not leave_input):
+        if pm is not None:
             # Pauli measurement to be removed
             if pm.axis == Axis.X:
                 if cmd.t_domain & non_pauli_node:  # cmd depend on non-Pauli measurement
diff --git a/tests/baseline/test_draw_graph_reference_False.png b/tests/baseline/test_draw_graph_reference_False.png
index 22999333f..bb7b73fc4 100644
Binary files a/tests/baseline/test_draw_graph_reference_False.png and b/tests/baseline/test_draw_graph_reference_False.png differ
diff --git a/tests/baseline/test_draw_graph_reference_True.png b/tests/baseline/test_draw_graph_reference_True.png
index 2f79065a6..7252879c3 100644
Binary files a/tests/baseline/test_draw_graph_reference_True.png and b/tests/baseline/test_draw_graph_reference_True.png differ
diff --git a/tests/test_gflow.py b/tests/test_gflow.py
index 87fbae81a..0d938be87 100644
--- a/tests/test_gflow.py
+++ b/tests/test_gflow.py
@@ -538,6 +538,7 @@ def test_rand_circ_gflow(self, fx_rng: Generator) -> None:
         pattern = circ.transpile().pattern
         pattern.standardize()
         pattern.shift_signals()
+        pattern.remove_input_nodes()
         pattern.perform_pauli_measurements()
         graph = pattern.extract_graph()
         input_ = set()
diff --git a/tests/test_optimization.py b/tests/test_optimization.py
index 4f3af777a..d1375fdda 100644
--- a/tests/test_optimization.py
+++ b/tests/test_optimization.py
@@ -61,6 +61,7 @@ def test_incorporate_pauli_results(fx_bg: PCG64, jumps: int) -> None:
     pattern = circuit.transpile().pattern
     pattern.standardize()
     pattern.shift_signals()
+    pattern.remove_input_nodes()
     pattern.perform_pauli_measurements()
     pattern2 = incorporate_pauli_results(pattern)
     state = pattern.simulate_pattern(rng=rng)
@@ -78,6 +79,7 @@ def test_flow_after_pauli_preprocessing(fx_bg: PCG64, jumps: int) -> None:
     pattern.standardize()
     pattern.shift_signals()
     # pattern.move_pauli_measurements_to_the_front()
+    pattern.remove_input_nodes()
     pattern.perform_pauli_measurements()
     pattern2 = incorporate_pauli_results(pattern)
     gflow = pattern2.extract_gflow()
@@ -93,6 +95,7 @@ def test_remove_useless_domains(fx_bg: PCG64, jumps: int) -> None:
     pattern = circuit.transpile().pattern
     pattern.standardize()
     pattern.shift_signals()
+    pattern.remove_input_nodes()
     pattern.perform_pauli_measurements()
     pattern2 = remove_useless_domains(pattern)
     state = pattern.simulate_pattern(rng=rng)
diff --git a/tests/test_parameter.py b/tests/test_parameter.py
index b57f9bba4..8b701f2ce 100644
--- a/tests/test_parameter.py
+++ b/tests/test_parameter.py
@@ -158,6 +158,7 @@ def test_random_circuit_with_parameters(fx_bg: PCG64, jumps: int, use_xreplace:
     pattern = circuit.transpile().pattern
     pattern.standardize()
     pattern.shift_signals()
+    pattern.remove_input_nodes()
     pattern.perform_pauli_measurements()
     pattern.minimize_space()
     assignment: dict[Parameter, float] = {alpha: rng.uniform(high=2), beta: rng.uniform(high=2)}
diff --git a/tests/test_pattern.py b/tests/test_pattern.py
index a1e600d3b..0b282cac9 100644
--- a/tests/test_pattern.py
+++ b/tests/test_pattern.py
@@ -97,6 +97,7 @@ def test_pauli_non_contiguous(self) -> None:
                 M(node=0, plane=Plane.XY, angle=0.0, s_domain=set(), t_domain=set()),
             ]
         )
+        pattern.remove_input_nodes()
         pattern.perform_pauli_measurements()
 
     @pytest.mark.parametrize("jumps", range(1, 11))
@@ -109,6 +110,7 @@ def test_minimize_space_with_gflow(self, fx_bg: PCG64, jumps: int) -> None:
         pattern = circuit.transpile().pattern
         pattern.standardize()
         pattern.shift_signals(method="mc")
+        pattern.remove_input_nodes()
         pattern.perform_pauli_measurements()
         pattern.minimize_space()
         state = circuit.simulate_statevector().statevec
@@ -186,6 +188,7 @@ def test_pauli_measurement_random_circuit(self, fx_bg: PCG64, jumps: int, backen
         pattern = circuit.transpile().pattern
         pattern.standardize()
         pattern.shift_signals(method="mc")
+        pattern.remove_input_nodes()
         pattern.perform_pauli_measurements()
         pattern.minimize_space()
         state = circuit.simulate_statevector().statevec
@@ -204,6 +207,7 @@ def test_pauli_measurement_random_circuit_all_paulis(
         pattern = circuit.transpile().pattern
         pattern.standardize()
         pattern.shift_signals(method="mc")
+        pattern.remove_input_nodes()
         pattern.perform_pauli_measurements(ignore_pauli_with_deps=ignore_pauli_with_deps)
         assert ignore_pauli_with_deps or not any(
             PauliMeasurement.try_from(cmd.plane, cmd.angle) for cmd in pattern if cmd.kind == CommandKind.M
@@ -216,26 +220,12 @@ def test_pauli_measurement_single(self, plane: Plane, angle: float) -> None:
         pattern.add(E(nodes=(0, 1)))
         pattern.add(M(node=0, plane=plane, angle=angle))
         pattern_ref = pattern.copy()
+        pattern.remove_input_nodes()
         pattern.perform_pauli_measurements()
         state = pattern.simulate_pattern()
         state_ref = pattern_ref.simulate_pattern(branch_selector=ConstBranchSelector(0))
         assert np.abs(np.dot(state.flatten().conjugate(), state_ref.flatten())) == pytest.approx(1)
 
-    @pytest.mark.parametrize("jumps", range(1, 11))
-    def test_pauli_measurement_leave_input_random_circuit(self, fx_bg: PCG64, jumps: int) -> None:
-        rng = Generator(fx_bg.jumped(jumps))
-        nqubits = 3
-        depth = 3
-        circuit = rand_circuit(nqubits, depth, rng)
-        pattern = circuit.transpile().pattern
-        pattern.standardize()
-        pattern.shift_signals(method="mc")
-        pattern.perform_pauli_measurements(leave_input=True)
-        pattern.minimize_space()
-        state = circuit.simulate_statevector().statevec
-        state_mbqc = pattern.simulate_pattern(rng=rng)
-        assert np.abs(np.dot(state_mbqc.flatten().conjugate(), state.flatten())) == pytest.approx(1)
-
     def test_pauli_measurement(self) -> None:
         # test pattern is obtained from 3-qubit QFT with pauli measurement
         circuit = Circuit(3)
@@ -243,7 +233,6 @@ def test_pauli_measurement(self) -> None:
             circuit.h(i)
         circuit.x(1)
         circuit.x(2)
-
         # QFT
         circuit.h(2)
         cp(circuit, ANGLE_PI / 4, 0, 2)
@@ -251,19 +240,26 @@ def test_pauli_measurement(self) -> None:
         circuit.h(1)
         cp(circuit, ANGLE_PI / 2, 0, 1)
         circuit.h(0)
-        swap(circuit, 0, 2)
-
+        circuit.swap(0, 2)
         pattern = circuit.transpile().pattern
         pattern.standardize()
         pattern.shift_signals(method="mc")
+        pattern.remove_input_nodes()
         pattern.perform_pauli_measurements()
-
         isolated_nodes = pattern.extract_isolated_nodes()
-        # 48-node is the isolated and output node.
-        isolated_nodes_ref = {48}
-
+        # 42-node is the isolated and output node.
+        isolated_nodes_ref = {42}
         assert isolated_nodes == isolated_nodes_ref
 
+    def test_pauli_measurement_error(self, fx_rng: Generator) -> None:
+        nqubits = 2
+        depth = 1
+        circuit = rand_circuit(nqubits, depth, fx_rng)
+        pattern = circuit.transpile().pattern
+        pattern.standardize()
+        with pytest.raises(ValueError):
+            pattern.perform_pauli_measurements()
+
     def test_pauli_measurement_leave_input(self) -> None:
         # test pattern is obtained from 3-qubit QFT with pauli measurement
         circuit = Circuit(3)
@@ -280,17 +276,62 @@ def test_pauli_measurement_leave_input(self) -> None:
         cp(circuit, ANGLE_PI / 2, 0, 1)
         circuit.h(0)
         swap(circuit, 0, 2)
+        pattern = circuit.transpile().pattern
+        pattern.standardize()
+        with pytest.raises(ValueError):
+            pattern.perform_pauli_measurements()
 
+    @pytest.mark.parametrize("jumps", range(1, 6))
+    @pytest.mark.parametrize("ignore_pauli_with_deps", [False, True])
+    def test_pauli_measured_against_nonmeasured(self, fx_bg: PCG64, jumps: int, ignore_pauli_with_deps: bool) -> None:
+        rng = Generator(fx_bg.jumped(jumps))
+        nqubits = 2
+        depth = 2
+        circuit = rand_circuit(nqubits, depth, rng)
         pattern = circuit.transpile().pattern
         pattern.standardize()
-        pattern.shift_signals(method="mc")
-        pattern.perform_pauli_measurements(leave_input=True)
+        pattern1 = copy.deepcopy(pattern)
+        pattern1.remove_input_nodes()
+        pattern1.perform_pauli_measurements(ignore_pauli_with_deps=ignore_pauli_with_deps)
+        state = pattern.simulate_pattern(rng=rng)
+        state1 = pattern1.simulate_pattern(rng=rng)
+        assert np.abs(np.dot(state.flatten().conjugate(), state1.flatten())) == pytest.approx(1)
 
-        isolated_nodes = pattern.extract_isolated_nodes()
-        # There is no isolated node.
-        isolated_nodes_ref: set[int] = set()
+    @pytest.mark.parametrize("jumps", range(1, 4))
+    def test_pauli_repeated_measurement(self, fx_bg: PCG64, jumps: int) -> None:
+        rng = Generator(fx_bg.jumped(jumps))
+        nqubits = 2
+        depth = 2
+        circuit = rand_circuit(nqubits, depth, rng, use_ccx=False)
+        pattern = circuit.transpile().pattern
+        pattern.remove_input_nodes()
+        assert not pattern.results
+        pattern.perform_pauli_measurements()
+        assert pattern.results
+        pattern.perform_pauli_measurements()
+        assert pattern.results
 
-        assert isolated_nodes == isolated_nodes_ref
+    @pytest.mark.parametrize("jumps", range(1, 4))
+    def test_pauli_repeated_measurement_compose(self, fx_bg: PCG64, jumps: int) -> None:
+        rng = Generator(fx_bg.jumped(jumps))
+        nqubits = 2
+        depth = 2
+        circuit = rand_circuit(nqubits, depth, rng, use_ccx=False)
+        circuit1 = rand_circuit(nqubits, depth, rng, use_ccx=False)
+        pattern = circuit.transpile().pattern
+        pattern1 = circuit1.transpile().pattern
+        composed_pattern, _ = pattern.compose(
+            pattern1, mapping=dict(zip(pattern1.input_nodes, pattern.output_nodes, strict=True)), preserve_mapping=True
+        )
+        pattern.remove_input_nodes()
+        pattern1.remove_input_nodes()
+        assert not pattern.results
+        assert not pattern1.results
+        pattern.perform_pauli_measurements()
+        pattern1.perform_pauli_measurements()
+        composed_pattern.remove_input_nodes()
+        composed_pattern.perform_pauli_measurements()
+        assert abs(len(composed_pattern.results) - len(pattern.results) - len(pattern1.results)) <= 2
 
     def test_get_meas_plane(self) -> None:
         preset_meas_plane = [
@@ -382,6 +423,7 @@ def test_pauli_measurement_then_standardize(self, fx_bg: PCG64, jumps: int) -> N
         depth = 3
         circuit = rand_circuit(nqubits, depth, rng)
         pattern = circuit.transpile().pattern
+        pattern.remove_input_nodes()
         pattern.perform_pauli_measurements()
         pattern.standardize()
         pattern.minimize_space()
@@ -594,41 +636,6 @@ def test_compose_3(self) -> None:
         assert p_1 == pc_1
         assert p_2 == pc_2
 
-    #  Pattern composition with Pauli preprocessing
-    def test_compose_4(self, fx_rng: Generator) -> None:
-        alpha = fx_rng.random()
-        i1 = [0]
-        o1 = [2]
-        cmds1: list[Command] = [N(1), N(2), E((0, 1)), E((1, 2)), M(0, angle=-alpha), M(1), X(2, {1}), Z(2, {0})]
-        p1 = Pattern(cmds=cmds1, input_nodes=i1, output_nodes=o1)
-        p2 = Pattern(cmds=cmds1, input_nodes=i1, output_nodes=o1)
-
-        p1.perform_pauli_measurements()
-        p2.perform_pauli_measurements()
-
-        mapping = {0: 2, 1: 3, 2: 4}
-        pc, mapping_complete = p1.compose(p2, mapping=mapping)
-
-        i = [0]
-        o = [4]
-        cmds: list[Command] = [
-            N(2),
-            E((0, 2)),
-            M(0, plane=Plane.YZ, angle=alpha),
-            Z(2, {0}),
-            X(2, {1}),
-            N(4),
-            E((2, 4)),
-            M(2, plane=Plane.YZ, angle=alpha),
-            Z(4, {2}),
-            X(4, {3}),
-        ]
-        p = Pattern(cmds=cmds, input_nodes=i, output_nodes=o)
-        p.results = {1: 0, 3: 0}
-
-        assert p == pc
-        assert mapping_complete == mapping
-
     # Equivalence between pattern and circuit composition
     def test_compose_5(self, fx_rng: Generator) -> None:
         circuit_1 = Circuit(1)
@@ -679,6 +686,7 @@ def test_compose_7(self, fx_rng: Generator) -> None:
         circuit_1.h(0)
         circuit_1.rz(0, alpha)
         p1 = circuit_1.transpile().pattern
+        p1.remove_input_nodes()
         p1.perform_pauli_measurements()
 
         circuit_2 = Circuit(1)
@@ -805,6 +813,7 @@ def test_extract_partial_order_layers_results(self) -> None:
         c = Circuit(1)
         c.rz(0, 0.2)
         p = c.transpile().pattern
+        p.remove_input_nodes()
         p.perform_pauli_measurements()
         assert p.extract_partial_order_layers() == (frozenset({2}), frozenset({0}))
 
@@ -930,6 +939,8 @@ def test_extract_causal_flow_rnd_circuit(self, fx_bg: PCG64, jumps: int) -> None
         p_ref = circuit_1.transpile().pattern
         p_test = p_ref.extract_causal_flow().to_corrections().to_pattern()
 
+        p_ref.remove_input_nodes()
+        p_test.remove_input_nodes()
         p_ref.perform_pauli_measurements()
         p_test.perform_pauli_measurements()
 
@@ -946,7 +957,8 @@ def test_extract_gflow_rnd_circuit(self, fx_bg: PCG64, jumps: int) -> None:
         circuit_1 = rand_circuit(nqubits, depth, rng, use_ccx=False)
         p_ref = circuit_1.transpile().pattern
         p_test = p_ref.extract_gflow().to_corrections().to_pattern()
-
+        p_ref.remove_input_nodes()
+        p_test.remove_input_nodes()
         p_ref.perform_pauli_measurements()
         p_test.perform_pauli_measurements()
 
@@ -1204,6 +1216,7 @@ def test_pauli_measurement_end_with_measure(self) -> None:
         p = Pattern(input_nodes=[0])
         p.add(N(node=1))
         p.add(M(node=1, plane=Plane.XY))
+        p.remove_input_nodes()
         p.perform_pauli_measurements()
 
     @pytest.mark.parametrize("backend", ["statevector", "densitymatrix"])
diff --git a/tests/test_pyzx.py b/tests/test_pyzx.py
index 646cd1552..d77217da8 100644
--- a/tests/test_pyzx.py
+++ b/tests/test_pyzx.py
@@ -85,9 +85,11 @@ def test_random_circuit(fx_bg: PCG64, jumps: int) -> None:
     zx_graph = to_pyzx_graph(opengraph)
     opengraph2 = from_pyzx_graph(zx_graph)
     pattern2 = opengraph2.to_pattern()
+    pattern.remove_input_nodes()
     pattern.perform_pauli_measurements()
     pattern.minimize_space()
     state = pattern.simulate_pattern()
+    pattern2.remove_input_nodes()
     pattern2.perform_pauli_measurements()
     pattern2.minimize_space()
     state2 = pattern2.simulate_pattern()
diff --git a/tests/test_tnsim.py b/tests/test_tnsim.py
index 2bfccc1f7..6129c3c05 100644
--- a/tests/test_tnsim.py
+++ b/tests/test_tnsim.py
@@ -334,6 +334,7 @@ def test_with_graphtrans(self, fx_bg: PCG64, jumps: int, fx_rng: Generator) -> N
         pattern = circuit.transpile().pattern
         pattern.standardize()
         pattern.shift_signals()
+        pattern.remove_input_nodes()
         pattern.perform_pauli_measurements()
         state = circuit.simulate_statevector().statevec
         tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
@@ -352,6 +353,7 @@ def test_with_graphtrans_sequential(self, fx_bg: PCG64, jumps: int, fx_rng: Gene
         pattern = circuit.transpile().pattern
         pattern.standardize()
         pattern.shift_signals()
+        pattern.remove_input_nodes()
         pattern.perform_pauli_measurements()
         state = circuit.simulate_statevector().statevec
         tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", graph_prep="sequential", rng=fx_rng)
@@ -401,6 +403,7 @@ def test_evolve(self, fx_bg: PCG64, jumps: int, fx_rng: Generator) -> None:
         pattern = circuit.transpile().pattern
         pattern.standardize()
         pattern.shift_signals()
+        pattern.remove_input_nodes()
         pattern.perform_pauli_measurements()
         state = circuit.simulate_statevector().statevec
         tn_mbqc = pattern.simulate_pattern(backend="tensornetwork", rng=fx_rng)
diff --git a/tests/test_visualization.py b/tests/test_visualization.py
index ccc0a1275..2c6722f79 100644
--- a/tests/test_visualization.py
+++ b/tests/test_visualization.py
@@ -133,6 +133,7 @@ def example_hadamard() -> Pattern:
 
 def example_local_clifford() -> Pattern:
     pattern = example_hadamard()
+    pattern.remove_input_nodes()
     pattern.perform_pauli_measurements()
     return pattern
 
@@ -230,9 +231,9 @@ def test_empty_pattern() -> None:
 # Compare with baseline/test_draw_graph_reference.png
 # Update baseline by running: pytest --mpl-generate-path=tests/baseline
 @pytest.mark.usefixtures("mock_plot")
-@pytest.mark.parametrize("flow_from_pattern", [False, True])
+@pytest.mark.parametrize("flow_and_not_pauli_presimulate", [False, True])
 @pytest.mark.mpl_image_compare
-def test_draw_graph_reference(flow_from_pattern: bool) -> Figure:
+def test_draw_graph_reference(flow_and_not_pauli_presimulate: bool) -> Figure:
     circuit = Circuit(3)
     circuit.cnot(0, 1)
     circuit.cnot(2, 1)
@@ -240,6 +241,9 @@ def test_draw_graph_reference(flow_from_pattern: bool) -> Figure:
     circuit.x(2)
     circuit.cnot(2, 1)
     pattern = circuit.transpile().pattern
-    pattern.perform_pauli_measurements(leave_input=True)
-    pattern.draw_graph(flow_from_pattern=flow_from_pattern, node_distance=(0.7, 0.6))
+    if not flow_and_not_pauli_presimulate:
+        pattern.remove_input_nodes()
+        pattern.perform_pauli_measurements()
+    pattern.standardize()
+    pattern.draw_graph(flow_from_pattern=flow_and_not_pauli_presimulate, node_distance=(0.7, 0.6))
     return plt.gcf()