diff --git a/pytket-mbqc-py/pytket_mbqc_py/__init__.py b/pytket-mbqc-py/pytket_mbqc_py/__init__.py
index 63b01d73..60022e55 100644
--- a/pytket-mbqc-py/pytket_mbqc_py/__init__.py
+++ b/pytket-mbqc-py/pytket_mbqc_py/__init__.py
@@ -5,3 +5,4 @@
 from .cnot_block import CNOTBlocksGraphCircuit
 from .graph_circuit import GraphCircuit
 from .qubit_manager import QubitManager
+from .random_register_manager import RandomRegisterManager
diff --git a/pytket-mbqc-py/pytket_mbqc_py/cnot_block.py b/pytket-mbqc-py/pytket_mbqc_py/cnot_block.py
index 75effc85..dcf4cf50 100644
--- a/pytket-mbqc-py/pytket_mbqc_py/cnot_block.py
+++ b/pytket-mbqc-py/pytket_mbqc_py/cnot_block.py
@@ -23,21 +23,25 @@ def __init__(
         n_physical_qubits: int,
         input_state: Tuple[int],
         n_layers: int,
+        n_registers: int,
     ) -> None:
         """Initialisation method.
 
         :param n_physical_qubits: The maximum number of physical qubits
             available. These qubits will be reused and so the total
             number of 'logical' qubits may be larger.
-        :type n_physical_qubits: int
         :param input_state: Integer tuple describing the input
             to the circuit. This will be a classical binary
             string and so the outcome is deterministic.
-        :type input_state: Tuple[int]
         :param n_layers: The number of layers of CNOT gates.
-        :type n_layers: int
+        :param n_registers: The number of classical registers to use
+            for state preparation information. Note that there should
+            be at least one register per logical qubit.
         """
 
+        # TODO: n_registers should be calculated from the given
+        # input variables.
+
         self.input_state = input_state
         self.n_layers = n_layers
 
@@ -56,7 +60,10 @@ def __init__(
             [[] for _ in range(n_rows)] for _ in range(n_layers)
         ]
 
-        super().__init__(n_physical_qubits=n_physical_qubits)
+        super().__init__(
+            n_physical_qubits=n_physical_qubits,
+            n_registers=n_registers,
+        )
 
         for layer in range(n_layers):
             # If this is the first layer then the control qubit of the first row needs
diff --git a/pytket-mbqc-py/pytket_mbqc_py/graph_circuit.py b/pytket-mbqc-py/pytket_mbqc_py/graph_circuit.py
index c071d92a..0817cf16 100644
--- a/pytket-mbqc-py/pytket_mbqc_py/graph_circuit.py
+++ b/pytket-mbqc-py/pytket_mbqc_py/graph_circuit.py
@@ -5,22 +5,21 @@
 """
 
 from functools import reduce
-from typing import Dict, List, Tuple
+from typing import Dict, List, Tuple, Union
 
 import networkx as nx  # type:ignore
 from pytket import Qubit
 from pytket.circuit.logic_exp import BitLogicExp
 from pytket.unit_id import BitRegister
 
-from pytket_mbqc_py.qubit_manager import QubitManager
+from pytket_mbqc_py.random_register_manager import RandomRegisterManager
 
 
-class GraphCircuit(QubitManager):
+class GraphCircuit(RandomRegisterManager):
     """Class for the automated construction of MBQC computations.
     In particular only graphs with valid flow can be constructed.
     Graph state construction and measurement corrections are added
-    automatically. A child of
-    :py:class:`~pytket_mbqc_py.qubit_manager.QubitManager`.
+    automatically.
 
     :ivar entanglement_graph: Graph detailing the graph state entanglement.
     :ivar flow_graph: Graph describing the flow dependencies of the graph state.
@@ -28,6 +27,10 @@ class GraphCircuit(QubitManager):
     :ivar vertex_measured: List indicating if vertex has been measured.
     :ivar vertex_x_corr_reg: List mapping vertex index to
         the classical register where the required X correction is stored.
+    :ivar vertex_init_reg: List mapping vertex to the register describing
+        the state it was initialised in. In particular this is a 3 bit register
+        with the 0th entry giving the T rotation, the 1st giving the
+        S rotation, and the 2nd giving the Z rotation.
     """
 
     entanglement_graph: nx.Graph
@@ -37,12 +40,16 @@ class GraphCircuit(QubitManager):
     def __init__(
         self,
         n_physical_qubits: int,
+        n_registers: int,
     ) -> None:
         """Initialisation method. Creates tools to track
         the graph state structure and the measurement corrections.
 
         :param n_physical_qubits: The number of physical qubits available.
-        :type n_physical_qubits: int
+        :param n_registers: The number of state initialisation registers
+            to generate. Each register describes the state that the logical
+            qubit is initialised in. Note that the number of such registers
+            should be at least the number of logical qubits.
         """
         super().__init__(n_physical_qubits=n_physical_qubits)
 
@@ -60,6 +67,16 @@ def __init__(
         # vertex has been measured.
         self.vertex_x_corr_reg: List[BitRegister] = []
 
+        # Generate one random register per vertex.
+        # When qubits are added they will be initialised in this
+        # random register. This in except for the case of input qubits
+        # which are initialised in the 0 state, and in which case this
+        # register is overwritten.
+        self.vertex_init_reg = list(
+            self.generate_random_registers(n_registers=n_registers)
+        )
+        self.add_barrier(units=self.qubits)
+
     def get_outputs(self) -> Dict[int, Qubit]:
         """Return the output qubits. Output qubits are those that
         are unmeasured, and which do not have a flow. This should
@@ -101,7 +118,29 @@ def get_outputs(self) -> Dict[int, Qubit]:
 
         # We need to correct all unmeasured qubits.
         for vertex in output_qubits.keys():
-            self._apply_x_correction(vertex=vertex)
+            # Corrections are first applied to invert the initialisation
+            self.T(self.vertex_qubit[vertex], condition=self.vertex_init_reg[vertex][0])
+            self.S(self.vertex_qubit[vertex], condition=self.vertex_init_reg[vertex][0])
+
+            self.S(self.vertex_qubit[vertex], condition=self.vertex_init_reg[vertex][1])
+
+            self.Z(
+                self.vertex_qubit[vertex],
+                condition=(
+                    self.vertex_init_reg[vertex][2]
+                    ^ self.vertex_init_reg[vertex][1]
+                    ^ self.vertex_init_reg[vertex][0]
+                ),
+            )
+
+            # Apply X correction according to correction register.
+            # These are corrections resulting from the measurements
+            self.X(
+                self.vertex_qubit[vertex],
+                condition=self.vertex_x_corr_reg[vertex][0],
+            )
+
+            # Apply Z corrections resulting from measurements.
             self._apply_z_correction(vertex=vertex)
 
         return output_qubits
@@ -114,6 +153,9 @@ def _add_vertex(self, qubit: Qubit) -> int:
 
         :param qubit: Qubit to be added.
         :return: The vertex in the graphs corresponding to this qubit
+
+        :raises Exception: Raised if an insufficient number of initialisation
+            registers were initialised.
         """
         self.vertex_qubit.append(qubit)
         self.vertex_measured.append(False)
@@ -126,6 +168,12 @@ def _add_vertex(self, qubit: Qubit) -> int:
         self.vertex_x_corr_reg.append(x_corr_reg)
         self.add_c_register(register=x_corr_reg)
 
+        if index >= len(self.vertex_init_reg):
+            raise Exception(
+                "An insufficient number of initialisation registers "
+                + "were initialised."
+            )
+
         return index
 
     def add_input_vertex(self) -> Tuple[Qubit, int]:
@@ -140,6 +188,10 @@ def add_input_vertex(self) -> Tuple[Qubit, int]:
         qubit = self.get_qubit()
         index = self._add_vertex(qubit=qubit)
 
+        # In the case of input qubits, the initialisation is not random.
+        # As such the initialisation register should be set to 0.
+        self.add_c_setreg(0, self.vertex_init_reg[index])
+
         return (qubit, index)
 
     def add_graph_vertex(self) -> int:
@@ -151,6 +203,12 @@ def add_graph_vertex(self) -> int:
         self.H(qubit)
         index = self._add_vertex(qubit=qubit)
 
+        # The graph state is randomly initialised based on the
+        # initialisation register.
+        self.T(qubit, condition=self.vertex_init_reg[index][0])
+        self.S(qubit, condition=self.vertex_init_reg[index][1])
+        self.Z(qubit, condition=self.vertex_init_reg[index][2])
+
         return index
 
     @property
@@ -240,8 +298,8 @@ def add_edge(self, vertex_one: int, vertex_two: int) -> None:
         ):
             raise Exception(
                 "This does not define a valid flow. "
-                f"In particular {vertex_two} is the flow of {self.flow_graph.predecessors(vertex_two)}, "
-                f"some of which are measured before {vertex_one}."
+                f"In particular {vertex_two} is the flow of {list(self.flow_graph.predecessors(vertex_two))}, "
+                f"some of which are measured after {vertex_one}."
             )
 
         # If this is the first future of vertex_one then it is taken to be its flow.
@@ -260,34 +318,26 @@ def add_edge(self, vertex_one: int, vertex_two: int) -> None:
             v_of_edge=vertex_two,
         )
 
-    def _apply_x_correction(self, vertex: int) -> None:
-        """Apply X correction. This correction is drawn from
-        the x correction register which is altered
-        by the corrected measure method as appropriate.
-
-        :param vertex: The vertex to be corrected.
-        """
-        self.X(
-            self.vertex_qubit[vertex],
-            condition=self.vertex_x_corr_reg[vertex][0],
-        )
-
-    def _get_z_correction_expression(self, vertex: int) -> BitLogicExp:
+    def _get_z_correction_expression(self, vertex: int) -> Union[None, BitLogicExp]:
         """Create logical expression by taking the parity of
         the X corrections that have to be applied to the neighbouring
-        qubits.
+        qubits. If there are no neighbours then None will be returned.
 
         :param vertex: Vertex to be corrected.
         :return: Logical expression calculating the parity
             of the neighbouring x correction registers.
         """
-        return reduce(
-            lambda a, b: a ^ b,
-            [
-                self.vertex_x_corr_reg[neighbour][0]
-                for neighbour in self.entanglement_graph.neighbors(n=vertex)
-            ],
-        )
+
+        neighbour_reg_list = [
+            self.vertex_x_corr_reg[neighbour][0]
+            for neighbour in self.entanglement_graph.neighbors(n=vertex)
+        ]
+
+        # This happens of this vertex has no neighbours.
+        if len(neighbour_reg_list) == 0:
+            return None
+
+        return reduce(lambda a, b: a ^ b, neighbour_reg_list)
 
     def _apply_z_correction(self, vertex: int) -> None:
         """Apply Z correction on qubit. This correction is calculated
@@ -296,10 +346,12 @@ def _apply_z_correction(self, vertex: int) -> None:
 
         :param vertex: Vertex to be corrected.
         """
-        self.Z(
-            self.vertex_qubit[vertex],
-            condition=self._get_z_correction_expression(vertex=vertex),
-        )
+        condition = self._get_z_correction_expression(vertex=vertex)
+        if condition is not None:
+            self.Z(
+                self.vertex_qubit[vertex],
+                condition=condition,
+            )
 
     def _apply_classical_z_correction(self, vertex: int) -> None:
         """Apply Z correction on measurement result. This correction is calculated
@@ -308,11 +360,13 @@ def _apply_classical_z_correction(self, vertex: int) -> None:
 
         :param vertex: Vertex to be corrected.
         """
-        self.add_classicalexpbox_bit(
-            expression=self.qubit_meas_reg[self.vertex_qubit[vertex]][0]
-            ^ self._get_z_correction_expression(vertex=vertex),
-            target=[self.qubit_meas_reg[self.vertex_qubit[vertex]][0]],
-        )
+        condition = self._get_z_correction_expression(vertex=vertex)
+        if condition is not None:
+            self.add_classicalexpbox_bit(
+                expression=self.qubit_meas_reg[self.vertex_qubit[vertex]][0]
+                ^ condition,
+                target=[self.qubit_meas_reg[self.vertex_qubit[vertex]][0]],
+            )
 
     def corrected_measure(self, vertex: int, t_multiple: int = 0) -> None:
         """Perform a measurement, applying the appropriate corrections.
@@ -321,7 +375,6 @@ def corrected_measure(self, vertex: int, t_multiple: int = 0) -> None:
         :param vertex: Vertex to be measured.
         :param t_multiple: The angle in which to measure, defaults to 0.
             This defines the rotated hadamard basis to measure in.
-        :type t_multiple: int, optional
         :raises Exception: Raised if this vertex has already been measured.
         :raises Exception: Raised if there are vertex in the past of this
             one which have not been measured.
@@ -342,11 +395,68 @@ def corrected_measure(self, vertex: int, t_multiple: int = 0) -> None:
                 + f"are in the future of {vertex} and have already been measured."
             )
 
-        # This is actually optional for graph vertices
-        # as the X correction commutes with the hadamard
-        # basis measurement.
-        self._apply_x_correction(vertex=vertex)
+        # Apply X correction according to correction register.
+        # This is to correct for measurement outcomes.
+        self.X(
+            self.vertex_qubit[vertex],
+            condition=self.vertex_x_corr_reg[vertex][0],
+        )
+
+        self.T(
+            self.vertex_qubit[vertex],
+            # Required to invert random T from initialisation.
+            condition=self.vertex_init_reg[vertex][0],
+        )
+        self.S(
+            self.vertex_qubit[vertex],
+            # Required to invert random T from initialisation.
+            # This additional term is required to account for the case where
+            # the correcting T is commuted through an X correction.
+            condition=(
+                self.vertex_init_reg[vertex][0] & self.vertex_x_corr_reg[vertex][0]
+            ),
+        )
+
+        self.S(
+            self.vertex_qubit[vertex],
+            # Required to invert random T from initialisation.
+            condition=self.vertex_init_reg[vertex][0],
+        )
+
+        self.S(
+            self.vertex_qubit[vertex],
+            # Required to invert random S from initialisation.
+            condition=self.vertex_init_reg[vertex][1],
+        )
+
+        self.Z(
+            self.vertex_qubit[vertex],
+            condition=(
+                # Required to invert random T from initialisation.
+                self.vertex_init_reg[vertex][0]
+                # Required to invert random S from initialisation.
+                ^ self.vertex_init_reg[vertex][1]
+                # Required to invert random Z from initialisation.
+                ^ self.vertex_init_reg[vertex][2]
+                # Required to invert random S from initialisation.
+                # This additional term is required to account for the case where
+                # the correcting S is commuted through an X correction.
+                ^ (self.vertex_init_reg[vertex][1] & self.vertex_x_corr_reg[vertex][0])
+                # Required to invert random T from initialisation.
+                # This additional term is required to account for the case where
+                # the correcting S is commuted through an X correction.
+                ^ (self.vertex_init_reg[vertex][0] & self.vertex_x_corr_reg[vertex][0])
+                # Required to invert random T from initialisation.
+                # This additional term is required to account for the case where
+                # the correcting T is commuted through an X correction.
+                ^ (self.vertex_init_reg[vertex][0] & self.vertex_x_corr_reg[vertex][0])
+            ),
+        )
 
+        # Rotate measurement basis.
+        # TODO: These measurements should be combined with the above
+        # so that the measurement angles are hidden by the initialisation
+        # angles.
         inverse_t_multiple = 8 - t_multiple
         inverse_t_multiple = inverse_t_multiple % 8
         if inverse_t_multiple // 4:
diff --git a/pytket-mbqc-py/pytket_mbqc_py/qubit_manager.py b/pytket-mbqc-py/pytket_mbqc_py/qubit_manager.py
index 8e904372..4ad92cce 100644
--- a/pytket-mbqc-py/pytket_mbqc_py/qubit_manager.py
+++ b/pytket-mbqc-py/pytket_mbqc_py/qubit_manager.py
@@ -14,8 +14,7 @@ class QubitManager(Circuit):
     Manages a collection of qubits. In particular maintains
     a list of qubits which are not in use. This can be added
     to by measuring qubits (making them available again)
-    and drawn from by initialising qubits. This is a child of
-    :py:class:`~pytket.circuit.Circuit`.
+    and drawn from by initialising qubits.
 
     :ivar qubit_list: List of available quits.
     :ivar qubit_initialised: The qubits which have been added
@@ -78,13 +77,26 @@ def physical_qubits_used(self) -> List[Qubit]:
             if initialised
         ]
 
+    @property
+    def initialised_qubits(self) -> List[Qubit]:
+        """Qubits which have been initialised."""
+        return [
+            qubit
+            for qubit, initialised in self.qubit_initialised.items()
+            if initialised
+        ]
+
     def managed_measure(self, qubit: Qubit) -> None:
         """Measure the given qubit, storing the result in the
         qubit's classical register. This will return the qubit to
         the list of available qubits.
 
         :param qubit: The qubit to be measured.
-        :type qubit: Qubit
+
+        :raises Exception: Raised if the qubit to be measured does
+            not belong to the circuit.
         """
+        if not self.qubit_initialised[qubit]:
+            raise Exception(f"The qubit {qubit} has not been initialised.")
         self.qubit_list.insert(0, qubit)
         self.Measure(qubit=qubit, bit=self.qubit_meas_reg[qubit][0])
diff --git a/pytket-mbqc-py/pytket_mbqc_py/random_register_manager.py b/pytket-mbqc-py/pytket_mbqc_py/random_register_manager.py
new file mode 100644
index 00000000..8fc6a694
--- /dev/null
+++ b/pytket-mbqc-py/pytket_mbqc_py/random_register_manager.py
@@ -0,0 +1,107 @@
+"""
+Tools for managing random bits. This include generating random bits
+to dedicated registers.
+"""
+
+from collections.abc import Iterator
+
+from pytket.unit_id import Bit, BitRegister
+
+from pytket_mbqc_py.qubit_manager import QubitManager
+
+
+class RandomRegisterManager(QubitManager):
+    """Class for generating random bits, and managing dedicated registers
+    where they are stored.
+    """
+
+    def generate_random_registers(
+        self,
+        n_registers: int,
+        n_bits_per_reg: int = 3,
+        max_n_randomness_qubits: int = 2,
+    ) -> Iterator[BitRegister]:
+        """Generate registers containing random bits. This is achieved
+        by initialising hadamard basis plus states and measuring them.
+
+        :param n_registers: The number of registers to generate.
+        :param n_bits_per_reg: The number of bits in each register.
+        :param max_n_randomness_qubits: The maximum number of qubits to use
+            to generate randomness. If a number of qubits less than this
+            number are actually available then the number available will
+            be used.
+        :raises Exception: Raised if there are no qubits left to use to
+            generate randomness.
+        :yield: Registers containing random bits.
+        """
+
+        if len(self.qubit_list) == 0:
+            raise Exception(
+                "There are no unused qubits "
+                + "which can be used to generate randomness."
+            )
+
+        # The number of qubits used is the smaller of the maximum number set
+        # by the user, or the number which is available.
+        n_randomness_qubits = min(len(self.qubit_list), max_n_randomness_qubits)
+
+        def get_random_bits(n_random_bits: int) -> Iterator[Bit]:
+            """An iterator over bits populated with random values.
+            These are created by initialising hadamard plus states and
+            measuring the state in the computation basis. At most
+            `n_randomness_qubits` qubit are created in the plus state, then
+            measured and reset as appropriate until the requested number
+            of bits have been created.
+
+            Note that the physical bits may be reused so should be used
+            before the next step of the iteration.
+
+            :param n_random_bits: The number of random bits to generate.
+            :yield: Bits containing random values.
+            """
+            # We repeatedly initialise and measure qubits in groups
+            # of size n_randomness_qubits. This allows randomness generation
+            # to be done in parallel where possible.
+            for _ in range(n_random_bits // n_randomness_qubits):
+                # Initialise all qubits.
+                qubit_list = [self.get_qubit() for _ in range(n_randomness_qubits)]
+
+                # For each qubit, initialise and measure.
+                for qubit in qubit_list:
+                    self.H(qubit=qubit)
+                    self.managed_measure(qubit=qubit)
+                    yield self.qubit_meas_reg[qubit][0]
+
+            # This may leave a number of random bits less than
+            # n_randomness_qubits to generate. We do this here.
+            qubit_list = [
+                self.get_qubit() for _ in range(n_random_bits % n_randomness_qubits)
+            ]
+
+            for qubit in qubit_list:
+                self.H(qubit=qubit)
+                self.managed_measure(qubit=qubit)
+                yield self.qubit_meas_reg[qubit][0]
+
+        # For each bit, copy it's value to a persistent
+        # register. This is done in groups of n_bits_per_reg.
+        for bit_index, bit in enumerate(
+            get_random_bits(n_random_bits=n_bits_per_reg * n_registers)
+        ):
+            # If this is the first bit to copy to a new register, create
+            # the new register.
+            if bit_index % n_bits_per_reg == 0:
+                reg = self.add_c_register(
+                    name=f"rand_{bit_index // n_bits_per_reg}",
+                    size=n_bits_per_reg,
+                )
+
+            # Copy the bit to the persistent register.
+            self.add_classicalexpbox_bit(
+                expression=reg[bit_index % n_bits_per_reg] | bit,
+                target=[reg[bit_index % n_bits_per_reg]],
+            )
+
+            # If this is the last bit to fill the register, yield the register.
+            if bit_index % n_bits_per_reg == n_bits_per_reg - 1:
+                yield reg
diff --git a/pytket-mbqc-py/source/graph_circuit.rst b/pytket-mbqc-py/source/graph_circuit.rst
index 05d998b5..79b97a8d 100644
--- a/pytket-mbqc-py/source/graph_circuit.rst
+++ b/pytket-mbqc-py/source/graph_circuit.rst
@@ -2,17 +2,6 @@ Graph Circuit
 =============
 
 .. automodule:: pytket_mbqc_py.graph_circuit
-
-.. autoclass:: pytket_mbqc_py.graph_circuit.GraphCircuit
-
-    .. automethod:: GraphCircuit.__init__
-
-    .. automethod:: GraphCircuit.add_input_vertex
-    
-    .. automethod:: GraphCircuit.add_graph_vertex
-
-    .. automethod:: GraphCircuit.add_edge
-
-    .. automethod:: GraphCircuit.corrected_measure
-
-    .. automethod:: GraphCircuit.get_outputs
\ No newline at end of file
+    :members:
+    :show-inheritance:
+    :special-members:
\ No newline at end of file
diff --git a/pytket-mbqc-py/source/index.rst b/pytket-mbqc-py/source/index.rst
index ee837f98..942df6a2 100644
--- a/pytket-mbqc-py/source/index.rst
+++ b/pytket-mbqc-py/source/index.rst
@@ -12,6 +12,7 @@ Welcome to pytket-mbqc's documentation!
 
    qubit_manager
    graph_circuit
+   random_register_manager
 
 
 
diff --git a/pytket-mbqc-py/source/qubit_manager.rst b/pytket-mbqc-py/source/qubit_manager.rst
index 3ac356aa..0ac80d25 100644
--- a/pytket-mbqc-py/source/qubit_manager.rst
+++ b/pytket-mbqc-py/source/qubit_manager.rst
@@ -2,13 +2,6 @@ Qubit Manager
 =============
 
 .. automodule:: pytket_mbqc_py.qubit_manager
-
-.. autoclass:: pytket_mbqc_py.qubit_manager.QubitManager
-
-    .. automethod:: QubitManager.__init__
-
-    .. automethod:: QubitManager.get_qubit
-
-    .. automethod:: QubitManager.managed_measure
-
-    .. autoproperty:: QubitManager.physical_qubits_used
\ No newline at end of file
+    :members:
+    :show-inheritance:
+    :special-members:
\ No newline at end of file
diff --git a/pytket-mbqc-py/source/random_register_manager.rst b/pytket-mbqc-py/source/random_register_manager.rst
new file mode 100644
index 00000000..4042191a
--- /dev/null
+++ b/pytket-mbqc-py/source/random_register_manager.rst
@@ -0,0 +1,7 @@
+Random Register Manager
+=======================
+
+.. automodule:: pytket_mbqc_py.random_register_manager
+    :members:
+    :show-inheritance:
+    :special-members:
\ No newline at end of file
diff --git a/pytket-mbqc-py/tests/test_graph_circuit.py b/pytket-mbqc-py/tests/test_graph_circuit.py
index 80cd267a..2e3eda4c 100644
--- a/pytket-mbqc-py/tests/test_graph_circuit.py
+++ b/pytket-mbqc-py/tests/test_graph_circuit.py
@@ -6,7 +6,7 @@
 
 
 def test_plus_state():
-    circuit = GraphCircuit(n_physical_qubits=2)
+    circuit = GraphCircuit(n_physical_qubits=2, n_registers=3)
 
     input_qubit, vertex_one = circuit.add_input_vertex()
     circuit.H(input_qubit)
@@ -36,7 +36,10 @@ def test_plus_state():
 
 
 def test_x_gate():
-    circuit = GraphCircuit(n_physical_qubits=3)
+    circuit = GraphCircuit(
+        n_physical_qubits=3,
+        n_registers=3,
+    )
 
     _, vertex_one = circuit.add_input_vertex()
 
@@ -74,7 +77,10 @@ def test_x_gate():
     [((0, 0), (0, 0)), ((0, 1), (0, 1)), ((1, 0), (1, 1)), ((1, 1), (1, 0))],
 )
 def test_cnot(input_state, output_state):
-    circuit = GraphCircuit(n_physical_qubits=5)
+    circuit = GraphCircuit(
+        n_physical_qubits=5,
+        n_registers=10,
+    )
 
     target_qubit, vertex_one = circuit.add_input_vertex()
     if input_state[1]:
@@ -149,6 +155,7 @@ def test_cnot_block(input_state, output_state, n_layers):
         n_physical_qubits=n_physical_qubits,
         input_state=input_state,
         n_layers=n_layers,
+        n_registers=40,
     )
 
     output_vertex_quibts = circuit.get_outputs()
@@ -181,6 +188,7 @@ def test_large_cnot_block():
         n_physical_qubits=n_physical_qubits,
         input_state=input_state,
         n_layers=n_layers,
+        n_registers=40,
     )
 
     output_vertex_quibts = circuit.get_outputs()
@@ -203,7 +211,7 @@ def test_large_cnot_block():
 
 
 def test_3_q_ghz():
-    graph_circuit = GraphCircuit(n_physical_qubits=5)
+    graph_circuit = GraphCircuit(n_physical_qubits=5, n_registers=5)
 
     input_quibt, input_vertex = graph_circuit.add_input_vertex()
 
@@ -255,7 +263,10 @@ def test_3_q_ghz():
     [((0, 0), (0, 0)), ((0, 1), (0, 1)), ((1, 0), (1, 1)), ((1, 1), (1, 0))],
 )
 def test_cnot_early_measure(input_state, output_state):
-    circuit = GraphCircuit(n_physical_qubits=3)
+    circuit = GraphCircuit(
+        n_physical_qubits=3,
+        n_registers=10,
+    )
 
     target_qubit, vertex_one = circuit.add_input_vertex()
     if input_state[1]:
@@ -314,3 +325,197 @@ def test_cnot_early_measure(input_state, output_state):
     )
 
     assert result.get_counts(output_reg)[output_state] == n_shots
+
+
+def test_2q_t_gate_example():
+    api_offline = QuantinuumAPIOffline()
+    backend = QuantinuumBackend(device_name="H1-1LE", api_handler=api_offline)
+
+    graph_circuit = GraphCircuit(
+        n_physical_qubits=6,
+        n_registers=20,
+    )
+
+    _, input_vertex_0 = graph_circuit.add_input_vertex()
+
+    # H[0]S[0]
+    graph_vertex_0_0 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(input_vertex_0, graph_vertex_0_0)
+    graph_circuit.corrected_measure(input_vertex_0, t_multiple=2)
+
+    # H[0]
+    graph_vertex_0_1 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(graph_vertex_0_0, graph_vertex_0_1)
+    graph_circuit.corrected_measure(graph_vertex_0_0, t_multiple=0)
+
+    _, input_vertex_1 = graph_circuit.add_input_vertex()
+
+    # CZ[0,1]H[1]H[0]
+    graph_vertex_0_0 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(graph_vertex_0_1, graph_vertex_0_0)
+    graph_circuit.corrected_measure(graph_vertex_0_1, t_multiple=0)
+
+    graph_vertex_1_0 = graph_circuit.add_graph_vertex()
+
+    graph_vertex_0_1 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(graph_vertex_0_0, graph_vertex_0_1)
+
+    graph_circuit.add_edge(input_vertex_1, graph_vertex_1_0)
+    graph_circuit.add_edge(graph_vertex_0_0, graph_vertex_1_0)
+
+    graph_circuit.corrected_measure(input_vertex_1, t_multiple=0)
+
+    graph_vertex_1_1 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(graph_vertex_1_0, graph_vertex_1_1)
+
+    # H[0]
+    graph_circuit.corrected_measure(graph_vertex_0_0, t_multiple=0)
+
+    # H[1]
+    graph_circuit.corrected_measure(graph_vertex_1_0, t_multiple=0)
+
+    # H[0]Z[0]
+    graph_vertex_0_0 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(graph_vertex_0_1, graph_vertex_0_0)
+    graph_circuit.corrected_measure(graph_vertex_0_1, t_multiple=4)
+
+    # H[0]
+    graph_vertex_0_1 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(graph_vertex_0_0, graph_vertex_0_1)
+
+    # CZ[0,1]H[0]H[1]
+    graph_vertex_0_2 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(graph_vertex_0_1, graph_vertex_0_2)
+
+    graph_vertex_1_0 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(graph_vertex_1_1, graph_vertex_1_0)
+    graph_circuit.corrected_measure(graph_vertex_1_1, t_multiple=0)
+
+    graph_circuit.corrected_measure(graph_vertex_0_0, t_multiple=0)
+    graph_circuit.corrected_measure(graph_vertex_0_1, t_multiple=0)
+
+    graph_vertex_0_1 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(graph_vertex_0_2, graph_vertex_0_1)
+
+    graph_vertex_1_1 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(graph_vertex_1_0, graph_vertex_1_1)
+
+    graph_circuit.add_edge(graph_vertex_0_2, graph_vertex_1_0)
+
+    # H[0]
+    graph_circuit.corrected_measure(graph_vertex_0_2, t_multiple=0)
+
+    # H[1]
+    graph_circuit.corrected_measure(graph_vertex_1_0, t_multiple=0)
+
+    # H[0]S[0]
+    graph_vertex_0_0 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(graph_vertex_0_1, graph_vertex_0_0)
+    graph_circuit.corrected_measure(graph_vertex_0_1, t_multiple=2)
+
+    # H[0]
+    graph_vertex_0_1 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(graph_vertex_0_0, graph_vertex_0_1)
+    graph_circuit.corrected_measure(graph_vertex_0_0, t_multiple=0)
+
+    outputs = graph_circuit.get_outputs()
+    out_meas_reg = graph_circuit.add_c_register(
+        name="output measure", size=len(outputs)
+    )
+    for qubit, bit in zip(outputs.values(), out_meas_reg):
+        graph_circuit.Measure(qubit=qubit, bit=bit)
+
+    copmiled_graph_circuit = backend.get_compiled_circuit(circuit=graph_circuit)
+    n_shots = 1000
+    result = backend.run_circuit(circuit=copmiled_graph_circuit, n_shots=n_shots)
+    assert result.get_counts(cbits=out_meas_reg)[(1, 0)] == n_shots
+
+
+def test_1q_t_gate_example():
+    ################################
+    # The following compiles to I
+
+    api_offline = QuantinuumAPIOffline()
+    backend = QuantinuumBackend(device_name="H1-1LE", api_handler=api_offline)
+
+    graph_circuit = GraphCircuit(
+        n_physical_qubits=2,
+        n_registers=5,
+    )
+
+    _, input_vertex_0 = graph_circuit.add_input_vertex()
+
+    # H[0]
+    graph_vertex_1 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(input_vertex_0, graph_vertex_1)
+    graph_circuit.corrected_measure(input_vertex_0, t_multiple=0)
+
+    # H[0]T[0]
+    graph_vertex_0 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(graph_vertex_1, graph_vertex_0)
+    graph_circuit.corrected_measure(graph_vertex_1, t_multiple=1)
+
+    # H[0]
+    graph_vertex_1 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(graph_vertex_0, graph_vertex_1)
+    graph_circuit.corrected_measure(graph_vertex_0, t_multiple=0)
+
+    # H[0]T[0]S[0]Z[0]
+    graph_vertex_0 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(graph_vertex_1, graph_vertex_0)
+    graph_circuit.corrected_measure(graph_vertex_1, t_multiple=7)
+
+    outputs = graph_circuit.get_outputs()
+    out_meas_reg = graph_circuit.add_c_register(
+        name="output measure", size=len(outputs)
+    )
+    for qubit, bit in zip(outputs.values(), out_meas_reg):
+        graph_circuit.Measure(qubit=qubit, bit=bit)
+
+    copmiled_graph_circuit = backend.get_compiled_circuit(circuit=graph_circuit)
+
+    n_shots = 1000
+    result = backend.run_circuit(circuit=copmiled_graph_circuit, n_shots=n_shots)
+    assert result.get_counts(cbits=out_meas_reg)[(0,)] == n_shots
+
+    ################################
+    # The following compiles to X
+
+    graph_circuit = GraphCircuit(
+        n_physical_qubits=2,
+        n_registers=5,
+    )
+
+    _, input_vertex_0 = graph_circuit.add_input_vertex()
+
+    # H[0]
+    graph_vertex_1 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(input_vertex_0, graph_vertex_1)
+    graph_circuit.corrected_measure(input_vertex_0, t_multiple=0)
+
+    # H[0]T[0]
+    graph_vertex_0 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(graph_vertex_1, graph_vertex_0)
+    graph_circuit.corrected_measure(graph_vertex_1, t_multiple=1)
+
+    # H[0]
+    graph_vertex_1 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(graph_vertex_0, graph_vertex_1)
+    graph_circuit.corrected_measure(graph_vertex_0, t_multiple=0)
+
+    # H[0]T[0]S[0]Z[0]
+    graph_vertex_0 = graph_circuit.add_graph_vertex()
+    graph_circuit.add_edge(graph_vertex_1, graph_vertex_0)
+    graph_circuit.corrected_measure(graph_vertex_1, t_multiple=3)
+
+    outputs = graph_circuit.get_outputs()
+    out_meas_reg = graph_circuit.add_c_register(
+        name="output measure", size=len(outputs)
+    )
+    for qubit, bit in zip(outputs.values(), out_meas_reg):
+        graph_circuit.Measure(qubit=qubit, bit=bit)
+
+    copmiled_graph_circuit = backend.get_compiled_circuit(circuit=graph_circuit)
+
+    result = backend.run_circuit(circuit=copmiled_graph_circuit, n_shots=n_shots)
+    assert result.get_counts(cbits=out_meas_reg)[(1,)] == n_shots
diff --git a/pytket-mbqc-py/tests/test_random_register_manager.py b/pytket-mbqc-py/tests/test_random_register_manager.py
new file mode 100644
index 00000000..e0867b88
--- /dev/null
+++ b/pytket-mbqc-py/tests/test_random_register_manager.py
@@ -0,0 +1,46 @@
+from itertools import product
+
+from pytket.circuit.display import render_circuit_jupyter
+from pytket.extensions.quantinuum import QuantinuumAPIOffline, QuantinuumBackend
+
+from pytket_mbqc_py import RandomRegisterManager
+
+
+def test_random_register_manager():
+    # Here we test a total of 15 random bits generated on 2 qubits.
+    # Note that this mean one of the qubits must be used more than the other.
+    rand_reg_mngr = RandomRegisterManager(n_physical_qubits=2)
+    n_bits_per_reg = 3
+    n_registers = 5
+    reg_list = list(
+        rand_reg_mngr.generate_random_registers(
+            n_registers=n_registers,
+            n_bits_per_reg=n_bits_per_reg,
+        )
+    )
+
+    # Check that there are 3 registers, and that each has 3 bits.
+    assert len(reg_list) == n_registers
+    assert all(reg.size == n_bits_per_reg for reg in reg_list)
+
+    api_offline = QuantinuumAPIOffline()
+    backend = QuantinuumBackend(device_name="H1-1LE", api_handler=api_offline)
+    n_shots = 100
+
+    compiled_circuit = backend.get_compiled_circuit(rand_reg_mngr)
+    result = backend.run_circuit(
+        circuit=compiled_circuit,
+        n_shots=n_shots,
+        seed=0,
+    )
+    result.get_counts(cbits=reg_list[0])
+
+    # Each bit string in each register should occur with equal probability.
+    # This may fail sometimes as we are performing a statistical check.
+    # With the correct seed it should pass every time.
+    for cbits in reg_list:
+        counts = result.get_counts(cbits=cbits)
+        assert all(
+            abs(counts[bit_string] - (n_shots / (2**n_bits_per_reg))) < (n_shots**0.5)
+            for bit_string in product([0, 1], repeat=n_bits_per_reg)
+        )