diff --git a/.gitignore b/.gitignore
index 5bc66443..b303af32 100644
--- a/.gitignore
+++ b/.gitignore
@@ -18,3 +18,5 @@ graphix/sim/graphix.code-workspace
 .vscode/settings.json
 .pre-commit-config.yaml
 graphix/_version.py
+.coverage
+coverage.xml
diff --git a/benchmarks/efficient_contraction_order_statevec.py b/benchmarks/efficient_contraction_order_statevec.py
new file mode 100644
index 00000000..9e8defec
--- /dev/null
+++ b/benchmarks/efficient_contraction_order_statevec.py
@@ -0,0 +1,147 @@
+"""
+Efficient contraction order statevector simulation of MBQC patterns
+===================================================================
+
+Here we benchmark our efficient contraction order statevector simulator for MBQC,
+which uses the TensorNetwork backend.
+
+The methods and modules we use are the followings:
+    1. :meth:`graphix.pattern.Pattern.simulate_pattern`
+        Pattern simulation with Statevector backend.
+        :meth:`graphix.pattern.Pattern.minimize_space` locally minimizes the space of the pattern.
+    2. :mod:`graphix.sim.tensornet`
+        Pattern simulation with TensorNetwork backend.
+        This enables the efficient contraction order simulation of the pattern.
+        Here we use the `cotengra` optimizer for the contraction order optimization.
+"""
+
+# %%
+# Firstly, let us import relevant modules:
+from __future__ import annotations
+
+from copy import deepcopy
+from time import perf_counter
+from typing import TYPE_CHECKING
+
+import cotengra as ctg
+import numpy as np
+import numpy.typing as npt
+import quimb as qu
+from numpy.random import PCG64, Generator
+
+from graphix.random_objects import get_rand_circuit
+
+if TYPE_CHECKING:
+    import graphix
+
+# %%
+# Next, set global seed and number of thread workers
+GLOBAL_SEED = 2
+qu.core._NUM_THREAD_WORKERS = 1
+
+# %%
+# We then run simulations.
+# Let's benchmark the simulation time of the statevector simulator and the efficient contraction order simulator.
+
+
+def prepare_pattern(n_qubit: int) -> tuple[graphix.pattern.Pattern, graphix.pattern.Pattern]:
+    rng = Generator(PCG64(GLOBAL_SEED))
+    circuit = get_rand_circuit(n_qubit, n_qubit, rng)
+    pattern = circuit.transpile().pattern
+    pattern.standardize()
+    pat_original = deepcopy(pattern)
+    pattern.standardize()
+    pat_original = deepcopy(pattern)
+    return pattern, pat_original
+
+
+def statevector_sim(n_qubit: int, sv_sim: list, pat_original: graphix.pattern.Pattern) -> None:
+    if n_qubit > 17:
+        sv_sim.append(-1)
+    else:
+        start = perf_counter()
+        pat_original.minimize_space()
+        max_sp = pat_original.max_space()
+        sv = pat_original.simulate_pattern(max_qubit_num=max_sp)
+        end = perf_counter()
+        sv_sim.append(end - start)
+        del sv
+
+
+def efficient_contraction_order_sim(eco_sim: list, pattern: graphix.pattern.Pattern) -> None:
+    tn = pattern.simulate_pattern("tensornetwork")
+    start = perf_counter()
+    tn_sv = tn.to_statevector(
+        backend="numpy",
+        optimize=ctg.HyperOptimizer(minimize="combo", max_time=600, progbar=True),
+    )
+    end = perf_counter()
+    eco_sim.append(end - start)
+    del tn_sv, tn
+
+
+def single_iter(n_qubit: int, sv_sim: list, eco_sim: list) -> None:
+    pattern, pat_original = prepare_pattern(n_qubit)
+    statevector_sim(n_qubit, sv_sim, pat_original)
+    efficient_contraction_order_sim(eco_sim, pattern)
+
+
+def benchmark(n_qubit_list: list) -> tuple[npt.NDArray, npt.NDArray]:
+    sv_sim = []
+    eco_sim = []
+    for n_qubit in n_qubit_list:
+        single_iter(n_qubit, sv_sim, eco_sim)
+    return np.array(sv_sim), np.array(eco_sim)
+
+
+# %%
+# Here, we calculate the simulation multiple times so that we can make error bars.
+n_qubit_list = list(range(2, 25))
+sv_sim_results: list[npt.NDArray] = []
+eco_sim_results: list[npt.NDArray] = []
+
+for _ in range(6):
+    sv_sim, eco_sim = benchmark(n_qubit_list)
+    sv_sim_results.append(sv_sim)
+    eco_sim_results.append(eco_sim)
+
+
+sv_sim_means = np.mean(sv_sim_results, axis=0)
+eco_sim_means = np.mean(eco_sim_results, axis=0)
+sv_sim_std = np.std(sv_sim_results, axis=0)
+eco_sim_std = np.std(eco_sim_results, axis=0)
+
+# %%
+# Finally, we save the results to a text file.
+with open("results.txt", "w") as f:
+    f.write("n_qubit, sv_sim_mean, sv_sim_std, eco_sim_mean, eco_sim_std\n")
+    for i in range(len(n_qubit_list)):
+        f.write(f"{n_qubit_list[i]}, {sv_sim_means[i]}, {sv_sim_std[i]}, {eco_sim_means[i]}, {eco_sim_std[i]}\n")
+
+# %%
+# We can now plot the simulation time results.
+import matplotlib.pyplot as plt  # noqa: E402
+
+data = np.loadtxt("results.txt", delimiter=",", skiprows=1)
+n_qubits = data[:, 0].astype(int)
+sv_sim_mean = data[:16, 1].astype(float)
+sv_sim_std = data[:16, 2].astype(float)
+eco_sim_mean = data[:, 3].astype(float)
+eco_sim_std = data[:, 4].astype(float)
+
+fig, ax = plt.subplots(figsize=(8, 5))
+color = "tab:red"
+ax.set_xlabel("Original Circuit Size [qubit]", fontsize=20)
+ax.set_ylabel("Simulation time [sec]", fontsize=20)
+ax.set_yscale("log")
+ax.errorbar(n_qubits[:16], sv_sim_mean, yerr=sv_sim_std, color=color, fmt="o", label="MBQC Statevector (minimizing sp)")
+ax.errorbar(n_qubits, eco_sim_mean, yerr=eco_sim_std, fmt="x", label="MBQC TN base")
+ax.tick_params(axis="y")
+ax.legend(loc="upper left")
+ax.set_title("Simulation time", fontsize=20)
+plt.grid(True, which="Major")
+plt.rcParams["svg.fonttype"] = "none"
+plt.savefig("simulation_time.png")
+plt.close()
+
+# %%
diff --git a/graphix/random_objects.py b/graphix/random_objects.py
index ca5bd62f..e74deb30 100644
--- a/graphix/random_objects.py
+++ b/graphix/random_objects.py
@@ -1,5 +1,8 @@
 from __future__ import annotations
 
+import functools
+from typing import TYPE_CHECKING
+
 import numpy as np
 import scipy.linalg
 from scipy.stats import unitary_group
@@ -7,6 +10,12 @@
 from graphix.channels import KrausChannel
 from graphix.ops import Ops
 from graphix.sim.density_matrix import DensityMatrix
+from graphix.transpiler import Circuit
+
+if TYPE_CHECKING:
+    from collections.abc import Iterable, Iterator
+
+    from numpy.random import Generator
 
 
 def rand_herm(sz: int):
@@ -173,3 +182,116 @@ def rand_Pauli_channel_kraus(dim: int, rank: int | None = None) -> KrausChannel:
     # think we don't really care
 
     return KrausChannel(data)
+
+
+def first_rotation(circuit: Circuit, nqubits: int, rng: Generator) -> None:
+    for qubit in range(nqubits):
+        circuit.rx(qubit, rng.random())
+
+
+def mid_rotation(circuit: Circuit, nqubits: int, rng: Generator) -> None:
+    for qubit in range(nqubits):
+        circuit.rx(qubit, rng.random())
+        circuit.rz(qubit, rng.random())
+
+
+def last_rotation(circuit: Circuit, nqubits: int, rng: Generator) -> None:
+    for qubit in range(nqubits):
+        circuit.rz(qubit, rng.random())
+
+
+def entangler(circuit: Circuit, pairs: Iterable[tuple[int, int]]) -> None:
+    for a, b in pairs:
+        circuit.cnot(a, b)
+
+
+def entangler_rzz(circuit: Circuit, pairs: Iterable[tuple[int, int]], rng: Generator) -> None:
+    for a, b in pairs:
+        circuit.rzz(a, b, rng.random())
+
+
+def generate_gate(
+    nqubits: int,
+    depth: int,
+    pairs: Iterable[tuple[int, int]],
+    rng: Generator,
+    *,
+    use_rzz: bool = False,
+) -> Circuit:
+    circuit = Circuit(nqubits)
+    first_rotation(circuit, nqubits, rng)
+    entangler(circuit, pairs)
+    for _ in range(depth - 1):
+        mid_rotation(circuit, nqubits, rng)
+        if use_rzz:
+            entangler_rzz(circuit, pairs, rng)
+        else:
+            entangler(circuit, pairs)
+    last_rotation(circuit, nqubits, rng)
+    return circuit
+
+
+def genpair(n_qubits: int, count: int, rng: Generator) -> Iterator[tuple[int, int]]:
+    choice = list(range(n_qubits))
+    for _ in range(count):
+        rng.shuffle(choice)
+        x, y = choice[:2]
+        yield (x, y)
+
+
+def gentriplet(n_qubits: int, count: int, rng: Generator) -> Iterator[tuple[int, int, int]]:
+    choice = list(range(n_qubits))
+    for _ in range(count):
+        rng.shuffle(choice)
+        x, y, z = choice[:3]
+        yield (x, y, z)
+
+
+def get_rand_circuit(
+    nqubits: int,
+    depth: int,
+    rng: Generator,
+    *,
+    use_rzz: bool = False,
+    use_ccx: bool = False,
+) -> Circuit:
+    """Generate a random circuit.
+
+    Parameters
+    ----------
+    nqubits : int
+        Number of qubits in the circuit. Must be at least 2.
+    depth : int
+        Number of layers in the circuit.
+    rng : Generator
+        Random number generator.
+    use_rzz : bool, optional
+        Whether to include RZZ gates, by default False.
+    use_ccx : bool, optional
+        Whether to include CCX gates, by default False.
+    """
+    circuit = Circuit(nqubits)
+    gate_choice = (
+        functools.partial(circuit.ry, angle=np.pi / 4),
+        functools.partial(circuit.rz, angle=-np.pi / 4),
+        functools.partial(circuit.rx, angle=-np.pi / 4),
+        circuit.h,
+        circuit.s,
+        circuit.x,
+        circuit.z,
+        circuit.y,
+    )
+    for _ in range(depth):
+        for j, k in genpair(nqubits, 2, rng):
+            circuit.cnot(j, k)
+        if use_rzz:
+            for j, k in genpair(nqubits, 2, rng):
+                circuit.rzz(j, k, np.pi / 4)
+        if use_ccx:
+            for j, k, l in gentriplet(nqubits, 2, rng):
+                circuit.ccx(j, k, l)
+        for j, k in genpair(nqubits, 4, rng):
+            circuit.swap(j, k)
+        for j in range(nqubits):
+            rng.choice(gate_choice)(j)
+    return circuit
diff --git a/graphix/sim/statevec.py b/graphix/sim/statevec.py
index f45bef90..e219cc99 100644
--- a/graphix/sim/statevec.py
+++ b/graphix/sim/statevec.py
@@ -41,6 +41,9 @@ def __init__(
         pr_calc : bool
             whether or not to compute the probability distribution before choosing the measurement result.
             if False, measurements yield results 0/1 with 50% probabilities each.
+        max_qubit_num : int, optional
+            optional argument specifying the maximum number of qubits
+            to be stored in the statevector at a time.
         """
         # check that pattern has output nodes configured
         # assert len(pattern.output_nodes) > 0
@@ -55,9 +58,7 @@ def __init__(
         self.Nqubit = 0
         self.to_trace = []
         self.to_trace_loc = []
-        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.")
+        self.max_qubit_num = _validate_max_qubit_num(max_qubit_num, pattern.max_space())
         super().__init__(pr_calc, rng)
 
         # initialize input qubits to desired init_state
@@ -309,8 +310,13 @@ def evolve(self, op: np.ndarray, qargs: list[int]):
             target qubits' indices
         """
         op_dim = int(np.log2(len(op)))
-        # TODO shape = (2,)* 2 * op_dim
-        shape = [2 for _ in range(2 * op_dim)]
+        shape = (
+            [
+                2,
+            ]
+            * 2
+            * op_dim
+        )
         op_tensor = op.reshape(shape)
         self.psi = np.tensordot(
             op_tensor,
@@ -496,7 +502,7 @@ def expectation_value(self, op, qargs):
         return np.dot(st2.psi.flatten().conjugate(), st1.psi.flatten())
 
 
-def _get_statevec_norm(psi):
+def _get_statevec_norm(psi) -> float:
     """returns norm of the state"""
     return np.sqrt(np.sum(psi.flatten().conj() * psi.flatten()))
 
@@ -514,3 +520,15 @@ def _get_statevec_norm(psi):
         Iterable[graphix.states.State],
         Iterable[numbers.Number],
     ]
+
+
+def _validate_max_qubit_num(max_qubit_num: int | None, max_space: int) -> int | None:
+    if max_qubit_num is None:
+        return
+    if not isinstance(max_qubit_num, int):
+        raise ValueError("max_qubit_num must be an integer")
+    if max_qubit_num < 1:
+        raise ValueError("max_qubit_num must be a positive integer")
+    if max_qubit_num < max_space:
+        raise ValueError("Pattern.max_space is larger than max_qubit_num. Increase max_qubit_num and try again")
+    return max_qubit_num
diff --git a/graphix/sim/tensornet.py b/graphix/sim/tensornet.py
index 03f93f31..1da39956 100644
--- a/graphix/sim/tensornet.py
+++ b/graphix/sim/tensornet.py
@@ -4,6 +4,7 @@
 from copy import deepcopy
 
 import numpy as np
+import numpy.typing as npt
 import quimb.tensor as qtn
 import typing_extensions
 from quimb.tensor import Tensor, TensorNetwork
@@ -467,7 +468,8 @@ def get_basis_coefficient(self, basis, normalize=True, indices=None, **kwagrs):
         normalize (optional): bool
             if True, normalize the coefficient by the norm of the entire state.
         indices (optional): list of int
-            target qubit indices to compute the coefficients, default is the MBQC output nodes (self.default_output_nodes).
+            target qubit indices to compute the coefficients,
+            default is the MBQC output nodes (self.default_output_nodes).
 
         Returns
         -------
@@ -522,28 +524,26 @@ def get_basis_amplitude(self, basis, **kwagrs):
         coef = self.get_basis_coefficient(basis, **kwagrs)
         return abs(coef) ** 2
 
-    def to_statevector(self, indices=None, **kwagrs):
-        """Retrieve the statevector from the tensornetwork.
-        This method tends to be slow however we plan to parallelize this.
-
-        Parameters
-        ----------
-        indices (optional): list of int
-            target qubit indices. Default is the MBQC output nodes (self.default_output_nodes).
+    def to_statevector(self, **kwagrs) -> npt.NDArray:
+        """Take outer product of the tensors in the network and return the statevector.
+        This method uses the contract method and full_simplify method of quimb.tensor.TensorNetwork,
+        which enables the efficient contraction order search and simplification of the tensor network.
+        This method maintains the same computational capabilities as the conventional statevector simulation,
+        as both calculate the same tensor network.
+        The difference is only in the contraction order.
+        This eco simulation is not always exactly optimal, but it is efficient enough for most cases.
+        The current implementation does NOT support intermediate measurements according to the probability distribution.
 
         Returns
         -------
-        numpy.ndarray :
+         : numpy.ndarray
             statevector
         """
-        if indices is None:
-            n_qubit = len(self.default_output_nodes)
-        else:
-            n_qubit = len(indices)
-        statevec = np.zeros(2**n_qubit, np.complex128)
-        for i in range(len(statevec)):
-            statevec[i] = self.get_basis_coefficient(i, normalize=False, indices=indices, **kwagrs)
-        return statevec / np.linalg.norm(statevec)
+        tn = self.copy()
+        tn_simplified = tn.full_simplify("ADCR")
+        output_inds = [self._dangling[str(index)] for index in self.default_output_nodes]
+        psi = tn_simplified.contract(output_inds=output_inds, **kwagrs)
+        return np.array([psi]) if isinstance(psi, (np.complex128, np.float64)) else psi.data.reshape(-1)
 
     def get_norm(self, **kwagrs):
         """Calculate the norm of the state.
@@ -588,7 +588,12 @@ def expectation_value(self, op, qubit_indices, output_node_indices=None, **kwagr
             target_nodes = [output_node_indices[ind] for ind in qubit_indices]
             out_inds = output_node_indices
         op_dim = len(qubit_indices)
-        op = op.reshape([2 for _ in range(2 * op_dim)])
+        op = op.reshape(
+            [
+                2,
+            ]
+            * (2 * op_dim)
+        )
         new_ind_left = [gen_str() for _ in range(op_dim)]
         new_ind_right = [gen_str() for _ in range(op_dim)]
         tn_cp_left = self.copy()
@@ -626,7 +631,8 @@ def evolve(self, operator, qubit_indices, decompose=True, **kwagrs):
             Applied positions of **logical** qubits.
         decompose : bool, optional
             default True
-            whether a given operator will be decomposed or not. If True, operator is decomposed into Matrix Product Operator(MPO)
+            whether a given operator will be decomposed or not.
+            If True, operator is decomposed into Matrix Product Operator(MPO)
         """
         if len(operator.shape) != len(qubit_indices) * 2:
             shape = [2 for _ in range(2 * len(qubit_indices))]
diff --git a/tests/conftest.py b/tests/conftest.py
index e0947f8e..69889d60 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -3,8 +3,8 @@
 import pytest
 from numpy.random import PCG64, Generator
 
+import graphix.random_objects as ro
 import graphix.transpiler
-import tests.random_circuit
 
 SEED = 25
 DEPTH = 1
@@ -34,7 +34,7 @@ def nqb(fx_rng: Generator) -> int:
 
 @pytest.fixture
 def rand_circ(nqb, fx_rng: Generator) -> graphix.transpiler.Circuit:
-    return tests.random_circuit.get_rand_circuit(nqb, DEPTH, fx_rng)
+    return ro.get_rand_circuit(nqb, DEPTH, fx_rng)
 
 
 @pytest.fixture
diff --git a/tests/random_circuit.py b/tests/random_circuit.py
deleted file mode 100644
index ae73579c..00000000
--- a/tests/random_circuit.py
+++ /dev/null
@@ -1,111 +0,0 @@
-from __future__ import annotations
-
-import functools
-from typing import TYPE_CHECKING
-
-import numpy as np
-
-from graphix.transpiler import Circuit
-
-if TYPE_CHECKING:
-    from collections.abc import Iterable, Iterator
-
-    from numpy.random import Generator
-
-
-def first_rotation(circuit: Circuit, nqubits: int, rng: Generator) -> None:
-    for qubit in range(nqubits):
-        circuit.rx(qubit, rng.random())
-
-
-def mid_rotation(circuit: Circuit, nqubits: int, rng: Generator) -> None:
-    for qubit in range(nqubits):
-        circuit.rx(qubit, rng.random())
-        circuit.rz(qubit, rng.random())
-
-
-def last_rotation(circuit: Circuit, nqubits: int, rng: Generator) -> None:
-    for qubit in range(nqubits):
-        circuit.rz(qubit, rng.random())
-
-
-def entangler(circuit: Circuit, pairs: Iterable[tuple[int, int]]) -> None:
-    for a, b in pairs:
-        circuit.cnot(a, b)
-
-
-def entangler_rzz(circuit: Circuit, pairs: Iterable[tuple[int, int]], rng: Generator) -> None:
-    for a, b in pairs:
-        circuit.rzz(a, b, rng.random())
-
-
-def generate_gate(
-    nqubits: int,
-    depth: int,
-    pairs: Iterable[tuple[int, int]],
-    rng: Generator,
-    *,
-    use_rzz: bool = False,
-) -> Circuit:
-    circuit = Circuit(nqubits)
-    first_rotation(circuit, nqubits, rng)
-    entangler(circuit, pairs)
-    for _ in range(depth - 1):
-        mid_rotation(circuit, nqubits, rng)
-        if use_rzz:
-            entangler_rzz(circuit, pairs, rng)
-        else:
-            entangler(circuit, pairs)
-    last_rotation(circuit, nqubits, rng)
-    return circuit
-
-
-def genpair(n_qubits: int, count: int, rng: Generator) -> Iterator[tuple[int, int]]:
-    choice = list(range(n_qubits))
-    for _ in range(count):
-        rng.shuffle(choice)
-        x, y = choice[:2]
-        yield (x, y)
-
-
-def gentriplet(n_qubits: int, count: int, rng: Generator) -> Iterator[tuple[int, int, int]]:
-    choice = list(range(n_qubits))
-    for _ in range(count):
-        rng.shuffle(choice)
-        x, y, z = choice[:3]
-        yield (x, y, z)
-
-
-def get_rand_circuit(
-    nqubits: int,
-    depth: int,
-    rng: Generator,
-    *,
-    use_rzz: bool = False,
-    use_ccx: bool = False,
-) -> Circuit:
-    circuit = Circuit(nqubits)
-    gate_choice = (
-        functools.partial(circuit.ry, angle=np.pi / 4),
-        functools.partial(circuit.rz, angle=-np.pi / 4),
-        functools.partial(circuit.rx, angle=-np.pi / 4),
-        circuit.h,
-        circuit.s,
-        circuit.x,
-        circuit.z,
-        circuit.y,
-    )
-    for _ in range(depth):
-        for j, k in genpair(nqubits, 2, rng):
-            circuit.cnot(j, k)
-        if use_rzz:
-            for j, k in genpair(nqubits, 2, rng):
-                circuit.rzz(j, k, np.pi / 4)
-        if use_ccx:
-            for j, k, l in gentriplet(nqubits, 2, rng):
-                circuit.ccx(j, k, l)
-        for j, k in genpair(nqubits, 4, rng):
-            circuit.swap(j, k)
-        for j in range(nqubits):
-            rng.choice(gate_choice)(j)
-    return circuit
diff --git a/tests/test_generator.py b/tests/test_generator.py
index c6bd5ce3..147c1fc3 100644
--- a/tests/test_generator.py
+++ b/tests/test_generator.py
@@ -7,7 +7,7 @@
 import pytest
 
 import graphix.pauli
-import tests.random_circuit as rc
+import graphix.random_objects as ro
 from graphix.generator import generate_from_graph
 
 if TYPE_CHECKING:
@@ -57,7 +57,7 @@ def test_pattern_generation_flow(self, fx_rng: Generator) -> None:
         nqubits = 3
         depth = 2
         pairs = [(0, 1), (1, 2)]
-        circuit = rc.generate_gate(nqubits, depth, pairs, fx_rng)
+        circuit = ro.generate_gate(nqubits, depth, pairs, fx_rng)
         # transpile into graph
         pattern = circuit.transpile().pattern
         pattern.standardize()
diff --git a/tests/test_gflow.py b/tests/test_gflow.py
index dbe30052..bd86f291 100644
--- a/tests/test_gflow.py
+++ b/tests/test_gflow.py
@@ -19,7 +19,7 @@
     verify_pauliflow,
 )
 from graphix.pattern import Pattern
-from tests.random_circuit import get_rand_circuit
+from graphix.random_objects import get_rand_circuit
 
 if TYPE_CHECKING:
     from collections.abc import Iterable, Iterator
diff --git a/tests/test_pattern.py b/tests/test_pattern.py
index 192a2eee..610f2c40 100644
--- a/tests/test_pattern.py
+++ b/tests/test_pattern.py
@@ -12,9 +12,9 @@
 import graphix.clifford
 import graphix.ops
 import graphix.pauli
+import graphix.random_objects as ro
 import graphix.sim.base_backend
 import graphix.states
-import tests.random_circuit as rc
 from graphix.command import E, M, N
 from graphix.pattern import CommandNode, Pattern
 from graphix.pauli import Plane
@@ -60,7 +60,7 @@ def test_manual_generation(self) -> None:
     def test_standardize(self, fx_rng: Generator) -> None:
         nqubits = 2
         depth = 1
-        circuit = rc.get_rand_circuit(nqubits, depth, fx_rng)
+        circuit = ro.get_rand_circuit(nqubits, depth, fx_rng)
         pattern = circuit.transpile().pattern
         pattern.standardize(method="global")
         assert pattern.is_standard()
@@ -71,7 +71,7 @@ def test_standardize(self, fx_rng: Generator) -> None:
     def test_minimize_space(self, fx_rng: Generator) -> None:
         nqubits = 5
         depth = 5
-        circuit = rc.get_rand_circuit(nqubits, depth, fx_rng)
+        circuit = ro.get_rand_circuit(nqubits, depth, fx_rng)
         pattern = circuit.transpile().pattern
         pattern.standardize(method="global")
         pattern.minimize_space()
@@ -85,7 +85,7 @@ def test_minimize_space_with_gflow(self, fx_bg: PCG64, jumps: int, use_rustworkx
         nqubits = 3
         depth = 3
         pairs = [(i, np.mod(i + 1, nqubits)) for i in range(nqubits)]
-        circuit = rc.generate_gate(nqubits, depth, pairs, rng)
+        circuit = ro.generate_gate(nqubits, depth, pairs, rng)
         pattern = circuit.transpile().pattern
         pattern.standardize(method="global")
         pattern.shift_signals(method="global")
@@ -121,7 +121,7 @@ def test_minimize_space_graph_maxspace_with_flow(self, fx_rng: Generator) -> Non
         for nqubits in range(2, max_qubits):
             depth = 5
             pairs = [(i, np.mod(i + 1, nqubits)) for i in range(nqubits)]
-            circuit = rc.generate_gate(nqubits, depth, pairs, fx_rng)
+            circuit = ro.generate_gate(nqubits, depth, pairs, fx_rng)
             pattern = circuit.transpile().pattern
             pattern.standardize(method="global")
             pattern.minimize_space()
@@ -130,7 +130,7 @@ def test_minimize_space_graph_maxspace_with_flow(self, fx_rng: Generator) -> Non
     def test_parallelize_pattern(self, fx_rng: Generator) -> None:
         nqubits = 2
         depth = 1
-        circuit = rc.get_rand_circuit(nqubits, depth, fx_rng)
+        circuit = ro.get_rand_circuit(nqubits, depth, fx_rng)
         pattern = circuit.transpile().pattern
         pattern.standardize(method="global")
         pattern.parallelize_pattern()
@@ -143,7 +143,7 @@ def test_shift_signals(self, fx_bg: PCG64, jumps: int) -> None:
         rng = Generator(fx_bg.jumped(jumps))
         nqubits = 2
         depth = 1
-        circuit = rc.get_rand_circuit(nqubits, depth, rng)
+        circuit = ro.get_rand_circuit(nqubits, depth, rng)
         pattern = circuit.transpile().pattern
         pattern.standardize(method="global")
         pattern.shift_signals(method="global")
@@ -161,7 +161,7 @@ def test_pauli_measurement_random_circuit(
         rng = Generator(fx_bg.jumped(jumps))
         nqubits = 3
         depth = 3
-        circuit = rc.get_rand_circuit(nqubits, depth, rng)
+        circuit = ro.get_rand_circuit(nqubits, depth, rng)
         pattern = circuit.transpile().pattern
         pattern.standardize(method="global")
         pattern.shift_signals(method="global")
@@ -178,7 +178,7 @@ def test_pauli_measurement_leave_input_random_circuit(
         rng = Generator(fx_bg.jumped(jumps))
         nqubits = 3
         depth = 3
-        circuit = rc.get_rand_circuit(nqubits, depth, rng)
+        circuit = ro.get_rand_circuit(nqubits, depth, rng)
         pattern = circuit.transpile().pattern
         pattern.standardize(method="global")
         pattern.shift_signals(method="global")
@@ -193,7 +193,7 @@ def test_pauli_measurement_opt_gate(self, fx_bg: PCG64, jumps: int, use_rustwork
         rng = Generator(fx_bg.jumped(jumps))
         nqubits = 3
         depth = 3
-        circuit = rc.get_rand_circuit(nqubits, depth, rng, use_rzz=True)
+        circuit = ro.get_rand_circuit(nqubits, depth, rng, use_rzz=True)
         pattern = circuit.transpile(opt=True).pattern
         pattern.standardize(method="global")
         pattern.shift_signals(method="global")
@@ -208,7 +208,7 @@ def test_pauli_measurement_opt_gate_transpiler(self, fx_bg: PCG64, jumps: int, u
         rng = Generator(fx_bg.jumped(jumps))
         nqubits = 3
         depth = 3
-        circuit = rc.get_rand_circuit(nqubits, depth, rng, use_rzz=True)
+        circuit = ro.get_rand_circuit(nqubits, depth, rng, use_rzz=True)
         pattern = circuit.standardize_and_transpile(opt=True).pattern
         pattern.standardize(method="global")
         pattern.shift_signals(method="global")
@@ -228,7 +228,7 @@ def test_pauli_measurement_opt_gate_transpiler_without_signalshift(
         rng = Generator(fx_bg.jumped(jumps))
         nqubits = 3
         depth = 3
-        circuit = rc.get_rand_circuit(nqubits, depth, rng, use_rzz=True)
+        circuit = ro.get_rand_circuit(nqubits, depth, rng, use_rzz=True)
         pattern = circuit.standardize_and_transpile(opt=True).pattern
         pattern.perform_pauli_measurements(use_rustworkx=use_rustworkx)
         pattern.minimize_space()
@@ -408,7 +408,7 @@ def test_get_graph(self, fx_rng: Generator) -> None:
         nqubits = 5
         depth = 4
         pairs = [(i, np.mod(i + 1, nqubits)) for i in range(nqubits)]
-        circuit = rc.generate_gate(nqubits, depth, pairs, fx_rng)
+        circuit = ro.generate_gate(nqubits, depth, pairs, fx_rng)
         pattern = circuit.transpile().pattern
         nodes_ref, edges_ref = pattern.get_graph()
 
@@ -450,7 +450,7 @@ def test_standardize(self, fx_bg: PCG64, jumps: int) -> None:
         rng = Generator(fx_bg.jumped(jumps))
         nqubits = 5
         depth = 4
-        circuit = rc.get_rand_circuit(nqubits, depth, rng)
+        circuit = ro.get_rand_circuit(nqubits, depth, rng)
         pattern = circuit.transpile().pattern
         localpattern = pattern.get_local_pattern()
         localpattern.standardize()
@@ -466,7 +466,7 @@ def test_shift_signals(self, fx_bg: PCG64, jumps: int) -> None:
         rng = Generator(fx_bg.jumped(jumps))
         nqubits = 5
         depth = 4
-        circuit = rc.get_rand_circuit(nqubits, depth, rng)
+        circuit = ro.get_rand_circuit(nqubits, depth, rng)
         pattern = circuit.transpile().pattern
         localpattern = pattern.get_local_pattern()
         localpattern.standardize()
@@ -483,7 +483,7 @@ def test_standardize_and_shift_signals(self, fx_bg: PCG64, jumps: int) -> None:
         rng = Generator(fx_bg.jumped(jumps))
         nqubits = 5
         depth = 4
-        circuit = rc.get_rand_circuit(nqubits, depth, rng)
+        circuit = ro.get_rand_circuit(nqubits, depth, rng)
         pattern = circuit.transpile().pattern
         pattern.standardize_and_shift_signals()
         assert pattern.is_standard()
@@ -507,7 +507,7 @@ def test_mixed_pattern_operations(self, fx_bg: PCG64, jumps: int) -> None:
         ]
         nqubits = 3
         depth = 2
-        circuit = rc.get_rand_circuit(nqubits, depth, rng)
+        circuit = ro.get_rand_circuit(nqubits, depth, rng)
         state_ref = circuit.simulate_statevector().statevec
         for process in processes:
             pattern = circuit.transpile().pattern
@@ -526,7 +526,7 @@ def test_opt_transpile_standardize(self, fx_bg: PCG64, jumps: int) -> None:
         rng = Generator(fx_bg.jumped(jumps))
         nqubits = 5
         depth = 4
-        circuit = rc.get_rand_circuit(nqubits, depth, rng)
+        circuit = ro.get_rand_circuit(nqubits, depth, rng)
         pattern = circuit.transpile(opt=True).pattern
         pattern.standardize(method="local")
         assert pattern.is_standard()
@@ -540,7 +540,7 @@ def test_opt_transpile_shift_signals(self, fx_bg: PCG64, jumps: int) -> None:
         rng = Generator(fx_bg.jumped(jumps))
         nqubits = 5
         depth = 4
-        circuit = rc.get_rand_circuit(nqubits, depth, rng)
+        circuit = ro.get_rand_circuit(nqubits, depth, rng)
         pattern = circuit.transpile(opt=True).pattern
         pattern.standardize(method="local")
         pattern.shift_signals(method="local")
@@ -570,7 +570,7 @@ def test_localpattern_is_standard(self, fx_bg: PCG64, jumps: int) -> None:
         rng = Generator(fx_bg.jumped(jumps))
         nqubits = 5
         depth = 4
-        circuit = rc.get_rand_circuit(nqubits, depth, rng)
+        circuit = ro.get_rand_circuit(nqubits, depth, rng)
         localpattern = circuit.transpile().pattern.get_local_pattern()
         result1 = localpattern.is_standard()
         localpattern.standardize()
diff --git a/tests/test_statevec_backend.py b/tests/test_statevec_backend.py
index f49fbc44..60b8ff8f 100644
--- a/tests/test_statevec_backend.py
+++ b/tests/test_statevec_backend.py
@@ -7,13 +7,26 @@
 import pytest
 
 import graphix.pauli
-from graphix.sim.statevec import Statevec, StatevectorBackend, meas_op
+from graphix.sim.statevec import Statevec, StatevectorBackend, _validate_max_qubit_num, meas_op
 from graphix.states import BasicStates, PlanarState
 
 if TYPE_CHECKING:
     from numpy.random import Generator
 
 
+class TestStatevectorBackend:
+    @pytest.mark.parametrize("max_qubit_num, max_space", [(1.0, 2), (1.1, 2), (0, 2), (1, 2), ("2", 1)])
+    def test_validate_max_qubit_num_fail(self, max_qubit_num: float | int, max_space: int):
+        with pytest.raises(ValueError):
+            _validate_max_qubit_num(max_qubit_num, max_space)
+
+    @pytest.mark.parametrize(
+        "max_qubit_num, max_space, expect", [(None, 2, None), (None, 3, None), (1, 1, 1), (10, 1, 10), (10, 2, 10)]
+    )
+    def test_validate_max_qubit_num_pass(self, max_qubit_num: int | None, max_space: int, expect: int | None):
+        assert _validate_max_qubit_num(max_qubit_num, max_space) == expect
+
+
 class TestStatevec:
     def test_remove_one_qubit(self) -> None:
         n = 10
diff --git a/tests/test_tnsim.py b/tests/test_tnsim.py
index 66ce5f25..81d84bda 100644
--- a/tests/test_tnsim.py
+++ b/tests/test_tnsim.py
@@ -8,7 +8,7 @@
 from numpy.random import PCG64, Generator
 from quimb.tensor import Tensor
 
-import tests.random_circuit as rc
+import graphix.random_objects as ro
 from graphix.clifford import CLIFFORD
 from graphix.command import C, E, X, Z
 from graphix.ops import Ops
@@ -302,7 +302,7 @@ def test_cnot(self, fx_rng: Generator) -> None:
     @pytest.mark.parametrize("jumps", range(1, 11))
     def test_ccx(self, fx_bg: PCG64, jumps: int) -> None:
         rng = Generator(fx_bg.jumped(jumps))
-        circuit = rc.get_rand_circuit(4, 6, rng)
+        circuit = ro.get_rand_circuit(4, 6, rng)
         circuit.ccx(0, 1, 2)
         pattern = circuit.transpile().pattern
         pattern.minimize_space()
@@ -316,7 +316,7 @@ def test_ccx(self, fx_bg: PCG64, jumps: int) -> None:
     @pytest.mark.parametrize("jumps", range(1, 11))
     def test_with_graphtrans(self, fx_bg: PCG64, jumps: int) -> None:
         rng = Generator(fx_bg.jumped(jumps))
-        circuit = rc.get_rand_circuit(4, 6, rng)
+        circuit = ro.get_rand_circuit(4, 6, rng)
         pattern = circuit.transpile().pattern
         pattern.standardize()
         pattern.shift_signals()
@@ -333,7 +333,7 @@ def test_with_graphtrans(self, fx_bg: PCG64, jumps: int) -> None:
     @pytest.mark.parametrize("jumps", range(1, 11))
     def test_with_graphtrans_sequential(self, fx_bg: PCG64, jumps: int) -> None:
         rng = Generator(fx_bg.jumped(jumps))
-        circuit = rc.get_rand_circuit(4, 6, rng)
+        circuit = ro.get_rand_circuit(4, 6, rng)
         pattern = circuit.transpile().pattern
         pattern = circuit.transpile().pattern
         pattern.standardize()
@@ -351,7 +351,7 @@ def test_with_graphtrans_sequential(self, fx_bg: PCG64, jumps: int) -> None:
     @pytest.mark.parametrize("jumps", range(1, 11))
     def test_coef_state(self, fx_bg: PCG64, jumps: int) -> None:
         rng = Generator(fx_bg.jumped(jumps))
-        circuit = rc.get_rand_circuit(4, 2, rng)
+        circuit = ro.get_rand_circuit(4, 2, rng)
 
         pattern = circuit.standardize_and_transpile().pattern
 
@@ -367,7 +367,7 @@ def test_coef_state(self, fx_bg: PCG64, jumps: int) -> None:
     @pytest.mark.parametrize(("nqubits", "jumps"), itertools.product(range(2, 6), range(1, 6)))
     def test_to_statevector(self, fx_bg: PCG64, nqubits: int, jumps: int) -> None:
         rng = Generator(fx_bg.jumped(jumps))
-        circuit = rc.get_rand_circuit(nqubits, 3, rng)
+        circuit = ro.get_rand_circuit(nqubits, 3, rng)
         pattern = circuit.standardize_and_transpile().pattern
         statevec_ref = circuit.simulate_statevector().statevec
 
@@ -380,7 +380,7 @@ def test_to_statevector(self, fx_bg: PCG64, nqubits: int, jumps: int) -> None:
     @pytest.mark.parametrize("jumps", range(1, 11))
     def test_evolve(self, fx_bg: PCG64, jumps: int) -> None:
         rng = Generator(fx_bg.jumped(jumps))
-        circuit = rc.get_rand_circuit(4, 6, rng)
+        circuit = ro.get_rand_circuit(4, 6, rng)
         pattern = circuit.transpile().pattern
         pattern.standardize()
         pattern.shift_signals()
diff --git a/tests/test_transpiler.py b/tests/test_transpiler.py
index d6647cfe..1af9cfc1 100644
--- a/tests/test_transpiler.py
+++ b/tests/test_transpiler.py
@@ -5,8 +5,8 @@
 from numpy.random import PCG64, Generator
 
 import graphix.pauli
+import graphix.random_objects as ro
 import graphix.simulator
-import tests.random_circuit as rc
 from graphix.transpiler import Circuit
 
 
@@ -99,7 +99,7 @@ def test_ccx(self, fx_bg: PCG64, jumps: int) -> None:
         rng = Generator(fx_bg.jumped(jumps))
         nqubits = 4
         depth = 6
-        circuit = rc.get_rand_circuit(nqubits, depth, rng, use_ccx=True)
+        circuit = ro.get_rand_circuit(nqubits, depth, rng, use_ccx=True)
         pattern = circuit.transpile().pattern
         pattern.minimize_space()
         state = circuit.simulate_statevector().statevec
@@ -113,7 +113,7 @@ def test_ccx_opt(self, fx_bg: PCG64, jumps: int) -> None:
         rng = Generator(fx_bg.jumped(jumps))
         nqubits = 4
         depth = 6
-        circuit = rc.get_rand_circuit(nqubits, depth, rng, use_ccx=True)
+        circuit = ro.get_rand_circuit(nqubits, depth, rng, use_ccx=True)
         circuit.ccx(0, 1, 2)
         pattern = circuit.transpile(opt=True).pattern
         pattern.minimize_space()
@@ -125,7 +125,7 @@ def test_transpile_opt(self, fx_rng: Generator) -> None:
         nqubits = 2
         depth = 1
         pairs = [(i, np.mod(i + 1, nqubits)) for i in range(nqubits)]
-        circuit = rc.generate_gate(nqubits, depth, pairs, fx_rng, use_rzz=True)
+        circuit = ro.generate_gate(nqubits, depth, pairs, fx_rng, use_rzz=True)
         pattern = circuit.transpile(opt=True).pattern
         state = circuit.simulate_statevector().statevec
         state_mbqc = pattern.simulate_pattern()
@@ -135,7 +135,7 @@ def test_standardize_and_transpile(self, fx_rng: Generator) -> None:
         nqubits = 3
         depth = 2
         pairs = [(i, np.mod(i + 1, nqubits)) for i in range(nqubits)]
-        circuit = rc.generate_gate(nqubits, depth, pairs, fx_rng, use_rzz=True)
+        circuit = ro.generate_gate(nqubits, depth, pairs, fx_rng, use_rzz=True)
         pattern = circuit.standardize_and_transpile().pattern
         state = circuit.simulate_statevector().statevec
         pattern.minimize_space()
@@ -146,7 +146,7 @@ def test_standardize_and_transpile_opt(self, fx_rng: Generator) -> None:
         nqubits = 3
         depth = 2
         pairs = [(i, np.mod(i + 1, nqubits)) for i in range(nqubits)]
-        circuit = rc.generate_gate(nqubits, depth, pairs, fx_rng, use_rzz=True)
+        circuit = ro.generate_gate(nqubits, depth, pairs, fx_rng, use_rzz=True)
         pattern = circuit.standardize_and_transpile(opt=True).pattern
         state = circuit.simulate_statevector().statevec
         pattern.minimize_space()